Iterative schemes for surfactant transport in porous media

In this work we consider the transport of a surfactant in a variably saturated porous media. The water flow is modelled by the Richards equations and it is fully coupled with the transport equation for the surfactant. Three linearization techniques are discussed: the Newton method, the modified Picard and the L-scheme. Based on these, monolithic and splitting schemes are proposed and their convergence is analyzed. The performance of these schemes is illustrated on four numerical examples. For these examples, the number of iterations and the condition numbers of the linear systems emerging in each iteration are presented.


Introduction
Many societal relevant problems are involving multiphase flow and multicomponent, reactive transport in porous media. Examples in this sense appear in the enhanced oil recovery, geological CO 2 storage, diffusion of medical agents into the human body, or water or soil pollution. In all these situations, experimental results are difficult to obtain, if not possible at all, and therefore numerical simulations become a key technology. The mathematical models for problems as mentioned above are (fully or partially) coupled, non-linear, possible degenerate partial differential equations. In most cases, finding explicit solutions is not possible, whereas developing appropriate algorithms for finding numerical solutions is a challenge in itself. Here we investigate robust and efficient methods for solving the nonlinear problems obtained after performing an implicit time discretization, the focus being on iterative, splitting or monolithic schemes for fully coupled flow and transport.
Of particular interest here is the special case of multiphase, reactive flow in porous media, namely the surfactant transport in soil [2,19,30,21,24,25]. Surfactants, which are usually organic compounds, are commonly used for actively combating soil and water pollution [16,35,11,40,12]. They contain both hydrophobic and hydrophilic groups and are dissolved in the water phase, being transported by diffusion and convection. Typically, the surfactants are employed in soil regions near the surface (vadose zone), where water and air are present in the pores. Consequently, the outcoming mathematical model accounts the transport of at least one species (the surfactant, but often also the contaminant) in a variably saturated porous medium. Whereas the dependence of the species transport on the flow is obvious, one can encounter the reverse dependece as well since if the surfactants are affecting the interfacial tension between water and air, leading to a dependency of the water flow on the concentration of surfactant. In other words, one has to cope with a fully coupled flow and transport problem, and not only with a one-way coupling, i.e. when only the transport depends on the flow, as mostly considered in reactive transport [32].
Whereas the surfactant transport is described by a reaction-diffusion-convection equation, water flow in variably saturated porous media is modelled by the Richards equation [7,18]. The main assumption in this case is that the air remains in contact with the atmosphere, having a constant pressure (the atmospheric pressure, here assumed zero). This allows reducing the flow model to one equation, the Richards' equation. In mathematical terms, this equation is degenerate parabolic, whose solution has typically low regularity [3].
From the above, and adopting the pressure head as the main unknown in the Richards' equation, we study here different linearization schemes for the model ∂ θ (Ψ, c) ∂t − ∇ · (K(θ (Ψ, c))∇(Ψ + z)) = H 1 (1) and holding for x ∈ Ω (z being the vertical coordinate of x, pointing against gravity) and t ∈ (0, T ]. Here Ω is a bounded, open domain in R d (d = 1, 2 or 3) having a Lipschitz continuous boundary ∂ Ω and T > 0 is the final time. Further, θ (·, ·) denotes the water content, and is a given function depending on the pressure head Ψ and of the surfactant concentration c. Also, K(·) is the hydraulic conductivity, D > 0 the diffusion/dispersion coefficient. Finally, u w is the water flux, R(·) the reaction term expressed as a function of the concentration c, and H 1 , H 2 are the external sinks/sources. Initial and boundary conditions, which are specified below, complete the system. We point out that the water content and the hydraulic conductivity, θ (·, ·) and K(·) are given non-linear functions. They are medium-and surfactant-dependent and are determined experimentally (see [18]). Specific choices are provided in Section 2.
To solve numerically the system (1) -(2) one needs to discretize it in time and space. We refer to [15] for a practical review of numerical methods for the Richards equation. Due to the low regularity of the solution and the need of relatively large time steps, the backward Euler method is a good candidate for the time discretization. Multiple spatial discretization techniques are available, such as the Galerkin Finite Element Method (FEM) [29,5,36], the Mixed Finite Element Method (MFEM) [4,23,33,41,42], the Multi-Point Flux Approximation (MPFA) [22,6,1] and the Finite Volume Method (FVM) [9,13,14].
Since the time discretization is not explicit, the outcome is a sequence of nonlinear problems, for which a linearization step has to be performed. Widely used linearization schemes are the quadratic, locally convergent Newton method and the modified Picard method [10]. For both, the convergence is guaranteed if the starting point is close to the solution. Since for evolution equations the initial gues is typically the solution at the previous time, this induces severe restrictions on the time step (see [34]). Among alternative approaches we mention the L-scheme (see [43,37,31,27] and the modified L-scheme [28], both being robust w.r.t. the mesh size, but converging linearly. In particular, the L-scheme converges for any starting point, and the restriction on the time step, if any, is very mild. The modified L-scheme makes explicit use of the choice of the starting point as the solution obtained at the previous time, and has an improved convergence behavior if the changes in the solutions at two successive times are controlled by the time step. Moreover, the robustness of the Newton method is significantly increased if one considers combinations of the Picard and the Newton methods [8], and in particular of the L-scheme and the Newton scheme [27]. We conclude this discussion by mentioning that in this paper we adopt the FEM and the MPFA, but the iterative schemes presented here can be applied in combination with any other spatial discretization method. The focus is on effectively solving the fully coupled flow and transport system (1) - (2), and in particular on the adequate treating of the coupling between the two model components (the flow and the reactive transport). The schemes are divided in three main categories: monolithic (MON), non-linear splitting (NonLinS) and alternate splitting (AltS). Subsequently, we denote e.g. by MON-NE, the monolithic scheme obtained by applying the Newton method as linearization. The nonlinear splitting schemes (NonLinS) should be understood as solving at each time step first the flow equation until convergence, by using the surfactant concentration from the last iteration, and then with the obtained flow solving the transport equation until convergence. The procedure can be continued iteratively, this being the usual or classical splitting method for transport problems. The convergence of NonLinS does not depend on the linearization approach used for each model component (Newton, Picard or L−scheme), because we assume that the nonlinear subproblems are solved exactly, i.e. until convergence. Finally, the alternate splitting methods (AltS) have a different philosophy. Instead solving each subproblem until convergence within each iteration, one performs only one step of the chosen linearization. For example, AltS-NE will perform one Newton step for each model component, and iterate. These schemes are illustrated in Figures 1, 2.
All the schemes can be analysed theoretically, and we do this exemplary for MON-LS, i.e. for the monolithic approach combined with the L-scheme. Based on comparative numerical tests performed for academic and benchmark problems, we see that the alternate methods can save substantial computational time, while maintaining the robustness of the L-scheme.
The remaining of the paper is organized as follows. In Section 2 we establish the mathematical model and the notation used and present the iterative monolithic and splitting schemes. In Section 3 we prove the convergence of the MON − LS scheme and briefly discuss the convergence of the other schemes. Section 4 presents four different numerical examples. They are inspired by the cases already studied in the literature [27,24]. Section 5 concludes this work.

Problem formulation, discretization and iterative schemes
We solve the fully coupled system (1)-(2), completed by homogeneous Dirichlet boundary conditions for both Ψ and c and the initial conditions: We use the van Genuchten-Mualem parameterization [17] where θ r and θ s represent the values of the residual and saturated water content, K s is the conductivity and α and n are model parameters depending on the soil. Observe that in the expression above for θ the influence of the surfactant on the water flow is neglected. As reported in [20,24,39], the surface tension between water and air does depend on the surfactant concentration c, implying the same for the function θ above. The following parametrization is proposed in [24] θ (Ψ, c) .
Here γ and γ 0 are the surface tensions at concentrations c and c 0 , the second being a reference concentration. The parameters a and b depend on the fluid and the medium. We refer to [38,39] for details about the scaling factor in (5).
This gives the following expressions for θ and K This shows that the flow component also depends on the reactive transport, implying that the model is coupled in both way.
In the following we proceed by discretizing the equations (1) and (2) in time and space. We will use common notations in functional analysis. We denote by L 2 (Ω) the space of real valued, squared integrable function defined on Ω and H 1 (Ω) its subspace, containing the functions having also the first order derivatives in L 2 (Ω). H 1 0 (Ω) is the space of functions belonging to H 1 (Ω) and vanishing on ∂ Ω. Further, we denote by < ·, · > the L 2 (Ω) scalar product (and by · the associated norm) or the pairing between H1 0 and its dual H −1 . Finally, by L 2 (0, T ; X) we mean the Bochner space of functions taking values in the Banach-space X, the extension to H 1 (0, T ; X) being straightforward.
With this we state the weak formulation of the problem related to (1) -(2): and hold for all Φ 1 , Φ 2 ∈ H 1 0 (Ω) and almost every t ∈ (0, T ]. We now combine the backward Euler method with linear Galerkin finite elements for the discretization of Problem P. We let N ∈ N be a strictly positive natural number and the time step τ := T /N. Correspondingly, the discrete times are t n . = nτ (n ∈ 0, 1, . . . , N). Further, we let T h be a regular decomposition of where P 1 (T ) denotes the space of linear polynomials on T and v h|T the restriction of v h to T . For the fully discrete counterpart of Problem P we let n ≥ 1 be fixed and assume that Ψ n−1 h , c n−1 h ∈ V h are given. The solution pair at time t n solves and e z denotes the unit vector in the direction opposite to gravity.

Remark 1
Observe that u n−1 w appears in the convective term in (12). This choice is made for the ease of presentation. Nevertheless, all calculations carried out in this paper were doubled by ones where u n w has replaced u n−1 w . The differences in the results were marginal.
Observe that Problem P n is a coupling system of two elliptic, nonlinear equations. In the following we discuss different iterative schemes for solving this system.

Iterative linearization schemes
We discuss monolithic and splitting approaches for solving Problem P n , combined with either the Newton-method, or the modified Picard [10] or the L-scheme [31,27]. In the following the index n always refers to the time step, whereas j denotes the iteration index. As a rule, the iterations start with the solution at the last time, t n−1 .
In the monolithic approach one solves the two equations of the system (11)-(12) at once, and combined with a linearization method. Formally, this becomes Problem PMon n, j+1 : Find Ψ n, j+1 and c n, j+1 such that F Lin k is a linearization of the expression F k (k = 1, 2) appearing in the system (11)- (12). Depending on the used linearization technique, one speaks about a monolithic-Newton scheme (Mon-Newton), or monolithic-Picard (Mon-Picard) or monolithic-L-scheme (Mon-LS). These three schemes will be presented in detail below.
In the iterative splitting approach one solves each equation separately and then iterates between these using the results obtained previously. We distinguish between two main splitting ways: the nonlinear slitting and the alternate splitting. This is schematized in Figure 1 and Figure 2 respectively. The former becomes Problem PNonLinS n, j+1 : Find Ψ n, j+1 and c n, j+1 such that For the linearization of F 1 and F 2 one can use one of the three linearization techniques mentioned before. In contrast, in the alternate splitting one performs only one linearization step per iteration, see also Figure 2. The alternate splitting scheme becomes Problem PAltS n, j+1 : Find Ψ n, j+1 and c n, j+1 such that Depending which linearization is used, one speaks about alternate splitting Newton (AltS-NE) or alternate splitting L-scheme (AltS-LS). Both schemes are presented in detail below.

The monolithic Newton method (Mon-Newton)
We recall that Newton scheme is quadratically, but only locally convergent. The monolithic Newton method applied to (11)- (12) gives Problem PMon-Newton n, j+1 : and

The monolithic Picard (Mon-Picard)
The modified Picard method was initially proposed by Celia [10] for the Richards equation. It is similar to the Newton method in dealing with the nonlinearity in the saturation, but not in the permeability. Being a modification of the Newton method, modified Picard method is only linearly convergent [34]. The monolithic Picard method applied to (11)-(12) becomes Problem PMon-Picard n, j+1 : and The equations (18) and (19) have been expressed as functions of only the unknown pressure Ψ n, j+1 h and concentration c n, j+1 h , respectively. To achieve this, in the former we used only the derivative of θ with respect to Ψ and only the derivative of θ with respect to c in the latter.
Alternatively, both of the partial derivatives can be involved in the formulation,

The monolithic L-scheme (Mon-LS)
The monolithic L-scheme for solving (8)(9) becomes L 1 and L 2 are free to be chosen parameters but should be large enough to ensure the convergence of the scheme, see Sec. 3. In practice, the values of L 1 , L 2 are The L-scheme does not involve the computations of derivatives, and the linear systems to be solved within each iteration are better conditioned compared to the ones given by Newton or Picard method (see [27]). Moreover, this scheme is (linearly) convergent for any initial guess for the iteration.

The non-linear splitting approach (NonLinS)
The non-linear splitting approach for solving (8)(9) becomes Problem PNonLinS n, j+1 : As for the monolithic schemes, one can apply the different linear iterative schemes to obtain fully linear versions of the splitting approach. This is done first to solve (23) and, once a solution to (23) is available, this is employed in the linearization of (24).

The alternate L-scheme (AltS-LS)
The alternate L-scheme for solving (8)(9) becomes hold true for all v h ∈ V h . Then, with Ψ n, j+1 h obtained above, find c n, j+1 Remark 2 (Stopping criterion) For either monolithic or splitting schemes, one stops the iteration process whenever where ε 1 , ε 2 are small numbers. Here we took 10 −07 or 10 −08 .

Convergence analysis
In this section we analyse the convergence of the monolithic L-scheme introduced through Problem PMon-LS n, j+1 . We restrict the analysis to this iteration, but mention that the convergence analysis for the other (monolithic or splitting) schemes introduced above, can be done in a similar fashion. We start by defining the errors e j+1 obtained at iteration j + 1. The scheme is convergent if both errors vanish when j → ∞. The convergence is obtained under the following assumptions: (A1) There exist α Ψ > 0 and α c ≥ 0 such that for any Ψ 1 , Ψ 2 ∈ R and c 1 , Furthermore, there exist two constants θ m ≥ 0 and θ M < ∞ such that θ m ≤ θ ≤ θ M .
(A2) The function K(θ (·, ·)) is Lipschitz continuous, with respect to both variables, and there exist two constants K m and K M such that 0 ≤ K m ≤ K ≤ K M < ∞.
Remark 3 (A2) is satisfied in most realistic situations. (A3) is a pure technical one, being satisfied when data is sufficiently regular, which is assumed to be the case for the present analysis. The inequality (32) in (A1) is a coercivity assumption.
It is in particular satisfied if Θ only depends on Ψ, and for common relationships Θ − −Ψ encountered in the engineering literature.

Remark 4
The convergence rate resulting from (44) do not depend on the spatial mesh size. Also observe that this convergence is obtained for any initial guess. Based on this, the method is globallly convergent, which is in contrast to the Newton or (modified) Picard schemes, converging only locally. It can be observed that, the larger the time step, the smaller the constants L 1 and L 2 are, resulting in a faster convergence. For small steps instead the convergence rate can approach 1. On the other hand, if the time step is small enough, one may reach the regime where the Newton scheme becomes convergent (see [34]). Alternatively, one may first perform a number of L-scheme iterations, and use the resulting as an initial guess for the Newton scheme (see [27]), or consider the modified L-scheme in [28]. In either situations, the convergence behavior was much improved.

Remark 5
The convergence of the modified Picard and Newton method applied to the Richards equation has been already proved in [34]. Such results can be extended to the coupled problems considered here.

Numerical examples
In this section we consider four test cases for the proposed linearization schemes, inspired by the literature [27,24]. The schemes have been implemented in the open source software package MRST [26], an open source toolbox based on Matlab, in which multiple solvers and models regarding flows in porous media are incorporated.

Example 1A: flow and transport in strictly unsaturated media
We start our numerical studies with a manufactured problem admitting an analytical solution [27]. The unit square Ω is divided into two sub-domains:   Table 1.
Further, for both Richards and transport equations, we have a source term, f (x, y) = .006 cos(4/3πy) sin(x), on Ω up . No external forces or sources, are defined in the lower region, i.e. f = 0 on Ω down . Finally, the initial condition for the concentration is given by c 0 = 1 and the boundary conditions by c |Γ D = 4.
We performed simulations using regular meshes, consisting of squares, whose sides were of length dx = [1/10, 1/20, 1/40, 1/80]. We consider also varying time steps of sizes τ = [1/10, 1/20, 1/40, 1/80]. In Fig. 3a we are plotting the pressure and concentration profiles at the final time T = 1. We point out that in this first example we are always in the strictly unsaturated regime, implying that Richards equation is a regular. All the proposed iterative schemes were converging for this example. In Fig. 4 is given the total number of iterations for the different schemes.
More details regarding the total number of iterations and the condition number of the linear systems are presented in Tables 2, 3. The condition number is computed at the first iteration of each algorithm and with respect to the Euclidean norm. In Table 2, we fixed a time step τ = 1/10 and we investigated different mesh sizes, precisely dx = [1/10, 1/20, 1/40, 1/80]. In Table 3 we use a constant dx = 1/40 and varying the time step sizes τ = [1/10, 1/20, 1/40, 1/80]. We point out that the alternate schemes are converging much faster than the classical ones. We also remark the high differences in the condition numbers, the L−scheme based algorithms being much better conditioned.

Example 1B: flow and transport in variably saturated porous media
For the second example we use the same domain, mesh sizes, boundary conditions and parameters, but we allow a saturated/unsaturated regime by changing the initial condition for the pressure. We consider a subdivision of p 0 between upper and lower regions, precisely: p 0 up = −2 and p 0 down = −y + 1/4. This new expression for p 0 down gives a positive pressure in the lower part of the domain (saturated region    At the iteration j + 1, the term R(c) is linearized in the following way: In Fig. 10a we show again the pressure and concentration profiles at the final step T = 1. The main differences to the previous example, i.e. Figure (3a) are in the values of the pressure. We can observe again a discontinuity in the pressure profile but, more importantly, it is evident a jump from negative to positive values. Such results were expected considering the initial pressure imposed on the domain.  In Fig. 6 are presented the total number of iterations. We remark that in this case only the L−scheme based algorithms are converging. It is also interesting to observe that the difference in the number of iterations between the more commonly used non-linear splitting approach (NonLinS) and alternate splitting (AltLinS) approach.
The alternate method appears to be a valid alternative to the common formulation. It produces equally accurate results, requiring fewer iterations.
As for the previous example we present in the Tables 4, 5 the precise numbers of iterations and condition numbers for each algorithm implemented for different mesh diameters and time steps. Each segment (−) in the tables below, implies that the method failed to converge for such particular combination of time step and space mesh. As already observed in Fig. 6 the L-scheme based solvers are the only ones converging in all cases. Moreover, the linear systems associated with the L-scheme are better conditioned than the ones for Picard or Newton methods. We finally remark that, as expected, for smaller time steps the Newton and Picard schemes converges, see Table 5.

Example 2A: well in unsaturated porous media
Our next example is inspired from [24]. We consider same domain (e.g. the unit square), boundary and initial condition and parameters as in the first numerical example (1A). The medium is again strictly unsaturated. We include now, in the upper part of the domain, a well and inject water with a specific concentration of the external component. No analytical solution is available for this case. Due to the higher complexity of the problem we use more refined meshes, precisely dx = [1/50, 1/100, 1/150, 1/200]. The pressure at the well is set to p W = −10 and the concentration of the surfactant to c W = 10.
In Fig. 7, we present the different profiles of pressure and concentration at the   Once more, in Fig. 8 we compare the different solving algorithms. We study the numbers of iterations and the conditions numbers of the linearized systems. As for the first example, the media being unsaturated, the Richards equation does not degenerate and all the schemes converge. We can observe, in the Tables 6, 7, that the monolithic Newton method is the fastest, in term of numbers of iterations. We remark that the alternate splitting approach (AltLinS), once more, requires fewer iterations than the classical splitting algorithm (NonLinS) for all of the linearization schemes. The linear systems resulting by applying the L-scheme based solvers are better conditioned compared with the other solvers.

Example 2B: well in variably saturated porous media
Our last numerical example is obtained by changing the initial condition for pressure in the example 2A. We use the same p 0 as in example 1B. The profiles of pressure and concentration at the beginning and end of the simulation, i.e. at T = 1 hour, are presented in Fig. 11. We can observe smaller changes, compared to the previous example, due to a smaller time interval (1 hour versus 1 day). In Fig. 11 we present the total number of iterations for the different schemes applied

Conclusions
In this paper we considered surfactant transport in variably saturated porous media. The water flow and the transport are in this case fully coupled. Three linearization techniques were considered: the Newton method, the modified Picard and the L-scheme. Based on these, monolithic and splitting schemes were designed, analyzed and numerically tested. We conclude that the only quadratic convergent scheme is the monolithic Newton, that the L-scheme based solvers are the most robust ones and produce well-conditioned linear systems and that the alternative schemes are faster than the classical splitting approaches.