Jacques Bodin
email: jacques.bodin@univ-poitiers.fr

Gilles Porel

Fred Delay

Fabrice Ubertosi

Stephane Bernard

Jean-Raynald De

Stéphane Bernard

Jean-Raynald De Dreuzy

Simulation and analysis of solute transport in 2D fracture/pipe networks: The SOLFRAC program

Keywords: numerical modelling, solute transport, fractured rocks

Introduction

In many geologic formations, fractures represent preferential pathways along which a dissolved contaminant can migrate rapidly. The understanding and quantification of solute transport in fractured rocks is therefore of considerable practical importance in terms of aquifer preservation. There are three main approaches to simulating solute transport in fractured rocks: discrete network simulations [START_REF] Andersson | Predicting mass transport in discrete fracture networks with the aid of geometrical field data[END_REF]Cacas et al. 1990a;[START_REF] Dverstorp | Discrete fracture network interpretation of field tracer migration in sparsely fractured rock[END_REF]Dershowitz et al. 1999;[START_REF] Huseby | Dispersion in three-dimensional fracture networks[END_REF], continuum approaches based on either equivalent-porous-medium (EPM) assumption [START_REF] Long | Porous media equivalents for networks of discontinuous fractures[END_REF][START_REF] Berkowitz | Continuum models for contaminant transport in fractured porous formations[END_REF][START_REF] Mckay | EPM modeling of a field-scale tritium tracer experiment in fractured, weathered shale[END_REF][START_REF] Gburek | Patterns of contaminant transport in a layered fractured aquifer[END_REF][START_REF] Jackson | Self-consistency of a heterogeneous continuum porous medium representation of a fractured medium[END_REF] or dual/triple-continuum approximations [START_REF] Bibby | Mass transport of solutes in dual-porosity media[END_REF]Huyakorn et al. 1983aHuyakorn et al. , 1983b;;[START_REF] Gerke | Macroscopic representation of structural geometry for simulating water and solute movement in dual-porosity media[END_REF][START_REF] Leo | A time-stepping finite element method for analysis of contaminant transport in fractured porous media[END_REF][START_REF] Kischinhevsky | Modelling and numerical simulations of contaminant transport in naturally fractured porous media[END_REF][START_REF] Arnold | Effective-porosity and dual-porosity approaches to solute transport in the saturated zone at Yucca Mountain: implications for repository performance assessment[END_REF][START_REF] Lichtner | Critique of dual continuum formulations of multicomponent reactive transport in fractured porous media[END_REF][START_REF] Hassan | On using particle tracking methods to simulate transport in single-continuum and dual continua porous media[END_REF][START_REF] Wu | A triple-continuum approach for modeling flow and transport processes in fractured rock[END_REF], and hybrid models that combine elements of both discrete network simulations and continuum approaches (Schwartz insu-00260744, version 1 -5 Mar 2008 and[START_REF] Schwartz | A continuum approach for modeling mass transport in fractured media[END_REF][START_REF] Grandi | Transport of radionuclides in isolated fractures in crystalline rocks[END_REF]National Research Council 1996;[START_REF] Tsang | Tracer transport in a stochastic continuum model of fractured media[END_REF][START_REF] Stafford | Hydrology influence of fracture truncation on dispersion: A dual permeability model[END_REF][START_REF] Lee | Hierachical modeling of flow in naturally fractured formations with multiple length scales[END_REF][START_REF] Ohman | A regional-scale particle-tracking method for nonstationary fractured media[END_REF].

As emphasized by [START_REF] Long | Porous media equivalents for networks of discontinuous fractures[END_REF] and [START_REF] Berkowitz | Continuum models for contaminant transport in fractured porous formations[END_REF], the validity of continuum approaches for fracture networks depends on several parameters, the most important being the existence of a Representative Elementary Volume (REV), i.e. a volume above which the flow and transport properties of the fractured rock can be considered as statistically homogeneous.

This assumption, however, does not hold for fractured rocks in general. Several authors have shown that most natural fracture networks exhibit fractal or scale-dependent properties [START_REF] Neuman | Universal scaling of hydraulic conductivities and dispersivities in geologic media[END_REF][START_REF] Velde | Fractal patterns of fractures in granites[END_REF][START_REF] Vignes-Adler | Fractal analysis of fracturing in two african regions, from satellite imagery to ground scale[END_REF][START_REF] Guimera | A discussion of scale effects on hydraulic conductivity at a granitic site (El Berrocal, Spain)[END_REF][START_REF] Cowie | Scaling laws for fault and fracture populations -Analyses and applications[END_REF][START_REF] Bour | Connectivity of random fault networks following a power law fault length distribution[END_REF][START_REF] Bodin | L'analyse d'images appliquée au traitement automatique de champs de fractures. Propriétés géométriques et lois d'échelles. [Use of image analysis techniques for the automatic processing of fracture maps. Geometrical properties and scaling laws[END_REF]. In that case, the homogenization scale of the network is difficult, or perhaps even impossible, to identify.

In the discrete framework, individual fractures are explicitly incorporated into the spatial domain. For 2D analyses, the fractures are represented by 1D-pipe elements [START_REF] Schwartz | A stochastic analysis of macroscopic dispersion in fractured media[END_REF][START_REF] Smith | An analysis of the influence of fracture geometry on mass transport in fractured media[END_REF][START_REF] Robinson | A laboratory and numerical investigation of solute transport in discontinuous fracture systems[END_REF]de Dreuzy et al. 2001ade Dreuzy et al. , 2001b[START_REF] De Dreuzy | Hydraulic properties of two-dimensional random fracture networks following power law distributions of length and aperture[END_REF][START_REF] De Dreuzy | Influence of spatial correlation of fracture centers on the permeability of two-dimensional fracture networks following a power law length distribution[END_REF]. For 3D problems, the fractures are represented by 2D-planar elements, which can be either disk-shaped [START_REF] Andersson | Conditional simulations of fluid flow in threedimensional networks of discrete fractures[END_REF][START_REF] Billaux | Three-dimensional statistical modeling of a fractured rock mass. An example from the Fanay-Augères Mine[END_REF]), elliptical [START_REF] De Dreuzy | Percolation parameter and percolationthreshold estimates for three-dimensional random ellipses with widely scattered distributions of eccentricity and size[END_REF], or polygonal [START_REF] Huseby | Dispersion in three-dimensional fracture networks[END_REF]. A good review on the art of simulating realistic 3D fracture networks can be found in National Research Council (1996). Once a 3D fracture network has been generated, the flow and transport calculations are performed on either 2D numerical grids constructed in each fracture plane (Nordqvist et al. 1992;[START_REF] Therrien | Three-dimensional analysis of variably-saturated flow and solute transport in discretely fractured porous media[END_REF][START_REF] Huseby | Dispersion in three-dimensional fracture networks[END_REF][START_REF] Reichenberger | A mixed-dimensional finite volume method for two-phase flow in fractured porous media[END_REF], or in a network of interconnected 1D pipes derived from the fracture network (Cacas et al. 1990a;Cacas et al. 1990b;[START_REF] Dverstorp | Discrete fracture network interpretation of field tracer migration in sparsely fractured rock[END_REF][START_REF] Moreno | Fluid flow and solute transport in a network of channels[END_REF]Dershowitz andinsu-00260744, version 1 -5 Mar 2008 Fidelibus 1999;[START_REF] Jourde | Flow behavior in a dual fracture network[END_REF]. These two alternatives are illustrated in Fig. 1.

Discretising fractures into regular or irregular grids is the more rigorous approach since aperture variations in each fracture plane are accounted for. Unfortunately, this is also very computationally demanding since flow and transport equations have to be solved for each grid element. Today's computer capabilities are far from sufficient to enable field-scale analyses using this method since a number of fractures of 10 5 -10 7 is often required to correctly depict the medium.

Modelling flow and transport in pipe networks is more computationally efficient as two-or three-dimensional problems are reduced to a series of one-dimensional ones. A physical justification of this conceptualization is that fluid flow through rock fractures is often concentrated along preferential pathways, corresponding to the paths with the lowest hydraulic resistance [START_REF] Tsang | Flow channeling in heterogeneous fractured rocks[END_REF]Bodin et al. 2003a). Of course, the use of the pipe network approach needs to address the problem of parameterisation. The number, location, and aperture distribution of flow channels have to be specified. For a real case problem, it is not clear whether such parameters may be accessed by current field investigation methods. Note that a similar problem exists with the gridding approach, because in-situ characterization of aperture distributions within fractures planes is virtually impossible. Also note that instead of considering gridded fracture models and pipe network models as rivals, these two approaches may be viewed as complementary. One could compute flow in a few gridded fractures in order to get statistics about the geometry of preferential pathways, and then simulate channel networks on a greater scale.

The aim of this paper is to present a program called SOLFRAC, which enables to simulate and analyze solute transport in complex fracture networks using a pipe network i.e. pipe networks interconnected in 2D space. Note that all the modelling concepts and numerical methods presented in this work could be easily adapted for simulating solute transport in pipe networks interconnected in 3D space (i.e. 3D fracture networks). In that case, most efforts should be put on the network simulator and not on the transport algorithms that are described in the present paper.

Input data files

The 2D fracture networks handled by the program are delimited by square domains of L x L in size. Fluid flow is established between the top and bottom faces by assigning constant values of hydraulic head to all fractures intersecting those faces. Impermeable boundaries are assigned to the remaining two faces. Each fracture is modelled as a rectangular pipe with constant aperture and width, the latter parameters possibly differing from one fracture to another. In such networks, the computation of hydraulic head at each node is performed by assuming that the flow in each bond obeys Darcy's law and by applying the principle of conservation of fluid mass at each node (i.e. Kirchhoff's law). This leads to a linear system of algebraic equations, which can be solved by either direct or iterative methods, see e.g. [START_REF] Gill | Numerical Linear Algebra and Optimization[END_REF].

The input data files required by the program are in text format. These files include the information describing the geometry of the fracture networks, the hydraulic conductivity of each fracture, and the hydraulic head values at each node. The file format used is very similar to that of the output files of the code developed by de Dreuzy et al. (2001a, 2001b[START_REF] Bernard | Le logiciel MODFRAC : simulation de réseaux de fractures 2D et de l'écoulement associé[END_REF], insu-00260744, version 1 -5 Mar 20082004) for analyzing the hydraulic properties of 2D random fracture networks following a power law length distribution. Further details about the structure of input files are given at the URL http://labo.univ-poitiers.fr/hydrasa/intranet/telechargement.htm. A program enabling to generate input files for SOLFRAC is also freely downloadable at this URL. This program, called MODFRAC, has been developed by [START_REF] Bernard | Le logiciel MODFRAC : simulation de réseaux de fractures 2D et de l'écoulement associé[END_REF] and [START_REF] Ubertosi | Résolution numérique de l'écoulement transitoire dans des réseaux de fractures 2D[END_REF]. It generates 2D fracture networks according to a wide range of statistical distribution functions for fracture density, orientation, aperture, and length. The hydraulic head value at each node can be computed by either direct methods (Cholesky or Gauss Jordan) or iterative ones (conjugate gradient or Gauss-Seidel).

Transport simulations

Transport equations

The transport simulations handled by the program are limited to scenarios for single-phase, isothermal flow conditions in which the solute concentration is diluted enough to neglect density effects. Hydraulic properties are assumed to remain constant in time, i.e. it is assumed that coupling between fluid pressure and rock stress is negligible and that chemical reactions (such as precipitation and/or dissolution) that could modify fracture openings do not occur.

The model accounts for the following transport processes: advection and hydrodynamic dispersion in the fractures, matrix diffusion, diffusion into stagnant zones within the fracture planes, sorption reactions onto the fracture walls and in the matrix, linear decay, and mass sharing at fracture intersections. Both the rock matrix and the stagnant zones are considered as immobile pore spaces in which the movement of tracers is due to molecular diffusion only (advection is neglected). A sketch illustrating a fracture-matrix system with stagnant zones in the fracture plane is given in Fig. 2. In the mathematical development of the transport equations, we make the following major assumptions: (1) the fluid velocity is constant insu-00260744, version 1 -5 Mar 2008

between two successive fracture intersections, (2) transverse diffusion and dispersion within each flow channel ensure complete mixing across its width, (3) dispersion in each flow channel is assumed to be Fickian or based on an expression making its first spatial derivative easy to approximate, (4) sorption reactions onto the fracture walls and in the matrix obey linear instantaneous equilibriums, (5) diffusion into stagnant zones within the fracture plane is perpendicular to the flow direction, (6) diffusion in the matrix is perpendicular to the fracture plane (diffusive fluxes parallel to the fracture plane are neglected), ( 7) the extent of solute diffusion within the rock matrix is small enough so that diffusion from adjacent fracture elements does not interact, (8) diffusion in the fracture plane and diffusion in the surrounding rock matrix are two independent processes that do not interact (i. e. diffusive fluxes between the stagnant zones and the rock matrix are neglected). With the above assumptions, the solute transport in each branch of the network (i.e. in each flowing bond between two nodes) is described by three coupled one-dimensional equations: one for the flow channel, one for the stagnant zones, and one for the porous matrix attached to the bond. The coupling is provided by the continuity of concentrations along the interface between the flow channel and the stagnant zones, and between the flow channel and the rock matrix. The transport equation in the fracture is written as: -3 ] is the solute concentration in the matrix; z [L] is the space coordinate in the matrix block, normal to the fracture plane; and R f [-] is a retardation factor accounting for the sorption of the solute onto the fracture walls. R f is defined as:

1 f f f f n f m e m f f z b f f f f y a c u c c c D D c c D t R x R x x a R y b R z λ = = ∂ ∂ ∂ ∂ ⎛ ⎞ ∂ ∂ + =- + + + ⎜ ⎟ ∂ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ (
1 f f K R b = + (2)
where K f [L] is the surface-sorption coefficient of the solute onto the fracture walls. The transport equation in the stagnant zones is written as:

2 ' ' 2 f nf m f f c c D c t R y λ ∂ ∂ + = ∂ ∂ (3)
The transport equation in the matrix is written as: Transport simulations are performed using the Time Domain Random Walk (TDRW) method, recently developed by [START_REF] Delay | Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks[END_REF] and Bodin et al. (2003c), from the original work by [START_REF] Banton | A new Time Domain Random Walk method for solute transport in 1-D heterogeneous media[END_REF]. The TDRW method is a Lagrangian method developed in the time domain. It allows for the one-step calculation of the particle residence time in each bond of the network and is thus very efficient in terms of computation costs, while preserving accuracy. The fundamentals of this method are summarized below.

2 2 m m m a c c c D t z λ ∂ ∂ + = ∂ ∂ (4)

insu-00260744, version 1 -5 Mar 2008

The TDRW method

For the sake of simplicity, the role of stagnant zones within fracture planes is neglected in this section and will be addressed later. The time needed for a solute particle to move between two fracture intersections can then be separated in two residence times: (1) residence time of the particle in the flow channel and (2) residence time of the particle in the rock matrix adjacent to the flow channel. The TDRW method allows solving these two residence times separately.

As a first step, the residence time of a particle in a bond is calculated assuming that the rock matrix is totally impervious, i.e. there is no matrix diffusion. In the Lagrangian framework, considering the change of variable:

f f f f f R u x c t x x c t c ∂ ∂ = ∂ ∂ ∂ ∂ = ∂ ∂ (6)
the transport equation in a bond can be rewritten as:

( )

f f f f f f f f f f f f f c D t u R c x D u t u R u R c x c 2 2 3 2 2 ∂ ∂ + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ - = + ∂ ∂ λ (7)
Note that expression ( 7) is an advection-dispersion equation written with the Fokker-Plank formalism. From the equivalence between this expression and the random-walk approach to transport, Bodin et al. (2003c) showed that the mean and variance of particle travel times could be easily identified. For an elementary displacement of length dx, these terms are written, respectively, as:

( ) 2 f f t f f R D dx u dx u x μ ∂ ⎛ ⎞ = + ⎜ ⎟ ∂ ⎝ ⎠ (8) ( ) 2 2 3 2 f t f f D dx R dx u σ = (9)
Because the particle motion along a bond of length L can be considered as a series of independent jumps, the means and variances are additive [START_REF] Rasmuson | Analysis of hydrodynamic dispersion in discrete fracture networks using the method of moments[END_REF]Bodin et al. insu-00260744, version 1 -5 Mar 20082003c). Thus, the mean and variance of the particle travel time distribution for a displacement of length L are written as:

2 2 0 0 L L f f f f t f f f f R D R D u d x u L d x u x u x μ ∂ ∂ ⎛ ⎞ ⎛ ⎞ = + = + ⎜ ⎟ ⎜ ⎟ ∂ ∂ ⎝ ⎠ ⎝ ⎠ ∫ ∫ (10) 2 2 3 0 2 L f t f f R D dx u σ = ∫ (11)
For Péclet numbers Pe = u f L/D f larger than 10, it can be shown that the travel time distribution is lognormal (Bodin et al. 2003c). Therefore, the stochastic calculation of travel times over a distance L is given by:

( ) ln ln ln f N t Z μ σ Δ = + (12) ( ) 2 2 ln 1 ln t t t μ σ μ μ + = (13) ( ) 2 2 2 ln 1 ln t t μ σ σ + = (14)
where μ t [T] and σ t 2 [T 2 ] are the mean and variance in ( 10) and ( 11), Δt f [T] the particle travel time for a travel distance L, Z N a random number drawn from a normal deviate, and μ ln [T] and σ ln 2 [T 2 ] the mean and variance of the log transform. Note that the TDRW method enables the scale-dependent dispersion coefficient to be dealt with, provided that the spatial derivative in (10) is calculable. If indexes n and n+1 refer to the upstream and downstream nodes of a bond of length L, the particle travel time in this bond can be written as:

( )

1 l n l n exp n n f N t t t Z μ σ + -= Δ = + (15)
In fracture networks, Péclet numbers can be locally less than 10 in very short bonds or in bonds with very low flow velocities. The assumption of a lognormal travel time distribution in these bonds can be flawed, and yield inaccurate results. Bodin et al. (2003c) propose an empirical correction of expression ( 13) to preserve accuracy of the TDRW method for Pe < insu-00260744, version 1 -5 Mar 2008

10. This correction consists in multiplying expression (13) by a factor β =1-1/(33Pe), which leads to:

' ln ln 2 2 1 1 l n 33 1 t t t Pe μ μ β μ σ μ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ = = - ⎜ ⎟ ⎜ ⎟ + ⎝ ⎠ ⎝ ⎠ (16)
Once the residence time of the particle in a bond has been determined, the next step is to calculate the time spent by diffusion and sorption in the matrix block attached to the bond.

Initial work by [START_REF] Delay | Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks[END_REF] only applied to non-reactive solutes. The TDRW method is here re-developed to handle matrix diffusion and sorption in the matrix. A set of analytical solutions for the transport problem described by Eq. (1-5) is provided by [START_REF] Tang | Contaminant transport in fractured porous media: Analytical solution for a single fracture[END_REF] for a continuous injection of constant concentration c 0 at the inlet of the fracture. In the case of both non-decaying solute and negligible dispersion in the fracture (i.e. λ = 0 and D f = 0), the concentration at the outlet of a fracture of length L is written as:

0 f c = for 0 0 t t ≤ < (17) 0 0 0 erfc f f t c c R t t ⎛ ⎞ Ω = ⎜ ⎟ ⎜ ⎟ - ⎝ ⎠ for 0 t t ≥ ( 18
)
where

0 f f L t R u = (19) ( ) 2 m m m e K D b θ ρ + Ω = (20)
Note that t 0 [T] is the particle residence time in the fracture, for pure advection delayed by sorption reactions onto the fracture walls. Since expressions ( 17) and ( 18) are the answer to a continuous injection, they have the significance of a cumulative probability density function for residence times by advection and matrix diffusion in a bond of length L. This distribution is written as:

insu-00260744, version 1 -5 Mar 2008

0 ad F = for 0 0 t t ≤ < (21) 0 0 erfc ad f t F R t t ⎛ ⎞ Ω = ⎜ ⎟ ⎜ ⎟ - ⎝ ⎠ for 0 t t ≥ (22)
Applying the so-called rejection method [START_REF] Yamashita | Particle-Tracking technique for nuclide decay chain transport in fractured porous media[END_REF][START_REF] Moreno | Fluid flow and solute transport in a network of channels[END_REF][START_REF] Delay | Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks[END_REF], one gets a stochastic expression of the particle residence time Δt m in the matrix:

( )

2 0 0 1 01 erfc m f t t t t R U - ⎛ ⎞ Ω Δ = -= ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ( 23
)
where U 01 is a random number drawn from a uniform distribution between 0 and 1. However, because Eq. ( 18) has been developed under the assumption of negligible dispersion in the fracture, the diffusion time Δt m cannot be merely added to the advection-dispersion time Δt f (15). [START_REF] Delay | Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks[END_REF] have shown that the mean and variance of the travel time distribution in the bond should be modified to account for the interaction between advectiondispersion and matrix diffusion. The modification is performed as follows:

* * * 2 * 0 L f f t f f R D u L dx x u μ ⎛ ⎞ ∂ = + ⎜ ⎟ ⎜ ⎟ ∂ ⎝ ⎠ ∫ (24) 2 * 2 * 3 0 2 L f t f f R D dx u σ = ∫ (
* f diff f u R u = (26) ( ) * * f f D f u = if dispersion is modelled as being proportional to fluid velocity (27) (e.g. * * f f D u α = with α [L] a dispersivity constant) insu-00260744, version 1 -5 Mar 2008 * f f D D = otherwise (28)
The coefficient R diff is a retardation factor that expresses the delay stemming from matrix diffusion as compared to pure advection. It can be viewed as the ratio t ad / t 0 , t ad being a characteristic time of advection-diffusion. [START_REF] Delay | Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks[END_REF] propose to calculate t ad as follows:

( )

( ) 0 0 0 0 t B ad t ad t B ad t f d t f d τ τ τ τ τ + + = ∫ ∫ (29)
where t 0 is defined by ( 19), B [T] is an integration boundary, and f ad is the (non-cumulative)

probability density function of residence times by advection-diffusion in a bond-matrix system, which can be derived from ( 22):

( ) ( ) ( ) ( ) 2 2 0 0 3 2 2 0 0 exp ad ad f f d F t t f t dt R t t R t t π ⎛ ⎞ Ω Ω = = - ⎜ ⎟ ⎜ ⎟ - - ⎝ ⎠ (30)
By experience, the best results are given for B = σ t /2 with σ t defined by (11). Introducing (30) into (29) and using B = σ t /2 yield the following expression for the coefficient R diff :

( )

( ) 2 exp 2 1 erfc t diff f R R ξ σ ξ π ξ ⎛ ⎞ - Ω ⎜ ⎟ = + - ⎜ ⎟ ⎝ ⎠ (31) where 0 2 2 f t f t t L R u ξ σ σ Ω Ω = = (32)
One can easily check that R diff tends to 1 when ξ tends to 0 (i.e. for negligible matrix diffusion, the advection velocity in the bond remains obviously unchanged). In summary, the total residence time of a particle in a bond-matrix system undergoing an advective-dispersive motion in the fracture, coupled to both sorption onto the fracture walls and diffusion-sorption into the matrix, can be calculated as follows:

insu-00260744, version 1 -5 Mar 2008 ( ) ( ) 2 ' 0 ln ln 1 01 exp erfc fm N f t t Z R U μ σ - ⎛ ⎞ Ω ⎜ ⎟ Δ = + + ⎜ ⎟ ⎝ ⎠ (33)
where expressions ( 24) and ( 25) have to be used in place of expressions ( 10) and ( 11) in the calculus of μ ln ' ( 16) and σ ln 2 ( 14). If the solute undergoes linear decay, the solute mass associated with each particle is decreased at each node according to:

1 exp n n f m mp mp t λ + ⎡ ⎤ = -Δ ⎣ ⎦ (36)
where mp n and mp n+1 [M] represent the mass of the particle at nodes n and n+1, respectively.

Solute diffusion in stagnant zones

In the preceding section, the rock matrix was considered to be the only part of the fractured system where the solute motion is purely diffusive. However, it is well known that fluid flow in natural fractures is often highly channelled, i.e. flow occurs in a relative small portion of the fracture plane [START_REF] Tsang | Flow channeling in heterogeneous fractured rocks[END_REF]Bodin et al. 2003a). The stagnant zones in the fracture plane therefore act as an additional "non-flowing" pore space available for solute diffusion. It is assumed that (1) diffusion in the fracture plane and diffusion in the surrounding rock matrix are two independent processes, and (2) diffusion in the stagnant zones is not influenced by the finite width of the fracture plane (unlimited diffusion). Using these assumptions, the transport problem described by equations ( 1)-( 5) can be simulated as described in the previous section, while adding a residence time Δt f ' of the particle in stagnant zones to the time Δt fm (expression. 33). By analogy with the term Δt m for matrix diffusion, the expression of Δt f ' can be written as:

( ) 2 0 ' 1 ' 01 ' erfc f f t t R U - ⎛ ⎞ Ω ⎜ ⎟ Δ = ⎜ ⎟ ⎝ ⎠ (35)
where insu-00260744, version 1 -5 Mar 2008

' 2 f m R D a Ω = (36)
Theoretically, for the same reasons as those evoked in the previous section, the diffusion time Δt f ' cannot be merely added to the time Δt fm . The mean and variance of the travel time distribution in the bond should be modified in order to account for the interaction between advection-dispersion in the flow channel and diffusion into the stagnant zones. However, the coupling error is expected to be negligible because the exchange surface between the flow channel and the stagnant zones is much more limited than the contact area with the matrix blocks. The complete expression for simulating advective-dispersive transport in a 1D-flow channel, coupled to both diffusion and sorption in adjacent rock matrix and stagnant zones in the fracture plane, can thus be written as:

( ) ( ) ( ) 2 ' 0 ' l n l n 1 1 ' 01 01 ' exp erfc erfc f f m N f t t Z R U U μ σ - - ⎡ ⎤ ⎛ ⎞ Ω Ω ⎢ ⎥ ⎜ ⎟ Δ = + + + ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ (37)
Note that the above developments could be generalized to any number of independent diffusion compartments.

Dispersion models

Solute dispersion in a fracture stems from the combined effects of molecular diffusion and heterogeneity of the fluid velocity field (hydrodynamic dispersion). In a (natural) variableaperture fracture, the heterogeneity of fluid velocities develops both along the fracture plane and across the fracture aperture. The part of hydrodynamic dispersion resulting from heterogeneity along the fracture plane is classically written as (Bodin et al. 2003b): x >> u f τ c before the Taylor-Aris dispersion regime is completely established. For t < τ c , the dispersion coefficient D L2 is time-dependent and its expression is given by [START_REF] Berkowitz | Reactive solute transport in a single fracture[END_REF]:

1 L f D u α = (
( ) ( ) 2 2 2 2 2 6 2 1 2 18 exp 105 f m L n m u b n D D t t D b n π π ∞ = ⎡ ⎤ ⎛ ⎞ = - - ⎢ ⎥ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ∑ (40)
Ippolito et al. (1994) and [START_REF] Detwiler | Solute transport in variable-aperture fractures: An investigation of the relative importance of Taylor dispersion and macrodispersion[END_REF] showed both experimentally and numerically that dispersions along the fracture plane and across its variable aperture add up. Thus, the "bulk" dispersion coefficient can be written as:

1 2 f m L L D D D D = + + (
f f m m m f f f D u D x D x D x x ub ub ub π π π π ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ∂ = - + - + - ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ∂ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ (42)
insu-00260744, version 1 -5 Mar 2008

The dispersivity α in D L1 has been considered as constant for the above derivation of D f . Note however that in the program, the value of α in each bond can be fixed as a user-defined constant, or as a function of the bond length:

L β α μ = ( 43
)
where μ and β are user-defined constants. The dispersivity α may also be considered as scaledependent, as suggested by numerous laboratory and field experiments, see e.g. [START_REF] Neuman | Universal scaling of hydraulic conductivities and dispersivities in geologic media[END_REF], [START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF], [START_REF] Schulze-Makuch | Longitudinal dispersivity data and implications for scaling behavior[END_REF]. Three models of scale-dependent dispersivity are handled by the program: power law model, asymptotic model, and exponential model. These models are written, respectively, as (Pickens and Grisak 1981):

( )

x x ξ α ϖ = (44) ( ) 1 x x η α η ⎛ ⎞ = Ψ - ⎜ ⎟ + ⎝ ⎠ (45)
( ) ( )

1 exp x x α κ = Φ - - ⎡ ⎤ ⎣ ⎦ (46

Solute injection

In the program, the injection of particles may be performed either along the inlet boundary, or at any point of the backbone (i.e. the flowing part) of the fracture network. In the first case, each particle is injected into one of the inlet bonds chosen randomly according to probability density either uniform or proportional to the flow rate in the bond. In the second case, the injection point is specified by the user via a graphic interface. Both short-term injection of a finite mass and continuous injection of constant concentration are handled by the program.

Short-term injection may be either instantaneous or exponential-decaying in time, and simulated by the release time T init assigned to each particle:

0 init T = if 0 γ = (47) ( ) 01 1 log init T U γ = - if 0 γ > (48)
where γ [T -1 ] is a decay coefficient defining the release-rate of the injected mass and U 01 is a random number drawn from a uniform deviate between 0 and 1. In the case where γ > 0, the inlet solute-mass flux

F i (t) [M•T -1 ] at each injection point i obeys: ( ) ( ) exp i i F t M t γ γ = - ( 49
)
where M i [M] is the total solute-mass injected in the network at the point i:

Σ i M i = M 0 , with
M 0 the total mass of solute injected into the network. If N is the total number of particles used in the simulation, the solute-mass mp i assigned to each particle for its entry in the network through the point i is equal to M 0 / N. In the case of a continuous injection of constant concentration C 0 [M•L -3 ], the mass mp i is computed as:

insu-00260744, version 1 -5 Mar 2008 ( )

0 i i i C Q t mp N δ = (50)
where

Q i [L 3 •T -1 ]
is the flow rate at the point i in the network, and N i [-] is the number of particles entering the network through this point (Σ i N i = N).

Identification of elementary paths

Once the injection points have been located, the program allows identifying all the elementary paths of the network. An elementary path is defined as a ranked series of connected bonds that a particle may experience to move from its injection point to the outlet boundary of the network. This computation step is optional because it is not required for the transport simulations. It enables (1) to visualize the set of bonds that are "available" for the solute motion, and (2) it may be useful for a breakthrough-curve computation based on a convolution of analytical solutions (see hereafter). However, in well-connected networks, the number of such paths may be much greater than the number of flowing bonds and the computational works needed to identify the paths is very time-consuming. As an example, several tens of hours were necessary with a 1 GHz PC to identify 8 x 10 6 elementary paths on the transport problem illustrated in Fig. 3.

Mass transfer through fracture intersections

Rules for particle transition through fracture intersections must be defined according to a model of solute mass partitioning. Three particle routing methods are implemented in the program. These methods are based on either (1) the "perfect-mixing" model [START_REF] Smith | An analysis of the influence of fracture geometry on mass transport in fractured media[END_REF], or (2) the "streamtube" model [START_REF] Endo | A model for investigating mechanical transport in fracture networks[END_REF], or (3) the "diffusionalmixing" model [START_REF] Park | Analytical solutions for solute transfer characteristics at continuous fracture junctions[END_REF]. The perfect-mixing model is the simplest mixing rule and has been widely used in discrete transport simulations (see e.g. [START_REF] Smith | An analysis of the influence of fracture geometry on mass transport in fractured media[END_REF], Cacas et al. (1990a), [START_REF] Bradbury | Effects of fracture density and anisotropy on delineation of wellhead-protection areas in fractured-rock aquifers[END_REF]). In this model, molecular insu-00260744, version 1 -5 Mar 2008

diffusion is assumed to ensure the homogenization of the solute mass fluxes at each fracture junction; therefore the same concentration value is observed at the entrance through the outlet bonds and mass sharing is proportional to the relative discharge flow rates. At the opposite limit, the streamtube model assumes that the solute molecules strictly follow the streamlines from the inflow bonds to the outflow bonds, with no mixing at the intersection. In this case, mass sharing depends on the configuration of inlet and outlet fluxes at the node. The streamtube model has been used by [START_REF] Robinson | A laboratory and numerical investigation of solute transport in discontinuous fracture systems[END_REF], [START_REF] Wels | Retardation of sorbing solutes in fractured media[END_REF], [START_REF] Parney | Fluid velocity and path length in fractured media[END_REF], among others. The diffusional-mixing model, recently developed by [START_REF] Park | Analytical solutions for solute transfer characteristics at continuous fracture junctions[END_REF], is an alternative mixing rule between the perfect mixing model and the streamtube model. The solute molecules are assumed to follow the streamlines from the inflow bonds to the outflow bonds, but are allowed to diffuse between the streamlines, leading somehow to "partial-mixing" at the fracture junction. The mixing rate depends on the relative importance of advective versus diffusive mass transfer, expressed as a node-Péclet number Pe, and defined as: of the solute in free water. [START_REF] Park | Analytical solutions for solute transfer characteristics at continuous fracture junctions[END_REF] have shown that the perfect-mixing and the streamtube models can be considered as two "end-member" cases of their diffusional-mixing model.

For the fracture networks handled by the SOLFRAC program, there is only one type of fracture junction for which the choice of the mixing model is important. This is the so-called "continuous junction" case, corresponding to the junction of two contiguous inlet bonds and insu-00260744, version 1 -5 Mar 2008 two contiguous outlet bonds (Fig. 4). For all the other types of fracture junction (one inlet and three outlet bonds, three inlet and one outlet bonds, or one or more non-flowing "dead-end" bonds), it can be easily shown that the streamline patterns leads to the same mass partitioning whatever the mixing model used. Note that the particular case of "discontinuous junction" (i.e. two inlet and two outlet bonds on the opposite sides of an intersection), studied by [START_REF] Hull | Streamline routing through fracture junctions[END_REF], may only occur if fluid sources or sinks are present within the flow domain, which is not provided for by the SOLFRAC program.

For all the "non-continuous" fracture junctions, particle-transition rules in the program are based on the perfect mixing model, and are written as:

∑ - = Q Q p j ij (52)
where p ij [-] is the probability of particle transition from an inlet bond i to an outlet bond j, Q j [L 3 •T -1 ] is the flow rate in the outlet bond j, and ΣQ -[L 3 •T -1 ] is the sum of the discharge flow rates over all the outlet bonds connected to the junction of interest.

For the continuous junctions, particle-transition rules depend on the mixing model chosen by the user. With the perfect-mixing model, the particle transition probabilities are calculated according to Eq. ( 52). Note that with this model, mass sharing is not influenced by the inlet fluxes. With the streamtube model, mass sharing depends on both the inlet and outlet fluid fluxes. The solute tends to flow preferentially into the outlet bond that is directly contiguous to the inlet where it comes from. Referring to the case illustrated in Fig. 4 where the solute particles come from the inlet bond i = 2, one can distinguish two possibilities [START_REF] Park | Analytical solutions for solute transfer characteristics at continuous fracture junctions[END_REF]:

⎩ ⎨ ⎧ = = ⇒ ≤ 0 1 24 23 3 2 p p Q Q (53) insu-00260744, version 1 -5 Mar 2008 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ - = = ⇒ > 2 3 2 24 2 3 23 3 2 Q Q Q p Q Q p Q Q (54)
Using the same flow pattern but with solute particles coming from the inlet bond i = 1, the transition probabilities are written as:

⎩ ⎨ ⎧ = = ⇒ ≤ 1 0 14 13 4 1 p p Q Q (55) ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = - = ⇒ > 1 4 14 1 4 1 13 4 1 Q Q p Q Q Q p Q Q (56)
The mathematical expressions for the diffusional-mixing model of [START_REF] Park | Analytical solutions for solute transfer characteristics at continuous fracture junctions[END_REF] are somewhat bigger, but do not entail any calculation problem. Referring to the above two cases, the transition probabilities are written as:

= + = + = C C C C C Q C Q C Q p C Q C Q C Q p C Q C Q C Q p C Q C Q C Q p (57)
with the expressions of C 3 and

C 4 for 1 4 Q Q ≤ : ( ) ( ) ( ) ( ) 2 2 22 2 2 0 3 0 22 22 2 2 32 exp exp erf erf C C C C δ δ δ π ⎡ ⎤ -Π - -Π -⎢ ⎥ = + +Π Π -Π Π Π ⎢ ⎥ ⎣ ⎦ (58) ( ) ( ) ( ) ( ) 2 2 11 2 1 0 4 0 11 11 2 2 41 exp exp erf erf C C C C δ δ δ π ⎡ ⎤ -Π - -Π -⎢ ⎥ = + +Π Π -Π Π Π ⎢ ⎥ ⎣ ⎦ (59) insu-00260744, version 1 -5 Mar 2008 while, for 1 4 Q Q > : ( ) ( ) ( ) ( ) 2 2 22 1 2 0 3 0 2 2 2 2 1 1 32 exp exp erf erf C C C C δ δ δ π ⎡ ⎤ -Π - -Π -⎢ ⎥ = + +Π Π -Π Π Π ⎢ ⎥ ⎣ ⎦ (60) ( ) ( ) ( ) ( ) 2 2 11 1 1 0 4 0 11 11 1 1 41 exp exp erf erf C C C C δ δ δ π ⎡ ⎤ -Π - -Π -⎢ ⎥ = + +Π Π -Π Π Π ⎢ ⎥ ⎣ ⎦ (61)
and where: Many attempts, either experimental, analytical, or numerical, have been undertaken by several authors to determine the best suited mixing-model for simulating solute mass partitioning at continuous fracture junctions [START_REF] Hull | Streamline routing through fracture junctions[END_REF][START_REF] Robinson | A laboratory and numerical investigation of solute transport in discontinuous fracture systems[END_REF][START_REF] Berkowitz | Mass transfer at fracture intersections: An evaluation of mixing models[END_REF][START_REF] Park | Analytical solutions for solute transfer characteristics at continuous fracture junctions[END_REF][START_REF] Mourzenko | Solute transport at fracture intersections[END_REF]. The most recent is that by [START_REF] Mourzenko | Solute transport at fracture intersections[END_REF], who performed three-dimensional particle tracking simulations over both parallel-plate and rough-walled fracture intersections. The following conclusions can be drawn from their results:

1 1 2 2 0 1 2 u C u C C u u + = + (62) 4 2 i ij m j m Q D u b Π = (63) 3 2 4 2 i m i m Q Q D u b δ - Π = (64) D m [L 2 •T -1 ]
insu-00260744, version 1 -5 Mar 2008

1. In the case of a parallel-plate fracture intersection, the streamtube model is suited for node-Péclet numbers Pe ≥ 30 (see Eq. 51). However, the sharp mass partitioning predictions in expressions ( 53) and ( 55) are unrealistic when the difference between the flow rates (Q 2 and Q 3 ) or (Q 1 and Q 4 ) is very small. The diffusional-mixing model is more universal since it applies fairly well for Pe ≥ 10. The perfect-mixing approximation is only valid for Pe values close to unity. For very low Péclet values, (i.e. Pe << 1), the solute transfer through the fracture junction becomes mostly diffusive and none of the above three mixing models is valid. Note that such eventuality occurs rarely in practice.

2. In the case of a rough-walled fracture intersection, the streamtube, diffusional mixing, and perfect mixing models remain qualitatively valid over Pe ranges similar to those specified above for a parallel-plate fracture intersection. However, the streamtube and diffusional mixing models slightly underestimate the solute mixing at fracture intersection for large Péclet numbers.

Breakthrough curve computation

Solute breakthrough curves are calculated from the residence time distribution of the particles in the network, according to the following algorithm:

1. Determination of the minimum and maximum residence times t min , t max of the whole set of where Np n is the number of particles leaving the network during the time step Δt n , and mp out (j) is the mass of the j th particle leaving the network during Δt n (see Eq. 34);

4. Calculation of the mean concentration at the outlet boundary of the network for each time step. In the case of short-term injection of finite mass M 0 (see above):

( ) ( )

1 n n out m t c t t Q Δ Δ = Δ ∑ (66)
where

ΣQ out [L 3 •T -1 ]
is the sum of the discharge flow rate over all the bonds connected to the outlet boundary of the network. In the case of a continuous injection of constant concentration C 0 (see above):

( ) ( ) min min 2 1 t n t n t c t c t dt + Δ Δ = ∫ (67)
In the program, the integral in expression ( 67) is numerically evaluated with an evolutive Simpson's rule (see [START_REF] Press | Numerical recipes in C: the art of scientific computing[END_REF] pp 132-134).

During the transport simulations, the history of each particle in the network is recorded as a series of data pairs describing the succession of the nodes encountered, and the arrival times of the particle at these nodes. These data may be used to compute the solute breakthrough curve at any observation node specified by the user, according to an algorithm similar to that described above. Note also that these data pairs further enable to analyze the solute-plume dispersion at the network scale (see hereafter).

The accuracy of the simulated breakthrough curves depends on both the number N of times intervals (i.e. the number of points c(Δt n ) on the graph), and the number of particles used in the simulation. The number of particles to be used increases with N because a minimum number of particles is required within each time interval in order to ensure the convergence of insu-00260744, version 1 -5 Mar 2008 c(Δt n ). Furthermore, the number of particles needed increases with the complexity of the network. From a practical point of view, it is difficult to estimate the optimum number of particles to be used in a given network. An empirical method consists in repeating simulations, increasing each time the number of particles, up to reach no significant difference between the simulated breakthrough curves. As an example, about 10 4 particles are required for simulating accurately the transport problem illustrated in Fig. 3, with 100 points on the breakthrough curve. Note that in a Monte Carlo framework involving a large number of simulations, the number of particles may be lower because the oscillations resulting from the discrete nature of the particles are random and cancel each other when ensemble statistics are computed (see e.g. [START_REF] Hassan | On using particle tracking methods to simulate transport in single-continuum and dual continua porous media[END_REF].

Verification problems

The purpose of this chapter is to verify the accuracy of the program through comparisons with analytical and semi-analytical solutions of solute transport in fractured media. Although not mentioned in previous papers, the SOLFRAC program has been used by [START_REF] Delay | Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks[END_REF] and Bodin et al. (2003c) for assessing the accuracy of the TDRW method through several basic transport problems. [START_REF] Delay | Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks[END_REF] first addressed the case of nonreactive solute transport in a single fracture-matrix system. The transport processes considered were advection-dispersion in the fracture plane, and matrix diffusion. Simulation results compared very well with the analytical solution of [START_REF] Sudicky | Contaminant transport in fractured porous media: Analytical solutions for a system of parallel fractures[END_REF]. This model verification was extended by Bodin et al. (2003c) to the case of reactive solute transport, including sorption on the fracture walls and radioactive decay. Bodin et al. (2003c) also showed the accuracy of the TDRW method (and of the SOLFRAC code) for simulating advection and hydrodynamic dispersion in a synthetic fracture network with sharp contrasts in insu-00260744, version 1 -5 Mar 2008 dispersion coefficients. The series of test problems presented in this section has been designed

to verify the accuracy of the program in transport scenarios that have not yet been addressed.

Tests SF1 and SF2: Scale-dependent dispersion in a single fracture

The purpose of the tests SF1 and SF2 is to investigate the ability of the program to deal with scale-dependent dispersion within the bonds of a fracture network. An advection-dispersion problem in a single fracture is addressed and dispersion is assumed to increase linearly (test SF1) or exponentially (test SF2) with distance along the fracture axis. A source of solute of constant strength is assumed to supply a continuous injection of constant concentration at the fracture inlet. Analytical solutions to the above-described problems were developed by [START_REF] Yates | An analytical solution for one-dimensional transport in heterogeneous porous media[END_REF][START_REF] Yates | An analytical solution for one-dimensional transport in porous media with an exponential dispersion function[END_REF]. The values of the parameters used for the calculation of both the analytical curves and SOLFRAC simulations are listed in Table 1. As shown in Fig. 5, simulations and analytical solutions are in very good agreement.

Test DFN1: Reactive solute transport in a fracture network

The test DFN1 performed over a synthetic discrete fracture network involves the following mechanisms: (1) advection and hydrodynamic dispersion in the fractures, (2) matrix diffusion,

(3) solute sorption on the fracture walls and in the matrix, and (4) mass sharing at fracture intersections. The network size is 350 x 350 m, and the constant hydraulic head values on the top and bottom boundaries are fixed to 100 m and 99 m, respectively. The flow takes place within two orthogonal sets of fractures, yielding 71 flowing bonds (Fig. 6). The fracture apertures and widths are constant over the network and set up to 2b = 2.5 x 10 -4 m and W = 1 m, respectively. A mass m 0 = 10 -2 g of solute is instantaneously injected at the inlet B 1 of the network (see Fig. 6). The values of transport parameters in the fractures and the rock matrix insu-00260744, version 1 -5 Mar 2008 are listed in Table 2. The theoretical breakthrough curve c out (t) at the outlet of the network is calculated as the summation of mass fluxes from each elementary flow path in the network (Bodin et al., 2003c):

( ) ( ) 1 EP N i i out out F t c t Q = = ∑ (
( ) ( ) ( ) ( ) 1 0 , 1 1 2 1 i i N i n n N n F t m p t t t - + = ⎛ ⎞ = ℑ * ℑ * ℑ ⎜ ⎟ ⎝ ⎠ ∏ (69)
where N i is the number of bonds in the elementary path i, and p n,n+1 is the fraction of solute mass in bond n that flows into bond n+1, this fraction depending on the flow configuration and on the mixing model at the fracture intersection (see expressions 52, 53-56, and 57).

Because the transport problem DFN1 involves both advection-dispersion in fractures and matrix diffusion, the convolution product in (69) cannot be simply replaced by an equivalent probability density function as the one developed by Bodin et al. (2003c). We rather propose to use the Laplace transform of ℑ n (t), which is easily derived from the work by [START_REF] Tang | Contaminant transport in fractured porous media: Analytical solution for a single fracture[END_REF]:

( ) ( ) ( ) ( )

1 2 1 2 2 exp exp 1 n s s L L s A ν ν β ⎡ ⎤ ⎧ ⎫ ⎛ ⎞ ⎪ ⎪ ⎢ ⎥ ℑ = - + + ⎨ ⎬ ⎜ ⎟ ⎢ ⎥ ⎪ ⎪ ⎝ ⎠ ⎩ ⎭ ⎣ ⎦ (
1 0 , 1 1 1 i i N N i n n n n n F s m p s - + = = ⎛ ⎞⎛ ⎞ = ℑ ⎜ ⎟⎜ ⎟ ⎝ ⎠⎝ ⎠ ∏ ∏ (75)
The substitution of the inverse transform of ( ) i F s for F i (t) in ( 68) yields the following semianalytical (SA) solution:

( ) ( ) { } ( ) ( ) 1 1 1 1 1 1 0 , 1 1 1 1 EP EP i i EP N N i N N i N i i out n n n i n n out out out L F s L F s m c t L p s Q Q Q - - - = - = + = = = ⎧ ⎫ ⎨ ⎬ ⎧ ⎫ ⎡ ⎤ ⎛ ⎞⎛ ⎞ ⎪ ⎪ ⎩ ⎭ = = = ℑ ⎨ ⎬ ⎢ ⎥ ⎜ ⎟⎜ ⎟ ⎝ ⎠⎝ ⎠ ⎪ ⎪ ⎣ ⎦ ⎩ ⎭ ∑ ∑ ∑ ∏ ∏ (76)
For the test DFN1, the inverse Laplace transform L -1 {} in (76) has been numerically evaluated using the MATLAB® routine "Invlap.m" written by [START_REF] Hollenbeck | INVLAP.M: A matlab function for numerical inversion of Laplace transforms by the de Hoog algorithm[END_REF] This spreading results from the cumulated effects of the flow-field heterogeneity on the network scale (controlled by the fracture network geometry), the mass sharing at fracture intersections, and the transport processes on the bond-scale. In the "ideal case" of a homogeneous 1D transport problem, the following relations hold: The concept of "dilution index" was developed by [START_REF] Kitanidis | The concept of the dilution index[END_REF] to better differentiate between spreading and dilution in heterogeneous media. The dilution index characterizes the volume of porous medium occupied by the solute at a given time. Mathematically, it is defined as: where V is the volume (or length) of the studied system, c(x, t) is the solute concentration at location x, and M(t) is the total mass of solute in the system at time t:

3 2 2 3 2 2 t F t F U D x D U L σ σ Δ = ⇔ = (77) 2 2 2 2 x x F F D t D t σ σ = ⇔ = (
T D T T σ = (82) ( ) ( ) ( ) 2 2 2 _ 2 2 1 2 2 obs x obs x obs x xx eff obs t T T t T t D T t t σ σ σ = + Δ - -Δ ∂ = = ∂ Δ (83) ( ) ( ) ( ) 2 2 2 _ 2 2 1 2 2 obs y y o b s y o b s yy eff obs t T T t T t D T t t σ σ σ = ∂ + Δ - -Δ = = ∂ Δ (
( ) ( ) ( ) ( ) ( ) , , exp ln V c t c t E t dV M t M t ⎡ ⎤ ⎛ ⎞ = - ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ∫ x x (
( ) ( ) , V M t c t dV = ∫ x (86)
As pointed out by [START_REF] Park | Effects of junction transfer characteristics on transport in fracture networks[END_REF], the dilution process in fracture networks is closely related to the solute mass partitioning among the multiple available flow paths. Today, it is well known that solute transport in natural fracture networks is often channelised, i.e. most of the mass transfer occurs in a limited number of flow paths (see e.g. [START_REF] Tsang | Flow channeling in heterogeneous fractured rocks[END_REF] and references therein). [START_REF] Park | Effects of junction transfer characteristics on transport in fracture networks[END_REF] suggested that the dilution index might be used as a quantitative measurement of the degree of channelised transport in fractured media.

The program enables to compute curves showing the evolution of the dilution index with respect to the mean position of the particle cloud. As suggested by [START_REF] Park | Effects of junction transfer characteristics on transport in fracture networks[END_REF], this type of plot is useful to analyze the dilution features in a fracture network. In the program, the dilution index is calculated as follows:

( ) ( ) ( ) where N B is the total number of flowing bonds in the fracture network (i.e. the number of pipes in the backbone of the network), NP n (t) is the number of particles in bond n at time t, and NP tot (t) the total number of particles in the whole network at time t.

1

Summary and conclusion

In 1993, [START_REF] Smith | Solute transport through fracture networks[END_REF] rightly argued that the computational constraints on the total number of fractures that can be included in a discrete network severely limited the applicability of this modelling approach to practical cases. This was justified a dozen years ago, but the computation capacities have strongly increased and today, representative discrete models become affordable to handle transport problems on the scale of a reservoir. However, insu-00260744, version 1 -5 Mar 2008

it is important to note that the question of parameterisation remains unanswered but still crucial since numerical simulations rely on available data (see e.g. [START_REF] Smith | The development of a safety assessment model of the geosphere for a repository sited in the crystalline basement of northern Switzerland[END_REF].

The software presented in this paper performs fast simulations of solute transport in complex

  1) where c f [M•L -3 ] is the solute concentration in the flow channel; λ [T -1 ] is the decay constant of the solute; t [T] is the time variable; x [L] is the space coordinate along the flow channel; u f [L•T -1 ] is the fluid velocity in the flow channel; D f [L 2 •T -1 ] is the hydrodynamic dispersion coefficient in the flow channel; D m [L 2 •T -1 ] is the molecular-diffusion coefficient of the solute; a [L] is the half-width of the flow channel; c nf [M•L -3 ] is the solute concentration in the non-flowing part of the fracture plane; y [L] is the space coordinate in the fracture plane, insu-00260744, version 1 -5 Mar 2008 perpendicular to the flow channel axis; D e [L 2 •T -1 ] is the effective-diffusion coefficient of the solute in the matrix; b [L] is the half-aperture of the fracture; c m [M•L

  where D a [L 2 •T -1 ] is the apparent-diffusion coefficient of the solute in the matrix, expressed as a function of θ m[-] the matrix porosity, ρ m [M•L -3 ] the bulk density of the matrix, and K m [L 3 •M -1 ] the volumetric-sorption coefficient of the solute in the matrix:

  are respectively the corrected mean and variance of the advectiondispersion travel time, based on an apparent fluid-velocity *

  38)where dispersivity α [L] is proportional to the correlation length of the flow field. Note that[START_REF] Gelhar | Three-dimensional stochastic analysis of macrodispersion in aquifers[END_REF] propose analytical expressions for the calculation of α, from the insu-00260744, version 1 -5 Mar 2008 correlation length and the standard deviation of the logarithm of the apertures in the fracture plane. Another part of solute dispersion results from the combination of molecular diffusion and velocity variations across the fracture aperture. This leads to the well-known Taylor-Aris dispersion, characterized by the dispersion coefficient[START_REF] Dewey | Longitudinal dispersion in flows that are homogeneous in the streamwise direction[END_REF]: is theoretically valid only beyond a critical travel time τ c , which corresponds to the minimum duration needed for a particle to experience the whole crosssectional parabolic profile of velocities across the fracture aperture. This critical time is proportional to a characteristic time of transverse diffusion: τ c ∝ b 2 /D m . In other words, the solute must travel over a minimum distance c

  41)In the program, the time-dependent Taylor-Aris dispersion of expression (40) is truncated to the third term. By applying the change of variable t = x/u f , the spatial derivative of expression (41) yields the derivative of D f used in (10

  ) where ϖ [-], ξ [-] and κ [L -1 ] are user-defined constants, Ψ [L] is the asymptotic dispersivity, η [L] is the migration distance for which dispersivity is half its asymptotic value, and Φ [L] is the maximal value of dispersivity in the exponential model. The derivative and/or integral of D f in the case where both D L1 and D L2 are scale-dependent are straightforward and not developed here. Note that when simulating advective-dispersive transport in the bonds and matrix-diffusion, the mean fluid velocity u f in expressions (38-42) is replaced by the corrected velocity u f * (Eq. 26).Transport simulations on the network scaleTransport simulations on the network scale are performed by tracking a set of particles injected into the flow field, until all the particles have left the network. For each particle, the TDRW algorithm is duplicated over the series of bonds experienced by the particle from its insu-00260744, version 1 -5 Mar 2008 injection point up to its exit point at the downstream boundary of the network. Solute breakthrough curves are then calculated from the residence time distribution of the particles in the network.

  u [L•T -1 ] is the mean fluid velocity in the inlet bonds connected to the junction, b [L] is the half mean-aperture of the inlet and outlet bonds, and D m [L 2 •T -1 ] is the molecular diffusion

  is the molecular diffusion of the solute, b m [L] the mean aperture of the bonds connected to the fracture junction, and u i [L•T -1 ] the fluid velocity in bond i. Note that the above expressions are those provided by Mourzenko at al. (2002) who corrected a few misprints of the initial work by Park and Lee (1999).

  of the time-interval [t min -t max ] in N time steps of duration Δt; 3. Summation of the mass associated with the particles leaving the network during each time step Δt n (n = 1 .. N):

  68) where N EP is the number of elementary paths starting from the injection point (N EP = 1800 in the current problem), F i (t) [M•T -1 ] is the solute mass flux in the i th elementary path, and Q out [L 3 •T -1 ] is the total flow rate at the outlet boundary of the network. The solute mass flux F i (t) can be expressed as a convolution product of the probability density functions ℑ n [T -1 ] of residence times in each bond of the elementary path i:

  s [-] is the Laplace parameter, L [L] is the length of the bond n, u n [L•T -1 ] and D n [L 2 •T -1 ] are the flow velocity and the dispersion coefficient in the bond, respectively, and R f [-] is a retardation coefficient due to solute sorption on the fracture walls. As stated above, the Laplace transform of a convolution is the product of the individual transforms. Thus, the Laplace transform of the solute mass flux in the i th elementary path is:

  , which is based on the De Hoog et al. (1982) algorithm. Both SA and SOLFRAC breakthrough curves have been computed using successively the perfect-mixing model, streamtube model, and diffusionalmixing model for mass partitioning at fracture junctions. As shown in Fig. 7, SA and SOLFRAC results coincide very well in each case.Macroscopic time-and scale-analysis of solute transport insu-00260744, version 1 -5 Mar 2008As stated earlier, the program can analyze various transport features on the network scale, either from the spatial distribution of particles at a given time or from the residence time distribution (RTD) of the particles in the network. These transport features are: (1) timeevolution of the mean position of the particle cloud, (2) time-evolution of the longitudinal and transverse variance of the particle cloud, (3) computation of the macroscopic longitudinal and transverse dispersion coefficients at a given time, (4) time-evolution of the macroscopic dispersion, (5) scale-evolution of the macroscopic dispersion according to the mean position of the particle cloud, and (6) evolution of the dilution index according to the mean position of the particle cloud. The term "macroscopic dispersion" is used here by analogy to transport in porous media, and describes the spreading of the solute/particle plume on the network scale.

  2 ] is the variance of the residence time distribution of the particles for a transportlength L, σ x 2 [T 2 ] is the variance of the spatial distribution of the particles at a given time t, D F [L 2 •T -1 ] is the (Fickian) dispersion coefficient, and U [L•T -1 ] is the mean flow velocity. The macroscopic dispersion (or "macrodispersion") coefficients computed in SOLFRAC rely on various definitions, which are all basically based on the relations above. Two expressions are based on the Residence Time Distribution (RTD) between injection and one location of the network: insu-00260744, version 1 -5 Mar 2008

  84)insu-00260744, version 1 -5 Mar 2008 where σ x 2 (T obs ) and σ y 2 (T obs ) are the spatial-variance of the particle cloud projected on the xdirection parallel to the flow gradient and y-direction perpendicular to x, at time T obs .Coefficients D xx_app and D yy_app are termed "apparent" because the analysis assumes spatial homogeneity along the travel distance, i.e. the definition is similar to that in (78) of a homogeneous Fickian dispersion. On the other hand, coefficients D xx_eff and D yy_eff are termed "effective" because their calculation only considers the temporal-increment in the plume development without any reference to homogeneity. In the program, the duration of Δt in (83) and (84) is that of the time-steps used for the computation of breakthrough curves (see Eq.

  66). SOLFRAC enables to analyze the time-evolution of macrodispersion by computing D xx_app , D yy_app , D xx_eff , and D yy_eff for several T obs values, evenly distributed between 0 and the first particle arrival-time at the network outlet. It is also possible to analyze the scaleevolution of macrodispersion by plotting D xx_app and D xx_eff values with respect to the x-mean position of the particle cloud for the same T obs values. Note that performing an injection over the whole inlet boundary of the network, which makes the particles to experience most of the flowing bonds, is preferable for analyzing the longitudinal macrodispersion D xx . On the other hand, a single injection point is preferable for analyzing the transverse macrodispersion D yy , in order to ensure that σ y 2 = 0 at t = 0. This injection point may be randomly chosen among the nodes located at the inlet boundary of the network.

  85) insu-00260744, version 1 -5 Mar 2008

Fig. 1 .Fig. 2 .

 12 Fig. 1. Modelling approaches for 3D fracture networks. A: discretisation of fracture planes,

Fig. 3 .Fig. 5 .Fig. 6 .Fig. 7 .

 3567 Fig. 3. Identification of elementary paths in the SOLFRAC program for an injection in (A).

  Δl io[L] is the longitudinal distance between injection and observation points in the network, t m [T] and σ t [T] are the mean and standard deviation of the RTD, respectively, and t 16 and t 84 are the time values for which the cumulative probability density function of the particle RTD is equal to 0.16 and 0.84, respectively. Expression (80) was first proposed by[START_REF] Fried | Dispersion in porous media[END_REF], and is similar to Eq. (79) except that the calculation of the dispersion coefficient relies on a truncated RTD. This truncation forces the relative influence of both rising and lowering parts of the breakthrough curve to balance in the computation of

		( ) io D l Δ =	2 3 io m l t Δ 2	σ	2 t	(79)
	( ) io D l Δ =	8 io l Δ	⎛ ⎜ ⎜ ⎝	1	-	16 16 t t t	m	-	1	-	84 84 t t t	m	⎞ ⎟ ⎟ ⎠	2	(80)
	where ( ) _ xx app obs D T	=	σ	2 2 x	( ) obs obs T T	(81)
	_ yy app	( ) obs					2 y 2	( ) obs obs

D, which may be useful when the time-evolution of concentrations exhibits an extensive tailing. The calculation of a macroscopic dispersion coefficient from the spatial distribution of particles at a given time may also be based on different approaches. Let us denote D xx and D yy the longitudinal and transverse macrodispersion coefficients, respectively. The program computes D xx and D yy from either an "apparent" or "effective" definition

[START_REF] Dagan | Flow and transport in porous formation[END_REF][START_REF] Jussel | Transport modeling in heterogeneous aquifers, 2, Three-dimensional transport model and stochastic numerical tracer experiments[END_REF]

):

Acknowledgments. We are grateful to the "French National Program for Research in Hydrology" (PNRH) for the financial support of this work. insu-00260744, version 1 -5 Mar 2008