**Multiscale, Generalised Stochastic Solute Transport Model in One Dimension**

### **6.1 Introduction**

In Chapter 3 and 4, we have developed a stochastic solute transport model in 1-D without rosorting to simplifying Fickian assumptions, but by using the idea that the fluctuations in velocity are influenced by the nature of porous medium. We model these fluctuations through the velocity covariance kernel. We have also estimated the dispersivity by taking the realisations of the solution of the SSTM and using them as the observations in the stochastic inverse method (SIM) based on the maximum likelihood estimation procedure for the stochastic partial differential equation obtained by adding a noise term to the advectiondispersion equation. We have confined the estimation of dispersitivities to a flow length of 1 m (i.e, *x* 0,1 ) except in Chapter 3, section 3.10, where we have estimated the dispersitivities up to 10 km using the SIM by simplifying the SSTM. This approach was proven to be computationally expensive and the approximation of the SSTM we have developed was based on the spatial average of the variance of the fluctuation term over the flow length. Further, the solution is based on a specific kernel. This development in Chapter 3 is inadequate to examine the scale dependence of the dispersitivity. Therefore, we set out to develop a dimensionless model for any given arbitrary flow length, *L* , in this Chapter for any given velocity kernel provided that we have the eigen functions in the form given by equation (4.2.3). Then we examine the dispersivities in relation to the flow lengths to understand the multi-scale behaviour of the SSTM.

The starting point of the development of the multi-scale SSTM is the Langevin equation for the SSTM, which is interpreted locally. From equation (4.9.1), the Langevin equation can be written as,

$$d\mathbb{C}\_{\mathbf{x}}(t) = -\alpha\_{\mathbf{x}}(\mathbb{C}\_{\mathbf{x}}(t), \overline{V}(\mathbf{x}, t), \mathbf{x})dt + \beta\_{\mathbf{x}}(\mathbb{C}\_{\mathbf{x}}(t), \frac{\partial \mathbb{C}\_{\mathbf{x}}}{\partial \mathbf{x}}, \frac{\partial^2 \mathbb{C}\_{\mathbf{x}}}{\partial \mathbf{x}^2}, \mathbf{x}) dw(t) \tag{6.1.1}$$

where the coefficients *x* and *<sup>x</sup>* are dependent on , () *<sup>x</sup> xC t* and *Vxt* ( ,) ; and 2 <sup>2</sup> ( ), , *x x x C C C t x x* and *x* , respectively. *dw t*( ) are the standard Wiener increments with zero-mean and *dt* variance. As discussed in Chapter 4, equation (6.1.1) has to be interpreted carefully to understand it better. Equation (6.1.1) is a SDE and also an Ito diffusion with the coefficients depending on the functions of space variables. It gives us the time evolution of the concentration of solute at a given point *x* which is denoted by subscript *x* . Obviously, the computation of *Cx* also depends on how the spatial

*j j j kj*

*P x g gx ge*

*h*

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

2

*j*

*p*

*k*

2 0 1

terms of discretized times, 0 1 <sup>1</sup> , ,..., , ,..., *i i t t tt* equation (6.1.1) can be written as,

account in developing numerical algorithms to solve it.

where the drift coefficient,

*<sup>x</sup>* , and the diffusion coefficient,

This restrictive nature of equation (6.1.14) in evaluating the coefficient has to be taken in to

*x x x x xi i i i*

1 , ,, ,

*V V dC t C t C t C t V x t t t*

*x xi i i*

*V V C t V xt d t*

, ,, ,

 

2

*i i*

*t t*

*x x*

*x x*

2 1

*<sup>x</sup>* , are evaluated at time *<sup>i</sup> t* .

(6.1.14)

2 2

*i i*

*t t*

,

2 2

2

*j j j kj*

Equation (6.1.1) to (6.1.13) constitute the Langevin form of the SSTM. It should be noted that the functions *Pij* are only valid for *x* 0,1. If we normalize the spatial variable *x* to remain with in 0,1 , then we can use the results in Chapter 4 to obtain *Pij* . We develop

One should note that the Langevin equation for any system reflect the role of external noise to the system under consideration (van Kampen, 1992). Even though we have derived equation (6.1.1) starting from the mass conservation of solute particles, the fluctuations associate with hydrodynamics dispersion are a result of dissipation of energy of particles due to momentum changes associated near to the surfaces of porous medium. For *a* physical ensemble of solute particles, porous medium through which it flows act as an external source of noise. From this point of review, the Langevin type equation for solute concentration is justified. As a SDE, equation (6.1.1) is a Wiener process with stochastic, at best nonlinear, time-dependent coefficients, and it is also an Ito diffusion which should be interpreted locally, i.e., for a given *x* and *t* , equation (6.1.1) is valid only for short time intervals beyond *t* . This naturally leads us to evaluate the associated spatial derivatives at the previous time, which is valid according to Ito's interpretation of stochastic integral. In

*<sup>h</sup> P x g g x ge*

*j*

*p r xs x ij kj kj kj k*

*j j*

*p p r xs r xs x kj kj kj kj kj*

Multiscale, Generalised Stochastic Solute Transport Model in One Dimension 179

*k k*

2

*j*

*p*

2

*h*

*j ij kj kj kj k*

1 0 1

the dimensionless Langevin form of the SSTM in section 6.2.

*P x g gr x s e*

0

and

2

2 2 , <sup>2</sup>

2 2

2

*kj kj*

*kj kj*

(6.1.13)

(6.1.11)

(6.1.12)

*kj kj kj kj*

2 . <sup>2</sup> *j*

*p r xs x*

*k*

2

*kj kj*

*r xs*

<sup>4</sup> <sup>2</sup> , <sup>2</sup>

*gr x s e gre*

*kj kj*

*r xs*

*g gr x s e*

<sup>2</sup>

derivatives of *Cx* are calculated. In that sense, equation (6.1.1) is a stochastic partial differential equation as the coefficients are functions of random quantities. But we avoid solving a SPDE by treating equation (6.1.1) as a SDE and interpreting it as an Ito integral which makes us to evaluate coefficients at the previous time point with respect to the current point of evaluation.

For simplicity, we will denote the coefficients as *x* and *<sup>x</sup>* . In Chapter 4, we have derived explicit function for *x* and *<sup>x</sup>* :

$$
\delta \mathbf{x} = \mathbf{C}\_{\mathbf{x}} \left( t \right) F\_{\mathbf{x},0} + \left( \frac{\partial \mathbf{C}}{\partial \mathbf{x}} \right)\_{\mathbf{x}} F\_{\mathbf{x},1} + \left( \frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2} \right)\_{\mathbf{x}} F\_{\mathbf{x},2}. \tag{6.1.2}
$$

Where

$$F\_{\mathbf{x},0} = \frac{\partial \overline{V}(\mathbf{x},t)}{\partial \mathbf{x}} + \frac{\mu\_{\mathbf{x}}}{2} \frac{\partial^2 \overline{V}(\mathbf{x},t)}{\partial \mathbf{x}^2},\tag{6.1.3}$$

$$F\_{\mathbf{x},1} = \overline{V}\left(\mathbf{x}, t\right) + h\_{\mathbf{x}} \frac{\partial \overline{V}\left(\mathbf{x}, t\right)}{\partial \mathbf{x}},\tag{6.1.4}$$

and,

$$F\_{\mathbf{x},3} = \frac{h\_{\mathbf{x}}}{2} \overline{V}(\mathbf{x}, t);\tag{6.1.5}$$

$$
\beta\_x = \left(\beta\_0^2 + \beta\_1^2 + \beta\_2^2\right)^{\frac{1}{2}},\tag{6.1.6}
$$

where,

$$
\mathcal{J}\_0 = \mathbb{C}\_{\times}(t) \sqrt{a\_{\infty}} \tag{6.1.7}
$$

$$\mathcal{J}\_1 = \left(\frac{\partial \mathbf{C}}{\partial \mathbf{x}}\right)\_\mathbf{x} \sqrt{a\_{11'}} \tag{6.1.8}$$

$$
\beta\_2 = \left(\frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2}\right)\_\times \sqrt{a\_{22}} \,\,\,\,\tag{6.1.9}
$$

and

$$a\_{\rm ii} = \sigma^2 \sum\_{j=1}^{m} \mathbb{X}\_j P\_{\rm ij}^2 \, \, \, \quad \qquad \left( \text{i}, 0, 1, 2 \right), \tag{6.1.10}$$

In equation (6.1.10), <sup>2</sup> is the variance of the covariance kernel, *<sup>j</sup>* are eigen functions, and for the domain of *x* 0,1,

$$\begin{split} P\_{0j}\left(\mathbf{x}\right) &= \left[ \mathcal{G}\_{il} - 2\sum\_{k=2}^{p\_j} \mathcal{G}\_{kj} r\_{kj} \left( \mathbf{x} - \mathbf{s}\_{kj} \right) e^{-\boldsymbol{\eta}\_j \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2} \right] \\ &+ \left( \frac{\boldsymbol{h}\_x}{2} \right) \left[ 4 \sum\_{k=2}^{p\_j} \mathcal{G}\_{kj} r\_{kj}^2 \left( \mathbf{x} - \mathbf{s}\_{kj} \right)^2 e^{-\boldsymbol{\eta}\_j \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2} - 2 \sum\_{k=2}^{p\_j} \mathcal{G}\_{kj} r\_{kj} e^{-\boldsymbol{\eta}\_j \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2} \right], \end{split} \tag{6.1.11}$$
 
$$P\_{1j}\left(\mathbf{x}\right) = \mathcal{G}\_{0j} + \mathcal{G}\_{1j}\mathbf{x} + \sum\_{k=2}^{p\_j} \mathcal{G}\_{kj} e^{-\boldsymbol{\eta}\_k \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2} \tag{6.1.12}$$

$$\begin{split} \mathbf{x} \mathbf{y} &= \mathbf{g}\_{0j} + \mathbf{g}\_{1j}\mathbf{x} + \sum\_{k=2} \mathbf{g}\_{kj} e^{-\mathbf{r}\_{kj}^{-\left(1-\mathbf{s}\_{kj}\right)}} \\ &+ \left(\frac{\mathbf{h}\_{\mathbf{x}}}{2}\right) \left[ \mathbf{2} \left( \mathbf{g}\_{ij} - 2 \sum\_{k=2}^{p\_j} \mathbf{g}\_{kj} r\_{kj} \left( \mathbf{x} - \mathbf{s}\_{kj} \right) e^{-\mathbf{r}\_{kj} \left( \mathbf{x} - \mathbf{s}\_{kj} \right)^2} \right) \right] \end{split} \tag{6.1.12}$$

and

Computational Modelling of Multi-Scale Non-Fickian Dispersion

(6.1.3)

*<sup>h</sup> F V xt* (6.1.5)

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

*Ct a <sup>x</sup>* (6.1.7)

(6.1.8)

(6.1.9)

*<sup>j</sup>* are eigen functions,

(6.1.10)

(6.1.4)

*<sup>x</sup>* . In Chapter 4, we have

(6.1.2)

178 in Porous Media - An Approach Based on Stochastic Calculus

derivatives of *Cx* are calculated. In that sense, equation (6.1.1) is a stochastic partial differential equation as the coefficients are functions of random quantities. But we avoid solving a SPDE by treating equation (6.1.1) as a SDE and interpreting it as an Ito integral which makes us to evaluate coefficients at the previous time point with respect to the

> *x* and

*x x*

*x x*

<sup>2</sup>

,3 , ; <sup>2</sup> *x*

 <sup>1</sup> <sup>222</sup> <sup>2</sup> <sup>012</sup> ,

<sup>0</sup> <sup>00</sup> ,

<sup>1</sup> <sup>11</sup> , *x*

*a*

*a*

, ,0,1,2 ,

*C*

*x*

2 <sup>2</sup> <sup>2</sup> <sup>22</sup> , *x*

*x*

*C*

2 2

is the variance of the covariance kernel,

, , , *<sup>x</sup> <sup>x</sup> V xt F V xt h*

, , , <sup>2</sup> *x*

*x*

current point of evaluation.

derived explicit function for

Where

and,

where,

and

In equation (6.1.10), <sup>2</sup>

and for the domain of *x* 0,1,

For simplicity, we will denote the coefficients as

*x* and

*x*

,1

*<sup>x</sup>* :

 <sup>2</sup> ,0 ,1 <sup>2</sup> ,2 . *x x <sup>x</sup> <sup>x</sup>*

*C C C tF F <sup>F</sup>*

,0 2

*x*

1

*aPi*

*m ii j ij j*

 

*V xt V xt <sup>h</sup> <sup>F</sup> x x*

 

$$P\_{2\_{j}}\left(\mathbf{x}\right) = \left(\frac{h\_{\mathbf{x}}}{2}\right) \left[ \mathbf{g}\_{0j} + \mathbf{g}\_{1j}\mathbf{x} + \sum\_{k=2}^{p\_{j}} \mathbf{g}\_{kj} e^{-r\_{\mathbf{y}}\left(\mathbf{x} - \mathbf{s}\_{k}\right)^{2}}\right].\tag{6.1.13}$$

Equation (6.1.1) to (6.1.13) constitute the Langevin form of the SSTM. It should be noted that the functions *Pij* are only valid for *x* 0,1. If we normalize the spatial variable *x* to remain with in 0,1 , then we can use the results in Chapter 4 to obtain *Pij* . We develop the dimensionless Langevin form of the SSTM in section 6.2. 

One should note that the Langevin equation for any system reflect the role of external noise to the system under consideration (van Kampen, 1992). Even though we have derived equation (6.1.1) starting from the mass conservation of solute particles, the fluctuations associate with hydrodynamics dispersion are a result of dissipation of energy of particles due to momentum changes associated near to the surfaces of porous medium. For *a* physical ensemble of solute particles, porous medium through which it flows act as an external source of noise. From this point of review, the Langevin type equation for solute concentration is justified. As a SDE, equation (6.1.1) is a Wiener process with stochastic, at best nonlinear, time-dependent coefficients, and it is also an Ito diffusion which should be interpreted locally, i.e., for a given *x* and *t* , equation (6.1.1) is valid only for short time intervals beyond *t* . This naturally leads us to evaluate the associated spatial derivatives at the previous time, which is valid according to Ito's interpretation of stochastic integral. In terms of discretized times, 0 1 <sup>1</sup> , ,..., , ,..., *i i t t tt* equation (6.1.1) can be written as,

$$\begin{split} d\mathbf{C}\_{x}(t) = \mathbf{C}\_{x}(t+1) - \mathbf{C}\_{x}(t) &= a\_{x} \Big( \mathbf{C}\_{x}(t\_{i}), \overline{V}(\mathbf{x}, t\_{i}), \Big( \frac{\partial \overline{V}}{\partial \mathbf{x}} \Big)\_{l\_{i}}, \Big( \frac{\partial^{2} \overline{V}}{\partial \mathbf{x}^{2}} \Big)\_{l\_{i}} \Big) (t\_{i+1} - t\_{i}) \\ &+ (\boldsymbol{\beta}\_{x}) \Big( \mathbf{C}\_{x}(t\_{i}), \overline{V}(\mathbf{x}, t\_{i}), \Big( \frac{\partial \overline{V}}{\partial \mathbf{x}} \Big)\_{l\_{i}}, \Big( \frac{\partial^{2} \overline{V}}{\partial \mathbf{x}^{2}} \Big)\_{l\_{i}} \Big) d o(t\_{i}) \end{split} \tag{6.1.14}$$

where the drift coefficient, *<sup>x</sup>* , and the diffusion coefficient, *<sup>x</sup>* , are evaluated at time *<sup>i</sup> t* . This restrictive nature of equation (6.1.14) in evaluating the coefficient has to be taken in to account in developing numerical algorithms to solve it.

### **6.2 Partially Dimensionless SSTM with Flow Length** *L*

We start the derivation of partially dimensionless SSTM by defining the dimensionless distance, *Z* , as:

$$\mathbf{Z} = \bigvee \mathbf{L} \tag{6.2.1}$$

Similarly,

where

Now we can write equation (6.1.1) in the following manner:

0

Therefore, the Langevin form of the generalized SSTM is given by

0

(6.2.14) first before discussing a completely dimensionless equation.

.

*Z C* 

2

0

and

*Z*

0 *Z <sup>Z</sup> C* 

where

for any given *L* .

moles within the matrix.

 <sup>0</sup> <sup>1</sup> <sup>11</sup> , *<sup>C</sup> a*

Multiscale, Generalised Stochastic Solute Transport Model in One Dimension 181

 <sup>2</sup> 0 2 2 2 22 *C*

*a*

 

,0 ,1 2 2 ,2

00 2 11 4 2 22

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

, 0 1. (6.2.14)

 

0 0 *<sup>Z</sup> <sup>d</sup> dt d t C C*

<sup>2</sup>

*<sup>Z</sup>* 1 1 *FF F C L L*

*d dt d t Z*

Using equation (6.2.14) we can compute the time course of the dimensionless concentration

The dimensionless/concentration, , varies from 0 to 1. *C t <sup>Z</sup>* is proportioned to the number of solute moles within a unit volume of porous/water matrix, and *C*0 is proportional to the maximum possible number of solute moles within the same matrix. Therefore, *Ct C <sup>Z</sup>* / 0 can be interpreted as the likelihood (probability) of finding solute

It should be noted that time, *t* , is not a dimensionless quantity and therefore, equation (6.2.14) is partially dimensionless equation. We will explore the dispersivity using equation

 

*<sup>Z</sup> <sup>Z</sup>* <sup>1</sup> <sup>1</sup> *a Z <sup>a</sup> <sup>a</sup>*

*C L L*

*L*

*L*

*d C* <sup>0</sup> 

<sup>0</sup> *C ta* <sup>0</sup> 00 (6.2.7)

and (6.2.8)

(6.2.9)

*<sup>Z</sup> dt d t* , (6.2.10)

, (6.2.12)

(6.2.13)

. (6.2.11)

where *L* is the total flow length.

When *x L* 0, , 0,1 .

If *C*0 is a constant concentration defined such a way that *C*<sup>0</sup> maximum of *C t <sup>x</sup>* for all *x* and *t* , then *C Ct* <sup>0</sup> *<sup>x</sup>* for any *t* and *x* . We can define dimensional concentration *t* as,

$$
\Gamma \begin{pmatrix} t \\ \end{pmatrix} = \frac{\mathbf{C}\_x}{\mathbf{C}\_o} \qquad . \tag{6.2.2}
$$

From equation (6.2.1),

$$\frac{\partial Z}{\partial \mathbf{x}} = \frac{1}{L} \,\mathrm{\,\,\,}\tag{6.2.3a}$$

$$\frac{\partial \mathbf{C}\_{\mathbf{x}}}{\partial \mathbf{x}} = \frac{\partial \{\mathbf{C}\_{0} \Gamma\}}{\partial \mathbf{Z}} \cdot \frac{\partial \mathbf{Z}}{\partial \mathbf{x}} = \frac{\mathbf{C}\_{0}}{\mathbf{L}} \frac{\partial \Gamma}{\partial \mathbf{Z}},\tag{6.2.3b}$$

$$\frac{\partial^2 \mathbb{C}\_x}{\partial \mathbf{x}^2} = \frac{\partial}{\partial \mathbf{x}} \left( \frac{\partial \mathbb{C}\_x}{\partial \mathbf{x}} \right) = \frac{\partial}{\partial \mathbf{x}} \left( \frac{\mathbb{C}\_0}{L} \frac{\partial \Gamma}{\partial \mathbf{Z}} \right) = \frac{\partial}{\partial \mathbf{Z}} \left( \frac{\mathbb{C}\_0}{L} \frac{\partial \Gamma}{\partial \mathbf{Z}} \right) \frac{\partial \mathbb{Z}}{\partial \mathbf{x}} = \frac{\mathbb{C}\_0}{L^2} \frac{\partial^2 \Gamma}{\partial \mathbf{Z}^2}. \tag{6.2.3c}$$

As the domain of *x* is the generalized SSTM is from 0 to 1, we can replace *x* with *Z* in the dimensionless generalized SSTM. For example, *Fx*,0 becomes *FZ*,0 .

$$\frac{\partial C\_c}{\partial x} = \frac{\partial (C\_c \overline{\alpha})}{\partial \overline{L}} \frac{\partial C\_c}{\partial x} = \frac{C\_c}{L} \frac{\partial L}{\partial Z},\tag{6.2.2b}$$

$$\frac{\partial^2 C\_c}{\partial x^2} = \frac{\partial}{\partial x} \left(\frac{C\_c}{\overline{\alpha}}\right) = \frac{\partial}{\overline{\alpha}} \left(\frac{C\_c}{L} \frac{\partial C}{\partial Z}\right) = \frac{\partial}{\overline{\alpha}Z} \left(\frac{C\_c}{L} \frac{\partial C}{\partial Z}\right) \frac{\partial Z}{\partial x} = \frac{C\_c}{L} \frac{\partial^2 C}{\partial Z^2},\tag{6.2.2c}$$

$$\text{for } x \text{ is the generalized STSM is from } 0 \text{ to } 1, \text{we can replace } x \text{ with } Z. \text{ in these generalized STSM for example, } F\_{L,\nu} \text{ becomes}$$

$$F\_{C,\nu} = \frac{\partial^2 \overline{F}(Z)}{\partial x} + \frac{\partial}{\overline{\alpha}} \frac{\partial^2 \overline{F}(Z)}{\partial z} = \frac{\partial \overline{F}(Z)}{\partial Z} \frac{\partial Z}{\partial x} + \frac{\partial}{\overline{\alpha}} \frac{\partial}{\partial z} \left(\frac{Z}{\overline{\alpha}} \frac{Z}{\overline{\alpha}}\right)$$

$$= \frac{1}{2} \frac{\partial \overline{V}}{\partial Z} + \frac{\hbar}{2} \frac{\partial}{\overline{\alpha}} \left(\frac{\partial \overline{V}}{\partial \overline{\alpha}}\right) \left(\frac{Z}{\overline{\alpha}}\right) = \frac{1}{2} \frac{\partial \overline{V}}{\partial Z} + \frac{\hbar}{2} \frac{\partial}{\overline{\alpha}} \frac{\partial^2 \overline{V}}{\partial Z}.$$

Similarly,

$$F\_{\mathbf{Z},1} = \overline{V}\left(\mathbf{Z}, \mathbf{t}\right) + \frac{h\_{\mathbf{Z}}}{\mathbf{L}} \frac{\partial \overline{V}}{\partial \mathbf{Z}}; \text{and} \tag{6.2.5}$$

$$F\_{\mathbf{z},2} = \frac{h\_{\mathbf{z}}}{2} \overline{V}(\mathbf{Z}, \mathbf{t}). \tag{6.2.6}$$

0 1 , *P P <sup>j</sup> <sup>j</sup>* and *P*<sup>2</sup> *<sup>j</sup>* are obtained by simply replacing *x* in 0 1 , *P xP x j j* and *P x* <sup>2</sup> *<sup>j</sup>* expressions by , because these expressions are derived for 0,1 domain.

Similarly,

Computational Modelling of Multi-Scale Non-Fickian Dispersion

*<sup>L</sup>* (6.2.1)

. (6.2.2)

(6.2.3a)

(6.2.4)

and (6.2.5)

(6.2.6)

(6.2.3b)

2 2 2

180 in Porous Media - An Approach Based on Stochastic Calculus

We start the derivation of partially dimensionless SSTM by defining the dimensionless

*x*

If *C*0 is a constant concentration defined such a way that *C*<sup>0</sup> maximum of *C t <sup>x</sup>* for all *x* and *t* , then *C Ct* <sup>0</sup> *<sup>x</sup>* for any *t* and *x* . We can define dimensional concentration

0 *Cx <sup>t</sup> C*

> <sup>1</sup> , *<sup>Z</sup> x L*

 <sup>0</sup> <sup>0</sup> , *<sup>x</sup> <sup>C</sup> <sup>C</sup> Z C x Z x LZ*

0 0 0

(6.2.3c)

2 2

<sup>2</sup> 2 2 . . *CCC <sup>x</sup> <sup>x</sup> C ZC x x x xLZ ZLZ x L*

As the domain of *x* is the generalized SSTM is from 0 to 1, we can replace *x* with *Z* in

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

,1 , ; *h V F Vt*

*L* 

,2 , . <sup>2</sup> *<sup>h</sup> F Vt*

0 1 , *P P <sup>j</sup> <sup>j</sup>* and *P*<sup>2</sup> *<sup>j</sup>* are obtained by simply replacing *x* in 0 1 , *P xP x j j* and

*P x* <sup>2</sup> *<sup>j</sup>* expressions by , because these expressions are derived for 0,1 domain.

 

*x x Z x x Zx*

the dimensionless generalized SSTM. For example, *Fx*,0 becomes *FZ*,0 .

,0 2

<sup>2</sup>

*L LZ Z x L L*

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

*Vh V Z V h V*

*Vt Vt Vt <sup>h</sup> Zh VZ <sup>F</sup>*

**6.2 Partially Dimensionless SSTM with Flow Length** *L*

distance, *Z* , as:

*t* as,

Similarly,

where *L* is the total flow length.

When *x L* 0, , 0,1 .

From equation (6.2.1),

$$\mathcal{B}\_0(Z) = \mathbb{C}\_0 \Gamma(t) \sqrt{a\_{\text{co}}(Z)} \tag{6.2.7}$$

$$\mathcal{B}\_1(Z) = \frac{\mathbb{C}\_0}{L} \frac{\partial \Gamma}{\partial Z} \sqrt{a\_{11}(Z)}, \quad \text{and} \tag{6.2.8}$$

$$\mathcal{J}\_2\left(\mathbf{Z}\right) = \frac{\mathbb{C}\_0}{\mathbf{L}^2} \frac{\partial^2 \Gamma}{\partial \mathbf{Z}^2} \sqrt{a\_{22}\left(\mathbf{Z}\right)}\tag{6.2.9}$$

Now we can write equation (6.1.1) in the following manner:

$$d\left(\mathbf{C}\_{0}\Gamma\right) = -\alpha\_{\mathbb{Z}}\left(\mathbf{Z}\right)dt + \beta\_{\mathbb{Z}}\left(\mathbf{Z}\right)do\left(t\right),\tag{6.2.10}$$

$$d\Gamma = \frac{-d\mathbf{Z}\left(\mathbf{Z}\right)}{\mathbf{C}\_0}dt + \frac{\beta\_\mathbf{z}\left(\mathbf{Z}\right)}{\mathbf{C}\_0}d\phi(t)\,. \tag{6.2.11}$$

where

$$\frac{\partial \sigma\_{\mathbf{z}}\left(\mathbf{Z}\right)}{\mathbf{C}\_{0}} = \Gamma F\_{\mathbf{z},0} + \frac{1}{\mathbf{L}} \frac{\partial \Gamma}{\partial \mathbf{Z}} F\_{\mathbf{z},1} + \frac{1}{\mathbf{L}^{2}} \frac{\partial^{2} \Gamma}{\partial \mathbf{Z}^{2}} F\_{\mathbf{z},2} \,\prime \,\tag{6.2.12}$$

$$\frac{\beta\_{\mathcal{Z}}(Z)}{\mathcal{C}\_{0}} = \left\{ \left( \Gamma^{2} a\_{00}(Z) \right) + \frac{1}{L^{2}} \left( \frac{\partial \Gamma}{\partial \mathcal{Z}} \right)^{2} a\_{11}(Z) + \frac{1}{L^{4}} \left( \frac{\partial^{2} \Gamma}{\partial \mathcal{Z}^{2}} \right)^{2} a\_{22}(Z) \right\}^{Y\_{2}} \tag{6.2.13}$$

Therefore, the Langevin form of the generalized SSTM is given by

$$d\Gamma = -\alpha\_{\mathbb{Z}}dt + \beta\_{\mathbb{Z}}do(t), \, 0 \le Z \le 1. \tag{6.2.14}$$

where 0 *Z <sup>Z</sup> C* and 0 *Z Z C* . 

Using equation (6.2.14) we can compute the time course of the dimensionless concentration for any given *L* .

The dimensionless/concentration, , varies from 0 to 1. *C t <sup>Z</sup>* is proportioned to the number of solute moles within a unit volume of porous/water matrix, and *C*0 is proportional to the maximum possible number of solute moles within the same matrix. Therefore, *Ct C <sup>Z</sup>* / 0 can be interpreted as the likelihood (probability) of finding solute moles within the matrix.

It should be noted that time, *t* , is not a dimensionless quantity and therefore, equation (6.2.14) is partially dimensionless equation. We will explore the dispersivity using equation (6.2.14) first before discussing a completely dimensionless equation.

(a)

(b)

(a)

*<sup>Z</sup>* at 0.5 *Z* when *L m* 1 , 0.1 *b* and

*<sup>Z</sup>* at 0.5 *Z* when *L m* 1 ,

Figure 6.1. (a) Realizations of

*b* 0.1 and *V m day* 0.5 / for different <sup>2</sup>

*V m day* 0.5 / for different <sup>2</sup>

values; (b) Realizations of

Multiscale, Generalised Stochastic Solute Transport Model in One Dimension 183

values.

### **6.3 Computational Exploration of the Langevin form of SSTM**

Equation (6.2.14) is not only an expression of how the solute disperses within a porous media but also an expression of nature of dispersion. Being a SDE, the drift coefficient *<sup>Z</sup>* portrays the dispersion due to the convective forces and the diffusive coefficient shows the dynamical behaviour of hydrodynamic dispersion. As *Z* has the range from 0 to 1 in equation (6.2.14), we can compute *<sup>Z</sup>* and values for a specific *Z* value and examine how they change over time. (We use 0 *C* 1.0 for computations, and therefore, *Z Z* and .) We have developed a finite difference algorithm to compute *<sup>Z</sup>* and adhering to the Ito integration as we have done before. Figure 6.1a and 6.1b show the time courses of *<sup>Z</sup>* and at 0.5 *Z* , respectively, for different <sup>2</sup> values when *L m* 1 (All times are given in days and *b* 0.1 . At low <sup>2</sup> values, *<sup>Z</sup>* behaves almost as a smooth deterministic function but at high <sup>2</sup> values it shows irregular behaviours. In these calculations, we have kept the mean velocity *V* at a constant value (0.5), therefore only fluctuating component affecting *<sup>Z</sup>* function is the solute concentration and its spatial derivatives. Further, Figure 6.1a and 6.1b only show a single realization for each <sup>2</sup> values. When we explore multiple realizations (not shown here), we see that randomness of *<sup>Z</sup>* and increases with higher <sup>2</sup> . One distinct feature of Figure 6.1b for *<sup>Z</sup>* is that is almost negligible for very small values of <sup>2</sup> but increases quite sharply for higher <sup>2</sup> values. *<sup>Z</sup>* does not behave in this manner. However, we can not ignore the effect of <sup>2</sup> at low values in computing *Z* , which has a follow-on affect on subsequent calculation. In other words, the affects of porous media, which <sup>2</sup> and the covariance kernel signify, can not be ignored as they affect the flow velocities significantly in making them stochastic. Figure 6.2a and 6.2b show *<sup>Z</sup>* and realization at 0.5 *Z* when 5 *L m* . The behaviours of *<sup>Z</sup>* and realizations are similar to those shown in Figures 6.1a and 6.1b. Figure 6.3a and Figure 6.3b show the similar trends for 10 *L m* . It should be noted that as *L* is increased, the time duration for the numerical solution of equation (6.2.14) should be increased. For example, when 10 *L m* , the model was run for 25 days to obtain Figures 6.3a and 6.3b. However, the order of magnitude for *<sup>Z</sup>* and *<sup>Z</sup>* has not changed as we change *L* in an order of magnitude.

*<sup>Z</sup>* and

adhering to the Ito integration as we have

*<sup>Z</sup>* and

values when *L m* 1 (All times are given in days and

*<sup>Z</sup>* behaves almost as a smooth deterministic function but

*<sup>Z</sup>* and

.) We have developed a finite

values. When we explore multiple

values.

and the covariance kernel signify, can not

*m*

*<sup>Z</sup>* and

realization at 0.5 *Z* when 5 *L m* . The

realizations are similar to those shown in Figures 6.1a and

values for a

at 0.5 *Z* ,

increases with

at low values in

*<sup>Z</sup>* has not

*<sup>Z</sup>* does not

is almost negligible for

182 in Porous Media - An Approach Based on Stochastic Calculus

Equation (6.2.14) is not only an expression of how the solute disperses within a porous media but also an expression of nature of dispersion. Being a SDE, the drift coefficient

*<sup>Z</sup>* portrays the dispersion due to the convective forces and the diffusive

specific *Z* value and examine how they change over time. (We use 0 *C* 1.0 for

 

values it shows irregular behaviours. In these calculations, we have kept the

*<sup>Z</sup>* is that

and

mean velocity *V* at a constant value (0.5), therefore only fluctuating component affecting

*<sup>Z</sup>* function is the solute concentration and its spatial derivatives. Further, Figure 6.1a

but increases quite sharply for higher <sup>2</sup>

computing *Z* , which has a follow-on affect on subsequent calculation. In other

be ignored as they affect the flow velocities significantly in making them stochastic.

6.1b. Figure 6.3a and Figure 6.3b show the similar trends for 10 *L m* . It should be noted that as *L* is increased, the time duration for the numerical solution of equation (6.2.14) should be increased. For example, when 10 *L m* , the model was run for 25 days to

shows the dynamical behaviour of hydrodynamic dispersion. As *Z* has

**6.3 Computational Exploration of the Langevin form of SSTM** 

the range from 0 to 1 in equation (6.2.14), we can compute

 *Z Z* 

done before. Figure 6.1a and 6.1b show the time courses of

realizations (not shown here), we see that randomness of

behave in this manner. However, we can not ignore the effect of <sup>2</sup>

*<sup>Z</sup>* and

obtain Figures 6.3a and 6.3b. However, the order of magnitude for

changed as we change *L* in an order of magnitude.

. One distinct feature of Figure 6.1b for

values,

and 6.1b only show a single realization for each <sup>2</sup>

words, the affects of porous media, which <sup>2</sup>

*<sup>Z</sup>* and

coefficient

computations, and therefore,

respectively, for different <sup>2</sup>

*b* 0.1 . At low <sup>2</sup>

very small values of <sup>2</sup>

behaviours of

Figure 6.2a and 6.2b show

*<sup>Z</sup>* and

at high <sup>2</sup> 

higher <sup>2</sup> 

difference algorithm to compute

Figure 6.1. (a) Realizations of *<sup>Z</sup>* at 0.5 *Z* when *L m* 1 , 0.1 *b* and *V m day* 0.5 / for different <sup>2</sup> values; (b) Realizations of *<sup>Z</sup>* at 0.5 *Z* when *L m* 1 , *b* 0.1 and *V m day* 0.5 / for different <sup>2</sup> values.

(a)

Multiscale, Generalised Stochastic Solute Transport Model in One Dimension 185

(b)

(c)

(d)

Figure 6.4. Realizations of *Z* at 0.5 *Z* when *b* = 0.1, *V* =0.5 for (a) *L* 1, (b) *L* 5,

(c) *L* 10 and (d) *L* 100

Figure 6.2. (a) Realizations of *<sup>Z</sup>* at 0.5 *Z* when *L m* 5 , 0.1 *b* and *V m day* 0.5 / for different <sup>2</sup> values; (b) Realizations of *<sup>Z</sup>* at 0.5 *Z* when *L m* 5 , *b* 0.1 and *V m day* 0.5 / for different <sup>2</sup> values.

Figure 6.3. (a) Realizations of *<sup>Z</sup>* at 0.5 *Z* when *L* 10 m, *b* = 0.1 and *V* =0.5 m/day for different <sup>2</sup> values; (b) Realizations of *<sup>Z</sup>* at 0.5 *Z* when *L* 10 m, *b* = 0.1 and *V* =0.5 m/day for different <sup>2</sup> values.

*<sup>Z</sup>* at 0.5 *Z* when *L m* 5 , 0.1 *b* and

*<sup>Z</sup>* at 0.5 *Z* when *L m* 5 ,

(b)

*<sup>Z</sup>* at 0.5 *Z* when *L* 10 m, *b* = 0.1 and *V* =0.5 m/day

*<sup>Z</sup>* at 0.5 *Z* when *L* 10 m, *b* = 0.1 and

184 in Porous Media - An Approach Based on Stochastic Calculus

(b)

values.

values; (b) Realizations of

(a)

Figure 6.2. (a) Realizations of

*b* 0.1 and *V m day* 0.5 / for different <sup>2</sup>

*V m day* 0.5 / for different <sup>2</sup>

Figure 6.3. (a) Realizations of

*V* =0.5 m/day for different <sup>2</sup>

for different <sup>2</sup>

values.

values; (b) Realizations of

Figure 6.4. Realizations of *Z* at 0.5 *Z* when *b* = 0.1, *V* =0.5 for (a) *L* 1, (b) *L* 5, (c) *L* 10 and (d) *L* 100

Multiscale, Generalised Stochastic Solute Transport Model in One Dimension 187

Figure 6.4a, 6.4b, 6.4c, and 6.4d show the realization of *Z* at 0.5 *Z* when *L* 1,5,10, and 100, respectively, for different values of <sup>2</sup> . (For all the calculations, we have used *b* = 0.1). When *L* 100 m, we computed *Z* values for 175 days and the affects of <sup>2</sup> on *Z* is quite dramatic, and this shows that equation (6.2.14) can display very complex behaviour patterns albeit its simplicity. It should be noted however that <sup>2</sup> plays major role in delimiting the nature of realizations; <sup>2</sup> values high than 0.25 in these situations produces highly irregular concentration realizations which could occur in highly heterogeneous porous formations such as fractured formations. that a

### **6.4 Dispersivities Based on the Langevin Form of SSTM for** *L* 10 **m**

One of the advantages of the partially dimensionless Langevin equation for the SSTM (equation 6.2.14) is that we can use it to compute the solute concentration profiles when the

travel length *L* is large. Equation (6.2.14) allows us to compute the dispersitivities using the stochastic inverse method (SIM) by estimating dispersivity for each realization of *Z* . For the SIM, we need to modify the deterministic-advection and dispersion equation into a partially dimensionless one. We start with the deterministic advection-dispersion equation with additive Gaussian noise,

$$\frac{\partial \mathbf{C}}{\partial t} = D\_{\perp} \frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2} - V\_{\times} \frac{\partial \mathbf{C}}{\partial \mathbf{x}} + \xi(\mathbf{x}, t), \tag{6.4.1}$$

where *DL* is the dispersion coefficient (dispersivity *Vx* ).

The partially dimensionless form of equation (6.4.1) is,

$$\frac{\partial \Gamma}{\partial t} = \frac{D\_L}{\mathcal{L}^2} \frac{\partial^2 \Gamma}{\partial \mathcal{Z}^2} - \frac{V\_\chi}{\mathcal{L}} \frac{\partial \Gamma}{\partial \mathcal{Z}} + \xi(\mathcal{Z}, t), \tag{6.4.2}$$

where 2 *DL L* is now estimated using SIM when *Vx* is known. Then the dispersivity value is (estimated 2 *DL <sup>L</sup>* ) <sup>2</sup> / . *L Vx*

Figure 6.5 show the scatter plots of dispersivity values estimated using the SIM for 1,5, *L* and 10 m. Each plot in Figure 6.5 gives 30 estimates of the dispersivity for a given value <sup>2</sup> . *Z* realizations were computed at 0.5 *Z* and 0.1 *b* for all plots. Table 6.1 summarizes the results giving the mean of each plot. We will compare these results with available data for dispersivities later in this chapter.

. (For all the calculations, we

values high than 0.25 in these

.

186 in Porous Media - An Approach Based on Stochastic Calculus

Figure 6.4a, 6.4b, 6.4c, and 6.4d show the realization of *Z* at 0.5 *Z* when

have used *b* = 0.1). When *L* 100 m, we computed *Z* values for 175 days and the

very complex behaviour patterns albeit its simplicity. It should be noted however that <sup>2</sup>

situations produces highly irregular concentration realizations which could occur in highly

One of the advantages of the partially dimensionless Langevin equation for the SSTM (equation 6.2.14) is that we can use it to compute the solute concentration profiles when the

travel length *L* is large. Equation (6.2.14) allows us to compute the dispersitivities using the stochastic inverse method (SIM) by estimating dispersivity for each realization of *Z* . For the SIM, we need to modify the deterministic-advection and dispersion equation into a partially dimensionless one. We start with the deterministic advection-dispersion equation

> 

*t L Z LZ*

Figure 6.5 show the scatter plots of dispersivity values estimated using the SIM for 1,5, *L* and 10 m. Each plot in Figure 6.5 gives 30 estimates of the dispersivity for a given value <sup>2</sup>

*Z* realizations were computed at 0.5 *Z* and 0.1 *b* for all plots. Table 6.1 summarizes the results giving the mean of each plot. We will compare these results with

 

 <sup>2</sup> <sup>2</sup> , , *<sup>L</sup> <sup>x</sup> C CC D V xt txx*

> <sup>2</sup> 2 2 , , *D V <sup>L</sup> <sup>x</sup> Z t*

*L* is now estimated using SIM when *Vx* is known. Then the dispersivity value is

(6.4.1)

the

(6.4.2)

on *Z* is quite dramatic, and this shows that equation (6.2.14) can display

*L* 1,5,10, and 100, respectively, for different values of <sup>2</sup>

plays major role in delimiting the nature of realizations; <sup>2</sup>

where *DL* is the dispersion coefficient (dispersivity *Vx* ).

The partially dimensionless form of equation (6.4.1) is,

available data for dispersivities later in this chapter.

heterogeneous porous formations such as fractured formations.

**6.4 Dispersivities Based on the Langevin Form of SSTM for** *L* 10 **m**

affects of <sup>2</sup>

with additive Gaussian noise,

where 2 *DL*

(estimated 2

*DL*

*<sup>L</sup>* ) <sup>2</sup> / . *L Vx*

As mentioned previously, the partially dimensionless equation (6.2.14) still requires us to compute for a large number of days when *L* is large. While the computational times are still manageable, we would like to develop a completely dimensionless Langevin equation for the SSTM. This could be especially useful and insightful when the mean velocity *V*

Multiscale, Generalised Stochastic Solute Transport Model in One Dimension 189

 ,., *<sup>t</sup> V Zt L*

where, *V Zt* , is mean velocity when 0 1, *Z* (m/day); *L* is travel length, m; and *t*

 . *Z dt d t* 

 *Z z*

and the variance of *<sup>L</sup> dt t*

, *<sup>Z</sup>*

 

*<sup>L</sup> Z d dt V* 

 , , *Z dd Z* 

, . *<sup>Z</sup>*

for the Ito integration to be accurate; therefore, we should have for maximum

*L V*

*z*

(6.5.1)

*V*

.

*<sup>z</sup>* (6.5.3)

variance, and

5 10 .

appropriately. Ideally

(6.5.2)

*V* 

<sup>4</sup> 10 0.5 , <sup>1000</sup>

and the range of

(6.5.4)

i.e, <sup>8</sup>

This allows us to

, as,

Therefore, if *V* 0.5 , 100 *L* and 0 200, *t* then, 0 1.0.

The completely dimentionless Langevin form of the SSTM is therefore,

are the Wiener increment with zero-mean and *<sup>L</sup> <sup>d</sup>*

could be considered as a constant.

We introduce dimensionless time,

compute *Z* realization for larger times.

*V* 

Equation (6.2.14) can be written as,

We can now change , *<sup>L</sup> dt d*

**6.5 Dimensionless Time** 

is time in days.

Therefore,

 ~ 0, *<sup>L</sup> <sup>d</sup> V* .

*<sup>L</sup>* . Suppose 0.5, 1000, *V L* then

Figure 6.5. Dispersivities estimated from SSTM for 1,5, *L* and 10.


As mentioned previously, the partially dimensionless equation (6.2.14) still requires us to compute for a large number of days when *L* is large. While the computational times are still manageable, we would like to develop a completely dimensionless Langevin equation for the SSTM. This could be especially useful and insightful when the mean velocity *V* could be considered as a constant.

### **6.5 Dimensionless Time**

Computational Modelling of Multi-Scale Non-Fickian Dispersion

188 in Porous Media - An Approach Based on Stochastic Calculus

Figure 6.5. Dispersivities estimated from SSTM for 1,5, *L* and 10.

Table 6.1. Mean dispersivities for the data in Figure 6.5.

Dispersivity

L=1 L=5 L=10

0.0001 0.050064 0.050112 0.0495545 0.001 0.0501232 0.05082 0.0511587 0.01 0.0520638 0.0604215 0.0778382 0.1 0.0669766 0.0832735 0.11899 0.25 0.0723413 0.111142 0.253195 0.4 0.0783754 0.142422 0.354335 0.6 0.0843219 0.170975 0.427603 0.8 0.0962623 0.225344 0.549473 1 0.110849 0.256348 0.609508

2 

We introduce dimensionless time, , as,

$$
\theta = \overline{V}(Z, \mathbf{t}) . \frac{\mathbf{t}}{L}' \tag{6.5.1}
$$

where, *V Zt* , is mean velocity when 0 1, *Z* (m/day); *L* is travel length, m; and *t* is time in days.

Therefore, if *V* 0.5 , 100 *L* and 0 200, *t* then, 0 1.0. This allows us to compute *Z* realization for larger times.

Equation (6.2.14) can be written as,

$$
\Gamma \left( Z \right) = -\alpha\_2 dt + \beta\_z d\alpha(t).
$$

We can now change , *<sup>L</sup> dt d V* and the variance of *<sup>L</sup> dt t V* .

Therefore,

$$
\Gamma\left(Z\right) = \frac{-\alpha\_2 L}{\overline{V}} d\theta + \beta\_z d\alpha(t) \,\tag{6.5.2}
$$

where *dw* ~ 0, *<sup>L</sup> <sup>d</sup> V* . ~

The completely dimentionless Langevin form of the SSTM is therefore,

$$
\Gamma(Z) = \alpha\_{z,\rho} d\theta + \beta\_z d\alpha(\theta),
\tag{6.5.3}
$$

where *d* are the Wiener increment with zero-mean and *<sup>L</sup> <sup>d</sup> V* variance, and

$$
\alpha\_{\mathcal{Z},\theta} = \frac{-\alpha\_{\mathcal{Z}}L}{\overline{V}}.\tag{6.5.4}
$$

To use equation (6.5.3), we need to choose and the range of appropriately. Ideally 0.0001 *<sup>L</sup> <sup>d</sup> V* for the Ito integration to be accurate; therefore, we should have for maximum as 0.0001*V <sup>L</sup>* . Suppose 0.5, 1000, *V L* then <sup>4</sup> 10 0.5 , <sup>1000</sup> i.e, <sup>8</sup> 5 10 . should 1000,

and Log10( <sup>2</sup>

0.94, respectively.

the boundary condition A

models for these significant relationships:

) for the most parts of the Log10 (Dispersivity) surface. Figure 6.7 shows the

, and (6.6.1)

*DCL <sup>s</sup>* , (6.6.2)

, and

(0.0001) and lower

) for different values of *L.*

linear relationship of Log10 (Dispersivity) vs Log10 (*L*) for different values of <sup>2</sup>

values of *L* (1 and 5). Therefore, we develop the following statistical nonlinear regression

Multiscale, Generalised Stochastic Solute Transport Model in One Dimension 191

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

where *Ds* is the dispersivity, and *C* <sup>1</sup> and *C* 2 are given in Tables 6.5 and 6.6 , respectively, along with *m1* and *m2* values. R-square values for equations (6.6.1) and (6.6.2) are 0.96 and

0.0001 0.0498 0.0500 0.0497 0.0498 0.0507 0.0686 0.001 0.0498 0.0499 0.0495 0.0477 0.0639 0.4982 0.01 0.0492 0.0510 0.0511 0.1642 0.5073 4.0672 0.1 0.0449 0.0592 0.1372 0.9309 2.9601 28.6151 0.25 0.0451 0.1123 0.2391 2.5441 6.1225 40.5301 0.4 0.0573 0.1340 0.3413 3.4365 8.1834 48.7567 0.6 0.0784 0.1824 0.4619 4.9440 10.9837 64.7589 0.8 0.0958 0.1987 0.7057 6.6800 14.9122 82.4423 1 0.1247 0.2159 0.8102 8.9878 19.9003 112.5246 L=1000 L=2000 L=4000 L=6000 L=8000 L=10000 0.0001 0.2697 0.7964 2.0630 4.1138 5.9939 8.1065 0.001 2.5154 7.2616 20.5460 32.6517 45.9978 69.0446 0.01 12.6500 30.0361 81.0270 155.2103 231.3154 324.3036 0.1 70.0564 156.8923 333.6665 523.4295 708.0212 903.6889 0.25 87.8303 185.7131 381.1019 569.9892 766.8287 978.5914 0.4 101.1441 203.0552 425.5467 625.6709 866.0189 1061.9651 0.6 131.0882 259.4990 528.0956 828.7496 1079.0040 1355.8468 0.8 173.1833 344.9935 691.7747 1070.9582 1399.5126 1771.3449 1 227.3204 453.3977 925.0844 1396.8663 1864.1378 2337.5588 Table 6.3. Longitudinal dispersivities (mean) for the range of *L* from 1 m to <sup>4</sup> 10 m under

*m*

 <sup>2</sup> 2

*m*

L=1 L=5 L=10 L=50 L=100 L=500

Figure 6.8 shows the same for Log10 (Dispersivity) vs Log10( <sup>2</sup>

σ2 Dispersivity

0.0499 0.0592 0.1340 0.1824 0.1987 7.2616

The gradient of the graphs are the same except for lower values of <sup>2</sup>

*D C <sup>s</sup>*

As we can see we may not gain much computational advantage with a completely dimensionless Langevin form of the SSTM.

### **6.6 Estimation of Field Scale Dispersivities**

We have estimated the longitudinal dispersivities using SSTM for two different boundary conditions:

$$\text{(A)}\quad \Gamma\_Z = 1 \quad \text{at} \quad Z = 0 \text{ and for } \quad t \ge 0; \quad \text{and} \quad \bar{\lambda}$$

(B) 1 *<sup>Z</sup>* at 0 *Z* and for 0 , *<sup>R</sup> t t* ; and 0 *<sup>Z</sup>* at *Z*=0 for *<sup>R</sup> t t* .

*Rt* is taken to be 1/3 of the total time (*T*) of the computational experiment. Table 6.1 and 6.2 show the dispersivity values for the boundary conditions A and B, respectively, when *L* 10 m based on 100 realisations for each of the boundary condition.


Table 6.2. Longitudinal dispersivities (mean) for the boundary condition A.

The values in Table 6.1 and 6.2 are similar for the similar values of <sup>2</sup> and *L* showing that (1) the SSTM procedure is robust in evaluating the dispersivities, and (2) the computed mean dispersivities do not depend on the boundary conditions, A and B. In these calculations, we have 0.5 *VZ* m/day.

We have also computed the dispersivities for larger scales up to 10,000 m, and Table 6.3 gives the mean values for the range of *L* from 1 m to <sup>4</sup> 10 m under the boundary condition A, and Table 6.4 gives the mean values for the range of *L* from 1 m to <sup>8</sup> 10 m for the boundary condition B. All mean values are calculated based on different sets of 100 realisations for each boundary condition. Except for the smallest <sup>2</sup> values (0.0001 and 0.001), the dispersivities have similar mean values for both boundary conditions, A and B. Therefore, it is quite reasonable to compute the dispersivities only for the boundary condition A for larger values of *L*. We can also hypothesise that the dispersivities are independent of the boundary conditions used to solve the SSTM. We have tested the SSTM for different values of *Rt* > (1/3) *T* when *L*>10 m. Figure 6.6 depicts the dispersivity plotted against <sup>2</sup> and *L* in Log10 scale, and Log10 (Dispersivity) is a linear function of Log10(*L*)

190 in Porous Media - An Approach Based on Stochastic Calculus

As we can see we may not gain much computational advantage with a completely

We have estimated the longitudinal dispersivities using SSTM for two different boundary

*Rt* is taken to be 1/3 of the total time (*T*) of the computational experiment. Table 6.1 and 6.2 show the dispersivity values for the boundary conditions A and B, respectively, when

0.0001 0.050013 0.050013 0.049828 0.001 0.050035 0.050223 0.050226 0.01 0.050646 0.055152 0.06112 0.1 0.055176 0.079403 0.136904 0.25 0.068846 0.108899 0.257902 0.4 0.083342 0.16346 0.333472 0.6 0.093185 0.191919 0.334818 0.8 0.109335 0.251033 0.54346 1 0.129395 0.331389 0.613823

(1) the SSTM procedure is robust in evaluating the dispersivities, and (2) the computed mean dispersivities do not depend on the boundary conditions, A and B. In these

We have also computed the dispersivities for larger scales up to 10,000 m, and Table 6.3 gives the mean values for the range of *L* from 1 m to <sup>4</sup> 10 m under the boundary condition A, and Table 6.4 gives the mean values for the range of *L* from 1 m to <sup>8</sup> 10 m for the boundary condition B. All mean values are calculated based on different sets of 100

0.001), the dispersivities have similar mean values for both boundary conditions, A and B. Therefore, it is quite reasonable to compute the dispersivities only for the boundary condition A for larger values of *L*. We can also hypothesise that the dispersivities are independent of the boundary conditions used to solve the SSTM. We have tested the SSTM for different values of *Rt* > (1/3) *T* when *L*>10 m. Figure 6.6 depicts the dispersivity plotted

and *L* in Log10 scale, and Log10 (Dispersivity) is a linear function of Log10(*L*)

Dispersivity L=1 L=5 L=10

and *L* showing that

values (0.0001 and

(B) 1 *<sup>Z</sup>* at 0 *Z* and for 0 , *<sup>R</sup> t t* ; and 0 *<sup>Z</sup>* at *Z*=0 for *<sup>R</sup> t t* .

*L* 10 m based on 100 realisations for each of the boundary condition.

Table 6.2. Longitudinal dispersivities (mean) for the boundary condition A.

The values in Table 6.1 and 6.2 are similar for the similar values of <sup>2</sup>

realisations for each boundary condition. Except for the smallest <sup>2</sup>

dimensionless Langevin form of the SSTM.

(A) 1 *<sup>Z</sup>* at 0 *Z* and for 0; *t* and

calculations, we have 0.5 *VZ* m/day.

against <sup>2</sup> 

conditions:

2 

**6.6 Estimation of Field Scale Dispersivities** 

and Log10( <sup>2</sup> ) for the most parts of the Log10 (Dispersivity) surface. Figure 6.7 shows the linear relationship of Log10 (Dispersivity) vs Log10 (*L*) for different values of <sup>2</sup> , and Figure 6.8 shows the same for Log10 (Dispersivity) vs Log10( <sup>2</sup> ) for different values of *L.* The gradient of the graphs are the same except for lower values of <sup>2</sup> (0.0001) and lower values of *L* (1 and 5). Therefore, we develop the following statistical nonlinear regression models for these significant relationships:

$$D\_s = \mathbb{C}\_{\,\_1} \left( \sigma^2 \right)^{m1} \text{, and } \tag{6.6.1}$$

$$D\_s = \mathbb{C}\_{\,\_2} \left( L \right)^{w2} \text{ \,\, \,} \tag{6.6.2}$$

where *Ds* is the dispersivity, and *C* <sup>1</sup> and *C* 2 are given in Tables 6.5 and 6.6 , respectively, along with *m1* and *m2* values. R-square values for equations (6.6.1) and (6.6.2) are 0.96 and 0.94, respectively.


Table 6.3. Longitudinal dispersivities (mean) for the range of *L* from 1 m to <sup>4</sup> 10 m under the boundary condition A

Figure 6.7. The linear relationship of Log10 (Dispersivity) vs Log10 (*L*) for different values of <sup>2</sup>

Multiscale, Generalised Stochastic Solute Transport Model in One Dimension 193

*L* (m) **1 5 10 50 100 500**  m1 0.039 0.125 0.311 0.605 0.677 0.704 C1 0.063 0.124 0.468 6.275 16.23 109.6 *L* (m) **1000 2000 4000 6000 8000 10000**  m1 0.690 0.642 0.605 0.578 0.567 0.552 C1 229.5 451.3 912.4 1368.7 1823.1 2281.4

) and Log10 (*L*)

Figure 6.8. The plot of Log10 (Dispersivity) vs Log10 ( <sup>2</sup>

Table 6.5. *m1* and *C*<sup>1</sup> values for different *L* for equation (6.6.1).

.


Table 6.4. Longitudinal dispersivities (mean) for the range of *L* from 1 m to <sup>8</sup> 10 m under the boundary condition B

Figure 6.6. The linear relationship of Log10 (Dispersivity) vs Log10 ( <sup>2</sup> ) for different values of L.

) for different values of L.

192 in Porous Media - An Approach Based on Stochastic Calculus

0.0001 0.0498 0.0500 0.0497 0.0498 0.0507 0.0686 0.1426 0.001 0.0498 0.0499 0.0495 0.0477 0.0639 0.4982 1.4690 0.01 0.0492 0.0510 0.0511 0.1642 0.5073 4.0672 12.0999 0.1 0.0449 0.0592 0.1372 0.9309 2.9601 28.6151 69.2489 0.25 0.0451 0.1123 0.2391 2.5441 6.1225 40.5301 87.0760 0.4 0.0573 0.1340 0.3413 3.4365 8.1834 48.7567 100.6075 0.6 0.0784 0.1824 0.4619 4.9440 10.9837 64.7589 132.1320 0.8 0.0958 0.1987 0.7057 6.6800 14.9122 82.4423 173.1823 1 0.1247 0.2159 0.8102 8.9878 19.9003 112.5246 221.6737

Table 6.4. Longitudinal dispersivities (mean) for the range of *L* from 1 m to <sup>8</sup> 10 m under

Figure 6.6. The linear relationship of Log10 (Dispersivity) vs Log10 ( <sup>2</sup>

Dispersivity L=1 L=5 L=10 L=50 L=100 L=500 L=1000

σ2

the boundary condition B

Figure 6.7. The linear relationship of Log10 (Dispersivity) vs Log10 (*L*) for different values of <sup>2</sup> .

Figure 6.8. The plot of Log10 (Dispersivity) vs Log10 ( <sup>2</sup> ) and Log10 (*L*)


Table 6.5. *m1* and *C*<sup>1</sup> values for different *L* for equation (6.6.1).

**The Stochastic Solute Transport** 

In Chapter 6, we developed the generalised Stochastic Solute Transport Model (SSTM) in 1 dimension and showed that it can model the hydrodynamic dispersion in porous media for the flow lengths ranging from 1 to 10000 m. For computational efficiency, we have employed one of the fastest converging kernels tested in Chapter 6 for illustrative purposes, but, in principle, the SSTM should provide scale independent behaviour for any other velocity covariance kernel. If the kernel is developed based on the field data, then the SSTM based on that particular kernel should give realistic outputs from the model for that particular porous medium. In the development of the SSTM, we assumed that the hydrodynamic dispersion is one dimensional but by its very nature, the dispersion lateral to

First, we solve the integral equation with the covariance kernel in two dimensions, and use the eigen values and functions thus obtained in developing the two dimensional stochastic solute transport model (SSTM2d). Then we solve the SSTM2d numerically using a finite difference scheme. In the last section of the chapter, we illustrate the behaviours of the

We consider the flow direction to be *x* and the coordinate perpendicular to *x* to be *y* in the 2 dimensional flow with in the porous matrix saturated with water. Then the distance

> 2 <sup>1122</sup> ( , , , ) exp *<sup>r</sup> qx y x y <sup>b</sup>*

2 12 12

( )( ) ( , , , ) exp ,

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

2 1 2 1 2

() () exp exp .

*x x y y b b*

*xx yy*

2

2 2

2 2

is the variance at a given point, i.e., when 1 2 *x x* and

between the points 1 1 (,) *x y* and 2 2 (,) *x y* , *r*, is given by 1/2 <sup>2</sup> <sup>2</sup>

the flow direction occurs. We intend to explore this aspect in this chapter.

SSTM2d graphically to show the robustness of the solution.

then define a velocity covariance kernel as follows:

1122

is a constant. <sup>2</sup>

1 2 *y y* . The covariance can be written as,

**7.2 Solving the Integral Equation** 

where <sup>2</sup> 

**7.1 Introduction** 

**Model in 2-Dimensions** 

12 12 ( )( ) *xx yy* . We can

, (7.2.1)

(7.2.2)


Table 6.6. *m2* and *C* <sup>2</sup> values for different <sup>2</sup> for equation (6.6.2).

Figure 6.9. Mean dispersivity from the SSTM and experimental dispersivity vs flow length (log10 scale).

We can estimate the approximate dispersivity values either from Figures 6.7 and 6.8, or from equations (6.6.1) and (6.6.2). It is quite logical to ask the question whether we can characterise the large scale aquifer dispersivities using a single value of <sup>2</sup> ? To answer this question, we resort to the published dispersivity values for aquifers. We use the dispersivity data first published by Gelhar et al. (1992) and reported to Batu (2006). We extracted the tracer tests data related to porous aquifers in 59 different locations characterised by different geologic materials. The longest flow length was less than 10000 m. We then plotted the experimental data and overlaid the plot with the dispersivity vs *L* curves from the SSTM for each <sup>2</sup> value. Figure 6.9 shows the plots, and <sup>2</sup> =0.1 best fit to the experimental data. In other words, by using one value of <sup>2</sup> , we can obtain the dispersivity for any length of the flow by using the SSTM. We can also assume that each experimental data point represents the mean dispersivity for any length of the flow by using the SSTM. We can also assume that each experimental data point represents the mean dispersivity at a particular flow length. If that is the case, Figure 6.9 can be interpreted as follows: by using the SSTM, we can obtain sufficiently large number of realisations for particular values of <sup>2</sup> and the mean flow velocity, and the mean values of the dispersivities estimated for those concentration realisations do represent the experimental dispersivities. <sup>2</sup> can be hypothesised to indicate the type of media (e.g. fractured, porous etc.). These findings support the hypothesis that the dimensionless SSTM is scale-independent, i.e., one value of <sup>2</sup> would be sufficient to characterise the dispersivity at different flow lengths. It is important to note that the role of the mean velocity in these calculations. We used 0.5 m/day to represent an indicative value in real aquifers, but the character of solutions do not change, if we assume a different value; only the specific values of <sup>2</sup> would be changed to represent a given flow situation. The
