Predicting solute transport in heterogeneous media from results obtained in homogeneous ones: an experimental approach

JOURNAL OF

Contaminant Hydrology ELSEVIER

Journal of Contaminant Hydrology 25 (1997) 63-84

Predicting solute transport in heterogeneous media from results obtained in homogeneous ones: an experimental approach Frrdrrick Delay a,*, Gilles Porel b, Ghislain de Marsily

a

a Laboratoire G~ologieAppliqu~e,

URA CNRS 1367, box 123, 4, Place Jussieu F-75252 Paris Cedex 05, France b Laboratoire Hydrog~ologie, URA CNRS 721, The University of Poitiers, Avenue du Recteur Pineau, F-86000 Poitiers, France

Received 13 November 1995; accepted 17 April 1996

Abstract The aim of this paper is to show that the hydrodispersive transport parameters of a solute in a heterogeneous medium can be predicted using parameters determined by several experiments with the different homogeneous media that constitute together the heterogeneous one. Elementary rules for predicting these parameters are defined. Three types of experiments were carried out. (1) Homogeneous medium transport experiments in 1-D columns. They were fitted with the help of a new numerical particle tracking method, including advection-dispersion and exchange between mobile and immobile water with first-order kinetics. (2) Layered heterogeneous media. The experiments were simulated by the numerical model without fitting, using the parameters of each of the layers of the homogeneous media included in the heterogeneous one. The results prove that, in such layered media, the global transport problem can be considered as the convolution of elementary transports in homogeneous blocks. We will therefore conclude that the reservoir simulation methods which represent the heterogeneous reservoir as a random set of elementary homogeneous blocks can be used to represent transport if each elementary block is given its appropriate hydrodispersive parameters. (3) Homogeneous mixtures of two pure components. The experiments produced three interesting results. The immobile water volume is the direct linear combination of its equivalent in the two components, weighted by the volumetric fraction of the component in the mixture. The dispersivity of the medium also evolves linearly, but is affected by an initial sharp rise due to a strong increase in local flow path heterogeneities whenever a "foreign body" is introduced into a pure component. The exchange rate between mobile and immobile water is also perturbed by

* Corresponding author. 0169-7722/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PH S 0 1 6 9 - 7 7 2 2 ( 9 6 ) 0 0 0 2 0 - 4

64

F. Delay et al. / Journal of Contaminant ltydrology 25 (1997) 63-84

mixing. In a binary mixture of one component with dual porosity and another one without, the exchange rate decreases with the increase of the volumetric fraction of the medium without dual porosity. It seems that, even if it does not modify the volumes of immobile and mobile water, it reduces the surface/volume ratio of the solid medium and consequently the exchange rate. Mixing rules could be adopted in order to modify the values of the hydrodispersive parameters that need to be introduced for the elementary blocks of heterogeneous reservoirs. © 1997 Elsevier Science B.V. Keywords: Advection-dispersion;Matrix diffusion; Dual porosity; Reservoir simulation

1. Introduction

Solute transport in heterogeneous media is a problem which, at present, still offers a promising future for numerical modeling. In this paper we will assume as a first approximation that a heterogeneous medium is a succession of " b l o c k s " of homogeneous media. Each homogeneous medium included in the heterogeneous one is assumed to have the following characteristics: - - a macro-porosity through which fluid flow occurs; -a micro-porosity (or matrix porosity) containing an immobile fluid capable of exchanging solute with the mobile phase. In the presence of an immobile fluid phase, the advection-dispersion equation for one-dimensional flow, for any homogeneous block, is: 52 C m c~u ~x ~

0C m -

u ~)x

0C m -

~)t

( L - w~) ~tCim +

w~

St

(l)

where C m = volumetric concentration in the mobile fluid [M L 3]; Ci m = volumetric concentration in the immobile fluid [M L-3 ]; c~ = longitudinal dispersivity [L]; u = mean pore velocity [L T -1 ]; t = time variable [T]; x = space variable [L]; wc = kinematic porosity of the porous medium (dimensionless); wt = total porosity of the medium; and L = dimensionless parameter taking a value between w~. and wt. L - wc represents the volumetric fraction of immobile fluid in the porous medium. It seems to us that two types of problems are the most important. The first one is mainly numerical. Although Eq. (1) is simple, solving it is not without problems. Methods such as those of finite differences and finite elements (Kinzelbach, 1986; Bobba, 1989) are difficult to use when solute exchange between the mobile and the immobile phase is added, as is the case here. The difficulties increase if one is dealing with a strongly heterogeneous medium with, for instance, sudden variations in the dispersive properties or in the "adsorption capacities" of the medium (Van Leer, 1979; Putti et al., 1990). Such constraints make it necessary to solve sparse matrices with a high variation coefficient which are numerically unstable. Classical methods of the " r a n d o m walk" type simulate the physical process of dispersion by a stochastic analysis of particles taken to represent the solute. They can be used on dual-porosity or fractured media (Ackerer and Kinzelbach, 1985; Ackerer, 1988) without creating any particular problem for writing of the algorithms. However, since these models manage the

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84

65

movement of each particle individually, they consume a great deal of computer time unless the total number of particles is severely limited. The second type of problem is added to the preceding whenever we are dealing with a heterogeneous medium and concerns the identification of transport parameters. Even if we have some "local data" on parameters that govern the transport over a short distance, how do we extend the numerical definition of these parameters to the entire domain to be modeled? One method is geostatistical conditional simulations (Journel and Huijbregts, 1978; Delhomme, 1979; Gomez-Hernandez and Srivastava, 1990; Gomez-Hernandez and Journel, 1992; de Marsily, 1993; Delay et al., 1995) which generate random fields consistent with the measured data. However, some questions remain unanswered. Often the numerical image of a heterogeneous reservoir is created by merely juxtaposing homogeneous elementary blocks. Can we be sure that this kind of representation does not alter the nature of the solute transport? In other words, can we solve the elementary transport equation in each successive homogeneous block, and thus predict the transport through several blocks in series without changing their parameters? This bring us to another question. Suppose that we have measured hydrodynamic and hydrodispersive parameters at the scale of one block for two typical homogeneous facies of a reservoir, e.g., pure clay and pure sand of an alluvial aquifer. If at a point in space, a numerical geostatistical simulation of the reservoir shows a high probability that the facies corresponds to clay, we know that, even if the simulation is close to reality, the true geological facies is probably not pure, but is a mixture where clays predominate. It would a priori be preferable to use the parameters of a clay-sand mixture in this reservoir block. Is it possible to predict these parameters from the data obtained at the same scale with pure components? The aim of this paper is to give some answers to the preceding questions on the basis of experimental one-dimensional (l-D) tracer tests and their interpretation by a numerical method for solving the transport Eq. (1). As we know that classical methods may give rise to numerical problems, we will use a new one (Delay et al., 1993, 1994) which has proved its worth in heterogeneous adsorbing media. For the sake of clarity, we begin by describing the calculation principles.

2. General principles of the calculation model The calculation can be described as a hybrid between discrete schemes and particle tracking. Transport is represented by particles moving in a grid. The particles are no longer managed individually, but as numerical variables representing their number allocated to each elementary cell of the grid. 2.1. Advection-dispersion

In one dimension, consider elementary cells placed side by side at a constant distance of A x and through which flow transports a solute. Convection is obtained by moving the particles from cell i - 1 to cell i, with a time step At = A x / u . The number of

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63 84

66

i" ~x-I a

I Convectioni-,-i ,~,

k

~,

,~.

~'~I

I ~'~

k

~.~, "~i

b Dispersion Fig. 1. Principle of the advective-dispersive jump in one-dimensional flow between times t and t + At involving the particles in cell i - 1. During At, all the particles in cell i 1 are moved by advection to cell i (top). Dispersion distributes them uniformly along a length 2D, centered in i. particles c o n t a i n e d in the free water in cell i - 1 is assigned to the n u m e r i c a l variable which represents the n u m b e r of particles in the free water of cell i (Fig. la). For the dispersion, a purely deterministic form of the results of U f f i n k ' s (1985) probabilistic approach is used, adding a dispersive displacement of the particles to the advective one. At each a d v e c t i v e - d i s p e r s i v e " j u m p " from a cell i - 1 to cell i, the Ni- 1 particles of cell i - 1 are uniformly distributed around a point centered on i, along a distance _+ D (Fig. lb), with:

D= +

(2)

where D = m a x i m u m length of the dispersive m o v e m e n t [L]; a = longitudinal dispersivity [L]; u = m e a n pore velocity [L T - j ]; A x = length of a cell [L]; and At = length of time step needed to cross A x by pure c o n v e c t i o n [T]. With a constant linear particle density, each cell entirely covered by s e g m e n t 2 D receives: (N,. i A x ) / 2 D particles; the two end cells receive: [ N , _ t { D - ( k + 0.5)A x } ] / 2 D particles, k is the integer n u m b e r of cells other than i which are entirely covered by D. Like U f f i n k ' s method, this one generates, after a certain n u m b e r of j u m p s , a n o r m a l distribution around the c o n v e c t i o n point in accordance with F i c k ' s law.

2.2. Exchange between the mobile and immobile phase For the sake of simplicity, the exchanges b e t w e e n the i m m o b i l e and m o b i l e fluid phases in the m e d i u m are represented by first-order kinetics. A s s u m e that, b e t w e e n time t' and time t' + At ( A t = A x / u ) , each cell is a closed system. D u r i n g the exchange, the total n u m b e r of particles in the cell remains constant:

U = N i m + N m = U i ° l +N~ ° -- Uim f + Nm~

(3)

where Ni','. , N',~ = n u m b e r of m o b i l e and i m m o b i l e particles, respectively, at time t" Nim, N~,~= their respective equivalents at time T (t' < T < t' + At); and Ni~~, Ni~ = their respective equivalents at time t' + At.

F. Delay et a L / Journal of Contaminant Hydrology 25 (1997) 63-84

67

Since the system is closed, the solution of the exchange problem can take two forms: - either the right-hand side of the advection-dispersion equation (1) is cancelled and one solves: aCm -

-

Ot'

L - w c OCim +

-

-

wc

-

-

-

0

(4)

Ot'

- - or one solves directly the equation of first-order kinetics: (L-

We) ~

(5)

= / 3 ( K p C m -- Cirri)

where L - w c = volumetric fraction of immobile fluid (dimensionless); Cim , C m = concentrations in the immobile and the mobile fluid phase, respectively [M L -3 ]; and /3 = exchange rate between the mobile and the immobile phase - - note that in reality, /3 is not a true first-order rate coefficient but rather the product of a mass transfer velocity (dimensions: [L T - I ] ) and the area/volume ratio of the granular medium (dimensions: [ L - l ] ) thus appearing as a first-order constant in [T 1]; and Kp = Cim//Cm ratio at equilibrium (dimensionless). The latter ratio is equal to 1 if the exchange only occurs between mobile and immobile fluid. But if there is also sorption on the solid, Kp can be globally greater than 1. Both methods lead to the same exact analytical solution (Delay et al., 1993): Nmf- RN°-Ni°mexp - / 3 L I+R

I+R wc

N L = U°m + Ui°m - N f

At

) + N°+Ni ° I+R (6)

where R=Kp[(L-Wc)/Wc]. With the analytical expressions obtained by globally managing the particles, the exchange simulation algorithm can only assign, to each cell, two variables representing, respectively, the number of particles in the mobile and immobile fluid.

2.3. Extension to heterogeneous media We want to make the main "hydrodispersive parameters" vary along the flow path. These parameters are: we , the kinematic porosity considered equal to the volumetric fraction of mobile fluid; a , the dispersivity of the medium; L - we, the volumetric fraction of immobile fluid; /3, the exchange rate; and Kv, the partition constant. The latter is not an adjustable parameter and its value is 1 if the exchange is between mobile and immobile water. Instead of using constant parameters for all cells, one specifies the parameter values for each cell. It is then easy to construct an 1-D "layered" heterogeneous medium. One only has to define the length of the sections in terms of the number of cells for which wc, a , L and /3 remain constant. This results in a medium with "dispersive character-

F. Delay et al./ Journal of Contaminant Hydrology 25 (1997) 63-84

68

istics" that evolve with the transfer distance. If A t represents the time step of an a d v e c t i v e - d i s p e r s i v e jump, we write:

Ax A xSwc ( min ) At-- - u(max) Q(t)

(7)

where Q(t) = flow rate in the system [L 3 T - I ] ; S = total cross-section of the flow [L2]; wc(min) = lowest value of kinematic porosity along the medium; and u(max) = highest mean pore velocity value. If A t remains constant throughout the column, the main difference with the homogeneous case is that the arrival point for pure advection will no longer coincide with the center of a cell (see Fig. l a ) and the algorithm has to be adapted so that it can calculate, for each cell of the column, the length covered by the segment of dispersion and define the indices of the cells beneath it. Note that we also have to add a correction of the advective j u m p to take into account the variations of the dispersion coefficient along the flow path (Ackerer and Kinzelbach, 1985; Ackerer, 1988; Uffink, 1990). The advective displacement now becomes:

x(t+±t)

-

x(t)

=

u+

a~

u--

+

~I

At

(8)

The first derivatives of the dispersivity c~ and of the mean pore velocity u with respect to the space dimension x are calculated by a centered difference approximation around a point situated at X ( t ) + uAt. Note that this enhanced particle tracking method was recently extended to two dimensions in Delay et al. (1996).

3. Experimental apparatus Most of the tracer experiments were made on a vertical Plexiglas column, 60 cm long and 5 cm in diameter. In some experiments, a column with the same diameter but 120 cm long was used. In order to avoid excessive differential compaction, the latter was installed horizontally. In all cases, the column is saturated and subjected to a constant flow rate by means of a peristaltic pump which feeds water into the system. All samples are taken at the outlet in order to systematically obtain a measure of the flowing concentrations. The transfer distance is varied by moving the injection point along the column. Short point injections, which are made with a syringe at the center of the section, introduce a concentrated solution of uranine into the medium. By adjusting as well as possible the injected mass, one makes the output concentrations, analysed by a Turner fluorine-meter, fall in the range of 1 - 5 0 0 0 p~g/1.

3.1. Homogeneous medium without matrix porosity This medium is made up of unfritted silica beads, 2 m m in diameter. The void volume of a pseudo-random arrangement of the beads can be measured by simply

69

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84

weighing it. One measures the weight difference between a given dry volume of silica beads and its equivalent saturated with water. The experiments give an intergranular porosity o f 41.4% which is taken as the reference value for the total porosity o f the medium. 3.2. Dual-porosity m e d i u m

This is also a granular medium, commercialized as catbox litter under the name " C a t s a n " . The material takes the shape o f roughly disk-shaped pellets, 4 m m in mean diameter and 1-1.5 m m thick. A granulometric analysis gives a quasi-linear distribution of the grain size with: minimum grain size, 0.8 mm, with a proportion in mass of 3%; and m a x i m u m grain size, 5 mm, with a proportion in mass of 30%. The visible intragranular porosity is considerable and the diameter of the pores in these disks is on the order of 0 . 1 - 0 . 5 mm. The total porosity is measured with the same method as that used for the silica beads. The total porosity can be taken as 73%. In order to estimate the matrix porosity, 150 g of water saturated material were kept at 105°C for 24 h and then weighed. The new weight was 98 g, so the evaporated mass of water was 42 g which, compared to the total volume, gives an approximate matrix porosity of 25%. The material was analysed by microprobe at 15 keV and is mainly composed of: Na, 0.24%; Mg, 0.15%; Ca, 19%; Si, 19%, in mass. X-ray diffraction by the D e b y e - S c h e r r e r method shows compositions o f quartz and tobermorite at 11 ,~ [stoichiometric formula A S T M sample No. 19.1364: C a s ( O H ) 2 S i 6 0 1 6 . 4 ( H 2 0 ) ] .

4. Tracer experiments on silica beads As there is no matrix porosity, the results of the tracer experiments can a priori be interpreted only on the basis of the advection-dispersion process. The characteristics of the experiments and their parameters to be fitted by the particle model are shown in Table 1. The four tracer experiments in Fig. 2 are fitted with the same kinematic porosity value wc = 40% and the same dispersivity value, c~ = 0.14 cm (Table 1). The model produces a good fit of the intrinsic dispersivity o f the medium which is independent o f variations in transfer distances and velocities. T h i s concerns, of course, variations of similar orders of magnitude.

Table 1 Tracer tests and fitting parameters in homogeneous segments of silica beads No.

x (cm)

Q (cm3/rain)

Mo (p,g)

c~ (cm)

wc (%)

$2 S1 $4 $3

54.5 26.5 54.5 26.5

20 20 40 40

200 200 200 200

0.14 0.14 0.14 0.14

40 40 40 40

wt (%)

70

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63 84 C (gg/1)

[ - - -

3500 3000

[

2000

Experimental curve[ Simulated curve Q = 20 cm3/min

S1

1500

2

1000

500

Time (min) 10

20

30

40

[ - - - Experimental curve / Simulated curve Q = 40 cm3/min

C (I.tg/l) 3500 3000

2000

$3

1500 1000 500 0

.

.

.

. 5

10

T i m e (min) ~ - i 15 20

Fig. 2. Fitting of tracer experiments on silica beads. For any hydrodynamic condition and transfer distance, the fittings are obtained with the same parameter set: dispersivity, a = 0 . 1 3 - 0 . 1 4 cm; kinematic porosity, w c = 0.4. The tracer transfer distances are 26.5 cm for S1 and $3, and 54.5 c m for $2 and $4. All injections are short point injections and the total cross-section of the column is 19.635 c m 2.

Note that the adjusted dispersivity does not depend on the mesh size of the model. From Eq. (2), it follows that a must be greater than A x / 2 4 for the algorithm to effectively disperse particles to other cells than their point of arrival by advection. Provided that the domain is finely discretized, it is therefore possible to handle low dispersivity values.

5. Tracer experiments in a dual-porosity homogeneous medium Highly asymmetrical output curves, where the tracer tails are noticeably flat, are often caused by the presence of an immobile fluid phase likely to exchange solute with the mobile phase (e.g., Villermaux, 1972; Porel, 1988). The exchange phenomenon generally produces a delay in the concentration peak as compared to the time of pure convection (Coats and Smith, 1964; James and Rubin, 1979; Zuber, 1986) and purely convective-dispersive models can then only fit an apparent dispersivity which varies with the transfer velocity and the transfer distance (Zuber, 1986). It would be presumptuous to argue that first-order kinetics, as used in our model, are capable of accounting for all the reversible exchange mechanisms between the mobile and the immobile phase or phases. However, in many of the complex and poorly understood cases, these kinetics constitute a simple and efficient means of integrating the noninstantaneous exchanges

F. Delay et aL / Journal of Contaminant Hydrology 25 (1997) 63-84

71

Table 2 Tracer tests and fitting parameters in homogeneous segments of Catsan pellets No.

x (cm)

Q (cm3/min)

Mo (Ixg)

ot (cm)

wc (%)

wt (%)

fl (min- 1)

C6 C5 C3 C4 C2 CI

54.5 26.5 26.5 54.5 54.5 26.5

46 46 24 24 16 16

1000 1000 1000 1000 1000 1000

0.65 0.65 0.65 0.65 0.68 0.7

47 47 47 47 47 47

72 70 70 72 72 72

0.055 0.052 0.040 0.035 0.022 0.023

into the transport equations (Gaudet et al., 1977; Legrand-Marc and Laudelout, 1985; Valocchi, 1986; Bahr, 1989; Maloszewsky and Zuber, 1991). The tracer experiments on the granular Catsan medium and the fittings on the particle model are shown in Table 2 and Fig. 3. The exchange phenomenon adds the immobile volumetric fraction L - wc and the exchange rate fl to the parameters ~ and w c of the a d v e c t i o n - d i s p e r s i o n phenomenon. The equilibrium constant Kp is set at l, assuming that the immobile phase is the same as the mobile one, i.e. water. The strength of this assumption was verified with batch experiments showing that the uranine was never sorbed on the solid (no decrease with time of the initial concentration C O in water in contact with Catsan pellets). A first analysis of the fittings (Table 2) leads to the following observations: wc and a are independent o f the distances and the transfer velocities; L is almost constant for all the tracer experiments and is in agreement with the measured value (73%) for the total porosity w t. One can therefore consider that the immobile volume fraction is constituted by the whole immobile fluid volume; - - As opposed to the immobile volumetric fraction, the pseudo-exchange rate /3 remains constant with the distance but increases with the mean pore velocity u (Table 2; Fig. 4). As stated above, /3 is the product of a mass transfer velocity and the a r e a / v o l u m e ratio of the solid. For a given mean pore velocity, /3 is constant with the distance because the a r e a / v o l u m e ratio is constant in a homogeneous colunm. However, the mass transfer velocity between immobile and mobile water is influenced by the mean pore velocity of the mobile water and results in the increase of /3 with u. As the exchange is likely to create a delay in the output curves, it is theoretically impossible to consider the modal time t m as equal to the pure convection time tc and to adjust first w c and then c~ as in a simple case of advection-dispersion. In practice, a close observation of the shape of the concentration increases in the output curves (examples, Fig. 3) shows that they do not vary very much. This means, in fact, that we admit that the exchange between the free and the bound phase only cuts the top off the concentration peak. The strength of this assumption was verified by a comparison between the values of ce and w~ obtained with the model and those that can be adjusted by fitting the first part o f the experimental curves on the analytical solution of a pure a d v e c t i v e - d i s p e r s i v e problem [Sauty's (1977) method of dimensionless tables]. More-

-

-

-

72

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84 C (~g/l)

-

-

-

5000 4000

]

....

l /

Simulated curve

J

Q= 16 cm3/min

J

~pe~e~tal

ooo. 2000

C2

1000

.

. )

100

120

0 0

20

40

60

C (gg/1)

80

- - - Experimental curve ] Simulated curve m l Q = 24 cm3/min

5000 4000

3000 2000

C4

1000

Time (min)

0

"

0

20

40

5000 C (gg/1)

,

60

80

- - - Experimental curve ] Simulated curve

l

Q = 46 cm3/min

4000

C5

30O0

1000

ox 10

20

30

40

Fig. 3. Fitting of the tracer experiments on the dual-porosity C a t s a n medium. O n l y the e x c h a n g e rate /3 varies with the flow rate. The other parameters are fairly constant with: dispersivity, ~ - 0 . 6 5 - 0 . 7 0 cm; kinematic porosity, wc = 0.47; a n d total porosity, L = w, = 0 . 7 0 - 0 . 7 2 . The tracer transfer distances are 26.5 c m for C I , C3 and C5, and 54.5 c m lbr C2, C4 and C6. All the injections are short point injections a n d the total cross-section of the c o l u m n is 19.635 c m 2.

13(min-1)

0,07 0,06 0,05 0,04 0,03 0,02

•

0,01

u (cm/min) 0

2

4

6

8

Fig. 4. Evolution o f the e x c h a n g e rate /3 of the first-order kinetics with the m e a n pore velocity u.

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84

73

over, we tested an automatic fitting of the experiments by minimizing a sum of squared differences as an objective function with a Broyden quasi-Newton method (in Dennis and More, 1977; Gill et al., 1991) and with a conjugate gradient method (Evans, 1967; Powell, 1984; Gilbert and Nocedal, 1992). The convergence down to a stable solution is obtained only for values of a and wc, that, give or take 15% for any experiment, remain constant at 0.65 cm and 47%, respectively. Note that the kinematic porosity value of 47% is almost the same as the experimental one measured by weighing. A dispersivity of 0.65 cm is a relatively high value compared to the mean grain size of the medium. Recall that in silica beads the dispersivity was about half the grain size, as usual for nonconsolidate media with a unique grain size. In Catsan, the grain size distribution is linear (see above) and the local heterogeneity of flow paths makes the apparent dispersivity to increase. Another explanation could be evoked and concerns the presence of intra granular porosity. The possible flow inside the grains, near the surface, may have a significant influence on dispersion. Given the relatively constant value of parameter L (Table 2) over all the experiments, the immobile volumetric fraction L - wc is constant in the fittings. This does not stem from a deliberate decision but from a characteristic imposed by the fitting of the

1

C/Cmax

a

" - " Experimental curve -Simulated curve

0,8

ct = 0.72 cm wt = 0.70 r 13a= i0.032 n - 1

0,6 0,4 0,2

_),J

0 0

1

,

20

40

60

80

C/Cmax

I It

0,6

I

0,4 0,2

'

~

~

Time (min)

60

"780

-)

o 0

20

40

0,8

. . -- - - -

: 100

~

0,2

Experimental curve Simulated curve

_

= 0.72 cm wt = 0.63 _ 13= 0.032 min-1

,

Time (min)

0

100

- - - Experimental curve I -- - -- Simulated curve I c~= 0.72 cm Ii wt = 0.63 ] 13= 0.022 min-1 I

t

~l

0,4

C I

0,8

b

C/Cmax

0,6

e (rain) ,

- - ~ - ~

1

0

1

20

40

60

~

.

-

-

-

Experimental curve ) - -

0,6

Saim=ul1.5~e'urve wt = 0.63 ] 13=°'°22"~i"_______A1 J

_ _

0,4

0

100

d

C/Cmax

0,8

0,2

80

Time (min)

3 20

40

60

80

100

Fig. 5. E x a m p l e o f the destabilizing o f a tracer e x p e r i m e n t fitting in a dual-porosity m e d i u m . E x p e r i m e n t on the Catsan m e d i u m over a distance o f 100 c m with a total cross-section o f 19.635 c m 2, flow rate o f the system 4 4 c m 3 / m i n , kinematic porosity 0.47. a. Best fit o f simulated c u r v e - e x p e r i m e n t a l curve. b. D e c r e a s e o f the tail o f the simulated output curve t h r o u g h the reduction o f the i m m o b i l e volumetric fraction. c. L e n g t h e n i n g o f the tail o f the output curve but n a r r o w i n g o f the concentration p e a k through a reduction o f the e x c h a n g e rate. d. S p r e a d i n g of the p e a k but shift o f the increase in concentration b e c a u s e o f an increase in dispersivity.

74

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84

experiments. Although the parameters L - w~ and /3 have similar effects on the kinetics (Eq. (5): if L or /3 increases, " t h e trapping" of solute in immobile water increases), their influences on the shape of the curves are different. Take the example of Fig. 5. A tracer experiment is correctly fitted (Fig. 5a) with w~ = 0.47, c~ = 0.72 cm, L = w t = 0.70 and /3 = 0.032 m i n - 1. If we destabilize the fitting by lowering w t to 0.63 (Fig. 5b), the simulated curve correctly adjusts the ascending and the first third of the descending portion of the experimental curve but the tail of the simulated output curve is much too short. While keeping w t = 0.63, one can try to lengthen the tail of the curve by reducing [3 in order to obtain slower diffusion out of the immobile water volume. This (Fig. 5c) again generates a flat tail but the peak is too narrow and shifts to the left as soon as the simulated curve starts to descend. In order to spread it out again, one has to increase the dispersivity (Fig. 5d), but the curve then reverts to type b (Fig. 5) besides having an added advance in time. This example shows how sensitive the model is to fluctuations in its parameters, thus providing evidence of the precision of the fittings.

6. Predicting tracer experiments in heterogeneous and mixed media The main objective is now to predict the transport in heterogeneous and mixed media on the basis of the different parameters defined for the two homogeneous media. Two sets of experiments were carried out. The first one deals with a heterogeneous medium in a vertical column containing alternating segments of silica beads and Catsan pellets disposed as shown in Fig. 6. The second set of experiments concerns a vertical column filled with homogeneous mixtures of silica beads and Catsan. 6.1. T r a c e r tests in a l a y e r e d h e t e r o g e n e o u s m e d i u m

Table 3 gives the tracer experiments and their parameters. First, a fitting was done by a classical trial-and-error method and by minimizing the sum of squared differences between simulations and experiments. Different parameters are fitted for each medium

Size(cm) ~ [ ' - - ' - 9 , 5

I

ll

Silicabeads

I ~

8

I

10

I

8

I

7---t

"Catsan"

Fig. 6. Column simulating a heterogeneous medium constructed by successive segments of glass beads and Catsan material. There are two injection points for transfer distances of 26.5 and 53.5 cm.

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84

75

Table 3 Tracer tests and fitting parameters in a heterogeneous medium composed of alternating segments of silica beads and Catsan pellets No.

SC1 SC2 SC3 SC4 SC5 SC6

x

Q

Mo

a (cm)

(cm)

(cm3/min)

(p~g)

silica

Catsan

silica

Catsan

26.5 53.5 26.5 53.5 26.5 53.5

20 20 46 46 70 70

500 500 500 500 500 500

0.14 0.14 0.14 0.14 0.14 0.14

0.7 0.7 0.7 0.68 0.7 0.7

40 40 40 40 40 40

47 47 47 47 47 47

wc (%)

(wt%)

/3 Catsan

(min 1) 72 72 72 72 70 70

0.034 0.034 0.052 0.052 0.063 0.062

(silica beads or Catsan). The fitting parameters were almost the same as their counterparts in the homogeneous medium: in the silica beads: constant dispersivity a = 0.14 cm, wc = 40%, no dual porosity, therefore no exchange; in Catsan: constant dispersivity a = 0.7 cm, wc = 46%, immobile volumetric fraction around 0.25, exchange rate varying from 0.06 rain-1 for u = 7.6 cm min-1 to 0.03 min -1 for u = 2.15 cm min -1. Then, simulated curves were computed without fitting by introducing into the model the same series of segments as in the true column and for each segment, a parameter set identical to that used when fitting the experiments with homogeneous media. If a comparison is made between the experimental and simulated curves for different flow rates and two transfer distances (26.5 and 53.5 cm, respectively, i.e. for different injection points; Fig. 6), the results are excellent (Fig. 7). The good agreement between simulations and experiments proves that this heterogeneous layered medium keeps the parameters of each of the constituent media. Note that, at least for the studied ranges, the variations in transfer distance and flow rates do not influence the values of these parameters. It must be emphasized that these results are only valid for conservative (i.e. non reactive) solutes. With radioactive decay or reactive processes between solutes and solid (Kp 4~ 1), the mass restituted and the resting time distribution depend on the mean transfer time and, of course, on the layer distribution along the column. This could therefore modify the transfer parameters of a layer according to the position in the column. In fact, at the scale of the global transport problem through the medium, there is no evidence of non-local dispersion. The non-local form of the dispersive flux term that was suggested for example by Koch and Brady (1987, 1988), Naff (1990), Cushman and Ginn (1993), Neuman (1993), and Sternberg and Greenkorn (1994) has not been observed here. This is probably due to the small difference in permeability between the Catsan and the silica beads (1.5.10 -2 and 1.2- 10 2 m s-1, respectively, for experimental measures in 20-cm-long columns). Of several suggested dispersion mechanisms (Greenkorn, 1983), it seems that the difference between dual-porosity and single-porosity dispersion does not generate non-local phenomena. In other words, for a 1-D system such as this, the transport problem can be -

-

-

-

76

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84 - - - Experimental. . . . 1 - -Simulatedcurve [ Q = 20 cm3/min [

C (p-g/l) 3000 2500

~

SC1

1500 -

SC2

1000

0

20

40

60

Experimental curve

3500 C (p-g/l) 3000 •

80

SC3

Q = 46 cm3/min

2500 20O0 1500

~

_

1000

500 0

Time (min) - -i 30 40

• 0

10

20

- - - Experimental curve] Simulated curve j Q = 70 cm3/min J

C (p,g/l)

35OO . 3000

~

2500

SC5

tl s

1500 1000 500

C6

~ 0

5

• 10

15

20

' 25

Fig. 7. Fitting of tracer experiments in a heterogeneous medium by global m o d e l l i n g ( o n l y one column

simulated). The fitting parameters are those of each one of the homogeneous media, i.e. a = 0.14 c m 1or the silica beads, a = 0.7 c m , w, = 0.7, /3 = 0.06 m i n i to 0.03 rain 1 b e t w e e n u = 7.6 cm rain 1 and u = 2.15 cm m i n - J for Catsan. The kinematic porosity is 0.40 for the silica beads and 0.47 for Catsan.

considered as the convolution of elementary homogeneous systems where the breakthrough curve at the outlet o f a system becomes the input of the next one. Take the example of modeling with a system in series and with the tracer experiment SC1 in Table 3. It was simulated numerically over a distance o f 53.5 c m with a f l o w rate of 46 c m 3 min -1 on 6 successive h o m o g e n e o u s columns, numbered from 1 to 6 in the input-output direction and with the m e d i u m structure shown in Fig. 6: - - columns 1, 3 and 5 with respective lengths of 9.5, 8 and 8 c m are given parameters specific to the silica beads, i.e. w c = 0.4, c~ = 0.14 cm, no i m m o b i l e water volume; - - columns 2, 4 and 6 with respective lengths of 11, 10 and 7 c m are given the parameters o f the Catsan medium, i.e. wc = 0.47, a = 0.7, w t = 0.72, /3 = 0.052 rain- 1, Kn=l.

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84

Silica beads

O1 12

~

O 2 I3

77

"Catsan"

0 3 14

O 4 15

0 5 I6

c ~g~) 14000 -

-

-

12000 - -

O1

.

Experimental curve Simulated curve Q = 46 cm3/min

10000 - -

8000 - i

6000-O2

03

4000 - -

2000.-

0

-J-0

2

4

6

8

10

12

14

16

18

20

F i g . 8. E x a m p l e o f the f i t t i n g o f a t r a c e r e x p e r i m e n t in a h e t e r o g e n e o u s m e d i u m b y a m o d e l w i t h c o l u m n s in series. F o r e a c h c o l u m n , the f i t t i n g p a r a m e t e r s are t h o s e o f the h o m o g e n e o u s m e d i u m : ot = 0 . 1 4 c m , w c = 0.4 in the g l a s s b e a d s , c~ = 0.7 c m , w c = 0.47, w t = 0 . 7 2 a n d /3 = 0 . 0 5 2 m i n -1 in the C a t s a n . T h e o u t p u t f u n c t i o n o f o n e c o l u m n b e c o m e s the i n p u t f u n c t i o n o f the n e x t o n e , e.g., 0 2 = I3. I n this e x p e r i m e n t , 0 6 c o r r e s p o n d s to the o u t p u t o f the s y s t e m a n d p e r f e c t l y fits t h e t r a c e r e x p e r i m e n t ( S C 4 , T a b l e 3).

The output function of the preceding column becomes the input function of the next one. The 6 output functions, numbered from O1 to 06, are shown in Fig. 8. It is thus possible to follow the evolution of the shape of the output curves as the tracer advances in the heterogeneous medium. As expected, one observes that, apart from a shift in time, the curves are only slightly deformed when the tracer passes through the silica beads (compare, e.g., 0 2 and 0 3 ; Fig. 8). 0 6 , which represents the output of the system, perfectly fits the experimental curve. It seems to us that these experimental results argue for the use of a geostatistical simulation technique to generate heterogeneous reservoirs when mass transfer simulations are necessary. The principle is the following: first, a spatial distribution of typical geological facies is calculated on a regular meshed domain; then, for each block, hydrodynamic and hydrodispersive properties are allocated according to the facies

78

F. Delay et a L / Journal of Contaminant Hydrology 25 (1997) 63-84

present in the block (de Marsily, 1993). Our 1-D experiments show that this principle is satisfactory because transport over the whole domain is, in effect, the convolution of a series of elementary transports in successive homogeneous blocks. In this sense, simulating a heterogeneous reservoir as the juxtaposition of elementary homogeneous blocks does not alter the nature of the transport process. It must be pointed out that our experiments deal with a flow perpendicular to the layers and the macrodispersion of the tracer would be different with a flow parallel to the layers. However, the convolution principle from one block to the other remains valid. In natural cases where the flow can, a priori, occur in all directions, the convolution principle is applied to each flow line, thus the global accuracy of the model depends on the velocity field resolution and on the size o f the blocks discretizing the reservoir. 6.2. Tracer tests in mixed media

Take the example of an alluvial aquifer consisting only of sand and clay. In fact, we know that the aquifer is not a layered system with areas of pure clay and areas of pure sand but rather a system of blocks where on some cases, clay predominates, and on others, sand predominates. Even if, at the scale of the block, there are data for the pure facies, how do we introduce this information into blocks where we know that the medium is mixed? The fourth series of experiments illustrates this problem and, although it does not allow us to propose a global solution, it suggests some general ideas that can be used to improve the definition of the parameters in each block. As mentioned above, three homogeneous mixtures of Catsan and silica beads were used with volumetric proportions of: 25% C a t s a n - 7 5 % silica beads, 50% C a t s a n - 5 0 % silica beads and 75% C a t s a n - 2 5 % silica beads. Henceforth, they will be noted $75C25, $50C50 and 825C75 , respectively. The tracer tests were carried out in vertical columns for a transfer distance of 53.5 cm and two flow rates: 20 and 40 cm 3 min -1 . The model of particle tracking including exchange between immobile and mobile water (see Section 2.2) was used to fit the experimental breakthrough curves (Fig. 9). Except for the exchange rate /3, the parameters, i.e. dispersivity, kinematic porosity and immobile fraction, remain the same whatever the flow rate, but vary from one mixture to another (Table 4): - - F o r the mixture 875C25 , the kinematic porosity wc is 35.5%, the total porosity w t is 42%, giving an immobile fraction of 6.5% and a dispersivity a of 0.4 cm; - - for the mixture $50C50, wc = 33%, w~ = 43% (immobile fraction = 10%) and a = 0.52 cm; - - for the mixture $25C75, wc = 47%, w t = 62% (immobile fraction = 15%) and a = 0.64 cm. The immobile fraction w t - w e can be considered as an approximately linear combination of the immobile fraction of each component (0% for silica beads, 2 2 - 2 5 % for Catsan) weighted by its volumetric fraction in the mixture. This result appears logical and shows that, a priori, the access to the microporosity of the Catsan pellets is not disturbed by the presence of silica beads. As to the dispersivity, although it increases with the fraction of Catsan and remains bounded between 0.13 cm (silica beads) and 0.7 cm (Catsan), it is not a simple linear

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84 1200 C ~g/l) MSC1 - Q = 42 1000800

~

- - . ~

10

Experimental curve Simulated curve

Mix : $75C25

20

800 - (2' (gg/1) 700 MSC3 - Q = 37 6005004OO

79

30

40

50

- - . Experimental curve Simulated curve Mix : $25C75

20

-

300200-

1®

i__j

0 .... 10 900

20

C ~g!l) MSC5

800 700

- Q - - 48

30

40

50

- - . Experimental curve ---Simulated curve [ _ _ Mix : $50C50

500 400. 300. 200

~

10o

0

10

~ 20

T ~ (mi~

3o

40

5O

Fig. 9. Fitting o f tracer experiments on mixtures o f silica b e a d s a n d C a t s a n pellets. All the injections are short point injections, the total cross-section o f the c o l u m n is 19.635 c m 2 and the transfer distance is 53.5 cm. The volumetric proportions o f silica b e a d s a n d C a t s a n pellets are 7 5 % : 2 5 % for M S C 1 a n d M S C 2 , 2 5 % : 7 5 % for M S C 3 and M S C 4 , a n d 5 0 % : 5 0 % for M S C 5 a n d M S C 6 . The f l o w rates given in the figure are expressed in cm 3 min- 1

Table 4 Tracer tests a n d fitting p a r a m e t e r s in h o m o g e n e o u s mixes c o m p o s e d of silica beads and C a t s a n pellets in different volumetric proportions No

MSCI MSC2 MSC3 MSC4 MSC5 MSC6

V o l u m e t r i c fraction

x

Q

Mo

o/

wc

Catsan

silica

25 25 75 75 50 50

75 75 25 25 50 50

wt

fl

(cm)

(cm3/min)

(ixg)

(cm)

54.5 54.5 54.5 54.5 54.5 54.5

42 20 37 20 48 20

175 175 200 200 175 175

0.40 0.35 0.67 0.61 0.52 0.56

(%)

(%)

( m i n - 1)

35.4 35.6 46.5 46.8 32.0 33.0

42.0 42.0 62.0 65.0 42.0 45.0

0.02 0.012 0.026 0.02 0.032 0.017

80

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84

0,7

Dispersivity (cm)

•

0,6 0,5 0,4

L

0,3 + 0,2

Volumetric fraction of Catsan (%)

0,1 0 0

20

40

60

80

100

Fig. 10. Evolution of the dispersivity a (cm) vs. the volumetric fraction of Catsan in a homogeneous mixture with silica beads.

combination of the values of the two pure components. The plot of c~ vs. the volumetric fraction of Catsan (Fig. 10) shows a logarithmic increase of a with the evolution of the Catsan fraction. The sharp increase near the origin may be interpreted as follows: if a small quantity of Catsan pellets is introduced among the silica beads, it generates, at the grain size scale, local heterogeneities in the geometric distribution of the flow paths in space. If one admits that the dispersivity is related to the tortuosity (Bear, 1979; de Marsily, 1986), one understands that a small fraction of Catsan can induce a strong increase in the dispersivity values. Beyond a 25% volumetric fraction of Catsan, the evolution of the dispersivity value is quasi-linear. Concerning the kinematic porosity of the different mixtures, there is no rule allowing their values to be predicted from those of the pure components. For two mixtures, $75Cz5 and $50C50 (Table 4), the values are lower than 40% which corresponds to pure silica beads, and for the 825C75 mixture the value of 47% is equal to the value of pure Catsan. W e therefore repeated the experiments several times with the three mixtures, emptying and filling the columns between each experiment. W e found that wc varied between 34% and 38% for $7sC25, between 32% and 38% for Ss0C50, and between 45% and 47% for $25Cv5 whereas the other parameters: /3, w t - w~ and a remained constant. In fact, only one valid conclusion could be drawn: wo varies between the theoretical values of the kinematic porosity for homogeneous granular media, 25.9% and 47.6% for tetrahedric and cubic stackings, respectively. If the exchange rate /3 is plotted vs. the mean pore velocity u (Fig. 11), /3 increases with u whatever the mixture but is always lower than the value obtained for the same velocity u in pure Catsan (Fig. 4). Like with pure Catsan, the evolution of /3 with u corresponds to the increase of the mass transfer coefficient [L T - 1] between the free and immobile water at the grain surface. The ratio s u r f a c e / v o l u m e of the granular medium was calculated as the linear combination of the measured value in Catsan (BET measure giving: 1.1 • 106 c m - 1 ) and the theoretical value in silica beads (4"rrR2/4"rrR3 with R = 0.1 cm, giving 30 cm ~) each value being weighted by the volumetric fraction of the component in the mixture. Then if the exchange rate /3 is divided by the s u r f a c e / v o l u m e ratio, one obtains an approximation of the mass transfer coefficient. The latter is plotted vs. the mean pore velocity u in Fig. 12 and results in a quasi-linear

F. Delay et al./ Journal of Contaminant Hydrology 25 (1997) 63-84 0,035

13(min-l)

81

--

0,030

.--"E3

0,025

*

0,020

$25C75

0,015

.--

0,010

.

..-"

. -"" .-

.-

+

-"" [] - " s50c50 .-:~

.--""

$75C25

0,005

u (cm/min)

0,000 0

,

,

,

2

4

6

- - , 8

Fig. 11. Evolution of the exchange rate /3 of the first-order kinetics with the mean pore velocity u and for different homogeneous mixtures of silica beads and Catsan pellets.

correlation. Furthermore, for a given value of u, the greater the fraction of silica beads, the more fl decreases (Fig. 11). This result agrees with the linear decrease of the surface/volume ratio in the mixture when silica beads are added to Catsan. Based on these results we offer the following comment: whereas intrinsically, the access to the immobile water (matrix porosity in Catsan) is not disturbed by the silica beads, i.e. the immobile fraction of a mixture is the linear combination of the immobile fraction of each component weighted by its volumetric fraction, and the mass transfer coefficient in the solid evolves linearly with the mean pore velocity, the rate of exchange decreases linearly as soon as silica beads are introduced among the Catsan pellets. This decrease is due to an approximately linear decrease of the surface/volume ratio of the solid medium. These results appear logical and confirm that fl is effectively the product of a mass transfer coefficient [L T -1] and a surface/volume ratio [L -l ] of the granular medium.

13 / Ss (10-8 cm/min)

8 7 6 5 4

•

3

•

•

2

u (cm/min)

1 0 0

2

4

I

I

6

8

Fig. 12. Variation of the mass transfer rate in the solid with the mean pore velocity. The mass transfer rate is calculated as the (exchange rate)/(specific surface) ratio, the specific surface of the granular medium being the linear combination of the surface/volume ratio of each pure component weigthed by its volumetric fraction in the mixture.

82

F. Delay et al. / Journal of Contaminant Hydrology 25 (1997) 63-84

7. Conclusions By comparing the model with the results of tracer experiments we were able to demonstrate that it is possible, in homogeneous media, to fit the experiments with the same limited parameter set, for any hydrodynamic conditions and transport distances. In a single-porosity medium, two parameters are fitted: the kinematic porosity, w~ and the dispersivity of the medium, a. In the presence of dual-porosity, four parameters are used: the kinematic porosity w~, the dispersivity of the medium c~, the volumetric fraction of immobile water L - w~, the exchange rate of first-order kinetics /3. Only the latter changes with the transport velocity. The tests made on tracer experiments with layered heterogeneous media show that the transport phenomenon of non-reactive solutes can easily be represented by a convoluted series of elementary transport problems. The breakthrough curve at the outlet of an elementary homogeneous system becomes the input of the following one. It can also be concluded that the difference in dispersion between dual and single porosity does, apparently, not generate non-local effects on the dispersion. The transport problem will then correctly be solved at the scale of a reservoir if, for each block of equivalent homogeneous medium within the reservoir, the parameters used in the transport model are those determined for that medium independently. Consequently, heterogeneous reservoirs created by geostatistical numerical simulations can be used in a straightforward manner in transport calculations. It must be noted that with reactive solutes the convolution principle should be re-examined. In fact, it would no longer be suitable for predicting transport, if the parameters that describe the reactions are dependent of the mass present in the medium or of the resting time distribution. By mixing two pure components in a column we tried to use experimental results to find some practical rules for improving the definition of hydrodispersive parameters for "homogeneous mixtures". Three results are important: the immobile fraction of a mixture is a linear combination of those observed for the pure components; the dispersivity of a mixture presents a strong modification whenever a "foreign b o d y " is introduced into a pure component and then evolves linearly; finally, the mixtures tend to reduce linearly the exchange rate of a first-order kinetics according to the decrease of the surface/volume ratio of the solid medium. This ratio corresponds approximately to the linear combination of the surface/volume value of each component weighted by its volumetric fraction in the mixture.

Acknowledgements The authors would like to thank M.L. Brusseau for helpful comments and discussions concerning an earlier draft of this paper.

References Ackerer, P., 1988. Random-walk method to simulate pollutant transport in alluvial aquifers or fractured rocks. In: Custodio et al. (Editors), Groundwater Flow and Quality Modelling. Reidel, Dordrecht, pp. 475-486.

F. Delay et aL / Journal of Contaminant Hydrology 25 (1997) 63-84

83

Ackerer, P. and Kinzelbach, W., 1985. Modrlisation du transport de contaminant par la m&hode de marche au hasard, influence des variations du champ d'rcoulement au cours du temps sur la dispersion. Proc. Int. Syrup. on Stochastic Approach of Sub-surface Flow, Montvillagene, pp. 446-458. Bahr, J.M., 1989. Analysis of nonequilibrium desorption of volatile organics during a field test of decontamination. J. Contam. Hydrol., 4: 205-222. Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York, NY, 567 pp. Bobba, A.G., 1989. Numerical model of contaminant transport through conduit-porous matrix systems. Math. Geol., 21: 861-890. Coats, K.H. and Smith, B.D., 1964. Dead-end pore volume and dispersion in porous media. Soc. Pet. Eng. J., 4: 78-84. Cushman, J.H. and Ginn, T.R., 1993. Non-local dispersion in media with continuously evolving scales of heterogeneity. Transp. Porous Media, 13: 123-138. Delay, F., Dzikowsky, M. and de Marsily, G., 1993. A new algorithm for representing transport in porous media in one dimension, including convection, dispersion and interaction with the immobile phase with first-order kinetics. Math. Geol., 25: 689-712. Delay, F., de Marsily, G. and Carlier, E., 1994. One-dimensional solution of the transport equation in porous media in transient state by a new numerical method for the management of particle track. Comput. Geosci., 20: 1169-1200. Delay, F., Housset-Resche, H. and Porel, G., 1995. Simulation numrrique de r~servoirs heterog~nes: apport dans le transfert de masse en milieu saturC Bull. Soc. Grol. Fr., 166: 627-637. Delay, F., Housset-Resche, H., Porel, G. and de Marsily, G., 1996. A new method of particle tracking to solve the transport equation in saturated porous media and in two dimensions. Math. Geol., 28: 45-71. Delhomme, J.F., 1979. Spatial variability and uncertainty in groundwater flow parameters: a geostatistical approach. Water Resour. Res., 15: 269-280. de Marsily, G., 1986. Quantitative Hydrogeology: Groundwater Hydrology for Engineers. Academic Press, San Diego, CA, 440 pp. de Marsily, G., 1993. Quelques mrthodes d'approche de la variabilit6 spatiale des rrservoirs souterrains. Hydrogrologie, 4: 259-268. Dennis, J.E. and More, J.J., 1977. Quasi-Newton methods, motivation and theory. SIAM Rev., 19: 46-89. Evans, J.D., 1967. The use of pre-conditioning in iterative methods for solving linear equations with symmetric positive definite matrices. J. Inst. Math. Appl., 4: 295-314. Gaudet, J.P., Jegat, H., Vachaux, G. and Wierenga, P.J., 1977. Solute transfer with exchange between mobile and stagnant water through unsaturated sand soil. Sci. Soc. Am. J., 41: 665-671. Gilbert, J.C. and Nocedal, J., 1992. Global convergence properties of conjugate gradient methods for optimization. SIAM J. Opt., 2: 21-42. Gill, P.E., Murray, W. and Wright, M.H., 1991. Numerical Linear Algebra and Optimization. Addison-Wesley, Redwood, CA, 423 pp. Gomez-Hernandez, J.J. and Journel, A.G., 1992. Joint sequential simulation of multigaussian fields. In: A. Soares (Editor), Geostatistics Troia. Kluwer, Dordrecht, pp. 85-94. Gomez-Hernandez, J.J. and Srivastava, R.M., 1990. ISIM3D: an ansi C three dimensional multiple indicator conditional simulation program. Comput. Geosci., 16: 395-440. Greenkorn, R.A., 1983. Flow Phenomena in Porous Media. Marcel Dekker, New York, NY, 550 pp. James, R.V. and Rubin J., 1979. Applicability of the local equilibrium assumption to transport through soils of solutions affected by ion exchange. In: E.A. Jenne (Editor), Chemical Modelling in Aqueous Systems. Am. Chem. Soc., Washington, DC, pp. 225-235. Journel, A.G. and Huijbregts, C.J., 1978. Mining Geostatistics. Academic Press, London, 600 pp. Kinzelbach, W., 1986. Groundwater Modelling, An Introduction with Sample Programs in Basic. Elsevier, Amsterdam, 333 pp. Koch, D.L. and Brady, J.F., 1987. A non-local description of advection-diffusion with application to dispersion in porous media. J. Fluid Mech., 130: 387-403. Koch, D.L. and Brady, J.F., 1988. Anomalous diffusion in heterogeneous porous media. Phys. Fluids, 31: 965-967. Legrand-Marc, C. and Laudelout, H., 1985. Longitudinal dispersion in a forest stream. J. Hydrol., 78: 317-324.

84

F. Delay et al./ Journal of Contaminant 14ydrology 25 (1997) 63-84

Maloszewsky, P. and Zuber, A., 1991. Influence of matrix diffusion and exchange reactions on radiocarbon ages in fissured carbonate aquifers. Water Resour. Res., 27: 1937-1945. Naff, R.L., 1990. On the nature of the dispersive flux in saturated heterogeneous porous media. Water Resour. Res., 26: 1013-1026. Neuman, S.P., 1993. Enlerian-Lagrangian theory of transport in space-time non stationary fields: exact non-local formalism by conditional moments and weak approximation. Water Resour. Res., 29: 633-645. Porel, G., 1988. Transfert de solute en aquifrre crayeux: causes de modifications des rEsultats de tra~ages. Th~se Doctorat Univ., UniversitE de Lille, Lille, 327 pp. Powell, M.J.D., 1984. Nonconvex minimization calculations and the conjugate gradient method. Lect. Notes Math., 1066: 122-141. Putti, M., Yeh, W.G. and Mulder, W.A., 1990. A triangular finite approach with high-resolution upwind terms for the solution of groundwater transport equations. Water Resour. Res., 26: 2865-2880. Sauty, J.P., 1977. Contribution h l'identification des param~tres de dispersion dans les aquif'eres par interpretation des experiences de tra~ages. Th~se Doct. lng., UniversitE de Grenoble, Grenoble, 157 pp. Sternberg, S.P.K. and Greenkorn, R.A., 1994. An experimental investigation of dispersion in layered porous media. Transp. Porous Media, 15: 10-30. Uffink, G.J.M., 1985. A random-walk method for simulation of macrodispersion in a stratified aquifer. In: I.U.G.G. 18th General Assembly Proceedings of the Hamburg Symposium. Int. Assoc. Hydrogeol. Sci. Publ., 146: 103-114. Uffink, G.J.M., 1990. Analysis of dispersion by random walk method. Ph.D. Thesis, Delft University, Delft, 150 pp. Valocchi, A.J., 1986. Effect of radial flow on deviations from local equilibrium during sorbing solute transport through homogeneous soils. Water Resour. Res., 22: 1693-1701. Van Leer, B., 1979. Towards ultimate conservative scheme: a second order sequel to Godunov's method. J. Comput. Phys., 32: 101-136. Villermaux, J., 1972. Analyse des processus chromatographiques linEaires ~t l'aide de modbtes phEnomEnologiques. Chem. Eng. Sci., 27: 1231-1243. Zuber, A., 1986. On the interpretation of tracer data in variable flow system. J. Hydrol., 86: 45-57.