Proposal for a probabilistic model of dispersion: A first validation

Available online at www.sciencedirect.com

SCIENCE ~__~DIRECTe ELSEVIER

MATHEMATICAL AND COMPUTER MODELLING

Mathematical and Computer Modelling 42 (2005) 1137-1144 www.elsevier.eom/locate/mcm

P r o p o s a l for a P r o b a b i l i s t i c M o d e l of Dispersion: A First Validation E. CARLIER* AND J. EL KHATTABI C e n t r e de Calcul et de Mod41isation de Lens Facult4 J e a n Perrin, Universit4 d ' A r t o i s S P 18, rue J e a n Souvraz, 62 307 Lens cedex, F r a n c e c a r l i e r ~ u n i v - a r t ois. fr

(Received and accepted May 2004) Abstract--The probabilistic approach is used to simulate the particle tracking for two types of porous media. The first one is sand grains with a single intergranular porosity. The particle tracking is carried out by advection and dispersion. The second one is chalk granulates with intergranular and matrix porosities. Sorption can occur with advection and dispersion during particle tracking. The particle tracking is simulated as the sum of elementary steps with independent random variables in the sand medium. An exponential distribution is obtained for each elementary step and shows that the whole process is Markovian. A gamma distribution or probability density function is then deduced. The relationship between dispersivity and the elementary step is given using the central limit theorem. The particle tracking in the chalky medium is a non-Markovian process. The probability density function depends of a power to the distance. Experimental simulation by dye tracer tests on a column have been performed for different distances and discharges. The probabilistic approach computations are in good agreement with the experimental data. The probabilistie computation seems an interesting and complementary approach to simulate transfer phenomena in porous media with respect to the traditional numerical methods. @ 2005 Elsevier Ltd. All rights reserved.

Keywords--Particle,

Tracking, Groundwater, Dispersion, Modelling.

1. I N T R O D U C T I O N The probabilistic approach of particle tracking in natural environment is not new. It was used successfully in sedimentary transport [1,2] in solutes transfer in porous medium without sorption phenomena [3,4] and in solutes transfer in fissured medium [5-7]. The modelling of particle tracking in a porous or fissured medium with sorption capacity is very difficult when the exchange coefficients between matrix and water are difficult to estimate. The modelling of such a process is carried out by the resolution of the dispersion equation by taking account of the sorption phenomena [8 17]. According to these authors, the models are very sensitive to the exchange coefficient values. The error on these coefficients must be very weak in order to keep the model reliable. Furthermore, particle tracking modelling using the classic solution of the *Author to whom all correspondence should be addressed. 0895-7177/05/$ - see front matter @ 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2004.05.014

Typeset by .AjMS-TEX

1138

E. CARLIER AND J. ELKHATTABI

dispersion equation leads to an increase in apparent dispersivity with the distance if adsorption phenomena occur but are not taken into account [18]. The aim of this paper is to derive the dispersion equations in a porous medium using the probabilistic approach. The starting point is to simulate the particle tracking in a single porosity medium without sorption. The particle tracking is considered as the result of the succession of independent random variables according to a Markovian process. The second point is to derive the dispersion equations in a medium with dual porosity by taking account of the former particle tracking. 2. T H E O R E T I C A L

ASPECT

According to the theorem of conditional probability, the probability of a particle to be in the interval {x, x + Ax} is

(

P(x
X<_x+-~->x

P(X>x).

(1)

If f is the probability density function (PDF) of the random variable x and F the corresponding cumulative distribution function (CDF), (1) can be expressed by

ff

+Ax

(2)

= F ( x + zXx) - F ( x )

and

P ( X > x) =

~x+°z f(o~) do~ = 1 - ~0x f(c~) dx = 1 -

F(x).

(3)

The properties of the solid medium and particularly its sorption capacity must be known to estimate P ( X <_ x + A x / X > x). The elementary steps of distance carried out by the particles can be regarded as independent random variables if the solid phase is nonadsorbent. Then P ( X < x + A x / X > x) must be proportional to the length of Ax:

P

(

X<_x+--ff->x

=bAx.

(4)

If sorption occurs in the medium with a random distribution of the sorption site, P ( X <_ x + A x / X > x) must depend on Ax but also on the distance X. Indeed, the probability of a particle to be in the interval {x, x + Ax} is stronger when X is larger, i.e., the probability for a particle to meet a site of sorption increases when the distance covered by the particle is longer:

P

(

Xx

= b A x x a.

(5)

Equation (4) was tested on siliceous material with a single intergranular porosity and equation (5) was tested on chalk granules medium with a dual porosity. 3. M E D I U M

WITHOUT

SITES

OF

SORPTION

By combining equations (1)-(4), the probability functions f and F can be expressed by

F ( z ) = l - e -bx, f ( x ) = b e -bx.

(6)

Equation (6) evidences the Markovian character of variable X. With hypothesis (4), result (6) follows straightforwardly; its derivation can be found in a classical textbook like [19]. The last state of the particle does not have any influence on its future state. If the total distance X covered

P r o p o s a l for a P r o b a b i l i s t i c M o d e l

1139

by a particle is the sum of n elementary steps x~, then X is the sum of N independent random variables. Its density of probability is the product of convolution of order N of the probability density function of the variable xi. Let X = X1 + X2, where X1 and X2 are two independent random variables (two elementary steps), the probability density function of X is

fx(x)=

fx2(x

x l ) f x ~ ( x l ) d x l = b 2 x e -b~.

~0a~

For three variates X1, X2, and X 3 with Y = X 1 q- X 2 , x

fx(x)

b2

j£o

fx~(x - y)fy(y) dy

b-~.x 2 e -b~.

For n variables, fx(

) = be

' (7)

fro+~° f x (x) dx = 1. F is the g a m m a function. Equation (7) is a probability density function, called the g a m m a distribution. The mean distance, variance and modal distance (distance to peak concentration at a fixed time) are n

2

cr

n

b2 ,

n--1 Xm°d -b

(8) '

where 1/b is the mean distance of an elementary step. If one assumes the immediate injection of a mass M of particles, the mass conservation equation gives

fro °° C(x, t ) S aodx = M,

(9)

where S (L 2) is the cross area of the medium and co the kinematic porosity. C(x,t) (ML -3) is the concentration at distance x and time t. The flow Q (L3T -1) is the product of the average velocity U (LT-1), the cross area and the kinematic porosity. The following equation checks (9): M

C(x, t) = -~ U b e -bx

(bx) n - i

(10)

Equation (10) evidences the distribution of concentration in space is not Gaussian. But, according to the central limit theorem [19,20], it is possible to relate the p a r a m e t e r b and the dispersivity a ( L ) [21]. The distance X is the sum of n independent r a n d o m variables of mean value 1/b and variance 1/b 2. According to the central limit theorem, the variate ( X - n/b)/v/7-/b 2 has a limiting cumulative distribution function which approaches a normal distribution with a zero m e a n and a variance equal to one. T h a t means that the variable X is normally distributed with a mean equal to n/b and a variance equal to n/b 2. The solution of the equation of dispersion, for an impulse input in one dimension, is C ( x , t ) - 5I U e_(x_Ut)2/4Dt ' (11)

Q 417-

1140

E. CARLIER AND J . E L K H A T T A B I

where D (L2T -1) is the dispersion coefficient and is equal to the product of the average velocity U and the dispersivity c~. According to (11), the concentration distribution in space is Gaussian with a mean equal to Ut and a variance equal to 2Dt. By analogy, one deduces Tt

- = Ut, b

(12)

n

2Dt = ~.

Replacing n by bUt and D by c~U in (12), the relationship between b and the dispersivity is 1 c~ = 2b"

(13)

If n is large, the average elementary step is equal to the double of dispersivity. According to (12), at time t, the number of steps carried out by the particle is equal to but. Equation (10) becomes M (bx) but-1 C ( z , t ) = - ~ U b e -b~ r(bUt) " (14)

4. M E D I U M

WITH

SITES

OF SORPTION

By combining equations (1)-(3) and (5), the functions f and F can be expressed by F ( x ) = 1 - e -b(z~+~)/(l+a),

(15) f ( x ) = bx ~ e -b(xl+~)/(l+~).

The random variable X does not present a Markovian character:

( X > x + ~ - Ax > x ) P ( X > x )

P(X>xAX>x+Ax)=P

=P(X>x+Ax).

If the cumulative distribution function is F ( x ) = 1 - e -b(~l+~)/(l+a), then Ax

P X>z+--X->Z

)

=

P(X

> x + A x ) = e_(b/(l+a))((x+Ax)l+~_xl+~),

P(X>z)

the future state of the particle depends on the former state, the random variates xi are dependent. The function f checks

~

+cc f ( x ) dx = 1,

Va.

(16)

According to the same procedure used to derive equation (10), concentration is C ( x , t) = @ Ubxae -b(xl+~)/(l+a).

(lr)

By solving the equation oc = 0, the modal distance for which the concentration is maximum is Xmod =

(b) 1/(l+a)

(18)

k bl/(l+a-------~

(19)

The distance Ut is unknown but is given by Ut-

Proposal for a Probabilistic Model

1141

Equation (19) allows us to compute the b value as a function of time:

b=

(20)

The k coefficient must be given according to the mass law conservation, by considering that a particle adsorbed by the matrix will be then restored with the fluid when the gradient of concentration between interstitial water and water in the matrix is reversed:

fo +~ C(z,t)Qdt

= M.

(21)

By replacing b in (20) and (17), the concentration equation is

C(2g, t) = M

( k ~ 1-t-a

(22)

While solving -~oc = 0, the average velocity U and the m a x i m u m concentration can be determined according to k, the modal time tin, the distance X , and the power number a:

kz U = t,~(1 + a)U(l+a) ' Cm~x = k ( l +

(23)

a) a/(l+a) Q eMl t r n

5. E X P E R I M E N T A L

•

(24)

METHOD

The column used for the tracer tests has a length of 1934 m m and a diameter of 90 ram. Two piezometers were installed in the column; the first piezometer was at a distance of 143 m m from the inlet, the second at a distance of 825 m m from the first one. The tracer was injected with a syringe in the first piezometer; the injection was followed by a flush of water in order to approximate a delta function input. The tracer tests were carried out using sand from the Seine River and in Senonian Chalk. The size of grains was between 2 and 3.15mm. A mass of 204kg of sand and of 10.7kg for chalk filled a volume of 12.21 in the colmnn. Using a value of 2 6 5 0 k g / m 3 for the density of sand and 1590 k g / m a for the density of chalk, the total porosity can be inferred. The inferred total porosity for sand is 37.2% and 45% for chalk.

6. D Y E

TRACING

IN

SAND

T h e t r a c e r t e s t s w e r e c a r r i e d o u t for t w o d i f f e r e n t d i s c h a r g e s a n d t w o m a s s e s o f f l u o r o s c e i n

(Table 1). The b p a r a m e t e r was determined using the concentration m a x i m u m Cmax

M

z

- -Q t-~

be_bx (bcc)b~-I C(bz)

- 0.

(25)

Table 1. Experimental and computed data for dye tracing in sand. Distance (mm)

Q (mrn3/min)

M (ing)

Concentration Max (mg/1)

t,~ (min)

c~(mm)

1 ~ (mm)

825

102600

5.0542

7.62

22.5

5.419

5.407

1791

102600

5.0542

5.6

43

5.964

5.958

825

355030

3.1322

5.02

6

5.632

5.619

1791

355030

3.1322

3.63

12

5.846

5.840

1142

E. EARLIER AND J. EL KHATTABI

The average velocity was estimated by the relationship between distance X and modal time tm. From (14), concentrations were computed to be compared with the experimental concentrations. Dispersivity a was calculated, according to (11), by

Q m

)2 z

tm--Cmax

a =

(26)

4w

C o m p u t e d and experimental concentrations are given in Figure 1 and show a good agreement between experimental and theoretical values. The dispersivity values are very close to one half of the step of tile distance (Table 1).

4 -

Concentration (mg/l)

2-

1-

,J I

/

!

I

'

I

12

8

16

4-

~

I

Time (min) 2 0

Figure 1. E x p e r i m e n t a l and c o m p u t e d concentrations versus time for tracer test in sand. Q = 355030 m m a / m i n ; X = 1791 mm.

7.

DYE

TRACING

IN

CHALK

The tracer tests were carried out for three different discharges and two masses of fluoroscein (Table 2). The k p a r a m e t e r and the a power number are estimated by an iterative method. An arbitrary value of k is chosen and the power number a is calculated by (23). These two values are introduced into (22). Equation (21) allows us to check the validity of the k and a values. C o m p u t e d and experimental concentrations are given in Figure 2 and show a rather good agreement between experimental and computed values. The k values range between 0.951 and 0.980. The a values range between 1.429 and 1.575. Table 2. E x p e r i m e n t a l and c o m p u t e d d a t a for 1791 m m (chalk).

Q (mm3/min)

M (rag)

Concentration tm (min)

k

a

U (mm/min) 18.29

Max (mg/1)

63160

3.434

67

0.526

0.98

1.575

89160

3.891

45

0.572

0.951

1.429

26.26

148800

3.892

25.8

0.62

0.965

1.485

46.44

Proposal for a Probabilistic Iviodel 0-% 1 Concentration

1143

rag/l)

0.5-

0.3-

0.1-

I

,

,

,

.

,

.

1610

Figure 2. Experimental and computed concentrations versus time for tracer test in chalk. Q = 148800 m m 3 / m i n ; X = 1791 mm. Table 3. Experimental and computed data for 825 mm (chalk).

Q (mm3/min)

Concentration

M (mg)

tm (min)

63160

3.434

27.6

89160

3.891

18

148800

3.892

11.2

k

a

1.18

0.957

1.453

19.84

1.48

0.964

1.483

30.06

1.438

0.968

1.494

49.43

Max (Ing/1)

U (mm/min)

8. C O N C L U S I O N The probabilistic model presented in this paper reproduces laboratory flow column experiments for granular materials with a single intergranular or a dual porosity. Unlike complex numerical models, the probabilistic model seems not to be very sensitive to the exchange coefficients used to calibrate the dispersion phenomena. However, our results do not prejudge of the validity of the model on the field scale. Nevertheless, our experiments suggest that this model could be an alternative to the traditional numerical models.

REFERENCES 1. H.A. Einstein, Bed load transport as a probability problem, Th~se, Fed. Inst. of Technol., Zurich, (1937). 2. C.T. Yang and W.W. Sayre, Stochastic model of sediments dispersion, J. Hydraul.~ Div. Amer. Soc. Civil. Eng. 97 (HYZ), 265 288, (1971). 3. P. Todorovie, A stochastic process of monotonous sample function., Mat. Vesn. 4 (19), 149-158, (1967). 4. P. Todorovic, On some problems involving a random number of random variables, A n n . Math. Statis. 41 (3), 1059-1063, (1970). 5. I. Neretnieks, Diffusion in the rocks matrix. An important factor in radionuelide retardation?, J. Geophys. Res. 85 (88), 4379 4397, (1980). 6. I. Neretnieks, Age dating of groundwater in fissured rocks: Influence of water volume in micropores, W a t e r Res. Res. 17 (2), 421-422, (1981). 7. I. Neretnieks, A stochastic multi-channel model for solute t r a n s p o r t - - A n a l y s i s of tracer tests in fractured rock, J. Contarn. I-~gdrol. 55, 175-211, (2002). 8. D.R. Cameron and A. Klute, Convective-dispersive solute transport with a combined equilibrium and kinectic adsorption model, W a t e r Res. Res. 13 (1), 183-189, (1977). 9. E.M. Sudicky and E.O. Frind, Carbon 14 dating of groundwater in confined aquifer: Implications of aquitard diffusion, W a t e r Res. Res. 17 (4), 106C~1064, (1981).

1144

E. CARLIER AND J. EL KHATTABI

10. E.M. Sudicky and E.O. Frind, Contaminant transport in fractured porous media: Analytical solutions for a system of parallel fractures, Water Res. Res. 18 (6), 1634-1642, (1982). l l . P. Maloszewski and A. Zuber, Determining the turnover of groundwater systems with the aid of environnemental tracers, 1. Models and their applicability, J. Hydrol. 57, 207 231, (1982). 12. P. Maloszewski and A. Zuber, On the theory of tracer experiments in fissured rocks with a porous matrix, J. Hydrol. 79, a33 358, (1985). 13. P. Maloszewski and A. Zuber, Mathematieals Inodels for interpreting tracer experiments in fissured aquifers~ In The Application of Isotopic Techniques in the Study of the Hydrology of Fractured and Fissured Rocks, International Atomic Energy Agency, Vienna, (1989). 14. P. Maloszewski and A. Zuber, Mathematical modeling of tracer behavior in short term experiments in fissured rocks, Water Res. Res. 26 (7), 1517-1528, (1990). 15. P. Maloszewski and A. Zuber, Influence of matrix diffusion and exchange reactions on radiocarbon ages in fissured carbonate aquifers, Water Res. Res. 27 (8), 1937-1945, (1991). i6. C. Zheng and G.D. Bennet, Applied Contaminant Transport Modelling. Theory and Practice, Van Nostrand Rheinold, (1995). 17. L. Pang, M. Goltz and M. Close, Application of the method of temporal moments to interpret solute transport with sorption and degradation, J. Contain. Hydrol. 60, 123-134, (2003). 18. E. Carlier and C. Boulemia, A method for the analysis of tracer tests in groundwater, The Quarterly Journal of Engineering Geology and Hydrogeology 35, 291-294, (2002). 19. W. Feller, An Introduction to Probability Theory and Its Applications, Wiley, New York, (1966). 20. O. Kallenberg, Foundations of Modern Probability, Springer-Verlag, New York, (1997). 21. P.G. Saffman, Dispersion due to molecular diffusion and macroscopic mixing in flow through a network of capillaries, J. Fluid Reeh. 2, 194-208, (1959).