Yannick Caillabet

Pierre Fabrie

Didier Lasseux

Michel Quintard

Computation of large-scale parameters for dispersion in fissured porous medium using finite-volume method

Keywords: fissured medium, joint element, large-scale dispersion, mass tranfer coefficient, porous medium, two-equation model

come

Computation of large-scale parameters for dispersion in fissured porous medium using finite-volume methods

Introduction

Mass transport phenomena in porous media is of major concern in a wide variety of domains ranging from chemical engineering to environmental issues. Most of the time, one is confronted with heterogeneous structures for which several different scales are involved, from the pore-scale to the column or field scale. This may result, for instance, to what is often referred to as anomalous dispersion at the field scale, or, in other words, to a non-fickian response of the system [START_REF] Bacri | Transient non-Gaussian tracer dispersion in porous media[END_REF][START_REF] Matheron | Is transport in porous media always diffusive? A counter example[END_REF][START_REF] Sposito | Fundamental problems in the stochastic convection-dispersion model of solute transport in aquifers and fields soils[END_REF].

In this work, we focus our attention to a special class of heterogeneous porous media called fissured porous media. Our definition of a fissured medium is classical, i.e., this is a class of material characterized by a porous matrix which contains most of the fluid in place and is embedded in a fissure network, which represents a tiny volume fraction but has a far greater permeability [START_REF] Reiss | Reservoir Engineering en Milieu Fissuré[END_REF]. The situation is schematically depicted in figure 1 where we have represented the pore-scale with the solid σ -phase, the liquid β-phase and their characteristic length scales l σ and l β , respectively. The next scale can be identified as the Darcy scale, the characteristic lengths being the matrix block size l ω and fissure length l η . Finally, L represents the column or field macroscopic size. In this paper we are interested in a large-scale description of the flow in the reservoir, i.e., the effective behavior of the large-scale volume V ∞ represented in figure 1. The one-phase flow process in fissured porous media has received a lot of attention in the literature. A review of this problem can be found in [START_REF] Chen | Transient flow of slightly compressible fluids through double-porosity, doublepermeability systems -a state-of-the-art review[END_REF]. Mixed models have been obtained by using homogenization theory (see for instance [START_REF] Douglas | Dual porosity models for flow in naturally fractured reservoirs[END_REF]). They corresponds to a large-scale equation for the effective fissured medium coupled with specific Darcy-scale problems for the flow in the blocks. Fully large-scale models have been widely used under the form of two large-scale equations for the effective fissured medium and matrix medium. The most popular forms are those proposed by Barenblatt and Zheltov [4] or Warren and Root [START_REF] Warren | The behavior of naturally fractured reservoirs[END_REF]. A general formulation of the two-equation model for fissured media has been proposed in [START_REF] Quintard | Transport in chemically and mechanically heterogeneous porous media, I: Theoretical development of region-averaged equations for slightly compressible single-phase flow[END_REF][START_REF] Quintard | Transport in chemically and mechanically heterogeneous porous media, II: Comparison with numerical experiments for slightly compressible single-phase flow[END_REF]. All effective properties in this latter theoretical development, i.e., the effective permeabilities and the mass exchange coefficient, are given explicitly by three 'closure problems' to be solved over representative unit cells of the fissured media under consideration. Numerical solutions of these closure problems have been obtained for unit cells with simple geometries using Cartesian grids [START_REF] Quintard | Transport in chemically and mechanically heterogeneous porous media, II: Comparison with numerical experiments for slightly compressible single-phase flow[END_REF]. In the case of more complex fissured media, numerical procedures based on a volume element formulation have been proposed in [START_REF] Caillabet | Implementation of finite-volume method for the determination of equivalent parameters in fissured porous medium[END_REF].

This paper deals specifically with the transport of solutes in fissured media, or, more generally, in the kind of highly heterogeneous medium represented in figure 2. The porous medium under consideration can be represented by a two-region system, identified as the matrix (ω-region) and the fissures (η-region), respectively. Moreover, it is assumed that the medium obeys a hierarchichal set of length-scales excluding media with strongly evolving heterogeneities or fractal properties [START_REF] Cushman | On unifying the concepts of scale, instrumentation and stochastics in the development of multiphase transport theory[END_REF]. We further consider that solute transport do not affect the physical properties of the fluid such as viscosity and density, and that adsorption do not occur on the solid phase. In two-region systems, we may have two different classes of systems. If the permeability of the two regions is highly contrasted, advection occurs mainly in the more permeable region, while diffusion takes place in the less permeable one. The time-scales associated with these different mechanisms leads to local mass non-equilibrium, and this may call for a twomedium description. The length and time scales constraints for local mass equilibrium are for instance discussed in [START_REF] Quintard | Transport in chemically and mechanically heterogeneous porous media, IV: Large-scale mass equilibrium for solute transport with adsorption[END_REF]. Two-equation models have already been used extensively in the literature to model such flows [START_REF] Brusseau | Modeling the transport of solutes influenced by multiprocess nonequilibrium[END_REF][START_REF] Coats | Dead-end pore volume and dispersion in porous media[END_REF][START_REF] Smedt | Mass transfer in porous media with immobile water[END_REF][START_REF] Goltz | Three-dimensional solutions for solute transport in an infinite medium whith mobile and immobile zones[END_REF][START_REF] Gvirtzam | Mass exchange between mobile freshwater and immobile saline water in the unsaturated zone[END_REF][START_REF] Gwo | Using a multiregion model to study the effects of advective and diffusive mass transfer on local physical nonequilibrium and solute mobility in a structured soil[END_REF][START_REF] Van Genuchten | Mass transfer studies in sorbing porous media I. Analytical solutions[END_REF][START_REF] Zurmühl | Modeling transient water and solute transport in a biporous soil[END_REF]. This two-equation concept has been generalized to a second class of system in which advection may also occur in the less permeable region [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF][START_REF] Dykhuizen | A new coupling term for dual-porosity models[END_REF][START_REF] Gerke | A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media[END_REF][START_REF] Gwo | Using a multiregion model to study the effects of advective and diffusive mass transfer on local physical nonequilibrium and solute mobility in a structured soil[END_REF][START_REF] Gwo | Mass transfer in structured porous media: Embedding mesoscale structure and microscale hydrodynamics in a two-region model[END_REF][START_REF] Skopp | Miscible displacement in structured soils: Two-region model with small interaction[END_REF]. This is often the case for moderate permeability ratio. Such a two-equation model has been derived using a volume averaging technique [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF], and in this theoretical development local problems are developed to estimate the effective coefficients in the large-scale equations. These 'closure problems' have been solved for simple unit cells, like for instance stratified media. Our objective in this paper is to compute the effective coefficients on complex heterogeneous media using the local problems proposed in [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF]. The proposed numerical techniques is inspired from the work in [START_REF] Caillabet | Implementation of finite-volume method for the determination of equivalent parameters in fissured porous medium[END_REF]. It is detailed below, and results are presented in subsequent sections.

Macro-scale model and related closure problems

The major objective of this work is to present an efficient numerical procedure to be able to compute the effective large-scale transport coefficients in the two-equation model. For this reason, we start with the large-scale two-equation model obtained from the averaging procedure described in [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF]. It is beyond the scope of this paper to provide all the details in this development. Only equations necessary to understand our work are provided. The large-scale equations obtained by Ahmadi et al. [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF] are listed below:

ε η φ η ∂C * η ∂t + ∇ φ η U * η C * η + • • • = ∇ D * * ηη • ∇C * η + D * * ηω • ∇C * ω -α * C * η -C * ω , ( 1
)
ε ω φ ω ∂C * ω ∂t + ∇ φ ω U * ω C * ω + • • • = ∇ D * * ωω • ∇C * ω + D * * ωη • ∇C * η -α * C * ω -C * η , ( 2
)
where ε α is the porosity in the α-region, φ α is the volume fraction of the α-region in the large-scale averaging volume, U * α and C * α are the average filtration velocity and intrinsic average concentration in the α-region.

These equations are obtained by establishing a direct while approximated link between the Darcy-scale and the large-scale concentration fields. Following Ahmadi et al. [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF], this link is expressed as

cη = C η -C * η = b ηη • ∇C * η + b ηω • ∇C * ω -r η C * ω -C * η , (3)
cω = b ωη • ∇C * η + b ωω • ∇C * ω + r ω C * ω -C * η , ( 4
)
where b ηη , b ηω , b ωη , b ωω , r η , and r ω are called closure variables. The closure variables are solution of three closure problems that are detailed later in this paper. These closure problems depend on the Darcy-scale velocity field V η and V ω . With the assumptions made, these fields can be obtained independently from the concentration fields by solving Darcy's law over the unit cell representative of the system under consideration. The problem of the pressure field at the region-averaged scale follows:

c η ∂P β βη ∂t = ∇ 1 µ β K βη • ∇P β βη in the η-region, (5)
B.C.1 n ηω K βη • ∇P β βη = n ηω K βω • ∇P β βω at A ηω , (6)
I.C. P β βη = f η (r), P β βω = f ω (r) at t = 0. (10)
The effective coefficients in equations ( 1) and (2) will be expressed as a function of these closure variables. The large-scale dispersion tensors D * * are calculated with the formulas below [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF],

D * * ηη = D * η • φ η I + 1 V ∞ A ηω n ηω b ηη dA - 1 V ∞ V η ṽβη b ηη dV , ( 12
)
D * * ηω = D * η • 1 V ∞ A ηω n ηω b ηω dA - 1 V ∞ V η ṽβη b ηω dV , ( 13
)
D * * ωη = D * ω • 1 V ∞ A ηω n ωη b ωη dA - 1 V ∞ V ω ṽβω b ωη dV , ( 14
)
D * * ωω = D * ω • φ ω I + 1 V ∞ A ωη n ωη b ωω dA - 1 V ∞ V ω ṽβω b ωω dV . ( 15
)

Closure problems

We have to solve the three closure problems written below in order to get largescale equations. The purpose of the first problem is to find the b ηη and b ωη vectors. It is given by

Problem I. ∇(V η b ηη ) + ṽη = ∇ D * η • ∇b ηη + ∇ D * η -φ -1 η c ηη , in the η-region, (16)
∇(V ω b ωη ) = ∇ D * ω • ∇b ωη -φ -1 ω c ωη , in the ω-region. ( 17
)
Interface conditions

n ηω • D * η • ∇b ηη + n ηω • D * η = n ηω • D * ω • ∇b ωη at A ηω , ( 18
) b ηη = b ωη at A ηω . ( 19
)
Boundary conditions

b ηη (r + l i ) = b ηη (r), b ωη (r + l i ) = b ωη (r), i = 1, 2, 3. (20)
Average conditions

{b ηη } η = 0, {b ωη } ω = 0, (21)
c ηη = -c ωη = - 1 V ∞ A ηω n ηω • V βη b ηη -D * η • ∇b ηη -D * η dA. ( 22
)
The velocity deviation is expressed as the deviation between the region-averaged velocity field and the Darcy's scale velocity field

V β = {V β } + ṽβ . ( 23
)
The second problem is similar, and the unknowns are the vectors b ηω and b ωω .

Problem II.

∇(V η b ηω ) = ∇ D * η • ∇b ηω -φ -1 η c ηω , in the η-region, (24)
∇(V ω b ωω ) + ṽω = ∇ D * ω • ∇b ωω + ∇ D * ω -φ -1 ω c ωω , in the ω-region. ( 25
)
Interface conditions

n ηω • D * η • ∇b ηω = n ηω • D * ω • ∇b ωω + n ηω • D * ω at A ηω , ( 26
) b ηω = b ωω at A ηω . ( 27
)
Boundary conditions b ηω (r

+ l i ) = b ηω (r), b ωω (r + l i ) = b ωω (r), i = 1, 2, 3. (28)
Average conditions

{b ηω } η = 0, {b ωω } ω = 0, (29)
c ηω = -c ωω = 1 V ∞ A ωη n ωη • V βω b ωω -D * ω • ∇b ωω -D * ω dA. ( 30
)
The last problem involves the scalar fields r η and r ω .

Problem III.

∇(V η r η ) = ∇ D * η • ∇r η -φ -1 η α * , in the η-region, (31)
∇(V ω r ω ) = ∇ D * ω • ∇r ω -φ -1 ω α * , in the ω-region. ( 32
)
Interface conditions

n ηω • D * η • ∇r η = n ηω • D * ω • ∇r ω at A ηω , (33) r η = r ω at A ηω . ( 34
)
Boundary conditions

r η (r + l i ) = r η (r), r ω (r + l i ) = r ω (r), i = 1, 2, 3. (35)
Average conditions

{r η } η = 0, {r ω } ω = 0. ( 36
)
The mass exchange coefficient is given by

α * = - 1 V ∞ A ηω n ηω • V η r η -D * η • ∇r η dA = 1 V ∞ A ωη n ωη • V ω r ω -D * ω • ∇r ω dA. ( 37
)
The equations in these three problems appear to be an integro-differential equation. However, this problem can be overcome by following a similar procedure as those proposed in [START_REF] Quintard | Transport in ordered and disordered porous media: Volume-averaged equations, closure problems, and comparison with experiment[END_REF]. This is explained below for the first closure problem.

The first problem is decomposed as indicated in equation (38), in which the unknowns are respectively a vector and a scalar field that must obey problems Ia and Ib described below,

b ηη = b Iη + B Iη c ωη ; b ωη = b Iω + B Iω c ωη . ( 38
)
Before giving problems Ia and Ib, it must be emphasized that we may introduce similar decompositions for problems II and III as given in equations ( 39) and (40):

b ηω = b IIη + B IIη c ωη , b ωω = b IIω + B IIω c ωη , (39) r η = α * r f , r ω = α * r m -1. ( 40
)
Only problems Ia and Ib are detailed as the other problems are deduced in a similar way.

Problem Ia.

∇(V η b Iη ) + ṽη = ∇ D * η • ∇b Iη + ∇ D * η , in the η-region, (41)
∇(V ω b Iω ) = ∇ D * ω • ∇b Iω , in the ω-region. ( 42
)
Interface conditions

n ηω • D * η • ∇b Iη + n ηω • D * η = n ηω • D * ω • ∇b Iω at A ηω , ( 43
) b Iη = b Iω at A ηω . ( 44
)
Boundary conditions b Iη (r

+ l i ) = b Iη (r), ( 45
) b Iω (r + l i ) = b Iω (r), i = 1, 2, 3. (46)
Average conditions

{b Iη } η + {b Iω } ω = 0. ( 47
)
Problem Ib.

∇(V η B Iη ) = ∇ D * η ∇B Iη -φ -1 η in the η-region, ( 48
)
∇(V ω B Iω ) = ∇ D * ω ∇B Iω -φ -1 ω in the ω-region. ( 49
)
Interface conditions

n ηω • D * η • ∇B Iη = n ηω • D * ω • ∇B Iω at A ηω , ( 50
) B Iη = B Iω at A ηω . ( 51
)
Boundary conditions

B Iη (r + l i ) = B Iη (r), ( 52
) B Iω (r + l i ) = B Iω (r), i = 1, 2, 3. ( 53
)
Average conditions

{B Iη } η + {B Iω } ω = 0, (54)
c ωη = - {b Iη } η {B Iη } η . ( 55
)

Large-scale parameters

From the solutions of the three closure problems, the large-scale dispersion tensors D * * are calculated, while the mass exchange coefficient is given by equation (37).

Problems of similar, while simpler, forms have been solved in [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: Numerical results for effective properties[END_REF]. The numerical solution was obtained by a finite difference scheme over Cartesian grids. Such grids are not suitable for achieving a good accuracy, especially with fissures of any orientation. In addition, to the natural complexity of the geometry is associated a difficulty due to the small thickness of the fractured regions. The objective of this paper is to propose a numerical solution over an unstructured mesh, with a special treatment of the elements in the fracture to account for the small thickness. This is presented in the next section.

Numerical study

In this section we detail the scheme we have used to solve the closure problems presented in the previous section.

The main idea is to solve the problems on a fissured porous medium using unstructured grids fitting closely the fissure network. One example is represented in figure 3. A finite volume method is used, and in order to reduce the numerical diffusion the closure problems are slightly modified. We write below the new problem Ia,

∂b Iη ∂t + ∇(V η b Iη ) + V η = ∇(D * η • ∇b Iη ) + ∇ D * η , in the η-region, ( 56
)
∂b Iω ∂t + ∇(V ω b Iω ) = ∇(D * ω • ∇b Iω ), in the ω-region. ( 57
)
Interface conditions

n ηω • D * η • ∇b Iη + n ηω • D * η = n ηω • D * ω • ∇b Iω at A ηω , (58) b Iη = b Iω at A ηω . ( 59
)
Figure 3. Fissured porous medium.

Boundary conditions

b Iη (r

+ l i ) = b Iη (r), ( 60
) b Iω (r + l i ) = b Iω (r), i = 1, 2, 3. (61)
Average conditions

{b Iη } η + {b Iω } ω = 0. ( 62
) Initial conditions b Iη (t = 0, r) = f (r); b Iω (t = 0, r) = g(r). ( 63
)
We have now an evolution problem, and, for long times, the asymptotic solution is the same as for original problem Ia. We notice that both diffusive and convective terms appears in the equations. We choose to use an explicit scheme to express the convection, and an implicit one to model the diffusion.

Another difficulty to overcome is the fact that the dispersion tensor D * ω may be anisotropic. In order to take care of the anisotropy which appears in the diffusive term, a 10-point finite volume scheme is built. It is a cell center scheme which provides a symmetric consistent approximation of the flux on each cell. In order to solve for the convection part, we use a scheme introduced by Eymard et al. [START_REF] Eymard | Finite Volume Methods[END_REF]. We need to know V β on the domain to be able to solve the problem. First, we have to calculate the pressure field given in equations [START_REF] Brusseau | Modeling the transport of solutes influenced by multiprocess nonequilibrium[END_REF]. However, we only have in this pressure equation a diffusive term similar with the one for the dispersion if we replace K β with D * . Therefore, we will use the same scheme for the pressure equation and we will only detail below the scheme for the dispersion problem. We also note that the Darcy-scale permeability tensor K βω may be anisotropic.

We have a fissured porous medium with a high contrast of transport properties between the matrix and the fissures. The porous matrix is meshed using triangular cells. In order to avoid the need for a very fine mesh in the fissure network, we adopt the joint element strategy which proved to be very efficient in solving the one-phase flow case. This is detailed in [START_REF] Caillabet | Implementation of finite-volume method for the determination of equivalent parameters in fissured porous medium[END_REF], and here we only recall what is a joint element.

A joint element is a small rectangle of thickness ε corresponding to the fissure thickness q and of length δx which is the average step of the triangular mesh in the porous matrix (figure 4). Some intersection nodes are added when several fissures cross or when the fissure turns.

On this mesh, we build a first order conservative finite volume scheme. On the cell equation (56) takes the form,

V ( ) b n+1 I -b n I δt + ∇ V β b n I + ṽβη d = ∇ D * • ∇b n+1 I + ∇ D * η d , ( 64
)
where δt is the time step. Using Green's formula, the equation becomes V ( )

b n+1 I -b n I δt + ∂ n• V β b n I d + ṽβη d = ∂ n• D * •∇b n+1 I d + ∇ D * η d .
(65) In this equation, the CFL condition has the following form,

L δt δx < 1, ( 66
)
where L is the domain length, and δx is the average mesh size. Under this condition, we model first the diffusive term, then the convection term. The treatment of the source terms is detailed later.

Diffusive term

Here, we build a first order, symmetric, conservative approximation of the flux on the cell side. The scheme is a 10 point finite volume method which takes into account the anisotropy. It is a hybrid scheme, i.e., the approximation is not the same on each cell but depends on the triangle geometry. With geometric considerations, the best scheme is applied. Unknowns will be affected in general to the medians intersection. However, this choice is not acceptable if this intersection is outside the triangle. These cases require a different treatment as indicated below.

Starting from equation (65), a first order approximation of the diffusion is provided by the following equation,

σ n • D * ω • ∇b Iω dσ = m(σ )n • D * ω • ∇b I ω (x) , ( 67
)
where σ is a side of a triangular porous cell, and x is the midpoint of this side. The scheme is cell center, thus we have to express the value at x using the values at the centers of the neighbouring cells.

The dispersion tensors may be anisotropic, so n • D * ω is generally not in the same direction as n which is the outwardly oriented normal from a side of a triangular porous matrix cell. We write then

n • D * ω = n • n • D * ω n + n • n • D * ω n . ( 68
)
Here, (n, n ) is an orthonormal base. Thus we split the interfacial flux in its normal and transversal components.

The normal component can be expressed in a similar way as in the isotropic case which is detailed in [START_REF] Caillabet | Implementation of finite-volume method for the determination of equivalent parameters in fissured porous medium[END_REF]. The transversal flux needs two points called γ 2 and γ 1 on the interface to get a first order discretization given below,

n • n • D * ω n • ∇b Iω (x 1 ) = n • n • D * ω b Iω (γ 2 ) -b Iω (γ 1 ) d(γ 2 , γ 1 ) . ( 69
)
Several geometric cases are described below. In these examples common notations have been done. The w i are the centers of the cells, and w is the center of the calculating cell.

The n i vectors are the normals of the sides of the porous cells.

First case: matrix-matrix interface without any neighbouring fissures

A typical situation is described in figure 5. The idea is to judiciously place the points γ 1 and γ 2 in order to eliminate them using the others unknowns. γ 2 is the midpoint of α 2 and β 2 , where α 2 and β 2 are respectively the intersection between the interface and n 2 and the intersection between the interface and n 5 . So α 2 , w, and w 2 are on the same

(α 2 ) is given below, b Iω (α 2 ) = d(w 2 , α 2 )b Iω (w) -d(w, α 2 )b Iω (w 2 ) d(w, w 2 ) , ( 70
) b Iω (γ 2 ) is given by b Iω (γ 2 ) = b Iω (α 2 ) + b Iω (β 2 ) 2 . ( 71
)
We have a similar expression for b Iω (γ 1 ), and the transversal flux is expressed on the interface. The resulting flux on the interface is given by

n • D * ω • ∇b I ω (x) = (n • n • D * ω ) b Iω (w 1 ) -b Iω (w) d(w, w 1 ) + n • n • D * ω b Iω (γ 2 ) -b Iω (γ 1 ) d(γ 2 , γ 1 ) . ( 72
)
Second case: matrix-matrix interface with a neighbouring fissured interface We suppose that there is a fissure on the interface 2 (figure 6). Thus we express b Iω (α 2 ) as a function of b Iω (w) and b Iη (F 2 ) and a source term appears. The gradient of the solution jumps at the fissure-matrix interface (equation ( 58)) (in x 2 on figure 6). We have the equations below,

b Iω (α 2 ) = b Iω (x 2 ) -d(x 2 , α 2 )n 2 • ∇b Iω (x 2 ). ( 73
)
Figure 6. Approximation of the flux with a neighbouring fissured interface.

We express the value of the solution at x 2 using the jump of the gradient

b Iω (w) = b Iω (x 2 ) -d(w, x 2 )n 2 • ∇b Iω (x 2 ), ( 74
) b Iη (F 2 ) = b Iη (x 2 ) + ε 2 n 2 • ∇b Iη (x 2 ), ( 75
)
where F 2 is the fissured joint element on side 2 of the triangular cell.

Then we use the expression for the jump of ∇b Iη in equation ( 58). We substract the two equations ( 74) and (75) in order to eliminate b Iω (x 2 ) in equation (76).

n 2 • n 2 • D * ω n 2 • ∇b 0 Iω (x 2 ) = 2(n 2 • n 2 • D * η )(n 2 • n 2 • D * ω )(b Iη (F 2 ) -b Iω (w)) 2(n 2 • n 2 • D * η )d(w, x 2 ) + (n 2 • n 2 • D * ω )ε + n 2 (n 2 • n 2 • D * η )(n 2 • n 2 • D * ω )ε 2(n 2 • n 2 • D * η )d(w, x 2 ) + (n 2 • n 2 • D * ω )ε . ( 76
)
Thus we are able to express the value at α 2 using the values at w and F 2 , thus obtaining a 10 point finite volume scheme. At this point, we have expressed the fluxes in the matrix but there are some geometrical cases which remain. That is the case of the rectangle triangles where a normal may not cross the interface. And on some triangles, the distance d(γ 2 , γ 1 ) may be far greater than the side of the cell. In both cases, we use an alternative scheme described below.

Third case: Alternative scheme -case 1

We consider a rectangle triangle (figure [START_REF] Chen | Transient flow of slightly compressible fluids through double-porosity, doublepermeability systems -a state-of-the-art review[END_REF]). In this case the normal n 2 does not cross the interface. Thus, a point on the interface must be chosen, and the value has to be expressed with the values of the neighboring centers. Here we say that α 2 = x. We can express the value at x using the values at w and w 1 .

Thus we have also a 10-point finite volume method, and we build a first order approximation of the transversal flux.

Fourth case: Alternative scheme -case 2

We consider a rectangle triangle (figure 8). Here the normals n 2 and n 3 cross at the same point on the interface, and this point is w.

In order to get a 10-point finite volume scheme α 2 , and α 1 are set at the nodes of the interface. We call G 2 the intersection of n 1 and (α 2 , x 2 ), so we obtain the expression below,

b Iω (α 2 ) = d(x 2 , α 2 )b Iω (G 2 ) -d(G 2 , α 2 )b Iω (x 2 ) d(G 2 , x 2 ) . ( 77
)
Thus the value at G 2 can be expressed using the values at w and w 1 with the equation below,

b Iω (G 2 ) = d(w 1 , G 2 )b Iω (w) -d(w, G 2 )b Iω (w 1 ) d(w, w 1 ) . ( 78
)
We have a similar expression for x Note that this scheme is also used when the distance between γ 2 and γ 1 is far greater than the side lenght of the cell. This case ends our study in the matrix cells, and we present our strategy in the fissured joint elements in the next cases.

Fifth case: Fissured joint element

This case is the same already studied in the isotropic case because the dispersion tensor in the fissure D * η is isotropic. So, these results can also be found in [START_REF] Caillabet | Implementation of finite-volume method for the determination of equivalent parameters in fissured porous medium[END_REF]. Here we study the approximation of the flux in a fissured zone (figure 9). Once again, we use the conservation law of the flux in a joint element. In this element, we have some exchange with the neighbouring fissures and the porous matrix. These two cases are detailed.

First, the flux in the fissure network (in x 3 and x 4 ) can be modeled thanks to first order Taylor's formulas because there is no source term on these boundaries. On boundary x 3 we have the equation

D * η n 3 • ∇b Iη (x 3 ) = D * η d(F, F 1 ) b Iη (F 1 ) -b Iη (F ) . ( 80
)
Next, we study the interface between the joint elements and the porous matrix. A source term appears at the interface (at x 1 for our example). With first order Taylor's formulas, we build the approximation of the flux in equations ( 81 Using the jump of ∇b Iη in equation ( 58), we eliminate b Iη (x 1 ) from equations ( 81) and ( 82), and we get the expression for D * η n 1 • ∇b Iη (x 1 ) in equation ( 83). This is exactly the same strategy already used for a fissured boundary in the porous matrix:

n 1 • n 1 • D * η n 1 • ∇b Iη (x 1 ) = 2(n 1 • n 1 • D * η )(n 1 • n 1 • D * ω )(b Iω (w 1 ) -b Iη (F )) 2(n 1 • n 1 • D * η )d(x 1 , w 1 ) + (n 1 • n 1 • D * ω )ε -n 1 2(n 1 • n 1 • D * η ) 2 d(x 1 , w 1 ) 2(n 1 • n 1 • D * η )d(x 1 , w 1 ) + (n 1 • n 1 • D * ω )ε . ( 83
)
Using this method we obtain an approximation of the flux in the most common case. However, we must be careful with the joint elements because we have some special cases to handle. These cases occur when a redirection of the flux is requested in the fissure ie at the intersection nodes. These cases have been treated by Caillabet et al. [START_REF] Caillabet | Implementation of finite-volume method for the determination of equivalent parameters in fissured porous medium[END_REF].

The approximation of the diffusive term is done, and we can note that though the approximation of the flux is symmetric, the resulting matrix which is to be reversed is not symmetric. In the next section, the convection term is treated.

Convection term

We have to model ∂ n • ( v β b n I ) d on each cell. We use a finite volume scheme developped by Eymard et al. [START_REF] Eymard | Finite Volume Methods[END_REF], and we build a first order approximation of this convective term. A typical case is described in figure 10.

We have the following approximation,

σ n • v β b n I dσ = m(σ ) n • v β b n I (σ ), ( 84
)
where b n I (σ ) is given by

b n I (σ ) = b n I (w), if n • v β > 0, b n I (w 1 ), if n • v β < 0. ( 85
)
Figure 10. Approximation of the convection term.

We have built an approximation for both the convective and the diffusive terms. In the next section, we study the source term, and especially the velocity deviation.

Source terms: velocity deviation

The expression of the velocity is given by

V β = - 1 µ K β • ∇ P β β . ( 86
)
Solving the pressure equation will provide approximations of the pressure at the center of each cell. On the contrary, the velocity is known at the boundary of each cell, which is exactly what we want in the expression of the convection term. However, the velocity deviation is needed at the center of each cell and is given by

V β = {V β } + ṽβ . ( 87
)
Thus, we interpolate the velocity at the center using the values at the boundary. A typical case is shown in figure 11.

The value at y 1 is expressed using the values at x 2 and x 3 .

v β (y 1 ) = d(y 1 , x 3 ) v β (x 2 ) + d(y 1 , x 2 ) v β (x 3 ) d(x 2 , x 3 ) . ( 88
)
Thus the following equation is obtained,

v β (w) = d(y 1 , w) v β (x 1 ) + d(w, x 1 ) v β (y 1 ) d(y 1 , x 1 ) . ( 89
)
In order to avoid a privileged direction, we average the velocity on each side of the cell. Then, using equation (87), we obtain the deviation velocity field.

At this point, we have modeled all the terms, and we can start the resolution of the evolutive problem.

Calculation of the equivalent dispersion tensors

Using the different approximation, we build a system on the domain from equation (65) which takes the following form,

A - V ( ) δt b n+1 I = B - V ( ) δt b n I + Source terms, ( 90
)
where A is the matrix built from the approximation of the diffusion part on each cell, and B is the matrix built from the approximation of the convection part. Therefore, we have a system to solve at each time step, and we can note that matrix A is not symmetric. We used a preconditionned Bi-Conjugate Gradient (BiCG) to solve the system. The matrix is symmetric if both dispersion tensors are isotropic. And the time iterations are done until we reach the permanent flow. Following these ideas, we can solve problems Ia and Ib, obtaining the vector b ηη which is the solution to problem I. Similarly, we can solve problem II.

The equivalent dispersion tensors can be determined via the formulas given below:

D * * ηη = D * η • φ η I + 1 V ∞ A ηω n ηω b ηη dA - 1 V ∞ V η ṽβη b ηη dV , (91)
D * * ηω = D * η • 1 V ∞ A ηω n ηω b ηω dA - 1 V ∞ V η ṽβη b ηω dV , (92)
D * * ωη = D * ω • 1 V ∞ A ηω n ωη b ωη dA - 1 V ∞ V ω ṽβω b ωη dV , ( 93
)
D * * ωω = D * ω • φ ω I + 1 V ∞ A ωη n ωη b ωω dA - 1 V ∞ V ω ṽβω b ωω dV . ( 94
)
In these equations, the problem is to determine the value of the integral term. We know the value of the closure problems at the center of the fissures and in the neighbouring triangle in the porous matrix. However, we need the value at the interface, so, again, we use the formulas used to establish the approximation of the flux, for instance from equation (83) we have

b ηη (x) = ε(n • n • D * ω )b ωη (K) + 2d(x, K)(n • n • D * η )b ηη (F ) ε(n • n • D * ω ) + 2d(x, K)(n • n • D * η ) -n ηω εd(x, F )(n • n • D * η ) ε(n • n • D * ω ) + 2d(x, F )(n • n • D * η ) . ( 95
)
This equation is valid for the interface between the fissure and the porous matrix. However at the intersections, and at the fissure stops, we have some special cases because there is no porous triangular cell attached to the ε-length interface. A typical case is summarized in figure 12. In this case, we have a matrix-fissure interface following n r , but in this direction we also have a joint element (here F 1 ). We express b ηη (x) with the known values in F and F 1 , which gives us

b ηη (x) = b ηη (F ) - ε 2 b ηη (F 1 ) -b ηη (F ) d(F, F 1 ) . ( 96
)
Such an expression is also valid for a fissure stop. Using these approximations, we are able to determine the large-scale dispersion tensors from the closure problem solutions.

Mass transfer coefficient

This coefficient is given by the third closure problem, and especially by the average conditions Average:

{r f } η = 0, {r m } ω = 1 α * . ( 97
)
The first equation allows us to fix the constant in order to determine the closure problem solution. Then, we obtain the mass transfer coefficient. Here we have average conditions and not interface conditions, so we can simply use the values given by the closure problems without any further operations unlike the case of the large-scale dispersion tensors.

This concludes the presentation of the numerical models. In the next section we give some results and a comparison with other available methods.

Results

In this part, we present some results on several fissured porous media. We first test the method on classical examples like the stratified or 'sugarbox' medium. Later on, more complex cases are presented with inclined fissures and with anisotropic dispersion tensors.

Stratified medium

We consider a stratified porous medium (figure 13). On this medium, we know an analytical solution of the closure problems. We compare it with the numerical results in order to validate the method. The analytical results are given in [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF], and take the following form:

D * * ηη yy = φ η D * η yy , (98)
D * * ηω yy = D * * ωη yy = 0, (99)
D * * ωω yy = φ ω D * ω yy , ( 100
)
α * = 12 (l η + l ω ) 2 (D * η ) xx (D * ω ) xx φ ω (D * η ) xx + φ η (D * ω ) xx . ( 101
)
First, we solve the pressure problem on this medium, where we suppose that a pressure gradient is applied on the side of the domain. The pressure and the velocity fields are displayed in figure 14. The flow is supposed to be in the layer direction (y-direction in figure 13). Then the dispersion closure problems are solved. In this example, the time and the space steps are fixed to the following values:

δt = 0.01; δx = 0.025. ( 102
)
The CFL condition is checked, and we note that we need about 20 time iterations to reach the convergence criteria. Each time iteration costs about 300 iterations in the BiCG process. So the computation is quite fast, and the numerical results are close to the analytical ones. For example, on the α * coefficient the error is roughly 0.1%. Hence, two decimal figures are guaranteed on the numerical results, which is in agreement with the mesh step and with the fact that we have a first order scheme.

The solutions of the closure problems I and III are displayed in figure 15.

The next graph (figure 16) shows the α * results when the dispersion tensors are anisotropic. Five different tensors are studied and displayed below. These tensors will also be used in the next cases. • Case 1 (isotropic):

D * ω = 1 0 0 1 . ( 103
)
• Case 2:

D * ω = 2 0 0 1 . ( 104
)
• Case 3:

D * ω = 1 0 0 2 . ( 105
)
• Case 4:

D * ω = 1 0.5 0.5 1 . ( 106
)
Once again, the numerical results are in good agreement with the analytical ones. This ends the stratified porous medium case. We have validate our scheme on such medium. In the next example, we study the sugarbox medium where we have no analytical solution, but this medium has been often studied [START_REF] Quintard | Transport in chemically and mechanically heterogeneous porous media, II: Comparison with numerical experiments for slightly compressible single-phase flow[END_REF]).

Sugarbox medium

This porous medium has two fissure directions which cross at 90 degrees (figure 17). On this porous medium, no analytical solution is available. We can compare our results with those obtained by Ahmadi et al. [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF]

These results have been obtained after the resolution of the pressure problem, where we have supposed that a pressure gradient is applied on the sides of the domain. The pressure and the velocity fields are shown in figure 18.

In figure 19, the solutions of the problems I and III are displayed. We see that the symmetry of the medium is respected by the closure problem solutions.

Figure 20 shows the variations of the mass transfer coefficient with the anisotropy. These two first examples have been used to validate the method. The results found shows only small differences with analytical solutions and others existing method. The computation time remains small too. In the next examples, we study more complex cases, and we fully use the advantages of the unstructured grids by studying fissured network with inclined fissures.

Nontrivial case: example 1

We consider a porous medium with a complex fissured network (figure 21). On this fissured medium, we first solve the pressure problem. The pressure and the velocity fields are displayed in figure 22.

With these values, we can start the resolution of the closure problems. The CFL condition is satisfied by using the following time and space steps:

δt = 0.01; δx = 0.025. ( 112
)
We solve the closure problems for different D * ω either isotropic or anisotropic. The mass transfer coefficient variations with anisotropy is displayed in figure 23.

Nontrivial case: example 2

We consider a porous medium with a complex fissured network (figure 24). On this fissured medium, we first solve the pressure problem. The pressure and the velocity fields are displayed in figure 25.

We solve the closure problems for different D * ω either isotropic or anisotropic. The mass transfer coefficient variations with anisotropy is displayed in figure 26.

Of course, we do not have results available to check the validity of our results. However, they illustrate the ability of the proposed numerical scheme to solve for very complex fissured networks.

Conclusion

In this paper we have presented a volume element numerical model to calculate the effective properties in the two-equation models proposed to model dispersion in fissured media. In the last part, several test cases have been treated. The method has been validated on the stratified medium case for which an analytical solution is available. The difference between the numerical and analytical solution is less than 1% for a typical problem. The proposed method is fast and accurate. The stratified medium has also provided a good test in order to determine how the method handles anisotropy effects. The obtained results show good agreement with the analytical solution. Here we must note that the numerical model size remains small because we do not finely mesh the fissures. Hence, the joint element method gives good results for a small computational cost.

Once validated, the method has been tested on more complex geometries. First, the classical sugarbox medium has been delt with. Our results are close to those found by Ahmadi et al. [START_REF] Ahmadi | Transport in chemically and mechanically heterogeneous porous media, V: Two-equation model for solute transport with adsorption[END_REF] using structured grids. Then the application of the proposed method has been illustrated on more complex fissure networks.

Two last cases show the advantages of using unstructured grids, because the triangular mesh is able to handle inclined fissures. Moreover, the method is tested for anisotropic tensors, and shows a good behaviour.

Nomenclature

A ηω = A ωη : area of the interface between the η and ω-region contained in the averaging volume V ∞ (m 2 ). g: gravitational acceleration (m 2 s -1 ). I: unit tensor. i, j, k: unit base vectors parallel to the x, y, z coordinate system. K η : dimensionless equilibrium coefficient for the η-region. K ω : dimensionless equilibrium coefficient for the ω-region. K βη : Darcy-scale permeability tensor in the η-region (m 2 ). K βω : Darcy-scale permeability tensor in the ω-region (m 2 ). l η : characteristic length for the η-region (m). l ω : characteristic length for the ω-region (m). l: length of a unit cell (m). n ηω = -n ωη : outwardly directed unit normal vector pointing from the η-region toward the ω-region. P β βη : intrinsic regional pressure for the η-region (Pa). P β βω : intrinsic regional pressure for the ω-region (Pa). r: position vector (m). r η : scalar field that maps C * ω -C * η onto cη . r ω : scalar field that maps C * η -C * ω onto cω . t: time, s. V βη : Darcy-scale, superficial average velocity in the η-region, m s -1 . {V βη } η : intrinsic regional average velocity in the η-region, m s -1 . ṽβη = V βη -{V βη } η : η-region spatial deviation velocity, m s -1 . V βω : Darcy-scale, superficial average velocity in the ω-region, m s -1 . {V βω } ω : intrinsic regional average velocity in the ω-region, m s -1 . ṽβω = V βω -{V βω } ω : ω-region spatial deviation velocity, m s -1 . V ∞ : large-scale averaging volume (m 3 ). V η : volume of the η-region contained within V ∞ (m 3 ).

V ω : volume of the ω-region contained within V ∞ (m 3 ). α * : mass exchange coefficient for the η-ω system, s -1 . ε: average thickness of the fissures (m). ε η : total porosity for the η-region. ε ω : total porosity for the ω-region. φ η : volume fraction of the η-region. φ ω : volume fraction of the ω-region. µ β : shear coefficient of viscosity (Ns m -2 ).

Figure 1 .

 1 Figure 1. Different observation scales.

Figure 2 .

 2 Figure 2. Porous medium description at the region-averaged scale.

Figure 4 .

 4 Figure 4. Joint elements in the fissure network.

Figure 5 .

 5 Figure 5. Approximation of the flux for a 10-point finite volume.

Figure 7 .

 7 Figure 7. Approximation of the flux for a rectangle triangle: case 1.

2 ,b

 2 Iω (x 2 ) = d(w 2 , x 2 )b Iω (w) + d(w, x 2 )b Iω (w 2 ) d(w, w 2 ) . (79) So we express b Iω (α 2 ) using b Iω (w), b Iω (w 1 ), and b Iω (w 2 ). And we still have a 10-point finite volume scheme.

Figure 8 .

 8 Figure 8. Approximation of the flux for a rectangle triangle: case 2.

Figure 9 .

 9 Figure 9. Typical fissured joint element with its neighbours.

Figure 11 .

 11 Figure 11. Approximation of the velocity.

Figure 12 .

 12 Figure 12. Estimation of the value at an interface.

Figure 13 .

 13 Figure 13. Stratified medium meshed.

Figure 14 .

 14 Figure 14. Pressure and velocity fields.

Figure 15 .

 15 Figure 15. Solutions of the closure problems I and III.

Figure 16 .

 16 Figure 16. Variations of α * with dispersion ratio and for various anisotropic properties.

1 .Figure 17 .

 117 Figure 17. Sugarbox medium.

Figure 18 .

 18 Figure 18. Pressure and velocity fields.

Figure 19 .

 19 Figure 19. Solution of the closure problems I and III.

Figure 20 .

 20 Figure 20. Variations of α * with dispersion ratio and for various anisotropic properties.

Figure 21 .

 21 Figure 21. Area meshed with 2122 triangles.

Figure 22 .

 22 Figure 22. Pressure and velocity fields.

Figure 23 .

 23 Figure 23. Variations of α * with dispersion ratio and for various anisotropic properties.

Figure 24 .

 24 Figure 24. Nontrivial case: example 2.

Figure 25 .

 25 Figure 25. Pressure and velocity fields.

Figure 26 .

 26 Figure 26. Variations of α * with dispersion ratio and for various anisotropic properties.

  b ηη : vector field that maps ∇C * η onto C η -C * η in the two-equation model (m). b ηω : vector field that maps ∇C * ω onto C η -C * η in the two-equation model (m). b ωη : vector field that maps ∇C * η onto C ω -C * ω in the two-equation model (m). b ωω : vector field that maps ∇C * ω onto C ω -C * ω in the two-equation model (m). b Iη : vector field used in the decomposion of b ηη (m). b Iω : vector field used in the decomposion of b ωη (m). b IIη : vector field used in the decomposion of b ηω (m). b IIω : vector field used in the decomposion of b ωω (m). B Iη : scalar field used in the decomposion of b ηη . B Iω : scalar field used in the decomposion of b ωη . B IIη : scalar field used in the decomposion of b ηω . B IIω : scalar field used in the decomposion of b ηη . c η : total compressibility of the η-region, Pa -1 . c ω : total compressibility of the ω-region, Pa -1 . C η : Darcy-scale intrinsic average concentration for the β σ system in the η-region, moles/m 3 . C ω : Darcy-scale intrinsic average concentration for the β σ system in the ω-region, moles/m 3 . C η : η-region superficial average concentration, moles/m 3 . C * η = φ -1 η C η : η-region intrinsic average concentration, moles/m 3 . cη = C η -C * η : spatial deviation concentration for the η-region, moles/m 3 . C ω : ω-region superficial average concentration, moles/m 3 . C * ω = φ -1 ω C ω : ω-region intrinsic average concentration, moles/m 3 . cω = C ω -C * ω : spatial deviation concentration for the ω-region, moles/m 3 . D * η : dispersion tensor for the β-σ system in the η-region, m 2 /s. D * ω : dispersion tensor for the β-σ system in the ω-region, m 2 /s. D * * ηη : dominant dispersion tensor for the η-region transport equation, m 2 /s. D * * ηω : coupling dispersion tensor for the η-region transport equation, m 2 /s. D * * ωω : dominant dispersion tensor for the ω-region transport equation, m 2 /s. D * * ωη : coupling dispersion tensor for the ω-region transport equation, m 2 /s. D 1 , D 2 , D 3 : lengths of the sides of a porous triangular cell. K 1 , K 2 , K 3 , K j , K: porous triangular cells used to model the flux.

Acknowledgement

Financial support for this study was provided by Laboratoire Energétique et Phénomènes de Transfert.