**A Stochastic Model for Hydrodynamic Dispersion**

### **3.1 Introduction**

Computational Modelling of Multi-Scale Non-Fickian Dispersion

 ) AAPE ()

<sup>2</sup> *R* ( 

and

 ) AAPE ( )

values

64 in Porous Media - An Approach Based on Stochastic Calculus

[1,2] [0.2,0.4] 0.35 0.47 0.53 0.39 0.2932 13.91 0.3891 12.98 [1,2] [1,2] 2 0.11 0.89 0.40 0.0006 18 0.6833 9.92 [0.8,1.6] [0.02,0.04] 0.15 0.89 0.11 0.42 0.7735 8.46 0.0108 17.94

To investigate the reason for largest differences in *R2* values for *α* and *β*, we change the

while keeping diffusion level an approximate constant. Table 2.4 shows that the bigger the contribution of a term containing a particular parameter (*Pα* or *Pβ*), the smaller the error (*AAPE*) and better the prediction (*R2*) for that parameter. Therefore, we conclude that the accuracy of a parameter in a nonlinear SDE is dependent on its term that contributes *pro rata*

In the data preparation stage, we use different time steps to solve SDEs and found 50 data points are sufficient to represent the realisation of SDEs. In addition, we emphasise the effect of the number of Wiener processes used to create training data sets. Increasing the number of Wiener processes boosts the performance of networks considerably and eliminates the over fitting problem. When over fitting occurs, the resulting network is accurate on the training set but perform poorly on the test set. When the number of Wiener processes used to generate training data sets is increased, the learning procedure finds common features amongst the training sets that enable the network to correctly estimate the parameter(s) in

In the ANN training procedure, we use early stopping to obtain the optimum test results. We also employ different MLP architectures, transfer functions, learning rates and momentums. However we find that these factors do not increase the performance of ANNs

The diffusion level in a SDE has a significant impact on the network performance. In the linear SDE, when the ratio of diffusion term and drift term is below 0.40, the network can estimate the parameter accurately ( <sup>2</sup> *R* >0.93). When the ratio reaches 0.67, the network estimates the parameter accurately only when Wiener processes in test sets and in training sets are similar. If the diffusion term is larger than the drift term, the network cannot predict the parameter(s) and only tends to give an average value of the parameters used for training datasets. For nonlinear SDEs, the estimation ability of a network is generally poorer than that for the linear SDEs. Furthermore, the accuracy of a parameter in a nonlinear SDE is

We can conclude that the classical neural networks method (MLP with backpropagation algorithm) provides a simple but robust parameter estimation approach for the SDEs that are under certain noisy conditions, but this estimation capability is limited for the SDEs having a high diffusion level. When the diffusion level is high (>10%-20%), the statistical

term in SDEs by altering parameters

γ Pα Pβ Pγ <sup>2</sup> *R* (

Table 2.4. Network performance for different parameters in the nonlinear SDE

Range of *α*

magnitudes of

to the drift term.

test data sets.

significantly.

Range of

term and

dependent on its term that contributes *pro rata* to the drift term.

methods also fail to estimate parameters accurately.

*β*

We have seen in chapter 1 that, in the derivation of advection–dispersion equation, also known as continuum transport model ( Rashidi et al. ,1999), the velocity fluctuations around the mean velocity enter into the calculation of solute flux at a given point through averaging theorems. The mean advective flux and the mean dispersive flux are then related to the concentration gradients through Fickian–type assumptions. These assumptions are instrumental in defining dispersivity as a measure of solute dispersion. Dispersivity is proven to be scale dependant.

To address the issue of scale dependence of dispersive fundamentally, it has been argued that a more realistic mathematical framework for modelling is to use stochastic calculus (Holden et al., 1996; Kulasiri and Verwoerd, 1999, 2002). Stochastic calculus deals with the uncertainty in the natural and other phenomena using nondifferentiable functions for which ordinary differentials do not exist (Klebaner, 1998). Stochastic calculus is based on the premise that the differentials of nondifferential functions can have meaning only through certain types of integrals such as Ito integrals which are rigorously developed in the literature. In addition, mathematically well-defined processes such as the Weiner process aid in formulating mathematical models of complex systems. Mathematical theories aside, one needs to question the validity of using stochastic calculus in each instance. In modelling the solute transport in porous media, we consider that the fluid velocity is fundamentally a random variable with respect to space and time and continuous but irregular, i.e., nondifferentiable. In many natural porous formations, geometrical structures are irregular and therefore, as fluid particles encounter porous structures, velocity changes are more likely to be irregular than regular. In many situations, we hardly have accurate information about the porous structure, which contributes to greater uncertainties. Hence, stochastic calculus provides a more sophisticated mathematical framework to model the advectiondispersion in porous media found in practical situations, especially involving natural porous formations. By using stochastic partial differential equations, for example, we could incorporate the uncertainty of the dispersion coefficient and hydraulic conductivity that are present in porous structures such as underground aquifers. The incorporation of the dispersivity as a random, irregular coefficient makes the solution of resulting partial differential equations an interesting area of study. However, the scale dependency of the dispersivity can not be addressed in this manner because the dispersivity itself is not a material property but a constant that depends on the scale of the experiment.

In this chapter we develop one dimensional model without resorting to Fickian assumptions and discuss the methods of estimating the parameters. As of many contracted description of

Let us assume that the terms up to second order would be sufficient to adequately represent the irregular behaviour of the flow, thus, the higher order derivatives of the flux greater

A Stochastic Model for Hydrodynamic Dispersion 67

2 2 <sup>1</sup> ( , ), <sup>2</sup>

*CJ J dx R x t*

*c*

(3.2.2)

(3.2.3)

(3.2.4)

(3.2.5)

(3.2.7)

*Jxt Vxt Cxt* ( , ) ( , ) ( , ). (3.2.6)

*x x*

2

*xx x*

2

*Jh J dC dt R x t dt x x* 

Equation (3.2.4) represents the mass conversation of the solute within the infinitesimal cylindrical volume shown in Figure 3.1. *Cxt* ( ,) describes the average solute concentration over the cylindrical volume and the smallest possible *x* would increase the accuracy of the stochastic model. However, if *x* is in the same order of magnitude of a typical grain size of the porous media, *Cxt* ( ,) would lose its meaning. Therefore, it is important to use a

Compared to the first term on the right hand side of equation (3.2.4), it is assumed that ( , ) 0. *R x t dt <sup>c</sup>* (However, this assumption needs to be tested for any given situation,

> *xx x Jh J dC dt x x*

The contaminant flux can be expressed in terms of the velocity in the *x* direction and the

*Vxt Vxt xt* ( , ) ( , ) ( , ),

( ) ( ,) ( ) *<sup>e</sup> K x <sup>p</sup> Vxt*

where ( , ) *Vxt* is the mean velocity of the flow that may be described by Darcy's law,

*K x*( ) = an average value of the hydraulic conductivity of the region,

*nx x* ,

*tx x* 

*xx x*

realistic *x* value for the scale under consideration for the computational solutions.

especially when porous media is extremely heterogeneous.) Under this assumption,

*C Jh J R xt*

<sup>2</sup> ( ,). <sup>2</sup>

<sup>2</sup> ( ,) . <sup>2</sup>

2 <sup>2</sup> . <sup>2</sup>

*c*

*c*

than three of the flux are negligible. Therefore, equation (3.2.1) can be written as,

*txx* 

3 4

*<sup>J</sup> <sup>J</sup> R x t dx dx R x x*

Substituting *<sup>x</sup> dx h* , we can write equation (3.2.2),

concentration of the contaminant is expressed as,

( ) *n x <sup>e</sup>* = porosity of the material,

The velocity is modelled as,

,)

3 4 <sup>1</sup> <sup>1</sup> ( ,) ( ) ( ) ( ) <sup>6</sup> <sup>24</sup> *x x*

2 3

.

Multiplying the both sides of equation (3.2.3) by *dt*, we obtain,

where

*c*

a natural phenomena the model presented in this chapter has its weaknesses. But we model the fluctuation of the solute velocity due to porous structure and incorporate the fluctuation in the mass conservation of solute. Then we need to characterise the fluctuations so that we can relate them to the porous structure

### **3.2 The Model Development**

The basic assumption on the formulation of this model is that the velocity of solute particles is fundamentally a stochastic variable with irregular but continuous realisations.

Let us consider an infinite cylindrical volume having a cross sectional area, *A*, within a onedimensional porous medium with the solute concentration of *C*(*x, t*) (Figure 3.1). Within the context of solute transport in heterogeneous groundwater systems, the velocity, *V*(*x, t*), the contaminant flux, *J*(*x, t*), and *C*(*x,t*) are stochastic functions of space and time, having irregular (may be highly irregular) and continuous realisations. Therefore, when formulating the mass conservation model for the solute transport, it is important to use higher order terms in the Taylor expansion. solute negligible. Therefore, assumption,

The mass balance for change of solute of the infinitesimal cylindrical volume for a small time interval, *t* , could be written as,

$$
\Delta \mathbf{C}(\mathbf{x}, t) n\_{\varepsilon} A \,\Delta \mathbf{x} = \left( f\_{\mathbf{x}}(\mathbf{x}, t) - f\_{\mathbf{x}}(\mathbf{x} + \Delta \mathbf{x}, t) \right) n\_{\varepsilon} A \,\Delta t \,\,\mathrm{d}\mathbf{x} \,, \,\mathrm{d}\mathbf{x} \,, \,\mathrm{d}\mathbf{x} = \left( \frac{f\_{\mathbf{x}}(\mathbf{x}, t) - f\_{\mathbf{x}}(\mathbf{x} + \Delta \mathbf{x}, t)}{\Delta \mathbf{x}} \right) \,\,\mathrm{d}\mathbf{x} \,\tag{3.2.1}
$$

$$
\left( \frac{\Delta \mathbf{C}}{\Delta t} \right)\_{\mathbf{x}, t} = \frac{\left( f\_{\mathbf{x}}(\mathbf{x}, t) - f\_{\mathbf{x}}(\mathbf{x} + \Delta \mathbf{x}, t) \right)}{\Delta \mathbf{x}} \,\,\mathrm{d}\mathbf{x} \,\tag{3.2.1}
$$

where *ne* is the effective porosity of the material, *Cxt* ( ,) and *t* are infinitesimal changes of solute concentration, *C*(*x*,*t*), and time, respectively.

For the convenience, let us indicate ( ,) *<sup>x</sup> J xt* as *<sup>x</sup> J* and ( ,) *<sup>x</sup> J x xt* as *x x J* .

The Taylor series expansion can be shown as

$$J\_{\chi\_{\mathbf{x}\leftarrow\mathbf{A}\mathbf{x}}} - J\_{\chi} = \frac{1}{1!} \frac{\partial f\_{\chi}}{\partial \mathbf{x}} \Delta \mathbf{x} + \frac{1}{2!} \frac{\partial^2 f\_{\chi}}{\partial \mathbf{x}^2} \left(\Delta \mathbf{x}\right)^2 + \frac{1}{3!} \frac{\partial^3 f\_{\chi}}{\partial \mathbf{x}^3} \left(\Delta \mathbf{x}\right)^3 + R(\boldsymbol{\varepsilon}) \text{ , } \mathbf{x}$$

where *R*( ) is the remainder of the series.

66 in Porous Media - An Approach Based on Stochastic Calculus

a natural phenomena the model presented in this chapter has its weaknesses. But we model the fluctuation of the solute velocity due to porous structure and incorporate the fluctuation in the mass conservation of solute. Then we need to characterise the fluctuations so that we

The basic assumption on the formulation of this model is that the velocity of solute particles

Let us consider an infinite cylindrical volume having a cross sectional area, *A*, within a onedimensional porous medium with the solute concentration of *C*(*x, t*) (Figure 3.1). Within the context of solute transport in heterogeneous groundwater systems, the velocity, *V*(*x, t*), the contaminant flux, *J*(*x, t*), and *C*(*x,t*) are stochastic functions of space and time, having irregular (may be highly irregular) and continuous realisations. Therefore, when formulating the mass conservation model for the solute transport, it is important to use

Figure 3.1. An infinitesimal cylindrical volume within a one-dimensional porous medium. The mass balance for change of solute of the infinitesimal cylindrical volume for a small

( ,) ( ( ,) ( ,) *Cxt nA x J xt J x xt nA t <sup>e</sup> <sup>x</sup> <sup>x</sup> <sup>e</sup>* ,

*C J xt J x xt*

where *ne* is the effective porosity of the material, *Cxt* ( ,) and *t* are infinitesimal

*J J <sup>J</sup> JJ x x x R x x x x*

 

*t x* 

For the convenience, let us indicate ( ,) *<sup>x</sup> J xt* as *<sup>x</sup> J* and ( ,) *<sup>x</sup> J x xt* as *x x J* .

( ( , ) ( , )) , *<sup>x</sup> <sup>x</sup>*

2 3

2 3 1 1 <sup>1</sup> ( ) ( ) () 1! 2! 3! *x x x*

,

2 3

(3.2.1)

,

*x t*

changes of solute concentration, *C*(*x*,*t*), and time, respectively.

is fundamentally a stochastic variable with irregular but continuous realisations.

can relate them to the porous structure

higher order terms in the Taylor expansion.

time interval, *t* , could be written as,

The Taylor series expansion can be shown as

*x x*

is the remainder of the series.

where *R*( ) 

**3.2 The Model Development** 

Let us assume that the terms up to second order would be sufficient to adequately represent the irregular behaviour of the flow, thus, the higher order derivatives of the flux greater than three of the flux are negligible. Therefore, equation (3.2.1) can be written as,

$$\frac{\partial \mathbf{C}}{\partial t} = -\frac{\partial f\_{\times}}{\partial \mathbf{x}} - \frac{1}{2} \frac{\partial^2 f\_{\times}}{\partial \mathbf{x}^2} d\mathbf{x} + R\_c(\mathbf{x}, t), \tag{3.2.2}$$

where 3 4 2 3 3 4 <sup>1</sup> <sup>1</sup> ( ,) ( ) ( ) ( ) <sup>6</sup> <sup>24</sup> *x x c <sup>J</sup> <sup>J</sup> R x t dx dx R x x* .

Substituting *<sup>x</sup> dx h* , we can write equation (3.2.2),

$$\frac{\partial \mathbb{C}}{\partial t} = -\frac{\partial \mathbb{J}\_x}{\partial \mathbf{x}} - \frac{h\_x}{2} \frac{\partial^2 \mathbb{J}\_x}{\partial \mathbf{x}^2} + R\_c(\mathbf{x}, t). \tag{3.2.3}$$

Multiplying the both sides of equation (3.2.3) by *dt*, we obtain,

$$d\mathbf{C} = -\left(\frac{\partial \mathbf{J}\_x}{\partial \mathbf{x}} + \frac{\mathbf{h}\_x}{2} \frac{\partial^2 \mathbf{J}\_x}{\partial \mathbf{x}^2}\right) dt + R\_c(\mathbf{x}, t) \, dt. \tag{3.2.4}$$

Equation (3.2.4) represents the mass conversation of the solute within the infinitesimal cylindrical volume shown in Figure 3.1. *Cxt* ( ,) describes the average solute concentration over the cylindrical volume and the smallest possible *x* would increase the accuracy of the stochastic model. However, if *x* is in the same order of magnitude of a typical grain size of the porous media, *Cxt* ( ,) would lose its meaning. Therefore, it is important to use a realistic *x* value for the scale under consideration for the computational solutions.

Compared to the first term on the right hand side of equation (3.2.4), it is assumed that ( , ) 0. *R x t dt <sup>c</sup>* (However, this assumption needs to be tested for any given situation, especially when porous media is extremely heterogeneous.) Under this assumption,

$$d\mathbf{C} = -\left(\frac{\partial \mathbf{J}\_{\mathbf{x}}}{\partial \mathbf{x}} + \frac{\mathbf{h}\_{\mathbf{x}}}{2} \frac{\partial^2 \mathbf{J}\_{\mathbf{x}}}{\partial \mathbf{x}^2}\right) dt. \tag{3.2.5}$$

The contaminant flux can be expressed in terms of the velocity in the *x* direction and the concentration of the contaminant is expressed as, ,)

$$\mathbf{J}(\mathbf{x},t) = \mathbf{V}(\mathbf{x},t) \, \mathbf{C}(\mathbf{x},t). \tag{3.2.6}$$

The velocity is modelled as,

$$V(\mathbf{x},t) = \overline{V}(\mathbf{x},t) + \xi(\mathbf{x},t),\tag{3.2.7}$$

where ( , ) *Vxt* is the mean velocity of the flow that may be described by Darcy's law,

$$
\overline{V}(\infty, t) = -\frac{K(\infty)}{n\_e(\infty)} \frac{\partial p}{\partial \mathfrak{x}}\,'\,'
$$

*K x*( ) = an average value of the hydraulic conductivity of the region,

( ) *n x <sup>e</sup>* = porosity of the material,

	- ( ,) *x t* = stochastic perturbation of the fluid velocity described by noise correlated in space and - correlated in time such that,

$$E\left[\mathcal{\xi}(\mathbf{x},t)\right] = 0,\tag{3.2.8}$$

Therefore, we can write equation (3.2.11) as,

The explanation on the transformation of

( ,) *x t dt* to *d t*

transformed by *d t*

The transformation of

Hilbert space (Hernandez, 1995). Equation (3.2.13) can be written as,

1 2

, where <sup>2</sup> 

*x x b*

2

 *e* 

vector field ( *d t*

(1970).

And

*d t* 

function,

*dC S V x t C x t dt S C x t x t dt* ( ,) ( ,) ( ,) ( ,) .

hand side need to be integrated as Ito integrals to obtain the concentration.

Equation (3.2.13) is a stochastic (partial) differential equation and both terms on the right

A Stochastic Model for Hydrodynamic Dispersion 69

fluctuating component of the travel length of a solute particle during of time interval *dt* . As we see later this fluctuating travel length component can be expressed as a random

ordinates depend on the covariance kernel 1 2 *q*(,) *x x* , and these are called basis functions in

*dC S V x t C x t dt S C x t d t* ( ,) ( ,) ( ,) (),

( ,) ( ,) ( ,) ( ,) ().

We call equation (3.2.15) a stochastic solute transport model (SSTM) in which *C*(*x,t*) is stochastic solute concentration at a given point in space and time; *Vxt* ( ,) is the expected value of stochastic velocity; *hx* approximates the spatial discretization interval, *dx*; and

In this model, instead of the Fickian assumption (equation (1.2.3)), a covariance kernel for the velocity fluctuations 1 2 *qx x* (,) is assumed. The emphasis on modelling of the velocity fluctuation as a random field of second order has its own strengths and weaknesses. The expectation is that by directly incorporating the velocity fluctuation in the solute mass conservation, we expect to reduce the scale dependency because the assumed covariance kernel is a function of the pore structure. The major weakness is that we usually do not have velocity data to construct the kernel for velocity fluctuations. On the other hand, different velocity kernels can be assumed based on plausible reasoning. In this chapter, for the illustration purposes, we assume that the covariance kernel to be exponentially decaying

In SSTM, the first term on the right hand side of equation (3.2.15) can be considered as the drift term of the stochastic partial differential equations (SPDE). The velocity, *Vxt* ( ,) , can be considered as the mean velocity of a local region but for the simplicity, we can assume that *Vxt* ( ,) is a constant in this chapter. This assumption is reasonable in practice as we

0 0

( ) can be considered as a "noise" term which models the velocity fluctuations.

*t t C x t S V x t C x t dt S C x t d t*

is the variance and *b* is the correlation length.

( ) which is a Wiener process increments in Hilbert space for a given *x*.

(3.2.15)

( ) can be understood if we interpret

( ,) *x t dt* to *d t*

( ) ) in an orthonormal set of co-ordinate space (Hilber space). The co-

(3.2.13)

(3.2.14)

( ) can be found in Jazwinski

( ,) *x t dt* can be

( ,) *x t dt* to be

$$E\left[\xi(\mathbf{x}\_1, t\_1)\xi(\mathbf{x}\_2, t\_2)\right] = q(\mathbf{x}\_1, \mathbf{x}\_2)\,\delta(t\_1 - t\_2) \,. \tag{3.2.9}$$

where 1 2 *qx x* (,) is the velocity covariance kernel in space, and 1 2 ( ) *t t* is the Dirac's delta function.

Since, 2 1 *tt t* , we can rewrite equation (3.2.9) as,

$$E\left[\xi(\mathbf{x}\_{1},t\_{1})\,\xi(\mathbf{x}\_{2},t\_{1}+\Delta t)\right] = q(\mathbf{x}\_{1},\mathbf{x}\_{2})\,\delta(\Delta t)\,. \tag{3.2.10}$$

In equation (3.2.6) the velocity, *Vxt* ( ,) , is modelled as a stochastic process with the fluctuating component ( ( ,) *x t* ) superimposed on the expectation of velocity (mean velocity). The expected value of the velocity *Vxt* ( ,) can be obtained from other considerations without involving Darcy's law. The use of Darcy's law is only one way of prescribing *Vxt* ( ,) within a given region of the porous medium. The fluctuating velocity ( ( ,) *x t* ) is assumed to be a zero-mean stochastic process as given by equation (3.2.7), and the two-time correlation function is expressed as a product of a function of space ( 1 2 *qx x* (,) ) and a function of time ( 1 2 ( ) *t t* ) in equation (3.2.9) and (3.2.10). (As ( ,) *x t* is a zeromean process, the two-time correlation function and the covariance are the same). The spatial function, 1 2 *qx x* (,) , reflects the contribution of the porous structure to the fluctuations and helps characterise the structural influence on the flow. The Dirac's delta function makes the stochastic process behave as a white noise process along the time line. For low velocity situations, such as in aquifers, this is not an unreasonable assumption, but the alternative approach of assembling a time correlation function along the time line would complicate the mathematics a fair deal, and unless the flow is turbulent, there is no justification to take that approach. law. structure

Substituting equations (3.2.5) and (3.2.6) into equation (3.2.4), we can obtain,

$$\begin{split} d\mathcal{C} &= -\frac{\partial}{\partial \mathbf{x}} \Big[ \overline{\mathcal{V}}(\mathbf{x},t) \mathbf{C}(\mathbf{x},t) + \mathbf{C}(\mathbf{x},t) \boldsymbol{\xi}(\mathbf{x},t) \Big] dt - \frac{h\_{\varepsilon}}{2} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} \Big[ \overline{\mathcal{V}}(\mathbf{x},t) \mathbf{C}(\mathbf{x},t) + \mathbf{C}(\mathbf{x},t) \boldsymbol{\xi}(\mathbf{x},t) \Big] dt, \\ &= -\left( \frac{\partial}{\partial \mathbf{x}} + \frac{h\_{\varepsilon}}{2} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} \right) \Big[ \overline{\mathcal{V}}(\mathbf{x},t) \mathbf{C}(\mathbf{x},t) + \mathbf{C}(\mathbf{x},t) \boldsymbol{\xi}(\mathbf{x},t) \Big] dt. \end{split} \tag{3.2.11}$$

An operator in space can be defined such that,

$$S = -\left(\frac{\partial}{\partial \mathbf{x}} + \frac{h\_x}{2} \frac{\partial^2}{\partial \mathbf{x}^2}\right) \tag{3.2.12}$$

for a given . *<sup>x</sup> h*

Therefore, we can write equation (3.2.11) as,

$$d\mathbb{C} = \mathcal{S}\{\bar{V}(\mathbf{x},t)\mathcal{C}(\mathbf{x},t)\}dt + \mathcal{S}\{\mathcal{C}(\mathbf{x},t)\mathcal{G}(\mathbf{x},t)\}dt.\tag{3.2.13}$$

Equation (3.2.13) is a stochastic (partial) differential equation and both terms on the right hand side need to be integrated as Ito integrals to obtain the concentration. ( ,) *x t dt* can be transformed by *d t* ( ) which is a Wiener process increments in Hilbert space for a given *x*. The explanation on the transformation of ( ,) *x t dt* to *d t* ( ) can be found in Jazwinski (1970).

The transformation of ( ,) *x t dt* to *d t* ( ) can be understood if we interpret ( ,) *x t dt* to be fluctuating component of the travel length of a solute particle during of time interval *dt* . As we see later this fluctuating travel length component can be expressed as a random vector field ( *d t* ( ) ) in an orthonormal set of co-ordinate space (Hilber space). The coordinates depend on the covariance kernel 1 2 *q*(,) *x x* , and these are called basis functions in Hilbert space (Hernandez, 1995).

Equation (3.2.13) can be written as,

$$d\mathbf{C} = \mathcal{S}\{V(\mathbf{x}, t)\mathcal{C}(\mathbf{x}, t)\}dt + \mathcal{S}\{\mathcal{C}(\mathbf{x}, t)d\mathcal{J}(t)\},\tag{3.2.14}$$

And

Computational Modelling of Multi-Scale Non-Fickian Dispersion

( , ) 0, (3.2.8)

( ,) *x t* ) superimposed on the expectation of velocity (mean

(3.2.9)

( ) *t t* is the Dirac's delta

. (3.2.10)

The

(3.2.12)

( ,) *x t* is a zero-

(3.2.11)

68 in Porous Media - An Approach Based on Stochastic Calculus


*E xt* 

11 21 1 2 *E xt xt t*

In equation (3.2.6) the velocity, *Vxt* ( ,) , is modelled as a stochastic process with the

velocity). The expected value of the velocity *Vxt* ( ,) can be obtained from other considerations without involving Darcy's law. The use of Darcy's law is only one way of prescribing *Vxt* ( ,) within a given region of the porous medium. The fluctuating velocity

( ,) *x t* ) is assumed to be a zero-mean stochastic process as given by equation (3.2.7), and the two-time correlation function is expressed as a product of a function of space ( 1 2 *qx x* (,) )

mean process, the two-time correlation function and the covariance are the same). The spatial function, 1 2 *qx x* (,) , reflects the contribution of the porous structure to the fluctuations and helps characterise the structural influence on the flow. The Dirac's delta function makes the stochastic process behave as a white noise process along the time line. For low velocity situations, such as in aquifers, this is not an unreasonable assumption, but the alternative approach of assembling a time correlation function along the time line would complicate the mathematics a fair deal, and unless the flow is turbulent, there is no

( ) *t t* ) in equation (3.2.9) and (3.2.10). (As

2 2

2 <sup>2</sup> 2

*x x* 

( ,) ( ,) ( ,) ( ,) ( ,) ( ,) ( ,) ( ,) , <sup>2</sup>

*x*

*<sup>h</sup> dC V x t C x t C x t x t dt V x t C x t C x t x t dt*

*<sup>x</sup> <sup>h</sup> <sup>S</sup>*

( , ) ( , ) ( , ) ( ), 11 22 1 2 1 2 *q xx t t*

( , ) ( , ) ( , )( ) *q xx t*

( ,) *x t* = stochastic perturbation of the fluid velocity described by noise correlated in

*p* = pressure head, and

space and

Since, 2 1 *tt t* , we can rewrite equation (3.2.9) as,

*E xt xt*

where 1 2 *qx x* (,) is the velocity covariance kernel in space, and 1 2

 

> 

Substituting equations (3.2.5) and (3.2.6) into equation (3.2.4), we can obtain,

<sup>2</sup> .

( ,) ( ,) ( ,) ( ,) <sup>2</sup>

*<sup>h</sup> V x t C x t C x t x t dt*

*x x*

function.

( 

fluctuating component (

and a function of time ( 1 2

justification to take that approach.

2

An operator in space can be defined such that,

*x*

*x x*

for a given . *<sup>x</sup> h*

$$\mathbf{C}(\mathbf{x},t) = \bigvee\_{\mathbf{0}}^{\mathsf{t}} \mathbf{S}\big(\overline{V}(\mathbf{x},t)\mathbf{C}(\mathbf{x},t)\big)dt + \bigvee\_{\mathbf{0}}^{\mathsf{t}} \mathbf{S}\big(\mathbf{C}(\mathbf{x},t)d\beta(t)\big). \tag{3.2.15}$$

We call equation (3.2.15) a stochastic solute transport model (SSTM) in which *C*(*x,t*) is stochastic solute concentration at a given point in space and time; *Vxt* ( ,) is the expected value of stochastic velocity; *hx* approximates the spatial discretization interval, *dx*; and *d t* ( ) can be considered as a "noise" term which models the velocity fluctuations.

In this model, instead of the Fickian assumption (equation (1.2.3)), a covariance kernel for the velocity fluctuations 1 2 *qx x* (,) is assumed. The emphasis on modelling of the velocity fluctuation as a random field of second order has its own strengths and weaknesses. The expectation is that by directly incorporating the velocity fluctuation in the solute mass conservation, we expect to reduce the scale dependency because the assumed covariance kernel is a function of the pore structure. The major weakness is that we usually do not have velocity data to construct the kernel for velocity fluctuations. On the other hand, different velocity kernels can be assumed based on plausible reasoning. In this chapter, for the illustration purposes, we assume that the covariance kernel to be exponentially decaying

function, 1 2 2 *x x b e* , where <sup>2</sup> is the variance and *b* is the correlation length.

In SSTM, the first term on the right hand side of equation (3.2.15) can be considered as the drift term of the stochastic partial differential equations (SPDE). The velocity, *Vxt* ( ,) , can be considered as the mean velocity of a local region but for the simplicity, we can assume that *Vxt* ( ,) is a constant in this chapter. This assumption is reasonable in practice as we

Now we can make use of Mercer's theorem (Hernandez, 1995) to express *q*(*x,y*) in terms of

A Stochastic Model for Hydrodynamic Dispersion 71

and convergence is uniform and absolute over *K × K.* Now let us define a set of (complex)

0, 1,2,... *EZ n <sup>n</sup>* , and

, 1,2,..... *EZ Zn m n nm* 

> 1 ( ): ( ) *n n <sup>n</sup>*

which converges in quadratic mean for any *x* in *K.* In addition, the following conditions

What Karhunen-Loeve (KL) expansion does is to model a random field by using two separate mathematical objects: continuous functions, which stem as the solutions of an integral equation (equation (3.3.1)) and random variables. If we assume the random variables ( *Zi* ) to be standard Gaussian (*N*(0,1)), then Karhunen-Loeve expansion takes the

and,

for any *x0* in *K.* 

( ) ( ), 1,2,..... *E xZ x n n nn*

1 ( ): ( ) *nn n n*

 *x x Z* 

*KL* expansion provides a way of modelling a random field of a single variable, which is considered to be the space variable, in terms of a set of orthonormal basis functions of a Hilbert space and a discrete set of standard Gaussian variables. We extend this development

symmetric and positive definite covariance function 1 2 *qx x* (,) , and *δ*-correlated in time (or white in time). As mentioned before this means that velocity changes in a porous medium are influenced by the porous matrix but behaves like white noise with respect to time, i.e., correlation in time is a Dirac delta function. Suppose that the mean velocity *Vxt* ( ,) of a region is 1 m per day, then within a second a solute particle travels 0.01 mm on the average, and within 0.001 days it travels 1 mm, which increases the probability of solute particle changing its course. We can assume that in slow moving fluids, the velocity fluctuation at

 *x xZ* 

 

<sup>2</sup> <sup>0</sup> lim ( ) ( ) 0

 

*Ex x* 

 

*n* .

, (3.3.5)

, (3.3.6)

. (3.3.7)

( ,) *x t* , which is spatially and

( , ) ( ) [0, ] *xt H L R T* where <sup>0</sup> *H* is a

( ,) *x t* to be spatially correlated through a

2

( ) *x* can be expressed as a

1 (,) () () *kk k <sup>k</sup> qxy x y* 

eigenfunctions (*ψi* ) and the corresponding eigenvalues (*λi*) ,

random variables {*Zi : i=*1, 2, 3*..*.} with the following moments:

Then the Karhunen-Loeve theorem states that random field

0

to model the velocity fluctuation of equation (3.2.6),

temporally correlated in such a way that <sup>0</sup>

separable Hilbert space in which we assume

*x x*

series expansion of Zi s,

hold,

following form,

often calculate the Darcian velocity over a region. As *C*(*x*,*t*) is a stochastic variable, at a given point in space, the integral can be evaluated as an Ito integral to simplify the numerical solution. This argument can also be applied to the dispersion term, the second integral of the right hand side of equation (3.2.15), but *d t* ( ) term has to be evaluated. (We use the term "dispersion" instead of "diffusion" to denote the second integral of the right hand side of equation (3.2.15) to allude to the physical nature that is being represented by the term. But in the mathematical literature, the term "diffusion" is used to denote the stochastic (Ito or Stratonovich) integral of analogous stochastic differential equations because it is a diffusion process in mathematical sense.) It should also be noted that *Vxt* ( ,) could also vary from location to location at a higher scale, and we can assume that we can identify regions of sufficiently large magnitudes within which *Vxt* ( ,) can be treated as a smoothly varying function if not a constant. However, this assumption is not necessary if we know the mean velocity profiles based on more aggregate properties. ( ,) *x t* , on the other hand, relaxes the assumption expressed by equation (1.2.3) and allows us to include a much more realistic assumption as to how the noise behaves. There is some experimental evidence to suggest that the normalised longitudinal velocity covariances for different flow rates had similar exponential behaviour with respect to time in a homogenous porous media (Moroni and Cushman, 2001). However, there is no experimental evidence to date as to how the covariance of velocity fluctuation around the mean velocity, ( ,) *x t* , behaves with respect to space, but the exponential decay in this situation is quite plausible. hand identify )

### **3.3 Construction of** *d t* ( ) **Random Fields Using Spectral Expansion**

We summarise the pertinent theoretical background here to explain the assumptions in our model of velocity fluctuations. Let *K* be a compact set in *Rd,* over which we define a second order random field, having a covariance function, *q*(*x,y*) which is assumed to be square integrable over *K × K* (see Hernandez,1995):

$$\iint\_{K \times K} \left| q(\mathbf{x}\_\prime \, y) \right|^2 d\mathbf{x} \, dy < +\infty \,\,\omega$$

Let *A: L2*(*K*)*→L2*(*K*) be the integral operator

$$A\,\phi(\mathbf{x}) = \left[ \_{\kappa}q(\mathbf{x},y)\,\phi(y)\,dy\right] \tag{3.3.1}$$

with the set of eigen values {*λi* : *i=* 1, 2*,* 3*.....*} and the set of orthonormal eigenfunctions {*ψi* : *i=* 1, 2, 3 *...*}. This means

$$A\,\psi\_{i} = \lambda\_{i}\psi\_{i}, \quad i = 1, 2, 3...,\tag{3.3.2}$$

$$\text{and} \quad \langle \boldsymbol{\psi}\_{\text{i}}, \boldsymbol{\psi}\_{\text{j}} \rangle = \delta\_{\boldsymbol{\psi}} \quad \text{i}, \mathbf{j} = \mathbf{1}, \mathbf{2}, \mathbf{3} \dots , \tag{3.3.3}$$

$$\text{where}\quad \left\langle \boldsymbol{\nu}\_{\boldsymbol{\cdot}}, \boldsymbol{\nu}\_{\boldsymbol{\cdot}} \right\rangle = \left\{ \boldsymbol{\nu}\_{\boldsymbol{\cdot}}(\mathbf{x}) \,\overline{\boldsymbol{\nu}\_{\boldsymbol{\cdot}}(\mathbf{x})} \, d\mathbf{x} \right\}. \tag{3.3.4}$$

If q(*x,y*) is continuous then , 0 *A i i* for 2 ( ) *<sup>i</sup> L K* .

( ) term has to be evaluated. (We use the term

( ,) *x t* , on the other hand, relaxes the

( ,) *x t* , behaves with respect to

*dy* (3.3.1)

, (3.3.2)

*i j* , (3.3.3)

*x x dx* . (3.3.4)

70 in Porous Media - An Approach Based on Stochastic Calculus

often calculate the Darcian velocity over a region. As *C*(*x*,*t*) is a stochastic variable, at a given point in space, the integral can be evaluated as an Ito integral to simplify the numerical solution. This argument can also be applied to the dispersion term, the second integral of the

"dispersion" instead of "diffusion" to denote the second integral of the right hand side of equation (3.2.15) to allude to the physical nature that is being represented by the term. But in the mathematical literature, the term "diffusion" is used to denote the stochastic (Ito or Stratonovich) integral of analogous stochastic differential equations because it is a diffusion process in mathematical sense.) It should also be noted that *Vxt* ( ,) could also vary from location to location at a higher scale, and we can assume that we can identify regions of sufficiently large magnitudes within which *Vxt* ( ,) can be treated as a smoothly varying function if not a constant. However, this assumption is not necessary if we know the mean

assumption expressed by equation (1.2.3) and allows us to include a much more realistic assumption as to how the noise behaves. There is some experimental evidence to suggest that the normalised longitudinal velocity covariances for different flow rates had similar exponential behaviour with respect to time in a homogenous porous media (Moroni and Cushman, 2001). However, there is no experimental evidence to date as to how the

( ) **Random Fields Using Spectral Expansion** We summarise the pertinent theoretical background here to explain the assumptions in our model of velocity fluctuations. Let *K* be a compact set in *Rd,* over which we define a second order random field, having a covariance function, *q*(*x,y*) which is assumed to be square

> 2 (,)

*q x y y*

with the set of eigen values {*λi* : *i=* 1, 2*,* 3*.....*} and the set of orthonormal eigenfunctions {*ψi* :

and , , 1,2,3..... *i j ij*

 

where , () () *ij i j K*

*i i* for 2 ( ) *<sup>i</sup> L K* .

.

*q x y dxdy*

right hand side of equation (3.2.15), but *d t*

velocity profiles based on more aggregate properties.

covariance of velocity fluctuation around the mean velocity,

integrable over *K × K* (see Hernandez,1995):

Let *A: L2*(*K*)*→L2*(*K*) be the integral operator

If q(*x,y*) is continuous then , 0 *A*

**3.3 Construction of** *d t*

*i=* 1, 2, 3 *...*}. This means

space, but the exponential decay in this situation is quite plausible.

*K K*

 () (,)() *A x <sup>K</sup>* 

 , 1,2,3... *A i i ii* 

> 

  Now we can make use of Mercer's theorem (Hernandez, 1995) to express *q*(*x,y*) in terms of eigenfunctions (*ψi* ) and the corresponding eigenvalues (*λi*) , *λi*

$$q(\mathbf{x}, \mathbf{y}) = \sum\_{k=1}^{\alpha} \mathbb{A}\_{k} \boldsymbol{\mu}\_{k}(\mathbf{x}) \overline{\boldsymbol{\mu}\_{k}(\mathbf{y})} \; , \tag{3.3.5}$$

and convergence is uniform and absolute over *K × K.* Now let us define a set of (complex) random variables {*Zi : i=*1, 2, 3*..*.} with the following moments:

$$\begin{array}{ll} \text{EZ}\_{n} = 0, & n = 1, 2, \dots & \text{ , and } \\\\ \text{EZ}\_{n} \overline{\text{Z}\_{m}} = \text{A}\_{n} \, \delta\_{nm}, & n = 1, 2, \dots & \text{ . } \end{array}$$

Then the Karhunen-Loeve theorem states that random field ( ) *x* can be expressed as a series expansion of Zi s,

$$\varphi(\mathbf{x}) \coloneqq \sum\_{n=1}^{\infty} \varphi\_n(\mathbf{x}) Z\_{n \text{ \textquotedblleft}} \tag{3.3.6}$$

which converges in quadratic mean for any *x* in *K.* In addition, the following conditions hold,

$$E\boldsymbol{\gamma}(\mathbf{x})\,\mathbf{Z}\_n = \mathbb{A}\_n \boldsymbol{\Psi}\_n(\mathbf{x}), \quad n = 1, 2, \dots \quad \text{and},$$

$$\lim\_{\mathbf{x} \to \mathbf{x}^0} E\left|\boldsymbol{\gamma}(\mathbf{x}) - \boldsymbol{\gamma}(\mathbf{x}^0)\right|^2 = 0 \qquad \text{for any} \quad \mathbf{x}^0 \text{ in } K.$$

What Karhunen-Loeve (KL) expansion does is to model a random field by using two separate mathematical objects: continuous functions, which stem as the solutions of an integral equation (equation (3.3.1)) and random variables. If we assume the random variables ( *Zi* ) to be standard Gaussian (*N*(0,1)), then Karhunen-Loeve expansion takes the following form,

$$\mathcal{Y}(\mathbf{x}) \coloneqq \sum\_{\boldsymbol{\pi} \in \mathbf{I}}^{\boldsymbol{\pi}} \sqrt{\lambda\_{\boldsymbol{\pi}}} \,\boldsymbol{\nu}\_{\boldsymbol{\pi}}(\mathbf{x}) \,\boldsymbol{Z}\_{\boldsymbol{\pi}} \,\,\, \,\tag{3.3.7}$$

*KL* expansion provides a way of modelling a random field of a single variable, which is considered to be the space variable, in terms of a set of orthonormal basis functions of a Hilbert space and a discrete set of standard Gaussian variables. We extend this development to model the velocity fluctuation of equation (3.2.6), ( ,) *x t* , which is spatially and temporally correlated in such a way that <sup>0</sup> 2 ( , ) ( ) [0, ] *xt H L R T* where <sup>0</sup> *H* is a separable Hilbert space in which we assume ( ,) *x t* to be spatially correlated through a symmetric and positive definite covariance function 1 2 *qx x* (,) , and *δ*-correlated in time (or white in time). As mentioned before this means that velocity changes in a porous medium are influenced by the porous matrix but behaves like white noise with respect to time, i.e., correlation in time is a Dirac delta function. Suppose that the mean velocity *Vxt* ( ,) of a region is 1 m per day, then within a second a solute particle travels 0.01 mm on the average, and within 0.001 days it travels 1 mm, which increases the probability of solute particle changing its course. We can assume that in slow moving fluids, the velocity fluctuation at respect

As the upper limit of integration is a constant, equation (3.3.10) is a homogenous Fredholm

A Stochastic Model for Hydrodynamic Dispersion 73

We now discuss the general solution of the integral equation generated by equation (3.3.10) considering the solutions of a general Fredholm integral equation of the form (Polyanin et

> ( ) 1 () ( ) *b xt <sup>a</sup> y xA e y t dt f x*

Then the function *y = y*(*x*) obeys the following second order linear non-homogeneous

<sup>2</sup> (2 1 ) ( ) ( ) *xx xx y*

The boundary conditions for this ordinary differential equation (ODE) have the form,

() () () () *<sup>x</sup> <sup>x</sup> y a ya f a f a*

() () () () *<sup>x</sup> <sup>x</sup> y b yb f b f b*

Polyanin et al. (1998) explain why the given ODE (equation (3.3.16)) under these boundary

(2 1 ) 0 *A* , the general solution of equation (3.3.13) is given by,

(2 1 ) 0 *A* , the general solution of equation (3.3.13) is given by

*k* 

2

, (3.3.18)

*<sup>a</sup> y x C Cx f x A xt f t dt* . (3.3.19)

2 1 () ( ) ( ) () [ ( )] ( ) *<sup>A</sup> <sup>x</sup> y x C Cosh k x C Sinh kx f x Sinh k x t f t dt <sup>k</sup>*

2 1 () ( ) ( ) () [ ( )] ( ) *<sup>A</sup> <sup>x</sup> y x C Cos k x C Sin kx f x Sin k x t f t dt*

*A* , the general solution of equation (3.3.13) is given by,

1 2 ( ) ( ) 4 1 ( ) () *<sup>x</sup>*

These results (Cases 1, 2, and 3) can be applied to the integral equation (3.3.11) by changing

<sup>2</sup> 12 1 (2 1 ) *<sup>A</sup>* <sup>0</sup> *b b* 

 

, (3.3.17)

conditions determines the solutions of the original integral equation (3.3.12).

1 2 0

( 2 1) *A* , and *C*1 and *C*2 are arbitrary constants.

<sup>1</sup> <sup>2</sup> <sup>0</sup>

The constants *C*<sup>1</sup> and *C*2 are determined by the boundary conditions.

, (3.3.13)

*A y f x f x* . (3.3.14)

,and (3.3.15)

. (3.3.16)

. (3.3.20)

integral equation of the second kind.

where *A1* and *θ* are constants.

differential equation with constant coefficients,

al., 1998),

**Case 1.** For

**Case 2.** For

where *k* 

where *k A* (2 1 ) .

**Case 3.** For 2 1 

notation, and observing that

 

 

one small interval in time is independent of that of within another interval in time. By incorporating appropriate discretisation time and space intervals in the solution process, the white-in-time assumption leads to a simpler but physically plausible model of velocity fluctuation. Therefore, we can express ( ,) *x t* as a second-order random field with the moments as given in equation (3.2.9).

Based on stochastic calculus considerations (Unny, 1985), ( ,) *x t dt* can be considered as a Wiener process increment *d t* ( ) in Hilbert space <sup>0</sup> *H* for a given *x;* and to obtain strong solutions to the stochastic partial differential equation (3.2.14) numerically, this representation is particularly useful. The orthonormal functions, {*ψi : i=*1, 2, 3.*..* }, provide a basis for <sup>0</sup> *H* , and then we can use the standard Wiener increments (*dbi*(*t*)) to replace random variables in *KL* expansion (equation (3.3.7)) to construct the random field *dβ*(*t*):

$$d\mathcal{J}(t) \coloneqq \sum\_{n=1}^{n} \sqrt{\mathcal{A}\_n} \, \psi\_n(\mathbf{x}) \, db\_n(t) = \mathcal{J}(\mathbf{x}, t)dt \,. \tag{3.3.8}$$

Intrinsic in this model is the fact that for any given element <sup>0</sup> *e H* , the inner product ( ), *t e* is a real Wiener process having the correlation,

$$E\left[<\beta(t), e\_1><\beta(t), e\_2>\right] = \iint\_{\mathbb{X}\times\mathbb{X}} q(\mathbf{x}\_1, \mathbf{x}\_2) e\_1(\mathbf{x}\_2) e\_2(\mathbf{x}\_1) d\mathbf{x}\_1 d\mathbf{x}\_2 \,. \tag{3.3.9}$$

An approximation for equation (3.3.8) can be constructed by considering the first *m* terms in the expansion which converges to a desired accuracy. As the increments in standard Wiener process have zero means and *dt* variances, care must be given to the choice of the time increments within the context of the problem. Smaller the time increments, higher the accuracy of KL expansion as in the case of solving stochastic differntial equations numerically (Kloeden and Platen, 1992).

Equation (3.3.8) gives a series expansion for the stochastic process *d t* ( ) but eigen functions and eigenvalues must be found for a given kernel. This task is not straight forward and often involves solving the Fredholm integral equations. Let us assume, based on previous discussion, an exponentially decaying kernel in the *x*-direction in this chapter,

$$q(\mathbf{x}\_1, \mathbf{x}\_2) = \sigma^2 e^{\frac{-\left|\mathbf{x}\_1 - \mathbf{x}\_2\right|}{\theta}}.\tag{3.3.10}$$

where *x1* and *x2* are two neighbouring points on *x*. The correlation length *b* signifies the extent to which the correlation of the stochastic process is decayed along the *x*-axis. *σ2* is the variance which acts as an amplitude factor. To obtain the orthornormal functions for the kernel in equation (3.3.10), the following integral equation must be solved within the domain [0, *a*], where *a* is the length (scale) of the experiment,

$$\int\_{0}^{\boldsymbol{\mu}} q(\mathbf{x}\_{1}, \mathbf{x}\_{2}) \, \, f\_{i}(\mathbf{x}\_{2}) \, \, d\mathbf{x}\_{2} = \mathcal{A}\_{i} f\_{i}(\mathbf{x}\_{1}) \,. \tag{3.3.11}$$

Equation (3.3.11) can be written as,

$$f\_i(\mathbf{x}\_1) = \frac{1}{\lambda\_i} f\_0^{\prime} \sigma^2 e^{\frac{-|\mathbf{x}\_1 - \mathbf{x}\_2|}{b}} \quad f\_i(\mathbf{x}\_2) \text{ } d\mathbf{x}\_2 \text{ }. \tag{3.3.12}$$

As the upper limit of integration is a constant, equation (3.3.10) is a homogenous Fredholm integral equation of the second kind.

We now discuss the general solution of the integral equation generated by equation (3.3.10) considering the solutions of a general Fredholm integral equation of the form (Polyanin et al., 1998), We

( ) 1 () ( ) *b xt <sup>a</sup> y xA e y t dt f x* , (3.3.13)

where *A1* and *θ* are constants.

Computational Modelling of Multi-Scale Non-Fickian Dispersion

( ,) *x t* as a second-order random field with the

( ,) *x t dt* can be considered as a

(3.3.10)

. (3.3.12)

*f x* . (3.3.11)

( ) but eigen

( ) in Hilbert space <sup>0</sup> *H* for a given *x;* and to obtain strong

. (3.3.8)

72 in Porous Media - An Approach Based on Stochastic Calculus

one small interval in time is independent of that of within another interval in time. By incorporating appropriate discretisation time and space intervals in the solution process, the white-in-time assumption leads to a simpler but physically plausible model of velocity

solutions to the stochastic partial differential equation (3.2.14) numerically, this representation is particularly useful. The orthonormal functions, {*ψi : i=*1, 2, 3.*..* }, provide a basis for <sup>0</sup> *H* , and then we can use the standard Wiener increments (*dbi*(*t*)) to replace random variables in *KL* expansion (equation (3.3.7)) to construct the random field *dβ*(*t*):

( ): ( ) ( ) ( , ) *nn n <sup>n</sup>*

Intrinsic in this model is the fact that for any given element <sup>0</sup> *e H* , the inner product

 <sup>1</sup> <sup>2</sup> 1 2122 1 1 2 ( ), ( ), ( , ) ( ) ( ) *K K E t e t e q x x e x e x dx dx*

An approximation for equation (3.3.8) can be constructed by considering the first *m* terms in the expansion which converges to a desired accuracy. As the increments in standard Wiener process have zero means and *dt* variances, care must be given to the choice of the time increments within the context of the problem. Smaller the time increments, higher the accuracy of KL expansion as in the case of solving stochastic differntial equations

functions and eigenvalues must be found for a given kernel. This task is not straight forward and often involves solving the Fredholm integral equations. Let us assume, based on previous discussion, an exponentially decaying kernel in the *x*-direction in this chapter,

> 2 1 2 (,) .

where *x1* and *x2* are two neighbouring points on *x*. The correlation length *b* signifies the extent to which the correlation of the stochastic process is decayed along the *x*-axis. *σ2* is the variance which acts as an amplitude factor. To obtain the orthornormal functions for the kernel in equation (3.3.10), the following integral equation must be solved within the

> <sup>0</sup> 12 2 2 <sup>1</sup> (,) () () *<sup>a</sup> <sup>i</sup> i i q x x f x dx*

2 1 0 2 2 <sup>1</sup> ( ) ( ) *x x a b i i i*

*f x e f x dx* 

1 2

*<sup>b</sup> qx x e* 

1 2

*x x*

 *x db t x t dt* 

. (3.3.9)

1

Equation (3.3.8) gives a series expansion for the stochastic process *d t*

 

fluctuation. Therefore, we can express

Based on stochastic calculus considerations (Unny, 1985),

*d t* 

( ), *t e* is a real Wiener process having the correlation,

 

domain [0, *a*], where *a* is the length (scale) of the experiment,

numerically (Kloeden and Platen, 1992).

Equation (3.3.11) can be written as,

moments as given in equation (3.2.9).

Wiener process increment *d t*

 

Then the function *y = y*(*x*) obeys the following second order linear non-homogeneous differential equation with constant coefficients,

$$
\mathcal{Y}\_{\text{xx}}'' + \theta(2A1 - \theta)\mathcal{Y} = f\_{\text{xx}}''(\mathbf{x}) - \theta^2 \; f(\mathbf{x})\,. \tag{3.3.14}
$$

The boundary conditions for this ordinary differential equation (ODE) have the form,

$$y\_x'(a) + \theta \, y(a) = f\_x'(a) + \theta \, f(a) \text{, and } \tag{3.3.15}$$

$$
\partial\_x \underline{v}'(b) - \theta \,\underline{v}(b) = f'\_x(b) - \theta \, f(b) \,. \tag{3.3.16}
$$

Polyanin et al. (1998) explain why the given ODE (equation (3.3.16)) under these boundary conditions determines the solutions of the original integral equation (3.3.12).

**Case 1.** For (2 1 ) 0 *A* , the general solution of equation (3.3.13) is given by,

$$y(\mathbf{x}) = \mathbb{C}\_1 \operatorname{Cash}(k\,\mathbf{x}) + \mathbb{C}\_2 \operatorname{Simh}(k\,\mathbf{x}) + f(\mathbf{x}) - \frac{2A1\,\theta}{k} \left[ \_0^x \operatorname{Simh}[k(\mathbf{x} - t)] f(t) \, dt \right. \tag{3.3.17}$$

where *k* ( 2 1) *A* , and *C*1 and *C*2 are arbitrary constants.

**Case 2.** For (2 1 ) 0 *A* , the general solution of equation (3.3.13) is given by

$$y(\mathbf{x}) = \mathbf{C}\_1 \cos(k\mathbf{x}) + \mathbf{C}\_2 \sin(k\mathbf{x}) + f(\mathbf{x}) - \frac{2A1\theta}{k} \int\_0^\mathbf{x} \text{Sin}[k(\mathbf{x} - t)] f(t) dt \,\tag{3.3.18}$$

where *k A* (2 1 ) .

**Case 3.** For 2 1 *A* , the general solution of equation (3.3.13) is given by,

$$f(\mathbf{x}) = \mathbf{C}\_1 + \mathbf{C}\_2 \mathbf{x} + f(\mathbf{x}) - 4A1^2 \int\_{\boldsymbol{\mu}}^{\boldsymbol{\mu}} (\mathbf{x} - t) f(t) dt \,. \tag{3.3.19}$$

The constants *C*<sup>1</sup> and *C*2 are determined by the boundary conditions.

These results (Cases 1, 2, and 3) can be applied to the integral equation (3.3.11) by changing notation, and observing that

$$\left(2A1-\theta\right) = \left(\frac{1}{b}\right)\left(\frac{2\sigma^2}{\mathcal{A}} + \frac{1}{b}\right) > 0\ .\tag{3.3.20}$$

The eigenvalues of equation (3.3.13) for constant <sup>2</sup> over [0, *a*] are given by

$$
\lambda\_{\parallel} = \frac{2\theta\sigma^2}{\alpha\_{\parallel}^2 + \theta^2},
\tag{3.3.21}
$$

Now for the sake of argument, let us assume that all the Wiener increments are the same, a situation which is never encountered in computation. Then we can take *dbj* s out of the

A Stochastic Model for Hydrodynamic Dispersion 75

By substituting for the differential operator *S* and for ( ) *<sup>j</sup> f x* , after symbolic manipulations,

2 2

*C C SCxt <sup>f</sup> <sup>x</sup> <sup>C</sup>*

*j j j j j*

1 1 1 1

*h h Sin x Cos x*

2 2 3

This exercise shows the complexity of the dispersion term and it is related to the local spatial

*C*(*x,t*). *<sup>j</sup>* , *j* and *j* coefficients are also dependent on *λi* s, *hx*, and ωi s, and they are analogous to Fourier series, i.e., the coefficients themselves have the form of spectral expansions. Recall that *hx* is a scale length introduced to retain the second order derivative

the discretization length in *x* direction. If we make *hx*=0 then *<sup>j</sup>* =0, which eliminates the second order spatial derivative of the concentration in equation (3.4.3). Even in a simplified setting of having a common standard Wiener increment (*dbi*(*t*)) for all *x*, and assuming a differentiable function for *C*(*x*,*t*), this analysis shows that the dispersion term depicts much more complicated behaviours that cannot be modelled by scale independent simplifying assumptions. Therefore, one should note that as time increases the effects of Wiener process increases, and it would be worthwhile to investigate the dependence of the random field

( ) for different scales of experiments. In other words, from equation (3.2.7), we see that,

*dX t V x t dt x t dt* () ( ,) ( ,)

*j jx j j x j j j j h h*

 

*x* 

of fluxes in the SSTM. We can assume that ( ) *<sup>x</sup> h x* , i.e., 0

where *X*(*t*) is travel distance of a solute particle in *x*-direction.

1 () ( ) *j jx <sup>j</sup> j j j x j*

 

( ) ( ) <sup>2</sup> <sup>2</sup>

*j j j j*

<sup>1</sup> () () <sup>2</sup> *j x j j x j j*

*m m m m*

*m*

*j*

*x x*

*Sin x h Cos x*

*Sin x Cos x*

, the second order spatial derivative

 

*x <sup>h</sup> Lim x* 

, (3.4.3)

( ( ,) ( ) )

.

*SCxt f x*

*j j*

; (3.4.3a)

; and (3.4.3b)

. (3.4.3c)

, as well as

where *x* is

2 2 *C x* 

0 *<sup>x</sup>*

, (3.4.4)

summation, and explore the behaviour of the term, 1

( ( ,) ( ) )

*N* 

*N* 

*N*  *h*

 we can express,

gradient of the concentration, *<sup>C</sup>*

where,

*d t* 

where <sup>1</sup> *b* , and

*ωi* s are the roots of the following equation,

$$\operatorname{Tan}(\phi\_{\!\!\! }a) = \frac{2\alpha\_{\!\!\! }\theta}{\phi\_{\!\!\! }^{2} - \theta^{2}} \,. \tag{3.3.22}$$

The orthonormal basis functions of the Hilbert space associated with the given exponential kernel (equation (3.3.10)) are the normalised eigenfunctions,

$$f\_i(\mathbf{x}) = \frac{1}{\sqrt{N}} \left( \text{Sin}(o\rho\mathbf{x}) + \frac{o\_i}{\theta} \text{Cos}(o\rho\mathbf{x}) \right)\_{\text{'}} \tag{3.3.23}$$

where,

$$N = \frac{1}{2}a \left(1 + \frac{a\rho\_{\parallel}^2}{\theta^2}\right) - \frac{1}{4a\rho\_{\parallel}} \left(1 - \frac{a\rho\_{\parallel}^2}{\theta^2}\right) \text{Sim}(2a\rho\_{\parallel}a) - \frac{1}{2\theta} (\text{Cov}(2a\rho a) + 1)\,. \tag{3.3.24}$$

As seen from equation (3.3.21), the eigenvalues are dependent on the correlation length *b* and the variance <sup>2</sup> directly and to the length of the experiment, *a*, indirectly through equation (3.3.22). Equation (3.3.22) is a transendental equation with roots *ωi*. Firstly, equation (3.3.22) has to be solved numerically for a given *a* and *b*, then the eigenvalues *λi* s and finally the basis functions of <sup>0</sup> *H* . Then the random field *d t* ( ) for a given number of terms, *m,* can be evaluated by using equation (3.3.8) in conjunction with the standard Wiener increments.

### **3.4 The Dispersion and Travel Length Fluctuations**

In this section, we explore the dispersion term in the SSTM to understand its behaviour in simplified settings; as we have seen previously, *SCxt d t* ( ( , ) ( )) can be expanded by using the spectral expansion,

$$\mathbf{S(C(x,t)d\beta(t)) = \mathbf{S(C(x,t)\sum\_{j=1}^{m} f\_j(x)\sqrt{\lambda\_j} \, db\_j(t))}\,. \tag{3.4.1}$$

Let us assume that we use the same realisation of the standard Wiener process across all *x* in evaluating this term, i.e., *dbj* is independent of x*.* Then we can express the dispersion term for this situation,

$$S(\mathbf{C}(\mathbf{x},t) \, d\beta(t)) = \sum\_{j=1}^{m} S(\mathbf{C}(\mathbf{x},t) \, f\_{\rangle} \sqrt{\mathcal{A}\_{j}} \, d\boldsymbol{b}\_{j}(t) \,. \tag{3.4.2}$$

Now for the sake of argument, let us assume that all the Wiener increments are the same, a situation which is never encountered in computation. Then we can take *dbj* s out of the summation, and explore the behaviour of the term, 1 ( ( ,) ( ) ) *m j j j SCxt f x* .

By substituting for the differential operator *S* and for ( ) *<sup>j</sup> f x* , after symbolic manipulations, we can express,

$$\sum\_{j=1}^{m} \mathbf{S}(\mathbf{C}(\mathbf{x}, t) \, f\_j(\mathbf{x}) \sqrt{\mathcal{A}\_j}) = \left[ \sum\_{j=1}^{m} \Phi\_j \right] \frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2} + \left[ \sum\_{j=1}^{m} \Psi\_j \right] \frac{\partial \mathbf{C}}{\partial \mathbf{x}} + \left[ \sum\_{j=1}^{m} \Gamma\_j \right] \mathbf{C} \, f \tag{3.4.3}$$

where,

Computational Modelling of Multi-Scale Non-Fickian Dispersion

over [0, *a*] are given by

, (3.3.21)

. (3.3.22)

, (3.3.23)

. (3.3.24)

( ) for a given number of

can be expanded by using

74 in Porous Media - An Approach Based on Stochastic Calculus

*i*

<sup>2</sup> ( ) *<sup>i</sup> i*

*i*

*Tan a*

*N*

2 2

*i*

 

and finally the basis functions of <sup>0</sup> *H* . Then the random field *d t*

simplified settings; as we have seen previously, *SCxt d t* ( ( , ) ( ))

*SCxt d t SCxt* 

*i i*

kernel (equation (3.3.10)) are the normalised eigenfunctions,

**3.4 The Dispersion and Travel Length Fluctuations** 

2 2 2 2

2 2

<sup>2</sup> <sup>2</sup>

  

*i i*

directly and to the length of the experiment, *a*, indirectly through

*j jj*

*jj j*

*db t*

 . (3.4.1)

. (3.4.2)

*i*

The orthonormal basis functions of the Hilbert space associated with the given exponential

<sup>1</sup> () ( ) ( ) *<sup>i</sup> <sup>i</sup> <sup>i</sup> <sup>i</sup> f x Sin x Cos x*

<sup>1</sup> <sup>1</sup> <sup>1</sup> <sup>1</sup> 1 (2 ) (2 ) 1 <sup>2</sup> <sup>4</sup> <sup>2</sup>

As seen from equation (3.3.21), the eigenvalues are dependent on the correlation length *b*

equation (3.3.22). Equation (3.3.22) is a transendental equation with roots *ωi*. Firstly, equation (3.3.22) has to be solved numerically for a given *a* and *b*, then the eigenvalues *λi* s

terms, *m,* can be evaluated by using equation (3.3.8) in conjunction with the standard

In this section, we explore the dispersion term in the SSTM to understand its behaviour in

( ( , ) ( )) ( ( , ) ( ) ( )) *m*

1 ( ( , ) ( )) ( ( , ) ) ( ) *m*

*f*

*j SCxt d t SCxt* 

1

*j*

Let us assume that we use the same realisation of the standard Wiener process across all *x* in evaluating this term, i.e., *dbj* is independent of x*.* Then we can express the dispersion term

*f x db t*

*N a Sin a Cos a*

The eigenvalues of equation (3.3.13) for constant <sup>2</sup>

*ωi* s are the roots of the following equation,

where <sup>1</sup>

where,

and the variance <sup>2</sup>

Wiener increments.

the spectral expansion,

for this situation,

*b* 

, and

$$\Phi\_{\rangle} = \frac{1}{2} \sqrt{\frac{\lambda\_{\rangle}}{N}} \left( -h\_{\text{x}} \text{Sim}(o\_{\rangle} \mathbf{x}) - \frac{h\_{\text{x}} \, o\_{\rangle}}{\theta} \text{Co}(o\_{\rangle} \mathbf{x}) \right);\tag{3.4.3a}$$

$$\Psi'\_{\neq} = \sqrt{\frac{\lambda\_j}{N}} \left( \left( \frac{o \rho\_j h\_x}{\theta} - 1 \right) \text{Sim}(o \rho\_j \mathbf{x}) - \left( o \rho\_j h\_x + \frac{o \rho\_j}{\theta} \right) \text{Co}(o \rho\_j \mathbf{x}) \right); \text{ and }\tag{3.4.3b}$$

$$\Gamma\_{j} = \sqrt{\frac{\lambda\_{j}}{N}} \left[ \left( \frac{\alpha \rho\_{j}^{2} \hbar\_{x}}{2} + \frac{\alpha \rho\_{j}^{2}}{\theta} \right) \text{Sim}(\alpha \rho\_{j} \mathbf{x}) - \left( \frac{\alpha \rho\_{j}^{3} \hbar\_{x}}{2 \theta} + \alpha \rho\_{j} \right) \text{Co}(\alpha \rho\_{j} \mathbf{x}) \right]. \tag{3.4.3c}$$

This exercise shows the complexity of the dispersion term and it is related to the local spatial gradient of the concentration, *<sup>C</sup> x* , the second order spatial derivative 2 2 *C x* , as well as *C*(*x,t*). *<sup>j</sup>* , *j* and *j* coefficients are also dependent on *λi* s, *hx*, and ωi s, and they are analogous to Fourier series, i.e., the coefficients themselves have the form of spectral expansions. Recall that *hx* is a scale length introduced to retain the second order derivative of fluxes in the SSTM. We can assume that ( ) *<sup>x</sup> h x* , i.e., 0 0 *<sup>x</sup> x <sup>h</sup> Lim x* where *x* is the discretization length in *x* direction. If we make *hx*=0 then *<sup>j</sup>* =0, which eliminates the second order spatial derivative of the concentration in equation (3.4.3). Even in a simplified setting of having a common standard Wiener increment (*dbi*(*t*)) for all *x*, and assuming a differentiable function for *C*(*x*,*t*), this analysis shows that the dispersion term depicts much more complicated behaviours that cannot be modelled by scale independent simplifying assumptions. Therefore, one should note that as time increases the effects of Wiener process increases, and it would be worthwhile to investigate the dependence of the random field *d t* ( ) for different scales of experiments. In other words, from equation (3.2.7), we see that,

$$dX(t) = \overline{V}(\mathbf{x}, t)dt + \xi(\mathbf{x}, t)dt,\tag{3.4.4}$$

where *X*(*t*) is travel distance of a solute particle in *x*-direction.

Figure 3.2 show the space-time grid for the explicit difference scheme that can be used to independently calculate the value at time *<sup>n</sup>* <sup>1</sup> *t* (denoted by ●) from the values at time *nt*

A Stochastic Model for Hydrodynamic Dispersion 77

Figure 3.2. An explicit space-time scheme used for the computational solution.

2

*k*

where *<sup>n</sup> Uk* = value of *U* at the grid point (*k, n*).

The second derivative can be given by,

The operator *S* can be written as,

The first derivative of a variable *U* can be described as (Morton and Mayers, 1994),

*k*

1

1 2

,

2 <sup>2</sup> 2

*x x* 

, (3.5.1)

. (3.5.2)

*n n n k k*

*U U U x x*

2 2 2 *<sup>n</sup> nn n kk k*

*Uh U <sup>x</sup> SU*

*U UUU x x* 

(denoted by ○).

As seen in equation (3.3.7), *d t x t dt* () ( ,) is the fluctuating part of the travel length of a solute within the time *dt* at a given *x*. Assuming *Vxt* ( ,) is a constant across the domain [0,*a*], by investigating the nature of the fluctuating component of the travel length as the scale of the experiments changes, we can understand the scale dependency of the dispersion term better. , )

*d t* ( ) is an irregular stochastic function of time with spatial component appearing in normalised eigen functions, ( ) *<sup>j</sup> f x* .A standard Wiener process increment (*dbi*) corresponding to a time increment, Δ*ti* , should be generated for each eigenfunction. *d t* ( ) can be approximated as a summation of *M* terms:

$$d\mathcal{J}(t) \approx \sum\_{i=1}^{M} \sqrt{\lambda\_i} \, f\_i(\mathbf{x}) db\_i(t) \,.$$

Let () () *<sup>i</sup> i i x fx* , then 1 ( ) () *M i i i d x db t* .

Recall ( ) *<sup>i</sup> db t* s are independent zero-mean Gaussian increments with Δ*ti* variance.

Therefore, [ ] ( ) 0 *E d E x db i i* , for a given *x*.

For a given *x*, <sup>2</sup> [ ] () [ ] *Var d x Var db <sup>i</sup> <sup>i</sup>* .

If we discretise the time axis equidistantly, 1 2 ... ... *<sup>i</sup> tt t t* , then,

$$\left[\operatorname{Var}[d\beta]\right] = \left(\sum \eta\_i^2(\mathbf{x})\right) \Delta t \,\,. \tag{3.4.5}$$

As seen from equation (3.4.5), *d* is a summation of independent Gaussian processes for a given *x* value making it a Gaussian process zero mean and variance proportional to Δ*t*. We will make use of equation (3.4.5) in the approximate numerical solution of SSTM for large scale experiments in chapter 4.

### **3.5 Numerical Solutions of the 1-D SSTM and Their Behaviours**

We solve equation (3.2.14) for strong solutions using a finite difference scheme which is based on Euler solution of Ito integral. One dimensional domain is discretised, and the basis

of numerical solutions to SPDEs as given by Gaines and Lyons (1997) is adopted. A constant mean velocity is assumed in the scheme. The differential operator *S* in equation (3.2.11) was expressed as a differential operator using a backward difference scheme. One dimensional spatial length, *a* ( 0*x a* ) on *x* axis was divided into (*k-*1) equidistant intervals of small lengths of *x* . The total model time, *T*, was divided into (*n-*1) equidistant small intervals of *t* . The space-time grid for the explicit difference scheme that can be used to independently calculate the concentration value at time level *<sup>n</sup>* <sup>1</sup> *t* from the concentration values at time intervals *nt* , thus preserving the non-antipating nature of Ito integral.

, )

is the fluctuating part of the travel length of a

*<sup>i</sup>* . (3.4.5)

is a summation of independent Gaussian processes for a

76 in Porous Media - An Approach Based on Stochastic Calculus

solute within the time *dt* at a given *x*. Assuming *Vxt* ( ,) is a constant across the domain [0,*a*], by investigating the nature of the fluctuating component of the travel length as the scale of the experiments changes, we can understand the scale dependency of the dispersion

( ) is an irregular stochastic function of time with spatial component appearing in normalised eigen functions, ( ) *<sup>j</sup> f x* .A standard Wiener process increment (*dbi*) corresponding to a time increment, Δ*ti* , should be generated for each eigenfunction.

() ( ) ()

 

 <sup>2</sup> [ ] () *Var d x t* 

given *x* value making it a Gaussian process zero mean and variance proportional to Δ*t*. We will make use of equation (3.4.5) in the approximate numerical solution of SSTM for large

We solve equation (3.2.14) for strong solutions using a finite difference scheme which is based on Euler solution of Ito integral. One dimensional domain is discretised, and the basis of numerical solutions to SPDEs as given by Gaines and Lyons (1997) is adopted. A constant mean velocity is assumed in the scheme. The differential operator *S* in equation (3.2.11) was expressed as a differential operator using a backward difference scheme. One dimensional spatial length, *a* ( 0*x a* ) on *x* axis was divided into (*k-*1) equidistant intervals of small lengths of *x* . The total model time, *T*, was divided into (*n-*1) equidistant small intervals of *t* . The space-time grid for the explicit difference scheme that can be used to independently calculate the concentration value at time level *<sup>n</sup>* <sup>1</sup> *t* from the concentration values at time

 

*ii i*

.

1

*i d t f x db t*

Recall ( ) *<sup>i</sup> db t* s are independent zero-mean Gaussian increments with Δ*ti* variance.

( ) ()

*i i*

If we discretise the time axis equidistantly, 1 2 ... ... *<sup>i</sup> tt t t* , then,

.

*M*

As seen in equation (3.3.7), *d t x t dt*

term better.

Let () () *i i i*

 

*x fx* , then 1

 

> 

For a given *x*, <sup>2</sup> [ ] () [ ] *Var d x Var db* 

Therefore, [ ] ( ) 0 *E d E x db* 

As seen from equation (3.4.5), *d*

scale experiments in chapter 4.

*d t* 

*d t* 

( ) can be approximated as a summation of *M* terms:

*M*

*i d x db t*

*i i* , for a given *x*.

*<sup>i</sup> <sup>i</sup>* .

**3.5 Numerical Solutions of the 1-D SSTM and Their Behaviours**

intervals *nt* , thus preserving the non-antipating nature of Ito integral.

 

() ( ,) 

Figure 3.2 show the space-time grid for the explicit difference scheme that can be used to independently calculate the value at time *<sup>n</sup>* <sup>1</sup> *t* (denoted by ●) from the values at time *nt* (denoted by ○).

Figure 3.2. An explicit space-time scheme used for the computational solution. The first derivative of a variable *U* can be described as (Morton and Mayers, 1994),

$$\left(\frac{\partial \mathcal{U}\mathcal{U}}{\partial \mathbf{x}}\right)\_k^n = \frac{\mathcal{U}\mathcal{U}\_k^n - \mathcal{U}\_{k-1}^n}{\Delta \mathbf{x}}\,,\tag{3.5.1}$$

where *<sup>n</sup> Uk* = value of *U* at the grid point (*k, n*).

The second derivative can be given by,

$$\left(\frac{\partial^2 \mathcal{U}}{\partial \mathbf{x}^2}\right)\_k^n = \frac{\mathcal{U}I\_k^n - 2\mathcal{U}I\_{k-1}^n + \mathcal{U}I\_{k-2}^n}{\Delta \mathbf{x}^2}.\tag{3.5.2}$$

The operator *S* can be written as,

$$SUI = -\left(\frac{\partial \mathcal{U}}{\partial \mathcal{X}} + \frac{h\_x}{2} \frac{\partial^2 \mathcal{U}}{\partial \mathcal{X}^2}\right) \prime$$

[0, 1], the first 30 values of the roots (*ωi* ) are generally sufficient. We generate the standard Wiener process increments for 0.001 day time intervals for total of three days. Then the

A Stochastic Model for Hydrodynamic Dispersion 79

to compute *d*(t) , the Hilbert space valued Wiener incremental processes, using the *KL* expansion (equation (3.3.7)). Then we calculate the concentration profile for the discretised values of spatial-temporal development for the mean velocity of 0.5 m/day. The numerical solution is implemented in a mathematical software package, Mathematica (Wolfram

We use a spatial grid length of 0.1 m for the numerical calculation. The initial concentration distribution profile of 1.0 unit at *x* = 0 is considered and it exponentially decreases through the rest of the domain according to the function, *e*-*5k x* , where *k* = 1, 2, …,10 and *x* = grid size. We begin the numerical scheme with very small numerical concentration values, rather than zero concentrations, to reduce the numerical errors at the beginning of the scheme. The concentration of 1.0 unit is maintained at the boundary of *x* = 0 for the whole time period of

To investigate the general behaviour of the SSTM, we obtain the temporal development of the concentration profiles at the mid point of the domain, *x* = 0.5 m, for various parameter

and constant mean velocity of 0.5 m/day are used for all the experiments, so that we would

Figure 3.3. Comparison of deterministic advection-dispersion (*D* = 0.01) and stochastic ( <sup>2</sup>

= 0.001 and *b* = 0.0001) model concentration profiles. An explicit space-time scheme used for

Figure 3.3 shows that the stochastic model can mimic the solution of the advectiondispersion equation with reasonable accuracy. The concentration breakthrough curves (the

the deterministic curve for the advection-dispersion equation for dispersion coefficient (*D*)

= 0.001 and *b* = 0.0001, and

time history of concentration at a fixed *x* ) for the SSTM for <sup>2</sup>

. The same realisation of the standard Wiener process increments

*<sup>n</sup>* , we calculate the basis functions using equation (3.3.23). These values are used

using equation (3.3.21). With these roots,

*n* are computed for the given <sup>2</sup>

the solution to mimic a continuous point source.

eigenvalues

Research, 1999).

combinations of *b* and <sup>2</sup>

not bias our comparisons.

the computational solution.

 and  Therefore, we can express the operator in the difference form

$$\left(\left(SL\right)\_{k}^{\boldsymbol{u}} = -\left(\left(\frac{\partial \boldsymbol{\mathcal{U}}}{\partial \boldsymbol{\mathfrak{x}}}\right)\_{k}^{\boldsymbol{u}} + \frac{h\_{\boldsymbol{x}}}{2} \left(\frac{\partial^{2} \boldsymbol{\mathcal{U}}}{\partial \boldsymbol{\mathfrak{x}}^{2}}\right)\_{k}^{\boldsymbol{u}}\right). \tag{3.5.3}$$

 Substituting the backward difference schemes from (3.5.1) and (3.5.3) and taking *<sup>x</sup> h x* ,

$$\left(\left(SLI\right)\_{k}^{u} = -\left(\frac{1}{2\Delta x}\right)\left[\Im \mathcal{U}\_{k}^{u} - \mathbf{4} \mathcal{U}\_{k-1}^{u} + \mathcal{U}\_{k-2}^{u}\right].\tag{3.5.4}$$

The first derivative of *U* with respect to time can be expressed using a forward difference scheme,

$$\frac{\partial \mathcal{L}I}{\partial t} = \frac{\mathcal{L}I\_k^{n+1} - \mathcal{U}\_k^n}{\Delta t} \,. \tag{3.5.5}$$

Applying equation (3.5.1) and (3.5.3) to (3.5.5) and considering the mean velocity for the region as a constant, *v*, we can obtain the following scheme:

$$\begin{split} \mathbf{C}\_{k}^{u+1} &= \mathbf{C}\_{k}^{u} - \left(\frac{\Delta t.\nu}{2\Delta x}\right) \mathbf{\*} \left[\mathfrak{R}\mathbf{C}\_{k}^{u} - \mathfrak{A}\mathbf{C}\_{k-1}^{u} + \mathbf{C}\_{k-2}^{u}\right] \\ &- \left(\frac{1}{2\Delta x}\right) \mathbf{\*} \left[\mathfrak{R}\mathbf{C}\_{k}^{u} d\mathcal{J}\boldsymbol{\theta}(\mathbf{t})\_{k}^{u} - \mathbf{4}\mathbf{C}\_{k-1}^{u} d\mathcal{J}\boldsymbol{\theta}(\mathbf{t})\_{k-1}^{u} + \mathbf{C}\_{k-2}^{u} d\mathcal{J}\boldsymbol{\theta}(\mathbf{t})\_{k-2}^{u}\right], \end{split} \tag{3.5.6}$$

where *d t* ( ) = Wiener process increments in Hilbert space for a given *x*.

The explicit difference model (3.5.6) gives the future values of a stochastic variable in terms of the past values preserving the properties of Ito definition of integration with respect to time. This scheme is stable and gives strong solutions to equation (3.2.14), as in many SPDEs, if <sup>4</sup> *t* 10 . For example, if *a=*1000 meters, and if we simulate the solute transport for at least 1000 days, taking *x* =0.01 m and *t* =0.0001 days for stability reasons, we need a grid of 5 7 10 10 *x* . In addition, the evaluation of *d t* ( ) for each *x* involves the summation of a large number of ( ) *<sup>i</sup> x dbi* terms. Because of the computational time it requires to solve the SPDE, we use the scheme given in equation (3.5.6) when *a<* 10 m, and for larger *a* values we approximate *d t* ( ) term as described later.

We can now investigate the behaviour of the stochastic solute transport model (SSTM) in one-dimension. The main parameters of the SSTM are correlation length, *<sup>b</sup>* and variance, 2 . As the statistical nature of the computational solution changes with different *b* and <sup>2</sup> , we would like to understand the effect of these parameters on the solution of the model. Furthermore, we attempt to understand these parameters in relation to the hydrodynamic dispersion.

The finite difference numerical schemes are used in the investigation taking the numerical convergence and stability into account for the domain [0, 1]. First we solve equation (3.3.20) for the given values of *b* and <sup>2</sup> **.** It is necessary to find and appropriate number of roots for equation (3.3.20) to produce the desired accuracy of the numerical solution. For the domain

. (3.5.4)

. (3.5.5)

. (3.5.3)

(3.5.6)

,

78 in Porous Media - An Approach Based on Stochastic Calculus

<sup>2</sup> 2 *n n*

*k k*

*x x* 

*x*

<sup>2</sup>

*n x*

*UhU SU*

Substituting the backward difference schemes from (3.5.1) and (3.5.3) and taking *<sup>x</sup> h x* ,

 1 2 <sup>1</sup> 3 4

*n nn n kk k <sup>k</sup> SU UUU*

The first derivative of *U* with respect to time can be expressed using a forward difference

Applying equation (3.5.1) and (3.5.3) to (3.5.5) and considering the mean velocity for the

*n n* <sup>1</sup> *U U U k k t t*

1 2

*Cd t C d t C d t*

*nnn n n n k kk k k k*

*x dbi* terms. Because of the computational time it requires to solve

days

**.** It is necessary to find and appropriate number of roots for

<sup>1</sup> \* 3 () 4 () () , <sup>2</sup>

The explicit difference model (3.5.6) gives the future values of a stochastic variable in terms of the past values preserving the properties of Ito definition of integration with respect to time. This scheme is stable and gives strong solutions to equation (3.2.14), as in many SPDEs, if <sup>4</sup> *t* 10 . For example, if *a=*1000 meters, and if we simulate the solute transport for at least 1000 days, taking *x* =0.01 m and *t* =0.0001 days for stability reasons, we need

the SPDE, we use the scheme given in equation (3.5.6) when *a<* 10 m, and for larger *a* values

We can now investigate the behaviour of the stochastic solute transport model (SSTM) in one-dimension. The main parameters of the SSTM are correlation length, *<sup>b</sup>* and variance, 2

. As the statistical nature of the computational solution changes with different *b* and <sup>2</sup>

we would like to understand the effect of these parameters on the solution of the model. Furthermore, we attempt to understand these parameters in relation to the hydrodynamic

The finite difference numerical schemes are used in the investigation taking the numerical convergence and stability into account for the domain [0, 1]. First we solve equation (3.3.20)

equation (3.3.20) to produce the desired accuracy of the numerical solution. For the domain

1 1 2 2

( ) for each *x* involves the summation

Therefore, we can express the operator in the difference form

region as a constant, *v*, we can obtain the following scheme:

2

*x*

a grid of 5 7 10 10 *x* . In addition, the evaluation of *d t*

( ) term as described later.

*i* 

1

*k*

2

. \*3 4

( ) = Wiener process increments in Hilbert space for a given *x*.

*n n nn n k k kk k*

*t v C C CCC x*

scheme,

where *d t* 

of a large number of ( )

for the given values of *b* and <sup>2</sup>

we approximate *d t*

dispersion.

[0, 1], the first 30 values of the roots (*ωi* ) are generally sufficient. We generate the standard Wiener process increments for 0.001 day time intervals for total of three days. Then the eigenvalues *n* are computed for the given <sup>2</sup> using equation (3.3.21). With these roots, and *<sup>n</sup>* , we calculate the basis functions using equation (3.3.23). These values are used to compute *d*(t) , the Hilbert space valued Wiener incremental processes, using the *KL* expansion (equation (3.3.7)). Then we calculate the concentration profile for the discretised values of spatial-temporal development for the mean velocity of 0.5 m/day. The numerical solution is implemented in a mathematical software package, Mathematica (Wolfram Research, 1999).

We use a spatial grid length of 0.1 m for the numerical calculation. The initial concentration distribution profile of 1.0 unit at *x* = 0 is considered and it exponentially decreases through the rest of the domain according to the function, *e*-*5k x* , where *k* = 1, 2, …,10 and *x* = grid size. We begin the numerical scheme with very small numerical concentration values, rather than zero concentrations, to reduce the numerical errors at the beginning of the scheme. The concentration of 1.0 unit is maintained at the boundary of *x* = 0 for the whole time period of the solution to mimic a continuous point source.

To investigate the general behaviour of the SSTM, we obtain the temporal development of the concentration profiles at the mid point of the domain, *x* = 0.5 m, for various parameter combinations of *b* and <sup>2</sup> . The same realisation of the standard Wiener process increments and constant mean velocity of 0.5 m/day are used for all the experiments, so that we would not bias our comparisons.

 Figure 3.3. Comparison of deterministic advection-dispersion (*D* = 0.01) and stochastic ( <sup>2</sup> = 0.001 and *b* = 0.0001) model concentration profiles. An explicit space-time scheme used for the computational solution.

Figure 3.3 shows that the stochastic model can mimic the solution of the advectiondispersion equation with reasonable accuracy. The concentration breakthrough curves (the time history of concentration at a fixed *x* ) for the SSTM for <sup>2</sup> = 0.001 and *b* = 0.0001, and the deterministic curve for the advection-dispersion equation for dispersion coefficient (*D*)

individual concentration profiles have worse fluctuatiing stochasticity, and there are significant differences between concentration values for different *b* values at a given time. The high values of the variance not only directly increase the unpredictable nature of the flow but also influence the ways in which *b* affects the flow. We also observe that with high *b* values the asymptotic values (sills) of the concentration profiles are lower than the

A Stochastic Model for Hydrodynamic Dispersion 81

In the note that when *b* is very small, the concentration profile is smooth, but when *b* is 0.1 m

Wiener process increments, we obtain the break through curves as shown in Figure 3.6. The flow tends to be significantly unsteady for larger correlation lengths and still shows smaller

0.5 <sup>1</sup> 1.5 <sup>2</sup> 2.5 <sup>3</sup> <sup>t</sup> (days)

0.5 1 1.5 2 2.5 3

= 0.01.

= 0.1.

by 10 times, to 0.01, and by using the same standard

t (days)

intensities the fluctuation

b0.25 b0.1 b0.01 b0.001 b0.0001

b0.25

b0.1 b0.01

b0.001 b0.0001

fluctuations smaller *b* values. Furthermore larger values of <sup>2</sup>

Figure 3.6. Concentration profiles at *x* = 0.5 m for <sup>2</sup>

Figure 3.7. Concentration profiles at *x* = 0.5 m for <sup>2</sup>

deterministic sill.

it lowers the sill. By increasing <sup>2</sup>

the effect of *b* significantly.

C

0.2

C

0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.4

0.6

0.8

of 0.01 m2/day are overlaid in Figure 3.3. We can always find a solution for the SSTM that reasonably represents the deterministic break through curve for any given dispersion coefficient using appropriate values for the parameters, <sup>2</sup> and *b*.

To study the influences of *b* and <sup>2</sup> on the solution of the problem, we keep one parameter constant and change the other within a reasonable range to examine the behavioural change of the concentration breakthrough curves. Figure 3.4 shows the concentration profile at *x* = 0.5 m for a small value of the variance, <sup>2</sup> = 0.0001, when the correlation length, *b* varies from 0.0001m to 0.25m. Although the range of *b* varies from 0.0001 to 0.25m (a change of 2500 times) the change of stochasticity (noise level) is negligible and the solutions of the SSTM are independent of *b* and behave like those of a deterministic model.

Figure 3.4. Concentration profiles at *x* = 0.5m for <sup>2</sup> = 0.0001.

Figure 3.5. Concentration profiles at *x* = 0.5 m for <sup>2</sup> = 0.001.

We gradually increase <sup>2</sup> and obtain the concentration profiles for the same regime of *b* (0.0001 m to 0.25 m) to examine the effect of <sup>2</sup> . With an increase in <sup>2</sup> by 10 times, Figure 3.5 shows that two types of changes have occurred in the concentration profiles;

and *b*.

on the solution of the problem, we keep one parameter

t (days)

t (days)

. With an increase in <sup>2</sup>

b 0.25 b 0.1 b 0.01 b 0.001 b 0.0001

b0.25 b0.1 b0.01 b0.001 b0.0001

by 10 times,

= 0.0001, when the correlation length, *b* varies

80 in Porous Media - An Approach Based on Stochastic Calculus

of 0.01 m2/day are overlaid in Figure 3.3. We can always find a solution for the SSTM that reasonably represents the deterministic break through curve for any given dispersion

constant and change the other within a reasonable range to examine the behavioural change of the concentration breakthrough curves. Figure 3.4 shows the concentration profile at *x* =

from 0.0001m to 0.25m. Although the range of *b* varies from 0.0001 to 0.25m (a change of 2500 times) the change of stochasticity (noise level) is negligible and the solutions of the

SSTM are independent of *b* and behave like those of a deterministic model.

0.5 1 1.5 2 2.5 3

0.5 1 1.5 2 2.5 3

Figure 3.5 shows that two types of changes have occurred in the concentration profiles;

= 0.001.

and obtain the concentration profiles for the same regime of *b*

= 0.0001.

Figure 3.4. Concentration profiles at *x* = 0.5m for <sup>2</sup>

Figure 3.5. Concentration profiles at *x* = 0.5 m for <sup>2</sup>

(0.0001 m to 0.25 m) to examine the effect of <sup>2</sup>

coefficient using appropriate values for the parameters, <sup>2</sup>

To study the influences of *b* and <sup>2</sup>

0.2

0.2

We gradually increase <sup>2</sup>

0.4

0.6

0.8

1 C

0.4

0.6

0.8

1 C

0.5 m for a small value of the variance, <sup>2</sup>

individual concentration profiles have worse fluctuatiing stochasticity, and there are significant differences between concentration values for different *b* values at a given time. The high values of the variance not only directly increase the unpredictable nature of the flow but also influence the ways in which *b* affects the flow. We also observe that with high *b* values the asymptotic values (sills) of the concentration profiles are lower than the deterministic sill.

In the note that when *b* is very small, the concentration profile is smooth, but when *b* is 0.1 m it lowers the sill. By increasing <sup>2</sup> by 10 times, to 0.01, and by using the same standard Wiener process increments, we obtain the break through curves as shown in Figure 3.6. The flow tends to be significantly unsteady for larger correlation lengths and still shows smaller fluctuations smaller *b* values. Furthermore larger values of <sup>2</sup> intensities the fluctuation the effect of *b* significantly.

Figure 3.6. Concentration profiles at *x* = 0.5 m for <sup>2</sup> = 0.01.

Figure 3.7. Concentration profiles at *x* = 0.5 m for <sup>2</sup> = 0.1.

We increase *b* by 10 times to obtain Figure 3.10, which shows considerable changes in the

A Stochastic Model for Hydrodynamic Dispersion 83

To understand the effects of different Wiener realizations on the concentration profiles by using 50 different Wiener realisations to calculate the 95% confidence intervals. They show

fluctuation regimes increase but the solutions remain stable. Obviously, the confidence

concentration profile are negligible, but for larger values (for example, <sup>2</sup>

0.5 1 1.5 2 2.5 3

for different model parameters of the correlation length, *b* and the variance, <sup>2</sup>

We can now extend the investigation to a 10 m flow length. Similar to the 1 m investigation, a simple setting of the one-dimensional case is used to explore the behaviour of the model

meaningful concentration profiles for the extended spatial domain, we also increase the time interval of the solution. The stochastic model is computationally solved for an average linear velocity of 0.5 m/day for 30 days, which allows sufficient time for the solute to travel the entire spatial domain. However, due to increases in the time length, a different realisation of the standard Wiener process increments was used for the 10 m length. This may raise the question of the validity in comparing of the model for two different spatial lengths. However, as we have shown earlier the model is reasonably stable for different Wiener process increments and therefore, it is reasonable to assume that the solutions of the model are not significantly affected by the different Wiener processes. A single realisation of the standard Wiener process is used throughout the investigation of the 10 m spatial domain. Similar to 1 m case, the 10 m scale shows that a SSTM could mimic the solution of the

Figure 3.10. Concentration profiles at *x* = 0.5 m for *b* = 0.001.

advection – dispersion equation for the same distance (Figure 3.11)

= 0.001, *b*= 0.01), the variations in the

t (days)

 2 0.25

. To obtain

 2 0.1

 2 0.01

 2 0.001

 2 0.0001

= 0.1, *b*= 0.1) the

breakthrough curves.

0.2

0.4

0.6

0.8

1 C

that, for smaller values of parameters (for example, <sup>2</sup>

intervals widen with the larger values of the parameters.

The concentration profiles for higher *b* regimes for the increased value of <sup>2</sup> (=0.1) are highly random as shown in Figure 3.7. With higher values of *b*, the concentration profiles become highly irregular making the numerical scheme unstable. Therefore, limit our experiments to smaller *b* values, that are less than 0.01 m. Figure 4.8 shows that the fluctuation invariably increase with the high <sup>2</sup> values and the behaviour of the model continues with the same trend that we noticed earlier, but with enhanced effects.

Figure 3.8. Concentration profiles at *x* = 0.5 m for <sup>2</sup> = 0.25.

Figure 3.9 shows the concentration profiles at *b* = 0.0001 for the range of <sup>2</sup> that varies from 0.0001 to 0.25. By comparing in Figure 3.4, when <sup>2</sup> was very small, the fluctuation are not distinguishable even for very high *b* values; however, in Figure 3.9, irrespective of smaller *<sup>b</sup>*, 2 influenced the behaviour of the flow. It is not possible to differentiate the concentration profiles for very small <sup>2</sup> , such as 0.0001 and 0.001. With the increase of <sup>2</sup> stochasticity increases rapidly, and <sup>2</sup> influences the behaviour of the flow more than *b* does.

Figure 3.9. Concentration profiles at *x* = 0.5 m for *b* = 0.0001.

t (days)

= 0.25.

, such as 0.0001 and 0.001. With the increase of <sup>2</sup>

influences the behaviour of the flow more than *b* does.

distinguishable even for very high *b* values; however, in Figure 3.9, irrespective of smaller *<sup>b</sup>*, 2

influenced the behaviour of the flow. It is not possible to differentiate the concentration

b0.005

that varies from

stochasticity

was very small, the fluctuation are not

t (days)

 2 0.25

 2 0.1

2 0.01

 2 0.001

 2 0.0001

b0.001

b0.0001

values and the behaviour of the model

(=0.1) are

82 in Porous Media - An Approach Based on Stochastic Calculus

highly random as shown in Figure 3.7. With higher values of *b*, the concentration profiles become highly irregular making the numerical scheme unstable. Therefore, limit our experiments to smaller *b* values, that are less than 0.01 m. Figure 4.8 shows that the

The concentration profiles for higher *b* regimes for the increased value of <sup>2</sup>

continues with the same trend that we noticed earlier, but with enhanced effects.

0.5 1 1.5 2 2.5 3

Figure 3.9 shows the concentration profiles at *b* = 0.0001 for the range of <sup>2</sup>

0.5 1 1.5 2 2.5 3

Figure 3.9. Concentration profiles at *x* = 0.5 m for *b* = 0.0001.

Figure 3.8. Concentration profiles at *x* = 0.5 m for <sup>2</sup>

0.0001 to 0.25. By comparing in Figure 3.4, when <sup>2</sup>

fluctuation invariably increase with the high <sup>2</sup>

0.2 0.3 0.4 0.5 0.6 0.7

0.2

profiles for very small <sup>2</sup>

increases rapidly, and <sup>2</sup>

0.2

0.4

0.6

0.8

1 C

C

We increase *b* by 10 times to obtain Figure 3.10, which shows considerable changes in the breakthrough curves.

To understand the effects of different Wiener realizations on the concentration profiles by using 50 different Wiener realisations to calculate the 95% confidence intervals. They show that, for smaller values of parameters (for example, <sup>2</sup> = 0.001, *b*= 0.01), the variations in the concentration profile are negligible, but for larger values (for example, <sup>2</sup> = 0.1, *b*= 0.1) the fluctuation regimes increase but the solutions remain stable. Obviously, the confidence intervals widen with the larger values of the parameters.

Figure 3.10. Concentration profiles at *x* = 0.5 m for *b* = 0.001.

We can now extend the investigation to a 10 m flow length. Similar to the 1 m investigation, a simple setting of the one-dimensional case is used to explore the behaviour of the model for different model parameters of the correlation length, *b* and the variance, <sup>2</sup> . To obtain meaningful concentration profiles for the extended spatial domain, we also increase the time interval of the solution. The stochastic model is computationally solved for an average linear velocity of 0.5 m/day for 30 days, which allows sufficient time for the solute to travel the entire spatial domain. However, due to increases in the time length, a different realisation of the standard Wiener process increments was used for the 10 m length. This may raise the question of the validity in comparing of the model for two different spatial lengths. However, as we have shown earlier the model is reasonably stable for different Wiener process increments and therefore, it is reasonable to assume that the solutions of the model are not significantly affected by the different Wiener processes. A single realisation of the standard Wiener process is used throughout the investigation of the 10 m spatial domain.

Similar to 1 m case, the 10 m scale shows that a SSTM could mimic the solution of the advection – dispersion equation for the same distance (Figure 3.11)

pores; and *b* is indicative of the size of pores, and, may be, geometric shapes of pores. <sup>2</sup>

A Stochastic Model for Hydrodynamic Dispersion 85

affects the profiles more dramatically, especially depressing asymptotic of the profiles, indicating that solute mass is dissolved in a larger volume of water. This alludes again to pore geometry (shapes and interconnecting paths) and if pore structure is heterogeneous

flexibility of defining the nature of solute dispersion, and the complex interaction between 2

The Lincoln University aquifer is 9.49 m long, 4.66 m wide and 2.6 m deep. As shown in Figure 3.13, constant head tanks are the boundaries of the aquifer at its upstream and downstream ends. A porous wall provides the hydraulic connection between the aquifer and the head tanks. A weir controls the water surface elevation in each head tank, and each weir can be adjusted to provide different hydraulic gradients. However, a uniform hydraulic gradient of 0.0018 (head loss along the aquifer / flow length = 0.017 m / 9.49 m) is

Figure 3.13. Schematic diagram of the artificial aquifer at Lincoln University, New Zealand

(Courtesy of Dr.John Bright, Lincoln Ventures Ltd).

maintained during the entire experiment along the longitudinal direction of the tank.

and *b* would help us to characterise the pore structures for a given velocity kernel.

as well as b to be high. <sup>2</sup>

**3.6 A Comparison of the SSTM with the Experiments Data** 

with high porosity, one could expect <sup>2</sup>

and *b* allow us more

Figure 3.11. Comparison of deterministic advection-dispersion (*D* = 0.035 m2/day) and the SSTM ( <sup>2</sup> = 0.001 and *b* = 0.0001) concentration profiles at *x* = 5m for 10 m domain.

The behaviours of the SSTM model for the 10-m flow length are quite similar to those for the 1-m flow length. The influence of the parameter *b* can be seen in Figure 3.12 for the fixed values of <sup>2</sup> . High *b* values increase the propensity of the flow to be more stochastic and decrease the asymptotic values of the concentration. For a given *<sup>b</sup>*, the effects of increasing 2 are more profound in comparison to those associated with increasing *b* for a fixed <sup>2</sup> . be

Figure 3.12. Comparison of deterministic advection-dispersion (*D* = 0.035 m2/day) and the SSTM ( <sup>2</sup> = 0.001 and *b* = 0.0001) concentration profiles at *x* = 5m for 10 m domain.

How do we relate the parameter, <sup>2</sup> and *b*, to the physical porous structure? The relationship need to be understood through the influences on the concentration profiles. *b* is the correlation length of the velocity Kernel 1 2 *qx x* (,) (see equation 3.3.9), and higher the *b* slower the rate at which 1 2 *qx x* (,) decays. This means that pore structure contains larger

.

84 in Porous Media - An Approach Based on Stochastic Calculus

Figure 3.11. Comparison of deterministic advection-dispersion (*D* = 0.035 m2/day) and the

The behaviours of the SSTM model for the 10-m flow length are quite similar to those for the 1-m flow length. The influence of the parameter *b* can be seen in Figure 3.12 for the fixed

decrease the asymptotic values of the concentration. For a given *<sup>b</sup>*, the effects of increasing 2

are more profound in comparison to those associated with increasing *b* for a fixed <sup>2</sup>

Figure 3.12. Comparison of deterministic advection-dispersion (*D* = 0.035 m2/day) and the

relationship need to be understood through the influences on the concentration profiles. *b* is the correlation length of the velocity Kernel 1 2 *qx x* (,) (see equation 3.3.9), and higher the *b* slower the rate at which 1 2 *qx x* (,) decays. This means that pore structure contains larger

= 0.001 and *b* = 0.0001) concentration profiles at *x* = 5m for 10 m domain.

and *b*, to the physical porous structure? The

= 0.001 and *b* = 0.0001) concentration profiles at *x* = 5m for 10 m domain.

. High *b* values increase the propensity of the flow to be more stochastic and

SSTM ( <sup>2</sup> 

values of <sup>2</sup>

SSTM ( <sup>2</sup> 

How do we relate the parameter, <sup>2</sup>

pores; and *b* is indicative of the size of pores, and, may be, geometric shapes of pores. <sup>2</sup> affects the profiles more dramatically, especially depressing asymptotic of the profiles, indicating that solute mass is dissolved in a larger volume of water. This alludes again to pore geometry (shapes and interconnecting paths) and if pore structure is heterogeneous with high porosity, one could expect <sup>2</sup> as well as b to be high. <sup>2</sup> and *b* allow us more flexibility of defining the nature of solute dispersion, and the complex interaction between 2 and *b* would help us to characterise the pore structures for a given velocity kernel.

### **3.6 A Comparison of the SSTM with the Experiments Data**

The Lincoln University aquifer is 9.49 m long, 4.66 m wide and 2.6 m deep. As shown in Figure 3.13, constant head tanks are the boundaries of the aquifer at its upstream and downstream ends. A porous wall provides the hydraulic connection between the aquifer and the head tanks. A weir controls the water surface elevation in each head tank, and each weir can be adjusted to provide different hydraulic gradients. However, a uniform hydraulic gradient of 0.0018 (head loss along the aquifer / flow length = 0.017 m / 9.49 m) is maintained during the entire experiment along the longitudinal direction of the tank.

Figure 3.13. Schematic diagram of the artificial aquifer at Lincoln University, New Zealand (Courtesy of Dr.John Bright, Lincoln Ventures Ltd).

Figure 3.15. Concentration profiles of deterministic advection-dispersion model with D = 0.15

A Stochastic Model for Hydrodynamic Dispersion 87

Subsequently, we develop a one-dimensional deterministic advection-dispersion model using the *D* obtained from a two-dimensional comparison. We then use a similar curve fitting technique with the 1-D deterministic model and 1-D stochastic model to find the most

Figures 3.15 show the best fitting curves for both the SSTM and the 1-D advectiondispersion model. Having determined the appropriate parameters of the SSTM ( <sup>2</sup>

and *b* = 0.01) that simulate the Lincoln University aquifer at the spatial location under consideration (Row 5 – well A) we investigated the robustness of the model for different Wiener processes and found that the SSTM is reasonably stable for the seven different

Even though the results show, as mentioned above, that the parameter combination of <sup>2</sup>

0.01 and *b* = 0.01 is a fairly accurate representation of the experimental aquifer for the given spatial point, we need to test these parameters other spatial locations. Figure 3.16 shows that the 2-D deterministic advection-dispersion model with the longitudinal dispersion coefficient of 0.15 m2/day fits the aquifer data reasonably well for the similar locations.

Figure 3.16. Concentration profiles of 2-D deterministic advection-dispersion model (*D* =

t (days)

Det

Aquifer

= 0.01

> =

m2/day and SSTM with <sup>2</sup>

Wiener realisations tested.

0.2

0.15 m2/day) and the experimental aquifer.

0.4

0.6

0.8

1 C

suitable <sup>2</sup> 

and *b* for the SSTM.

= 0.01 and *b* = 0.01.

2 4 6 8 10

Multi-port monitoring wells are laid out on a 1m x 1m grid. The computer controlled peristaltic pumps enable fully automated, simultaneous solute water samples to be collected from sample points that are uniformly distributed throughout the aquifer (four sample points for each grid point at 0.4m, 1.0m, 1.6m and 2.2m depth from the top surface of the aquifer). The tracer used is Rhodamine WT (RWT) dye with an initial concentration of 200 parts per million and then allowed to decrease exponentially. Tracer is injected at the middle of the header tank using an injection box (dimensions of 50 cm length, 10 cm width and 20 cm depth). This tracer is rapidly mixed in the upstream header tank and, thus, infiltrates across the whole of the upstream face of the aquifer. This particular experiment described here lasted 432 hours, and two samples were taken at four-hour intervals from the wells.

Since, STTM described in this chapter is a one-dimensional model, we experiment in directly relating the one-dimension solute concentration profiles of the aquifer. However, as one can assume, the actual aquifer is subjected to transverse dispersion, and consideration of only a one-dimensional flow is not sufficiently accurate. Hence, we employ the following methodology to approximate the aquifer parameters.

Solute concentration values for the artificial aquifer are available for a large number of spatial points at different temporal intervals. The data are available mainly for header tank, row 1, row 3, row 5, row 7 and row 9 (see Figure 3.13) at all levels. Initially, we select a few spatial coordinates at row 5 of well A – level YE. We then develop a two-dimensional deterministic advection-dispersion transport model and obtain corresponding concentration values of the model that are similar to the selected spatial locations of the aquifer. As past studies show, we approximate that the transverse dispersion coefficient is 10% of the longitudinal dispersion (Fetter, 1999). The mean velocity is 0.5 m/day. The profiles of both the aquifer and the deterministic model are plotted in one axis system, to compare their similarities. This trial-and-error curve fitting technique is carried out to determine the most accurate dispersion coefficients of the deterministic model. In this procedure concentration values of the aquifer are normalised (i.e., the values vary from 0 to 1).

By trial and error, we find that the closest fit is given by the longitudinal dispersion coefficient of 0.15 m2/day, i.e., transverse dispersion is 0.015 m2/day, (Figure 3.14) for the aquifer.

Figure 3.14. Concentration profile of trial and error curve fit for D = 0.15 m2/day of the advection dispersion model with row 5 of the aquifer data.

86 in Porous Media - An Approach Based on Stochastic Calculus

Multi-port monitoring wells are laid out on a 1m x 1m grid. The computer controlled peristaltic pumps enable fully automated, simultaneous solute water samples to be collected from sample points that are uniformly distributed throughout the aquifer (four sample points for each grid point at 0.4m, 1.0m, 1.6m and 2.2m depth from the top surface of the aquifer). The tracer used is Rhodamine WT (RWT) dye with an initial concentration of 200 parts per million and then allowed to decrease exponentially. Tracer is injected at the middle of the header tank using an injection box (dimensions of 50 cm length, 10 cm width and 20 cm depth). This tracer is rapidly mixed in the upstream header tank and, thus, infiltrates across the whole of the upstream face of the aquifer. This particular experiment described here lasted 432 hours, and two samples were taken at four-hour intervals from the wells.

Since, STTM described in this chapter is a one-dimensional model, we experiment in directly relating the one-dimension solute concentration profiles of the aquifer. However, as one can assume, the actual aquifer is subjected to transverse dispersion, and consideration of only a one-dimensional flow is not sufficiently accurate. Hence, we employ the following

Solute concentration values for the artificial aquifer are available for a large number of spatial points at different temporal intervals. The data are available mainly for header tank, row 1, row 3, row 5, row 7 and row 9 (see Figure 3.13) at all levels. Initially, we select a few spatial coordinates at row 5 of well A – level YE. We then develop a two-dimensional deterministic advection-dispersion transport model and obtain corresponding concentration values of the model that are similar to the selected spatial locations of the aquifer. As past studies show, we approximate that the transverse dispersion coefficient is 10% of the longitudinal dispersion (Fetter, 1999). The mean velocity is 0.5 m/day. The profiles of both the aquifer and the deterministic model are plotted in one axis system, to compare their similarities. This trial-and-error curve fitting technique is carried out to determine the most accurate dispersion coefficients of the deterministic model. In this procedure concentration

By trial and error, we find that the closest fit is given by the longitudinal dispersion coefficient

Figure 3.14. Concentration profile of trial and error curve fit for D = 0.15 m2/day of the

of 0.15 m2/day, i.e., transverse dispersion is 0.015 m2/day, (Figure 3.14) for the aquifer.

methodology to approximate the aquifer parameters.

values of the aquifer are normalised (i.e., the values vary from 0 to 1).

advection dispersion model with row 5 of the aquifer data.

Figure 3.15. Concentration profiles of deterministic advection-dispersion model with D = 0.15 m2/day and SSTM with <sup>2</sup> = 0.01 and *b* = 0.01.

Subsequently, we develop a one-dimensional deterministic advection-dispersion model using the *D* obtained from a two-dimensional comparison. We then use a similar curve fitting technique with the 1-D deterministic model and 1-D stochastic model to find the most suitable <sup>2</sup> and *b* for the SSTM.

Figures 3.15 show the best fitting curves for both the SSTM and the 1-D advectiondispersion model. Having determined the appropriate parameters of the SSTM ( <sup>2</sup> = 0.01 and *b* = 0.01) that simulate the Lincoln University aquifer at the spatial location under consideration (Row 5 – well A) we investigated the robustness of the model for different Wiener processes and found that the SSTM is reasonably stable for the seven different Wiener realisations tested.

Even though the results show, as mentioned above, that the parameter combination of <sup>2</sup> = 0.01 and *b* = 0.01 is a fairly accurate representation of the experimental aquifer for the given spatial point, we need to test these parameters other spatial locations. Figure 3.16 shows that the 2-D deterministic advection-dispersion model with the longitudinal dispersion coefficient of 0.15 m2/day fits the aquifer data reasonably well for the similar locations.

Figure 3.16. Concentration profiles of 2-D deterministic advection-dispersion model (*D* = 0.15 m2/day) and the experimental aquifer.

much about the uniqueness of these set of parameters. A correlation length (*b*) of 0.01 m seems to be reasonable for a porous medium containing sand. But this has not been ascertained through the experiments in a laboratory. Therefore, one of the ways of determining the uniqueness of the parameters set is to develop the methodologies to estimate them using all the observations we have. We have used the maximum likelihood estimation (MLE) to determine the parameters of the deterministic advection and dispersion model for the aquifer data. However, MLE can not be used to estimate parameters of SSTM as the SSTM contains the stochastic integral for dispersion which can not be built into likelihood functions. Therefore, we make use of ANN to estimate the parameters of the

A Stochastic Model for Hydrodynamic Dispersion 89

t (days)

, (3.7.1)

Figure 3.19. Concentration profiles of the deterministic advection-dispersion model (*D* =

In chapter 1, we discuss the parameter estimation of the partial differential equations when the additive noise terms are present. The development of the likelihood functions are based on the theory of parameter estimation for stochastic process (Kutoyants, 1984; Lipster and Shirayer, 1977; Basawa and Prakasa Rao, 1980) and we use this theory to estimate

Let us consider a one-dimensional stochastic advection-dispersion equation with additive

where *DL* is the longitudinal hydrodynamic dispersion coefficient, in m2/day, and *<sup>x</sup> v* is

The two parameters to be estimated are *DL* and *<sup>x</sup> v* , and we can express the right hand

<sup>2</sup> ( ,) *<sup>L</sup> <sup>x</sup>*

*D v x t*

2

side of equation (3.7.1) as a linear function of the parameters (see section 1.4)

*C C C*

 

*t x x*

= 0.01 and *b* = 0.01 for row 7 well A.

SSTM

DET

SSTM for the experimental aquifer in section 3.6.

0.2

parameters.

fluctuations:

0.15 m2/day) and SSTM with <sup>2</sup>

the velocity of water flow, in m/day.

0.4

0.6

0.6

0.8

1 C

2 4 6 8 10 12 14

**3.7 Parameter Estimation using Maximum Likelihood Method** 

We compare the curves of the one-dimensional deterministic advection-dispersion model with D= 0.15 and of the SSTM with <sup>2</sup> = 0.01 and *b* = 0.01. Figure 3.17 illustrates that the curves for the deterministic model and the SSTM are in agreement.

Figure 3.17. Concentration profiles of the deterministic advection-dispersion model (*D* = 0.15 m2/day) and SSTM with <sup>2</sup> = 0.01 and *b* = 0.01 for row 3 well A.

Figures 3.18 and 3.19 show that the same parameter combination is valid for the data from well A of row 7.

Figure 3.18. Concentration profiles of the 2-D deterministic advection-dispersion model (*D* = 0.15 m2/day) and the experimental aquifer at row 7.

The comparative study of the SSTM and the experimental data described here is to show that the SSTM would faithfully mimic the concentration breakthrough curves in a real experiment situation. However, we do not attempt to validate the SSTM using the trial-anderror curve fitting techniques for all the data sets as if is too labour intensive. As the experimental aquifer consists of sand of similar particle size, one would expect that one set of <sup>2</sup> and b values should be able to characterise the whole aquifer. It appears that 2 =0.01 and b=0.01 would be an acceptable set of parameters. However, we can not say

= 0.01 and *b* = 0.01. Figure 3.17 illustrates that the

SSTM

DET

t (days)

88 in Porous Media - An Approach Based on Stochastic Calculus

We compare the curves of the one-dimensional deterministic advection-dispersion model

Figure 3.17. Concentration profiles of the deterministic advection-dispersion model (*D* =

Figures 3.18 and 3.19 show that the same parameter combination is valid for the data from

Figure 3.18. Concentration profiles of the 2-D deterministic advection-dispersion model (*D* =

The comparative study of the SSTM and the experimental data described here is to show that the SSTM would faithfully mimic the concentration breakthrough curves in a real experiment situation. However, we do not attempt to validate the SSTM using the trial-anderror curve fitting techniques for all the data sets as if is too labour intensive. As the experimental aquifer consists of sand of similar particle size, one would expect that one set

and b values should be able to characterise the whole aquifer. It appears that 2

=0.01 and b=0.01 would be an acceptable set of parameters. However, we can not say

= 0.01 and *b* = 0.01 for row 3 well A.

curves for the deterministic model and the SSTM are in agreement.

2 4 6 8 10

0.15 m2/day) and the experimental aquifer at row 7.

with D= 0.15 and of the SSTM with <sup>2</sup>

0.2

well A of row 7.

of <sup>2</sup> 

0.15 m2/day) and SSTM with <sup>2</sup>

0.4

0.6

0.8

1 C much about the uniqueness of these set of parameters. A correlation length (*b*) of 0.01 m seems to be reasonable for a porous medium containing sand. But this has not been ascertained through the experiments in a laboratory. Therefore, one of the ways of determining the uniqueness of the parameters set is to develop the methodologies to estimate them using all the observations we have. We have used the maximum likelihood estimation (MLE) to determine the parameters of the deterministic advection and dispersion model for the aquifer data. However, MLE can not be used to estimate parameters of SSTM as the SSTM contains the stochastic integral for dispersion which can not be built into likelihood functions. Therefore, we make use of ANN to estimate the parameters of the SSTM for the experimental aquifer in section 3.6.

Figure 3.19. Concentration profiles of the deterministic advection-dispersion model (*D* = 0.15 m2/day) and SSTM with <sup>2</sup> = 0.01 and *b* = 0.01 for row 7 well A.

### **3.7 Parameter Estimation using Maximum Likelihood Method**

In chapter 1, we discuss the parameter estimation of the partial differential equations when the additive noise terms are present. The development of the likelihood functions are based on the theory of parameter estimation for stochastic process (Kutoyants, 1984; Lipster and Shirayer, 1977; Basawa and Prakasa Rao, 1980) and we use this theory to estimate parameters.

Let us consider a one-dimensional stochastic advection-dispersion equation with additive fluctuations:

$$\frac{\partial \mathbf{C}}{\partial t} = D\_{\perp} \left( \frac{\partial^{2} \mathbf{C}}{\partial \mathbf{x}^{2}} \right) - \left. \upsilon\_{\mathbf{x}} \right| \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right) + \left. \xi \langle \mathbf{x}, t \rangle \right|,\tag{3.7.1}$$

where *DL* is the longitudinal hydrodynamic dispersion coefficient, in m2/day, and *<sup>x</sup> v* is the velocity of water flow, in m/day.

The two parameters to be estimated are *DL* and *<sup>x</sup> v* , and we can express the right hand side of equation (3.7.1) as a linear function of the parameters (see section 1.4)

*C* = solute concentration, mg/l,

*ne* = effective porosity.

*dh*

Depth (m)

aquifer could be neglected.

*K* = hydraulic conductivity, m/day,

*dl* = hydraulic gradient, m/m, and

longitudinal hydrodynamic dispersion, *DL* (m2/day).

We have used the experimental data from the artificial aquifer (Figure 3.13) to estimate *DL* and *<sup>x</sup> v* . Table 3.1 shows the estimated and experimental values for *DL* and *K* (calculated

A Stochastic Model for Hydrodynamic Dispersion 91

0.4 203.2 137 0.167 0.1596 1.0 210.6 137 0.143 0.1596 1.6 208.9 137 0.134 0.1596 2.2 262.3 137 0.242 0.1596

Table 3.1. Estimated and experimental parameters hydraulic conductivity, *K* (m/day) and

The results show the accuracy of *DL* estimates is better than that of *K*. We estimate the *vx* by using equation (3.7.6) for and for simplicity, we assume that the hydraulic gradient, *dh*

parameters to less

and the effective porosity, *ne* to be constants, or in other words, their spatial distributions are homogeneous. However, distribution of these values in the aquifer may be slightly heterogeneous. This assumption may have influenced the accuracy of less precise *K* estimates. However, the *DL* values are estimated directly and not affected by such

Another possible phenomenon that can be present in solute transport is adsorption and the occurrence of short circuits which are, assumed to be included in the random component,

( ,) *x t* , of the governing equation, equation (3.7.1). However, we assumed that in the experiments the tracer is mixed in the upstream header tank, hence, adsorption in the

As mentioned earlier, in the process of calibration to obtain parameters, it is assumed that the experimental aquifer is homogeneous. However, Figures 3.20, 3.21 and 3.22 show that the aquifer did not behave as expected (please refer to Figure 3.13 for notation). Figure 3.10 shows that the concentration values at a well which is closer to the middle of the artificial aquifer (Row 5 – Well B). This well is approximately 3 m away from the header tank as shown in Figure 3.13. Other Figures, 3.21 and 3.22, show that the concentration profiles of

assumptions, and the maximum likelihood estimates are closer to the experimental.

Hydraulic conductivity, *K* (m/day) Longitudinal hydrodynamic

Estimated Experimental Estimated Experimental

dispersion, *DL* (m2/day)

*dl*

from equation (3.7.6)) at different depths of the aquifer from the top surface.

$$f(t, \mathbb{C}, \theta) = a\_0(\mathbb{C}, \ t) + \theta\_1 a\_1(\mathbb{C}, t) + \theta\_2 a\_2(\mathbb{C}, \ t). \tag{3.7.2}$$

where

$$a\_0 \mathbf{C}\_{\iota'} \mathbf{t} = \mathbf{0}; \quad a\_1 \mathbf{C}\_{\iota'} \mathbf{t} \mathbf{t} = \left(\frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2}\right)\_{\iota}; \quad a\_2 \mathbf{C}\_{\iota'} \mathbf{t} \mathbf{t} = -\left(\frac{\partial \mathbf{C}}{\partial \mathbf{x}}\right)\_{\iota}; \; \boldsymbol{\theta}\_1 = \boldsymbol{D}\_{\iota}; \quad \boldsymbol{\theta}\_2 = \boldsymbol{v}\_{\iota} \dots$$

Then the log-likelihood function can be written as,

$$\begin{split} l(\boldsymbol{\theta}\_{1}, \boldsymbol{\theta}\_{2}) &= \sum\_{i=1}^{M} \Big[ \left\{ a\_{0}(\mathbf{C}\_{i}, t) + \theta\_{1} a\_{1}(\mathbf{C}\_{i}, t) + \theta\_{2} a\_{2}(\mathbf{C}\_{i}, t) \right\} d\mathbf{C}\_{i}(t) \\ &- \frac{1}{2} \sum\_{i=1}^{M} \Big[ \left\{ a\_{0}(\mathbf{C}\_{i}, t) + \theta\_{1} a\_{1}(\mathbf{C}\_{i}, t) + \theta\_{2} a\_{2}(\mathbf{C}\_{i}, t) \right\}^{2} dt \end{split} \tag{3.7.3}$$

Differentiating equation (3.7.3) with respect to 1 and <sup>2</sup> , respectively, we obtain the following two simultaneous equations:

$$\sum\_{i=1}^{M} \sum\_{0}^{T} \left\{ a\_{1}(\mathbf{C}\_{i}, t) \right\} d\mathbf{C}\_{i}(t) \cdot \sum\_{i=1}^{M} \left\{ \left\{ a\_{0}(\mathbf{C}\_{i}, t) + \theta\_{1} a\_{1}(\mathbf{C}\_{i}, t) + \theta\_{2} a\_{2}(\mathbf{C}\_{i}, t) \right\} \left\{ a\_{1}(\mathbf{C}\_{i}, t) \right\} dt = 0 \right\}.$$

$$\sum\_{i=1}^{M} \left[ \left\{ a\_{2}(\mathbf{C}\_{i}, t) \right\} d\mathbf{C}\_{i}(t) \cdot \sum\_{i=1}^{M} \left\{ a\_{0}(\mathbf{C}\_{i}, t) + \theta\_{1} a\_{1}(\mathbf{C}\_{i}, t) + \theta\_{2} a\_{2}(\mathbf{C}\_{i}, t) \right\} \left\{ a\_{2}(\mathbf{C}\_{i}, t) \right\} dt = 0 \right. \tag{3.7.4}$$

Suppose we have the observations of solute concentration, *Ci* at *M* independent space coordinates along the *x*-axis, where 1  *i M,* at different time intervals, *t* (where 0  *t T* , and *T* is an integer that represents the last reading taken at unit intervals on *t*-axis). In other words, we have *M* number of *Ci* observations for each time step. Hence, there are, altogether, ((*T+1*)*M*) *Ci* observations. We use these observations to estimate the parameter 1 and <sup>2</sup> ,

We substitute 0 ( ,) *<sup>i</sup> aCt* , 1( ,) *<sup>i</sup> aCt* , 2 ( ,) *<sup>i</sup> aCt* , 1 and 2 in equations (3.7.4) to obtain the following set of equations,

$$\sum\_{i=1}^{M} \sum\_{t=0}^{T} \left\{ \frac{\partial^2 \mathbf{C}\_i}{\partial \mathbf{x}^2} \right\} \Delta \mathbf{C}\_i - D\_L \sum\_{i=1}^{M} \sum\_{t=0}^{T} \left\{ \frac{\partial^2 \mathbf{C}\_i}{\partial \mathbf{x}^2} \right\}^2 \Delta t - \left. \upsilon\_x \sum\_{i=1}^{M} \sum\_{t=0}^{T} \left\{ \frac{\partial \mathbf{C}\_i}{\partial \mathbf{x}} \right\} \right\} \left\{ \frac{\partial^2 \mathbf{C}\_i}{\partial \mathbf{x}^2} \right\} \Delta t = 0$$
 
$$\sum\_{i=1}^{M} \sum\_{t=0}^{T} \left\{ \frac{\partial \mathbf{C}\_i}{\partial \mathbf{x}} \right\} \Delta \mathbf{C}\_i - D\_L \sum\_{i=1}^{M} \sum\_{t=0}^{T} \left\{ \frac{\partial \mathbf{C}\_i}{\partial \mathbf{x}} \right\} \left\{ \frac{\partial^2 \mathbf{C}\_i}{\partial \mathbf{x}^2} \right\} \Delta t - \left. \upsilon\_x \sum\_{i=1}^{M} \sum\_{t=0}^{T} \left\{ \frac{\partial \mathbf{C}\_i}{\partial \mathbf{x}} \right\}^2 \Delta t = 0. \tag{3.7.5}$$

Therefore, *DL* and *<sup>x</sup> v* values can be obtained by solving these two simultaneous equations. Once *<sup>x</sup> v* is known, the hydraulic conductivity (*K*) can be obtained from,

$$w\_x = \frac{K}{n\_\epsilon} \left(\frac{dl}{dl}\right) \tag{3.7.6} = \text{average linear velocity, m/day} \tag{3.7.6}$$


*i*

*a C t dt* <sup>2</sup> ( ,) 0 *<sup>i</sup>* . (3.7.4)

 <sup>1</sup> ; 

2

*x*  (3.7.2)

(3.7.3)

<sup>2</sup> , respectively, we obtain the

2 in equations (3.7.4) to obtain the

 *DL* 2 *<sup>x</sup> v* .

 

90 in Porous Media - An Approach Based on Stochastic Calculus

<sup>0</sup> 1 1 2 2 *f* ( , , ) ( , ) ( , ) ( , ). *tC a C t a Ct a C t*

<sup>2</sup> ( ,) ; *<sup>i</sup>*

( , ) ( , ) ( , ) ( , ) (t)

<sup>1</sup> <sup>0</sup> 1 1 2 2 <sup>1</sup>

*i i i i i i*

*a C t dC t a C t a C t a C t a C t dt* 

Suppose we have the observations of solute concentration, *Ci* at *M* independent space coordinates along the *x*-axis, where 1  *i M,* at different time intervals, *t* (where 0  *t T* , and *T* is an integer that represents the last reading taken at unit intervals on *t*-axis). In other words, we have *M* number of *Ci* observations for each time step. Hence, there are, altogether, ((*T+1*)*M*) *Ci* observations. We use these observations to estimate the parameter

> 1 and

<sup>2</sup> <sup>2</sup> <sup>2</sup> <sup>2</sup> 2 2 2

*i L x*

*x x x x*

1 0 1 0 1 0

*x x x x*

*i t i t i t*

Once *<sup>x</sup> v* is known, the hydraulic conductivity (*K*) can be obtained from,

*M T M T M T*

*i i i i*

*i i i i i L x*

Therefore, *DL* and *<sup>x</sup> v* values can be obtained by solving these two simultaneous equations.

*<sup>C</sup> C C <sup>C</sup> C D t v <sup>t</sup>*

*<sup>C</sup> <sup>C</sup> C C C D t v <sup>t</sup>*

2

.

*l a C t a C t a C t dC*

( , ) ( )- ( , ) ( , ) ( , ) ( , ) 0

{ ( ,) ( ,) ( ,)

*aCt aCt aCt* 

*i i i*

*i i i*

*a C t a C t a C t dt*

1 and

<sup>1</sup> ( ,) ( ,) ( ,) 2

0 1 1 2 2

*<sup>C</sup> aCt*

*i i i i*

 

> 

> >

 



(3.7.5)

2 2

= average linear velocity, m/day, (3.7.6)

 

*<sup>C</sup> aCt*

1 0

*i M T*

*M T*

1 0

*M T*

*i*

*i*

Then the log-likelihood function can be written as,

Differentiating equation (3.7.3) with respect to

1 0 1 0

*i i*

2 1 0

*M T*

*i*

<sup>2</sup> ,

following set of equations,

1 and { ( , )} ( )-

*a C t dC t*

*i i*

We substitute 0 ( ,) *<sup>i</sup> aCt* , 1( ,) *<sup>i</sup> aCt* , 2 ( ,) *<sup>i</sup> aCt* ,

*x*

*v*

*e K dh*

*n dl* 

*M T M T*

 

following two simultaneous equations:

2 <sup>1</sup> <sup>2</sup> ( ,) ; *<sup>i</sup>*

*x* 

*i*

1 2 0 1 1 2 2

<sup>0</sup> 1 1 2 2 <sup>1</sup> <sup>0</sup>

1 0 1 0 1 0

*i t i t i t*

*M T M T M T*

where

<sup>0</sup> ( , ) 0; *<sup>i</sup> aCt*

We have used the experimental data from the artificial aquifer (Figure 3.13) to estimate *DL* and *<sup>x</sup> v* . Table 3.1 shows the estimated and experimental values for *DL* and *K* (calculated from equation (3.7.6)) at different depths of the aquifer from the top surface.


Table 3.1. Estimated and experimental parameters hydraulic conductivity, *K* (m/day) and longitudinal hydrodynamic dispersion, *DL* (m2/day).

The results show the accuracy of *DL* estimates is better than that of *K*. We estimate the *vx* by using equation (3.7.6) for and for simplicity, we assume that the hydraulic gradient, *dh dl* and the effective porosity, *ne* to be constants, or in other words, their spatial distributions are homogeneous. However, distribution of these values in the aquifer may be slightly heterogeneous. This assumption may have influenced the accuracy of less precise *K* estimates. However, the *DL* values are estimated directly and not affected by such assumptions, and the maximum likelihood estimates are closer to the experimental.

Another possible phenomenon that can be present in solute transport is adsorption and the occurrence of short circuits which are, assumed to be included in the random component, ( ,) *x t* , of the governing equation, equation (3.7.1). However, we assumed that in the experiments the tracer is mixed in the upstream header tank, hence, adsorption in the aquifer could be neglected.

As mentioned earlier, in the process of calibration to obtain parameters, it is assumed that the experimental aquifer is homogeneous. However, Figures 3.20, 3.21 and 3.22 show that the aquifer did not behave as expected (please refer to Figure 3.13 for notation). Figure 3.10 shows that the concentration values at a well which is closer to the middle of the artificial aquifer (Row 5 – Well B). This well is approximately 3 m away from the header tank as shown in Figure 3.13. Other Figures, 3.21 and 3.22, show that the concentration profiles of

The reason that the artificial aquifer does not behave homogeneously may be due to the method of construction. The aquifer was constructed using sand blocks that were laid layer by layer. We assume that even though material used in the aquifer is uniform, joints in the blocks can create diverse flow patterns and different flow lengths. Besides, due to the high pressure on the bottom layers (from the top layers), they may be more compacted and,

A Stochastic Model for Hydrodynamic Dispersion 93

We can extend the parameter estimation procedure to determine parameters of a twodimensional groundwater problem. In the two-dimensional case, an advancing solute front will also tend to spread in the directions normal to the direction of flow because at the pore scale the flow paths can diverge. This results in mixing in the directions normal to the flow path, which is called transverse dispersion. Considering the transverse dispersion, the two-

dimensional advection-dispersion equation can be written as (Fetter, 1999),

2 2 *L* 2 *T* 2 *x C C <sup>C</sup> C D D dt v dt x y x*

 *DL =* hydrodynamic dispersion coefficient parallel to the principal direction of

 *DT =* hydrodynamic dispersion coefficient perpendicular to the principal direction

, (3.7.7)

therefore, behave differently.

Figure 3.22. Concentrations profiles at Row 7 – Well A.

where *C* = solute concentration (mg/l),

flow (longitudinal) (m2/day),

of flow (transverse) (m2/day), and

*vx* = average linear velocity (m/day).

t = time (day),

wells at Well B at Row 3 and Well A at Row 7, respectively. These figures demonstrate that the concentration values are not the same at all the depths and, hence, the behaviour of the aquifer is not similar throughout. The plots of other wells also exhibit heterogeneous behaviour. Therefore, we can state that the aquifer is not behaving homogeneously, meaning that the aquifer parameters, such as hydraulic gradient and effective porosity are not uniformly distributed throughout the system. The variables used to calibrate the aquifer parameters are subjected to randomness and the accuracy of the results could be affected, considerably.

Figure 3.20. Concentration profiles at Row 5 – Well B.

Figure 3.21. Concentrations profiles at Row 3 – Well B.

92 in Porous Media - An Approach Based on Stochastic Calculus

wells at Well B at Row 3 and Well A at Row 7, respectively. These figures demonstrate that the concentration values are not the same at all the depths and, hence, the behaviour of the aquifer is not similar throughout. The plots of other wells also exhibit heterogeneous behaviour. Therefore, we can state that the aquifer is not behaving homogeneously, meaning that the aquifer parameters, such as hydraulic gradient and effective porosity are not uniformly distributed throughout the system. The variables used to calibrate the aquifer parameters are subjected to randomness and the accuracy of the results could be affected,

considerably.

Figure 3.20. Concentration profiles at Row 5 – Well B.

Figure 3.21. Concentrations profiles at Row 3 – Well B.

The reason that the artificial aquifer does not behave homogeneously may be due to the method of construction. The aquifer was constructed using sand blocks that were laid layer by layer. We assume that even though material used in the aquifer is uniform, joints in the blocks can create diverse flow patterns and different flow lengths. Besides, due to the high pressure on the bottom layers (from the top layers), they may be more compacted and, therefore, behave differently.

Figure 3.22. Concentrations profiles at Row 7 – Well A.

We can extend the parameter estimation procedure to determine parameters of a twodimensional groundwater problem. In the two-dimensional case, an advancing solute front will also tend to spread in the directions normal to the direction of flow because at the pore scale the flow paths can diverge. This results in mixing in the directions normal to the flow path, which is called transverse dispersion. Considering the transverse dispersion, the twodimensional advection-dispersion equation can be written as (Fetter, 1999),

$$\partial \mathbf{C} = \left\{ D\_{\perp} \left( \frac{\partial^{2} \mathbf{C}}{\partial \mathbf{x}^{2}} \right) + D\_{\mathbf{r}} \left( \frac{\partial^{2} \mathbf{C}}{\partial y^{2}} \right) \right\} dt - \upsilon\_{\mathbf{x}} \left( \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right) dt \,, \tag{3.7.7}$$

where *C* = solute concentration (mg/l),

t = time (day),

	- *DT =* hydrodynamic dispersion coefficient perpendicular to the principal direction

of flow (transverse) (m2/day), and

*vx* = average linear velocity (m/day).

We use the observations for *C*(*x, y, t*) at *M* discrete points in (*x, y*) coordinate space for a period of time *t* (where 0 ) *t T* . Then we obtain the estimates for two unknown

A Stochastic Model for Hydrodynamic Dispersion 95

2 2 2 2 2 2 2 2

*C C K dh C C C <sup>D</sup> dt x y n dl x x y*

*<sup>C</sup> C C K dh C C dC D dt*

1 1 1

<sup>2</sup> <sup>1</sup> <sup>2</sup> <sup>2</sup> 1 2 2 1

*C C <sup>k</sup> t t x y*

0.1 *M n j j*

1 2 2 1

0.1 *M n j j*

2 2 2 1

, and

<sup>1</sup> <sup>1</sup>

.

*C CC <sup>k</sup> t t xxy*

*M n j j j i i i*

,

,

*j j i i i*

*j j*

*C C m C C x y* 

*CC C <sup>l</sup> t t xx y*

 ,

*M n j j j i i i*

*i i*

,

*K dh Dk l m n dl K dh Dk l m n dl*

2 2 2


. (3.7.12)

0,

0,

1

*i i*

<sup>1</sup> <sup>2</sup> <sup>2</sup>

0.1

*j j i i*

<sup>1</sup> 2 2

<sup>1</sup> <sup>2</sup> <sup>2</sup>

*j j*

*j j*

*j j*

,

(3.7.13)


2 2 2 2

*x x y n dl x x*

*e*

*e*

1 2 2

<sup>2</sup> <sup>1</sup> 2 1 1 0

*C m C C x* 

*<sup>C</sup> <sup>l</sup> t t x*

parameters as the solution to the following simultaneous equations:

*i*

We can simplify equation (3.7.11) to

*L*

1 0 1 0

We can rewrite equations (3.7.12) as

*M T M T*

1 0

*i M T*

*M T*

1 0

where

2 2 2 2

0.1

*i e*

*i L*

*i i e*

*L*

*L*

1 0

1 0

1 0

1 0

*M n j i*

*M n j*

*i j*

*i j*

2

1 0

*i j*

*i j*

*i j*

*i j*

*C C dC x y*

The randomness of heterogeneous groundwater systems can be accounted for by adding a noise component to equation (3.7.7), and it can be given by

$$
\partial \mathbf{C} = \left\{ D\_L \left( \frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2} \right) + D\_T \left( \frac{\partial^2 \mathbf{C}}{\partial y^2} \right) \right\} dt - v\_x \left( \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right) dt + \xi(\mathbf{x}, t) \, dt \,, \tag{3.7.8}
$$

where ( ,) *x t* is assumed to be a zero-mean stochastic process.

We multiply equation (3.6.8) by *dt* throughout and replace ( ,) *x t dt* by <sup>2</sup> *dB t*( ) (Jazwinski, 1970), where <sup>2</sup> is the amplitude of the Wiener increments, *dB t*( ) , to obtain equation (3.7.9). We can now obtain the stochastic partial differential equation as follows,

$$\mathcal{C}\mathbf{C} = \left| D\_{\perp} \left( \frac{\partial^{2} \mathbf{C}}{\partial \mathbf{x}^{2}} \right) + D\_{\Gamma} \left( \frac{\partial^{2} \mathbf{C}}{\partial \mathbf{y}^{2}} \right) \right| dt - v\_{\mathbf{x}} \left( \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right) dt + \sigma^{2} dB(\mathbf{f}) \,. \tag{3.7.9}$$

As we described in equation (3.7.6) average linear velocity can be expressed by

$$v\_x = \frac{K}{n\_e} \left(\frac{dh}{dl}\right)\_{\prime}$$

Hence, equation (3.7.9) becomes,

$$\partial \mathbf{C} = \left\{ D\_{\perp} \left( \frac{\partial^{2} \mathbf{C}}{\partial \mathbf{x}^{2}} \right) + D\_{r} \left( \frac{\partial^{2} \mathbf{C}}{\partial y^{2}} \right) \right\} dt - \frac{\mathbf{K}}{n\_{e}} \left( \frac{dh}{dl} \right) \left( \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right) dt + \sigma^{2} d\mathbf{B}(t) \,. \tag{3.7.10}$$

We assume that transverse dispersion ( *DT* ) can be approximated to 10% of the longitudinal dispersion *DL* , i.e., *DT* = 0.1 *DL* (Felter,1999).

Then equation (3.7.10) becomes,

2 2 2 <sup>2</sup> <sup>2</sup> 0.1 ( ) *<sup>L</sup> e C C K dh C C D dt dt dB t x y n dl x* . (3.7.11) 

Equation (3.7.11) can be written in the following form:

$$f(t, \mathbb{C}, \theta) \,=\, a\_0(\mathbb{C}, t) + \theta\_1 a\_1(\mathbb{C}, t) + \theta\_2 a\_2(\mathbb{C}, t),$$

where,

$$a\_0 \begin{pmatrix} \mathbf{C}\_{\vee} & t \end{pmatrix} = \mathbf{0}; \qquad a\_1 \mathbf{C}\_{\vee} t \mathbf{} = \left\{ \frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2} + \mathbf{0}. \mathbf{1} \left( \frac{\partial^2 \mathbf{C}}{\partial y^2} \right) \right\}; \quad a\_2 \mathbf{(C}\_{\vee} t \mathbf{)} = -\left( \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right); \quad \mathbf{0}$$

$$\boldsymbol{\theta}\_1 = \boldsymbol{D}\_{\perp}; \qquad \boldsymbol{\theta}\_2 = \boldsymbol{\upsilon}\_x = \frac{\mathbf{K}}{n\_\epsilon} \left( \frac{d\mathbf{t}}{dl} \right).$$

We use the observations for *C*(*x, y, t*) at *M* discrete points in (*x, y*) coordinate space for a period of time *t* (where 0 ) *t T* . Then we obtain the estimates for two unknown parameters as the solution to the following simultaneous equations:

We can simplify equation (3.7.11) to

$$\begin{split} &\sum\_{i=1}^{M} \Biggl[ \left| \frac{\partial^{2} \mathbb{C}}{\partial \mathbf{x}^{2}} + \mathbf{0}.1 \left( \frac{\partial^{2} \mathbb{C}}{\partial \mathbf{y}^{2}} \right) \right| \mathbf{d} \mathbf{C}\_{i} \\ &\sum\_{i=1}^{M} \prod\_{j=1}^{T} \Bigg[ D\_{L} \left\{ \frac{\partial^{2} \mathbb{C}}{\partial \mathbf{x}^{2}} + \mathbf{0}.1 \left( \frac{\partial^{2} \mathbb{C}}{\partial \mathbf{y}^{2}} \right) \right\} + \frac{K}{n\_{\epsilon}} \left( \frac{d \mathbb{h}}{d \mathbb{I}} \right) \left( \frac{\partial \mathbb{C}}{\partial \mathbf{x}} \right) \Bigg] \Bigg\{ \frac{\partial^{2} \mathbb{C}}{\partial \mathbf{x}^{2}} + \mathbf{0}.1 \left( \frac{\partial^{2} \mathbb{C}}{\partial \mathbf{y}^{2}} \right) \Bigg] d\mathbf{t} = \mathbf{0} \end{split}$$

$$\sum\_{i=1}^{M} \Bigg[ \left\{ \frac{\partial \mathbb{C}}{\partial \mathbf{x}} \right\} d \mathbf{C}\_{i} - \sum\_{i=1}^{M} \int\_{0}^{t} \Bigg\{ D\_{L} \left\{ \frac{\partial^{2} \mathbb{C}}{\partial \mathbf{x}^{2}} + \mathbf{0}.1 \left( \frac{\partial^{2} \mathbb{C}}{\partial \mathbf{y}^{2}} \right) \Bigg\} + \frac{K}{n\_{\epsilon}} \left( \frac{d \mathbb{h}}{d \mathbb{I}} \right) \Bigg\} \Bigg\} \Bigg\{ \frac{\partial \mathbb{C}}{\partial \mathbf{x}} \Bigg\} d \mathbf{t} = \mathbf{0} \end{split} \tag{3.7.12}$$

We can rewrite equations (3.7.12) as

$$\begin{aligned} D\_L k\_1 + \frac{K}{n\_e} \left( \frac{d\mathbf{l}}{dl} \right) l\_1 + m\_1 &= \mathbf{0}, \\ D\_L k\_2 + \frac{K}{n\_e} \left( \frac{d\mathbf{l}}{dl} \right) l\_2 + m\_2 &= \mathbf{0}, \end{aligned} \tag{3.7.13}$$

where

Computational Modelling of Multi-Scale Non-Fickian Dispersion

is the amplitude of the Wiener increments, *dB t*( ) , to obtain

2

2

2

; 2 ( ,) ; *<sup>i</sup>*

*<sup>C</sup> aCt*

*i*

*x* 

, (3.7.8)

*dB t*( )

( ,) *x t dt* by <sup>2</sup>

. (3.7.9)

. (3.7.10)

. (3.7.11)

94 in Porous Media - An Approach Based on Stochastic Calculus

The randomness of heterogeneous groundwater systems can be accounted for by adding a

<sup>2</sup> <sup>2</sup> ( ,) *<sup>L</sup> <sup>T</sup> <sup>x</sup> C C <sup>C</sup> C D D dt v dt x t dt x y x*

equation (3.7.9). We can now obtain the stochastic partial differential equation as follows,

<sup>2</sup> <sup>2</sup> ( ) *<sup>L</sup> <sup>T</sup> <sup>x</sup> C C <sup>C</sup> C D D dt v dt dB t x y x*

noise component to equation (3.7.7), and it can be given by

where

(Jazwinski, 1970), where <sup>2</sup>

Hence, equation (3.7.9) becomes,

Then equation (3.7.10) becomes,

where,

2 2

( ,) *x t* is assumed to be a zero-mean stochastic process.

2 2

2 2

2 2

<sup>1</sup> ; 

dispersion *DL* , i.e., *DT* = 0.1 *DL* (Felter,1999).

Equation (3.7.11) can be written in the following form:

<sup>0</sup> ( , ) 0; *<sup>i</sup> aC t*

As we described in equation (3.7.6) average linear velocity can be expressed by

*x*

*v*

*e K dh*

*n dl* ,

<sup>2</sup> <sup>2</sup> ( ) *<sup>L</sup> <sup>T</sup> e*

We assume that transverse dispersion ( *DT* ) can be approximated to 10% of the longitudinal

<sup>2</sup> <sup>2</sup> 0.1 ( ) *<sup>L</sup> e*

<sup>0</sup> 1 1 2 2 *f* ( , , ) ( , ) ( , ) ( , ), *tC a Ct a Ct a Ct*

2 2

*x y* 

*v*

 

*e K dh*

*n dl*

.

 

<sup>1</sup> <sup>2</sup> <sup>2</sup> ( , ) 0.1 *C C a Ct*

*DL* 2 *<sup>x</sup>*

*C C K dh C C D dt dt dB t x y n dl x*

*C C K dh C C D D dt dt dB t x y n dl x*

We multiply equation (3.6.8) by *dt* throughout and replace

 <sup>2</sup> <sup>1</sup> <sup>2</sup> <sup>2</sup> 1 2 2 1 1 0 0.1 *M n j j i i j j i j C C <sup>k</sup> t t x y* , <sup>1</sup> <sup>2</sup> <sup>2</sup> 1 2 2 1 1 0 0.1 *M n j j j i i i j j i j CC C <sup>l</sup> t t xx y* , <sup>1</sup> <sup>2</sup> <sup>2</sup> 1 1 2 2 1 0 0.1 *M n j j j j i i i i i j C C m C C x y* , <sup>1</sup> 2 2 2 2 2 1 1 0 *M n j j j i i i j j i j C CC <sup>k</sup> t t xxy* , <sup>2</sup> <sup>1</sup> 2 1 1 0 *M n j i j j i j <sup>C</sup> <sup>l</sup> t t x* , and <sup>1</sup> <sup>1</sup> 2 1 0 *M n j j j i i i i j C m C C x* . 1 0

The deterministic solute concentration values are generated for a 10 m x 5 m 2-D aquifer using equation (3.6.7). Eight hundred data examples (patterns) are generated for different hydraulic conductivity, *K,* values that ranged from 40 to 240 m/day. It is assumed that all other parameters, control variables and subsidiary conditions are fixed. An initial concentration value of 100 ppm is considered as a point source at the middle of the header boundary of the aquifer and the same source is maintained at the boundary throughout the 10 day period considered. Exponentially distributed concentration values of the point source (at the middle of the header boundary) are considered along with the longitudinal and lateral directions as the initial conditions for the other spatial coordinates. We gather 50 input values for each example. These input values represent solute concentration values at 10 spatial locations (Figure 6.1) at five different time intervals; *t* = 1, *t* = 3, *t* = 5, *t* = 7, *t* = 10 day. We examine the possibility of amalgamating the time as an independent variable into the concentration input data. However, it is difficult to meaningfully integrate them into presently available ANN architectures and innovative model structures need to be

A Stochastic Model for Hydrodynamic Dispersion 97

Figure 3.23. Spatial coordinates of concentration observations of the 2-D aquifer.

It is noted, in many studies, that determination of an appropriate network architecture is one of the most important but also one of the most difficult tasks in the model building process (Maier and Dandy, 2000). The network architecture determines the number of connection weights (free parameters) and the way information flows through the network. It is common practice to fix the number of hidden layers in the network and then to choose the number of nodes in each of these layers. Initially, it may be worthwhile to consider a network with simple properties. Such networks have higher processing speeds and can be implemented on hardware more economically (Towell et al., 1991; Bebis and Georgiopoulos,

developed.

We can develop the finite difference approximations for *<sup>i</sup> i C x* , 2 2 *i i C x* and 2 2 *i i C y* , to obtain values of *k1, l1, m1*, *k2, l2,* and *m2*. Furthermore, by substituting average values for the hydraulic gradient, *dh dl* , and the effective porosity, *ne* , two simultaneous equations (3.7.13) can be solved to obtain *DL* and *K,* for a two-dimensional groundwater system.


### **3.8 Parameter Estimation using ANN**

The two-dimensional deterministic advection–dispersion equation (Fetter, 1999), is used as the governing equation for this section. It is important to mention that other possible phenomenon that can be present in the solute transport, such as adsorption and the occurrence of short circuits, are neglected in the governing equation on the assumption that the introduction of noise into the solute concentration values used to estimate the parameters would compensate for them.

2 2 *i i C x* 

Longitudinal hydrodynamic dispersion,

and

2 2 *i i C y* 

, to obtain

*i C x* ,

, and the effective porosity, *ne* , two simultaneous equations

*DL* (m2/day)

96 in Porous Media - An Approach Based on Stochastic Calculus

values of *k1, l1, m1*, *k2, l2,* and *m2*. Furthermore, by substituting average values for the

We estimate the parameters of the artificial aquifer at Lincoln University using the procedure developed. The observations of solute concentration are obtained at 1 m grid intervals at four different levels (0.4, 1.0, 1.6, and 2.2 m from the surface of the aquifer).

Table 3.2 shows the estimated parameters for the two-dimensional aquifer and the estimates have come closer to the experimental values. However, as the experimental values may not represent the real values within the aquifer, further computational experiments by changing

Estimated Experimental Estimated Experimental

Table 3.2. Estimated and experimental parameters, hydraulic conductivity, *K* (m/day), and

In estimating parameters with ANN, first we employ a deterministic 2-D advectiondispersion transport numerical model to generate synthetic data. Afterwards, ANN are trained to learn the complex excitation and response relationship of the generated data. This is done by training the network sufficiently to minimise the error between the actual and network response while retaining the generalising capabilities of the network. We then estimate the associated parameters using noisy concentration data that represent real world aquifer systems. We also test the ability of the model to estimate hydraulic conductivity of

The two-dimensional deterministic advection–dispersion equation (Fetter, 1999), is used as the governing equation for this section. It is important to mention that other possible phenomenon that can be present in the solute transport, such as adsorption and the occurrence of short circuits, are neglected in the governing equation on the assumption that the introduction of noise into the solute concentration values used to estimate the

(3.7.13) can be solved to obtain *DL* and *K,* for a two-dimensional groundwater system.

Hence, we rearrange the dataset to be four two-dimensional datasets for each level.

0.4 198.2 137 0.165 0.1596 1.0 166.7 137 0.162 0.1596 1.6 171.5 137 0.143 0.1596 2.2 231.1 137 0.197 0.1596

longitudinal hydrodynamic dispersion, *DL* (m2/day), for the aquifer.

We can develop the finite difference approximations for *<sup>i</sup>*

*dl* 

values of Experimental

the lateral dispersion would not improve the estimates.

Hydraulic conductivity, *K*

(m/day)

**3.8 Parameter Estimation using ANN** 

an artificial experimental aquifer.

parameters would compensate for them.

hydraulic gradient, *dh*

Depth (m)

The deterministic solute concentration values are generated for a 10 m x 5 m 2-D aquifer using equation (3.6.7). Eight hundred data examples (patterns) are generated for different hydraulic conductivity, *K,* values that ranged from 40 to 240 m/day. It is assumed that all other parameters, control variables and subsidiary conditions are fixed. An initial concentration value of 100 ppm is considered as a point source at the middle of the header boundary of the aquifer and the same source is maintained at the boundary throughout the 10 day period considered. Exponentially distributed concentration values of the point source (at the middle of the header boundary) are considered along with the longitudinal and lateral directions as the initial conditions for the other spatial coordinates. We gather 50 input values for each example. These input values represent solute concentration values at 10 spatial locations (Figure 6.1) at five different time intervals; *t* = 1, *t* = 3, *t* = 5, *t* = 7, *t* = 10 day. We examine the possibility of amalgamating the time as an independent variable into the concentration input data. However, it is difficult to meaningfully integrate them into presently available ANN architectures and innovative model structures need to be developed.

Figure 3.23. Spatial coordinates of concentration observations of the 2-D aquifer.

It is noted, in many studies, that determination of an appropriate network architecture is one of the most important but also one of the most difficult tasks in the model building process (Maier and Dandy, 2000). The network architecture determines the number of connection weights (free parameters) and the way information flows through the network. It is common practice to fix the number of hidden layers in the network and then to choose the number of nodes in each of these layers. Initially, it may be worthwhile to consider a network with simple properties. Such networks have higher processing speeds and can be implemented on hardware more economically (Towell et al., 1991; Bebis and Georgiopoulos,

smaller permissible parameter regimes of *K*; (i) 40-90, (ii) 90-140, (iii) 140-190, and (iv) 190-240

A Stochastic Model for Hydrodynamic Dispersion 99

Figure 3.24. Absolute error of estimated parameter, *K,* when considering the whole range,

The maximum error in the 190-240 range has been reduced by about 90% (Figure 3.25 and Table 3.3). Therefore, it is reasonable to assume that if we can gather prior information about the system under consideration, it is possible to obtain more accurate estimates. However, in the real world problems the prior knowledge of the system is limited. Later, we discuss a method to identify the range of parameters by using Self-Organising Maps (SOM). However, before using SOM, we explore the robustness of the ANN estimation models.

> Max =2.98

Figure 3.25. Absolute Error of estimated parameter *K,* only for the range 190 – 240 m/day.

40–240 m/day.

m/day. Table 3.3 shows that accuracy of the estimates has improved considerably.

1994; Castellano et al., 1997). Cybenko (1989) showed that only one hidden layer was required to approximate any continuous function, given that sufficient degrees of freedom (i.e. connection weights) are provided. However, in practice many functions are difficult to approximate with one hidden layer, requiring a prohibitive number of hidden layer nodes (Flood and Kartam, 1994). The use of more than one hidden layer provides greater flexibility and enables the approximation of complex functions with fewer connection weights in many situations (Sarle, 1994).

Despite numerous studies, no systematic approach has been developed for the selection of an optimal network architecture and its geometry. Hence, we employ a trial and error exercise to determine the network architecture (number of hidden layers and number of nodes for each layer). Initially, a simple three layer MLP network (only one hidden layer) is used to train the network to build the complex relationship of output, *K*, and the associated concentration values. Hecht-Nielsen (1987) suggested that an optimum number of hidden nodes for a single hidden layer network can be selected from following relationship, Number if hidden nodes = 2 ( Number of input nodes ) + 1.

This suggests we should have 101 hidden nodes for 50 inputs. Nevertheless, after a number of trial and error tests, it was found that the optimum results can be achieved by 20 hidden neurons. As reported elsewhere, Abrahart et al. (1999) presented a method based on genetic algorithms to identify the best number of suitable hidden nodes. Chakilam (1998) used principal component analysis to determine the optimal structure of the multi-layered feed forward neural network for a time series forecasting problems, thus reducing the generalisation error and overcoming the over fitting problems.

The dataset is divided into two categories with a random selection of 80% used for training and the rest for testing. The maximum and minimum values of the training network are set by selecting the values from both training (and testing) and estimating dataset, to prevent the ANN from extrapolating beyond its range. We apply scale functions of none, logistic and logistic for input, hidden and output layers, respectively. The default network parameters of NeuroShell2 (neural computing software package) are used; learning rate = 0.1, momentum = 0.1, initial weight = 0.3. The network reached the stopping criterion of average error on test set, fixed at 0.000002, in less than two minutes in a 1GHz personal computer with performance measurements of the coefficient of multiple determination, R2 = 0.9999 and the square of the correlation coefficient, r2 = 0.9999. The network that produces the best results on the test set is the one most capable of generalising, so this is saved as the best network. hidden nodes. network with *K*is Max

Having completed the successful training, another dataset is employed to test the performance of the trained network. We make use of the same model to generate 800 new data values, however, the initial concentration is randomly changed by up to 5%, and up to 5% noise is arbitrarily added to all concentration input values. The reason for adding the noise is to simulate the real world problem of erratic behaviour of aquifers. The estimation error of each *K* value is given in Figure 3.4.4, which shows that the error increases with *K*.

Table 3.3 illustrates that the statistical measurements of error with mean square error (MSE) of 45.25, an average absolute percentage error (AAPE) of 5.63%, and a maximum error of 22.45 m/day. Such high error values may not be acceptable in the most practical cases. Since the objective range of parameters is fairly large (40 –240 m/day), the accuracy of the approximation tends to decrease (Figure 3.24). Therefore, we conduct the same estimation procedure with four

98 in Porous Media - An Approach Based on Stochastic Calculus

1994; Castellano et al., 1997). Cybenko (1989) showed that only one hidden layer was required to approximate any continuous function, given that sufficient degrees of freedom (i.e. connection weights) are provided. However, in practice many functions are difficult to approximate with one hidden layer, requiring a prohibitive number of hidden layer nodes (Flood and Kartam, 1994). The use of more than one hidden layer provides greater flexibility and enables the approximation of complex functions with fewer connection weights in

Despite numerous studies, no systematic approach has been developed for the selection of an optimal network architecture and its geometry. Hence, we employ a trial and error exercise to determine the network architecture (number of hidden layers and number of nodes for each layer). Initially, a simple three layer MLP network (only one hidden layer) is used to train the network to build the complex relationship of output, *K*, and the associated concentration values. Hecht-Nielsen (1987) suggested that an optimum number of hidden nodes for a single hidden layer network can be selected from following relationship,

This suggests we should have 101 hidden nodes for 50 inputs. Nevertheless, after a number of trial and error tests, it was found that the optimum results can be achieved by 20 hidden neurons. As reported elsewhere, Abrahart et al. (1999) presented a method based on genetic algorithms to identify the best number of suitable hidden nodes. Chakilam (1998) used principal component analysis to determine the optimal structure of the multi-layered feed forward neural network for a time series forecasting problems, thus reducing the

The dataset is divided into two categories with a random selection of 80% used for training and the rest for testing. The maximum and minimum values of the training network are set by selecting the values from both training (and testing) and estimating dataset, to prevent the ANN from extrapolating beyond its range. We apply scale functions of none, logistic and logistic for input, hidden and output layers, respectively. The default network parameters of NeuroShell2 (neural computing software package) are used; learning rate = 0.1, momentum = 0.1, initial weight = 0.3. The network reached the stopping criterion of average error on test set, fixed at 0.000002, in less than two minutes in a 1GHz personal computer with performance measurements of the coefficient of multiple determination, R2 = 0.9999 and the square of the correlation coefficient, r2 = 0.9999. The network that produces the best results on the test set is

Having completed the successful training, another dataset is employed to test the performance of the trained network. We make use of the same model to generate 800 new data values, however, the initial concentration is randomly changed by up to 5%, and up to 5% noise is arbitrarily added to all concentration input values. The reason for adding the noise is to simulate the real world problem of erratic behaviour of aquifers. The estimation error of each *K* value is given in Figure 3.4.4, which shows that the error increases with *K*. Table 3.3 illustrates that the statistical measurements of error with mean square error (MSE) of 45.25, an average absolute percentage error (AAPE) of 5.63%, and a maximum error of 22.45 m/day. Such high error values may not be acceptable in the most practical cases. Since the objective range of parameters is fairly large (40 –240 m/day), the accuracy of the approximation tends to decrease (Figure 3.24). Therefore, we conduct the same estimation procedure with four

many situations (Sarle, 1994).

Number if hidden nodes = 2 ( Number of input nodes ) + 1.

generalisation error and overcoming the over fitting problems.

the one most capable of generalising, so this is saved as the best network.

smaller permissible parameter regimes of *K*; (i) 40-90, (ii) 90-140, (iii) 140-190, and (iv) 190-240 m/day. Table 3.3 shows that accuracy of the estimates has improved considerably.

Figure 3.24. Absolute error of estimated parameter, *K,* when considering the whole range, 40–240 m/day.

The maximum error in the 190-240 range has been reduced by about 90% (Figure 3.25 and Table 3.3). Therefore, it is reasonable to assume that if we can gather prior information about the system under consideration, it is possible to obtain more accurate estimates. However, in the real world problems the prior knowledge of the system is limited. Later, we discuss a method to identify the range of parameters by using Self-Organising Maps (SOM). However, before using SOM, we explore the robustness of the ANN estimation models. a 

Figure 3.25. Absolute Error of estimated parameter *K,* only for the range 190 – 240 m/day.

values of all time intervals. Table 3.4 shows the statistics of the estimates. The estimates exhibit a direct relationship to the noise; however, the most of the results are dependable

A Stochastic Model for Hydrodynamic Dispersion 101

Random boundary conditions and an irregular porous structure can result in an erratic distribution of flow paths. Therefore, solute concentration spreads could be highly stochastic. We address this issue by extending the investigation of the robustness by adding different level of randomness to the concentration values. First, the data is generated using the deterministic solutions of equation (3.6.7) for each case and then noise is added randomly to each deterministic concentration value to generate a noisy dataset. For example, to generate up to a ±10% noise component to a deterministic value, *d*, two random

Table 3.5 demonstrates the statistics of the estimates obtained for noisy concentration data.

Max Mean StDev MSE AAPE

Estimates show that the ANN model is stable even for highly stochastic systems.

10 2.29 0.85 1.10 1.98 0.68 20 3.54 1.05 1.19 2.04 0.76 30 5.46 1.74 1.22 2.25 0.81 40 5.88 1.96 1.35 2.31 0.92 50 6.00 1.99 1.39 2.39 0.93

Table 3.5. Statistics of estimates for noisy data for the *K* range, 190-240 m/day.

The ANN model gives more accurate estimates when the target parameter range is small. However, in real world heterogeneous aquifers, it may not be a trivial task to identify the accurate parameter range without reliable prior information. Self Organising Maps (SOM ) has the ability to sort items into categories of similar objects by nonlinearly projecting the data onto a lower dimensional display and by clustering the data (Kohonen, 1990). We use this power of SOM and develop a method to identify the parameter range for given solute concentration values. We employ SOM to cluster an 800 x 50 dimension noisy dataset used before with a parameter range of 40 –240 m/day into four different categories. The "Supervised Kohonen" network architecture of NeroShell2 successfully categorised four different groups with 201, 200, 197 and 202 data patterns in each cluster, respectively. The SOM put data into categories with a high accuracy, with few exceptions which can be expected with noisy data, at the boundaries of the parameter ranges. To test the accuracy of the prediction capability of the trained model, we then create and feed 12 different test

random function 1 generate a random number between 0-1 (say *n*)

Error of Estimate *K* (m/day)

even at higher noise levels.

functions are used as follows,

%

added noise

random function 2 generate either + or -.

Therefore, noisy data = *d* (1 ± 10% *n*).


Table 3.3. Statistics of estimated error for different ranges of *K* with up to a ±5% difference in initial value and up to a ±5% noise in observations.


Table 3.4. Statistics of estimated error for different initial values with up to a ±5% error for the *K* range, 190-240.

As discussed in chapter 1, real world aquifer systems are subject to numerous random effects. One of them may be an initial value problem. First, we investigate in the range of *K* between 190 – 240 m/day for the stability of the model for different initial values. The point source value of the initial concentration are changed from –50% to 50%. Since, as explained above, an exponentially distributed pattern was used to determine the initial concentration values of other spatial points, change of initial value at the point source also resulted in changing the initial values of every spatial location. Furthermore, to illustrate the heterogeneity of the aquifers, up to 5% extra noise is added to the generated concentration subject used Table 3.5. Statistics of estimates for noisy data for the range, 190-240 m/day.

100 in Porous Media - An Approach Based on Stochastic Calculus

Table 3.3. Statistics of estimated error for different ranges of *K* with up to a ±5% difference in

Point C value Noise Maximum Mean Stdev MSE AAPE 50 -50% 9.72 3.48 1.68 12.66 0.95 60 -40% 7.55 2.06 1.57 3.36 0.94 70 -30% 6.18 2.01 1.42 3.01 0.92 80 -20% 4.36 1.59 0.97 2.58 0.77 90 -10% 2.84 0.97 0.72 1.46 0.43

110 +10% 2.92 1.04 0.86 1.46 0.44 120 +20% 4.57 1.68 1.05 2.95 0.84 130 +30% 6.49 2.11 1.48 3.26 1.06 140 +40% 7.58 2.08 1.57 3.47 1.07 150 +50% 10.14 3.67 1.73 13.81 1.07

Table 3.4. Statistics of estimated error for different initial values with up to a ±5% error for

As discussed in chapter 1, real world aquifer systems are subject to numerous random effects. One of them may be an initial value problem. First, we investigate in the range of *K* between 190 – 240 m/day for the stability of the model for different initial values. The point source value of the initial concentration are changed from –50% to 50%. Since, as explained above, an exponentially distributed pattern was used to determine the initial concentration values of other spatial points, change of initial value at the point source also resulted in changing the initial values of every spatial location. Furthermore, to illustrate the heterogeneity of the aquifers, up to 5% extra noise is added to the generated concentration

Max 22.45 1.88 2.23 2.99 2.98 Mean 8.03 0.27 0.38 0.36 0.39 StDev 5.15 0.32 0.41 0.49 0.47 MSE 45.25 0.11 0.14 0.20 0.19 AAPE(%) 5.63 0.11 0.12 0.18 0.18

40-240 40–90 90–140 140–190 190-240

*K* Range (m/day)

initial value and up to a ±5% noise in observations.

Initial condition Error of *K* (m/day)

Error

Trained value

the *K* range, 190-240.

(100)

values of all time intervals. Table 3.4 shows the statistics of the estimates. The estimates exhibit a direct relationship to the noise; however, the most of the results are dependable even at higher noise levels.

Random boundary conditions and an irregular porous structure can result in an erratic distribution of flow paths. Therefore, solute concentration spreads could be highly stochastic. We address this issue by extending the investigation of the robustness by adding different level of randomness to the concentration values. First, the data is generated using the deterministic solutions of equation (3.6.7) for each case and then noise is added randomly to each deterministic concentration value to generate a noisy dataset. For example, to generate up to a ±10% noise component to a deterministic value, *d*, two random functions are used as follows,


The ANN model gives more accurate estimates when the target parameter range is small. However, in real world heterogeneous aquifers, it may not be a trivial task to identify the accurate parameter range without reliable prior information. Self Organising Maps (SOM ) has the ability to sort items into categories of similar objects by nonlinearly projecting the data onto a lower dimensional display and by clustering the data (Kohonen, 1990). We use this power of SOM and develop a method to identify the parameter range for given solute concentration values. We employ SOM to cluster an 800 x 50 dimension noisy dataset used before with a parameter range of 40 –240 m/day into four different categories. The "Supervised Kohonen" network architecture of NeroShell2 successfully categorised four different groups with 201, 200, 197 and 202 data patterns in each cluster, respectively. The SOM put data into categories with a high accuracy, with few exceptions which can be expected with noisy data, at the boundaries of the parameter ranges. To test the accuracy of the prediction capability of the trained model, we then create and feed 12 different test

errors, it is fair to state that the estimate obtained from the ANN model is reasonable and

A Stochastic Model for Hydrodynamic Dispersion 103

Our investigation emphasise that the importance of modelling a sufficiently true representation of the physical system and subsidiary conditions to obtain accurate parameter values. As Minns et al. (1996) pointed out, the ANN is susceptible to becoming "a prisoner of its training data". Therefore, prior information, such as type of contaminant source, boundary conditions and subsidiary conditions is crucial to modelling the system accurately. If we could gain such prior information and model the system with ANN, it would be capable of solving the inverse problem with greater accuracy, even with highly

In this section, we apply the ANN hybrid approach presented in the previous section to estimate parameters of the Stochastic Solute Transport Model (SSTM) developed in this chapter. Since SSTM consists of two parameters ; variance ( σ<sup>2</sup> ) and correlation length (*b*),

In section 3.8 we showed that the accuracy of the estimates was inversely proportional to the size of the objective range of the output (maximum error of the estimate of parameter *K* was reduced by about 90% for smaller output ranges). In addition, the accuracy of the estimates may reduce when two parameters are estimated simultaneously. To limit the diminution of accuracy of the results imposed by the above mentioned performance characteristics of the method and, as this research is conducted in a general personal computer, we select a smaller permissible output ranges. Additionally, it has shown that the higher parameter values of SSTM (especially σ2 around 0.25) represent greater heterogeneous flow systems. Thus, we limit the parameter range for both parameters, σ2 and *b*, to be between 0.0001

The main objective of this section is to estimate parameters of SSTM using a hybrid ANN approach developed in the above section. We then extend the exercise using a case study to validate the method. As used in the previous section, the results from the artificial aquifer is used for the validation. However, we intend to achieve an auxiliary objective in the validation process by comparing the artificial aquifer parameters with the estimates

As the first step of the implementation process, we use SSTM to simulate a one-dimensional aquifer of 10 m in length. Eight hundred data patterns for different combinations of σ2 and *b* are generated. Each parameter ranges from 0.0001 and 0.2. Every data pattern consists of 200 inputs for 10 various spatial locations of the aquifer for 20 distinct time intervals. The same standard Wiener process is used for generating all data sets. The initial condition of the concentration value of unit 1.0 at *x* = 0.0 and exponentially distributed values for other spatial locations are considered. Throughout the simulation the same concentration (unit 1.0) is maintained at the upper end boundary. It is assumed that the mean velocity of the solute is 0.5 m/day. As mentioned above, it is very important to limit the objective range of the parameters for smaller regimes to attain accurate approximations. Thus, Kohonen's Self

Organising Map (SOM) is employed to cluster the data set into four categories.

noisy data, as well as different system input values.

we estimate both parameters simultaneously.

obtained for the same aquifer using a curve fitting technique.

**3.9 Parameter Estimation of SSTM** 

acceptable.

and 0.2.

datasets with the same number of input variables (50) into the trained SOM and it accurately identified the correct parameter range for all the datasets.

We extend the hybrid methodology to solve the groundwater inverse problem in the case of two unknown system parameters. We simulate the same aquifer as used above. The two parameters to be estimated are hydraulic conductivity, *K* (m/day) and longitudinal dispersion coefficient, *DL* (m2/day). In line with earlier work, we fed 50 concentration values and two actual outputs (*K* and *DL*) to train the network.

We use a simple three layer network and produce R2 = 0.9999 and r2 = 0.9999 for both outputs in 2 min and 50 sec. We then feed a different dataset, which has not been seen by the trained network before. The new dataset consisted of randomly varying (up to 5%) initial conditions and added noise to replicate a natural system. We explore two different levels of noise; up to 5% and 50%. The parameter ranges are; *K* between 190 – 240 m/day, *DL* between 0.03 – 0.08 m2/day. The ANN model produces reasonable estimates for both parameters and the summary of estimates is given in Table 3.6.


Table 3.6. Statistics of estimates for two parameter case

We apply the hybrid ANN inverse approach presented in the above sections to estimate parameters of the artificial aquifer at Lincoln University. Although, initial conditions, other parameters and the subsidiary conditions of the aquifer are somewhat known, we have to conduct a fairly tiresome, "trial and error" exercise to replicate the aquifer. Eight hundred data patterns are generated for the hydraulic conductivity range of 80 to 280 m/day. Each pattern consists of 100 concentration input variables for 10 distinct spatial locations for 10 different time intervals. We then use Kohonen's SOM (80% data for training and 20% for testing) to classify the input values into four clusters. Next we feed the actual aquifer data into the trained network and the selected subrange; it is established that the aquifer parameter should be within the second cluster (130 – 180 m/day). Based on this information, we generate a separate dataset for the specified range and train an MLP network with the associated *K* values. subrange; 137 <sup>a</sup>

The estimate of *K* given by the trained ANN is 152.86 m/day. The experimental value of hydraulic conductivity, *K*, is found to be 137 m/day, which is calculated by calibration tests conducted by aquifer testing staff. In these experiments they have assumed that the aquifer is homogeneous. The difference between two estimates is only 10.37 %. Considering the assumptions of homogeneity made by the aquifer researchers and other possible human

Error of Estimate Max Mean MSE AAPE

5 2.48 0.99 2.65 0.81 50 6.78 2.35 3.18 1.12

5 0.00341 0.00092 0.0014 0.0005 50 0.00875 0.00247 0.0029 0.0010

102 in Porous Media - An Approach Based on Stochastic Calculus

datasets with the same number of input variables (50) into the trained SOM and it accurately

We extend the hybrid methodology to solve the groundwater inverse problem in the case of two unknown system parameters. We simulate the same aquifer as used above. The two parameters to be estimated are hydraulic conductivity, *K* (m/day) and longitudinal dispersion coefficient, *DL* (m2/day). In line with earlier work, we fed 50 concentration values

We use a simple three layer network and produce R2 = 0.9999 and r2 = 0.9999 for both outputs in 2 min and 50 sec. We then feed a different dataset, which has not been seen by the trained network before. The new dataset consisted of randomly varying (up to 5%) initial conditions and added noise to replicate a natural system. We explore two different levels of noise; up to 5% and 50%. The parameter ranges are; *K* between 190 – 240 m/day, *DL* between 0.03 – 0.08 m2/day. The ANN model produces reasonable estimates for both

We apply the hybrid ANN inverse approach presented in the above sections to estimate parameters of the artificial aquifer at Lincoln University. Although, initial conditions, other parameters and the subsidiary conditions of the aquifer are somewhat known, we have to conduct a fairly tiresome, "trial and error" exercise to replicate the aquifer. Eight hundred data patterns are generated for the hydraulic conductivity range of 80 to 280 m/day. Each pattern consists of 100 concentration input variables for 10 distinct spatial locations for 10 different time intervals. We then use Kohonen's SOM (80% data for training and 20% for testing) to classify the input values into four clusters. Next we feed the actual aquifer data into the trained network and the selected subrange; it is established that the aquifer parameter should be within the second cluster (130 – 180 m/day). Based on this information, we generate a separate dataset for the specified range and train an MLP

The estimate of *K* given by the trained ANN is 152.86 m/day. The experimental value of hydraulic conductivity, *K*, is found to be 137 m/day, which is calculated by calibration tests conducted by aquifer testing staff. In these experiments they have assumed that the aquifer is homogeneous. The difference between two estimates is only 10.37 %. Considering the assumptions of homogeneity made by the aquifer researchers and other possible human

identified the correct parameter range for all the datasets.

and two actual outputs (*K* and *DL*) to train the network.

parameters and the summary of estimates is given in Table 3.6.

Parameter Actual Range % noise

Table 3.6. Statistics of estimates for two parameter case

K 190-240

*DL* 0.03-0.08

network with the associated *K* values.

errors, it is fair to state that the estimate obtained from the ANN model is reasonable and acceptable.

Our investigation emphasise that the importance of modelling a sufficiently true representation of the physical system and subsidiary conditions to obtain accurate parameter values. As Minns et al. (1996) pointed out, the ANN is susceptible to becoming "a prisoner of its training data". Therefore, prior information, such as type of contaminant source, boundary conditions and subsidiary conditions is crucial to modelling the system accurately. If we could gain such prior information and model the system with ANN, it would be capable of solving the inverse problem with greater accuracy, even with highly noisy data, as well as different system input values.

### **3.9 Parameter Estimation of SSTM**

In this section, we apply the ANN hybrid approach presented in the previous section to estimate parameters of the Stochastic Solute Transport Model (SSTM) developed in this chapter. Since SSTM consists of two parameters ; variance ( σ<sup>2</sup> ) and correlation length (*b*), we estimate both parameters simultaneously.

In section 3.8 we showed that the accuracy of the estimates was inversely proportional to the size of the objective range of the output (maximum error of the estimate of parameter *K* was reduced by about 90% for smaller output ranges). In addition, the accuracy of the estimates may reduce when two parameters are estimated simultaneously. To limit the diminution of accuracy of the results imposed by the above mentioned performance characteristics of the method and, as this research is conducted in a general personal computer, we select a smaller permissible output ranges. Additionally, it has shown that the higher parameter values of SSTM (especially σ2 around 0.25) represent greater heterogeneous flow systems. Thus, we limit the parameter range for both parameters, σ2 and *b*, to be between 0.0001 and 0.2.

The main objective of this section is to estimate parameters of SSTM using a hybrid ANN approach developed in the above section. We then extend the exercise using a case study to validate the method. As used in the previous section, the results from the artificial aquifer is used for the validation. However, we intend to achieve an auxiliary objective in the validation process by comparing the artificial aquifer parameters with the estimates obtained for the same aquifer using a curve fitting technique.

As the first step of the implementation process, we use SSTM to simulate a one-dimensional aquifer of 10 m in length. Eight hundred data patterns for different combinations of σ2 and *b* are generated. Each parameter ranges from 0.0001 and 0.2. Every data pattern consists of 200 inputs for 10 various spatial locations of the aquifer for 20 distinct time intervals. The same standard Wiener process is used for generating all data sets. The initial condition of the concentration value of unit 1.0 at *x* = 0.0 and exponentially distributed values for other spatial locations are considered. Throughout the simulation the same concentration (unit 1.0) is maintained at the upper end boundary. It is assumed that the mean velocity of the solute is 0.5 m/day. As mentioned above, it is very important to limit the objective range of the parameters for smaller regimes to attain accurate approximations. Thus, Kohonen's Self Organising Map (SOM) is employed to cluster the data set into four categories. assumptions 

Gaussian and logistic are used for layers of input, hidden (3 layers) and output, respectively. The default network parameters (NeuroShell2) are employed; learning rate = 0.1, momentum = 0.1, initial weight = 0.3. The stopping criterion is set to a minimum error of 0.000001. All four networks reach the stopping condition in about 30 minutes in a 1 GHz

A Stochastic Model for Hydrodynamic Dispersion 105

After the completion of successful training for each model, separate datasets are generated to test the prediction capability of each model. The same SSTM is employed to produce another dataset of 441 data patterns for each parameter range. However, different standard Wiener process increments are used. In addition, initial and boundary conditions are adjusted up to ±5% by adding random values. Input data values of each dataset are then fed

Square of the correlation

Mean absolute

error

coefficient, r2

σ<sup>2</sup> b σ<sup>2</sup> b σ<sup>2</sup> b

personal computer with the performance measurements shown in Table 3.7.

into the corresponding trained network and processed to obtain model predictions.

0.0001 – 0.05 0.9912 0.9876 0. 9911 0.9876 0.0 0.0 0.05 – 0.1 0.9911 0.9876 0. 9899 0.9876 0.0 0.0 0.1 – 0.15 0.9898 0.9870 0. 9872 0.9868 0.0 0.0 0.15 – 0.2 0.9728 0.9774 0. 9721 0.9661 0.0001 0.0001 Table 3.7. Performance measurements of trained ANN model for four different parameter

Figure 3.26 illustrates the absolute error of estimated parameter σ<sup>2</sup> *,* for the range of 0.0001 – 0.05. It shows that the ANN model prediction is extremely satisfactory and that the average absolute error is approximately 0.04%. Figure 3.27 shows that prediction for the other

Figure 3.26. Absolute error of estimated parameter σ<sup>2</sup> *,* for the range of 0.0001 – 0.05.

Coefficient of multiple

parameter, *b* also met with similar accuracy for the same range.

determination, R2

Parameter range

ranges.

Since the data set represents the stochastic behaviour of the flow, the time needed to classify the data into separate groups is much more than for the similar case, in the deterministic advection – dispersion data, in the previous section. Randomly selected 80% of data were used for training the network and the rest for the validation. Notwithstanding the random nature of the current dataset, the SOM has clustered the data to an adequate degree of accuracy that may be sufficient for the problem at hand.

To test the performance of the trained network model, we use eight different data patterns. Two datasets are generated for each group of clustered network. Moreover, each dataset is produced using different standard Wiener process. New data patterns are fed into the trained model and it accurately identified the parameter range it fitted into.

Having substantiated that the SOM could be successfully used to cluster large dataset that represent heterogeneous data, such as given by SSTM, we create four separate specialised network models for smaller parameter regimes of both parameters in the following ranges;

$$0.0001 - 0.05 \text{\AA} 0.05 - 0.1 \text{\AA} 0.1 - 0.15 \text{\AA} \text{ and } 0.15 - 0.2 \text{\AA}$$

Each dataset comprises 441 data patterns of different combinations of σ2 and *b* at intervals of 0.0025*.* For instance, for the range of 0.05 – 0.1, there are 21 different values of each parameter; 0.05, 0.0525, 0.055, 0.0575, … , 0.1 (21 x 21 = 441 data patterns). Same SSTM model used above for SOM cluster distribution is also used for data generation. Nevertheless, in this case each data pattern contains not only 200 inputs but also two corresponding output parameters. The number of training patterns has considerable influence on the performance of the ANN model (Flood and Kartam, 1994). Increasing the number of data patterns provides more information about the shape of the solution surface and, thus, improves the accuracy of the model prediction. However, in most real world applications, numerous logistical issues impose limitations on the amount of data available and, consequently, the size of the training set. Hence, in developing a method for practical applications, it is important to test the robustness of the method for such data limitations. For that reason, we limit the parameter range of each model to 21 values (between 0.0001 and 0.2 at intervals of 0.0025). the 1>,

Maier and Dandy (2000) showed that it was important to select a suitable network architecture and model validation method in the development of ANN models to achieve optimum results. In addition, it may be necessary to select the most suitable model for handling highly random data such as SSTM data. Therefore, we conduct a few trial and error exercises to choose the appropriate model structure, and training and testing procedure. The following are the MLP ANN models that are considered for different combinations of hidden layers and various grouping of activation functions: three layer standard connections; four layer standard connections; five layer standard connections; one hidden layer with two parallel slabs with different activation functions; one hidden layer with three hidden slabs with different activation functions; one hidden layer of two parallel slabs, different activation functions and jump connection; three layer jump connections; four layer jump connections; and five layer jump connections.

After numerous attempts, it is found that the network model of five layer standard connections could produce the best trained model in the least time. In the selected model, each hidden layer consists of 30 neurons. Activation functions of linear <0, 1>, logistic, tanh,

104 in Porous Media - An Approach Based on Stochastic Calculus

Since the data set represents the stochastic behaviour of the flow, the time needed to classify the data into separate groups is much more than for the similar case, in the deterministic advection – dispersion data, in the previous section. Randomly selected 80% of data were used for training the network and the rest for the validation. Notwithstanding the random nature of the current dataset, the SOM has clustered the data to an adequate degree of

To test the performance of the trained network model, we use eight different data patterns. Two datasets are generated for each group of clustered network. Moreover, each dataset is produced using different standard Wiener process. New data patterns are fed into the

Having substantiated that the SOM could be successfully used to cluster large dataset that represent heterogeneous data, such as given by SSTM, we create four separate specialised network models for smaller parameter regimes of both parameters in the following ranges; 0.0001 – 0.05;0.05 – 0.1;0.1 – 0.15; and 0.15 – 0.2. Each dataset comprises 441 data patterns of different combinations of σ2 and *b* at intervals of 0.0025*.* For instance, for the range of 0.05 – 0.1, there are 21 different values of each parameter; 0.05, 0.0525, 0.055, 0.0575, … , 0.1 (21 x 21 = 441 data patterns). Same SSTM model used above for SOM cluster distribution is also used for data generation. Nevertheless, in this case each data pattern contains not only 200 inputs but also two corresponding output parameters. The number of training patterns has considerable influence on the performance of the ANN model (Flood and Kartam, 1994). Increasing the number of data patterns provides more information about the shape of the solution surface and, thus, improves the accuracy of the model prediction. However, in most real world applications, numerous logistical issues impose limitations on the amount of data available and, consequently, the size of the training set. Hence, in developing a method for practical applications, it is important to test the robustness of the method for such data limitations. For that reason, we limit the parameter range of each model to 21 values (between 0.0001

Maier and Dandy (2000) showed that it was important to select a suitable network architecture and model validation method in the development of ANN models to achieve optimum results. In addition, it may be necessary to select the most suitable model for handling highly random data such as SSTM data. Therefore, we conduct a few trial and error exercises to choose the appropriate model structure, and training and testing procedure. The following are the MLP ANN models that are considered for different combinations of hidden layers and various grouping of activation functions: three layer standard connections; four layer standard connections; five layer standard connections; one hidden layer with two parallel slabs with different activation functions; one hidden layer with three hidden slabs with different activation functions; one hidden layer of two parallel slabs, different activation functions and jump connection; three layer jump connections; four

After numerous attempts, it is found that the network model of five layer standard connections could produce the best trained model in the least time. In the selected model, each hidden layer consists of 30 neurons. Activation functions of linear <0, 1>, logistic, tanh,

trained model and it accurately identified the parameter range it fitted into.

accuracy that may be sufficient for the problem at hand.

and 0.2 at intervals of 0.0025).

layer jump connections; and five layer jump connections.

Gaussian and logistic are used for layers of input, hidden (3 layers) and output, respectively. The default network parameters (NeuroShell2) are employed; learning rate = 0.1, momentum = 0.1, initial weight = 0.3. The stopping criterion is set to a minimum error of 0.000001. All four networks reach the stopping condition in about 30 minutes in a 1 GHz personal computer with the performance measurements shown in Table 3.7.

After the completion of successful training for each model, separate datasets are generated to test the prediction capability of each model. The same SSTM is employed to produce another dataset of 441 data patterns for each parameter range. However, different standard Wiener process increments are used. In addition, initial and boundary conditions are adjusted up to ±5% by adding random values. Input data values of each dataset are then fed into the corresponding trained network and processed to obtain model predictions.


Table 3.7. Performance measurements of trained ANN model for four different parameter ranges.

Figure 3.26 illustrates the absolute error of estimated parameter σ<sup>2</sup> *,* for the range of 0.0001 – 0.05. It shows that the ANN model prediction is extremely satisfactory and that the average absolute error is approximately 0.04%. Figure 3.27 shows that prediction for the other parameter, *b* also met with similar accuracy for the same range.

Figure 3.26. Absolute error of estimated parameter σ<sup>2</sup> *,* for the range of 0.0001 – 0.05.

We generate two separate datasets for two extreme cases:

0.0001 – 0.05.

absolute error is approximately 4%.

and *b* range of 0.15 – 0.2.

range of 0.0001 – 0.05.

where σ2 is smaller and *b* is higher - σ2 ranges of 0.0001 – 0.05 and *b* ranges of 0.15 – 0.2; and where σ2 is higher and *b* is smaller - σ2 ranges of 0.15 – 0.2 and *b* ranges of

A Stochastic Model for Hydrodynamic Dispersion 107

A similar method to that which was used for earlier investigations is employed to gauge the capability of the ANN model. Figures 3.29 and 3.30 reveal that the trained network has predicted the estimates with reasonable precision. In both cases the percentage average

Figure 3.29. Absolute error of estimated parameter σ<sup>2</sup> *,* for the σ2 range of 0.0001 – 0.05

Figure 3.30. Absolute error of estimated parameter b*,* for the σ2 range of 0.15 – 0.2 and *b*

We use the artificial aquifer data to validate the developed hybrid ANN method. We make use of the same aquifer dataset in this section to estimate the parameters of SSTM for the

A similar approach is also applied to other parameter ranges. The precision of the estimates given by ANN models shrinks with highly heterogeneous data. As larger values of parameters indicate excessive stochastic flows, we can expect the accuracy of the prediction to diminish for highly stochastic flows. Nonetheless, the average absolute error for the estimates for a range of 0.15 to 0.2 is approximately 5.5% (for example, Figure 3.28 illustrates the error of estimated parameter σ2 for parameter range of 0.15 to 0.2), which may be acceptable for the most of practical applications.

The above prediction accuracy analyses of the ANN models are based on similar ranges for both parameters. However, in real world applications we may have to deal with extremely different values for two parameters. Therefore, the robustness of the ANN method for different values of parameter regimes for two parameters need to be assessed.

Figure 3.28. Absolute error of estimated parameter σ<sup>2</sup> *,* for the range of 0.15– 0.2.

106 in Porous Media - An Approach Based on Stochastic Calculus

A similar approach is also applied to other parameter ranges. The precision of the estimates given by ANN models shrinks with highly heterogeneous data. As larger values of parameters indicate excessive stochastic flows, we can expect the accuracy of the prediction to diminish for highly stochastic flows. Nonetheless, the average absolute error for the estimates for a range of 0.15 to 0.2 is approximately 5.5% (for example, Figure 3.28 illustrates the error of estimated parameter σ2 for parameter range of 0.15 to 0.2), which may be

The above prediction accuracy analyses of the ANN models are based on similar ranges for both parameters. However, in real world applications we may have to deal with extremely different values for two parameters. Therefore, the robustness of the ANN method for

different values of parameter regimes for two parameters need to be assessed.

Figure 3.27. Absolute error of estimated parameter *b,* for the range of 0.0001 – 0.05.

Figure 3.28. Absolute error of estimated parameter σ<sup>2</sup> *,* for the range of 0.15– 0.2.

acceptable for the most of practical applications.

We generate two separate datasets for two extreme cases:

where σ2 is smaller and *b* is higher - σ2 ranges of 0.0001 – 0.05 and *b* ranges of 0.15 – 0.2; and where σ2 is higher and *b* is smaller - σ2 ranges of 0.15 – 0.2 and *b* ranges of 0.0001 – 0.05.

A similar method to that which was used for earlier investigations is employed to gauge the capability of the ANN model. Figures 3.29 and 3.30 reveal that the trained network has predicted the estimates with reasonable precision. In both cases the percentage average absolute error is approximately 4%.

Figure 3.29. Absolute error of estimated parameter σ<sup>2</sup> *,* for the σ2 range of 0.0001 – 0.05 and *b* range of 0.15 – 0.2.

Figure 3.30. Absolute error of estimated parameter b*,* for the σ2 range of 0.15 – 0.2 and *b* range of 0.0001 – 0.05.

We use the artificial aquifer data to validate the developed hybrid ANN method. We make use of the same aquifer dataset in this section to estimate the parameters of SSTM for the method.

(*D*) using the maximum likelihood estimation procedure for the 1-D advection-dispersion equation for a given velocity using each of the realization; and (3) take the mean of the

A Stochastic Model for Hydrodynamic Dispersion 109

We have demonstrated that the SSTM can be used to characterise an experimental homogeneous sand aquifer of 5 m width x 10 m length x 2.7 height using a single set of values of ( *σ2* = 0.01 and *b* = 0.01) quite satisfactorily. It has also been shown that the numerical solutions of SSTM, in conjunction with the parameter estimation methods such as maximum likelihood method and artificial neural networks, can be used to estimate reliable effective dispersion coefficients, therefore, effective dispersivity, for different scale

We use this procedure, which we call stochastic inversive method (SIM) to estimate *D* for

and a different kernel would give different D values. In addition, SIM is not accurate for very noisy realizations because the theory of estimation is strictly valid when noise is small

0.0001 0.0001 0.01634 0.03502 0.05804 0.07328 0.09447 0.12493 0.0001 0.001 0.01942 0.03738 0.06126 0.07526 0.09930 0.13006 0.0001 0.01 0.03844 0.05287 0.08469 0.10502 0.13591 0.18294 0.0001 0.05 0.04758 0.07799 0.12986 0.16064 0.20645 0.27009 0.0001 0.1 0.06786 0.08690 0.14214 0.17914 0.23342 0.31205 0.0001 0.15 0.06968 0.09060 0.15044 0.18790 0.24022 0.31327 0.0001 0.2 0.07047 0.09259 0.15717 0.20209 0.25756 0.33937 0.0001 0.25 0.07188 0.09382 0.15905 0.19747 0.25542 0.33289 0.0001 0.3 0.07258 0.09466 0.16131 0.19895 0.26026 0.34742 0.001 0.0001 0.01917 0.03738 0.05519 0.07807 0.09628 0.12984 0.001 0.001 0.03749 0.05289 0.08020 0.11085 0.13335 0.18273 0.001 0.01 0.06739 0.08698 0.13530 0.18702 0.23005 0.31770 0.001 0.05 0.07424 0.09650 0.15306 0.21656 0.26236 0.36391 0.001 0.1 0.08492 0.09910 0.15522 0.21291 0.25431 0.35222

1 m 10 m 20 m 30 m 50 m 100 m

dispersion coefficients (*D*) for the range of scales; 1, 10, 20, 30, 50 and 100 m.

<sup>b</sup>σ<sup>2</sup> Estimated *<sup>D</sup>*

We need to remember that the SSTM is based on the velocity covariance kernel,

and *b* for different flow lengths. Table 3.8 exhibits the estimated

1 2

small

,

*x x b*

2

 *e* 

estimates of *D* as the dispersivity for the given set of parameters

experiments up to 10 meters.

and Gaussian (Kutoyants, 1984).

other combinations of <sup>2</sup>

aquifer. Additionally, we use the same data to partially validate the stochastic model (SSTM) previously using a curve fitting technique to approximate the aquifer parameters. It was found that the approximate artificial aquifer parameters are variance, σ2 = 0.01 and correlation length*, b* = 0.01.

In the present validation, first, we used the known conditions such as initial and boundary conditions, and hydraulic gradient to simulate the aquifer using the SSTM. The initial concentration at *x* = 0 is 1.0 unit and it is reduced exponentially with time. Initial values of other spatial points are considered as zero. 1681 data patterns are generated for different combinations of parameters, σ2 and *b*. Single standard Wiener process increments are retained for every simulation run. Both parameters varied between 0.0001 and 0.2.

The 1681 generated data patterns are fed into Kohonen's self-organising map architecture to cluster them into four different groups. After classifying them with reasonable accuracy the aquifer dataset is fed into the SOM model to identify relevant groups that the data resemble the most. In this case, we make an assumption that the effect of the transverse dispersion of the flow is reflected in the stochastic flow described by SSTM. As shown in Figure 3.13, the data have been collected along five wells at four levels (A to E wells at levels BR, BL, RE and YE). We consider the data collected along each well at a certain level as a one-dimensional flow path. Hence, we reproduce 20 different one-dimensional datasets.

Concentration values of the aquifer are normalised to enable to weigh them against the normalised SSTM data.

Initially, a dataset closer to the middle of the aquifer, well C – level RE is chosen for the estimation of parameters. An aquifer dataset has to be interpolated to produce missing data and fabricate uniform spatial and temporal grids. Having constructed an exact number of data for similar spatial and time intervals as for the original ANN model, the aquifer dataset is fed into the trained model. The selected dataset is then separated from the larger set and trained for the smaller range selected. Based on the findings described previously, a five layer standard connections is used with the same activation functions, initial weights, momentum and learning rates.

After completing sufficient training, the artificial aquifer dataset is fed into the ANN model to estimate parameters. The estimates that are produced by the model are σ2 = 0.01364 and *b* = 0.01665. The estimates produced by the ANN model show close resemblance to the values given by the curve fitting technique. Since we avoid considering lateral dispersions for estimation by ANN model, the results may be subject to slight errors.

We extend the estimation procedure to determine the parameters of the Lincoln University aquifer for other flow lines. Since the earlier dataset is closer to the middle of the aquifer, we choose the next dataset that is nearer to the boundary of the aquifer. Estimates obtained using well A – level YE are σ2 = 0.01483 and *b* = 0.00912. These estimates are similar to earlier values.

### **3.10 Dispersivity Based on the SSTM**

To estimate the parameter of the SSTM, we develop a procedure consisting of the following steps: (1) we generate a large number of realizations of concentration (usually 100) for a particular set of values of σ2 and b using the SSTM; (2) estimate the diffusion coefficient

108 in Porous Media - An Approach Based on Stochastic Calculus

aquifer. Additionally, we use the same data to partially validate the stochastic model (SSTM) previously using a curve fitting technique to approximate the aquifer parameters. It was found that the approximate artificial aquifer parameters are variance, σ2 = 0.01 and

In the present validation, first, we used the known conditions such as initial and boundary conditions, and hydraulic gradient to simulate the aquifer using the SSTM. The initial concentration at *x* = 0 is 1.0 unit and it is reduced exponentially with time. Initial values of other spatial points are considered as zero. 1681 data patterns are generated for different combinations of parameters, σ2 and *b*. Single standard Wiener process increments are

The 1681 generated data patterns are fed into Kohonen's self-organising map architecture to cluster them into four different groups. After classifying them with reasonable accuracy the aquifer dataset is fed into the SOM model to identify relevant groups that the data resemble the most. In this case, we make an assumption that the effect of the transverse dispersion of the flow is reflected in the stochastic flow described by SSTM. As shown in Figure 3.13, the data have been collected along five wells at four levels (A to E wells at levels BR, BL, RE and YE). We consider the data collected along each well at a certain level as a one-dimensional

Concentration values of the aquifer are normalised to enable to weigh them against the

Initially, a dataset closer to the middle of the aquifer, well C – level RE is chosen for the estimation of parameters. An aquifer dataset has to be interpolated to produce missing data and fabricate uniform spatial and temporal grids. Having constructed an exact number of data for similar spatial and time intervals as for the original ANN model, the aquifer dataset is fed into the trained model. The selected dataset is then separated from the larger set and trained for the smaller range selected. Based on the findings described previously, a five layer standard connections is used with the same activation functions, initial weights,

After completing sufficient training, the artificial aquifer dataset is fed into the ANN model to estimate parameters. The estimates that are produced by the model are σ2 = 0.01364 and *b* = 0.01665. The estimates produced by the ANN model show close resemblance to the values given by the curve fitting technique. Since we avoid considering lateral dispersions

We extend the estimation procedure to determine the parameters of the Lincoln University aquifer for other flow lines. Since the earlier dataset is closer to the middle of the aquifer, we choose the next dataset that is nearer to the boundary of the aquifer. Estimates obtained using well A – level YE are σ2 = 0.01483 and *b* = 0.00912. These estimates are similar to earlier values.

To estimate the parameter of the SSTM, we develop a procedure consisting of the following steps: (1) we generate a large number of realizations of concentration (usually 100) for a particular set of values of σ2 and b using the SSTM; (2) estimate the diffusion coefficient

retained for every simulation run. Both parameters varied between 0.0001 and 0.2.

flow path. Hence, we reproduce 20 different one-dimensional datasets.

for estimation by ANN model, the results may be subject to slight errors.

correlation length*, b* = 0.01.

normalised SSTM data.

momentum and learning rates.

**3.10 Dispersivity Based on the SSTM** 

(*D*) using the maximum likelihood estimation procedure for the 1-D advection-dispersion equation for a given velocity using each of the realization; and (3) take the mean of the estimates of *D* as the dispersivity for the given set of parameters

We have demonstrated that the SSTM can be used to characterise an experimental homogeneous sand aquifer of 5 m width x 10 m length x 2.7 height using a single set of values of ( *σ2* = 0.01 and *b* = 0.01) quite satisfactorily. It has also been shown that the numerical solutions of SSTM, in conjunction with the parameter estimation methods such as maximum likelihood method and artificial neural networks, can be used to estimate reliable effective dispersion coefficients, therefore, effective dispersivity, for different scale experiments up to 10 meters.

We use this procedure, which we call stochastic inversive method (SIM) to estimate *D* for other combinations of <sup>2</sup> and *b* for different flow lengths. Table 3.8 exhibits the estimated dispersion coefficients (*D*) for the range of scales; 1, 10, 20, 30, 50 and 100 m.

We need to remember that the SSTM is based on the velocity covariance kernel, 1 2 2 *x x b e* , and a different kernel would give different D values. In addition, SIM is not accurate for very noisy realizations because the theory of estimation is strictly valid when noise is small and Gaussian (Kutoyants, 1984).


Table 3.8. Estimates of *D* obtained by using a stochastic inverse method for different combinations of parameters of SSTM for different flow lengths (velocity = 0.5 m/day).

0.1 0.3 0.07678 0.29791 0.48068 0.67575 0.91363 1.41871 0.15 0.0001 0.01731 0.08784 0.14149 0.20423 0.28678 0.45299 0.15 0.001 0.02570 0.10921 0.18047 0.25898 0.35905 0.56336 0.15 0.01 0.05454 0.16218 0.26117 0.37775 0.52802 0.82934 0.15 0.05 0.06574 0.24498 0.40463 0.57946 0.81625 1.27853 0.15 0.1 0.07627 0.24830 0.40532 0.58494 0.82098 1.28586 0.15 0.15 0.07722 0.24896 0.40327 0.58331 0.82088 1.29136 0.15 0.2 0.07785 0.25016 0.40435 0.58084 0.81027 1.27489 0.15 0.25 0.07807 0.25278 0.41587 0.59926 0.83683 1.32106 0.2 0.3 0.07527 0.27394 0.44577 0.65306 0.93236 1.47560 0.25 0.0001 0.01694 0.08642 0.14037 0.20553 0.29769 0.47547 0.25 0.001 0.02292 0.11779 0.19345 0.28363 0.41296 0.65983 0.25 0.01 0.04829 0.13433 0.22093 0.32259 0.46790 0.74928 0.25 0.05 0.06187 0.18828 0.30661 0.45919 0.67190 1.08123 0.25 0.1 0.07470 0.23226 0.37711 0.55761 0.81495 1.30630 0.25 0.15 0.07537 0.25020 0.41219 0.60827 0.88655 1.42797 025 0.2 0.07559 0.26583 0.42875 0.63902 0.92617 1.48670 0.25 0.25 0.07567 0.28102 0.46074 0.68281 0.99764 1.60438 0.25 0.3 0.07593 0.29162 0.46948 0.69110 1.00010 1.61016 0.3 0.0001 0.01785 0.08229 0.13906 0.21167 0.31537 0.51395 0.3 0.001 0.02169 0.11016 0.18535 0.28209 0.42416 0.69259 0.3 0.01 0.04787 0.14860 0.25288 0.38448 0.58112 0.95251 0.3 0.05 0.06316 0.22140 0.38149 0.57637 0.87271 1.42535 0.3 0.1 0.07527 0.25171 0.42778 0.64452 0.97782 1.58726 0.3 0.15 0.07728 0.27513 0.47773 0.73238 1.10681 1.80695 0.3 0.2 0.07799 0.28755 0.49003 0.73771 1.11325 1.81448 0.3 0.25 0.07824 0.29803 0.50666 0.77144 1.17052 1.91032 0.3 0.3 0.07911 0.30699 0.53747 0.80802 1.22297 1.98500

0.25 0.001 0.02292 0.11779 0.19345 0.28363 0.41296 0.65983

A Stochastic Model for Hydrodynamic Dispersion 111


110 in Porous Media - An Approach Based on Stochastic Calculus

0.001 0.15 0.08671 0.10782 0.17443 0.24632 0.29578 0.40749 0.001 0.2 0.09482 0.11886 0.18858 0.26626 0.32559 0.44440 0.001 0.25 0.10701 0.12933 0.20057 0.28264 0.33817 0.46265 0.001 0.3 0.11008 0.13860 0.21898 0.30438 0.36830 0.50826 0.01 0.0001 0.02541 0.05302 0.08035 0.11327 0.14376 0.20657 0.01 0.001 0.05697 0.08773 0.13444 0.19613 0.24574 0.34753 0.01 0.01 0.07906 0.10003 0.15332 0.22578 0.28589 0.40472 0.01 0.05 0.09457 0.16952 0.26730 0.38284 0.48728 0.69328 0.01 0.1 0.11483 0.21334 0.33450 0.48658 0.61670 0.87920 0.01 0.15 0.13745 0.24065 0.37173 0.54195 0.67787 0.96664 0.01 0.2 0.15574 0.25994 0.40478 0.58780 0.73802 1.05102 0.01 0.25 0.18468 0.27120 0.41747 0.59966 0.75083 1.06491 0.01 0.3 0.18994 0.27751 0.43013 0.62369 0.78586 1.11591 0.05 0.0001 0.01874 0.07828 0.12467 0.17543 0.22998 0.34209 0.05 0.001 0.03559 0.10128 0.16245 0.22834 0.29946 0.44681 0.05 0.01 0.06957 0.15994 0.25211 0.34326 0.45352 0.68013 0.05 0.05 0.07651 0.26271 0.41319 0.56701 0.75557 1.13054 0.05 0.1 0.08450 0.29209 0.45872 0.63193 0.83988 1.26490 0.05 0.15 0.08725 0.29659 0.47462 0.66040 0.88362 1.33137 0.05 0.2 0.08987 0.29700 0.46840 0.64492 0.85205 1.27735 0.05 0.25 0.09219 0.29480 0.46693 0.65051 0.87047 1.30991 0.05 0.3 0.09294 0.29137 0.45748 0.63153 0.83935 1.26428 0.1 0.0001 0.01797 0.08546 0.14242 0.19805 0.26901 0.42102 0.1 0.001 0.02971 0.10536 0.16856 0.23555 0.31622 0.48560 0.1 0.01 0.05925 0.17795 0.29101 0.41507 0.56815 0.88488 0.1 0.05 0.06452 0.28179 0.45334 0.64813 0.87743 1.35846 0.1 0.1 0.07232 0.29698 0.48594 0.68641 0.93357 1.45105 0.1 0.15 0.07485 0.29776 0.47662 0.67992 0.91941 1.42114 0.1 0.2 0.07587 0.29876 0.48105 0.68887 0.93379 1.44814 0.1 0.25 0.07608 0.29928 0.48535 0.68808 0.93290 1.45029


Table 3.8. Estimates of *D* obtained by using a stochastic inverse method for different combinations of parameters of SSTM for different flow lengths (velocity = 0.5 m/day).

Figure 3.33. *D* for the combination of b = 0.3 and σ2 = 0.0001.

A Stochastic Model for Hydrodynamic Dispersion 113

Figure 3.34. Field measured values of longitudinal dispersivity as a function of the scale of measurement. The largest circles represent the most reliable data. The estimated dispersivity

Figures 3.35 – 3.37 demonstrate that corresponding *D* values obtained for SSTM parameters

estimated *D* s are in reasonable agreement with the most of the reliable field measurements

*<sup>L</sup> = 0.1 L.* However, for mid and larger ranges of *b*,

from the SSTM are given by the squares. (Source: Gelhar (1986).)

do not agree with the relationship of

We can summarise some of the estimates in the plots in Figures 3.31, 3.32 and 3.33 as functions of the scales of the experiments. What these plots show is that, for a given set of parameters, the SSTM would give the estimates of dispersivity increasing with the flow length.

As discussed in chapter 1, Pickens and Grisak (1981), and Lallemand-Barres and Peaudecerf (1978, cited in Fetter, 1999) showed that the scale dependency of *<sup>L</sup>* has a linear relationship of *<sup>L</sup> = 0.1 L*, where *L* is the mean travel distance. However, Pickens and Grisak (1981) recognised that the linear increase of dispersivity with the mean travel distance was unlikely for large travel distances. It was expected that tracer migration between aquifer layers could cause a reduction in the magnitude of the proportionality constants, since the transverse migration would tend to reduce the spreading effect caused by the stratification. Field measurements obtained by Gelhar (1986) illustrate that the scale dependence relationship between *<sup>L</sup>* and the flow length is non-linear (Figure 3.34).

To evaluate the comparative estimates of *D* obtained from the inverse method for the SSTM parameters and the field measurements observed by Gelhar (1986), we plot them on the same graph (Figure 3.35 – 3.37). Only reliable observations of Figure 3.34 (indicated by larger symbols) are considered. Since the parameter estimated from the inverse approach is *D*, *<sup>L</sup>* values of Figure 3.34 are converted to *D* ( *D v <sup>L</sup>* ). Furthermore, we plot the relationship of *<sup>L</sup> = 0.1 L* in the same graph to assess our estimates. Three different ranges of *b* are chosen. Figure 3.35 shows the estimates for smaller *b*, 0.0001 m, for four values of σ2 (0.0001, 0.05, 0.2 and 0.3). Figures 3.36 and 3.37 illustrate the similar σ2 values for a mid range value *b*, 0.01 m, and larger *b*, 0.3, respectively.

Figure 3.31. *D* for the parameter combination of *b* = 0.0001 and <sup>2</sup> = 0.0001.

Figure 3.32 *D* for the parameter combination of *b* = 0.001 and <sup>2</sup> = 0.0001.

*<sup>L</sup>* ). Furthermore, we plot the

*<sup>L</sup>* has a linear

112 in Porous Media - An Approach Based on Stochastic Calculus

We can summarise some of the estimates in the plots in Figures 3.31, 3.32 and 3.33 as functions of the scales of the experiments. What these plots show is that, for a given set of parameters,

As discussed in chapter 1, Pickens and Grisak (1981), and Lallemand-Barres and Peaudecerf

Grisak (1981) recognised that the linear increase of dispersivity with the mean travel distance was unlikely for large travel distances. It was expected that tracer migration between aquifer layers could cause a reduction in the magnitude of the proportionality constants, since the transverse migration would tend to reduce the spreading effect caused by the stratification. Field measurements obtained by Gelhar (1986) illustrate that the scale

To evaluate the comparative estimates of *D* obtained from the inverse method for the SSTM parameters and the field measurements observed by Gelhar (1986), we plot them on the same graph (Figure 3.35 – 3.37). Only reliable observations of Figure 3.34 (indicated by larger symbols) are considered. Since the parameter estimated from the inverse approach is

of *b* are chosen. Figure 3.35 shows the estimates for smaller *b*, 0.0001 m, for four values of σ2 (0.0001, 0.05, 0.2 and 0.3). Figures 3.36 and 3.37 illustrate the similar σ2 values for a

*<sup>L</sup> = 0.1 L*, where *L* is the mean travel distance. However, Pickens and

*<sup>L</sup>* and the flow length is non-linear (Figure 3.34).

= 0.0001.

= 0.0001.

*<sup>L</sup> = 0.1 L* in the same graph to assess our estimates. Three different ranges

the SSTM would give the estimates of dispersivity increasing with the flow length.

(1978, cited in Fetter, 1999) showed that the scale dependency of

*<sup>L</sup>* values of Figure 3.34 are converted to *D* ( *D v*

Figure 3.31. *D* for the parameter combination of *b* = 0.0001 and <sup>2</sup>

Figure 3.32 *D* for the parameter combination of *b* = 0.001 and <sup>2</sup>

mid range value *b*, 0.01 m, and larger *b*, 0.3, respectively.

relationship of

*D*, 

relationship of

dependence relationship between

Figure 3.33. *D* for the combination of b = 0.3 and σ2 = 0.0001.

Figure 3.34. Field measured values of longitudinal dispersivity as a function of the scale of measurement. The largest circles represent the most reliable data. The estimated dispersivity from the SSTM are given by the squares. (Source: Gelhar (1986).)

Figures 3.35 – 3.37 demonstrate that corresponding *D* values obtained for SSTM parameters do not agree with the relationship of *<sup>L</sup> = 0.1 L.* However, for mid and larger ranges of *b*, estimated *D* s are in reasonable agreement with the most of the reliable field measurements

observed by Gelhar (1986). Data in Figure 3.35 represent a solute transport system with very low stochasticity, which may not be realistic in real world aquifers. Figures 3.38 and 3.39

A Stochastic Model for Hydrodynamic Dispersion 115

were obtained from many sites around the world. The estimated data that agrees with field

To simplify the computational burden, we make use of the fact that the expected value of

 <sup>2</sup> <sup>2</sup> 0 1 0 1

0

To illustrate this approximation, the number of eigenvalues required for each *a* (*M*) is determined by calculating a number of roots that are sufficient for the *M*th eigenvalue to reach a value , where was set to be 5% of the value of the first eigenvalue. The relationship between *a* and the number of roots required is approximately 150 \* *a*, for *a* = 1

{150,750,1500,7500,15000,75000,150000,300000,450000}, then for a =

425.0663, 850.1131, 1710.6345, and 2684.8313}. This shows the extent of the computational

using equation (3.10.1) for different values of *a* with *b=*0.01*, <sup>t</sup>* =0.001, 0.01 *<sup>x</sup>* , 2

 =1.0, and 150\* *a* number of eigen values. For each *a* , the computation was done using a set of routines written in the Python™ language with the *Numeric* extension for fast array

expensive operation with the computing times increasing exponentially with *a*. As an example, for *a* = 1000 meters, it takes approximately 200 minutes, with *a* = 2000 taking approximately 760 minutes on a 1.8 GHZ computer. For *a*

{0.8653, 0.8647,0.8630,0.8615,0.8611,0.8610,0.8610,0.8655, 0.8983}, respectively, and we can expect that for an infinite number of eigen functions it will approach 1.0. Therefore, the

={1,5,10,50,100,500,1000,2000,3000}, the corresponding values for *E Var d t* [ ( ( ))]

it can be shown that approximately, <sup>2</sup>

Therefore, <sup>2</sup> *E Var d t t* [ ( ( ))] 

*a*

*a M a M <sup>i</sup> i i <sup>i</sup> <sup>i</sup>*

*E Var d t t*

( ) ( ( ))

*a a*

( ) is zero and the variance is given by equation (3.4.5). Instead of calculating

*t dx f x dx*

() 1

 

1 *M <sup>i</sup> <sup>M</sup> <sup>i</sup> Limit a* 

= {0.8762, 4.2920, 8.5443, 42.5518, 85, 04312,

( ) . We have computed ( *E Var d t* [ ( ( ))]

/ *t* ) for larger values of *a* is a computationally

over the length *a,* 

*<sup>f</sup> x dx <sup>i</sup>* , and, (3.10.2)

. (3.10.3)

. (3.10.4)

=1.0 and b=0.01, *M* =

/ *t* )

/ *t* are

*<sup>L</sup>* values of Figure 3.36

(3.10.1)

may be a better representation of an actual aquifer. Note that the

( ) for each *x* within *a,* we take the mean of *Var d*[ ]

to 10 000 metres. For example, when <sup>2</sup>

{1,5,10,50,100,500,1000,2000,3000}. 1

operations. Computing ( *E Var d t* [ ( ( ))]

problem without the approximation for *d t*

[ ( ( ))]

mean *d t* 

*d t* 

measurements may be obtained from a variety of soil and heterogeneity.

But by orthogonality, <sup>2</sup>

*M i i* 

Figure 3.35. Estimated *D* values for *b* = 0.0001 m for range of σ2 (0.0001 – 0.3), *D= 0.1 L v,*  and reliable field measurements observed by Gelhar (1986) for different flow lengths.

Figure 3.36. Estimated *D* values for *b* = 0.01 m for range of σ2 (0.0001 – 0.3), *D= 0.1Lv,* and reliable field measurements observed by Gelhar (1986) for different flow lengths.

Figure 3.37 Estimated *D* values for *b* = 0.1 m for range of <sup>2</sup> (0.0001 – 0.3), *D= 0.1Lv,* and reliable field measurements observed by Gelhar (1986) for different flow lengths.

114 in Porous Media - An Approach Based on Stochastic Calculus

Figure 3.35. Estimated *D* values for *b* = 0.0001 m for range of σ2 (0.0001 – 0.3), *D= 0.1 L v,* 

Figure 3.36. Estimated *D* values for *b* = 0.01 m for range of σ2 (0.0001 – 0.3), *D= 0.1Lv,* and

(0.0001 – 0.3), *D= 0.1Lv,* and

reliable field measurements observed by Gelhar (1986) for different flow lengths.

reliable field measurements observed by Gelhar (1986) for different flow lengths.

Figure 3.37 Estimated *D* values for *b* = 0.1 m for range of <sup>2</sup>

and reliable field measurements observed by Gelhar (1986) for different flow lengths.

observed by Gelhar (1986). Data in Figure 3.35 represent a solute transport system with very low stochasticity, which may not be realistic in real world aquifers. Figures 3.38 and 3.39

may be a better representation of an actual aquifer. Note that the *<sup>L</sup>* values of Figure 3.36 were obtained from many sites around the world. The estimated data that agrees with field measurements may be obtained from a variety of soil and heterogeneity.

To simplify the computational burden, we make use of the fact that the expected value of mean *d t* ( ) is zero and the variance is given by equation (3.4.5). Instead of calculating *d t* ( ) for each *x* within *a,* we take the mean of *Var d*[ ] over the length *a,* 

$$E[Var(d\mathcal{J}(t))] = \frac{\int\_0^t \left(\sum\_{i=1}^M (\eta\_i^2)\,\Delta t\right) d\mathbf{x}}{a} = \left(\frac{\int\_0^t \sum\_{i=1}^M (\lambda\_i f\_i^2(\mathbf{x})) \,d\mathbf{x}}{a}\right) \Delta t \tag{3.10.1}$$

$$\text{But by orthogonality, }\begin{array}{c} \stackrel{a}{\underset{f}{\underset{f}{\underset{1}}}}(x)dx = 1 \text{, and,} \\ 0 \end{array} \tag{3.10.2}$$

it can be shown that approximately, <sup>2</sup> 1 *M <sup>i</sup> <sup>M</sup> <sup>i</sup> Limit a* . (3.10.3)

$$\text{Therefore,} \quad \text{E}[\text{Var}(d\beta(t))] = \sigma^2 \,\Delta t \,\,. \tag{3.10.4}$$

To illustrate this approximation, the number of eigenvalues required for each *a* (*M*) is determined by calculating a number of roots that are sufficient for the *M*th eigenvalue to reach a value , where was set to be 5% of the value of the first eigenvalue. The relationship between *a* and the number of roots required is approximately 150 \* *a*, for *a* = 1 to 10 000 metres. For example, when <sup>2</sup> =1.0 and b=0.01, *M* = {150,750,1500,7500,15000,75000,150000,300000,450000}, then for a = {1,5,10,50,100,500,1000,2000,3000}. 1 *M i i* = {0.8762, 4.2920, 8.5443, 42.5518, 85, 04312, 425.0663, 850.1131, 1710.6345, and 2684.8313}. This shows the extent of the computational problem without the approximation for *d t* ( ) . We have computed ( *E Var d t* [ ( ( ))] / *t* ) using equation (3.10.1) for different values of *a* with *b=*0.01*, <sup>t</sup>* =0.001, 0.01 *<sup>x</sup>* , 2 =1.0, and 150\* *a* number of eigen values. For each *a* , the computation was done using a set of routines written in the Python™ language with the *Numeric* extension for fast array operations. Computing ( *E Var d t* [ ( ( ))] / *t* ) for larger values of *a* is a computationally expensive operation with the computing times increasing exponentially with *a*. As an example, for *a* = 1000 meters, it takes approximately 200 minutes, with *a* = 2000 taking approximately 760 minutes on a 1.8 GHZ computer. For *a* ={1,5,10,50,100,500,1000,2000,3000}, the corresponding values for *E Var d t* [ ( ( ))] / *t* are {0.8653, 0.8647,0.8630,0.8615,0.8611,0.8610,0.8610,0.8655, 0.8983}, respectively, and we can expect that for an infinite number of eigen functions it will approach 1.0. Therefore, the

**A Generalized Mathematical** 

In the previous chapter we derived a stochastic solute transport model (equation (3.2.14)); we developed the methods to estimate its parameters, and investigated its behaviour numerically. We see some promise to characterise the solute dispersion at different flow lengths, and there are some indications that equation (3.2.14) produce the behaviours that would be interpreted as capturing the scale-dependency of dispersivity. However, there are weaknesses in the model as evident from Chapter 3. These weaknesses, which are discussed in the next section, are stemming from the very assumptions we made in the development of the model. One could argue that by relaxing the Fickian assumptions, we are actually complicating the problem quite unnecessarily. But as we see in Chapter 3 and in this chapter, we develop a new mathematical and computational machinery at a more

We see that equation (3.2.14) is based on assuming a covariance kernel for the velocity fluctuations, and the solution is dependent on solving an integral equation (see equation (3.3.11)). In Chapter 3, the integral equation is solved analytically for the covariance kernel given by equation (3.3.10) to obtain the eigen values and eigen functions, but analytical solutions of integral equations can not be easily derived for any arbitrary covariance kernel. This limits the flexibility of the SSTM in employing a suitable covariance kernel independent of the ability to solve relevant integral equations. Further, we need to solve the SSTM in a much more computationally efficient manner, and estimating dispersivity by always relating to the deterministic advection-dispersion equation is not quite satisfactory. Therefore, we seek to develop a more general form of equation (3.2.14) in this chapter.

*dC S V x t C x t dt S C x t d t* , , ,

where

We use the same notations and symbols as in Chapter 3. In equation (4.2.2), *d t*

 

*d t* 

associated independent Wiener process increment ( *db t <sup>j</sup>* ).

1 *m m jj j j*

*f db t*

. (4.2.2)

*jj j f db t* ), and for each eigen function, *<sup>j</sup> f* , there is an

*<sup>m</sup>* , (4.2.1)

*<sup>m</sup>* is

fundamental level for the hydrodynamic dispersion in saturated porous media.

**4.2 The Development of the Generalized Model**  We restate equation (3.2.14) in the differential form:

calculated by summing *m* terms of (

**4.1 Introduction** 

**Model in One-Dimension** 

spatially averaged *d t* ( ) is a Gaussian random variable with zero-mean and <sup>2</sup> *t* variance, which can be readly incorporated into the numerical solution scheme for SSTM.

Because of the high computational times involved, instead of using the SIM procedure, we estimate the dispersivity for a computational experiment by limiting ourselves to the SSTM parameters obtained for experimental sand aquifer (i.e., <sup>2</sup> =0.01; *b=*0.01 m ) in the following way:

(a) Set the initial and boundary conditions as,

$$\begin{aligned} \text{C}(\infty,0) &= 0, \quad \text{x} \ge 0\\ \text{C}(0,t) &= 1.0, \; t \ge 0 \quad \text{and} \\ \text{C}(\infty,t) &= 0, \quad t \ge 0 \end{aligned}$$

(b) Solve equation (21) assuming mean velocity to be 1.0 m/day,

(c) Use a realization of *C*(*x,t*) at x=*a* to estimate the dispersivity, *<sup>L</sup>* , using equation (3.10.5) (Fetter, 1999) given below using nonlinear regression: ( Mathematica® was used for this purpose.)

$$\mathbf{C}(\mathbf{x},t) = 0.5 \left( \text{erfc}\left(\frac{a-t}{2\sqrt{a\_\perp t}}\right) + \exp\left(\frac{a}{a\_\perp}\right) \text{erfc}\left(\frac{a+t}{2\sqrt{a\_\perp t}}\right) \right). \tag{3.10.5}$$

Equation (3.10.5) is the analytical solution for the one-dimensional advection-dispersion equation for the initial and boundary conditions given in (a) above. To get a reliable estimate of *<sup>L</sup>* for given *a*, we need to have *x-*axis length of 1.5 *a* meters and should have *C*(*x,t*) realization upto 2 *a* days. For example, to obtain an estimate of *<sup>L</sup>* for *a* = 3000 meters, we need to run the simulation of a domain of 4500 meters for 6000 days. However, to reduce the computational time, one can use higher *x* and *t* values than ideally suitable in solving SPDEs thereby sacrificing the reliability. ( *x* =0.01 m and *t* =0.00001 days would give very good solutions to SSTM.) However, this approximate procedure is only valid for the velocity covariance kernel used in this chapter. This procedure is not as reliable as the SIM.

We overlay some representative values of dispersivities from the SSTM on a graph of the field measurements obtained by Gelhar (1986) from different experiments in Figure 3.34, which shows that SSTM could model the multi-scale dispersion with a single set of parameters, <sup>2</sup> =0.01; *b=*0.01 m, that would give rise to similar non-linear scale dependency of "deterministic" dispersivity evaluated using the realisations of SSTM. evaluated

In Figure 3.34, the larger (hollow) circles depict the most reliable experimental data whereas the smaller (filled) circles give the data with lesser experimental accuracy. The dotted lines show the bounded region of experimental data. The estimated dispersivity values (filled squares) from SSTM are within the bounded region and follow similar trends to the most reliable experimental data. We estimated the dispersivities from SSTM using only a limited number of computational experiments, and each computational experiment produces a random realisation. Therefore, the estimated dispersivities are stochastic quantities just like the experimental values, and it is reasonable to expect discrepancies within the bounded region.
