Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods

Advances in Water Resources 32 (2009) 340–352

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods Anis Younes *, Marwan Fahs, Selim Ahmed Institut de Mécanique des Fluides et des Solides, UMR 7507, 2 rue Boussingault, 67000 Strasbourg, France

a r t i c l e

i n f o

Article history: Received 9 July 2008 Received in revised form 12 November 2008 Accepted 19 November 2008 Available online 30 November 2008 Keywords: Density driven flow Mixed finite elements Discontinuous finite elements MPFA Method of lines ODE BDF Time error

a b s t r a c t Modelling density driven flow problems requires an excessive computational time and/or heavy equipments due to the non-linear coupling between flow and transport equations. In this work, we develop a robust numerical model with efficient advanced approximations for both spatial and temporal discretizations in order to reduce the excessive computational requirement while maintaining accuracy. The spatial discretization is based upon the lumped formulation of the mixed finite element method for the flow equation, the discontinuous Galerkin method for advection and the multipoint flux approximation method for dispersion. The method of lines (MOL) is used to allow higher-order temporal discretization of the coupled flow-transport system. The developed model adapts in both the order of approximation and time step to provide the necessary accuracy. Standard numerical test cases as well as a laboratory scale experiment with a high viscosity contrast between injected and displaced fluids are simulated to show drawbacks of the standard approach and to highlight the efficiency, accuracy and robustness of the developed model for density driven flow problems. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Density driven flow plays a significant role in water resources management and engineering [53,55,69]. The need for accurate models to simulate the transport of salt water continues to increase due to a large number of environmental problems such as saltwater intrusion in exploited coastal aquifers [52] or aquifers overlying salt formations, leakage from landfills, and disposal of radioactive or toxic wastes in salt rock formations [45]. Density effects can also alter tracer tests’ interpretations which can lead to erroneous estimates of transport parameters [8,68]. Modeling density driven flow problems requires a coupled flow-transport numerical model where the Darcy’s flow equation and the advection dispersion transport equation are coupled by the state equations linking density and viscosity variations to mass fraction. Due to these strong nonlinearities, the simulation of density driven flow problems requires meshes of several millions of nodes and heavy computational equipments [42] even for simulations at the laboratory scale [35]. In this work, we develop a robust numerical model with efficient advanced approximations for both spatial and temporal discretizations in order to reduce the large computational requirement while maintaining the accuracy.

* Corresponding author. Fax: +33 388 614 300. E-mail address: [email protected] (A. Younes). 0309-1708/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2008.11.003

For the spatial discretization of flow and transport equations, we use numerical methods that are specifically suited to achieve high accuracy for each kind of equation. Indeed, the mixed finite element (MFE) method is used for the flow equation since it is locally conservative and produces accurate and consistent velocity field even for highly heterogeneous domains [18]. For the transport equation, the discontinuous Galerkin (DG) method is used to discretize the advection equation and combined with the multipoint flux approximation (MPFA) method for the discretization of the dispersion equation [72]. The DG method combines the advantageous of the finite element (FE) and finite volume (FV) methods to obtain a robust and accurate numerical scheme for problems involving sharp fronts [59]. Compared to other standard Eulerian schemes, the DG method allows handling complicated geometries, defining strategies for grid refinement or un-refinement and changing the degree of approximation from one element to the other. On the other hand, the MPFA method is locally conservative and handle general irregular grids on anisotropic heterogeneous domains [1,2,19,39,67]. The MPFA method uses the same type of unknowns (average value per element) than the DG method and therefore, both discretizations can be gathered into one system matrix [72]. The spatial discretization (MFE_DG_MPFA) based upon the combination of MFE, DG and MPFA methods was shown to be robust and accurate for modelling density driven flow problems [4]. On the other hand, the temporal discretization has received relatively little attention in the literature of density driven flow problems and was often limited to the first order backward Euler time

341

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

discretization. The linearization technique is often based upon the Picard iteration method since the full Newton Raphson method can be complex for such applications [56]. Moreover, the time step management is usually based upon the number of iterations required to reach the convergence and does not explicitly measure temporal truncation errors of the solution which may lead to unreliable and/or inaccurate results. In this work, we use the method of lines (MOL) for the coupled flow-transport system. The MOL is a technique for improving temporal accuracy and efficiency when solving partial differential equations (PDEs) by the use of adaptive higher-order time discretizations with formal error control [22]. With the MOL, we discretize first spatial derivatives and then integrate in time the semidiscrete problem as a system of ordinary differential equations (ODEs) or differential algebraic equations (DAEs). The temporal accuracy can be specified by the user, therefore the error checking, robustness, order selection and time step adaptativity features available in sophisticated ODE/DAE codes can be applied to the time integration of the PDEs [62]. In the context of flow in porous media, the MOL was used to solve saturated and unsaturated flow problems [22,36,47,48,62] and was shown to be very effective for the resolution of the nonlinear Richards’ equation [23,24,36,43,47,48,62,63]. In this paper, we show how the MOL can be used to solve density driven flow problems. Our objective is to formulate a numerical solution based upon the MFE_DG_MPFA spatial discretization and a higher-order method for time integration. The scheme is developed for unstructured triangular meshes. These meshes are suitable for practical problems with complex geometry and local mesh refinement. Accuracy and robustness of the scheme are evaluated with standard numerical test cases and with a laboratory scale experiment performed with a high viscosity contrast between injected and displaced fluids. 2. Mathematical models The most common mathematical models of brine transport in porous media are based upon the mass conservation equations and Darcy’s law. Following Huyakorn et al. [34], the flow system can be written in terms of equivalent fresh water head:

@h @ q @C þe þ qr  q ¼ 0; @t @C@t  qg q  q0 rz ; q ¼  0 k rh þ

qS

l

q0

ð1Þ ð2Þ

where q is the fluid density [ML3], S the specific mass storativity related to head changes [L1], h the equivalent freshwater head [L], t the time [T], e the porosity [–], C the solute mass fraction [M salt/M fluid], q the Darcy’s velocity [L T1], q0 the density of the displaced fluid [ML3], g the gravity acceleration [L T2], l the fluid dynamic viscosity [ML1 T1], k the permeability tensor [L2] and z the depth [L]. The solute mass conservation is written in term of mass fraction:

@ðeqCÞ þ r  ðqCq  qD  rCÞ ¼ 0; @t

ð3Þ

q ¼ q0 þ ðq1  q0 ÞC and l ¼ l0

ð4Þ

with aL and aT the longitudinal and transverse dispersivities [L], Dm the pore water diffusion coefficient [L2T1] and I the unit tensor. The associated boundary conditions of the flow-transport system (1)–(3) are of Dirichlet, Neuman or mixed type.





l1 C ; l0

ð5Þ

where q1 and l1 are respectively, density and viscosity of the injected fluid and l0 the viscosity of the displaced fluid. Different state equations may be used for density or viscosity [16]. 3. Spatial discretization To solve the set of Eqs. (1)–(5), we favour advanced numerical methods which allow handling complicated geometries with unstructured grids and which lead to accurate and consistent evaluation of flow fluxes and solute mass fraction. We implement the spatial discretization based upon the MFE method for the flow equation, the upwind DG method for advection and the MPFA method for dispersion. 3.1. Spatial discretization of the flow The idea of the MFE method is to approximate simultaneously the equivalent freshwater head h and the velocity q. The method is locally conservative and leads to accurate and consistent evaluation of the velocity field. In highly heterogeneous domains, MFE results’ were shown to be more accurate than conforming finite element [32,49] and control volume finite element results [18]. The hybridization technique is often used with MFE to reduce the unknowns to the pressure at the edges or faces and to obtain a positive definite linear system [10,11]. The MFE method has been successfully combined with the MOL in [22,23]. Note that this combination is much more involved than the combination of the MOL with FDs. Indeed, the hybridization process requires a discrete approximation in time [22] whereas the time derivative should be continuous with the MOL. To avoid these difficulties and to remove unphysical oscillations when the MFE method is used with small time steps, we use the mass lumping procedure introduced in [74]. The lumped formulation of MFE was applied to density driven flow problems in [4]. In the following, we recall the main stages in order to use this formulation within a MOL context. With MFE, the velocity inside each triangle E is approximated with linear vectorial basis functions

q¼

3 X

Q Ej wEj ;

ð6Þ

j¼1

where Q Ej is the flux across the edge j of the element E and wEj are the Raviart-Thomas basis function given by

wEj ¼

1 2jEj



x  xE;i z  zE;i

 ;

ð7Þ

where |E| is the area of the element E, xE,i and zE,i are the coordinates of the vertex i of E (opposed to the edge i). The variational formulation of the Darcy’s law (2) using test function wEi leads to (see [4] for details):

Z

where the dispersion tensor D is given by

D ¼ Dm I þ ðaL  aT Þq  q=jqj þ aT jqjI

Flow and transport equations are coupled by state equations linking density and viscosity to mass fraction. We use a linear model for density and a power formulation for viscosity:

E

1

ðkE qÞ:wEi ¼





q0 g ðq  q 0 Þ E h  Thi þ E ðzE  zEi Þ ; lE E q0

ð8Þ

where hE is the average head in the element E, zE the z-coordinate of E the centre of E, Thi the average head on the edge i of E, zEi the z-coordinate of the midpoint of the edge i. Viscosity and density are approximated by constant values over E.

342

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

Combining (6) and (8) leads to the following matrix form: 3 X

qg BEij Q Ej ¼ 0 lE j¼1

    ðq  q0 Þ ðq  q0 Þ E E hE þ E zE  Thi þ E zi ;

q0

q0

xk ð9Þ

R T where the local matrix BEij ¼ E wEi ðkE Þ1 wEj can be expressed analytically for triangles (see [73]). With the mass lumping procedure, in a first step, we approxi E across the edge i of the element E under steady mate the flux Q i state conditions and without sink/source terms. By inverting (9) and using (1), we obtain (see [4] for details),

 E ¼ q0 g Q i

" X

lE

xj

#

E NEij Thj

j

q  q0 X E E þ E Nij zj : q0 j

ð10Þ

EÞ T With ¼  detðk r i k r j (see [73]), where ri is the edge vector face jEj to the vertex i of the element E. In the second step of the mass lumping procedure, sink/source and accumulation terms are distributed over edges as in the case of edge centred FV formulation. The total flux Q Ei across the edge i of E is given by

E jEjSE @Thi

 E  jEj eE @ q @C E  Q Ei ¼ Q i 3 qE @C @t 3

@t

ð11Þ

;

where eE, SE, CE are respectively the porosity, the storage coefficient and the concentration of the element E. Using (10), the general expression of the flux becomes,

"

Q Ei

#

q g X E E qE  q0 X E E jEj eE @ q @C E ¼ 0 N Th þ Nij zj  3 qE @C @t lE j ij j q0 j E



jEjSE @Thi : 3 @t

ð12Þ

The last step in the mass lumping procedure consists in writing con0 E E0 tinuities of pressure heads ðThi ¼ Thi Þ and fluxes ðQ Ei þ Q Ei ¼ 0Þ between elements E and E’ having a common edge i. The continuity is written including storage and sink/source terms distributed over each element edge. In order to use the MOL, terms with time derivative are assembled in the left hand side as following:

     0  E jEjSE jE0 jSE0 @Thi jEj eE @ q @C E jE j eE0 @ q @C E0 þ þ þ 3 qE @C @t 3 qE0 @C @t 3 3 @t " # X X qg q  q0 E ¼ 0 NEij Thj þ E NEij zEj

lE

þ

q0

j

" q0 g X

lE0

j

0

E

0

NEij Thj þ

j

qE  q0 X 0

q0

0

QD1 ,i Fig. 1. Triangle splitting into three subcells and linear concentration approximation on each subcell.

As stated in [6], the DG method for hyperbolic systems has been proven to be clearly superior to the already existing FE methods. In this work, the DG method, where fluxes are upwinded using a Riemann solver is used to solve the advection equation and combined with the multipoint flux approximation (MPFA) method for the dispersion equation. MPFA method has properties similar to the MFE method. Indeed, both are locally conservative and handle general irregular grids on anisotropic heterogeneous domains. The MPFA method is preferred to MFE since it can be combined with the DG method without the time splitting procedure [72]. Indeed, MPFA uses the same type of unknowns (average value per element) than DG and therefore, both discretizations can be gathered into one system matrix. Several papers have been devoted to MPFA methods recently [1,2,19,39] and to the link between MPFA and Raviart Thomas and Brezzi–Douglas–Marini MFE [40,46,65,70,71]. In the following, the main stages for the DG-MPFA spatial discretization are recalled in order to use this formulation within the MOL context. Substituting the mass conservation of the fluid in the transport equation allows to reduce the coupling between flow and transport equations without loss of accuracy [51]. The transport equation (3) can then be written in the following mixed form:

xj

# 0

NEij zEj :

QD2 ,i

xi

xi1

1

N Eij

xi2

x

ð13Þ

j

E5

This equation is written for each edge of the mesh and constitutes the first part of the final ODE system. The second part is obtained from the spatial discretization of the transport equation.

E1

3.2. Spatial discretization of the transport For the transport equation, it is known that when advection is dominant, standard numerical methods, such as FE or FV, generate solution with numerical diffusion and/or non-physical oscillations. These problems can be avoided with DG [60]. The DG method leads to a high-resolution scheme for advection that maintains the local conservation of FV methods but allows higher order approximations to enter through a variational formulation rather than by some hybridised difference or functional reconstruction [38]. The method has been used to diffusion-advection problems in [7,12,33] and several strategies have been used to adapt the DG method to elliptic problems [5,6]. A detailed review of DG methods can be found in [6,13,14].

E4

xi

xk

QD1 ,i

E2 E3

Fig. 2. The interaction region sharing the vertex i.

343

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

where oE is the boundary edge of the element E and C* the upstream mass fraction on oE. The expression of each term of (18) will be detailed in the following. P @C E R The first term in (18): j eE @tj E /Ej /Ei . This term corresponds to the mass integral, using (17), it leads to the following mass matrix:

Table 1 Parameters and boundary conditions for the modified Henry problem. kx = ky = 1.0204  109 m2 e = 0.35 aL = aT = 0 m Dm = 18.86  106 m2 s1

Permeability Porosity Dispersivity Molecular diffusion coefficient Boundary conditions for flow

Hydrostatic pressure at the right hand side Constant flux at the inland boundary: Q = 3.3  105 m2/s No flow along the top and bottom q0 = 1000 kg/m3 on the left boundary

Boundary conditions for transport

0 X

eE

@C Ej

j

Z

@t

E

@C E1 B @t

C B EC /Ej /Ei ¼ ½M 1 B @C 2 C: @ @t A

(

ð14Þ

qD ¼ D  rC:

The dispersive flux qD is assumed to vary linearly inside the element E, therefore,

Z

r:qD ¼

Z

E

qD  n@E ¼

@E

X

Q ED;i

)

r  qD ¼

i

1 X E Q D;i ; jEj i

ð15Þ

where Q ED;i is the dispersive flux across the edge i. We use the P1 DG method where the approximate solution Ch(x, t) is expressed with linear basis functions /Ei on each element E as follows:

C h ðx; tÞjE ¼

3 X

C Ei ðtÞ; /Ei ðxÞ;

ð16Þ

i¼1

where C Ei ðtÞ are the three unknown coefficients corresponding to the degrees of freedom. The three unknowns for each element are the average value of E Þ and its the mass fraction defined at the triangle centroid ð xE ; y deviations in each space direction [13] with the corresponding interpolation functions:

/E1 ðx; yÞ

¼ 1;

/E2 ðx; yÞ

¼ x  xE ;

/E3 ðx; yÞ

E : ¼yy

j

eE

@C Ej @t

Z E

/Ej /Ei 

XZ j

þ

E

Z

C @E

C Ej /Ej q:r/Ei  

/Ei q

 g@E þ

Z E

XZ j

E

C Ej /Ej /Ei r  q

1 X E E Q / ¼ 0; jEj @Ej D;j i

jEj 0 B ½M 1  ¼ eE @ 0 Ixx 0 Ixy

0

R

E

/Ej /Ei

1

C Ixy A

ð20Þ

Iyy

R R R  Þ2 and Ixy ¼ E ðx  xE Þðy  y E Þ. with Ixx ¼ E ðx  xE Þ2 ; Iyy ¼ E ðy  y P R E EE E The second term in (18): j E C j /j q:r/i . This term is written in a non-standard way, the Raviart-Thomas approximation of the flow velocity (6) is used to obtain a matrix form where unknowns are the water fluxes across edges

2

XZ E

j

Q E1

3

6 7 C Ej /Ej q:r/Ei ¼ ½M 2 4 Q E2 5;

ð21Þ

Q E3

where the terms of the matrix [M2] are given by M2 ði; jÞ ¼  R P3 E E E E k¼1 C k /k wj r/i . E P R The third term in (18): j E C Ej /Ej /Ei r  q. Again, this term is written using Eq. (6) and properties of the Raviart-Thomas vectorial basis functions (7) in the following matrix form:

2

XZ j

ð17Þ

The variational formulation of (14) over the element E using /Ei as test functions gives (see [72]),

X

The terms of the mass matrix are given by M 1 ði; jÞ ¼ eE which leads to

0

e @C þ q  rC þ r  qD ¼ 0; @t

ð19Þ

@C E3 @t

q1 = 1025 kg/m3 on the right boundary Zero concentration gradient along the top and bottom

1

E

Q E1

3

6 7 C Ej /Ej /Ei r:q ¼ ½M 3 4 Q E2 5;

ð22Þ

Q E3

where the terms of the matrix [M3] are given by M3 ði; jÞ ¼  R P3 E E E 1 k¼1 C k /k /i . jEj E R The fourth term in (18): @E C  /Ei q  g@E . This term corresponds to the boundary integral over the three edges of the element E

Z ð18Þ @E

C  /Ei q:g@E ¼

Z 3 X Q Ej C @Ej /Ei ; j@Ejj @Ej j¼1

Fig. 3. Comparison between semi-analytical (h) and numerical (–––) solutions for the modified Henry problem on a coarse mesh (400 triangles).

ð23Þ

344

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

Table 2 Number of time steps (Nb-dt), number of evaluations of the Jacobian matrix (Nb_Jac) and CPU time with the five codes for the Henry problem.

Nb_dt Nb_Jac CPU

Std_1

Std_2

Std_3

Mol_1

Mol_V

96 0 55.73

96 0 54.46

100 0 80.5

1301 77 40.2

291 42 15.8

the mass fraction values at respectively midpoint edges x1i and x2i  of the element E (Fig. 1). Therefore, the dispersive and the centre x  ; q1i Þ and flux qD ¼ D  rC is constant inside the subcell ðxi ; q2i ; q 1 2 sub-edge (half-edge) dispersive fluxes Q D;i and Q D;i (see Fig. 1), taken positive for outflow, are given by

Q 1D;i Q 2D;i

where goEj is the unit outward normal vector to the edge oEj of length |oEj| and Q Ej the water flux across oEj given by (12). C @Ej is the concentration over oEj, defined using the Riemann solver [64] which corresponds to the upstream concentration value:

  C @Ej ¼ kE@Ej C E@Ej þ 1  kE@Ej C Ej @Ej ;

ð24Þ

where Ej is the adjacent element to E such that oEj is the common edge of E and Ej. At each edge, we define,

( kE@Ej

¼

1 if

q:g@Ej  0;

0

q:g@Ej < 0:

if

ð25Þ

When written in a matrix form, (23) leads to

2

Z

C



/Ei q:g@E

@E

3 Q E1 6 E7 ¼ ½M4 4 Q 2 5

ð26Þ

Q E3

Due to the Riemann solver, the terms in the matrix [M4] depend on the mass fractions of the considered element E and those of its adjacent elements E1, E2, E3 following:

M 4 ði; jÞ ¼

kE@Ej j@Ejj

Z

C E@Ej /Ei þ

@Ej

R

1 E jEj

ð1  kE@Ej Þ j@Ejj

Z @Ej

E C Ej @Ej /i :

P

E E @Ej Q D;j /i .

The last term in (18): Using (17), this term can be written in the following form:

Z E

2P E 3 Q D;j X 6 @Ej 7 1 7 Q ED;j /Ei ¼ 6 4 0 5: jEj @Ej 0

ð27Þ

Dispersive fluxes Q ED;j across edges are approximated using the MPFA method. Contrarily to standard FV schemes, the MPFA method allows to treat rigorously full discontinuous dispersion tensor on unstructured meshes. The basic idea of the MPFA method [1] is to divide each triangle ; x1i Þ, we asinto three subcells (Fig. 1). Inside the subcell ðxi ; x2i ; x sume a linear variation of the mass fraction between c1i ; c2i and C E1

Table 3 Parameters, initial and boundary conditions for the Elder problem. Permeability Porosity Dispersivity Molecular diffusion coefficient Dynamic viscosity State equation Initial conditions Boundary conditions

kx = kz = 4.845 1013 m2 e = 0.1 aL = aT = 0 m Dm = 3.565106 m2/s

l = 103kg/m/s q = q0+200 C p(x, z) = 0, (x, z) e [0, 300]  [0, 150] C(x, z) = 0, (x, z) e [0, 300]  [0, 150] pðx; zÞ ¼ 0; x 2 ½0; 1; z ¼ 150 nablap ¼ 0; elsewhere Cðx; zÞ ¼ 1; x 2 ½150; 300; z ¼ 150 Zero concentration Cðx; zÞ ¼ 1; x 2 ½0; 300; z ¼ 0 gradient elsewhere

!

! Þ? ðx1i  xi Þ? Dðx   x1i Þ? ðx1i  xi Þ? Dðx2i  x ¼ A 2jT xx1 x2 j ðxi  x2i Þ? Dðx2i  x Þ? ðxi  x2i Þ? Dðx   x1i Þ? Gi i i |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ! c1i  C E1 ; ð28Þ  c2i  C E1 1

where jT xx1 x2 j is equal to the area of the triangle spanned by the i i , x1i and x2i and the vector ðx1i  xi Þ? is obtained by a p/2 points x rotation of the vector x1i  xi . All subcells sharing the vertex i create an interaction volume (see Fig. 2). The discretization is accomplished by assuming continuous dispersive fluxes across each of the subedges of the interaction region, and a weak continuity condition of the mass fraction across the same subedges. From these assumptions, an explicit discrete flux can be found after the resolution of a local linear system and eliminating the edge mass fraction for each subedge of the interaction volume. Each subedge dispersive flux can then be written explicitly as a weighted sum of the cell concentrations of elements forming the interaction volume. For example, for the Fig. 2 we obtain

Q 1D;i ¼

5 X

E

t ki C 1k ;

ð29Þ

k¼1

where t ki are transmissibility coefficients. Eq. (29) is written for the six half fluxes of the element E. The sum of these fluxes gives the exchange of mass by dispersion between E and its adjacent elements. The sum of dispersive fluxes is then replaced in (27). Note that contrary to advection terms (the second, third and fourth terms of (18)), the dispersion term (the last term of (18)) cannot be expressed as a linear function of the water fluxes across edges. Indeed, the dispersion tensor D used in (28) is nonlinear dependent on the velocity vector (see (4)). Moreover, the expression of the transmissibility coefficients tki on an unstructured triangular mesh cannot be obtained analytically since it requires the inversion of a local system which depends on the number of elements sharing the node i. For theses reasons and because usually in density driven flow problems, the transfer of solute is mainly driven by advection, the transmissibility coefficients tki are computed using dispersion tensor performed with the flow velocity at the old time level. For example, for Fig. 2, the dispersive flux at the new time level is written as,

ðQ 1D;i Þnþ1 ¼

5 X E ðtki Þn ðC 1k Þnþ1 ;

ð30Þ

k¼1

where n is the time level. 3.3. Forming the global ODE system In this section, we write the flow and transport equations in a ¼ Gðt; yÞ where A(t, y) is a single linearly implicit form Aðt; yÞ dy dt square matrix and y the vector of unknowns which is formed by E both flow unknowns (the pressure head Thi at each edge) and transport unknowns (the three mass fractions ðC E1 ; C E2 ; C E3 Þ for each element).To this aim, we proceed as following:

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

Fig. 4. Simulation results of the Elder test case at 20 years with mesh Levels 4 and 5.

345

346

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

1. the transport equation (18) is written in the following form:

0

E 1

@C 1 2 E 3 2 P QE 3 Q1 D;j B @tE C 7 B C 6 7 6 @Ej 7 ½M 1 B @C 2 C ¼ ½M4 Q E2 5 þ 6 4 0 5; @ @t A E @C E3 Q3 0

ð31Þ

@t

where the matrix [M] is given by [M] = [M2] + [M3] + [M4]. 2. The flow fluxes in (12) are written in the following matrix form:

2

Q E1

2

3

E

Th1

3

2

zE

3

1 7 qE  q0 6 E 7 q0 g E 6 6 7 E þ ½N 6 g½NE 4 zE2 5 4 Q2 5 ¼ Th2 7 4 5 lE lE E E E

Q3

Th3 0



@C E1 B @t

1

0

@ThE1 B @t

z3 1

C jEj eE @ q B @C E C C jEjSE B @ThE2 C B @t1 C  B @t C: A 3 qE @C @ A 3 @ @C E1

@ThE3

@t

@t

ð32Þ

3. The flux expressions (32) are substituted into (31), which yields

0

@ThE1 B @t

1

0

@C E1 B @t

1

0

@C E1 B @t

1

C C C jEjSE B EC B E C jEj eE @ q B EC ½MB @Th2 C þ ½M1 B @C 2 C þ ½MB @C1 C @t @t @t A @ @ A 3 qE @C @ A 3 @ThE3 @t

2

@C E3 @t

@C E1 @t

3 2 E 3 2 P QE 3 E Th1 D;j z1 6 7 6 @Ej 7 qg qE  q0 E 6 E7 E7 6 7: ¼ 0 ½M½N E 6 þ þ g½M½N  z 4 5 Th 2 4 5 4 5 2 0 lE lE E E z Th3 3 0 ð33Þ Therefore the final matrix A is obtained by written the flow equation (13) for each edge of the mesh, and then, for each element, the transport equation (33). The final matrix A is a square matrix with a rank equal to the number of edges plus three times the number of elements of the mesh. 5. Temporal discretization Higher-order methods in time are an advantage if they advance the approximate solution in time for less computational effort than lower-order methods. This situation can occur if higher-order methods lead to longer time steps and/or less effort in the linear or nonlinear solver compared to the lowest order methods [62]. In order to investigate the advantage of using higher-order methods in time, one can use mature and sophisticated time integration packages which are available for large-scale systems of ODE or DAEs. In this work, we implement the recent publicly available version of DLSODIS developed 18 November 2003 which is briefly described in this section. DLSODIS is a set of general-purpose Fortran routines solver for the ODE systems given in the linearly implicit form A(t, y)dy/ dt = G(t, y). The matrix A is a square matrix and can be singular, in which case the system is a DAE system. For the coupled flow transport system formed by equations (13) and (33), the matrix A is singular if the storage coefficient is zero. The solver is a double precision version of LSODI (Livermore Solver for Ordinary Differential equations-implicit form [27–30,54]) where the matrices involved (mass and Jacobian matrices) are all assumed to be sparse. DLSODIS uses a fixed-coefficient implementation of the implicit Adams and BDF methods. The maximum order allowed is 12 for Adams and 5 for BDF. The BDF approach is chosen during this wok because it works well on stiff problems and has good

stability properties [36]. With BDF, the solver uses a modified Newton iteration (quasi-Newton method). Thus, the value of the Jacobian matrix and therefore the iteration matrix are approximated from a previous iterate and/or time level whenever possible. Reusing the Jacobian matrix can result in significant computational savings but reduces the formal convergence rate of the nonlinear solver from quadratic to superlinear [37]. Moreover, matrix iteration may be altered if the step size and/ or method order is changed. In this case, the code updates the matrix and performs its LU-decomposition in order to minimize convergence failures caused by the inaccurate matrix iteration. Note that in the DLSODIS code, the same iteration matrix is used for a maximum number of 20 steps, after which it is revaluated [57]. With DSLOSIS, the Jacobian matrix can be provide analytically by the user or calculated numerically by the solver using finite difference approximation. In this case, the number of evaluations of the residue is optimized using the column grouping technique of Curtis et al., [15,28] and the variables are perturbed by group. DLSODIS determines the sparsity structure on its own or optionally accepts this information from the user to enable efficient solution. It uses the CDRV solver from of the Yale Sparse Matrix Package (YSMP) to solve the linear systems that arise [20,21]. The CDRV direct solver for nonsymmetric systems performs the following four steps: (i) reordering, (ii) symbolic factorization, (iii) numerical factorization and (iv) linear system solution. To obtain significant computational savings, the two first steps are performed only once (for the first iteration of the simulation). Moreover, when the matrix is unchanged, only the last step (iv) is performed. Note that periodically, the DSLODIS solver attempts to change the step size and/or the method order to minimize computational work while maintaining prescribed accuracy. The local accuracy is specified to the code by both a relative er and absolute ea local error tolerance of the form si 6 ei = er|yi|+ea where si is the local temporal truncation error and yi is component of the solution vector. Finally, to initialize the integration, the solver requires the value 0 of y (0) as well as y(0). If the matrix A is non singular, the user can 0 provide only y(0) and the solver can find a consistent value of y (0). 0 Whereas, both values of y (0) and y(0) should be provided by the user if the matrix A is singular. For the studied coupled flow-transport system formed by equations (13) and (33), only y(0) is known 0 and the value of y (0) is computed by the solver. To this aim, the matrix A is made non singular by setting the storage coefficient to 1012 if it should be zero. 6. Numerical experiments In this part, we investigate the benefit of using the MOL for density driven flow problems. Numerical experiments are performed for the simulation of two standard numerical test cases and for a laboratory scale experiment performed with a high viscosity contrast between injected and displaced fluids.

Table 4 Number of time steps (Nb-dt), number of evaluations of the Jacobian matrix (Nb_Jac) and CPU time with the five codes for the Elder problem on the mesh level 5.

Nb_dt Nb_Jac CPU

Std_1

Std_2

Std_3

Mol_1

Mol_V

593 0 1135

590 0 1058

668 0 752

1626 93 155

794 73 94

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

Fig. 5. Results of the Elder problem with different initial time steps.

347

348

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

Five numerical codes are used for the simulations:  Std_1: the standard iterative method based upon the Picard linearization where (i) the flow system (13) is solved to obtain fluxes across edges, (ii) these values are then used for both advection and dispersion to solve the transport system (31) in order to obtain the three unknown mass fractions at each element, (iii) the Riemann problem (25) is solved at each iteration and (iv) densities and viscosities are updated before solving the flow at the new iteration. The temporal approximation is based on the first order backward Euler for both flow and transport. The global time step management during the simulation is of heuristic type. The time step Dtn is changed depending on the number of iterations k necessary to reach the convergence in the following way:

8 > if k < 5 Dt nþ1 ¼ 1:05Dt n ; > < if 5 k 10 Dt nþ1 ¼ 1:0Dt n ; > > : if k > 10 Dt nþ1 ¼ 0:8Dt n :











ð34Þ

The convergence is reached if the maximum change on the pressure head and concentration is less than a fixed tolerance. If the convergence is not reached after 20 iterations, the computation is restarted with a smaller time step (the time step is divided by 2). Std_2: This code is similar to Std_1. However, in order to test the validity of the modifications made previously for the implementation of the MOL, the Riemann problem (25) and the dispersion tensor (4) are not updated at each iteration but only at each time step. Std_3: This code is similar to Std_1 but imitate the well known IMPES (IMplicit Pressure-Explicit Saturation) scheme. Therefore, the flow is solved implicitly whereas, the transport is solved by combining an explicit DG discretization for advection which now requires a slope limiting procedure and an implicit MPFA discretization for dispersion. We use the slope limiter given in [31] for two-dimensional DG method on triangular elements. MOL_1: The global system A(t, y)dy/dt = G(t, y) formed by both flow equations (13) and transport equations (33) is solved with variable order and variable step-size discretization using the DLSODIS solver. MOL_1: In order to show the effect of the temporal discretization, we use the DLSODIS solver with the order of the integration method fixed at one, the code in this configuration only adapts the time step size. The 5 previous codes are used to simulate two standard test cases, the Henry saltwater intrusion problem and the Elder salt convection problem as well as a laboratory scale experiment.

Remarks 1 1. Simulations with the five codes with varying tolerance values between 102 and 1010 show that the CPU cost highly

increases for very small values (<106) and accuracy can deteriorate for large values (>103). A compromise is obtained where the error tolerance is fixed to 105 for standard codes. The relative er and absolute ea local errors are also fixed to 105 for the DSLODIS solver. 2. The direct CDRV solver used in DSLODIS is also implemented in the three first codes to solve separately flow and transport systems. All codes are optimized to take advantage of the direct solver. For example, the factorization step is performed only when required (when the matrix has been changed). 3. The code Std_1 can also be used with the Jacobi linearization method. In this case flow and transport are performed using values of the previous iteration. This scheme is less efficient (it requires very small time steps) than with the Picard linearization where the transport uses the fluxes obtained at the new iteration. 4. For all studied problems, DLSODIS is found to be less efficient when the Jacobian matrix is provided analytically than when only its structure is provided and the Jacobian matrix is calculated numerically by the solver using finite difference approximation. Indeed, the number of evaluations of the residue is kept to a minimum by a column grouping technique where the variables are perturbed by group (see [15] for more details). Therefore, for all simulations, the Jacobian is evaluated numerically. 6.1. Simulation of the modified Henry problem The Henry problem is one of the most widely used test cases for benchmarking density driven flow models. The problem deals with salt water intrusion which remains an important issue. It describes the advance of a saltwater front in a homogeneous confined aquifer initially free of salt. It is a classic test because of the existence of a semi-analytic steady-state solution developed by Henry [26] and revaluated by Segol [58]. The worthiness of the Henry problem has been improved in Simpson and Clement [61] by halving the freshwater recharge to increase density effects as compared to the boundary forcing. Results of simulations of this modified problem are compared to the semi-analytic results obtained by Simpson and Clement [61] using 203-term truncation for the semi-analytical solution. All simulations are conducted for 280 min for which the position of the 25%, 50% and 75% isochlors become steady. The parameters and boundary conditions are given in Table 1. The simulations are performed on a very coarse mesh of only 400 triangles and on a finer mesh of 1600 triangles used in [61]. Similar results are obtained with both meshes. The five codes give the same results which show that the accuracy is not altered when

Fig. 6. Results of Elder problem with fine meshes using the Mol_V code.

349

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

the Riemann problem (25) and the dispersion tensor (4) are not updated at each iteration but only at each time step. The coarse mesh results of the MOL_V code are plotted in Fig. 3. Even with the coarse mesh, the numerical results are very close to the semi-analytical solution given by Simpson and Clement [61]. Efficiency of the different codes for the fine mesh simulation is shown in Table 2. The Std_1 and Std_2 codes use the same number of time steps and require similar CPU times. Both require an average of four iterations per time step to reach the convergence criteria. Due to the CFL constraint, the explicit scheme Std_3 requires more CPU time. MOL_1 allows a reduction of the CPU time of about 27% compared to Std_1 and Std_2 codes. Note that the MOL_1 code is based upon the quasi Newton method whereas the three first codes are based upon the fixed point (Picard) scheme. MOL_1 uses 1301 time steps and requires only 77 evaluations of the Jacobian matrix (Table 2). The MOL_V code where both the order of approximation and the time step size are varying during the simulation is the more efficient code. It allows larger time steps. Indeed, the simulation is performed using 291 time steps and only 42 evaluations of the Jacobian matrix (Table 2). The CPU time with MOL_V is 3.5 times less than with Std_1 and Std_2 codes. 6.2. Simulation of the elder problem The Elder problem is a free convection problem where fluid flow is driven purely by fluid density differences in a 2D domain of 600 m  150 m. The problem involves total density variations of 20% which makes this a strongly coupled flow case. The viscosity is assumed to be constant. The problem being symmetrical at x = 300 m, only one half of the domain is simulated. The parameters, initial and boundary conditions for the problem are given in Table 3. A triangular mesh, obtained by dividing square elements in 2, is used. The number of squares is obtained following [25], i.e. 2  4l quadrilaterals for the half domain. The orientation of the diagonals is optimized to avoid a non-symmetric mesh and possible effects on the numerical results [9].

Volumetric pump y

Input tank

x z

Z1 Z2

1

2

3

4

Z3 Z4 Z5 Porous filter

Fig. 7. Experimental setup.

Table 5 Porous media characteristics for Loggia experiments (H: thickness, d: average diameter, k: permeability). Layer 1 2 3 4

Hi (m)

ki (m2)

di (m) 2

1.22  10 1.08  102 1.08  102 1.22  102

6

85  10 155  106 194  106 238  106

7.1  1012 23.7  1012 37.2  1012 55.8  1012

Numerous published results [3,9,25,41,50,66] discussed the existence of an upwelling or downwelling flow in the central part of the domain at 20 years simulation time. Generally, downwelling flow has been found with coarse meshes [17,25,66] and upwelling flow with fine meshes [4,9,41,50,51]. 6.2.1. Results for coarse meshes The five codes are used to simulate the Elder problem on the mesh levels 4 and 5 (l = 4 and 5 respectively). The three first codes start with a small time step (fixed to 1 day) which increases during the simulation following the procedure (34). The results of the simulations at t = 7300 days are plotted in Fig. 4. All simulations give a centred upwelling. As previously, Std_1 and Std_2 give similar results which confirm that the accuracy is not altered when the Riemann problem (25) and the dispersion tensor (4) are updated at each time step instead of at each iteration. Efficiency of the different codes for the simulation with the mesh level 5 is shown in Table 4. The Std_1 code requires more CPU time than Std_2 since the Riemann problem is solved at each iteration. Both require an average of five iterations per time step to reach the convergence criteria. The explicit scheme Std_3 gives similar results as the implicit schemes (Fig. 4) but requires less CPU time. The Mol_1 code which uses the quasi Newton method allows a strong reduction (between 5 and 6) of the CPU time compared to the Std_1 code. MOL_1 uses 1626 time steps and requires only 93 evaluations of the Jacobian matrix (Table 4). The MOL_V code where both the order of approximation and the time step size are varying during the simulation is the more efficient code. The gain in CPU time is more important than for the Henry problem since the Elder problem is highly nonlinear. With MOL_V, the simulation is performed using only 494 time steps and 73 evaluations of the Jacobian matrix. The MOL_V code requires 10 times less CPU time than the Std_1 code. 6.2.2. Temporal error With the MOL, the time step size is adapted using a detailed algorithm which ensures that the local truncation error is met. On the other hand, the time step control with standard codes is based upon the effort required for the convergence as given by (34). This procedure gives a rough measure of the actual temporal error which may lead to unreliable and/or inaccurate results. Fig. 5 shows results obtained with the Std_1 and Std_3 codes s when starting with a small time step of Dt init ¼ 1 day or with a large time step of Dt init ¼ 20 day for mesh levels 4 and 5. This figure shows that depending on the time step size used at the beginning of the simulation, we can obtain different results with an upwelling or a downwelling flow in the central part of the domain at 20 years simulation time. This difference exists even if the tolerance value (fixed to 105) is decreased to 1012. This shows that the temporal error with standard codes may remain significant even if the convergence is reached. This problem can also be observed with finer meshes. Fig. 5 shows results of Levels 6 and 7 with the Std_3 code when starting with a small time step of Dtinit ¼ 1 day and with a very small time step of Dt init ¼ 0:01 day. Different results are obtained with a central downwelling if Dtinit ¼ 1 day and a central upwelling if Dt init ¼ 0:01 day. This

350

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

shows that for fine meshes, the temporal error may remain significant with small time steps (even if the convergence is reached) and very small time steps have to be used to reduce this error. This problem is avoided with the MOL since the method adapt in both the order of approximation and time step to provide the necessary accuracy, efficiency and robustness.

0.8

C/C0

6.2.3. Efficiency of the MOL on fine meshes To show the efficiency of the MOL, results of Levels 6 and 7 with the Mol_V code are given in Fig. 6. These results show an upwelling in the central part of the domain at 20 years simulation time and are similar to the results obtained with the Std_3 code when used with very small initial time steps (see Fig. 5). The CPU time required by the Mol_V code is about 30 times less than with the Std_3 code for the very fine mesh Level 7. This shows the high efficiency of the Mol_V code which uses variable-order, variable step-size time discretization to maximize the size of the time step while maintaining a specified temporal integration error.

1.0

Z= 4.5 Z= 9.7 Z=14.9 Z=20.1 Z=25.3

0.6 0.4 0.2 0.0 0

10000 Time (s)

5000

15000

20000

Fig. 9. Measured (symbol) and computed breakthrough curves (exp. 2).

6.3. Simulation of Loggia’s experiments In this part, we use the laboratory experiments run by Loggia [44] as a benchmark to study stable flow occurring in a layered medium with density and viscosity contrasts. Mixture of water-sucrose and water-glycerin are used to obtain large variations of viscosity where the injected fluid can be seven times more viscous than the displaced one. This increases the nonlinearity of the problem since the viscosity term is given in the denominator of expressions (13) and (33). The laboratory experiments were run on a tank of x = 4.5 cm, y = 4.5 cm and z = 30 cm (Fig. 7), filled with glass beads of various sizes distributed in layers. The characteristics of each layer are given in Table 5. The longitudinal dispersivity is equal to the diameter of the beads and the transverse dispersivity is assumed to be 10 times smaller. The same porosity (e = 0.4) is used for the four layers. An input tank ensures homogeneous pressure and concentration at the upstream boundary of the layers. All experiments are done at a fixed temperature of 30 °C. A volumetric pump maintains a fixed flow rate at the entry of the domain. Measured concentrations are based on acoustic process. It gives the average of concentration at five fixed altitudes (Fig. 7). Two experiments were simulated, one with the injection of a tracer, the second with a fluid of high density (q1  q0 = 53 kg/ m3) and viscosity (l1/l0 = 7.54).

1.0

C/C0

0.8 0.6

Z= 4.5 Z= 9.7 Z=14.9 Z=20.1 Z=25.3

0.4

Due to symmetry, the system is modeled in a 2D cross section perpendicular to the layers (X–Z section). The measured permeability is k ¼ 5:1012 m2 and porosity is equal to 0.2 for the porous plate. Longitudinal and transversal dispersivities are assumed to be equal to 85  106 m and 8.5  106 m. The triangular discretization is obtained by subdivision of a 24  142 rectangular mesh. The first experiment is a tracer test, where density and viscosity of the injected fluid and the displaced one are very close. The numerical results obtained with the five codes are similar. Although no model calibration has been performed, the match between simulated and measured concentrations is very good (Fig. 8). In the second experiment, the injected fluid has a viscosity which is about 7.5 times higher than the displaced fluid. All codes give similar results. The match between simulated and measured mass fractions is satisfactory (Fig. 9). The effects of the heterogeneities are smoothed out by the viscosity contrast. The CPU time required by the five codes is given in Table 6. Due to the small opening at the output, the explicit Std_3 code requires very small time steps in order to satisfy the CFL constraint ðCFL 1=2Þ. On the other hand, the implicit schemes (Std_1 and Std_2) are not subject to the stability constraint. Results of Table 6 show that Mol_V is highly efficient, it requires about 15 times less CPU time than the Std_1 code. The results obtained for the three test cases point out the high efficiency of the developed scheme for 2D density driven flow problems. Note that for 3D problems, direct methods (as used in the DLSODIS solver) can be very restrictive (because of memory requirement) and iterative methods should be used. In this case, the gain can be less important than with direct solvers. Indeed, with direct methods the iteration matrix is computed and factored, and is then used over as many time steps as possible. However, higher-order time integration schemes with formal error estimation and control and adaptative time stepping combined with iterative methods can still be highly efficient compared to standard schemes as shown in [36] where the preconditioned Krylov method is used with higher-order integration in time for two-phase flow problems.

0.2 0.0 0

5000

10000 Time (s)

15000

20000

Fig. 8. Measured (symbol) and computed breakthrough curves (exp. 1).

Table 6 CPU time with the different codes for the Loggia experiment with high viscosity contrast.

l1 l0 ’ 7:54

Std_1

Std_2

Std_3

Mol_1

Mol_V

5058

5099

8026

856

349

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

7. Conclusion Due to nonlinearities, modeling density driven flow problems requires a coupled flow-transport numerical model which requires a lot of computational time and/or heavy equipments. We developed in this work a robust numerical model with efficient advanced approximations for both spatial and temporal discretizations in order to reduce the excessive computational requirement while maintaining the accuracy. The model is developed for a general triangular mesh and uses a spatial discretization based upon the lumped formulation of MFE for the flow equation, the DG method for advection and the MPFA method for dispersion. For the temporal discretization, the model uses higher-order time integration methods with formal error estimation and control and adaptative time stepping. The developed model is used for the simulation of the Henry and the Elder test cases as well as a laboratory scale experiment with a high viscosity contrast between injected and displaced fluids. Numerical results show that standard numerical models where the time step management is based upon the number of iterations required to reach the convergence and does not explicitly measure temporal truncation errors of the solution may lead to inaccurate results. This problem is avoided with the developed model since the scheme adapt in both the order of approximation and time step to provide the necessary accuracy, efficiency and robustness. The numerical results highlights the efficiency of the developed scheme for density driven flow problems since it can be more than 20 times more efficient than the standard methods. Acknowledgements This work was partly supported by the GdR MoMas CNRS-2439 sponsored by ANDRA, BRGM, CEA and EDF whose support is gratefully acknowledged. References [1] Aavatsmark I. An introduction to multipoint flux approximations for quadrilateral grids. Comput Geosci 2002;6:404–32. [2] Aavatsmark I, Barkve T, Bøe Ø, Mannseth T. Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media. J Comput Phys 1996;127:2–14. [3] Ackerer Ph, Younes A, Mosé R. Modeling variable density flow and solute transport in porous medium: 1. Numerical model and verification. Transport Porous Med 1999;35(3):345–73. [4] Ackerer P, Younes A. Efficient approximations for the simulation of density driven flow in porous media. Adv Water Res 2008;31:15–27. [5] Aizinger V, Dawson C, Cockburn B, Castillo P. The local discontinuous Galerkin method for contaminant transport. Adv Water Res 2001;24:73–87. [6] Arnold DN, Brezzi F, Cockburn B, Marini LD. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal. 2002;5:1749–79. [7] Baumann CE, Oden JT. A discontinuous hp finite element method for convection–diffusion problems. Comput Methods Appl Mech Eng 1999;175:311–41. [8] Beinhorn M, Dietrich P, Kolditz O. 3D numerical evaluation of density effects on tracer tests. J Contam Hydrol 2005;81:89–105. [9] Boufadel MC, Suidan MT, Venosa AD. Numerical modelling of water flow below dry salt lakes: effect of capillarity and viscosity. J Hydrol 1999:55–74. [10] Brezzi F, Fortin M. Mixed and hybrid finite element methods. New York: Springer; 1991. [11] Chavent G, Roberts JE. A unified physical presentation of mixed, mixed-hybrid finite elements and standard finite difference approximations for the determination of velocities in waterflow problems. Adv Water Res 1991;14:329–48. [12] Cockburn B, Hou S, Shu CW. TVB Runge Kutta local projection discontinuous Galerkin finite element method for conservative laws: III. One dimensional systems. J Comput Phys 1989;84:90–113. [13] Cockburn B, Shu CW. The Runge-Kutta discontinuous Galerkin method for conservative laws: V. Multidimentional systems. J Comput Phys 1998;141: 199–224. [14] Cockburn B, Karniadakis GE, Shu CW, editors. Discontinuous Galerkin methods: theory, computation and applications. Lecture notes in computational science and engineering, vol. 11. Springer-Verlag; 2000.

351

[15] Curtis AR, Powell MJD, Reid JK. On the estimation of sparse Jacobian matrices. J Inst Math Appl 1974;13:117–9. [16] Diersch HJ, Kolditz O. Coupled groundwater flow and transport: 2. Thermoline and 3D convection systems. Adv Water Res 1998;21:401–25. [17] Diersch HJ, Kolditz O. Variable-density flow and transport in porous media: approaches and challenges. Adv Water Res 2002;25:899–944. [18] Durlofsky L. Accuracy of mixed and control volume finite element approximations to Darcy velocity and related quantities. Water Resour Res 1994;30:965–73. [19] Edwards MG, Rogers CF. Finite volume discretization with imposed flux continuity for the general tensor pressure equation. Comput Geosci 1998;2: 259–90. [20] Eisenstat SC, Gursky MC, Schultz MH, Sherman AH. Yale Sparse Matrix Package: II. The Nonsymmetric Codes. Research Report No. 114, Department of Computer Sciences, Yale University, 1977. [21] Eisenstat SC, Gursky MC, Schultz MH, Sherman AH. Yale sparse matrix package: I. The symmetric codes. Int J Numer Methods Eng 1982;18:1145–51. [22] Farthing MW, Kees CE, Miller CT. Mixed finite element methods and higherorder temporal approximations. Adv Water Res 2002;25:85–101. [23] Farthing MW, Kees CE, Miller CT. Mixed finite element methods and higherorder temporal approximations for variably saturated groundwater flow. Adv Water Res 2003;26:373–94. [24] Farthing MW, Kees CE, Coffey TS, Kelley CT, Miller CT. Efficient steady-state solution techniques for variably saturated groundwater flow. Adv Water Res 2003;26:833–49. [25] Frolkovic P, De Schepper H. Numerical modelling of convection dominated transport coupled with density driven flow in porous media. Adv Water Res 2001;24:63–72. [26] Henry HR. Interfaces between salt water and fresh water in coastal aquifers, US Geological Survey Water-Supply Paper 1613-C, Sea Water in Coastal Aquifers, 1964. p. C35–70. [27] Hindmarsh AC. LSODE and LSODI, two new initial value ordinary differential equation solvers. SIGNUM Newslett 1980;15:10–1. [28] Hindmarsh AC. Large ordinary differential equation systems and software. IEEE Cont Sys Mag 1982;2:24–30. [29] Hindmarsh AC. ODEPACK*: a systematized collection of ODE solvers. In: Stepleman RS et al., editors. Scientific computing. Amsterdam, The Netherlands: North-Holland; 1983. p. 55–64. [30] Hindmarsh AC. The ODEPACK solvers. In: Aiken RC, editor. Stiff computation. Oxford University Press; 1985. p. 167–74. [31] Hoteit H, Ackerer Ph, Mosé R, Erhel J, Philippe B. New two-dimensional slope limiters for discontinuous Galerkin methods on arbitrary meshes. Int J Numer Methods Eng 2004;61:2566–93. [32] Huber R, Helmig R. Multiphase flow in heterogeneous porous media: a classical finite element method versus an implicit pressure-explicit saturationbased mixed finite element-finite volume approach. Int J Numer Methods Fluids 1999;29:899–920. [33] Hughes TJR, Masud A, Wan J. A stabilized mixed discontinuous Galerkin method for Darcy flow. Comput Methods Appl Mech Eng 2006;195(25-28):3347–81. [34] Huyakorn P, Anderson P, Mercer J, White H. Saltwater intrusion in aquifers: development and testing of a three-dimensional finite element model. Water Resour Res 1987;23:293–312. [35] Johannsen K, Kinzelbach W, Oswald S, Wittum G. The salt-pool benchmark problem-numerical simulation of salt water upcoming in a porous medium. Adv Water Res 2002;25:335–48. [36] Kees CE, Miller CT. Higher order time integration methods for two-phase flow. Adv Water Res 2002;25:159–77. [37] Kelly CT. Iterative methods for linear and nonlinear equations. Philadelphia: SIAM; 1995. [38] Kirby R. A posteriori error estimates and local time-stepping for flow and transport problems in porous media. PhD thesis, University of Texas at Austin, 2000. [39] Klausen RA, Russell TF. Relationships among some locally conservative discretization methods which handle discontinuous coefficients. J Comput Geosci 2004;8(4):1–37. [40] Klausen RA, Winther R. Convergence of multipoint flux approximations on quadrilateral grids. Numer Methods Partial Differ Equat 2006;22(6):1438–54. [41] Kolditz O, Ratke R, Diersch HJ, Zielke W. Coupled groundwater flow and transport: 1. Verification of variable-density flow and transport models. Adv Water Res 1998;21:27–46. [42] Lang S, Wittum G. Large-scale density-driven flow simulations using parallel unstructured grid adaptation and local multigrid methods. Concurr Comput: Pract Exp 2005;17:1415–40. [43] Li H, Farthing MW, Dawson CN, Miller CT. Local discontinuous Galerkin approximations to Richards’ equation. Adv Water Res 2007;30:555–75. [44] Loggia D, Rakotomalala N, Salin D, Yortsos YC. Phase diagram of stable miscible displacements in layered porous media. Europhys Lett 1996;36:105–11. [45] Ludwig R, Schelkes K, Vogel P, Wollrath J. Implications of large-scale heterogeneities for hydraulic model studies at the potential site of a radioactive waste repository at Gorleben, Germany. Eng Geol 2001;61(2-3):119–30. [46] Mayur Pal, Edwards MG, Lamb AR. Convergence study of a family of fluxcontinuous, finite-volume schemes for the general tensor pressure equation. Int J Numer Methods Fluids 2006;51:1177–203. [47] Miller CT, Williams GA, Kelly CT, Tocci MD. Robust solution of Richards’equation for nonuniform porous media. Water Resour Res 1998;34: 2599–610.

352

A. Younes et al. / Advances in Water Resources 32 (2009) 340–352

[48] Miller CT, Abhishek C, Farthing M. A spatially and temporally adaptive solution of Richards’ equation. Adv Water Res 2006;29:525–45. [49] Mosé R, Siegel P, Ackerer P. Application of the mixed hybrid finite element approximation in a groundwater flow model: luxury or necessity? Water Resour Res 1994;30:3001–12. [50] Oldenburg C, Pruess K. Dispersive transport dynamics in a strongly coupled groundwater-brine flow system. Water Resour Res 1995;31:289–302. [51] Oltean C, Buès MA. Coupled groundwater flow and transport in porous media. A conservative or non-conservative form? Transport Porous Med 2001;44(2): 219–46. [52] Oude Essink GHP. Improving fresh groundwater supply-problems and solutions. Ocean Coastal Manage 2001;44(5-6):429–49. [53] Oude Essink GHP. Salt water intrusion in a three-dimensional groundwater system in The Netherlands: a numerical study. Transport Porous Med 2001;43:137–58. [54] Painter JF. Solving the Navier–Stokes equations with LSODI and the method of lines. LLNL Report UCID-19262, 1981. [55] Paniconi C, Khlaifi I, Lecca G, Giacomelli A, Tarhouni J. A modelling study of seawater intrusion in the korba coastal plain, Tunisia. Phys Chem Earth (B) 2001;26(4):345–51. [56] Putti M, Paniconi C. Picard and Newton linearization for the coupled model of saltwater intrusion in aquifers. Adv Water Resour 1995;18(3):159–70. [57] Radhakrishnan K, Hindmarsh AC. Description and Use of LSODE, the Livermore solver for ordinary differential equations. LLNL report UCRL-ID-113855, December 1993. [58] Segol G. Classic groundwater simulations: proving and improving numerical codes. Englewood Cliffs, NJ: Prentice-Hall; 1994. [59] Shuangzhang T, Shahrouz A. A slope limiting procedure in Discontinuous Galerkin finite element method for gas dynamics applications. Int J Numer Anal Model 2005;2:163–78. [60] Siegel P, Mosé R, Ackerer P, Jaffre J. Solution of the advection-diffusion equation using a combination of discontinuous and mixed finite elements. Int J Numer Methods Fluids 1997;24(6):595–613. [61] Simpson MJ, Clement TP. Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models. Water Resour Res 2004;40:W01504. doi:10.1029/2003WR00219.

[62] Tocci MD, Kelly CT, Miller CT. Accurate and economical solution of the pressure-head form of Richards’ equation by the method of lines. Adv Water Res 1997;20:1–14. [63] Tocci MD, Kelly CT, Miller CT, Kees CE. Inexact Newton methods and the method of lines for solving Richards’ equation in two space dimensions. Comput Geosci 1998;2:291–309. [64] Riemann Toro E. Solvers and numerical methods for fluid dynamics. Berlin: Springer; 1997. ´ [65] Vohralik M. Equivalence between lowest-order mixed finite element and multi-point finite volume methods on simplicial meshes. Math Model Numer Anal 2006;40(2):367–91. [66] Voss C, Souza W. Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater–saltwater transition zone. Water Resour Res 1987;26:2097–106. [67] Wheeler MF, Yotov I. A multipoint flux mixed finite element method. SIAM 2006;44(5):2082–106. [68] Woumeni R, Vauclin M. A field study of the coupled effects of aquifer stratification, fluid density, and groundwater fluctuations on dispersivity assessments. Adv Water Res 2006;29:1037–55. [69] Xue Y, Xie C, Wu J, Lie P, Wang J, Jiang Q. A three-dimensional miscible transport model for seawater intrusion in China. Water Resour Res 1995;31(4):903–12. [70] Younes A, Fontaine V. Hybrid and multi point formulations of the lowest order mixed methods for Darcy’s flow on triangles. Int J Numer Methods Fluids 2008. doi:10.1002/fld.178. [71] Younes A, Fontaine V. Efficiency of mixed hybrid finite element and multi point flux approximation methods on quadrangular grids and highly anisotropic media. Int J Numer Methods Eng 2008. doi:10.1002/nme.232. [72] Younes A, Ackerer P. Solving the advection–dispersion equation with discontinuous Galerkin and multipoint flux approximation methods on unstructured meshes. Int J Numer Methods Fluids 2008. doi:10.1002/fld.178. [73] Younes A, Ackerer P, Chavent G. From mixed finite elements to finite volumes for elliptic PDEs in two and three dimensions. Int J Numer Methods Eng 2004;59:365–88. [74] Younes A, Ackerer P, Lehmann F. A new mass lumping scheme for the mixed hybrid finite element method. Int J Numer Methods Eng 2006;67(1):89–107.