An efficient numerical scheme for fully coupled flow and reactive transport in variably saturated porous media including dynamic capillary effects

by   Davide Illiano, et al.

In this paper, we study a model for the transport of an external component, e.g., a surfactant, in variably saturated porous media. We discretize the model in time and space by combining a backward Euler method with the linear Galerkin finite elements. The Newton method and the L-Scheme are employed for the linearization and the performance of these schemes is studied numerically. A special focus is set on the effects of dynamic capillarity on the transport equation.



There are no comments yet.


page 1

page 2

page 3

page 4


Efficient Solvers for Nonstandard Models for Flow and Transport in Unsaturated Porous Media

We study several iterative methods for fully coupled flow and reactive t...

A mathematical model for thermal single-phase flow and reactive transport in fractured porous media

In this paper we present a mathematical model and a numerical workflow f...

A Locally Conservative Mixed Finite Element Framework for Coupled Hydro-Mechanical-Chemical Processes in Heterogeneous Porous Media

This paper presents a mixed finite element framework for coupled hydro-m...

Reactive flow in fractured porous media

In this work we present a model reduction procedure to derive a hybrid-d...

A multi-layer reactive transport model for fractured porous media

An accurate modeling of reactive flows in fractured porous media is a ke...

Space-time upscaling of reactive transport in porous media

Reactive transport in saturated/unsaturated porous media is numerically ...

Accelerated reactive transport simulations in heterogeneous porous media using Reaktoro and Firedrake

This work investigates the performance of the on-demand machine learning...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In this work, we concentrate on efficiently solving reactive transport models in saturated/unsaturated porous media [8, 10]. Such media are observable in the section of the soil closer to the surface where, in the upper part of the domain, we have a coexistence of both water and air phases while, below the water table, the soil becomes fully saturated.

In particular, our model includes dynamic capillarity effects. The capillary pressure is commonly defined as the difference between the pressures of the two phases, in our case, the air and the water. Note that, in the Richards model, the air pressure is set to be equal to zero.

Typically, the capillary pressure is assumed to be a nonlinear decreasing function depending on the water saturation. However, numerous studies are showing that such formulation is often too simplistic and that dynamic effects, due to the changes in time of the water phase, should also be included [2, 3, 5, 11, 13]. Based on this, we consider here the system:


The first equation is the Richards equation, whereas the second is an ordinary differential equation used to include the non-equilibrium effects in the capillary pressure/water content relation. Equilibrium models are obtained for

. Furthermore, the third equation is the reactive transport equation. Here, is the water content, the pressure head, the concentration of the chemical component, the conductivity,

the unit vector in the direction opposite to gravity,

the diffusion/dispersion coefficient, the water flux, the reaction term and finally and are any source terms or external forces involved in the process. Note that where is a nonlinear function depending on and . In the van Genuchten model [4] one has:


is the saturated conductivity and is a soil dependent parameter.

The system (1) is completed by boundary conditions for and , and initial conditions for and .

The rest of the paper is organized as follows: in Section 2 the equations are discretized and linearized. Section 3 includes a numerical example, based on the literature [6], which allows us to compare the different numerical schemes. Finally, Section 4 will conclude this paper with our final remarks.

2 The Numerical Schemes

Applying a Euler implicit time-stepping to (1) gives a sequence of time discrete nonlinear equations. To solve them we apply different linearization schemes: the Newton method, the L-Scheme and a combination of the two [7, 9]. They are here compared thanks to a numerical example inspired by reactive models.

The equations in (1) are fully coupled due to the double dependency of the capillary pressure of both the water content and the concentration . In general, is a function of only , e.g., as presented in [4]. Anyhow, it has been observed [12] that, if an external component is involved, the surface tension becomes a function of the concentration and thus, the capillary pressure itself is influenced by this, i.e. .

In the following, we use the standard notations of functional analysis. The domain , or , is bounded, open and has a Lipschitz continuous boundary . We denote by the space of real-valued, square-integrable functions defined on and its subspace containing the functions having also the first order weak partial derivatives in . is the space of functions belonging to , having zero values on the boundary . We denote by the scalar product and by the associated norm. Finally, assume that is continuous and increasing, is decreasing and .

We now combine the backward Euler method with linear Galerkin finite elements for the discretization of the problem (1). Let be a strictly positive natural number, define the time step size and . Furthermore, is a regular decomposition of , , into -dimensional simplices, with denoting the maximal mesh diameter. The finite element space is defined by


where denotes the space of the afine polynomials on .

The fully discrete Galerkin formulation of the system (1) can be written as:

Problem P(n) Let be fixed. Given , find such that there holds




for all .

Remark 1 We use for the convective term in the transport equation, for simplicity reasons. Nevertheless, all the simulations presented in this paper have also been performed with instead of and the results were almost identical.

In the following, we propose different solving strategies for the system of equations presented above. These strategies are built on the ones discussed in [7], extending them to the case of dynamic capillary pressure (). They are either a monolithic solver of the full system, or a splitting approach obtained by solving first the flow component, using a previously computed concentration, then updating the transport equation, using the newly computed pressure and water content. In both cases, one has to iterate. Each iteration requires solving a non-linear problem, for which, either the Newton methods or the L-Scheme [7, 9, 10] are considered. These strategies are then named: monolithic-Newton scheme (MON-Newton), monolithic-L-Scheme (MON-LS), nonlinear splitting-Newton (NonLinS-Newton) and nonlinear splitting-L-Scheme (NonLinS-LS).

The index denotes the iteration index. As a rule, the iterations start with the solution obtained at the previous time step, for example . This is not necessary for the L-Scheme, which is globally convergent, but it appears to be a natural choice.

2.1 The monolithic Newton method (MON-NEWTON)

The Newton method is a well-known linearization scheme, which is quadratic but only locally convergent. Applying the monolithic Newton method to (4)-(6) leads to

Problem MN(n,j+1) Let be given, find such that




hold true for all .

2.2 The monolithic -scheme (MON-LS)

The monolithic -scheme for solving (4)-(6) reads as

Problem ML(n,j+1) Let be given,
, big enough.
Find such that




hold true for all .

The L-Scheme does not involve the computations of derivatives, the linear systems to be solved within each iteration are better conditioned, compared to the ones given by the Newton method [7, 9], and it is globally (linearly) convergent. The convergence of the scheme has been proved, for the equilibrium model () in [7], and can be easily extended to the non-equilibrium formulation given by the system (10)-(12).

2.3 The splitting approach (NonLinS)

The splitting approach for solving (4)-(6) reads as

Problem S(n,j+1) Let be given, find such that


hold true for all .

Then, with and obtained from the equations above, find such that


holds true for all .

The three equations above can be then linearised using either the Newton method (NonLinS-Newton) or the L-Scheme (NonLinS-LScheme).

2.4 The mixed linearization scheme

It has been already observed, for a different set of equations [9], that combining the Newton method and the L-Scheme can improve the convergence of the scheme. The Newton method is quadratically but only locally convergent and it can produce badly conditioned linearized systems. Moreover, the time step is subject to severe restrictions for guaranteeing the convergence of the scheme, and this has also been observed in numerical examples [1, 7, 9].

Contrarily, the L-Scheme is globally convergent and the linear systems to be solved within each iteration are better conditioned, however, it has only a linear rate of convergence.

The mixed formulation, obtained combining the two schemes, appears to be the best approach and shows practically both global and quadratic convergence. The Newton method commonly fails to converge, if the initial guess is too far from the actual solution. Since this guess is usually the solution at the previous time, this can force restriction on the time step. Instead of reducing the time step one can obtain a better approximation of the initial guess, for the Newton method, by performing few L-Scheme iterations. In the numerical simulation here presented, up to 5 iterations were sufficient to reach a good initial guess for the Newton iteration, which ensured its convergence.

3 Numerical examples

In this section, we use a benchmark problem, from [6], to compare the different linearization schemes and solving algorithms defined above. It describes the recharge of a two-dimensional underground reservoir , in the interval of time . The boundary of the domain and the Dirichlet boundary conditions are defined below.

Furthermore, no flow conditions are imposed on . The initial conditions are given by and . The capillary pressure is defined as , the conductivity is given by (2) and . Finally, the parameters implemented are: , , and the iterations stop whenever all the error norms, and , are below .

Figure 1: Total numbers of iterations for different solvers

We performed the simulations using regular meshes, consisting of squares, with sides . We considered two fixed time steps and .

In Figure 1, we can observe the total numbers of iterations required by the different linearization schemes and solving algorithms. Next to the name of each scheme we report, between parenthesis, which time step has been used.

We can observe, as the Newton method in the monolithic formulation, converges only for coarse meshes, for . For the smaller time step, , it converges for all of the tested meshes.

The L Scheme converges for both time steps, but, since it is linearly convergent, for would require more iterations than the Newton method.

The results obtained thanks to the mixed formulation are particularly interesting. We can observe that this scheme, both in the monolithic and splitting formulation, converges for all the tested meshes also in case of a large time step. Moreover, thanks to the Newton iterations, it appears to be faster than the classical L Scheme. It is as robust as the L Scheme and as fast as the Newton method. For more details regarding the mixed scheme, we refer to [9].

4 Conclusions

In this paper, we considered multiphase flow coupled with a one-component reactive transport in variably saturated porous media, including also the dynamic effects in the capillary pressure. The resulting model is nonlinear and for this reason, three different linearization schemes are investigated: the L-Scheme, the Newton method and a combination of the two. We also studied both monolithic solvers and splitting ones.

The tests show that, for this particular set of equations, the best linearization scheme is the one obtained combining the Newton method and the L-Scheme. Such scheme appears to be both quadratically and globally convergent.

The research of D. Illiano was funded by VISTA, a collaboration between the Norwegian Academy of Science and Letters and Equinor, project number 6367, project name: adaptive model and solver simulation of enhanced oil recovery. The research of I.S. Pop was supported by the Research Foundation-Flanders (FWO), Belgium through the Odysseus programme (project G0G1316N) and Equinor through the Akademia grant.
We thank the members of the Sintef research group and in particular to Dr. Olav Moyner for the assistance with the implementation of the numerical examples in MRST, the toolbox based on Matlab developed at Sintef itself.


  • [1] Cao, X., Pop, I.S.: Uniqueness of weak solutions for a pseudo-parabolic equation modeling two phase flow in porous media, Applied Mathematics Letters, Volume 46, Pages 25-30, (2015).
  • [2] Di Carlo, D.: Experimental measurements of saturation overshoot on infiltration, Water Resources Research, Volume 40, Issue 4, (2004).
  • [3] Fucik, R., Mikyska, J., Sakaki, T., Benes, M., Illangasekare, T.H.: Significance of Dynamic Effect in Capillarity during Drainage Experiments in Layered Porous Media, Vadose Zone Journal 9, Volume 3, (2010).
  • [4] van Genuchten, M.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Science Society of America Journal, Volume 44, Issue 5, Pages 892-898, (1980).
  • [5] Hassanizadeh, S.M., Celia, M.A., Dahle, H.K.: Dynamic Effect in the Capillary Pressure Saturation Relationship and its Impacts on Unsaturated Flow, Vadose Zone Journal, Volume 1, Issue 1, Pages 38-57, (2002).
  • [6] Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P.J., Vachaud, G.: A comparison of numerical simulation models for one-dimensional infiltration, Soil Science Society of America Journal, Volume 41, Pages 285-294, (1977).
  • [7] Illiano, D., Pop, I.S., Radu, F.A.: Iterative schemes for surfactant transport in porous media, arXiv preprint arXiv:1906.00224, (2019).
  • [8] Knabner, P.: Finite element simulation of saturated-unsaturated flow through porous media, LSSC 7, Pages 83-93, (1987).
  • [9] List, F., Radu, F.A.: A study on iterative methods for solving Richards’ equation, Computational Geoscience, Volume 20, Issue 2, Pages 341-353, (2016).
  • [10] Pop, I.S., Radu, F.A., Knabner, P.: Mixed finite elements for the Richards’ equation: linearization procedure, Journal of computational and applied mathematics, Volume 168, Issue 1, Pages 365-373, (2004).
  • [11] Shubao, T., Lei, G., Shun-li, H., Yang, L.: Dynamic effect of capillary pressure in low permeability reservoirs, Petroleum Exploration and Development, Volume 39, Issue 3, Pages 405-411, (2012).
  • [12] Smith, J., Gillham, R.: Effects of solute concentration-dependent surface tension on unsaturated flow: Laboratory sand column experiments, Water Resource Research, Volume 35, Issue 4, Pages 973-982, (1999).
  • [13] Zhuang, L., van Duijn, C.J., Hassanizadeh, S.M.: The effect of dynamic capillarity in modeling saturation overshoot during infiltration, Vadose Zone Journal, Volume 18, Issue 1, Pages 1-14, (2019).