Anis Younes
email: younes@unistra.fr

Behshad Koohbor

Marwan Fahs

Hussein Hoteit

An efficient Discontinuous Galerkin -Mixed finite element model for variable density flow in fractured porous media

Keywords: Fractured porous media, variable density flow, advection-dispersion transport, mixed finite element, discontinuous Galerkin method, method of lines

come

An efficient Discontinuous Galerkin -Mixed finite element model for variable density flow in fractured porous media

Introduction

Variable-density flow (VDF) processes have been largely investigated in unfractured domains, but studies in fractured domains remain relatively scarce. VDF in fractured domains is encountered in several applications such as fractured coastal aquifers contamination by saltwater intrusion [START_REF] Grillo | Simulation of density-driven flow in fractured porous media[END_REF][START_REF] Werner | Seawater intrusion processes, investigation and management: Recent advances and future challenges[END_REF][START_REF] Sebben | Seawater intrusion in fractured coastal aquifers: A preliminary numerical investigation using a fractured Henry problem[END_REF][START_REF] Koohbor | Uncertainty analysis for seawater intrusion in fractured coastal aquifers: Effects of fracture location, aperture, density and hydrodynamic parameters[END_REF][START_REF] Hosseini | Numerical modeling of density-driven solute transport in fractured porous media with the extended finite element method[END_REF][START_REF] Etsias | Laboratory and numerical investigation of saline intrusion in fractured coastal aquifers[END_REF], nuclear waste management [START_REF] Follin | A transmissivity model for deformation zones in fractured crystalline rock and its possible correlation to in situ stress at the proposed high-level nuclear waste repository site at Forsmark, Sweden[END_REF], CO2 storage [START_REF] Jafari Raad | Onset of density-driven instabilities in fractured aquifers[END_REF][START_REF] Kim | Density-Driven Convection in a Fractured Porous Media: Implications for Geological CO 2 Storage[END_REF][START_REF] Shafabakhsh | Convectivereactive transport of dissolved CO2 in fractured-geological formations[END_REF][11][12][13][14][15], thermal convection [16,17] or oil and gas production [18].

The propagation of dense plumes in a porous medium is known to be highly affected by heterogeneity [19]. This is even more significant in the case of fractured porous media where the plume migration can be highly impacted by the characteristics of the fractures such as location, size, density, orientation and aperture [20][21][22][23]. Modeling flow and transport in fractured porous media can be performed using either the implicit or the explicit representation of the fractures [24]. In the implicit approach (e.g. equivalent porous medium and dual porosity approaches), the effects of fractures are considered using single or multicontinuum models. In the single continuum model, an equivalent bulk permeability for the porous medium is employed which depends on the properties of the fractures and the fracture network [25,26]. In multi-continuum models, two or more superimposed media with different flow and/or transport equations are considered [27][28][29].

The explicit approach uses a discrete representation of the fractures. Among different variants of this approach, the discrete fracture network (DFN) model considers the matrix continuum as completely impermeable. Therefore, flow and transport processes only occur through interconnected fractures [30,31]. In this work, we consider the discrete fracture matrix (DFM) model where flow and transport processes occur in both fractures and matrix continua [24,32,33]. Fractures are considered thin regions filled with porous material with different properties from those of the porous medium enclosing them [START_REF] Grillo | Simulation of density-driven flow in fractured porous media[END_REF]. They are assumed to be codimension one with respect to the dimension of the surrounding matrix [34][35][36]. The main goal of this paper is to develop an accurate and efficient DFM model for the simulation of VDF in fractured domains by combining advanced methods for space discretization with high-order time integration schemes. Notice that, the DFM models can only handle a limited number of fractures (of the order of thousands) due to the limitation of computational resources [37]. One approach to overcome this limitation is to consider only the major fractures and use an upscaling technique to take into account the effect of the minor fractures [38]. Furthermore, DFM models require an unstructured mesh conform to the geometry of fractures which may complicate the development of the gridding algorithm. To avoid gridding problems, Embedded Discrete Fracture Model (EDFM) and projection-based EDFM (pEDFM) have been recently developed [39][40][41] where fractures are coupled with the matrix domain by modifying the transmissibility factors to reflect the geometry and properties of the fractures [42].

Because of the high permeability contrast between the matrix and the fractures, appropriate numerical methods should be employed for an accurate simulation of the flow in fractured domains. Further, the employed numerical method should fulfill the local (elementwise) mass conservation property for problems involving the transport of contaminants [43]. Cellcentered methods, like finite volumes (FV), finite differences (FD), Discontinuous Galerkin (DG) and mixed finite element (MFE) methods ensure local mass balance and are well suited for modeling VDF in highly heterogeneous domains such as fractured porous media.

However, these methods require computing water and contaminant matrix-fracture exchange fluxes which may require extra variables [33]. A practical way to overcome this difficulty is to employ the fracture cross-flow equilibrium assumption where the state unknowns (i.e. pressure and concentration) in the fracture and in the adjacent grid-cells are supposed to be equal [44,45]. This assumption requires the matrix grid cells next to the fractures to be sufficiently small which can hamper the efficiency of DFM models [46]. For instance, Hoteit and Firoozabadi [45] found that the matrix grid cell size next to the fractures should be of the order of ten centimeters in a fractured reservoir of length-scale in kilometers. This constraint on the mesh size is avoided by the hybrid formulation of the MFE method, which has the traces of pressures at edges as unknowns [47][48][49]. Hence, with this formulation, the cross-flow equilibrium is only assumed across edges by imposing the matrix pressure (respectively concentration) at an edge and the pressure (respectively concentration) at the fracture which coincides with that edge to be the same. This powerful feature renders the hybrid MFE method an excellent candidate for the discretization of the flow in fractured porous media.

The MFE has many other interesting features, notably: (i) the method is well adapted for high heterogeneous domains with full permeability tensors, (ii) it ensures an exact local mass balance and (iii) it yields more accurate fluxes than the conventional finite element (FE) and FV methods, particularly for unstructured grids [50,51]. Indeed, the classical cell-entered FV methods use the two-point flux approximation (TPFA) which may not converge unless the grid is K-orthogonal [52]. The MFE was also shown to be more efficient than the Multi-Point Flux Approximation method for anisotropic media and unstructured meshes [53]. Further, the MFE method was shown to be superior to the Local Discontinuous Galerkin (LDG) method for non-smooth grids [37]. The main drawback of the MFE method is the size of the linear system which is equal to the number of the mesh edges. A lumped formulation of the MFE method with improved monotonicity has been developed by Younes et al. [54]. Unlike the standard method, the lumped formulation has the advantage of keeping the time derivative continuous and hence, allows the employment of high-order time integration methods via the method of lines (MOL). A fully MFE discretization has been proposed for both flow and transport in unsaturated fractured porous media in Younes et al. [55].

In this work, we develop a MFE-DG model for VDF in a low permeable rock matrix with a limited number of embedded highly-conductive fractures where the fluid front movement is much faster than in the matrix continuum. The model is developed for 2D space and can be extended to 3D space. The fractures grid cells are (D-1)-dimensional grid-cells in a Ddimensional domain, thus for 3D space, fractures are 2D (i.e. planar) interfaces [46].

The lumped formulation of MFE is applied for the discretization of the flow equation both in the fracture and the matrix continua. The advection-dispersion transport in the porous matrix is discretized using a combination of the lumped MFE method and an upwind edge-centered FV scheme to produce a monotonic solution. In the fracture network, water and contaminants travel at much greater speeds than in the surrounding porous matrix [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF] making the transport process advection-dominated and its differential equation hyperbolic [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. Hyperbolic equations are characterized by moving sharp fronts that classical numerical methods (such as standard FE or FV methods) fail to capture, yielding solutions with non-physical oscillations and/or strong numerical diffusion [START_REF] Huyakorn | Computational Methods in Subsurface Flow[END_REF]. To avoid these phenomena, the transport in the fracture network is solved in this work with the DG finite element method. The DG method conserves mass at the element level and produces low numerical dispersion compared to classical methods. The method was shown to be clearly superior to existing finite element methods for hyperbolic systems [START_REF] Siegel | solution of the advection-diffusion equation using a combination of discontinuous and mixed finite elements[END_REF][START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF]. Time integration of the nonlinear MFE-DG system is performed with high-order methods using the MOL [START_REF] Younes | Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods[END_REF]. The main idea of this method is to convert partial differential equations into a system of ordinary differential equations by discretizing the space derivatives while keeping the time derivative continuous. The resulting system of ordinary differential equations is then solved with the DASPK time solver. This solver optimizes the time step length and the order of the temporal integration to reduce the computational time while maintaining accuracy. DASPK has been successfully employed for VDF problems in unfractured porous media [START_REF] Younes | Modeling variable-density flow in saturated-unsaturated porous media: An advanced numerical model[END_REF].

In this work, we show how we construct an efficient and accurate DFM model for VDF in fractured porous media with cross-flow equilibrium assumed only across fractures using:

-The lumped hybrid MFE method for the spatial discretization of the flow in the matrix and in the fracture continua; -The upwind MFE scheme for the spatial discretization of the transport in the matrix continuum;

-The DG method for the spatial discretization of the advection-dominated transport in the fractures; -High-order time integration techniques and efficient adaptive time stepping for the obtained coupled nonlinear flow transport system.

The paper is organized as follows: First, the partial differential equations describing the flow and the transport with density and viscosity variations in fractured porous media are presented. Second, the new MFE-DG numerical model is described. The lumped MFE method is detailed for the discretization of the flow in the porous matrix and in the fractures.

The upwind hybrid MFE method is presented for the discretization of the advectiondispersion transport in the porous matrix and the DG method is detailed for the advectiondominated transport in the fractures. Finally, we present numerical examples in unfractured and fractured media to demonstrate the validity, efficiency and accuracy of the developed MFE-DG model by comparison against semi-analytical and standard FE solutions.

Governing Equations

The fractures are assumed infilled by a porous medium having different properties from those of the porous matrix. Assuming the Boussinesq approximation is valid, VDF in the matrix and in the fracture network can be ruled by the following mass conservation equation of the fluid (in the absence of source/sink term):

0 h S t       q (1)
and Darcy law:

0 0 hz            qK (2)
where S is the mass storativity related to head changes [L -1 ], The transport of contaminant in the porous matrix and through the fracture branches is ruled by the following advection-dispersion equation (in the absence of reaction terms):

  0 C CC t           qD (3)
where  is the porosity, C is the relative concentration [-] and D is the dispersion tensor given by: The flow and transport equations are coupled by density and viscosity relationships as function of concentration. Although, when salt is added to water, the volume of the solution may change due to the contraction of solutes, the fluid density remains almost proportional to the solute concentration [START_REF] Guo | Density Slopes in Variable Density Flow Modeling[END_REF]. Assuming isothermal conditions and incompressible fluid, density is often represented by the linear relationship:

  D I q q q q I m L T T D/         (4)
  0 1 0 C        (5)
where 1  is the density of the injected contaminant.

Note that there are also nonlinear forms of approximation for the relationship between fluid density and solution concentration [64]. The viscosity of the mixture,

  C  
, is mostly represented by a nonlinear function of the concentration [64].

The fractures can be considered as open empty channels without a porous medium. For such a case, the permeability of the fracture is replaced by (see [START_REF] Adler | Fractured Porous Media[END_REF]), which is often referred as the cubic law [START_REF] Witherspoon | Validity of Cubic Law for fluid flow in a deformable rock fracture[END_REF]. In this case, dispersion in the fractures reduces to numerical diffusion

( 0 L   ).
In the following, the nonlinear system (1)-( 5) is solved using a combination of MFE and DG methods on 2D triangular elements for the porous matrix and on a set of 1D linear elements for the fracture network.

The new MFE-DG model

In this section, we provide details on how MFE and DG methods can be efficiently combined to solve coupled flow and transport in a fractured porous medium. More specifically, original details in space discretization are developed to show how (i) the lumped MFE method, developed in Koohbor et al. [START_REF] Koohbor | An advanced discrete fracture model for variably saturated flow in fractured porous media[END_REF] for unsaturated fractured flow, can be extended to VDF in fractured porous media, (ii) the upwind MFE scheme, recently developed in Younes et al. [START_REF] Younes | A robust upwind mixed hybrid finite element method for transport in variably saturated porous media[END_REF], can be employed for transport in the rock matrix and combined with the DG method for advection dominated transport in fractures, (iii) the upstream DG concentrations can be calculated at the intersection of several fractures and (iv) the continuity of the dispersive flux can be ensured at the nodes of the intersection of several fractures with different apertures.

The plan of this section is as follows: First, the MFE is developed for the spatial discretization of the flow equation in the matrix continuum and in the fracture continuum. Then, the transport equation in the matrix continuum is discretized using the upwind hybrid MFE scheme. Finally, the DG method is used for the spatial discretization of the advection dominated transport in the fractures.

MFE discretization of the flow in the porous matrix

We use the lowest order Raviart-Thomas space (RT0) approximation for the Darcy velocity q inside each triangular element E as:

  1 0 0 E E E E hz                E E E i i i
K q w w w [START_REF] Follin | A transmissivity model for deformation zones in fractured crystalline rock and its possible correlation to in situ stress at the proposed high-level nuclear waste repository site at Forsmark, Sweden[END_REF] which becomes [START_REF] Younes | From mixed finite elements to finite volumes for elliptic PDEs in two and three dimensions[END_REF]  

0 0 E E E E E i , j j E i E i j B Q h T h z T z         (8)
where E

B

is the elemental matrix of terms  

E E E , E E EE i i E E i , j j j j Q h z B T h T z                          (9)
where

1 E E , i i , j j B    
.

The lumped formulation of MFE [54] is then used as follows  Combination of Eq. ( 1) and Eq. ( 9) provides the steady state flux

E i Q
as function of head traces at edges,

1 0 0 EE E ij E , E E E i , j j j E i j Q B T h T z                 (10)
where EE i i     The mass conservation Eq. ( 1) is written over the lumped region i R (blue area in x , as following (see Figure 1 for notations), 0 33

EE E E E ' E ' im E E ' ij ik m j m k E E ' d T h d T h S S Q Q Q Q d t d t       (11)
where

  1 3 E E E ij j i Q Q Q 
is the interior flux edge evaluated using Eq [START_REF] Etsias | Laboratory and numerical investigation of saline intrusion in fractured coastal aquifers[END_REF].  The mass conservation Eq. ( 11) can be rewritten as an equation of continuity of the water flux between E and E  as following,

1 0 0 0 3 E E ' im EE E ij E E , E E Ei i i , j j j E E j QQ E d T h Q B T h T z S dt                        (12)
The Eq. ( 12) is written for all mesh edges which are not coincident with a fracture. The continuity of head traces  

EE i i m T h T h T h  
is implicitly imposed at each edge sharing the two adjacent elements E and E'.

Note that with the lumped MFE formulation Eq. ( 12), the time derivative remains continuous which allows employment of efficient high order time integration methods via the MOL.

Spatial discretization of the flow in the presence of a fracture

We consider a one-dimensional flow in a fracture l located at the interface between the elements E and E  (see Figure 2).

  0, xL  l .
The mass conservation Eq. ( 1) integrated over the fracture l writes

12 0 EE im d T h e L S q q Q Q dt       l ll lll ( 14
)
with Th l the mean head in the fracture l and

S l the storage coefficient in l , E i Q and E m Q  are
the matrix-fracture exchange fluxes, respectively, from the elements E and E'.

Using Eq. ( 13), the variational formulation of the Darcy law in the fracture l writes

  0 0 0 K L f k j i j i i j q w w T h N h T z N z e             l l l l l l l l l l (15)
where Applying the trapezoidal rule for the left term in Eq. (15) (see [START_REF] Koohbor | An advanced discrete fracture model for variably saturated flow in fractured porous media[END_REF]) yields

  0 0 2 f i i i K q e T h N h T z N z L           l l l l l l l l l (16) which yields   1 2 1 2 2 2 f K q q e T h N h N h L     l l l l l l l l ( 17
)
Using the cross-flow equilibrium assumption in the fracture by imposing equality of the head in the fracture and in the matrix edge  

EE im T h T h T h   l
, the mass conservation Eq. ( 14) becomes   11 00 00 12 2 20 33

E E E ' E ' i j m j E , E E E ', E ' E ' E E ' i , j j j m , j j j E E ' jj E f E i EE B T h T z B T h T z K EE d T h e T h N h N h S S e L S L d t                                                             l ll l l l l l l ( 18
)
To close the system (18), mass conservation is written at each intersection of fracture branches. At the node i of head where * i Q is the pumped/injected quantity at the node i of the fracture network.

Using Eq (16), we obtain

  0 0 2 f * i i i K e T h N h T z N z Q L                  l ll l l l l (20)
The head i Nh at the intersection node i writes

  0 0 2 2 f * ii i f K e T h T z N z Q L Nh K e L                       l ll l l l l l l l l (21)
The final system to solve is obtained by substituting Eq. ( 21) for (18). This system is solved for the head traces at all mesh edges (except Dirichlet boundary edges). Hence, the flow system in the presence of fractures has the same DOF than for unfractured domains.

Discretization of the transport equation in the porous matrix

The RT0-MFE method is employed for the solution of the advection-dispersion transport equation in the porous matrix. An upwind scheme is developed to provide a monotonic solution. Upwind schemes deliver stable solutions but are known to introduce some artificial numerical diffusion. In the matrix continuum, this drawback is limited since the transport is generally diffusive due to the slow flow caused by the low permeability of the matrix compared to that of fractures.

The transport equation ( 3) is written in the following form:

0 C C t C                  % % qq qD (22)
The concentration is considered to be piecewise constant and the dispersion vector q% is piecewise linear. Inside the triangular element E , q% is given by:

EE jj j Q qw   % % ( 23
)
where

j EE j d j E Q.    % %
q η is the dispersive flux across the edge As with the flow equation ( 9), the dispersive flux can be written as

1 E E E , E i i E i , j j j Q C B T C     %% % (
  1 E i , j E j E B    % EE i

D w w

The upwind hybrid MFE scheme [START_REF] Younes | A robust upwind mixed hybrid finite element method for transport in variably saturated porous media[END_REF] is employed following the main steps below  Similarly to Eq. ( 10), the steady state dispersive flux

E i Q % writes 1 EE ij E E E , E i j i , j j E jj Q T C B T C      %% % % % (25) with ,1 , EE i i j j B     % %
and

EE i i    %% .
 The integration of the advection dispersion equation Eq. ( 22) over the lumped region i R (see Figure 1) gives:

    0 3 E E E E E E E E E i E ij ij ik ik ij ik i i E E E TC Q T C Q T C Q Q T C Q ... t               % ( 26
)
where

  1 3 E E E ij j i Q Q Q 
is the interior water flux between edges

  1 E E E E E ij ij i ij j T C T C T C     (27)
with

1 E ij   for an outward flux   0 E ij Q  , else 0 E ij   .
Denoting as ± E i Qt , the total (advective+dispersive) flux at the edge i E  of the element E, the mass conservation Eq. ( 26) is expressed as the continuity of the total flux between the two adjacent elements E and ' E , as

± ± ±         1 0 11 3 E E ' im EE E E ij E , E E E E E E E E E i i , j j ij ij i j ik ik i k E i E j Q t Q t E TC Q t B T C Q T C T C Q T C T C t                               %% % % (28)

Spatial discretization of the transport equation in the fracture network

Because of the high-permeability of fractures, the flow across the fracture network can be much faster than in the matrix. As a consequence, the transport through fracture branches can be highly advection dominated and its equation tends to be hyperbolic. Hyperbolic equations can present solutions involving moving sharp fronts that classical methods fail to capture. Indeed, for advection dominated transport, standard numerical methods (e.g. FV, FD or FE methods) generally provide solutions with significant unphysical oscillations and/or strong numerical diffusion [START_REF] Huyakorn | Computational Methods in Subsurface Flow[END_REF]. For such situations, the DG method is clearly superior to existing FE methods [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF] as it is well adapted for hyperbolic systems. Indeed, the DG method yields a high-resolution scheme which is strictly conservative at the element level and can accurately reproduce solutions involving sharp fronts [START_REF] Siegel | solution of the advection-diffusion equation using a combination of discontinuous and mixed finite elements[END_REF]. The DG method is detailed hereafter for the transport of dense contaminants in the fractures network. An implicit time discretization is employed to avoid time step limitation caused by the CFL condition of explicit schemes.

The integration of the transport equation ( 22) over the fracture element l of aperture e l and length L l , taking into account matrix-fracture exchanges writes

± ±   1 EE im C C Q t Q t t e L            ll % qq (29)
where the right hand side term

± ±   1 EE im Q t Q t eL   ll
designates the contaminant flux per unit area exchanged between the fracture l and the matrix elements E and E  .

For the facility of the equation explanation, x is denoted as the local coordinate from the first to the second node of the fracture element l . A linear discontinuous approximation is assumed for the concentration in the fracture l as

      12 ,| x C x t x T C x C   l l l l l l ( 30
)
in which TC l and x C l are the two degrees of freedom (DOF) corresponding to the concentration at the fracture centre x l and its deviation along the fracture direction with the following interpolation functions:

    12 1, x x x x     ll l (31)
The variational formulation of Eq (29) writes

      ± ±   0 0 0 0 0 0 EE L L L L L L im i i i i i i Q t Q t C e e C e C e C e tL                              l l l l l l l l l l l l l % q q q q (32)
The calculations of terms of Eq. (32) are developed in the following.

The mass term in Eq. (32) writes

1 3 0 0 0 12 L i x M TC L t C ee L t C t                        l l l l l l l l l 1 4 4 2 4 43 (33)
The second term in Eq. (32) corresponds to the element boundary term:

  1 1 2 2 0 0 L ,* ,* i i i x x L e . C q C q C         l l l l l l l q (34)
where ,* i C l is the upstream concentration at the extremity i of l , expressed by

  1 ,* i ,o u t i i i ,in i C N C N C     l l l l l (35) with i ,in i x N C T C C   l l l
l the interior concentration at the node i of the fracture l , using the following notations:

  12 1 , 2 10 ,, 00
22 i i i i if q LL if q              l l l ll (36)
Thus, in the case of an outflow at the node i of l , the nodal concentration is

,* i i ,in C N C  ll ,
whereas, in the case of an inflow, the nodal concentration is noted , we use the mass conservation equation at the node i, shared by n fractures (see Figure 3) which writes:

1 0 n f f ,* ii f qC    (37)
Substituting Eq. (35) into Eq. (37) gives

1 f f f i ,o u t i i i ,in ff l ii f N C q N C q      (38) where f f f f i ,in i x N C T C C   .
Finally, Eq. (34) can be written in the following form

          2 12 1 1 2 2 1 1 1 2 2 2 22 0 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 12
11 11 The continuity of the dispersive fluxes of all fracture branches sharing the node i (see Figure 3) writes 2 0

I L i x M ,o u t VV q q q q TC eC C q q q q q N C q N                                                     l l l l l l l l l l
ff f i f f eD T C N C L      (45)
Which yields

f f f f f i f f f f e D T C L NC e D L    (46)
Substituting Eq. ( 46) into Eq. (43) yields

2 f f f f f i f f f f e D T C L eD q T C e L D L          ll ll l % (47)
Finally, the right hand side of Eq. (32) writes

± ±   ± ± 0 0 EE EE L im im i Q t Q t Q t Q t L            l l (48)
and the final local system for the transport in the fracture l writes

        ± ± 12 12 1 2 3 4 1 2 0 0 EE im ,o u t ,o u t x x TC TC q q Q t Q t t M M V N C V N C C C t                               l l ll ll l l %% (49)
where the local matrices 1 M and 2 3 4

2 3 4 I M M M M   
are defined in Eq. (33), Eq. (39), Eq.

(40) and Eq. ( 41). The local vectors 1 V and 2 V are defined in Eq. (39).

The final transport system is obtained by substituting the total fluxes ± E i Qt and ± E m Qt  expressed by Eq. (28) into the system (49) and by imposing the concentration of the fracture l to be equal to that of element edges coincident with

l   EE im T C T C T C   l .
Hence, the unknowns of the transport system are the traces of concentration at all the edges of the mesh plus the concentration deviation along the fracture direction at all fracture branches.

The temporal discretization of the nonlinear flow-transport system

One of the main advantages of the developed formulation is that the transient terms in the discretized flow and transport systems remain continuous which allows the employment of efficient high-order methods for time integration. Thus, the flow system formed by Eq. ( 12) or Eq. ( 18), in which we substitute Eq. ( 21), is combined with the transport system formed by Eq. ( 28) and Eq. ( 49) to form a single nonlinear implicit system of ordinary differential equations (ODEs) of the form

  0 F t , ,   yy (50) Where       11 1 i i x i ,.

.,n b _ e d g e s i ,..,n b _ e d g e s ,..,n b _ fr a c tu r e s T h , T C , C

      l l y is the vector of unknowns formed by (i) the head traces at all the edges of the mesh (except boundary edges with prescribed head), (ii) the concentration traces at all edges of the mesh (except boundary edges with prescribed concentration) and (iii) the concentration deviation in all fracture cells. Thus, the number of unknowns for the flow-transport system is approximately twice the number of edges, plus the number of fractures.

The time integration of the nonlinear system of Eq. ( 50) is performed with higher-order methods using the DASPK time solver [START_REF] Van Keken | DASPK: A new high order and adaptive time-integration technique with applications to mantle convection with strongly temperature-and pressure-dependent rheology[END_REF]. Higher-order methods use larger time steps and require less effort in the nonlinear solver than lowest-order methods [START_REF] Farthing | Mixed finite element methods and higher-order temporal approximations[END_REF]. Further, higherorder methods are often combined with an efficient automatic time stepping scheme that allows to select the optimized time-step size that maintains a given temporal discretization error [START_REF] Farthing | Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow[END_REF][START_REF] Tocci | Accurate and economical solution of the pressurehead form of Richards' equation by the method of lines[END_REF][START_REF] Kavetski | Adaptive backward Euler time stepping with truncation error control for numerical modelling of unsaturated fluid flow[END_REF]. DASPK is a sophisticated time integration solver for ODE and DAE systems. It uses a variable order time integration method (up to fifth order) based on the Fixed Leading Coefficient Backward Difference Formulas (FLCBDF) which has good stability properties [START_REF] Kees | Higher order time integration methods for two-phase flow[END_REF]. With the FLCBDF method, the Newton iteration matrix can be reused for more steps than in a fully variable-step approach [START_REF] Van Keken | DASPK: A new high order and adaptive time-integration technique with applications to mantle convection with strongly temperature-and pressure-dependent rheology[END_REF][START_REF] Hindmarsh | Algorithms and software for ordinary differential equations and differential-algebraic equations[END_REF].

DASPK involves a prediction Lagrange polynomial of order k using points   0 ,..,

,

n i n i ik t   y .
A corrector Lagrange polynomial of order k is then formed by interpolating the predictor polynomial at k equally spaced steps before (see [START_REF] Brenan | The numerical solution of initial value problems in differential-algebraic equations[END_REF] for more details).

Differentiating the previous system with respect to With DASPK, the linear systems arising at each time step can be solved with iterative or direct methods. In this work, we choose the preconditioned Krylov iterative method which is more adapted for large linear systems. A preconditioner matrix, based on the ILU (Incomplete Lower Upper) decomposition of the Jacobian matrix J , is computed and used over as many time steps as possible.

To reduce the calculation costs of the evaluation of the Jacobian matrix J , which can be very high, the terms of J are calculated using the finite difference approximation. Furthermore, we implement in DASPK the column grouping technique (see [START_REF] Curtis | On the Estimation of Sparse Jacobian Matrices[END_REF][START_REF] Hindmarsh | Large ordinary differential equation systems and software[END_REF]) to reduce the number of evaluations of the residue. With this technique, variables are perturbed by group, knowing the structure and sparsity of the Jacobian matrix.

DASPK, automatically adapts both the order of integration and the time step length to minimize the computational effort while keeping the temporal truncation error small. In this work, the maximum order is fixed to 5 and both relative and absolute time errors are maintained under a fixed tolerance of

Numerical experiments

In this section, three test problems are investigated to assess the validity, efficiency and

Flow and transport in a 2D domain with a ''+''-shaped barrier/fracture network

A homogeneous 2D problem is considered (see Figure 4a), where a ''+''-shaped barrier/fracture network is located in the middle of a 9m x 9m square domain. No-flow boundary conditions are applied at the top and bottom. At the left vertical side, a no-flow boundary occurs, except between y=4m and y=5m, where the pressure is fixed to P=1 and the concentration is set to C=1. The right vertical side has a fixed pressure P=0 and a zero diffusive boundary condition for the transport. The rock matrix has a permeability Intel Xeon E-2246G processor and 32 GB memory. Thus, the new 1D-2D model is around 17 times more efficient than the 2D-2D model. This highlights the great benefit of robust and accurate 1D-2D models to simulate flow and transport processes in fractured porous media.

The Henry saltwater intrusion problem

The Henry problem [START_REF] Hr | Interfaces between salt water and fresh water in coastal aquifers, US Geological Survey Water-Supply Paper 1613-C[END_REF] is a popular hypothetical saltwater intrusion problem [START_REF] Fahs | A Generalized Semi-Analytical Solution for the Dispersive Henry Problem: Effect of Stratification and Anisotropy on Seawater Intrusion[END_REF]. It represents a vertical cross-section of a confined coastal aquifer in which the inland freshwater coming from the left boundary is in equilibrium with seawater intruded from the right side due to its higher density (Figure 5). Sebben et al. [START_REF] Sebben | Seawater intrusion in fractured coastal aquifers: A preliminary numerical investigation using a fractured Henry problem[END_REF] -zero concentration gradient along the top and bottom Table 1. Parameters and boundary conditions for the original Henry problem.

The Henry problem is simulated using the parameters and boundary conditions described in Table 1. The analytical solution under this set of parameters has been developed by Henry [START_REF] Hr | Interfaces between salt water and fresh water in coastal aquifers, US Geological Survey Water-Supply Paper 1613-C[END_REF] using the Fourier-Galerkin (FG) method. The FG solution yields the steady-state isochlors positions by expanding the salt concentration and the stream function in double Fourier series. Ségol et al. [START_REF] Ségol | Classic groundwater simulations: proving and improving numerical models[END_REF] showed that the initial semi-analytical solution developed by Henry [START_REF] Hr | Interfaces between salt water and fresh water in coastal aquifers, US Geological Survey Water-Supply Paper 1613-C[END_REF] was not accurate because it did not contain enough terms (only 78) in the Fourier series. An accurate semi-analytical solution was developed later using more coefficients (138) in the Fourier series by Segol [START_REF] Ségol | Classic groundwater simulations: proving and improving numerical models[END_REF] and Zidane et al. [START_REF] Zidane | The Henry semianalytical solution for saltwater intrusion with reduced dispersion: Henry semianalytical solution with reduced dispersion[END_REF]. However, as pointed out by Kolditz et al. [START_REF] Kolditz | Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models[END_REF], the standard Henry problem (Table 1) involves an unrealistic large amount of dispersion and, as thus, cannot be used to check whether a model can represent field situations with relatively narrow transition zones. A new semi-analytical solution has been developed by Younes and Fahs [START_REF] Younes | A semi-analytical solution for saltwater intrusion with a very narrow transition zone[END_REF] for the challenging case of saltwater intrusion with a narrow transition zone. The semi-analytical solution was obtained for a very small diffusion coefficient that is 20 times less than that of the original Henry problem. The new semi-analytical solution required around 6330 terms in the Fourier series to accurately represent the solution. Younes and Fahs [START_REF] Younes | A semi-analytical solution for saltwater intrusion with a very narrow transition zone[END_REF] showed that, compared to the standard Henry problem, the new semi-analytical solution is more suited for benchmarking density-driven flow codes. In the following, the original and the new semi-analytical solutions are used for the validation of the MFE-DG numerical model in the case of high and low dispersion. For both cases, two configurations are explored. In the first configuration, the domain does not contain any fracture. Thus, flow and transport occur in the matrix continuum which has a permeability 0 m kk  . In the second configuration, the matrix is considered completely impermeable, thus, density-driven flow occurs only within the fracture network.

The results of the first configuration (unfractured domain) with the MFE-DG numerical model on a general triangular mesh formed by 1442 triangles are depicted in Figure 6. In this figure, saltwater intrudes from the right side and reaches equilibrium with the inland freshwater flow. The intrusion is more pronounced near the bottom because of density effects.

The transition zone is more significant for the standard Henry problem (Figure 6a) than for the low-dispersion problem (Figure 6b). On the other hand, density-dependent effects (and as a consequence, the amount of intruded saltwater) are more significant for the low-dispersion case (Figure 6b) than for the original Henry problem (Figure 6a). In the litterature, the comparison between the semi-analytical and numerical results are often limited to the position of the 0.25, 0.5 and 0.75 isochlors. A more complete and sensitive comparison is performed here by investigating the agreement between analyical and numerical concentrations on two concentration profiles: a horizontal profile at y=0.2m and a vertical profile at x=1.5m (dashed lines in the Figure 6). The semi-analytical and numerical concentration profiles are plotted in the Figures 7a and7b for the standard Henry problem and in the Figures 7c and7d for the low-dispersion problem. As expected, the concentration profiles are sharper for the low-dispersion problem than for the standard Henry problem. A very good agreement is observed between the horizontal numerical and semi-analytical profiles (Figure 7a) as well as between the vertical numerical and semi-analytical profiles (Figure 7b). This validates the results of the MFE-DG code in the case of high dispersive density driven flow in unfractured domains. Likewise, a very good agreement is also observed, for the low-dispersion Henry problem, between the horizontal numerical and semianalytical profiles (Figure 7c) as well as between the vertical numerical and semi-analytical profiles (Figure 7d) which validates the MFE-DG results in the case of low-dispersion density-driven flow in unfractured domains. 0,0 0,5

MFE-DG solution Semi-analytical solution 0,0 0,2 0,4

MFE-DG solution Semi-analytical solution 0,0 0,5 1,0

MFE-DG solution Semi-analytical solution 0,0 0,2 0,4 0,6 0,8

MFE-DG solution Semi-analytical solution In Figure 8, the concentration profile at y=0.2m and at x=1.5m do not cross the impermeable matrix and remain in the fracture network. The semi-analytical and numerical horizontal (at y=0.2m) and vertical (at x=1.5m) concentration profiles are plotted in Figures 9a and9b for the standard Henry problem and in Figures 9c and9d for the low-dispersion problem. For the standard Henry problem, a very good agreement is observed between the horizontal numerical and semi-analytical profiles (Figure 9a) but a small discrepancy is observed between the vertical numerical and semi-analytical profiles (Figure 9b). For the low dispersion problem, characterized by a sharp concentration front, the discrepancy between semi-analytical and numerical results becomes more significant and is observed for both horizontal (Figure 9c) and vertical (Figure 9d) profiles. The mismatch between numerical and semi-analytical results observed in Figure 9 is related to the fact that the semi-analytical solution is developed for a continuous 2D domain whereas

(a) (b) (c) (d)
the numerical solution is obtained on a set of 1D fractures. Thus, the fracture network could be not sufficiently dense to reproduce the processes occurring in the continuous 2D domain.

To avoid this problem, both the original and low-dispersion Henry problems are now simulated on a dense fracture network using 50 nf 

. The results of the simulations are depicted in Figure 10. A very good agreement can be observed between the numerical and semi-analytical horizontal (Figure 10a) and vertical (Figure 10b) concentration profiles for the original Henry problem. Likewise, a very good agreement is observed, for the low-dispersion Henry problem, between the numerical and semi-analytical horizontal (Figure 10c) and vertical (Figure 10d) concentration profiles. The results show that in the pure diffusion case, (i) the semi-analytical and DFN models yield similar steady-state salt distributions provided the fracture network is dense enough and (ii) a higher fracture density is required with the DFN model for the low diffusion case to approximate well the semi-analytical salt distribution. These results validate the MFE-DG numerical model for both high and lowdispersion density-driven flow in the fracture network. to obtain a stable flow with high density and viscosity variations [START_REF] Loggia | Phase diagram of stable miscible displacements in layered porous media[END_REF]. In this case, the density and the viscosity of the mixture are given by: The test problem is simulated with the MFE-DG model using a coarse mesh of 1950 triangles and 165 fracture edges (Figure 11a). The total flow-transport system has 6165 unknowns formed by 3000 pressure traces (pressure at mesh edges), 3000 concentration traces (concentration at mesh edges) and 165 concentration deviations along fractures (due to the P1-DG approximation in fractures). The simulation is performed for a final time of 8000s and the output breakthrough curve is compared to a standard FE solution obtained with COMSOL Multiphysics®. The results of the two models for the tracer configuration are depicted in Figure 12 shows that the output concentration calculated with the MFE-DG model reaches a first plateau at around 1200s due to the fast arrival of the front of concentration from the right high-permeable fracture. A second plateau is observed at around 3000s due to the arrival of the front of concentration from the low-permeable left vertical fracture. Note that horizontal fractures remain free of contaminant and there is no transfer between the two vertical fractures as shown in Figure 11b.

When used with the same coarse spatial discretization, the breakthrough curve of the standard FE model shows a significant delay. Further, the concentration distribution shows strong unphysical oscillations (with negative concentrations as low as -0.4) at regions of high concentration gradients (Figure 13a). Although they do not disappear completely, these oscillations are strongly reduced when the simulation is performed with a very fine mesh formed by 6800 elements (Figure 13b). Further, the breakthrough curve obtained with the standard FE model on the very fine mesh is in good agreement with the MFE-DG coarsemesh solution (Figure 12). These results show the efficiency and accuracy of the MFE-DG model as compared to the FE model for the simulation of advection-dominated transport in the heterogeneous fractured porous media. The second configuration of this test problem is highly nonlinear because of the significant density and viscosity contrasts between the injected and the displaced fluids. For this configuration, the displacement is unstable due to density contrast (gravity number >0) and stable due to viscosity contrast (Mobility ratio < 1). For this configuration, the breakthrough curve obtained with the MFE-DG model depicts a single diffusive front, thus one cannot distinguish the arrival of the first front (traveling in the right high-permeable fracture) from the second (traveling in the left low-permeable fracture). This phenomenon is due to the lateral transfer between the left low-permeable and the right high-permeable vertical fractures through the horizontal fractures as can be seen in Figure 11c. This transfer is induced by the high viscosity of the injected fluid compared to that of the displaced one which implies a reduction of the hydraulic conductivity when the injected fluid advances (the hydraulic conductivity is proportional to the ratio of permeability to viscosity).

As a consequence, in the first stages, the velocity inside the high-permeable fracture reduces more than in the low-permeable fracture as the injected fluid advances further. Thus, a lateral transfer takes place to the low-permeability fracture where the front is delayed. At later stages, both fractures are filled with the injected fluid and the lateral transfer between the two vertical fractures vanishes.

The breakthrough curves obtained with the standard FE and the MFE-DG models for this highly nonlinear problem are compared in Figure 15. As for the tracer case (Figure 12), a significant difference is observed between FE and MFE-DG concentration profiles when used on the same coarse mesh. When used on a finer mesh of around 3.5 times more elements, FE yields a breakthrough curve in good agreement with that obtained by the MFE-DG model on the coarse mesh. These results validate the developed MFE-DG model and point out its efficiency and accuracy for solving flow and transport with high density and viscosity variations in heterogeneous fractured porous media.

Conclusion

In The new model uses the lumped formulation of the MFE method to solve the flow on both 2D-matrix elements and 1D-fracture elements. The transport in the matrix continuum is discretized using the lumped MFE method combined with an upwind edge centered scheme in order to improve its stability. Since the flow in the fracture network can be fast and the transport can be strongly dominated by advection, the DG method, which is well adapted to hyperbolic equations, is used to solve the transport in the fracture network. The time discretization of the coupled nonlinear flow-transport system is performed with high-order time integration methods via the MOL using the DASPK time solver.



  freshwater head [L], P the pressure [Pa], z is the upward vertical coordinate [L], t is the time [T], q is the Darcy velocity [LT -1 ], Kk is the hydraulic conductivity tensor [LT -1 ], 0  is the density of the displaced fluid [ML -3 ], g is the gravity acceleration [LT -2 ],  is the fluid dynamic viscosity [ML -1 T -1 ], k is the permeability tensor of the porous medium [L 2 ],  is the fluid density [ML -3 ]. The meaning of the two densities will be clarified below.



  are the longitudinal and transverse dispersivities [L], m D is the pore water diffusion coefficient [L 2 T -1 ] and I is the unit tensor.

  the node j faced to the edge j E  of the element E and E the area of E . The variational formulation of the Darcy's law (2) writes

Figure 1 )

 1 Figure 1), associated to the common edge

Figure 1 :

 1 Figure 1: The lumping region i R (blue area) associated to the interface of the elements E

Figure 2 :

 2 Figure2: The one-dimensional fracture l between adjacent elements E and E  .

  the head and the elevation of the extremity i of l .

  the mean concentration at the element E and at the edge i E  and % B the elemental matrix of terms

  corresponds to the upstream concentration:

  used for the discretization of the dispersive fluxC    qD %through the fracture l . Thus, the last term in the left-hand side of Eq. (32) becomes the longitudinal dispersivity, the molecular diffusion and the mean velocity in the fracture l .

  accuracy of the newly developed MFE-DG code for modeling flow and transport in fractured porous media. The first problem deals with linear flow and transport in a 2D domain involving a ''+''-shaped barrier/fracture network. The results of the developed 1D-2D model are compared to the results of a 2D-2D model where both fractures and matrix continua are discretized with small 2D mesh elements. The second problem is the Henry saltwater intrusion problem which has a semi-analytical solution that is commonly used for benchmarking codes dealing with VDF. The Henry problem is first used for the validation of the new MFE-DG code in the case of unfractured domains. Then, the semi-analytical solution of the Henry problem is used for the validation of the MFE-DG code in the case of fractured domains. The third test problem, inspired by the laboratory experiments of Loggia [80], deals with a stable flow with high density and viscosity contrasts in a fractured heterogeneous porous medium. This test case is simulated to investigate the validity and efficiency of the MFE-DG model by comparison against standard FE solutions obtained with the COMSOL multiphysics model.

2 . 2 .

 22 The horizontal fracture has 8m length, 0.04m aperture and a has the same length and aperture as the horizontal fracture, but a The porosity for the rock matrix, the horizontal fracture and the vertical barrier is fixed to 0.5. For the developed 1D-2D model, the domain is discretized with 6594 triangles for the matrix continuum, 46 lines for the horizontal fracture and 46 lines for the vertical barrier (Figure4a). For the 2D-2D model, a highly unstructured triangular mesh including 13753 elements is used with local mesh refinements to correctly discretize flow and transport processes in the fractures and barriers with 0.04m aperture (Figure4b).

Figure 4 :

 4 Figure 4: The ''+''-shaped barrier/fracture network problem: (a) the 2D-2D and (b) the 1D-

  simulated the velocity-dependent dispersion Henry problem in a fractured aquifer with the DFN and the equivalent porous medium (EPM) models. They showed that, the EPM model yields a good estimate of the steady-state flow field but is unable to reproduce the salt distribution of the DFN model, irrespective of fracture density. In this work, we consider the diffusive Henry saltwater intrusion problem for which a semi-analytical solution exists. High and low diffusion cases are investigated with the new MFE-DG model and the results with different fracture densities are compared to the semi-analytical solutions. We consider a fractured aquifer with an orthogonal fracture network as in[START_REF] Koohbor | Uncertainty analysis for seawater intrusion in fractured coastal aquifers: Effects of fracture location, aperture, density and hydrodynamic parameters[END_REF][START_REF] Sebben | Seawater intrusion in fractured coastal aquifers: A preliminary numerical investigation using a fractured Henry problem[END_REF]. The fracture network is formed by   1 f n  horizontal fractures of 2m length crossing   21 f n vertical fractures of 1m length (Figure5).

Figure 5 :

 5 Figure 5: Conceptual model of the fractured Henry problem with an orthogonal fracture network   5

Figure 6 :

 6 Figure 6: Steady-state solution for the original (a) and the low-dispersion (b) Henry problems.

Figure 7 :

 7 Figure 7: Comparison between semi-analytical and numerical results: (a) horizontal (at y=0.2m) and (b) vertical (at x=1.5m) concentration profiles for the standard Henry problem,

Figure 8 :

 8 Figure 8: Steady state solution for the original (a) and the low-dispersion (b) Henry problems. Flow and transport processes occur only in the fracture network with 20 nf  .

Figure 9 :

 9 Figure 9: Comparison between semi-analytical and numerical results: (a) horizontal (at y=0.2m) and (b) vertical (at x=1.5m) concentration profiles for the standard Henry problem, (c) horizontal (at y=0.2m) and (d) vertical (at x=1.5m) concentration profiles for the lowdispersion Henry problem. Flow and transport processes occur only in the fracture network corresponding to 20 nf  .

Figure 10 :. 4 . 1

 1041 Figure 10: Comparison between semi-analytical and numerical results: (a) horizontal (at y=0.2m) and (b) vertical (at x=1.5m) concentration profiles for the standard Henry problem, (c) horizontal (at y=0.2m) and (d) vertical (at x=1.5m) concentration profiles for the lowdispersion Henry problem. Flow and transport processes occur only in the fracture network corresponding to 50 nf  .

Figure 11 :Table 2 .

 112 Figure 11: (a) description of the flow and transport problem in fractured porous media with

  Figure 12.

Figure 12 .

 12 Figure 12. Breakthrough curves for the tracer simulation with the MFE-DG model using a

Figure 14 :

 14 Figure 14: Breakthrough curves for the tracer simulation with MFE-DG and MFE-FV models

Figure 15 :

 15 Figure 15: Breakthrough curves for transport with density and viscosity contrasts with MFE-

  this work, a new model is developed for density and viscosity coupled fluid flow and mass transport in heterogeneous fractured porous media. The model is based on the DFM approach and uses an explicit description of the fracture network. Fractures are considered infilled by a porous medium having different properties from those of the porous matrix. Both advective and dispersive processes are considered in the matrix and the fracture continua. Cell-based numerical methods which ensure local mass conservation are employed. Matrix-fracture transfers are implicitly taken into account without any transfer function. The cross-flow equilibrium is assumed only across the fractures which avoids constraints on the size of the matrix grid-cells next to fractures. Continuity of variables (pressure and concentration) as well as fluxes (fluid flux, advective and dispersive contaminant fluxes) are ensured at the matrixfracture interface as well as at the intersection of several fractures.

  Three test problems are investigated to assess the validity and efficiency of the new MFE-DG model. The first test problem deals with linear flow and transport in a 2D domain containing a ''+''-shaped barrier/fracture network. This problem is simulated with the new 1D-2D model as well as with a 2D-2D model with a very fine mesh to correctly take into account flow and transport processes in the fracture/barrier network. The results of the two models are in very good agreement for pressure and concentration distributions. The new 1D-2D model is around 17 times more efficient than the 2D-2D model. The second test problem deals with the Henry saltwater intrusion problem. Two configurations are investigated. The first configuration deals with an unfractured domain. The results for this configuration show a very good agreement between numerical and semi-analytical breakthrough curves for the original Henry problem as well as for the challenging case of saltwater intrusion with a narrow transition zone involving 20 times less diffusion. In the second configuration, saltwater intrusion occurs only through the fracture network. The results, in the case of a dense fracture network, show a very good agreement between the numerical and semi-analytical solutions for both high and lowdispersion problems. These results validated the MFE-DG numerical model for density-driven flow in the porous matrix as well as in the fracture network.The third test problem deals with stable flow involving high density and viscosity contrasts in a fractured heterogeneous porous medium containing two vertical fractures of different permeabilities. The test case is simulated with both MFE-DG and standard FE models. When used on the same coarse mesh, the results of the two models are different. Further, contrarily to MFE-DG, FE yields inaccurate results with strong unphysical oscillations. When FE is used on a finer mesh of around 3.5 times more elements, it yields a concentration profile in

It is worth noticing that the good results obtained by the MFE-DG code on the coarse mesh are largely due to the high performance of the DG method for advection-dominated transport problems. To illustrate this point, the test case is simulated using the upwind FV method instead of the DG method in the fracture network. To this aim, a variant of the developed numerical model, named MFE-FV, is obtained by reducing the spatial order of the DG method from piecewise linear to constant approximation. The obtained scheme corresponds to the standard upwind FV scheme for advection in the fracture network. The breakthrough curve obtained with MFE-FV is compared to that obatained with MFE-DG in Figure 14.

good agreement with the coarse-mesh MFE-DG solution. This validates the developed MFE-DG model and points out its efficiency and accuracy for solving flow and transport with high density and viscosity variations in heterogeneous fractured porous media.