**NonFickian Solute Transport**

### **1.1 Models in Solute Transport in Porous Media**

VIII Preface

This work is funded by the Foundation for Science and Technology of New Zealand (FoRST) through Lincoln Ventures Ltd. (LVL), Lincoln University. I am grateful to the Chief Scientist of LVL, my colleague, Dr. Ian Woodhead for overseeing the contractual matters to facilitate the work with a sense of humour. I also acknowledge Dr. John

Finally I am grateful to my wife Professor Sandhya Samarasinghe for understanding the value of this work. Her advice on neural networks helped in the computational methods developed in this work. Sandhya's love and patience remained intact during this piece

Centre for Advanced Computational Solutions (C-fACS)

**Don Kulasiri** Professor

Lincoln University, New Zealand

Bright of Aqualinc Ltd. for managing the project for many years.

of work. To that love and patience, I dedicate this monograph.

This research monograph presents the modelling of solute transport in the saturated porous media using novel stochastic and computational approaches. Our previous book published in the North-Holland series of Applied Mathematics and Mechanics (Kulasiri and Verwoerd, 2002) covers some of our research in an introductory manner; this book can be considered as a sequel to it, but we include most of the basic concepts succinctly here, suitably placed in the main body so that the reader who does not have the access to the previous book is not disadvantaged to follow the material presented.

NonFickian Solute Transport 1

The motivation of this work has been to explain the dispersion in saturated porous media at different scales in underground aquifers (i.e., subsurface groundwater flow), based on the theories in stochastic calculus. Underground aquifers render unique challenges in determining the nature of solute dispersion within them. Often the structure of porous formations is unknown and they are sometimes notoriously heterogeneous without any recognizable patterns. This element of uncertainty is the over-arching factor which shapes the nature of solute transport in aquifers. Therefore, it is reasonable to review briefly the work already done in that area in the pertinent literature when and where it is necessary. These interludes of previous work should provide us with necessary continuity of thinking in this work.

There is monumental amount of research work done related to the groundwater flow since 1950s. During the last five to six decades major changes to the size and demographics of human populations occurred; as a result, an unprecedented use of the hydrogeological resources of the earth makes contamination of groundwater a scientific, socio-economic and, in many localities, a political issue. What is less obvious in terms of importance is the way a contaminant, a solute, disperses itself within the geological formations of the aquifers. Experimentation with real aquifers is expensive; hence the need for mathematical and computational models of solute transport. People have developed many types of models over the years to understand the dynamics of aquifers, such as physical scale models, analogy models and mathematical models (Wang and Anderson, 1982; Anderson and Woessner, 1992; Fetter, 2001; Batu, 2006). All these types of models serve different purposes.

Physical scale models are helpful to understand the salient features of groundwater flow and measure the variables such as solute concentrations at different locations of an artificial aquifer. A good example of this type of model is the two artificial aquifers at Lincoln University, New Zealand, a brief description of which appears in the monograph by Kulasiri and Verwoerd (2002). Apart from understanding the physical and chemical processes that occur in the aquifers, the measured variables can be used to partially validate the mathematical models. Inadequacy of these physical models is that their flow lengths are

There are three distinct processes that contribute to the transport of solute in groundwater: convection, dispersion, and diffusion. Convection or advective transport refers to the dissolved solid transport due to the average bulk flow of the ground water. The quantity of solute being transported, in advection, depends on the concentration and quantity of ground water flowing. Different pore sizes, different flow lengths and friction in pores cause ground water to move at rates that are both greater and lesser than the average linear velocity. Due to these multitude of non-uniform non-parallel flow paths within which water moves at different velocities, mixing occurs in flowing ground water. The mixing that occurs in parallel to the flow direction is called hydrodynamic longitudinal dispersion; the word "hydrodynamic" signifies the momentum transfers among the fluid molecules. Likewise, the hydrodynamic transverse dispersion is the mixing that occurs in directions normal to the direction of flow. Diffusion refers to the spreading of the pollutant due to its concentration gradients, i.e., a solute in water will move from an area of greater concentration towards an area where it is less concentrated. Diffusion, unlike dispersion will occur even when the fluid has a zero mean velocity. Due to the tortuosity of the pores, the rate of diffusion in an aquifer is lower than the rate in water alone, and is usually considered negligible in aquifer flow when compared to convection and dispersion (Fetter, 2001). (Tortuosity is a measure of the effect of the shape of the flow path followed by water molecules in a porous media). The latter two processes are often lumped under the term hydrodynamic dispersion. Each of the three transport processes can dominate under different circumstances, depending on the

NonFickian Solute Transport 3

The combination of these three processes can be expressed by the advection – dispersion equation (Bear, 1979; Fetter, 1999; Anderson and Woessner, 1992; Spitz and Moreno, 1996; Fetter, 2001). Other possible phenomenon that can present in solute transport such as adsorption and the occurrence of short circuits are assumed negligible in this case. Derivation of the advection-dispersion equation is given by Ogata (1970), Bear (1972), and Freeze and Cherry (1979). Solutions of the advection-dispersion equation are generally based on a few working assumptions such as: the porous medium is homogeneous, isotropic and saturated with fluid, and flow conditions are such that Darcy's law is valid (Bear, 1972; Fetter, 1999). The two-dimensional deterministic advection – dispersion

> 2 2 *L* 2 *T* 2 *x C C CC DDv tx y x*

where *C* is the solute concentration (*M/L3*), *t* is time (*T*), *DL* is the hydrodynamic dispersion coefficient parallel to the principal direction of flow (longitudinal) (*L2/T*), *DT* is the hydrodynamic dispersion coefficient perpendicular to the principal direction of flow (transverse) (*L2/T*), and *<sup>x</sup> v* is the average linear velocity (*L/T*) in the direction of flow.

It is usually assumed that the hydrodynamic dispersion coefficients will have Gaussian distributions that is described by the mean and variance; therefore we express them as

, (1.1.1)

water an

rate of fluid flow and the nature of the medium (Bear, 1972).

equation can be written as (Fetter, 1999),

follows:

fixed (in the case of Lincoln aquifers, flow length is 10 m), and the porous structure cannot be changed, and therefore a study involving multi-scale general behaviour of solute transport in saturated porous media may not be feasible. Analog models, as the name suggests, are used to study analogues of real aquifers by using electrical flow through conductors. While worthwhile insights can be obtained from these models, the development of and experimentation on these models can be expensive, in addition to being cumbersome and time consuming.These factors may have contributed to the popular use of mathematical and computational models in recent decades (Bear, 1979; Spitz and Moreno, 1996; Fetter, 2001).

A mathematical model consists of a set of differential equations that describe the governing principles of the physical processes of groundwater flow and mass transport of solutes. These time-dependent models have been solved analytically as well as numerically (Wang and Anderson, 1982; Anderson and Woessner, 1992; Fetter, 2001). Analytical solutions are often based on simpler formulations of the problems, for example, using the assumptions on homogeneity and isotropy of the medium; however, they are rich in providing the insights into the untested regimes of behaviour. They also reduce the complexity of the problem (Spitz and Moreno, 1996), and in practice, for example, the analytical solutions are commonly used in the parameter estimation problems using the pumping tests (Kruseman and Ridder, 1970). Analytical solutions also find wide applications in describing the onedimensional and two-dimensional steady state flows in homogeneous flow systems (Walton, 1979). However, in transport problems, the solutions of mathematical models are often intractable; despite this difficulty there are number of models in the literature that could be useful in many situations: Ogata and Banks' (1961) model on one-dimensional longitudinal transport is such a model. A one-dimensional solution for transverse spreading (Harleman and Rumer (1963)) and other related solutions are quite useful (see Bear (1972); Freeze and Cherry (1979)).

Numerical models are widely used when there are complex boundary conditions or where the coefficients are nonlinear within the domain of the model or both situations occur simultaneously (Zheng and Bennett, 1995). Rapid developments in digital computers enable the solutions of complex groundwater problems with numerical models to be efficient and faster. Since numerical models provide the most versatile approach to hydrology problems, they have outclassed all other types of models in many ways; especially in the scale of the problem and heterogeneity. The well-earned popularity of numerical models, however, may lead to over-rating their potential because groundwater systems are complicated beyond our capability to evaluate them in detail. Therefore, a modeller should pay great attention to the implications of simplifying assumptions, which may otherwise become a misrepresentation of the real system (Spitz and Moreno, 1996).

Having discussed the context within which this work is done, we now focus on the core problem, the solute transport in porous media. We are only concerned with the porous media saturated with water, and it is reasonable to assume that the density of the solute in water is similar to that of water. Further we assume that the solute is chemically inert with respect to the porous material. While these can be included in the mathematical developments, they tend to mask the key problem that is being addressed.

2 in Porous Media - An Approach Based on Stochastic Calculus

fixed (in the case of Lincoln aquifers, flow length is 10 m), and the porous structure cannot be changed, and therefore a study involving multi-scale general behaviour of solute transport in saturated porous media may not be feasible. Analog models, as the name suggests, are used to study analogues of real aquifers by using electrical flow through conductors. While worthwhile insights can be obtained from these models, the development of and experimentation on these models can be expensive, in addition to being cumbersome and time consuming.These factors may have contributed to the popular use of mathematical and computational models in recent decades (Bear, 1979; Spitz and Moreno, 1996; Fetter,

A mathematical model consists of a set of differential equations that describe the governing principles of the physical processes of groundwater flow and mass transport of solutes. These time-dependent models have been solved analytically as well as numerically (Wang and Anderson, 1982; Anderson and Woessner, 1992; Fetter, 2001). Analytical solutions are often based on simpler formulations of the problems, for example, using the assumptions on homogeneity and isotropy of the medium; however, they are rich in providing the insights into the untested regimes of behaviour. They also reduce the complexity of the problem (Spitz and Moreno, 1996), and in practice, for example, the analytical solutions are commonly used in the parameter estimation problems using the pumping tests (Kruseman and Ridder, 1970). Analytical solutions also find wide applications in describing the onedimensional and two-dimensional steady state flows in homogeneous flow systems (Walton, 1979). However, in transport problems, the solutions of mathematical models are often intractable; despite this difficulty there are number of models in the literature that could be useful in many situations: Ogata and Banks' (1961) model on one-dimensional longitudinal transport is such a model. A one-dimensional solution for transverse spreading (Harleman and Rumer (1963)) and other related solutions are quite useful (see Bear (1972);

Numerical models are widely used when there are complex boundary conditions or where the coefficients are nonlinear within the domain of the model or both situations occur simultaneously (Zheng and Bennett, 1995). Rapid developments in digital computers enable the solutions of complex groundwater problems with numerical models to be efficient and faster. Since numerical models provide the most versatile approach to hydrology problems, they have outclassed all other types of models in many ways; especially in the scale of the problem and heterogeneity. The well-earned popularity of numerical models, however, may lead to over-rating their potential because groundwater systems are complicated beyond our capability to evaluate them in detail. Therefore, a modeller should pay great attention to the implications of simplifying assumptions, which may otherwise become a

Having discussed the context within which this work is done, we now focus on the core problem, the solute transport in porous media. We are only concerned with the porous media saturated with water, and it is reasonable to assume that the density of the solute in water is similar to that of water. Further we assume that the solute is chemically inert with respect to the porous material. While these can be included in the mathematical

misrepresentation of the real system (Spitz and Moreno, 1996).

developments, they tend to mask the key problem that is being addressed.

2001).

Freeze and Cherry (1979)).

There are three distinct processes that contribute to the transport of solute in groundwater: convection, dispersion, and diffusion. Convection or advective transport refers to the dissolved solid transport due to the average bulk flow of the ground water. The quantity of solute being transported, in advection, depends on the concentration and quantity of ground water flowing. Different pore sizes, different flow lengths and friction in pores cause ground water to move at rates that are both greater and lesser than the average linear velocity. Due to these multitude of non-uniform non-parallel flow paths within which water moves at different velocities, mixing occurs in flowing ground water. The mixing that occurs in parallel to the flow direction is called hydrodynamic longitudinal dispersion; the word "hydrodynamic" signifies the momentum transfers among the fluid molecules. Likewise, the hydrodynamic transverse dispersion is the mixing that occurs in directions normal to the direction of flow. Diffusion refers to the spreading of the pollutant due to its concentration gradients, i.e., a solute in water will move from an area of greater concentration towards an area where it is less concentrated. Diffusion, unlike dispersion will occur even when the fluid has a zero mean velocity. Due to the tortuosity of the pores, the rate of diffusion in an aquifer is lower than the rate in water alone, and is usually considered negligible in aquifer flow when compared to convection and dispersion (Fetter, 2001). (Tortuosity is a measure of the effect of the shape of the flow path followed by water molecules in a porous media). The latter two processes are often lumped under the term hydrodynamic dispersion. Each of the three transport processes can dominate under different circumstances, depending on the rate of fluid flow and the nature of the medium (Bear, 1972).

The combination of these three processes can be expressed by the advection – dispersion equation (Bear, 1979; Fetter, 1999; Anderson and Woessner, 1992; Spitz and Moreno, 1996; Fetter, 2001). Other possible phenomenon that can present in solute transport such as adsorption and the occurrence of short circuits are assumed negligible in this case. Derivation of the advection-dispersion equation is given by Ogata (1970), Bear (1972), and Freeze and Cherry (1979). Solutions of the advection-dispersion equation are generally based on a few working assumptions such as: the porous medium is homogeneous, isotropic and saturated with fluid, and flow conditions are such that Darcy's law is valid (Bear, 1972; Fetter, 1999). The two-dimensional deterministic advection – dispersion equation can be written as (Fetter, 1999),

$$\frac{\partial \mathbf{C}}{\partial t} = 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) - v\_x \left( \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right)\_{\prime} \tag{1.1.1}$$

where *C* is the solute concentration (*M/L3*), *t* is time (*T*), *DL* is the hydrodynamic dispersion coefficient parallel to the principal direction of flow (longitudinal) (*L2/T*), *DT* is the hydrodynamic dispersion coefficient perpendicular to the principal direction of flow (transverse) (*L2/T*), and *<sup>x</sup> v* is the average linear velocity (*L/T*) in the direction of flow.

It is usually assumed that the hydrodynamic dispersion coefficients will have Gaussian distributions that is described by the mean and variance; therefore we express them as follows:

Longitudinal hydrodynamic dispersion coefficient,

$$D\_{\mu} = \frac{\sigma\_{\ll}^{2}}{2t} \quad \text{and} \tag{1.1.2}$$

matrix through which solute flows. To understand how equation (1.1.1), which is a working model of dispersion, came about, it is important to understand its derivation better and the

NonFickian Solute Transport 5

There is much work done in this area using the deterministic description of mass conservation. In the derivation of advection–dispersion equation, also known as continuum transport model, (see Rashidi et al. (1999)), one takes the velocity fluctuations around the mean velocity to calculate the solute flux at a given point using the averaging theorems. The solute flux can be divided into two parts: mean advective flux which stems from the mean velocity and the mean concentration at a given point in space; and the mean dispersive flux which results from the averaging of the product of the fluctuating velocity component and the fluctuating concentration component. These fluctuations are at the scale of the particle sizes, and these fluctuations give rise to hydrodynamic dispersion over time along the porous medium in which solute is dispersed. If we track a single particle with time along one dimensional direction, the velocity fluctuation of the solute particle along that direction is a function of the pressure differential across the medium and the geometrical shapes of the particles, consequently the shapes of the pore spaces. These factors get themselves incorporated into the advection-dispersion equation through the assumptions which are

To understand where the dispersion terms originate, it is worthwhile to review briefly the continuum model for the advection and dispersion in a porous medium (see Rashidi et al. (1999)). The mass conservation has been applied to a neutral solute assuming that the porosity of the region in which the mass is conserved does not change abruptly, i.e., changes in porosity would be continuous. This essentially means that the fluctuations which exist at the pore scale get smoothened out at the scale in which the continuum model is derived. However, the pore scale fluctuations give rise to hydrodynamic dispersion in the first place, and we can expect that the continuum model is more appropriate for homogeneous media. Consider the one dimensional problem of advection and dispersion in a porous medium without transverse dispersion. Assuming that the porous matrix is saturated with water of density, *ρ*, the local flow velocity with respect to pore structure and the local concentration are denoted by *v*(*x,t*) and *c*(*x,t*) at a given point *x*, respectively. These variables are interpreted as intrinsic volume average quantities over a representative elementary volume (Thompson and Gray, 1986). Because the solute flux is transient, conservation of solute mass is expressed by the time-dependent equation of continuity, a form of which is given below:

 

*<sup>c</sup> v c <sup>J</sup> <sup>c</sup> <sup>D</sup> t x xx x*

 

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

 

*A B C*

where *vx* is the mean velocity in the *x*- direction, *c* is the intrinsic volume average concentration, *φ* is the porosity, *Jx* and *τx* are the macroscopic dispersive flux and diffusive tortuosity, respectively. They are approximated by using constitutive relationships for the

*m x*

  , (1.2.1)

 

 

assumptions underpinning the development of the model.

**1.2 Deterministic Models of Dispersion** 

similar to the Fick's law in physics.

medium.

 

 

> 

transverse hydrodynamic dispersion coefficient,

$$D\_r = \frac{\sigma\_r^2}{2t} \tag{1.1.3}$$

where <sup>2</sup> *<sup>L</sup>* is the variance of the longitudinal spreading of the plume, and <sup>2</sup> *<sup>T</sup>* is the variance of the transverse spreading of the plume.

The dispersion coefficients can be thought of having two components: the first measure would reflect the hydrodynamic effects and the other component would indicate the molecular diffusion. For example, for the longitudinal dispersion coefficient,

$$\mathbf{D}\_{\rm L} = \mathbf{a}\_{\rm L} \mathbf{v}\_{\rm L} + \mathbf{D}^{\rm L} \tag{1.1.4}$$

where *<sup>L</sup>* is the longitudinal dynamic dispersivity, *<sup>L</sup> v* is the average linear velocity in longitudinal direction, and \* *D* is the effective diffusion coefficient.

A similar equation can be written for the transverse dispersion as well. Equation (1.1.4) introduces a measure of dispersivity, *<sup>L</sup>* , which has the length dimension, and it can be considered as the average length a solute disperses when mean velocity of solute is unity. Usually in aquifers, diffusion can be neglected compared to the convective flow. Therefore, if velocity is written as a derivative of travel length with respect to time, the simplified version of equation (1.1.4) ( *D v L Li* ) shows a similar relationship as Fick's law in physics.

(Fick's first law expresses that the mass of fluid diffusing is proportional to the concentration gradient. In one dimension, Fick's first law can be expressed as:

$$F = -D\_d \frac{dC}{d\mathbf{x}}\,'$$

where *F* is the mass flux of solute per unit area per unit time (*M/ L2/T*), *Dd* is the diffusion coefficient (*L2/T*), *C* is the solute concentration (*M/L3*), and *dC dx* is the concentration gradient (*M/L3*/*L*).

Fick's second law gives, in one dimension,

$$\frac{dC}{dt} = D\_d \frac{\partial^2 C}{\partial \mathbf{x}^2} \cdot \mathbf{\hat{\beta}}$$

In general, dispersivity is considered as a property of a porous medium. Within equation (1.1.1) hydrodynamic dispersion coefficients represent the average dispersion for each direction for the entire domain of flow, and they mainly allude to and help quantifying the fingering effects on dispersing solute due to granular and irregular nature of the porous

Computational Modelling of Multi-Scale Non-Fickian Dispersion

, and (1.1.2)

(1.1.3)

, (1.1.4)

*<sup>L</sup>* , which has the length dimension, and it can be

(

) shows a similar relationship as Fick's law in physics.

*<sup>T</sup>* is the

*dx*

is the

4 in Porous Media - An Approach Based on Stochastic Calculus

2

2 *L*

*t* 

2

2 *T*

*t* 

*L D*

> *T D*

*<sup>L</sup>* is the variance of the longitudinal spreading of the plume, and <sup>2</sup>

The dispersion coefficients can be thought of having two components: the first measure would reflect the hydrodynamic effects and the other component would indicate the

> \* *D vD L LL*

A similar equation can be written for the transverse dispersion as well. Equation (1.1.4)

considered as the average length a solute disperses when mean velocity of solute is unity. Usually in aquifers, diffusion can be neglected compared to the convective flow. Therefore, if velocity is written as a derivative of travel length with respect to time, the simplified

(Fick's first law expresses that the mass of fluid diffusing is proportional to the

*d dC F D dx* ,

where *F* is the mass flux of solute per unit area per unit time (*M/ L2/T*), *Dd* is the

2 *d* 2 *dC C <sup>D</sup> dt <sup>x</sup>* . )

In general, dispersivity is considered as a property of a porous medium. Within equation (1.1.1) hydrodynamic dispersion coefficients represent the average dispersion for each direction for the entire domain of flow, and they mainly allude to and help quantifying the fingering effects on dispersing solute due to granular and irregular nature of the porous

diffusion coefficient (*L2/T*), *C* is the solute concentration (*M/L3*), and *dC*

*<sup>L</sup>* is the longitudinal dynamic dispersivity, *<sup>L</sup> v* is the average linear velocity in

molecular diffusion. For example, for the longitudinal dispersion coefficient,

concentration gradient. In one dimension, Fick's first law can be expressed as:

longitudinal direction, and \* *D* is the effective diffusion coefficient.

Longitudinal hydrodynamic dispersion coefficient,

transverse hydrodynamic dispersion coefficient,

variance of the transverse spreading of the plume.

introduces a measure of dispersivity,

version of equation (1.1.4) ( *D v L Li*

concentration gradient (*M/L3*/*L*).

Fick's second law gives, in one dimension,

where <sup>2</sup> 

where

matrix through which solute flows. To understand how equation (1.1.1), which is a working model of dispersion, came about, it is important to understand its derivation better and the assumptions underpinning the development of the model.

### **1.2 Deterministic Models of Dispersion**

There is much work done in this area using the deterministic description of mass conservation. In the derivation of advection–dispersion equation, also known as continuum transport model, (see Rashidi et al. (1999)), one takes the velocity fluctuations around the mean velocity to calculate the solute flux at a given point using the averaging theorems. The solute flux can be divided into two parts: mean advective flux which stems from the mean velocity and the mean concentration at a given point in space; and the mean dispersive flux which results from the averaging of the product of the fluctuating velocity component and the fluctuating concentration component. These fluctuations are at the scale of the particle sizes, and these fluctuations give rise to hydrodynamic dispersion over time along the porous medium in which solute is dispersed. If we track a single particle with time along one dimensional direction, the velocity fluctuation of the solute particle along that direction is a function of the pressure differential across the medium and the geometrical shapes of the particles, consequently the shapes of the pore spaces. These factors get themselves incorporated into the advection-dispersion equation through the assumptions which are similar to the Fick's law in physics.

To understand where the dispersion terms originate, it is worthwhile to review briefly the continuum model for the advection and dispersion in a porous medium (see Rashidi et al. (1999)). The mass conservation has been applied to a neutral solute assuming that the porosity of the region in which the mass is conserved does not change abruptly, i.e., changes in porosity would be continuous. This essentially means that the fluctuations which exist at the pore scale get smoothened out at the scale in which the continuum model is derived. However, the pore scale fluctuations give rise to hydrodynamic dispersion in the first place, and we can expect that the continuum model is more appropriate for homogeneous media.

Consider the one dimensional problem of advection and dispersion in a porous medium without transverse dispersion. Assuming that the porous matrix is saturated with water of density, *ρ*, the local flow velocity with respect to pore structure and the local concentration are denoted by *v*(*x,t*) and *c*(*x,t*) at a given point *x*, respectively. These variables are interpreted as intrinsic volume average quantities over a representative elementary volume (Thompson and Gray, 1986). Because the solute flux is transient, conservation of solute mass is expressed by the time-dependent equation of continuity, a form of which is given below:

$$
\frac{\partial(\rho \overline{\boldsymbol{\varepsilon}})}{\partial t} + \frac{\partial \widehat{\left(\rho \overline{\boldsymbol{\upsilon}} \overline{\boldsymbol{\varepsilon}} \right)}}{\partial \mathbf{x}} + \frac{\partial \widehat{\left(\rho \overline{\boldsymbol{I}} \overline{\boldsymbol{\varepsilon}} \right)}}{\partial \mathbf{x}} - \frac{\partial}{\partial \mathbf{x}} \widehat{\left(\rho \overline{\boldsymbol{I}} \overline{\boldsymbol{I}} \right)}} - \frac{\partial}{\partial \mathbf{x}} \widehat{\left(\rho \overline{\boldsymbol{D}\_{m}} \frac{\partial \overline{\boldsymbol{\varepsilon}} \overline{\boldsymbol{\varepsilon}}}{\partial \mathbf{x}} + \boldsymbol{\tau}\_{\times} \right)} = 0,\tag{1.2.1}
$$

where *vx* is the mean velocity in the *x*- direction, *c* is the intrinsic volume average concentration, *φ* is the porosity, *Jx* and *τx* are the macroscopic dispersive flux and diffusive tortuosity, respectively. They are approximated by using constitutive relationships for the medium.

By substituting equations (1.2.3), (1.2.4) and (1.2.5) into equation (1.2.1), the working model

NonFickian Solute Transport 7

( ) ( ) (1 ) 0. *v c <sup>c</sup> <sup>x</sup> <sup>c</sup> DD G t xx m x*

() ( ) (1 ) 0. *<sup>x</sup>*

*c vc <sup>c</sup> DD G t xx x*

many cases, *D>>Dm* , therefore, *DH ≈ D.* We simply refer to *D* as the dispersion coefficient. For a flow with a constant mean velocity through a porous matrix having a constant

In his pioneering work, Taylor (1953) used an equation analogous to equation (1.2.6) to study the dispersion of a soluble substance in a slow moving fluid in a small diameter tube, and he primarily focused on modelling the molecular diffusion coefficient using concentration profiles along a tube for large time. Following that work, Gill and Sankarasubramanian (1970) developed an exact solution for the local concentration for the fully developed laminar flow in a tube for all time. Their work shows that the timedependent dimensionless dispersion coefficient approaches an asymptotic value for larger time proving that Taylor's analysis is adequate for steady-state diffusion through tubes. Even though the above analyses are primarily concerned with the diffusive flow in smalldiameter tubes, as a porous medium can be modelled as a pack of tubes, we could expect similar insights from the advection-dispersion models derived for porous media flow.

The assumptions described by equations (1.2.3) and (1.2.5) above are similar in form to Fick's first law, and therefore, we refer to equations (1.2.3) and (1.2.5) as Fickian assumptions. In particular, equation (1.2.3) defines the dispersivity and dispersion coefficient, which have become so integral to the modelling of dispersion in the literature. As we have briefly explained, dispersivity can be expected to be dependent on the scale of the experiment. This means that, in equations (1.1.1) and (1.2.6), the dispersion coefficient depends on the total length of the flow; mathematically, dispersion coefficient is not only a function of the distance variable *x*, but also a function of the total length. To circumvent the problems associated with solving the mathematical problem, the usual practice is to develop statistical relationships of dispersivity as a function of the total flow length. We discuss some of the relevant research related to ground water flow addressing the scale dependency

flow;

The differences between longitudinal dispersion observed in the field experiments and to the those conducted in the laboratory may be a result of the wide distribution of permeabilities and consequently the velocities found within a real aquifer (Theis 1962, 1963). Fried (1972) presented a few longitudinal dispersivity observations for several sites which were within the range of 0.1 to 0.6 m for the local (aquifer stratum) scale, and within 5 to 11

 

> 

*m*

is called the coefficient of hydrodynamic dispersion. In

(1.2.6)

for solute transport in porous media can be expressed as,

 

> 

 

*x* 

porosity, we see that equation (1.2.6) becomes equation (1.1.1).

The sum (1 ) *<sup>H</sup> <sup>m</sup>*

problem in the next section.

**1.3 A Short Literature Review of Scale Dependency** 

 

*<sup>c</sup> D DD G*

  In equation (1.2.1), the rate of change of the intrinsic volume average concentration is balanced by the spatial gradients of *A0*, *B0*, and *C0* terms, respectively. *A0* represents the average volumetric flux of the solute transported by the average flow of fluid in the *x*direction at a given point in the porous matrix, *x*. However, the fluctuating component of the flux due to the velocity fluctuations around the mean velocity is captured through the term *Jx*(*x,t*) in *B0*,

$$J\_{\mathbf{x}}(\mathbf{x},t) = \overline{\xi\_{\mathbf{x}}^{\prime} \mathbf{c}^{\prime}} \,, \tag{1.2.2}$$

where *ξx* and *c* are the "noise" or perturbation terms of the solute velocity and the concentration about their means, respectively. *C0* denotes the diffusive flux where *Dm* is the fundamental solute diffusivity.

The mean advective flux (*A0*) and the mean dispersive flux (*B0*) can be thought of as representations of the masses of solute carried away by the mean velocity and the fluctuating components of velocity. Further, we do not often know the behaviour of the fluctuating velocity component, and the following assumption, which relates the fluctuating component of the flux to the mean velocity and the spatial gradient of the mean concentration, is used to describe the dispersive flux,

$$J\_{\mathbf{x}}(\mathbf{x},t) \approx -\alpha\_{\mathbf{L}} \overline{v}\_{\mathbf{x}} \frac{\partial \overline{c}}{\partial \mathbf{x}}.\tag{1.2.3}$$

The plausible reasoning behind this assumption is as follows: dispersive flux is proportional to the mean velocity and also proportional to the spatial gradient of the mean concentration. The proportionality constant, *αL*, called the dispersivity, and the subscript *L* indicates the longitudinal direction. Higher the mean velocity, the pore-scale fluctuations are higher but they are subjected to the effects induced by the geometry of the pore structure. This is also true for the dispersive flux component induced by the concentration gradient. Therefore, the dispersivity can be expected to be a material property but its dependency on the spatial concentration gradient makes it vulnerable to the fluctuations in the concentration as so often seen in the experimental situations. The concentration gradients become weaker as the solute plume disperses through a bed of porous medium, and therefore, the mean dispersivity across the bed could be expected to be dependent on the scale of the experiment. This assumption (equation (1.2.3)) therefore, while making mathematical modelling simpler, adds another dimension to the problem: the scale dependency of the dispersivity; and therefore, the scale dependency of the dispersion coefficient, which is obtained by multiplying dispersivity by the mean velocity.

The dispersion coefficient can be expressed as,

$$D = \alpha\_{\perp} \overline{v}\_{\times} \,. \tag{1.2.4}$$

The diffusive tortuosity is typically approximated by a diffusion model of the form,

$$\pi\_{\chi}(\mathbf{x},t) \approx -G \frac{\partial \overline{\mathcal{L}}}{\partial \mathbf{x}} \,' \tag{1.2.5}$$

where *G* is a material coefficient bounded by 0 and 1.

Computational Modelling of Multi-Scale Non-Fickian Dispersion

, (1.2.2)

. (1.2.3)

*L x* . (1.2.4)

, (1.2.5)

6 in Porous Media - An Approach Based on Stochastic Calculus

In equation (1.2.1), the rate of change of the intrinsic volume average concentration is balanced by the spatial gradients of *A0*, *B0*, and *C0* terms, respectively. *A0* represents the average volumetric flux of the solute transported by the average flow of fluid in the *x*direction at a given point in the porous matrix, *x*. However, the fluctuating component of the flux due to the velocity fluctuations around the mean velocity is captured through the

where *ξx* and *c* are the "noise" or perturbation terms of the solute velocity and the concentration about their means, respectively. *C0* denotes the diffusive flux where *Dm* is the

The mean advective flux (*A0*) and the mean dispersive flux (*B0*) can be thought of as representations of the masses of solute carried away by the mean velocity and the fluctuating components of velocity. Further, we do not often know the behaviour of the fluctuating velocity component, and the following assumption, which relates the fluctuating component of the flux to the mean velocity and the spatial gradient of the mean

> ( ,) *<sup>x</sup> L x <sup>c</sup> J xt v*

The plausible reasoning behind this assumption is as follows: dispersive flux is proportional to the mean velocity and also proportional to the spatial gradient of the mean concentration. The proportionality constant, *αL*, called the dispersivity, and the subscript *L* indicates the longitudinal direction. Higher the mean velocity, the pore-scale fluctuations are higher but they are subjected to the effects induced by the geometry of the pore structure. This is also true for the dispersive flux component induced by the concentration gradient. Therefore, the dispersivity can be expected to be a material property but its dependency on the spatial concentration gradient makes it vulnerable to the fluctuations in the concentration as so often seen in the experimental situations. The concentration gradients become weaker as the solute plume disperses through a bed of porous medium, and therefore, the mean dispersivity across the bed could be expected to be dependent on the scale of the experiment. This assumption (equation (1.2.3)) therefore, while making mathematical modelling simpler, adds another dimension to the problem: the scale dependency of the dispersivity; and therefore, the scale dependency of the dispersion coefficient, which is

> *D v*

( ,) *<sup>c</sup> xt G x x*

The diffusive tortuosity is typically approximated by a diffusion model of the form,

*x*

( ,) *<sup>x</sup> <sup>x</sup> J xt c*

term *Jx*(*x,t*) in *B0*,

fundamental solute diffusivity.

concentration, is used to describe the dispersive flux,

obtained by multiplying dispersivity by the mean velocity.

The dispersion coefficient can be expressed as,

where *G* is a material coefficient bounded by 0 and 1.

By substituting equations (1.2.3), (1.2.4) and (1.2.5) into equation (1.2.1), the working model for solute transport in porous media can be expressed as,

$$\frac{\partial \left(\rho \overline{\mathbf{C}}\right)}{\partial \mathbf{t}} + \frac{\partial \left(\rho \overline{\mathbf{v}} \frac{\overline{\mathbf{c}}}{\mathbf{x}} \overline{\mathbf{c}}\right)}{\partial \mathbf{x}} - \frac{\partial}{\partial \mathbf{x}} \left(\rho \mathbf{D} + \rho \mathbf{D}\_{\text{m}} \mathbf{(1-G)} \frac{\overline{\mathbf{c}} \overline{\mathbf{c}}}{\overline{\mathbf{c}} \mathbf{x}}\right) = \mathbf{0}.\tag{1.2.6}$$
 
$$\mathbf{0} = \mathbf{0}, \qquad \mathbf{x} = \mathbf{0}, \qquad \mathbf{x} = \mathbf{0}$$

$$\frac{\partial \left(\boldsymbol{\varrho} \,\, \overline{\boldsymbol{\varepsilon}}\right)}{\partial t} + \frac{\partial \left(\boldsymbol{\varrho} \,\, \overline{\boldsymbol{v}}\_{\times} \, \overline{\boldsymbol{\varepsilon}}\right)}{\partial \mathbf{x}} - \frac{\partial}{\partial \mathbf{x}} \bigg(\boldsymbol{\varrho} \,\, \boldsymbol{D} + \boldsymbol{\varrho} \boldsymbol{D}\_{\boldsymbol{m}} (1 - \boldsymbol{G}) \frac{\boldsymbol{\mathcal{E}} \overline{\boldsymbol{\varepsilon}}}{\partial \mathbf{x}}\bigg) = \mathbf{0}.$$

The sum (1 ) *<sup>H</sup> <sup>m</sup> <sup>c</sup> D DD G x* is called the coefficient of hydrodynamic dispersion. In many cases, *D>>Dm* , therefore, *DH ≈ D.* We simply refer to *D* as the dispersion coefficient. For a flow with a constant mean velocity through a porous matrix having a constant porosity, we see that equation (1.2.6) becomes equation (1.1.1).

In his pioneering work, Taylor (1953) used an equation analogous to equation (1.2.6) to study the dispersion of a soluble substance in a slow moving fluid in a small diameter tube, and he primarily focused on modelling the molecular diffusion coefficient using concentration profiles along a tube for large time. Following that work, Gill and Sankarasubramanian (1970) developed an exact solution for the local concentration for the fully developed laminar flow in a tube for all time. Their work shows that the timedependent dimensionless dispersion coefficient approaches an asymptotic value for larger time proving that Taylor's analysis is adequate for steady-state diffusion through tubes. Even though the above analyses are primarily concerned with the diffusive flow in smalldiameter tubes, as a porous medium can be modelled as a pack of tubes, we could expect similar insights from the advection-dispersion models derived for porous media flow.

The assumptions described by equations (1.2.3) and (1.2.5) above are similar in form to Fick's first law, and therefore, we refer to equations (1.2.3) and (1.2.5) as Fickian assumptions. In particular, equation (1.2.3) defines the dispersivity and dispersion coefficient, which have become so integral to the modelling of dispersion in the literature. As we have briefly explained, dispersivity can be expected to be dependent on the scale of the experiment. This means that, in equations (1.1.1) and (1.2.6), the dispersion coefficient depends on the total length of the flow; mathematically, dispersion coefficient is not only a function of the distance variable *x*, but also a function of the total length. To circumvent the problems associated with solving the mathematical problem, the usual practice is to develop statistical relationships of dispersivity as a function of the total flow length. We discuss some of the relevant research related to ground water flow addressing the scale dependency problem in the next section. bounded dependency

### **1.3 A Short Literature Review of Scale Dependency**

The differences between longitudinal dispersion observed in the field experiments and to the those conducted in the laboratory may be a result of the wide distribution of permeabilities and consequently the velocities found within a real aquifer (Theis 1962, 1963). Fried (1972) presented a few longitudinal dispersivity observations for several sites which were within the range of 0.1 to 0.6 m for the local (aquifer stratum) scale, and within 5 to 11

mean that this deterministic formalism is adequate to capture the reality of solute transport within, often unknown, porous structures. While the deterministic formalisms provide tractable and useful solutions for practical purposes, they may deviate from the reality they represent, in some situations, to unacceptable levels. One could argue that any contracted description of the behaviour of physical ensemble of moving particles must be mechanistic as well as statistical (Keizer, 1987); this may be one of the plausible reasons why there are many stochastic models of groundwater flow. Other plausible reasons are: formations of real world groundwater aquifers are highly heterogeneous, boundaries of the system are multifaceted, inputs are highly erratic, and other subsidiary conditions can be subject to variation as well. Heterogeneous underground formations pose major challenges of developing contracted descriptions of solute transport within them. This was illustrated by injecting a colour liquid into a body of porous rock material with irregular permeability (Øksendal, 1998). These experiments showed that the resulting highly scattered

NonFickian Solute Transport 9

distributions of the liquid were not diffusing according to the deterministic models.

mathematical models of complex systems.

**1.4 Stochastic Models** 

To address the issue of scale dependence of dispersivity and dispersion coefficient fundamentally, it has been argued that a more realistic approach to 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). This well established branch of applied mathematics 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 Weiner processes aid in formulating

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 advection-dispersion 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 it depends on the scale of the experiment.

In respect therefore, uncertainties. hydraulic

The last three decades have seen rapid developments in theoretical research treating groundwater flow and transport problems in a probabilistic framework. The models that are developed under such a theoretical basis are called stochastic models, in which statistical

m for the global (aquifer thickness) scale. These values show the differences in magnitude of the dispersivities. Fried (1975) revisited and redefined these scales in terms of 'mean travelled distance' of the tracer or contaminant as local scale (total flow length between 2 and 4 m), global scale 1 (flow length between 4 and 20 m), global scale 2 (flow length between 20 and 100 m), and regional scale (greater than 100 m; usually several kilometres).

When tested for transverse dispersion, Fried (1972) found no scale effect on the transverse dispersivity and thought that its value could be obtained from the laboratory results. However, Klotz et al. (1980) illustrated from a field tracer test that the width of the tracer plume increased linearly with the travel distance. Oakes and Edworthy (1977) conducted the two-well pulse and the radial injection experiments in a sandstone aquifer and showed that the dispersivity readings for the fully penetrated depth to be 2 to 4 times the values for discrete layers. These results are inconclusive about the lateral dispersivity, and it is very much dependent on the flow length as well as the characteristics of porous matrix subjected to the testing.

Pickens and Grisak (1981), by conducting the laboratory column and field tracer tests, reported that the average longitudinal dispersivity, *<sup>L</sup>* , was 0.035 cm for three laboratory tracer tests with a repacked column of sand when the flow length was 30 cm. For a stratified sand aquifer, by analysing the withdrawal phase concentration histories of a single–well test of an injection withdrawal well, they showed *<sup>L</sup>* were 3 cm and 9 cm for flow lengths of 3.13 m and 4.99 m, respectively. Further, they obtained 50 cm dispersivity in a two-well recirculating withdrawal–injection tracer test with wells located 8 m apart. All these tests were conducted in the same site. Pickens and Grisak (1981) showed that the scale dependency of *L* for the study site has a relationship of *<sup>L</sup> = 0.1 L*, where *L* is the mean travel distance. Lallemand-Barres and Peaudecerf (1978, cited in Fetter, 1999) plotted the field measured *<sup>L</sup>* against the flow length on a log-log graph which strengthened the finding of Pickens and Grisak (1981) and suggested that *<sup>L</sup>* could be estimated to be about 0.1 of the flow length. Gelhar (1986) published a similar representation of the scale of dependency*<sup>L</sup>* using the data from many sites around the world, and according to that study, *<sup>L</sup>* in the range of 1 to 10 m would be reasonable for a site of dimension in the order of 1 km. However, the relationship of *<sup>L</sup>* and the flow length is more complex and not as simple as shown by Pickens and Grisak (1981), and Lallemand-Barres and Peaudecerf (1978, cited in Fetter, 1999). Several other studies on the scale dependency of dispersivity can be found in Peaudecef and Sauty (1978), Sudicky and Cherry (1979), Merritt et al. (1979), Chapman (1979), Lee et al. (1980), Huang et al. (1996b), Scheibe and Yabusaki (1998), Klenk and Grathwohl (2002), and Vanderborght and Vereecken (2002). These empirical relationships influenced the way models developed subsequently. For example, Huang et al. (1996a) developed an analytical solution for solute transport in heterogeneous porous media with scale dependent dispersion. In this model, dispersivity was assumed to increase linearly with flow length until some distance and reaches an asymptotic value. 

Scale dependency of dispersivity shows that the contracted description of the deterministic model has inherent problems that need to be addressed using other forms of contracted descriptions. The Fickian assumptions, for example, help to develop a description which would absorb the fluctuations into a deterministic formalism. But this does not necessarily

8 in Porous Media - An Approach Based on Stochastic Calculus

m for the global (aquifer thickness) scale. These values show the differences in magnitude of the dispersivities. Fried (1975) revisited and redefined these scales in terms of 'mean travelled distance' of the tracer or contaminant as local scale (total flow length between 2 and 4 m), global scale 1 (flow length between 4 and 20 m), global scale 2 (flow length between 20 and 100 m), and regional scale (greater than 100 m; usually several kilometres). When tested for transverse dispersion, Fried (1972) found no scale effect on the transverse dispersivity and thought that its value could be obtained from the laboratory results. However, Klotz et al. (1980) illustrated from a field tracer test that the width of the tracer plume increased linearly with the travel distance. Oakes and Edworthy (1977) conducted the two-well pulse and the radial injection experiments in a sandstone aquifer and showed that the dispersivity readings for the fully penetrated depth to be 2 to 4 times the values for discrete layers. These results are inconclusive about the lateral dispersivity, and it is very much dependent on the flow length as well as the characteristics of porous matrix subjected

Pickens and Grisak (1981), by conducting the laboratory column and field tracer tests,

tracer tests with a repacked column of sand when the flow length was 30 cm. For a stratified sand aquifer, by analysing the withdrawal phase concentration histories of a single–well test

3.13 m and 4.99 m, respectively. Further, they obtained 50 cm dispersivity in a two-well recirculating withdrawal–injection tracer test with wells located 8 m apart. All these tests were conducted in the same site. Pickens and Grisak (1981) showed that the scale

travel distance. Lallemand-Barres and Peaudecerf (1978, cited in Fetter, 1999) plotted the

0.1 of the flow length. Gelhar (1986) published a similar representation of the scale of

simple as shown by Pickens and Grisak (1981), and Lallemand-Barres and Peaudecerf (1978, cited in Fetter, 1999). Several other studies on the scale dependency of dispersivity can be found in Peaudecef and Sauty (1978), Sudicky and Cherry (1979), Merritt et al. (1979), Chapman (1979), Lee et al. (1980), Huang et al. (1996b), Scheibe and Yabusaki (1998), Klenk and Grathwohl (2002), and Vanderborght and Vereecken (2002). These empirical relationships influenced the way models developed subsequently. For example, Huang et al. (1996a) developed an analytical solution for solute transport in heterogeneous porous media with scale dependent dispersion. In this model, dispersivity was assumed to increase

Scale dependency of dispersivity shows that the contracted description of the deterministic model has inherent problems that need to be addressed using other forms of contracted descriptions. The Fickian assumptions, for example, help to develop a description which would absorb the fluctuations into a deterministic formalism. But this does not necessarily

*L* for the study site has a relationship of

linearly with flow length until some distance and reaches an asymptotic value.

*<sup>L</sup>* and the flow length is more complex and not as

*<sup>L</sup>* against the flow length on a log-log graph which strengthened the

*<sup>L</sup>* using the data from many sites around the world, and according to that

*<sup>L</sup>* in the range of 1 to 10 m would be reasonable for a site of dimension in the order

*<sup>L</sup>* , was 0.035 cm for three laboratory

*<sup>L</sup> = 0.1 L*, where *L* is the mean

*<sup>L</sup>* could be estimated to be about

*<sup>L</sup>* were 3 cm and 9 cm for flow lengths of

to the testing.

dependency of

field measured

dependency

study,

reported that the average longitudinal dispersivity,

finding of Pickens and Grisak (1981) and suggested that

of an injection withdrawal well, they showed

of 1 km. However, the relationship of

mean that this deterministic formalism is adequate to capture the reality of solute transport within, often unknown, porous structures. While the deterministic formalisms provide tractable and useful solutions for practical purposes, they may deviate from the reality they represent, in some situations, to unacceptable levels. One could argue that any contracted description of the behaviour of physical ensemble of moving particles must be mechanistic as well as statistical (Keizer, 1987); this may be one of the plausible reasons why there are many stochastic models of groundwater flow. Other plausible reasons are: formations of real world groundwater aquifers are highly heterogeneous, boundaries of the system are multifaceted, inputs are highly erratic, and other subsidiary conditions can be subject to variation as well. Heterogeneous underground formations pose major challenges of developing contracted descriptions of solute transport within them. This was illustrated by injecting a colour liquid into a body of porous rock material with irregular permeability (Øksendal, 1998). These experiments showed that the resulting highly scattered distributions of the liquid were not diffusing according to the deterministic models.

To address the issue of scale dependence of dispersivity and dispersion coefficient fundamentally, it has been argued that a more realistic approach to 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). This well established branch of applied mathematics 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 Weiner processes 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 advection-dispersion 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 it depends on the scale of the experiment. Weiner differential the

### **1.4 Stochastic Models**

The last three decades have seen rapid developments in theoretical research treating groundwater flow and transport problems in a probabilistic framework. The models that are developed under such a theoretical basis are called stochastic models, in which statistical

Welty and Gelhar (1992) studied that the density and fluid viscosity as a function of concentration in heterogeneous aquifers. The spatial and temporal behaviour of the solute front resulting from variable macrodispersion were investigated using analytical results and numerical simulations. The uncertainty in the mass flux for the solute advection in heterogeneous porous media was the research focus of Dagan et al. (1992) and Cvetkovic et al. (1992). Rubin and Dagan (1992) developed a procedure for the characterisation of the head and velocity fields in heterogeneous, statistically anisotropic formations. The velocity field was characterised through a series of spatial covariances as well as the velocity-head and velocity-log conductivity. Other important contributions of stochastic studies in subsurface hydrology can be found in Painter (1996), Yang et al. (1996), Miralles-Wilhelm and Gelhar (1996), Harter and Yeh (1996), Koutsoyiannis (1999), Koutsoyiannis (2000), Zhang and Sun (2000), Foussereau et al. (2000), Leeuwen et al. (2000), Loll and Moldrup (2000), Foussereau et al. (2001) and, Painter and Cvetkovic (2001). In additional to that, Farrell (1999), Farrell (2002a), and Farrell (2002b) made important contributions to the

NonFickian Solute Transport 11

Kulasiri (1997) developed a preliminary stochastic model that describes the solute dispersion in a porous medium saturated with water and considers velocity of the solute as a fundamental stochastic variable. The main feature of this model is it eliminates the use of the hydrodynamic dispersion coefficient, which is subjected to scale effects and based on Fickian assumptions that were discussed in section 1.2. The model drives the mass

In the process of developing the differential equations of any model, we introduce the parameters, which we consider the attributes or properties of the system. In the case of groundwater flow, for example, the parameters such as hydraulic conductivity, transmissivity and porosity are constant within the differential equations, and it is often necessary to assign numerical values to these parameters. There are a few generally accepted direct parameter measurement methods such as the pumping tests, the permeameter tests and grain size analysis (details on these tests can be found in Bear et al. (1968) and Bear (1979)). The values of the parameters obtained from the laboratory experiments and/or the field scale experiments, may not represent the often complex patterns across a large geographical area, hence limiting the validity and credibility of a model. The inaccuracies of the laboratory tests are due to the scale differences of the actual aquifer and the laboratory sample. The heterogeneous porous media is, most of the time, laterally smaller than the longitudinal scale of the flow; in laboratory experiments, due to practical limitations, we deal with proportionally larger lateral dimensions. Hence, the parameter values obtained from the laboratory tests are not directly usable in the models, and generally need to be upscaled using often subjective techniques. This difficulty is recognised as a major impediment to wider use of the groundwater models and their full utilisation (Frind and Pinder, 1973). For this reason, Freeze (1972) stated that the estimation

of

Often we are interested in modelling the quantities such as the depth of water table and solute concentration, which are relevant to environmental decision making, and we measure these variables regularly and the measuring techniques tend to be relatively inexpensive. In

conservation for solute transport based on the theories of stochastic calculus.

of the parameters is the 'Achilles' heel' of groundwater modelling.

stochastic theory in uncertain flows.

**1.5 Inverse Problems of the Models** 

uncertainty of a natural phenomenon, such as solute transport, is expressed within the stochastic governing equations rather than based on deterministic formulations. The probabilistic nature of this outcome is due to the fact that there is a heterogeneous distribution of the underlying aquifer parameters such as hydraulic conductivity and porosity (Freeze, 1975).

The researchers in the field of hydrology have paid more attention to the scale and variability of aquifers over the two past decades. It is apparent that we need to deal with larger scales more than ever to study the groundwater contaminant problems, which are becoming serious environmental concerns. The scale of the aquifer has a direct proportional relationship to the variability. Hence, the potential role of modelling in addressing these challenges is very much dependent on spatial distribution. When working with deterministic models, if we could measure the hydrogeologic parameters at very close spatial intervals (which is prohibitively expensive), the distribution of aquifer properties would have a high degree of detail. Therefore, the solution of the deterministic model would yield results with a high degree of reliability. However, as the knowledge of finegrained hydrogeologic parameters are limited in practice, the stochastic models are used to understand dynamics of aquifers thus recognising the inherent probabilistic nature of the hydrodynamic dispersion.

Early research on stochastic modelling can be categorised in terms of three possible sources of uncertainties: (i) those caused by measurement errors in the input parameters, (ii) those caused by spatial averaging of input parameters, and (iii) those associated with an inherent stochastic description of heterogeneity porous media (Freeze, 1975). Bibby and Sunada (1971) utilised the Monte Carlo numerical simulation model to investigate the effect on the solution of normally distributed measurement errors in initial head, boundary heads, pumping rate, aquifer thickness, hydraulic conductivity, and storage coefficient of transient flow to a well in a confined aquifer. Sagar and Kisiel (1972) conducted an error propagation study to understand the influence of errors in the initial head, transmissibility, and storage coefficient on the drawdown pattern predicted by the Theis equation. We can find that some aspects of the flow in heterogeneous formations have been investigated even in the early 1960s (Warren and Price, 1961; McMillan, 1966). However, concerted efforts began only in 1975, with the pioneering work of Freeze (1975).

Freeze (1975) showed that all soils and geologic formations, even those that are homogeneous, are non-uniform. Therefore, the most realistic representation of a nonuniform porous medium is a stochastic set of macroscopic elements in which the three basic hydrologic parameters (hydraulic conductivity, compressibility and porosity) are assumed to come from the frequency distributions. Gelhar et al. (1979) discussed the stochastic microdispersion in a stratified aquifer, and Gelhar and Axness (1983) addressed the issue of three-dimensional stochastic macro dispersion in aquifers. Dagan (1984) analysed the solute transport in heterogeneous porous media in a stochastic framework, and Gelgar (1986) demonstrated that the necessity of the use of theoretical knowledge of stochastic subsurface hydrology in real world applications. Other major contributions to stochastic groundwater modelling in the decade of 1980 can be found in Dagan (1986), Dagan (1988) and Neuman et al. (1987).

10 in Porous Media - An Approach Based on Stochastic Calculus

uncertainty of a natural phenomenon, such as solute transport, is expressed within the stochastic governing equations rather than based on deterministic formulations. The probabilistic nature of this outcome is due to the fact that there is a heterogeneous distribution of the underlying aquifer parameters such as hydraulic conductivity and

The researchers in the field of hydrology have paid more attention to the scale and variability of aquifers over the two past decades. It is apparent that we need to deal with larger scales more than ever to study the groundwater contaminant problems, which are becoming serious environmental concerns. The scale of the aquifer has a direct proportional relationship to the variability. Hence, the potential role of modelling in addressing these challenges is very much dependent on spatial distribution. When working with deterministic models, if we could measure the hydrogeologic parameters at very close spatial intervals (which is prohibitively expensive), the distribution of aquifer properties would have a high degree of detail. Therefore, the solution of the deterministic model would yield results with a high degree of reliability. However, as the knowledge of finegrained hydrogeologic parameters are limited in practice, the stochastic models are used to understand dynamics of aquifers thus recognising the inherent probabilistic nature of the

Early research on stochastic modelling can be categorised in terms of three possible sources of uncertainties: (i) those caused by measurement errors in the input parameters, (ii) those caused by spatial averaging of input parameters, and (iii) those associated with an inherent stochastic description of heterogeneity porous media (Freeze, 1975). Bibby and Sunada (1971) utilised the Monte Carlo numerical simulation model to investigate the effect on the solution of normally distributed measurement errors in initial head, boundary heads, pumping rate, aquifer thickness, hydraulic conductivity, and storage coefficient of transient flow to a well in a confined aquifer. Sagar and Kisiel (1972) conducted an error propagation study to understand the influence of errors in the initial head, transmissibility, and storage coefficient on the drawdown pattern predicted by the Theis equation. We can find that some aspects of the flow in heterogeneous formations have been investigated even in the early 1960s (Warren and Price, 1961; McMillan, 1966). However, concerted efforts began only in

Freeze (1975) showed that all soils and geologic formations, even those that are homogeneous, are non-uniform. Therefore, the most realistic representation of a nonuniform porous medium is a stochastic set of macroscopic elements in which the three basic hydrologic parameters (hydraulic conductivity, compressibility and porosity) are assumed to come from the frequency distributions. Gelhar et al. (1979) discussed the stochastic microdispersion in a stratified aquifer, and Gelhar and Axness (1983) addressed the issue of three-dimensional stochastic macro dispersion in aquifers. Dagan (1984) analysed the solute transport in heterogeneous porous media in a stochastic framework, and Gelgar (1986) demonstrated that the necessity of the use of theoretical knowledge of stochastic subsurface hydrology in real world applications. Other major contributions to stochastic groundwater modelling in the decade of 1980 can be found in Dagan (1986), Dagan (1988) and Neuman et

porosity (Freeze, 1975).

hydrodynamic dispersion.

al. (1987).

1975, with the pioneering work of Freeze (1975).

Welty and Gelhar (1992) studied that the density and fluid viscosity as a function of concentration in heterogeneous aquifers. The spatial and temporal behaviour of the solute front resulting from variable macrodispersion were investigated using analytical results and numerical simulations. The uncertainty in the mass flux for the solute advection in heterogeneous porous media was the research focus of Dagan et al. (1992) and Cvetkovic et al. (1992). Rubin and Dagan (1992) developed a procedure for the characterisation of the head and velocity fields in heterogeneous, statistically anisotropic formations. The velocity field was characterised through a series of spatial covariances as well as the velocity-head and velocity-log conductivity. Other important contributions of stochastic studies in subsurface hydrology can be found in Painter (1996), Yang et al. (1996), Miralles-Wilhelm and Gelhar (1996), Harter and Yeh (1996), Koutsoyiannis (1999), Koutsoyiannis (2000), Zhang and Sun (2000), Foussereau et al. (2000), Leeuwen et al. (2000), Loll and Moldrup

Kulasiri (1997) developed a preliminary stochastic model that describes the solute dispersion in a porous medium saturated with water and considers velocity of the solute as a fundamental stochastic variable. The main feature of this model is it eliminates the use of the hydrodynamic dispersion coefficient, which is subjected to scale effects and based on Fickian assumptions that were discussed in section 1.2. The model drives the mass conservation for solute transport based on the theories of stochastic calculus.

(2000), Foussereau et al. (2001) and, Painter and Cvetkovic (2001). In additional to that, Farrell (1999), Farrell (2002a), and Farrell (2002b) made important contributions to the

### **1.5 Inverse Problems of the Models**

stochastic theory in uncertain flows.

In the process of developing the differential equations of any model, we introduce the parameters, which we consider the attributes or properties of the system. In the case of groundwater flow, for example, the parameters such as hydraulic conductivity, transmissivity and porosity are constant within the differential equations, and it is often necessary to assign numerical values to these parameters. There are a few generally accepted direct parameter measurement methods such as the pumping tests, the permeameter tests and grain size analysis (details on these tests can be found in Bear et al. (1968) and Bear (1979)). The values of the parameters obtained from the laboratory experiments and/or the field scale experiments, may not represent the often complex patterns across a large geographical area, hence limiting the validity and credibility of a model. The inaccuracies of the laboratory tests are due to the scale differences of the actual aquifer and the laboratory sample. The heterogeneous porous media is, most of the time, laterally smaller than the longitudinal scale of the flow; in laboratory experiments, due to practical limitations, we deal with proportionally larger lateral dimensions. Hence, the parameter values obtained from the laboratory tests are not directly usable in the models, and generally need to be upscaled using often subjective techniques. This difficulty is recognised as a major impediment to wider use of the groundwater models and their full utilisation (Frind and Pinder, 1973). For this reason, Freeze (1972) stated that the estimation of the parameters is the 'Achilles' heel' of groundwater modelling.

Often we are interested in modelling the quantities such as the depth of water table and solute concentration, which are relevant to environmental decision making, and we measure these variables regularly and the measuring techniques tend to be relatively inexpensive. In

should be considered. Even though there are many advantages of this method such as not having to solve an ill-posed inverse problem, this is a rather tedious way of finding parameters when the model is a large one, and subjective judgements of experts may play a

NonFickian Solute Transport 13

The indirect method transfers the inverse problem into an optimisation problem, still using the forward solutions. Steps such as a criterion to decide the better parameters between previous and present values, and also a stopping condition, can be replaced with the computer-aided algorithms (Neuman, 1973; Sun, 1994). One draw back is that this method tends to converge towards local minima rather than global minima of objective functions

The direct method is another optimisation approach to the inverse problem. If the state variables and their spatial and temporal derivatives are known over the entire region, and if the measurement and mass balance errors are negligible, the flow equation becomes a first order partial differential equation in terms of the unknown aquifer parameters. Using numerical methods, the linear partial differential equations can be reduced to a linear system of equations, which can be solved directly for the unknown aquifer parameters, and

The above three methods (trial and error, indirect, and direct) are well established and a large number of advanced techniques have been added. The algorithms to use in these methods can be found in any numerical recipes (for example, Press, 1992). Even though we change the parameter estimation problem for an optimisation problem, the ill-posedness of the inverse problems do still exist. The non-uniqueness of the inverse solution strongly displays itself in the indirect method through the existence of many local minima (Keidser and Rosbjerg, 1991). In the direct method the solution is often unstable (Kuiper, 1986). To overcome the ill-posedness, it is necessary to have supplementary information, or as often referred to as prior information, which is independent of the measurement of state variables. This can be designated parameter values at some specific time and space points or reliable information about the system to limit the admissible range of possible parameters to a narrower range or to assume that an unknown parameter is piecewise constant (Sun, 1994).

The above described optimisation methods are limited to producing the best estimates and can only assess a residual uncertainty. Usually, output is an estimate of the confidence interval of each parameter after a post-calibration sensitivity study. This approach is deemed insufficient to characterise the uncertainty after calibration (Zimmerman et al., 1998). Moreover, these inverse methods are not suitable enough to provide an accurate representation of larger scales. For that reason, the necessity of having statistically sound methods that are capable of producing reasonable distribution of data (parameters) throughout larger regions was identified. As a result, a large number of geostatisticallybased inverse methods have been developed to estimate groundwater parameters (Keidser and Rosbjerg, 1991; Zimmerman et al., 1998). A theoretical underpinning for new geostatistical inverse methods and discussion of geostatistical estimation approach can be found in many publications (Kitanidis and Vomvoris, 1983; Hoeksema and Kitanidis, 1984; Kitanidis, 1985; Carrera, 1988; Gutjahr and Wilson, 1989; Carrera and Glorioso, 1991;

role in determining the parameters (Keidser and Rosbjerg, 1991).

hence the method is named "direct method" (Neuman, 1973; Sun, 1994).

(Yew, 1986; Kuiper, 1986; Keidser and Rosbjerg, 1991).

**1.8 Geostatistical Approach to the Inverse Problem** 

Cressie, 1993; Gomez-Hernandez et al., 1997; Kitanidis, 1997).

addition, we can continuously monitor these decision (output) variables in many situations. Therefore, it is reasonable to assume that these observations of the output variables represent the current status of the system and measurement errors. If the dynamics of the system can be reliably modelled using relevant differential equations, we can expect the parameters estimated, based on the observations, may give us more reliable representative values than those obtained from the laboratory tests and literature. The observations often contain noise from two different sources: experimental errors and noisy system dynamics. Noise in the system dynamics may be due to the factors such as heterogeneity of the media, random nature of inputs (rainfall) and variable boundary conditions. Hence, the question of estimating the parameters from the observations should involve the models that consist of plausible representation of "noises".

### **1.6 Inherent Ill-Posedness**

A well-posed mathematical problem derived from a physical system must satisfy the existence, uniqueness and stability conditions, and if any one of these conditions is not satisfied the problem is ill-posed. But in a physical system itself, these conditions do not necessarily have specific meanings because, regardless of their mathematical descriptions, the physical system would respond to any situation. As different combinations of hydrological factors would produce almost similar results, it may be impossible to determine a unique set of parameters for a given set of mathematical equations. So this lack of uniqueness could only be remedied by searching a large enough parameter space to find a set of parameters that would explain the dynamics of the maximum possible number, if not all, of the state variables satisfactorily. However, these parameter searches guarantee neither uniqueness nor stability in the inverse problems associated with the groundwater problems (Yew, 1986; Carrera, 1987; Sun, 1994; Kuiper, 1986; Ginn and Cushman, 1990; Keidser and Rosbjerg, 1991). The general consensus among groundwater modellers is that the inverse problem may at times result in meaningless solutions (Carrera and Neuman, 1986b). There are even those who argue that the inverse problem is hopelessly ill-posed and as such, intrinsically unsolvable (Carrera and Neuman, 1986b). This view aside, it has been established that a well-posed inverse problem can, in practice, yield an acceptable solution (McLauglin and Townley, 1996). We adopt a positive view point that a mixture of techniques smartly deployed would render us the sets of effective parameters under the regimes of behaviours of the system which we are interested in. Given this stance, we would like to briefly discuss a number of techniques we found useful in the parameter estimation of the models we describe in this monograph. This discussion does not do justice to the methods mentioned and therefore we include the references for further study. We attempt to describe a couple of methods, which we use in this work, inmore detail, but the reader may find the discussion inadequate; therefore, it is essential to follow up the references to understand the techniques thoroughly.

### **1.7 Methods in Parameter Estimation**

The trial and error method is the most simple but laborious for solving the inverse problems to estimate the parameters. In this method, we use a model that represents the aquifer system with some observed data of state variables. It is important, however, to have an expert who is familiar with the system available, i.e., a specific aquifer (Sun, 1994). Candidate parameter values are tried out until satisfactory outputs are obtained. However, if a satisfactory parameter fitting cannot be found, the modification of the model structure

12 in Porous Media - An Approach Based on Stochastic Calculus

addition, we can continuously monitor these decision (output) variables in many situations. Therefore, it is reasonable to assume that these observations of the output variables represent the current status of the system and measurement errors. If the dynamics of the system can be reliably modelled using relevant differential equations, we can expect the parameters estimated, based on the observations, may give us more reliable representative values than those obtained from the laboratory tests and literature. The observations often contain noise from two different sources: experimental errors and noisy system dynamics. Noise in the system dynamics may be due to the factors such as heterogeneity of the media, random nature of inputs (rainfall) and variable boundary conditions. Hence, the question of estimating the parameters from the observations should involve the models that consist of

A well-posed mathematical problem derived from a physical system must satisfy the existence, uniqueness and stability conditions, and if any one of these conditions is not satisfied the problem is ill-posed. But in a physical system itself, these conditions do not necessarily have specific meanings because, regardless of their mathematical descriptions, the physical system would respond to any situation. As different combinations of hydrological factors would produce almost similar results, it may be impossible to determine a unique set of parameters for a given set of mathematical equations. So this lack of uniqueness could only be remedied by searching a large enough parameter space to find a set of parameters that would explain the dynamics of the maximum possible number, if not all, of the state variables satisfactorily. However, these parameter searches guarantee neither uniqueness nor stability in the inverse problems associated with the groundwater problems (Yew, 1986; Carrera, 1987; Sun, 1994; Kuiper, 1986; Ginn and Cushman, 1990; Keidser and Rosbjerg, 1991). The general consensus among groundwater modellers is that the inverse problem may at times result in meaningless solutions (Carrera and Neuman, 1986b). There are even those who argue that the inverse problem is hopelessly ill-posed and as such, intrinsically unsolvable (Carrera and Neuman, 1986b). This view aside, it has been established that a well-posed inverse problem can, in practice, yield an acceptable solution (McLauglin and Townley, 1996). We adopt a positive view point that a mixture of techniques smartly deployed would render us the sets of effective parameters under the regimes of behaviours of the system which we are interested in. Given this stance, we would like to briefly discuss a number of techniques we found useful in the parameter estimation of the models we describe in this monograph. This discussion does not do justice to the methods mentioned and therefore we include the references for further study. We attempt to describe a couple of methods, which we use in this work, inmore detail, but the reader may find the discussion inadequate; therefore, it is

essential to follow up the references to understand the techniques thoroughly.

The trial and error method is the most simple but laborious for solving the inverse problems to estimate the parameters. In this method, we use a model that represents the aquifer system with some observed data of state variables. It is important, however, to have an expert who is familiar with the system available, i.e., a specific aquifer (Sun, 1994). Candidate parameter values are tried out until satisfactory outputs are obtained. However, if a satisfactory parameter fitting cannot be found, the modification of the model structure

plausible representation of "noises".

**1.7 Methods in Parameter Estimation** 

**1.6 Inherent Ill-Posedness** 

should be considered. Even though there are many advantages of this method such as not having to solve an ill-posed inverse problem, this is a rather tedious way of finding parameters when the model is a large one, and subjective judgements of experts may play a role in determining the parameters (Keidser and Rosbjerg, 1991).

The indirect method transfers the inverse problem into an optimisation problem, still using the forward solutions. Steps such as a criterion to decide the better parameters between previous and present values, and also a stopping condition, can be replaced with the computer-aided algorithms (Neuman, 1973; Sun, 1994). One draw back is that this method tends to converge towards local minima rather than global minima of objective functions (Yew, 1986; Kuiper, 1986; Keidser and Rosbjerg, 1991).

The direct method is another optimisation approach to the inverse problem. If the state variables and their spatial and temporal derivatives are known over the entire region, and if the measurement and mass balance errors are negligible, the flow equation becomes a first order partial differential equation in terms of the unknown aquifer parameters. Using numerical methods, the linear partial differential equations can be reduced to a linear system of equations, which can be solved directly for the unknown aquifer parameters, and hence the method is named "direct method" (Neuman, 1973; Sun, 1994).

The above three methods (trial and error, indirect, and direct) are well established and a large number of advanced techniques have been added. The algorithms to use in these methods can be found in any numerical recipes (for example, Press, 1992). Even though we change the parameter estimation problem for an optimisation problem, the ill-posedness of the inverse problems do still exist. The non-uniqueness of the inverse solution strongly displays itself in the indirect method through the existence of many local minima (Keidser and Rosbjerg, 1991). In the direct method the solution is often unstable (Kuiper, 1986). To overcome the ill-posedness, it is necessary to have supplementary information, or as often referred to as prior information, which is independent of the measurement of state variables. This can be designated parameter values at some specific time and space points or reliable information about the system to limit the admissible range of possible parameters to a narrower range or to assume that an unknown parameter is piecewise constant (Sun, 1994). Ginn existence

### **1.8 Geostatistical Approach to the Inverse Problem**

The above described optimisation methods are limited to producing the best estimates and can only assess a residual uncertainty. Usually, output is an estimate of the confidence interval of each parameter after a post-calibration sensitivity study. This approach is deemed insufficient to characterise the uncertainty after calibration (Zimmerman et al., 1998). Moreover, these inverse methods are not suitable enough to provide an accurate representation of larger scales. For that reason, the necessity of having statistically sound methods that are capable of producing reasonable distribution of data (parameters) throughout larger regions was identified. As a result, a large number of geostatisticallybased inverse methods have been developed to estimate groundwater parameters (Keidser and Rosbjerg, 1991; Zimmerman et al., 1998). A theoretical underpinning for new geostatistical inverse methods and discussion of geostatistical estimation approach can be found in many publications (Kitanidis and Vomvoris, 1983; Hoeksema and Kitanidis, 1984; Kitanidis, 1985; Carrera, 1988; Gutjahr and Wilson, 1989; Carrera and Glorioso, 1991; Cressie, 1993; Gomez-Hernandez et al., 1997; Kitanidis, 1997).

The explanation on the transformation of

For 0 *stT* ,

Note that *d t*

The estimate ˆ

function, *l*(

paths,

) *= ln L*(

solution to the equation

parameter

*V t*( ) 

( ) and ( ) *V t* 

> *L*(

Maximising the likelihood function *L*( )

of the groundwater system. The estimate ˆ

given by (Basawa and Prakasa Rao, 1980):

that depends continuously on *t T* 0, and satisfies the following:

NonFickian Solute Transport 15

(1970), and we develop this approach further in the later chapters. A standard Wiener process (often called a Brownian motion) on the interval 0,*T* is a random variable *W t*( )

*W t W s t sN* ( ) ( ) (0,1),

 ) = exp <sup>2</sup> 0 0

*T T*

*θ* can be obtained as the solution to the equation,

0 0

*T T*

Taking log on both sides of equation (1.9.6) we obtain,

*l*(

The parameter is estimated as the solution to the equation

where *N*(0,1) is a random variable generated with zero mean and unit variance.

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

*<sup>W</sup>*(0) 0, (1.9.5)

using the maximum likelihood approach using all the available observations

*θ* of 

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

 <sup>0</sup> *L* 

( ) <sup>0</sup> *<sup>l</sup> θ*

 ) = <sup>2</sup> 0 0

The parameters can be estimated from equation (1.9.10), based on a single sample path. Let us now consider the case when *M* independent sample paths are being observed. The likelihood-function becomes the product of the likelihood functions for *M* individual sample

*<sup>t</sup>* <sup>0</sup>

*S t, V, θ dV - S t, V, θ S t, V, θ dt =*

, , ( )- , , . *<sup>T</sup> <sup>T</sup> <sup>1</sup> S t V dV t S t V dt 2*

); hence, the maximum likelihood estimate can also be obtained as a

*S t V dV t S t V dt*

are defined on the same event space. We estimate the

. (1.9.6)

. (1.9.7)

is equivalent to maximising the log-likelihood

*<sup>θ</sup>* . (1.9.8)

(1.9.9)

*<sup>θ</sup> <sup>θ</sup>* . (1.9.10)

( ) can be found in Jazwinski

maximises the likelihood functions

### **1.9 Parameter Estimation by Stochastic Partial Differential Equations**

The geostatistical approaches mentioned briefly above estimate the distribution of the parameter space based on a few direct measurements and the geological formation of the spatial domain. Therefore, the accuracy of each method is largely dependent on direct measurements that, as mentioned above, are subject to randomness, numerical errors, and the methods of measurements tend to be expensive. Unny (1989) developed an approach based on the theory of stochastic partial differential equations to estimate groundwater parameters of a one-dimensional aquifer fed by rainfall by considering the water table depth as the output variable to identify the current state of the system. The approach inversely estimates the parameters by using stochastic partial differential equations that model the state variables of the system dynamics. Theory of the parameter estimation of stochastic processes can be found in Kutoyants (1984), Lipster and Shirayev (1977), and Basawa and Prakasa Rao (1980). We summarise this approach in some detail as we use this approach to estimate the parameters in our models in this monograph. stochastic standard

Let *V t*( ) denote a stochastic process having many realisations. We define the parameter set of a probability space which is given by a stochastic process ( ) *V t* , based on a set of realisations { ( ) *V t* ; *0tT* }. Let the evolution of the family of stochastic processes { () *V t* ; *t T* ;} be described by a stochastic partial differential equation (SPDE),

$$
\left\| V^{\theta} \right\| \left( t \right) = A \left| V^{\theta} \right| dt + \left\| \left( \mathbf{x}, t \right) \right\| dt \,, \tag{1.9.1}
$$

where *A* is a partial differential operator in space, and ( ,) *x t dt* is the stochastic process to represent a space- and time- correlated noise process.

The stochastic process ( ) *V t* forms infinitely many sub event spaces with increasing times. We can describe the stochastic process *VttT* ( ); ; , and *AV* as a known function of the system,

$$AV^{\theta} = \mathcal{S}\{\mathbf{t}, V, \theta\}. \tag{1.9.2}$$

Therefore, the stochastic process ( ) *V t* can be represented as the solution of the stochastic differential equation (SDE),

$$
\hat{\varepsilon}V''(t) = \mathbf{S}\{t, V, \theta\}dt + \tilde{\xi}\{\mathbf{x}, t\}dt,\tag{1.9.3}
$$

where *S*(.) is a given function.

We can transform the noise process by a Hilbert space valued standard Wiener process increments, ( )*t* . (A Hilbert space is an inner product space that is complete with respect to the norm defined by the inner product; and a separable Hilbert space should contain a complete orthonormal sequence (Young, 1988).) Therefore,

$$
\hat{\sigma}V^{\theta}(t) = \mathcal{S}\left(t, V, \theta\right)dt + d\mathcal{J}(t). \tag{1.9.4}
$$

, (1.9.1)

can be represented as the solution of the stochastic

. (1.9.2)

(1.9.3)

(1.9.4)

( ,) *x t dt* is the stochastic process to

, based on a set of

as a known function

14 in Porous Media - An Approach Based on Stochastic Calculus

The geostatistical approaches mentioned briefly above estimate the distribution of the parameter space based on a few direct measurements and the geological formation of the spatial domain. Therefore, the accuracy of each method is largely dependent on direct measurements that, as mentioned above, are subject to randomness, numerical errors, and the methods of measurements tend to be expensive. Unny (1989) developed an approach based on the theory of stochastic partial differential equations to estimate groundwater parameters of a one-dimensional aquifer fed by rainfall by considering the water table depth as the output variable to identify the current state of the system. The approach inversely estimates the parameters by using stochastic partial differential equations that model the state variables of the system dynamics. Theory of the parameter estimation of stochastic processes can be found in Kutoyants (1984), Lipster and Shirayev (1977), and Basawa and Prakasa Rao (1980). We summarise this approach in some detail as we use this approach to

Let *V t*( ) denote a stochastic process having many realisations. We define the parameter set

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

*AV S t V* , , 

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

We can transform the noise process by a Hilbert space valued standard Wiener process

the norm defined by the inner product; and a separable Hilbert space should contain a

*V t S t V dt d t* ( ) , , ( ).

} be described by a stochastic partial differential equation (SPDE),

; *0tT* }. Let the evolution of the family of stochastic processes

 

( )*t* . (A Hilbert space is an inner product space that is complete with respect to

 

forms infinitely many sub event spaces with increasing times.

, and *AV*

**1.9 Parameter Estimation by Stochastic Partial Differential Equations** 

of a probability space which is given by a stochastic process ( ) *V t*

estimate the parameters in our models in this monograph.

where *A* is a partial differential operator in space, and

We can describe the stochastic process *VttT* ( ); ;

complete orthonormal sequence (Young, 1988).) Therefore,

represent a space- and time- correlated noise process.

{ () *V t* 

realisations { ( ) *V t*

; *t T* ;

of the system,

increments,

The stochastic process ( ) *V t*

differential equation (SDE),

Therefore, the stochastic process ( ) *V t*

where *S*(.) is a given function.

The explanation on the transformation of ( ,) *x t* to *d t* ( ) can be found in Jazwinski (1970), and we develop this approach further in the later chapters. A standard Wiener process (often called a Brownian motion) on the interval 0,*T* is a random variable *W t*( ) that depends continuously on *t T* 0, and satisfies the following:

$$\mathcal{W}(0) = 0,\tag{1.9.5}$$

For 0 *stT* ,

$$\mathcal{W}(t) - \mathcal{W}(s) \approx \sqrt{t - s} \text{ N}(0, 1)\_{\mathcal{H}}$$

where *N*(0,1) is a random variable generated with zero mean and unit variance.

Note that *d t* ( ) and ( ) *V t* are defined on the same event space. We estimate the parameter using the maximum likelihood approach using all the available observations of the groundwater system. The estimate ˆ *θ* of maximises the likelihood functions *V t*( ) given by (Basawa and Prakasa Rao, 1980):

$$L(\theta) = \exp\left\{ \int\_0^T S(t, \ V, \ \theta) \, dV(t) - \frac{1}{2} \Big| \prescript{T}{S}^2(t, \ V, \ \theta) \, dt \right\}. \tag{1.9.6}$$

The estimate ˆ *θ* can be obtained as the solution to the equation,

 <sup>0</sup> *L* . (1.9.7)

Maximising the likelihood function *L*( ) is equivalent to maximising the log-likelihood function, *l*( ) *= ln L*( ); hence, the maximum likelihood estimate can also be obtained as a solution to the equation

$$\frac{\partial l(\theta)}{\partial \theta} = 0 \quad . \tag{1.9.8}$$

Taking log on both sides of equation (1.9.6) we obtain,

$$d(\theta \,) = \int\_0^T S(t, \,' \, V, \, \theta) \, dV(t) \text{-} \frac{1}{2} \Big| \, \overset{\circ}{S}^2 \{t, \,' \, V, \, \theta\} dt. \tag{1.9.9}$$

The parameter is estimated as the solution to the equation

$$\int\_{0}^{\tau} \frac{\partial}{\partial \theta} S(t, V, \theta) \, dV(t) \cdot \int\_{0}^{\tau} S(t, V, \theta) \frac{\partial}{\partial \theta} S(t, V, \theta) \, dt \,\,\,\,= 0\,\,. \tag{1.9.10}$$

The parameters can be estimated from equation (1.9.10), based on a single sample path. Let us now consider the case when *M* independent sample paths are being observed. The likelihood-function becomes the product of the likelihood functions for *M* individual sample paths,

The log-likelihood function from equation (1.9.13) is,

*i M T*

 

1 0

*M T*

1 0

Differentiating the above two expressions with respect to

*i*

obtain the following two simultaneous equations:

1 0 1 0

1 0 1 0

1 and

*i i*

(1996), and in many other excellent texts.

We obtain the values for

*M T M T*

*i i*

and

*M T M T*

*l a V t a V t a V t dV t*

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

1 2 0 1 1 2 2

NonFickian Solute Transport 17

*i i i*

( , ) ( )- ( , ) ( , ) ( , ) { ( , )} 0

{ ( , )} ( )- ( , ) ( , ) ( , ) { ( , )} 0.

Over the past decades, Artificial Neural Networks (ANN) have become increasingly popular in many disciplines as a problem solving tool in data rich areas (Samarasinghe, 2006). ANN's flexible structure is capable of approximating almost any input-output relationship. Their application areas are almost limitless but fall into categories such as

ANNs are a massively parallel-distributed information processing system that has certain performance characteristics resembling biological neural networks of the human brain (Samarasinghe, 2006, Haykin, 1994). We discuss only a few of main ANN techniques that are used in this work. General detail descriptions of ANN can be found in Samarasinghe (2006), Maren et al. (1990), Hertz et al. (1991), Hegazy et al. (1994), Hassoun (1995), Rojas

Back propagation may be the most popular algorithm for training ANN in a multi-layer perceptron (MLP), which is one of many different types of neural networks. MLP comprises a number of active 'neurons' connected together to form a network. The 'strengths' or 'weights' of these links between the neurons are where the functionality of the network

Rumelhart et al. (1986) developed the standard back propagation algorithm. Since then it has undergone many modifications to overcome the limitations; and the back propagation is essentially a gradient descent technique that minimises the network error function between the output vector and the target vector. Each input pattern of the training data set is passed through the network from the input layer to the output layer. The network output is compared with the described target output, and an error is computed based on the error

*i i i i i i*

*a V t dV t a V t a V t a V t a V t dt* 

*a V t a V t a V t dt*

*i i i i*

 

> 

> >

> > >

2 as the solutions to these two equations.

, (1.9.21)

<sup>1</sup> , , ,. <sup>2</sup>

0 1 1 2 2

*i i i i i i*

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

classification, forecasting and data modelling (Maren et al., 1990; Hassoun, 1995).

**1.10 Use of Artificial Neural Networks in Parameter Estimation**

resides (NeuralWare, 1998). Its basic structure is shown in Figure 1.1.

*a V t dV t a V t a V t a V t a V t dt* 

, , , , ()

2

1 and (1.9.20)

<sup>2</sup> , respectively, we can

(1.9.22)

$$L(\theta) = L(\theta, V\_1)L(\theta, V\_2)...L(\theta, V\_M). \tag{1.9.11}$$

Taking the log on both sides of equation (1.9.11) we have the log-likelihood function,

$$l(\theta) = l(\theta, V\_1) + l(\theta, V\_2) + \dots + l(\theta, V\_M). \tag{1.9.12}$$

Using equation (1.9.10) and (1.9.12)

$$d\left(\theta\right) = \sum\_{i=1}^{M} \int\_{0}^{T} \mathbf{S}\begin{pmatrix} t\_{\prime} \ V\_{i\prime} & \theta \end{pmatrix} \,dV\_{i}\left(t\right) - \frac{1}{2} \sum\_{i=1}^{M} \int\_{0}^{T} \mathbf{S}^{2}\begin{pmatrix} t\_{\prime} \ V\_{i\prime} & \theta \end{pmatrix} dt.\tag{1.9.13}$$

Now the parameter estimate is obtained as the solution to

$$\sum\_{i=1}^{M} \left[ \frac{\partial \mathcal{S}\left(t\_{\prime} \; V\_{\prime \prime} \; \theta\right)}{\partial \theta} \; dV\_{i}(t) - \sum\_{i=1}^{M} \Big[ \mathcal{S}\left(t\_{\prime} \; V\_{\prime \prime} \; \theta\right) \frac{\partial \mathcal{S}\left(t\_{\prime} \; V\_{\prime \prime} \; \theta\right)}{\partial \theta} dt = 0. \tag{1.9.14}$$

Lets consider two particular examples in which the drift term *S*(*t, V,*  ) depends linearly on its parameters *.* 

### **Example 1**

We define the problem of estimating a single parameter as follows,

$$S(t, \ V, \ \theta) = a\_0 \begin{pmatrix} V, t \end{pmatrix} + \theta\_1 a\_1 \begin{pmatrix} V, t \end{pmatrix}; \qquad \theta\_1 \in \theta \tag{1.9.15}$$

The log-likelihood function from equation (1.9.13) is

$$I(\theta\_1) = \sum\_{i=1}^{M} \left\{ \left\langle a\_0 \left( V\_{i'}, t \right) + \theta\_1 a\_1 \left( V\_{i'}, t \right) \right\rangle dV\_i(t) - \frac{1}{2} \sum\_{i=1}^{M} \left\langle \left\langle a\_0 \left( V\_{i'}, t \right) + \theta\_1 a\_1 \left( V\_{i'}, t \right) \right\rangle \left\langle a\_1 \left( V\_{i'}, t \right) \right\rangle^2 dt = 0. \tag{1.9.16} \right\}$$

The estimate ˆ *θ* is obtained as a solution to the equation,

$$\sum\_{i=1}^{M} \int \left\{ a\_1 \left( V\_i, t \right) \right\} dV\_i(t) - \sum\_{i=1}^{M} \int \left\{ a\_0 \left( V\_i, t \right) + \theta\_1 a\_1 \left( V\_i, t \right) \right\} \left\{ a\_1 \left( V\_i, t \right) \right\} dt = 0. \tag{1.9.17}$$

Hence the estimate is given by

$$\hat{\theta} = \frac{\sum\_{i=1}^{M} \int \{a\_1(V\_{i'}, t)\} \, dV\_i(t) \text{--} \sum\_{i=1}^{M} \int \{a\_0(V\_{i'}, t)\} \left\{a\_1(V\_{i'}, t)\right\} dt}{\sum\_{i=1}^{M} \int \{a\_1^2(V\_{i'}, t)\} \, dt} \,. \tag{1.9.18}$$

### **Example 2**

When there are two unknown parameters to be estimated,

$$S(t, \ V, \ \theta) = a\_0 \begin{pmatrix} V, t \end{pmatrix} + \theta\_1 a\_1 \begin{pmatrix} V, t \end{pmatrix} + \theta\_2 a\_2 \begin{pmatrix} V, t \end{pmatrix}; \quad \theta\_1, \theta\_2 \in \theta \ . \tag{1.9.19}$$

The log-likelihood function from equation (1.9.13) is,

$$\begin{split} I\left(\theta\_{1},\theta\_{2}\right) &= \sum\_{i=1}^{M} \Big[ \left\{ a\_{0}\left(V\_{i},t\right) + \theta\_{1}a\_{1}\left(V\_{i},t\right) + \theta\_{2}a\_{2}\left(V\_{i},t\right) \right\} dV\_{i}(t) \\ &- \frac{1}{2} \sum\_{i=1}^{M} \Big[ \left\{ a\_{0}\left(V\_{i},t\right) + \theta\_{1}a\_{1}\left(V\_{i},t\right) + \theta\_{2}a\_{2}\left(V\_{i},t\right) \right\}^{2} dt. \end{split} \tag{1.9.20}$$

Differentiating the above two expressions with respect to 1 and <sup>2</sup> , respectively, we can obtain the following two simultaneous equations:

$$\sum\_{i=1}^{M} \Big| \{a\_1(V\_i, t)\} \, dV\_i(t) - \sum\_{i=1}^{M} \Big| \{a\_0(V\_i, t) + \theta\_1 a\_1(V\_i, t) + \theta\_2 a\_2(V\_i, t)\} \{a\_1(V\_i, t)\} \, dt = 0 \,, \tag{1.9.21}$$

and

Computational Modelling of Multi-Scale Non-Fickian Dispersion

) depends linearly on

(1.9.15)

. (1.9.18)

*.* (1.9.19)

 

 *<sup>M</sup>* (1.9.12)

*Lθ θ LL L , V1 2 θ, V* ...... . *θ, VM* (1.9.11)

(1.9.13)

16 in Porous Media - An Approach Based on Stochastic Calculus

Taking the log on both sides of equation (1.9.11) we have the log-likelihood function,

 , , ...... , .<sup>1</sup> <sup>2</sup>

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

<sup>0</sup> 1 1 <sup>1</sup> *St V a Vt a Vt* ( , , ) , , ;

<sup>2</sup>

*i i i i i i*

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

*i i i i i*

*a V t dV t a V t a V t a V t dt*

*aVt aVt Vt*

( ,) - ( ,) ( ,)

*i 1*

*dV (t) a dt*

*i i i*

*2 1*

( ,)

*a dt*

   

*V t*

*i*

(1.9.16)

, () , , , 0.

(1.9.17)

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

*l a V t a V t dV t a V t a V t a V t dt*

1 0 1 0 1 0

*M T*

1 0

<sup>0</sup> 1 1 2 2 1 2 *St V a Vt a Vt a Vt* ( , , ) , , , ; ,

 

*i*

*i i i i*

*St V St V dV t S t V dt*

, , , , ( ) , , 0.

*l S t V dV t S t V dt*

<sup>1</sup> , , , , . <sup>2</sup>

(1.9.14)

*i i i*

*l lVlV lV*

*M T M T*

*i i*

Using equation (1.9.10) and (1.9.12)

*.* 

its parameters

**Example 1** 

The estimate ˆ

**Example 2** 

Hence the estimate is given by

Now the parameter estimate is obtained as the solution to

1 0 1 0

Lets consider two particular examples in which the drift term *S*(*t, V,* 

We define the problem of estimating a single parameter as follows,

1 0 1 1 0 1 1 1

*θ* is obtained as a solution to the equation,

*M T M T*

*i i*

*i i*

The log-likelihood function from equation (1.9.13) is

1 0 1 0

1 0 1 0

When there are two unknown parameters to be estimated,

*i i*

*M T M T*

*i i*

ˆ *θ* =

*M T M T*

*M T M T*

$$\sum\_{i=1}^{M} \int \{a\_2(V\_i, t)\} dV\_i(t) \cdot \sum\_{i=1}^{M} \Big[ \{a\_0(V\_i, t) + \theta\_1 a\_1(V\_i, t) + \theta\_2 a\_2(V\_i, t)\} \{a\_2(V\_i, t)\} dt = 0. \tag{1.9.22}$$

We obtain the values for 1 and 2 as the solutions to these two equations.

### **1.10 Use of Artificial Neural Networks in Parameter Estimation**

Over the past decades, Artificial Neural Networks (ANN) have become increasingly popular in many disciplines as a problem solving tool in data rich areas (Samarasinghe, 2006). ANN's flexible structure is capable of approximating almost any input-output relationship. Their application areas are almost limitless but fall into categories such as classification, forecasting and data modelling (Maren et al., 1990; Hassoun, 1995).

ANNs are a massively parallel-distributed information processing system that has certain performance characteristics resembling biological neural networks of the human brain (Samarasinghe, 2006, Haykin, 1994). We discuss only a few of main ANN techniques that are used in this work. General detail descriptions of ANN can be found in Samarasinghe (2006), Maren et al. (1990), Hertz et al. (1991), Hegazy et al. (1994), Hassoun (1995), Rojas (1996), and in many other excellent texts.

Back propagation may be the most popular algorithm for training ANN in a multi-layer perceptron (MLP), which is one of many different types of neural networks. MLP comprises a number of active 'neurons' connected together to form a network. The 'strengths' or 'weights' of these links between the neurons are where the functionality of the network resides (NeuralWare, 1998). Its basic structure is shown in Figure 1.1.

Rumelhart et al. (1986) developed the standard back propagation algorithm. Since then it has undergone many modifications to overcome the limitations; and the back propagation is essentially a gradient descent technique that minimises the network error function between the output vector and the target vector. Each input pattern of the training data set is passed through the network from the input layer to the output layer. The network output is compared with the described target output, and an error is computed based on the error

ANN applications in groundwater problems are limited when compared to other disciplines in hydrology. A few of applications relevant to our work are reviewed here. Ranjithan et al. (1993) successfully used ANNs to simulate the pumping index for hydraulic conductivity realisation to remediate groundwater under uncertainty. In the process of designing a reliable groundwater remediation strategy, clear identification of heterogeneous spatial variability of the hydrology parameters is an important issue. The association of hydraulic conductivity patterns and the level of criticalness need to be understood sufficiently for efficient screening. ANNs have been used to recognize and classify the variable patterns (Ranjithan et al., 1993). Similar work has been conducted by Rogers and Dowla (1994) to simulate a regulatory index for multiple pumping realizations at a contaminated site. In this study the supervised learning algorithm of back propagation has been used to train a network. The conjugate gradient method and weight elimination procedures have been employed to speed up the convergence and improve the performance, respectively. After training the networks, the ANN begins a search through various realizations of pumping patterns to determine matching patterns. Rogers et al. (1995) took another step forward to simulate the regulatory index, remedial index and cost index by using ANN for groundwater remediation. This research contributed towards addressing the issue of

NonFickian Solute Transport 19

Zhu (2000) used ANN to develop an approach to populate a soil similarity model that was designed to represent soil landscape as spatial continua for hydrological modelling at watershed of mesoscale size. Coulibaly et al. (2001) modelled the water table depth fluctuations by using three types of functionally different ANN models: Input Delay Neural Network (IDNN), Recurrent Neural Network (RNN) and Radial Basis Function Network (RBFN). This type of study has significant implications for groundwater management in the areas with inadequate groundwater monitoring networks (Maier and Dandy, 2000). Hong and Rosen (2001) demonstrated that the unsupervised self-organising map was an efficient tool for diagnosing the effect of the storm water infiltration on the groundwater quality variables. In addition, they showed that SOM could also be useful in extracting the

Balkhair (2002) presented a method for estimating the aquifer parameters in large diameter wells using ANN. The designed network was trained to learn the underlying complex relationship between input and output patterns of the normalized draw down data generated from an analytical solution and its corresponding transmissivity values. The ANN was trained with a fixed number of input draw down data points obtained from the analytical solution for a pre-specified ranges of aquifer parameter values and time-series data. The trained network was capable of producing aquifer parameter values for any given input pattern of normalized draw down data and well diameter size. The values of aquifer parameters obtained using this approach were in a good agreement with those obtained by other published results. Prior knowledge about the aquifer parameter values has served as a

Rudnitskaya et al. (2001) developed a methodology to monitor groundwater quality using an array of non-specific potentiometric chemical sensors with data processing by ANN. Lischeid (2001) studied the impact of long-lasting non-point emissions on groundwater and stream water in remote watersheds using a neural network approach. Scarlatos (2001) used ANN method to identify the sources, distribution and fate of fecal coliform populations in

dependencies between the variables in a given groundwater quality dataset.

valuable piece of information in this ANN approach.

escalating costs of environmental cleanup.

function. This error is propagated backward through the network to each node, and correspondingly the connection weights are adjusted.

Figure 1.1. Basic structure of a multi-layer perceptron network.

The Self-Organizing Map (SOM) was developed by Kohonen (1982) and arose from the attempts to model the topographically organized maps found in the cortices of the more developed animal brains. The underlying basis behind the development of the SOM was that topologically correct maps can be formed in an *n*-dimensional array of processing elements that did not have this initial ordering to begin with. In this way, input stimuli, which may have many dimensions, can cluster to be represented by a one or twodimensional vector which preserves the order of the higher dimensional data (NeuralWare, 1998)*.* The SOM employs a type of learning commonly referred to as competitive, unsupervised or self-organizing, in which adjacent cells within the network are able to interact and adaptively evolved into the detectors of a specific input pattern (Kohonen, 1990). The SOM can be considered to be *"neural"* because the results have indicated that the adaptive processes utilized in the SOM may be similar to the processes at work within the brain (Kohonen, 1990). The SOM has the potential for extending its capability beyond the original purpose of modelling biological phenomena. Sorting items into categories of similar objects is a challenging, yet frequent task. The SOM achieves this task by nonlinearly projecting the data onto a lower dimensional display and by clustering the data (Kohonen, 1990). This attribute has been used in a wide number of applications ranging from engineering (including image and signal processing, image recognition, telecommunication, process monitoring and control, and robotics) to natural sciences, medicine, humanities, economics and mathematics (Kaski et al., 1998). cluster represented which to the medicine, the extracting on

### **1.11 ANN Applications in Hydrology**

It has been shown that ANN's flexible structure can provide simple and reasonable solutions to various problems in hydrology. Since the beginning of the last decade, ANN have been successfully employed in hydrology research such as rainfall-runoff modelling, stream flow forecasting, precipitation forecasting, groundwater modelling, water quality and management modelling (Morshed and Kaluarachchi, 1998; ASCE Task Committee on Application of ANN in Hydrology, 2000a, b; Maier and Dandy, 2000).

18 in Porous Media - An Approach Based on Stochastic Calculus

function. This error is propagated backward through the network to each node, and

The Self-Organizing Map (SOM) was developed by Kohonen (1982) and arose from the attempts to model the topographically organized maps found in the cortices of the more developed animal brains. The underlying basis behind the development of the SOM was that topologically correct maps can be formed in an *n*-dimensional array of processing elements that did not have this initial ordering to begin with. In this way, input stimuli, which may have many dimensions, can cluster to be represented by a one or twodimensional vector which preserves the order of the higher dimensional data (NeuralWare, 1998)*.* The SOM employs a type of learning commonly referred to as competitive, unsupervised or self-organizing, in which adjacent cells within the network are able to interact and adaptively evolved into the detectors of a specific input pattern (Kohonen, 1990). The SOM can be considered to be *"neural"* because the results have indicated that the adaptive processes utilized in the SOM may be similar to the processes at work within the brain (Kohonen, 1990). The SOM has the potential for extending its capability beyond the original purpose of modelling biological phenomena. Sorting items into categories of similar objects is a challenging, yet frequent task. The SOM achieves this task by nonlinearly projecting the data onto a lower dimensional display and by clustering the data (Kohonen, 1990). This attribute has been used in a wide number of applications ranging from engineering (including image and signal processing, image recognition, telecommunication, process monitoring and control, and robotics) to natural sciences, medicine, humanities,

It has been shown that ANN's flexible structure can provide simple and reasonable solutions to various problems in hydrology. Since the beginning of the last decade, ANN have been successfully employed in hydrology research such as rainfall-runoff modelling, stream flow forecasting, precipitation forecasting, groundwater modelling, water quality and management modelling (Morshed and Kaluarachchi, 1998; ASCE Task Committee on

Application of ANN in Hydrology, 2000a, b; Maier and Dandy, 2000).

correspondingly the connection weights are adjusted.

Figure 1.1. Basic structure of a multi-layer perceptron network.

economics and mathematics (Kaski et al., 1998).

**1.11 ANN Applications in Hydrology** 

ANN applications in groundwater problems are limited when compared to other disciplines in hydrology. A few of applications relevant to our work are reviewed here. Ranjithan et al. (1993) successfully used ANNs to simulate the pumping index for hydraulic conductivity realisation to remediate groundwater under uncertainty. In the process of designing a reliable groundwater remediation strategy, clear identification of heterogeneous spatial variability of the hydrology parameters is an important issue. The association of hydraulic conductivity patterns and the level of criticalness need to be understood sufficiently for efficient screening. ANNs have been used to recognize and classify the variable patterns (Ranjithan et al., 1993). Similar work has been conducted by Rogers and Dowla (1994) to simulate a regulatory index for multiple pumping realizations at a contaminated site. In this study the supervised learning algorithm of back propagation has been used to train a network. The conjugate gradient method and weight elimination procedures have been employed to speed up the convergence and improve the performance, respectively. After training the networks, the ANN begins a search through various realizations of pumping patterns to determine matching patterns. Rogers et al. (1995) took another step forward to simulate the regulatory index, remedial index and cost index by using ANN for groundwater remediation. This research contributed towards addressing the issue of escalating costs of environmental cleanup. higher specific decade, other on

Zhu (2000) used ANN to develop an approach to populate a soil similarity model that was designed to represent soil landscape as spatial continua for hydrological modelling at watershed of mesoscale size. Coulibaly et al. (2001) modelled the water table depth fluctuations by using three types of functionally different ANN models: Input Delay Neural Network (IDNN), Recurrent Neural Network (RNN) and Radial Basis Function Network (RBFN). This type of study has significant implications for groundwater management in the areas with inadequate groundwater monitoring networks (Maier and Dandy, 2000). Hong and Rosen (2001) demonstrated that the unsupervised self-organising map was an efficient tool for diagnosing the effect of the storm water infiltration on the groundwater quality variables. In addition, they showed that SOM could also be useful in extracting the dependencies between the variables in a given groundwater quality dataset.

Balkhair (2002) presented a method for estimating the aquifer parameters in large diameter wells using ANN. The designed network was trained to learn the underlying complex relationship between input and output patterns of the normalized draw down data generated from an analytical solution and its corresponding transmissivity values. The ANN was trained with a fixed number of input draw down data points obtained from the analytical solution for a pre-specified ranges of aquifer parameter values and time-series data. The trained network was capable of producing aquifer parameter values for any given input pattern of normalized draw down data and well diameter size. The values of aquifer parameters obtained using this approach were in a good agreement with those obtained by other published results. Prior knowledge about the aquifer parameter values has served as a valuable piece of information in this ANN approach. from by

Rudnitskaya et al. (2001) developed a methodology to monitor groundwater quality using an array of non-specific potentiometric chemical sensors with data processing by ANN. Lischeid (2001) studied the impact of long-lasting non-point emissions on groundwater and stream water in remote watersheds using a neural network approach. Scarlatos (2001) used ANN method to identify the sources, distribution and fate of fecal coliform populations in

to

to denote the set of all such values.

is observed for *Y*, this is called an

may be a finite set of

). By using the standard

values. For this to

(the complementary event), or that a value

*, F, P )*" .

**Stochastic Differential Equations** 

As we have discussed in chapter 1, the deterministic mathematical formulation of solute transport through a porous medium introduces the dispersivity, which is a measure of the distance a solute tracer would travel when the mean velocity is normalized to be one. One would expect such a measure to be a mechanical property of the porous medium under consideration, but the evidence are there to show that dispersivity is dependent on the scale of the experiment for a given porous medium. One of the challenges in modelling the phenomena is to discard the Fickian assumptions, through which dispersivity is defined, and develop a mathematical discription containing the fluctuations associated with the mean velocity of a physical ensemble of solute particles. To this end, we require a sophisticated mathematical framework, and the theory of stochastic processes and differential equations is a natural mathematical setting. In this chapter we review some essential concepts in stochastic processes and stochastic differential equations in order to

A deterministic variable expressed as a function of time uniquely determines the value of the variable at a given time. A stochastic variable *Y*, on the other hand, is one that does not have a unique value; it can have any one out of a set of values. We assign a unique label

discrete numbers, and when *Y* is the instantaneous position of a fluid particle, it may be a

event *F*. In fact, this is only the simplest prototype of an event; other possibilities might be

Even though the outcome of a particular observation of *Y* is unpredictable, the probability of

methods of probability calculus, this implies that a probability *P*(*F*) can also be assigned to

work, *F* must satisfy the criteria that for any event *F* in its complement *Fc* must also belong to *F*, and that for any subset of *F'*s the union of these must also belong to *F*. The explanation above of what it means to call *Y* a stochastic variable, is encapsulated in formal

In describing physical systems, deterministic variables usually depend on additional parameters such as time. Similarly, a stochastic variable may depend on an additional

must be determined by a probability function *P*(

compound events *F* e.g. by appropriate summation or integration over

mathematical language by saying "*Y* is defined on a probability space *(*

values is observed. The set of all possible events is denoted by *F.* 

**2.1 Concepts in Stochastic Calculus** 

understand the stochastic calculus in a more applied context.

When *Y* represents, for example the outcome of throwing dice,

each possible value of the stochastic variable, and set

continuous range of real numbers. If a particular value *y*

that the value of *Y* is observed not to be *y*

within a certain range of

observing *y*

**and Related Inverse Problems** 

the North Fork of the New River that flows through the City of Fort Lauderdale, Florida, USA and how the storm water drainage from sewers affects the groundwater. Other ANN applications in water resources can be found in Aly and Peralta (1999), Mukhopadhyay (1999), Freeze and Gorelick (2000), Johnson and Rogers (2000), Hassan and Hamed (2001), Beaudeau et al. (2001), and Lindsay et al. (2002). ANN variable
