**A Generalized Mathematical Model in One-Dimension**

### **4.1 Introduction**

Computational Modelling of Multi-Scale Non-Fickian Dispersion

*t*

*t* -

=0.01; *b=*0.01 m ) in the

*<sup>L</sup>* , using equation (3.10.5)

. (3.10.5)

*<sup>L</sup>* for *a* = 3000 meters, we

116 in Porous Media - An Approach Based on Stochastic Calculus

variance, which can be readly incorporated into the numerical solution scheme for SSTM. Because of the high computational times involved, instead of using the SIM procedure, we estimate the dispersivity for a computational experiment by limiting ourselves to the

> ( ,0) 0, 0 (0, ) 1.0, 0 ( , ) 0, 0

 

(Fetter, 1999) given below using nonlinear regression: ( Mathematica® was used for this

Equation (3.10.5) is the analytical solution for the one-dimensional advection-dispersion equation for the initial and boundary conditions given in (a) above. To get a reliable estimate

need to run the simulation of a domain of 4500 meters for 6000 days. However, to reduce the computational time, one can use higher *x* and *t* values than ideally suitable in solving SPDEs thereby sacrificing the reliability. ( *x* =0.01 m and *t* =0.00001 days would give very good solutions to SSTM.) However, this approximate procedure is only valid for the velocity

We overlay some representative values of dispersivities from the SSTM on a graph of the field measurements obtained by Gelhar (1986) from different experiments in Figure 3.34, which shows that SSTM could model the multi-scale dispersion with a single set of

In Figure 3.34, the larger (hollow) circles depict the most reliable experimental data whereas the smaller (filled) circles give the data with lesser experimental accuracy. The dotted lines show the bounded region of experimental data. The estimated dispersivity values (filled squares) from SSTM are within the bounded region and follow similar trends to the most reliable experimental data. We estimated the dispersivities from SSTM using only a limited number of computational experiments, and each computational experiment produces a random realisation. Therefore, the estimated dispersivities are stochastic quantities just like the experimental values, and it is reasonable to expect discrepancies within the bounded region.

=0.01; *b=*0.01 m, that would give rise to similar non-linear scale dependency

*<sup>L</sup>* for given *a*, we need to have *x-*axis length of 1.5 *a* meters and should have *C*(*x,t*)

*Cx x Ct t Ct t*

( , ) 0.5 exp <sup>2</sup> *<sup>L</sup>* <sup>2</sup> *<sup>L</sup> <sup>L</sup> a t <sup>a</sup> a t C x t erfc erfc t*

solution

covariance kernel used in this chapter. This procedure is not as reliable as the SIM.

kernel used the

of "deterministic" dispersivity evaluated using the realisations of SSTM.

SSTM parameters obtained for experimental sand aquifer (i.e., <sup>2</sup>

(b) Solve equation (21) assuming mean velocity to be 1.0 m/day, (c) Use a realization of *C*(*x,t*) at x=*a* to estimate the dispersivity,

realization upto 2 *a* days. For example, to obtain an estimate of

( ) is a Gaussian random variable with zero-mean and <sup>2</sup>

and

spatially averaged *d t*

following way:

purpose.)

parameters, <sup>2</sup>

of  (a) Set the initial and boundary conditions as,

In the previous chapter we derived a stochastic solute transport model (equation (3.2.14)); we developed the methods to estimate its parameters, and investigated its behaviour numerically. We see some promise to characterise the solute dispersion at different flow lengths, and there are some indications that equation (3.2.14) produce the behaviours that would be interpreted as capturing the scale-dependency of dispersivity. However, there are weaknesses in the model as evident from Chapter 3. These weaknesses, which are discussed in the next section, are stemming from the very assumptions we made in the development of the model. One could argue that by relaxing the Fickian assumptions, we are actually complicating the problem quite unnecessarily. But as we see in Chapter 3 and in this chapter, we develop a new mathematical and computational machinery at a more fundamental level for the hydrodynamic dispersion in saturated porous media.

We see that equation (3.2.14) is based on assuming a covariance kernel for the velocity fluctuations, and the solution is dependent on solving an integral equation (see equation (3.3.11)). In Chapter 3, the integral equation is solved analytically for the covariance kernel given by equation (3.3.10) to obtain the eigen values and eigen functions, but analytical solutions of integral equations can not be easily derived for any arbitrary covariance kernel. This limits the flexibility of the SSTM in employing a suitable covariance kernel independent of the ability to solve relevant integral equations. Further, we need to solve the SSTM in a much more computationally efficient manner, and estimating dispersivity by always relating to the deterministic advection-dispersion equation is not quite satisfactory. Therefore, we seek to develop a more general form of equation (3.2.14) in this chapter.

### **4.2 The Development of the Generalized Model**

We restate equation (3.2.14) in the differential form:

$$d\mathbf{C} = \mathbf{S}\left(\overline{V}\left(\mathbf{x},t\right)\mathbf{C}\left(\mathbf{x},t\right)\right)dt + \mathbf{S}\left(\mathbf{C}\left(\mathbf{x},t\right)d\boldsymbol{\beta}\_{w}\left(t\right)\right),\tag{4.2.1}$$

$$\text{where}\quad d\beta\_m(t) = \sigma \sum\_{j=1}^m \sqrt{\lambda\_j} f\_j db\_j(t). \tag{4.2.2}$$

We use the same notations and symbols as in Chapter 3. In equation (4.2.2), *d t <sup>m</sup>* is calculated by summing *m* terms of ( *jj j f db t* ), and for each eigen function, *<sup>j</sup> f* , there is an associated independent Wiener process increment ( *db t <sup>j</sup>* ).

 

*j j kj k*

*j j kj*

2

*g C x t g xC x t g C x t e*

*j j j kj*

Now the derivative within the summation in equation (4.2.6) is evaluated.

Substituting this derivation back to equation (4.2.6) and defining,

2

*j j ij kj kj kj*

*j*

*p*

*k*

*x*

, 2 , .

*C xt rC xt x s <sup>e</sup>*

*r xs*

*kj kj*

*x t C xt g gr x s e x x*

, , <sup>2</sup>

2

4 2 ,

0 1 2 1

*kj kj kj kj kj*

*C xt C xt C xt*

*x x x*

1 0 1

2

*j kj kj kj j j kj*

*g g r x s e C xt g g x g e*

*j*

*p*

0 1

*C xt g g x g e*

0 1 1

*C xt C xt*

0 1

*kj kj*

*j j j kj*

1

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

2

*k*

*j*

*xt C xt g g x g e*

2 ,

2

*j kj kj kj k p p*

*j*

 

*p*

2 2

*k k p*

*j j*

0 1

, , ,

 

,

*x*

function *j*,

0 1

*x*

2

Then taking the derivative of

*j*

*p*

*k*

*x*

*j*

*k*

*g g x g C xt g C xt e x x x*

 

2

, , , ,

*j*

*p*

*kj kj*

A Generalized Mathematical Model in One-Dimension 119

*r xs*

*kj kj kj kj kj kj*

*C xt C xt e C xt r x s e e x x*

*r xs r xs r xs kj kj*

*kj kj*

*r xs*

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

2

*kj kj*

*r xs*

*j*

*p*

*k*

 

2

*<sup>j</sup> x t*, with respect to x,

2 2

*kj kj*

*r xs*

2

*kj kj kj*

 <sup>2</sup> 2

*k*

*j*

 

*p*

*r xs r*

, <sup>2</sup>

<sup>2</sup>

*j*

*p*

*k*

*g g x gC xt g r x s C xt e*

*j j kj j kj kj kj*

*C xt g g x ge g gr x s e x*

*r xs r xs*

*kj kj kj kj*

2

*gr x s e g re C xt*

*kj kj*

*r xs*

2

, , , 2

2 2 2

2

, , , , .

2

*j*

*p*

*k*

2

2

*kj kj*

*r xs*

2

(4.2.6)

*kj kj*

*<sup>j</sup> x t*, for eigen

*kj kj*

, ,

*C xt x*

(4.2.7)

*r xs*

, .

*C xt x*

*kj*

*x s*

*r xs*

This Wiener process increment is -correlated in time only, and for a particular point in time, *<sup>i</sup> t* , the Wiener process *B f* , generates the Wiener increments *db t <sup>j</sup> <sup>i</sup>* for each of the eigen functions. At another time, *<sup>k</sup> t* , *B f* , generates *db t <sup>j</sup> <sup>k</sup>* for each of the eigen function. However, as the Wiener increments are Gaussian with *t* variance (see Chapter 2), i.e., *<sup>j</sup> <sup>i</sup> db t* and *db t <sup>j</sup> <sup>k</sup>* are essentially the same stochastic variable.

The number of terms that need to be summed up, *m*, has to be determined by considering the contribution each pair of eigen value *<sup>i</sup>* and corresponding eigen function *<sup>j</sup> f* make to approximate the covariance.

Instead of solving the integral equation to obtain eigen function, we assume the following function to be a generalized eigen function. (We will discuss the method to obtain this function for any given kernel in section 4.3.)

We define,

$$f\_{\uparrow}(\mathbf{x}) = \mathbf{g}\_{0\uparrow} + \mathbf{g}\_{1\uparrow}\mathbf{x} + \sum\_{k=2}^{p\_{\downarrow}} \mathbf{g}\_{k\downarrow} e^{-\imath\_{\downarrow} \left(\mathbf{x} - \imath\_{\downarrow}\right)^{2}},\tag{4.2.3}$$

where *fj x* is the *j*th eigen function of a given covariance kernel, 0 *<sup>j</sup> g* , <sup>1</sup> *<sup>j</sup> g* and *kj g* are numerical coefficients, *jk r* and *jk s* are the coefficients in the exponent where *k* is an integer index ranging from 2 to *<sup>j</sup> p* . *<sup>j</sup> p* is determined by the computational method in section 4.3.

The differential operator is defined by,

$$S = -\left(\frac{\hbar\_x}{\mathcal{D}} \frac{\partial^2}{\partial \boldsymbol{\alpha}^2} + \frac{\partial}{\partial \boldsymbol{\alpha}}\right) \boldsymbol{\nu}$$

and the second term on the right hand side of equation (4.2.1) can be written as,

$$S\left(\mathbf{C}\left(\mathbf{x},t\right)\sum\_{j=1}^{m}\sqrt{\lambda\_{j}}f\_{j}db\_{j}\left(t\right)\right)=\sigma\sum\_{j=1}^{m}\sqrt{\lambda\_{j}}S\left(\mathbf{C}\left(\mathbf{x},t\right)f\_{j}\right)db\_{j}\left(t\right),\tag{4.2.4}$$

after substituting equation (4.2.2) into the second terms of equation (4.2.1) and taking *<sup>j</sup>* and *db t <sup>j</sup>* out of the operand. We can now expand the differential,

$$\mathbf{S}\left(\mathbf{C}\left(\mathbf{x},t\right)f\_{/}\right) = \mathbf{S}\left(\mathbf{C}\left(\mathbf{x},t\right)\left(\mathbf{g}\_{0/} + \mathbf{g}\_{1/}\mathbf{x} + \sum\_{k=2}^{p\_{j}} \mathbf{g}\_{k\not\!\!/} \mathbf{e}^{-\mathbf{e}\_{\mathbb{Q}}\left(\mathbf{x}-\mathbf{e}\_{\mathbb{Q}}\right)^{2}}\right)\right). \tag{4.2.5}$$

Let us evaluate the derivatives separately;


, generates *db t <sup>j</sup> <sup>k</sup>* for each of the eigen

*<sup>i</sup>* and corresponding eigen function *<sup>j</sup> f* make to

, generates the Wiener increments *db t <sup>j</sup> <sup>i</sup>* for each of

118 in Porous Media - An Approach Based on Stochastic Calculus

function. However, as the Wiener increments are Gaussian with *t* variance (see Chapter

The number of terms that need to be summed up, *m*, has to be determined by considering

Instead of solving the integral equation to obtain eigen function, we assume the following function to be a generalized eigen function. (We will discuss the method to obtain this

<sup>2</sup>

*f x g gx ge*

2 *j*

*p*

*k*

2 <sup>2</sup> 2

after substituting equation (4.2.2) into the second terms of equation (4.2.1) and taking

 <sup>2</sup> 0 1

*SC xt f S C xt g g x g e*

*j j j kj*

*x x* 

*jj j j j j*

 

where *fj x* is the *j*th eigen function of a given covariance kernel, 0 *<sup>j</sup> g* , <sup>1</sup> *<sup>j</sup> g* and *kj g* are numerical coefficients, *jk r* and *jk s* are the coefficients in the exponent where *k* is an integer index ranging from 2 to *<sup>j</sup> p* . *<sup>j</sup> p* is determined by the computational method in

*kj kj*

, (4.2.3)

*r xs*

,

, (4.2.4)

2

*kj kj*

. (4.2.5)

*r xs*

*j*

*p*

*k*

*S C x t f db t*

0 1

*j j j kj*

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

and the second term on the right hand side of equation (4.2.1) can be written as,

1 1 , , *m m*

*<sup>j</sup>* and *db t <sup>j</sup>* out of the operand. We can now expand the differential,

, ,

*j j*

*S C x t f db t* 

2), i.e., *<sup>j</sup> <sup>i</sup> db t* and *db t <sup>j</sup> <sup>k</sup>* are essentially the same stochastic variable.

This Wiener process increment is

time, *<sup>i</sup> t* , the Wiener process *B f*

the contribution each pair of eigen value

function for any given kernel in section 4.3.)

The differential operator is defined by,

Let us evaluate the derivatives separately;

approximate the covariance.

We define,

section 4.3.

the eigen functions. At another time, *<sup>k</sup> t* , *B f*

$$\begin{split} &\frac{\partial}{\partial \mathbf{x}} \Big( \mathbf{C}(\mathbf{x},t) \Big( \mathbf{g}\_{0j} + \mathbf{g}\_{1j}\mathbf{x} + \sum\_{k=2}^{p\_{j}} \mathbf{g}\_{k} e^{-\mathbf{e}\_{0}\cdot(\mathbf{x}-\mathbf{e}\_{0})} \Big)^{2} \Big) \\ &= \frac{\partial}{\partial \mathbf{x}} \Big[ \mathbf{g}\_{0j}\mathbf{C}(\mathbf{x},t) + \mathbf{g}\_{1j}\mathbf{x}\mathbf{C}(\mathbf{x},t) + \sum\_{k=2}^{p\_{j}} \mathbf{g}\_{kj}\mathbf{C}(\mathbf{x},t) e^{-\mathbf{e}\_{0}\cdot(\mathbf{x}-\mathbf{e}\_{0})} \Big] \\ &= \mathbf{g}\_{0j}\frac{\partial \mathbf{C}(\mathbf{x},t)}{\partial \mathbf{x}} + \mathbf{g}\_{1j}\mathbf{x} \frac{\partial \mathbf{C}(\mathbf{x},t)}{\partial \mathbf{x}} + \mathbf{g}\_{1j}\mathbf{C}(\mathbf{x},t) + \sum\_{k=2}^{p\_{j}} \mathbf{g}\_{kj} \frac{\partial}{\partial \mathbf{x}} \Big( \mathbf{C}(\mathbf{x},t) e^{-\mathbf{e}\_{0}\cdot(\mathbf{x}-\mathbf{e}\_{0})^{2}} \Big). \end{split} \tag{4.2.6}$$

Now the derivative within the summation in equation (4.2.6) is evaluated.

$$\begin{split} &\frac{\partial}{\partial\mathbf{x}}\bigg(\mathsf{C}\left(\mathbf{x},t\right)e^{-\imath\_{\boldsymbol{\eta}}\left(\mathbf{x}-\mathbf{s}\_{\boldsymbol{\eta}}\right)^{2}}\bigg) = \mathsf{C}\left(\mathbf{x},t\right)\bigg(-2\mathsf{r}\_{\boldsymbol{\eta}}\left(\mathbf{x}-\mathbf{s}\_{\boldsymbol{\eta}}\right)e^{-\imath\_{\boldsymbol{\eta}}\left(\mathbf{x}-\mathbf{s}\_{\boldsymbol{\eta}}\right)^{2}}\bigg) + e^{-\imath\_{\boldsymbol{\eta}}\left(\mathbf{x}-\mathbf{s}\_{\boldsymbol{\eta}}\right)^{2}}\frac{\partial\mathsf{C}\left(\mathbf{x},t\right)}{\partial\mathbf{x}} \\ &= \left[-2\mathsf{r}\_{\boldsymbol{\eta}}\mathsf{C}\left(\mathbf{x},t\right)\bigg(\mathbf{x}-\mathbf{s}\_{\boldsymbol{\eta}}\right) + \frac{\partial\mathsf{C}\left(\mathbf{x},t\right)}{\partial\mathbf{x}}\right]e^{-\imath\_{\boldsymbol{\eta}}\left(\mathbf{x}-\mathbf{s}\_{\boldsymbol{\eta}}\right)^{2}}.\end{split}$$

Substituting this derivation back to equation (4.2.6) and defining, *<sup>j</sup> x t*, for eigen function *j*,

 2 2 0 1 2 0 1 2 1 0 1 2 , , , , , , , <sup>2</sup> , , 2 , *j kj kj j kj kj j kj kj p r xs j j j kj k p r xs j j ij kj kj kj k p r xs j kj kj kj j j kj k xt C xt g g x g e x C xt C xt C xt g g x gC xt g r x s C xt e x x x g g r x s e C xt g g x g e* <sup>2</sup> 2 , . *j kj kj p r xs k C xt x* (4.2.7) *kjj* 

Then taking the derivative of *<sup>j</sup> x t*, with respect to x,

$$\begin{split} &\frac{\partial\mathcal{J}\boldsymbol{\nu}\left(\mathbf{x},t\right)}{\partial\mathbf{x}} = \left\{\mathcal{g}\_{1j} - 2\sum\_{k=2}^{p\_{j}} \mathcal{g}\_{kj}\boldsymbol{r}\_{kj}\left(\mathbf{x} - \mathbf{s}\_{kj}\right)e^{-\boldsymbol{\eta}\_{j}\left(\mathbf{x} - \mathbf{s}\_{kj}\right)^{2}}\right\} \frac{\partial\mathbf{C}\left(\mathbf{x},t\right)}{\partial\mathbf{x}} \\ &+ \left\{\mathbf{4}\sum\_{k=2}^{p\_{j}} \mathcal{g}\_{kj}\boldsymbol{r}\_{kj}^{2}\left(\mathbf{x} - \mathbf{s}\_{kj}\right)^{2}e^{-\boldsymbol{\eta}\_{j}\left(\mathbf{x} - \mathbf{s}\_{kj}\right)^{2}} - 2\sum\_{k=2}^{p\_{j}} \mathcal{g}\_{kj}\boldsymbol{r}\_{kj}e^{-\boldsymbol{\eta}\_{j}\left(\mathbf{x} - \mathbf{s}\_{kj}\right)^{2}}\right\} \mathbf{C}\left(\mathbf{x},t\right) \\ &+ \left\{\mathbf{g}\_{nj} + \mathbf{g}\_{1j}\mathbf{x} + \sum\_{k=2}^{p\_{j}} \mathcal{g}\_{kj}e^{-\boldsymbol{\eta}\_{j}\left(\mathbf{x} - \mathbf{s}\_{kj}\right)^{2}}\right\} \frac{\partial^{2}\mathbf{C}\left(\mathbf{x},t\right)}{\partial\mathbf{x}^{2}} + \left\{\mathbf{g}\_{1j} - 2\sum\_{k=2}^{p\_{j}} \mathcal{g}\_{kj}\boldsymbol{r}\_{kj}\left(\mathbf{x} - \mathbf{s}\_{kj}\right)e^{-\boldsymbol{\eta}\_{j}\left(\mathbf{x} - \mathbf{s}\_{kj}\right)^{2}}\right\} \frac{\partial\mathbf{C}\left(\mathbf{x},t\right)}{\partial\mathbf{x}},\end{split}$$

The expression of stochastic partial differential equation (equation (4.2.1)) as an Ito diffusion in time would give us the mathematical justification in solving it using Ito calculus. In our development, the eigen functions, *fj x* are continuous, differentiable functions so are the coefficients, *Pxi ij* 0,1,2 . Therefore, the Ito stochastic product rule is the same as the product rule in standard calculus (Klebaner, 1998), and we employ this fact in the previous derivations. Further, we assume that the mean velocity *V xt* , is a continuous, differentiable function which is a reasonable assumption given that the average velocity in aquifer situation is based on the hydraulic conductivity, porosity and the pressure gradient

A Generalized Mathematical Model in One-Dimension 121

across a large enough domain within a much larger total flow length. Therefore, we can write the drift term of equation (4.2.1) as follows:

, , . <sup>2</sup>

*x*

*h C xt V xt*

*x*

*x*

1

 

*j*

and then,

where,

*m*

2

2

2

*x*

2

By substituting equations (4.2.14) and (4.2.9) into equation (4.2.1),

, , <sup>2</sup>

*h C xt V xt*

2

0 1 2 2

, , , <sup>2</sup>

*x*

*x x*

*j j j j j*

2

*V xt C xt SV xtC xt C xt V xt*

, , , , , ,

2 2

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

*<sup>h</sup> V xt C xt V xt C xt C xt <sup>V</sup> x xx x*

2 2

 

*x x*

, , , ,

<sup>2</sup>

, , , , , *C xt C xt dC x t C x t dI dI dI*

 <sup>2</sup> 0 2 0

*V xt V xt <sup>h</sup> dI x t dt P db t*

*<sup>x</sup> <sup>x</sup> x x dC x t*

2

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

*x x*

2

*C xt C xt PC xt P <sup>P</sup> db t x x*

*x*

 

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

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

0 1 2 2

(4.2.15)

1

, (4.2.16)

 

*j*

*jj j*

*m*

*x x*

*V xt V xt <sup>h</sup> <sup>h</sup> V xt C xt C xt V xt*

2 2

(4.2.14)

*x x*

$$\begin{split} \mathbf{x} &= 2\left\{ \mathbf{g}\_{1j} - 2\sum\_{k=2}^{p\_j} \mathbf{g}\_{kj} r\_{lj} \left( \mathbf{x} - \mathbf{s}\_{lj} \right) e^{-\boldsymbol{\eta}\_{lj} \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2} \right\} \frac{\partial \mathbf{C}(\mathbf{x}, t)}{\partial \mathbf{x}} + \left\{ \mathbf{g}\_{0j} + \mathbf{g}\_{1j} \mathbf{x} + \sum\_{k=2}^{p\_j} \mathbf{g}\_{kj} e^{-\boldsymbol{\eta}\_{lj} \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2} \right\} \frac{\partial^2 \mathbf{C}(\mathbf{x}, t)}{\partial \mathbf{x}^2} \\ &+ \left\{ 4 \sum\_{k=2}^{p\_j} \mathbf{g}\_{kj} r\_{lj}^2 \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2 e^{-\boldsymbol{\eta}\_{lj} \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2} - 2 \sum\_{k=2}^{p\_j} \mathbf{g}\_{kj} r\_{lj} e^{-\boldsymbol{\eta}\_{lj} \left( \mathbf{x} - \mathbf{s}\_{lj} \right)^2} \right\} \mathbf{C}(\mathbf{x}, t) \end{split} \tag{4.2.8}$$

Therefore, equation (4.2.5) can be expressed as,

$$\begin{split}-S\left(\mathbf{C}\left(\mathbf{x},t\right)f\_{/}\right) &= \frac{h\_{\mathbf{x}}}{2} \left[\frac{\hat{\sigma}^{2}\left(\mathbf{C}\left(\mathbf{x},t\right)f\_{/}\right)}{\hat{\sigma}\mathbf{x}^{2}}\right] + \left[\frac{\partial\left(\mathbf{C}\left(\mathbf{x},t\right)f\_{/}\right)}{\hat{\sigma}\mathbf{x}}\right], \\ &= \nu\_{/}\left(\mathbf{x},t\right) + \frac{h\_{\mathbf{x}}}{2} \frac{\partial\left(\mathbf{w}\_{/}\left(\mathbf{x},t\right)\right)}{\partial\mathbf{x}}, \\ &= P\_{0/}\left(\mathbf{x}\right)\mathbf{C}\left(\mathbf{x},t\right) + P\_{1/}\frac{\partial\mathbf{C}\left(\mathbf{x},t\right)}{\partial\mathbf{x}} + P\_{2/}\frac{\partial^{2}\mathbf{C}\left(\mathbf{x},t\right)}{\partial\mathbf{x}^{2}}, \end{split} \tag{4.2.9}$$

Where

 2 2 2 0 2 <sup>2</sup> <sup>2</sup> 2 2 2 <sup>4</sup> <sup>2</sup> , <sup>2</sup> *j kj kj j j kj kj kj kj p r xs j ij kj kj kj k p p r xs r xs x kj kj kj kj kj k k P x g gr x s e h gr x s e gre* (4.2.10)

$$P\_{1}\left(\mathbf{x}\right) = \mathbf{g}\_{0\rangle} + \mathbf{g}\_{1\rangle}\mathbf{x} + \sum\_{k=2}^{p\_{\parallel}} \mathbf{g}\_{k\rangle} e^{-s\_{\parallel}\left(\mathbf{x} - s\_{\parallel}\right)^{2}} + \left(\frac{\hbar\_{\times}}{2}\right) \left[ 2\left(\mathbf{g}\_{\neq} - 2\sum\_{k=2}^{p\_{\parallel}} \mathbf{g}\_{k\rangle} r\_{kj} \left(\mathbf{x} - s\_{\neq}\right) e^{-s\_{\parallel}\left(\mathbf{x} - s\_{\neq}\right)^{2}}\right) \right],\tag{4.2.11}$$

And

$$P\_{2\_{\!/\!}}\left(\mathbf{x}\right) = \left(\frac{\hbar\_{\mathbf{x}}}{2}\right) \left[\mathbf{g}\_{0\!/\!} + \mathbf{g}\_{1\!/\!}\mathbf{x} + \sum\_{k=2}^{p\_{\!/\!}} \mathbf{g}\_{k\!/\!} \mathbf{e}^{-\mathbf{r}\_{\!/\!} \left(\mathbf{x} - \mathbf{s}\_{\!/\!}\right)^{2}}\right].\tag{4.2.12}$$

We see that, in equation (4.2.9), the coefficients, *Px i ij* 0,1,2 are functions of *x* only.

If we recall the premise on which stochastic calculus is discussed in Chapter 2, *C xt* , is a continuous, non-differentiable function, and equation (4.2.1) should be interpreted as an Ito stochastic differential, which essentially mean equation (4.2.1) has to be understood only in the following form:

$$\mathbf{C}\left(\mathbf{x},t\right) = \int A\left(\mathbf{C}\left(\mathbf{x},t\right),t,\mathbf{x}\right)dt + \sum\_{i} \int B\_{i}\left(\mathbf{C}\left(\mathbf{x},t\right),t,\mathbf{x}\right)dw\_{i}d\mathbf{w}\_{i} \tag{4.2.13}$$

where *A C xt tx* , ,, is the drift coefficient and *B C xt tx <sup>i</sup>* , ,, are diffusion coefficients of the Ito diffusion, equation (4.2.13); here *dw t <sup>i</sup>* are increments of the standard Wiener process. Ito diffusion are an interesting class of stochastic integrals and has many advantageous properties of practical importance (Klebaner, 1998).

*C xt C xt*

*x x*

2

2 2

2 2 , <sup>2</sup>

*kj kj*

(4.2.12)

(4.2.11)

*kj kj kj kj*

2

(4.2.9)

(4.2.10)

(4.2.8)

120 in Porous Media - An Approach Based on Stochastic Calculus

1 0 1 2

, , , ,

*x x*

*kj kj kj kj*

*r xs r xs x*

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

*C x t A C x t t x dt B C x t t x dw* , (4.2.13)

*p r xs x*

*k*

2

0 1 2 2

*j j j*

*x C xt C xt P xC xt P <sup>P</sup>*

 

, , 2 2

2

*kj gr* 

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

2 2

<sup>2</sup>

We see that, in equation (4.2.9), the coefficients, *Px i ij* 0,1,2 are functions of *x* only.

 , , ,, *<sup>i</sup>* , ,, *<sup>i</sup> i*

where *A C xt tx* , ,, is the drift coefficient and *B C xt tx <sup>i</sup>* , ,, are diffusion coefficients of the Ito diffusion, equation (4.2.13); here *dw t <sup>i</sup>* are increments of the standard Wiener process. Ito diffusion are an interesting class of stochastic integrals and has many

If we recall the premise on which stochastic calculus is discussed in Chapter 2, *C xt* , is a continuous, non-differentiable function, and equation (4.2.1) should be interpreted as an Ito stochastic differential, which essentially mean equation (4.2.1) has to be understood only in

*j j*

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

*k k*

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

2 2

*j j j kj*

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

*j j*

*p p*

*j j j kj ij kj kj kj k k <sup>h</sup> P x g gx ge g gr x s e* 

2 0 1

advantageous properties of practical importance (Klebaner, 1998).

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

*gr x s e gre*

*kj kj*

*r xs*

2

*j*

*p*

2

*h*

*j ij kj kj kj k*

*P x g gr x s e*

*r xs r xs*

*kj kj kj kj*

4 2 ,

*gr x s e g re C xt*

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

0

1 0 1

the following form:

2 2

Therefore, equation (4.2.5) can be expressed as,

*kj kj kj kj kj*

*j*

*j*

*x t*

*j j*

*p p*

*k k*

2 2

*<sup>h</sup> C xt f C xt f SC xt f*

2

*h x t*

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

*j j x*

*x x*

, , , <sup>2</sup>

*j x*

2 2

*j j*

*p p*

*j kj kj kj j j kj k k*

*g gr x s e g g x ge*

*kj kj kj kj*

*r xs r xs*

2 2

Where

And

The expression of stochastic partial differential equation (equation (4.2.1)) as an Ito diffusion in time would give us the mathematical justification in solving it using Ito calculus. In our development, the eigen functions, *fj x* are continuous, differentiable functions so are the coefficients, *Pxi ij* 0,1,2 . Therefore, the Ito stochastic product rule is the same as the product rule in standard calculus (Klebaner, 1998), and we employ this fact in the previous derivations. Further, we assume that the mean velocity *V xt* , is a continuous, differentiable function which is a reasonable assumption given that the average velocity in aquifer situation is based on the hydraulic conductivity, porosity and the pressure gradient across a large enough domain within a much larger total flow length.

Therefore, we can write the drift term of equation (4.2.1) as follows:

 2 2 2 2 2 2 2 2 , , , , , , , , , , , 2 <sup>2</sup> , , , , , ,2 <sup>2</sup> <sup>2</sup> , , . <sup>2</sup> *x x x x V xt C xt SV xtC xt C xt V xt x x <sup>h</sup> V xt C xt V xt C xt C xt <sup>V</sup> x xx x V xt V xt <sup>h</sup> <sup>h</sup> V xt C xt C xt V xt x x x x h C xt V xt x* (4.2.14)

By substituting equations (4.2.14) and (4.2.9) into equation (4.2.1),

$$\begin{split} d\mathbb{C}(\mathbf{x},t) &= -\left[ \left( \frac{\partial \overline{V}(\mathbf{x},t)}{\partial \mathbf{x}} + \frac{\hbar\_{\mathbf{r}}}{2} \frac{\partial^{2} \overline{V}(\mathbf{x},t)}{\partial \mathbf{x}^{2}} \right) \mathbb{C}(\mathbf{x},t) + \left( \overline{V}(\mathbf{x},t) + \frac{\hbar\_{\mathbf{r}}}{2} 2 \frac{\partial \overline{V}(\mathbf{x},t)}{\partial \mathbf{x}} \right) \frac{\partial \mathbb{C}(\mathbf{x},t)}{\partial \mathbf{x}^{2}} \right] \\ &+ \left( \frac{\hbar\_{\mathbf{r}}}{2} \overline{V}(\mathbf{x},t) \right) \frac{\partial^{2} \mathbb{C}(\mathbf{x},t)}{\partial \mathbf{x}^{2}} \\ &- \sigma \sum\_{j=1}^{m} \sqrt{\lambda\_{j}} \left( P\_{0j} \mathbb{C}(\mathbf{x},t) + P\_{1j} \frac{\partial \mathbb{C}(\mathbf{x},t)}{\partial \mathbf{x}} + P\_{2j} \frac{\partial^{2} \mathbb{C}(\mathbf{x},t)}{\partial \mathbf{x}^{2}} \right) db\_{j}(\mathbf{t}), \end{split}$$

and then,

$$d\mathbf{C}(\mathbf{x},t) = -\mathbf{C}(\mathbf{x},t)d\mathbf{l}\_0 - \frac{\partial \mathbf{C}(\mathbf{x},t)}{\partial \mathbf{x}}d\mathbf{l}\_1 - \frac{\partial^2 \mathbf{C}(\mathbf{x},t)}{\partial \mathbf{x}^2}d\mathbf{l}\_{2'} \tag{4.2.15}$$

where,

$$dI\_0\left(\mathbf{x},t\right) = \left(\frac{\partial \overline{V}\left(\mathbf{x},t\right)}{\partial \mathbf{x}} + \frac{h\_x}{2} \frac{\partial^2 \overline{V}\left(\mathbf{x},t\right)}{\partial \mathbf{x}^2}\right)dt + \sigma \sum\_{j=1}^{m} \sqrt{\lambda\_j} P\_{0\_j} db\_j\left(t\right),\tag{4.2.16}$$

develop an empirical SSTM based on equation (4.2.19). We explore the behaviours of *I xt <sup>i</sup>* , in section 4.5. Next we discuss the derivation of the generalized eigen function, the

A Generalized Mathematical Model in One-Dimension 123

We discuss the approach in this section to obtain the eigen functions for any given kernel in the form given by equation (4.2.3). We calculate the covariance kernel matrix (COV) for any given kernel function and COV can be decomposed in to eigen values and the corresponding eigen vectors using singular value decomposition method or principle component analysis. This can easily be done using mathematical software. Then we use the

2

*y*

<sup>1</sup> , for 0,1, 2, , *<sup>k</sup> x kx k n* , and (4.3.2)

<sup>2</sup> , for 0,1, 2, , *<sup>j</sup> x jx j n* . (4.3.3)

. (4.3.4)

, (4.3.5)

, (4.3.1)

form of which is given by equation (4.2.3).

where 1 2 *y x x* ,

 2 

*b* is the correlation length , and

using the Karhunen-Loève theorem as,

and 2 *x* can be displayed as:

is the variance when 1 2 *x x* .

**4.3 A Computational Approach for Eigen Functions** 

eigen vectors to develop eigen functions using neural networks.

Suppose that we already have an exponential covariance kernel as given by,

1 2 (,)

*<sup>b</sup> q xx e* 

In terms of 1 *x* and 2 *x* , both of them have the domain of [0,*L*]; we equally divide this range into (*n*) equidistant intervals of *x* for both variables. Thus, the particular position for 1 *x*

In equation (4.3.4), COV is the covariance matrix which contains all the variances and covariances, and it is a symmetric matrix with size *n n* where n is the number of intervals. The diagonals of COV represent variances where 1 *x* is equal to 2 *x* and off-

After the covariance matrix is defined, we can transform the COV matrix into a new matrix with new scaled variables according to Karhunen-Loève (KL) theorem. In this new matrix, all the variables are independent of each other having their own variances. i.e, the covariance between any two new variables is zero. The new matrix can be represented by

> 

<sup>2</sup> , for , 1, 1 *kj x*

*e kj n*

By substituting equations (4.3.2) and (4.3.3) into equation (4.3.1), we can obtain,

diagonals represent covariances between any two different discrete 1 *x* and 2 *x* .

1 1 () () *<sup>n</sup> jj j <sup>j</sup> COV x x* 

*<sup>b</sup> COV k j* 

( )

$$dI\_1(\mathbf{x}, t) = \left(\overline{V}(\mathbf{x}, t) + \frac{h\_x}{2} \frac{2\overline{\mathcal{E}}\overline{V}(\mathbf{x}, t)}{\partial \mathbf{x}}\right) dt + \sigma \sum\_{j=1}^{w} \sqrt{\lambda\_j} P\_{1\_j} \, db\_j(t) \text{, and} \tag{4.2.17}$$

$$dI\_2\left(\mathbf{x}, t\right) = \left(\frac{h\_x}{2}\overline{V}\left(\mathbf{x}, t\right)\right)dt + \sigma\sum\_{j=1}^{m}\sqrt{\lambda\_j}P\_{2\_j}\,db\_j\left(\mathbf{t}\right).\tag{4.2.18}$$

In equation (4.2.15), 0 1 *dI dI* , and 2 *dI* are Ito stochastic differentials with respect to *t* as well, and we can rewrite a generalized SSTM given by equation (4.2.1) in terms of another set of stochastic differentials,

$$\mathbf{C}\left(\mathbf{x},t\right) = -\int\_{t} \mathbf{C}\left(\mathbf{x},t\right)dI\_{0} - \int\_{t} \frac{\partial \mathbf{C}\left(\mathbf{x},t\right)}{\partial \mathbf{x}}dI\_{1} - \int\_{t} \frac{\partial^{2} \mathbf{C}\left(\mathbf{x},t\right)}{\partial \mathbf{x}^{2}}dI\_{2} + \mathbf{C}\left(\mathbf{x},\phi\right). \tag{4.2.19}$$

We note that 0 1 *dI dI* , and 2 *dI* are only functions of the mean velocity, *V xt* , , and 0,1,2 *Pxi ij* ; and *C xt* , is separated out in equation (4.2.19), which can be interpreted as *C xt* , and its spatial derivatives are modulating *dI i <sup>i</sup>* 0,1,2 s. *<sup>i</sup> I* are Ito stochastic integrals of the form given by equation (4.2.13), and can be expressed as,

$$I\_i(\mathbf{x}, t) = \int F\_i(\mathbf{x}, t)dt + \sigma \sum\_{j=1}^{m} \int G\_{\neq}(\mathbf{x}, t)db\_j(t), \quad \text{(i.0.1.2) and (j.1.2)}\tag{4.2.20}$$

where,

$$F\_0(\mathbf{x}, t) = \left(\frac{\partial \overline{V}(\mathbf{x}, t)}{\partial \mathbf{x}} + \frac{h\_x}{2} \frac{\partial^2 \overline{V}(\mathbf{x}, t)}{\partial \mathbf{x}^2}\right) \tag{4.2.21}$$

$$F\_1(\mathbf{x}, t) = \left(\overline{V}(\mathbf{x}, t) + \frac{\hbar\_\chi}{2} 2 \frac{\partial \overline{V}(\mathbf{x}, t)}{\partial \mathbf{x}}\right) \tag{4.2.22}$$

$$F\_2(\mathbf{x}, t) = \left(\frac{h\_x}{2}\overline{V}(\mathbf{x}, t)\right)\_{\prime} \tag{4.2.23}$$

$$\mathbf{G}\_{0j} = \sqrt{\mathcal{k}\_j} P\_{0j} \begin{pmatrix} \mathbf{x} \\ \end{pmatrix},\tag{4.2.24}$$

$$\mathcal{G}\_{1/} = \sqrt{\mathcal{k}\_{/}} P\_{1/}(\mathbf{x}),\tag{4.2.25}$$

$$\mathbf{j} \text{ and } \ G\_{2/} = \sqrt{\lambda\_{\uparrow}} P\_{2/} (\mathbf{x}). \tag{4.2.26}$$

As the stochastic integrals, *I xt <sup>i</sup>* , are dependent only on the behaviours of *V xt* , , the eigen values of the velocity covariance kernel and *P x ij* functions, i.e., *I xt <sup>i</sup>* , s are only dependent on the velocity fields within the porous media. A corollary to that is if we know the velocity fields and characterize them as stochastic differentials, we can then

develop an empirical SSTM based on equation (4.2.19). We explore the behaviours of *I xt <sup>i</sup>* , in section 4.5. Next we discuss the derivation of the generalized eigen function, the form of which is given by equation (4.2.3).

### **4.3 A Computational Approach for Eigen Functions**

We discuss the approach in this section to obtain the eigen functions for any given kernel in the form given by equation (4.2.3). We calculate the covariance kernel matrix (COV) for any given kernel function and COV can be decomposed in to eigen values and the corresponding eigen vectors using singular value decomposition method or principle component analysis. This can easily be done using mathematical software. Then we use the eigen vectors to develop eigen functions using neural networks.

Suppose that we already have an exponential covariance kernel as given by,

$$q(\mathbf{x}\_1, \mathbf{x}\_2) = \sigma^2 \quad e^{-\frac{\mathbf{y}}{\theta}} \quad , \tag{4.3.1}$$

where 1 2 *y x x* ,

Computational Modelling of Multi-Scale Non-Fickian Dispersion

*jj j*

(4.2.21)

(4.2.22)

1

 

(4.2.18)

. (4.2.19)

(*i*,0,1,2) and (*j*,1,*m*), (4.2.20)

, and (4.2.17)

*j*

*jj j*

*m*

122 in Porous Media - An Approach Based on Stochastic Calculus

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

1

 

In equation (4.2.15), 0 1 *dI dI* , and 2 *dI* are Ito stochastic differentials with respect to *t* as well, and we can rewrite a generalized SSTM given by equation (4.2.1) in terms of another

<sup>2</sup>

We note that 0 1 *dI dI* , and 2 *dI* are only functions of the mean velocity, *V xt* , , and 0,1,2 *Pxi ij* ; and *C xt* , is separated out in equation (4.2.19), which can be interpreted as *C xt* , and its spatial derivatives are modulating *dI i <sup>i</sup>* 0,1,2 s. *<sup>i</sup> I* are Ito stochastic

, , , , ,

*C xt C xt C x t C x t dI dI dI C x*

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

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

*x x* 

 

As the stochastic integrals, *I xt <sup>i</sup>* , are dependent only on the behaviours of *V xt* , , the eigen values of the velocity covariance kernel and *P x ij* functions, i.e., *I xt <sup>i</sup>* , s are only dependent on the velocity fields within the porous media. A corollary to that is if we know the velocity fields and characterize them as stochastic differentials, we can then

*hx V xt*

*x*

(4.2.23)

(4.2.24)

(4.2.25)

(4.2.26)

0 1 2 2

*x x*

*j*

*m*

*<sup>h</sup> V xt dI x t V x t dt P db t x*

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

*x*

set of stochastic differentials,

where,

*x*

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

*<sup>h</sup> dI x t V x t dt P db t*

, , . <sup>2</sup>

*t t t*

integrals of the form given by equation (4.2.13), and can be expressed as,

 1 , , , , *m*

*F xt V xt*

<sup>1</sup> <sup>1</sup> , *G Px j jj* 

and 2 <sup>2</sup> . *G Px j jj* 

<sup>0</sup> <sup>0</sup> , *G Px j jj*

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

*i i ij j j I x t F x t dt G x t db t* 

*b* is the correlation length , and

 2 is the variance when 1 2 *x x* .

In terms of 1 *x* and 2 *x* , both of them have the domain of [0,*L*]; we equally divide this range into (*n*) equidistant intervals of *x* for both variables. Thus, the particular position for 1 *x* and 2 *x* can be displayed as:

$$\mathbf{x}\_{1k} = k \,\Delta \mathbf{x} \,, \quad \text{for } k = 0, 1, 2, \cdots, n \text{ } \text{ and } \tag{4.3.2}$$

$$\text{fix}\_{2\_{\parallel}} = j \,\text{\AA} \,\text{\mathsf{x}} \,, \quad \text{for } j = 0, 1, 2, \cdots, n \,\text{ }. \tag{4.3.3}$$

By substituting equations (4.3.2) and (4.3.3) into equation (4.3.1), we can obtain,

$$\text{Cov}\left(k,j\right) = \sigma^2 \quad e^{\frac{\left[\left(k-j\right)\cdot\text{Ar}\right]}{\hbar}} \quad \text{for} \quad k,j \in \left[\begin{array}{c} 1, \ n-1 \end{array}\right].\tag{4.3.4}$$

In equation (4.3.4), COV is the covariance matrix which contains all the variances and covariances, and it is a symmetric matrix with size *n n* where n is the number of intervals. The diagonals of COV represent variances where 1 *x* is equal to 2 *x* and offdiagonals represent covariances between any two different discrete 1 *x* and 2 *x* .

After the covariance matrix is defined, we can transform the COV matrix into a new matrix with new scaled variables according to Karhunen-Loève (KL) theorem. In this new matrix, all the variables are independent of each other having their own variances. i.e, the covariance between any two new variables is zero. The new matrix can be represented by using the Karhunen-Loève theorem as,

$$COV1 = \sum\_{j=1}^{n} \lambda\_{/} \phi\_{/}(\mathbf{x}) \,\overline{\phi\_{/}(\mathbf{x})}\,\,\tag{4.3.5}$$

eigenvector. Now let us have the following function the form of which is previously given in

A Generalized Mathematical Model in One-Dimension 125

1 1 ( )2 0 1 0 1

2 2 ( ) (,) *k k p p r xs*

*k k*

where 0 *g* is the bias weight of the network, 1 *g* and *<sup>k</sup> g* is the 1*st* and *th k* weight of the network, and *p* is the total number of neurons in this RBF network (the reason for using *p*+1 for the summation is that *k* starts at 2). Figure 4.1 displays the architecture of a neural

After we decide the input-output mapping and architecture of deterministic neural networks, the next step is to choose the learning algorithm, i.e., the method used to update weights and other parameters of networks. The backpropagation algorithm is used as the learning algorithm in this work. The backpropagation algorithm is used to minimize the network's global error between the actual network outputs and their corresponding desired outputs. The backpropagation leaning method is based on gradient descent that updates weights and other parameters through partial derivatives of the network's global error with respect to the weights and parameters. A stable approach is to change the weights and parameters in the network after the whole sample has been presented and to repeat this process iteratively until the desired minimum error level is reached. This is called batch (or

Figure 4.1. Architecture of the one-dimensional RBF network given by equation (4.3.9).

Their values are based on the summation over all training examples of the partial derivative of the network's global error with respect to the weights and parameters in the whole

 *x g g x g e g g x g Gs r* 

*k k kk*

*x* is used to represent the output of

(4.3.9)

equation (4.2.3) to define the RBF network,

network for the case of one-dimensional input and

the neural network given by equation (4.3.9).

epoch based) learning.

sample.

where *n* is the total number of variables in the new matrix; *<sup>j</sup>* represents the variance of the *th j* rescaled variable and *<sup>j</sup>* is also the *th j* eigen values of the *COV*1 matrix; and ( ) *<sup>j</sup> x* is the eigen vectors of the *COV*1 matrix. The number of eigen vectors depends on the number of discrete intervals. This decomposition of the *COV*1 matrix which is called the singular value decomposition method can easily be done by using mathematical or statistical software.

Once we have the eigen vectors, the next step is to develop suitable neural networks to represent or mimic these eigen vectors. In fact, it is not necessary to simulate all the eigen vectors. The number of neural networks is decided by the number of eigen values which are significant in the KL representation. In some situations, for example, we may have 100 eigen values in the KL representation but only 4 significant eigen values. Then, we just create four networks to simulate these four eigenvectors which correspond to the most significant eigen values. The way we decide on the number of significant eigen values is based on the following equation: eigen the In

$$R[i] = \frac{\lambda\_i}{\sum\_{j=1}^n \lambda\_j} \quad \text{for} \quad i \in \{1, n\} \,\tag{4.3.6}$$

$$\text{Th} = \sum\_{i}^{k} R[i] \text{ .} \tag{4.3.7}$$

where *R i*[ ] represents the contribution of the *th i* eigen value in capturing the total original variance, *k* represents the number of the significant eigen values, and *Th* is the contribution of all significant eigen values in capturing the total original variance. In this chapter, *Th* is chosen to vary between 0.95 and 1. This means that if the total number of eigen values is 100 from the KL expansion and the contribution of the first 4 eigen values takes up more than 95% of the original variance, there are only four individual neural networks that need to be developed.

The main factors that need to be decided in the development of neural networks are the number of neurons needed, the structure of neural networks and the learning algorithm. The number of neurons in neural networks is case-dependent. It is difficult to define the number of neurons before the learning stage. In general, the number of neurons is adjusted during training until the network output converges on the actual output based on least square error minimization. A neural network with an optimum number of neurons will reach the desired minimum error level more quickly than other networks with more complex structures. The proposed neural network is a Radial Basis Function (RBF) Network. The approximation function is the Gaussian Function given below:

$$G(s,r) = e^{-r \cdot (x-s)^2} \, , \tag{4.3.8}$$

where *r* and *s* are constants. In this symmetric function, *s* defines the centre of symmetry and *r* defines the sharpness of Gaussian function.

Based on numerical values of the significant eigen vectors, several RBF networks with one input ( *x* ) and one output (eigen vector) are developed to approximate each significant

Computational Modelling of Multi-Scale Non-Fickian Dispersion

*j* eigen values of the *COV*1 matrix; and ( ) *<sup>j</sup>*

*<sup>j</sup>* represents the variance of the

, (4.3.6)

*Th R i* , (4.3.7)

<sup>2</sup> ( ) (,) *r xs Gsr e* , (4.3.8)

*i* eigen value in capturing the total original

*x* is

124 in Porous Media - An Approach Based on Stochastic Calculus

the eigen vectors of the *COV*1 matrix. The number of eigen vectors depends on the number of discrete intervals. This decomposition of the *COV*1 matrix which is called the singular value decomposition method can easily be done by using mathematical or statistical

Once we have the eigen vectors, the next step is to develop suitable neural networks to represent or mimic these eigen vectors. In fact, it is not necessary to simulate all the eigen vectors. The number of neural networks is decided by the number of eigen values which are significant in the KL representation. In some situations, for example, we may have 100 eigen values in the KL representation but only 4 significant eigen values. Then, we just create four networks to simulate these four eigenvectors which correspond to the most significant eigen values. The way we decide on the number of significant eigen values is based on the

[ ] for [1, ] *<sup>i</sup>*

[ ] *<sup>k</sup> i*

*R i i n*

variance, *k* represents the number of the significant eigen values, and *Th* is the contribution of all significant eigen values in capturing the total original variance. In this chapter, *Th* is chosen to vary between 0.95 and 1. This means that if the total number of eigen values is 100 from the KL expansion and the contribution of the first 4 eigen values takes up more than 95% of the original variance, there are only four individual neural

The main factors that need to be decided in the development of neural networks are the number of neurons needed, the structure of neural networks and the learning algorithm. The number of neurons in neural networks is case-dependent. It is difficult to define the number of neurons before the learning stage. In general, the number of neurons is adjusted during training until the network output converges on the actual output based on least square error minimization. A neural network with an optimum number of neurons will reach the desired minimum error level more quickly than other networks with more complex structures. The proposed neural network is a Radial Basis Function (RBF) Network.

where *r* and *s* are constants. In this symmetric function, *s* defines the centre of

Based on numerical values of the significant eigen vectors, several RBF networks with one input ( *x* ) and one output (eigen vector) are developed to approximate each significant

The approximation function is the Gaussian Function given below:

symmetry and *r* defines the sharpness of Gaussian function.

1

*n <sup>j</sup> <sup>j</sup>*

where *n* is the total number of variables in the new matrix;

*j* is also the *th*

where *R i*[ ] represents the contribution of the *th*

networks that need to be developed.

*th*

software.

following equation:

*j* rescaled variable and

eigenvector. Now let us have the following function the form of which is previously given in equation (4.2.3) to define the RBF network,

$$\boldsymbol{\phi}(\mathbf{x}) = \mathbf{g}\_0 + \mathbf{g}\_1 \,\,\mathbf{x} + \sum\_{k=2}^{p+1} \mathbf{g}\_k \,\,\mathbf{e}^{-\mathbf{e}\_1 \,(\mathbf{x} - \mathbf{e}\_1)2} = \mathbf{g}\_0 + \mathbf{g}\_1 \,\,\mathbf{x} + \sum\_{k=2}^{p+1} \mathbf{g}\_k \,\,\mathbf{G}(\mathbf{s}\_k, r\_k) \tag{4.3.9}$$

where 0 *g* is the bias weight of the network, 1 *g* and *<sup>k</sup> g* is the 1*st* and *th k* weight of the network, and *p* is the total number of neurons in this RBF network (the reason for using *p*+1 for the summation is that *k* starts at 2). Figure 4.1 displays the architecture of a neural network for the case of one-dimensional input and *x* is used to represent the output of the neural network given by equation (4.3.9).

After we decide the input-output mapping and architecture of deterministic neural networks, the next step is to choose the learning algorithm, i.e., the method used to update weights and other parameters of networks. The backpropagation algorithm is used as the learning algorithm in this work. The backpropagation algorithm is used to minimize the network's global error between the actual network outputs and their corresponding desired outputs. The backpropagation leaning method is based on gradient descent that updates weights and other parameters through partial derivatives of the network's global error with respect to the weights and parameters. A stable approach is to change the weights and parameters in the network after the whole sample has been presented and to repeat this process iteratively until the desired minimum error level is reached. This is called batch (or epoch based) learning.

Figure 4.1. Architecture of the one-dimensional RBF network given by equation (4.3.9).

Their values are based on the summation over all training examples of the partial derivative of the network's global error with respect to the weights and parameters in the whole sample.

Now let us assume that the actual network output is *T* ( ( ) *x* ), the desired network output is *Z* (( ) *x* ). The network's global error between the network output and actual output is

$$E = \frac{1}{2} \sum\_{i}^{N} \left( T\_i - Z\_i \right)^2,\tag{4.3.10}$$

Therefore, we need to focus on the significant eigen values as well as their corresponding eigen functions. There are 6 significant eigen values whose contribution takes up 99.9035% of the original variance and table 4.1 shows the value of each significant eigen value and the

A Generalized Mathematical Model in One-Dimension 127

Thus, the six eigenvectors corresponding to the significant eigenvalues are simulated by the individual RBF network. Although we use a different RBF network to approximate each of these six eigenvectors, the structure of RBF network is the same but the weights and

Figure 4.2. The covariance matrix calculated by the given covariance kernel (equation 4.3.12)

The number of eigen values Values Contribution as a propotion

Table 4.1. Significant eigen values obtained from the KL expansion of the covariance kernel given by equation (4.3.12) (six significant eigen values take up 99.9035% of total variance)

<sup>1</sup> 48.128 0.477

<sup>2</sup> 30.636 0.303

<sup>3</sup> 14.692 0.145

<sup>4</sup> 5.446 0.054

<sup>5</sup> 1.61 0.016

<sup>6</sup> 0.392 0.004

proportion of variance captured by the corresponding eigen value.

parameters inside the individual RBF network are different.

when 2 

1and 0.1 *b* .

where *Ti* and *Zi* are the actual output and the network output for the *th i* training pattern, and *N* is the total number of training patterns. The multiplication by 1 / 2 is a mathematical convenience (Samarasinghe, 2006).

The method of modifying a weight or a parameter is the same for all weights and parameters so we show the change to an arbitrary weight as an example. The change to a single weight of a connection between neuron *j* and neuron *i* in the RBF network based on batch learning can be defined as,

$$\Delta \mathbf{g}\_{\boldsymbol{\beta}} = \eta \sum\_{p=1}^{k} \left( \frac{dE}{d\mathbf{g}\_{\boldsymbol{\beta}}} \right)\_p \prime \tag{4.3.11}$$

where is called the learning rate with a constant value between 0 and 1. It controls the step size and the speed of weight adjustments. *k* is the total number of input vectors. The process that propagates the error information backwards into the network and updates weights and the parameters of network is repeated until the network minimizes the global error between the actual network outputs and their corresponding desired outputs. In the learning process, the weights and the parameters of the network converge on the optimal values.

To illustrate the computational approach, we give some examples here. The first covariance kernel is chosen to be

$$\mathbf{g(x\_1, x\_2)} = \sigma^2 \ e^{-\frac{y}{b}} \, \, \, \, \, \, \, \tag{4.3.12}$$

where <sup>2</sup> 1 2 *yxx* ,

*b* is constant , and

 <sup>2</sup> is variance when 1 2 *x x* .

This covariance kernel needs a relatively lower number of significant eigen values to capture 95% or more of the total variance; therefore we choose to work with this kernel and later we use two other forms of kernels: one discussed previously in Chapter 3 and other one is empirically based.

Figure 4.2 displays the covariance matrix based on the covariance kernel (equation (4.3.12)) when 2 1 and 0.1 *b* .

Table 4.1 reports all eigen values in the KL representation of the covariance matrix. The most of the eigen values is equal to zero and these eigen values can not affect the covariance matrix and just a few are significant and capture the total variance in the original data.

*<sup>N</sup>* , (4.3.10)

, (4.3.11)

, (4.3.12)

( ) *x* ), the desired network output

*i* training

126 in Porous Media - An Approach Based on Stochastic Calculus

2

where *Ti* and *Zi* are the actual output and the network output for the *th*

( ) *x* ). The network's global error between the network output and actual output is

*N*

*i E TZ*

pattern, and *N* is the total number of training patterns. The multiplication by 1 / 2 is a

The method of modifying a weight or a parameter is the same for all weights and parameters so we show the change to an arbitrary weight as an example. The change to a single weight of a connection between neuron *j* and neuron *i* in the RBF network based on

> 1 *k*

step size and the speed of weight adjustments. *k* is the total number of input vectors. The process that propagates the error information backwards into the network and updates weights and the parameters of network is repeated until the network minimizes the global error between the actual network outputs and their corresponding desired outputs. In the learning process, the weights and the parameters of the network converge on the optimal

To illustrate the computational approach, we give some examples here. The first covariance

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

This covariance kernel needs a relatively lower number of significant eigen values to capture 95% or more of the total variance; therefore we choose to work with this kernel and later we use two other forms of kernels: one discussed previously in Chapter 3 and other one is

Figure 4.2 displays the covariance matrix based on the covariance kernel (equation (4.3.12))

Table 4.1 reports all eigen values in the KL representation of the covariance matrix. The most of the eigen values is equal to zero and these eigen values can not affect the covariance matrix and just a few are significant and capture the total variance in the original data.

1 2 (,)

2

*y*

*<sup>g</sup> dg* 

*<sup>p</sup> ji <sup>p</sup> dE*

is called the learning rate with a constant value between 0 and 1. It controls the

*ji*

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

*i i*

Now let us assume that the actual network output is *T* (

mathematical convenience (Samarasinghe, 2006).

batch learning can be defined as,

is *Z* (

where

values.

 <sup>2</sup> 

when 2 

kernel is chosen to be

where <sup>2</sup>

empirically based.

1 2 *yxx* ,

1 and 0.1 *b* .

is variance when 1 2 *x x* .

*b* is constant , and

Therefore, we need to focus on the significant eigen values as well as their corresponding eigen functions. There are 6 significant eigen values whose contribution takes up 99.9035% of the original variance and table 4.1 shows the value of each significant eigen value and the proportion of variance captured by the corresponding eigen value.

Thus, the six eigenvectors corresponding to the significant eigenvalues are simulated by the individual RBF network. Although we use a different RBF network to approximate each of these six eigenvectors, the structure of RBF network is the same but the weights and parameters inside the individual RBF network are different.

Figure 4.2. The covariance matrix calculated by the given covariance kernel (equation 4.3.12) when 2 1and 0.1 *b* . 4.3.12)


Table 4.1. Significant eigen values obtained from the KL expansion of the covariance kernel given by equation (4.3.12) (six significant eigen values take up 99.9035% of total variance)

Eigen values Values The analytical forms of eigen functions

0.745057

9.47326

88.8399

significant eigen vector and their corresponding eigen values.

Figure 4.4. The covariance matrix of equation (4.3.12) when 2

<sup>1</sup> 48.128 <sup>2</sup> <sup>12</sup> 2.34049 ( 0.5) 0.0456723 1.61025 10 0.170974 *<sup>x</sup> x e*

*e*

2

*x e*

6.6592 ( 0.683601) 9.56556 ( 0.0999343)

14.6016 ( 0.306358) 0.385895 ( 0.699959)

2 2

*x e*

*x x*

*x x x e*

8.41938 ( 0.5) 5.67596 ( 0.5)

*x x*

5.57986 ( 0.511944) 5.3708 ( 0.46363)

1 and *b* = 0.2.

*x x*

*x e*

2.56631 ( 0.335246) 0.689211 0.560102 0.360116

0.00132933 0.0185462 3.7566 3.7712 0.183884

0.282871 12.8992

0.0403405 3.13566 10 9.47326

*e e*

8.95446 ( 0.272435)

0.318516 0.548883 127.699 429.099 237.989

8.18627 ( 0.433595)

*x*

Table 4.2. The final formula from the developed RBF networks to approximate each

*x*

19.0735 31.1501

 

*e*

*e*

*e e*

*x e e e*

2

2

*e e*

*x*

 

A Generalized Mathematical Model in One-Dimension 129

 

10.4236 5.90955 0.3771

 

<sup>2</sup> 30.636

<sup>3</sup> 14.692

<sup>4</sup> 5.446

<sup>5</sup> 1.61

<sup>6</sup> 0.392

*i*

5.37562 ( 0.290915)

6.62157 ( 0.703168)

*x*

2

7.08944 ( 0.658744)

*x*

*x*

2

*x*

2 2

*x*

2 2 13.4883 ( 0.68965)


2 2

( ) *x* for *x* 0,1

2

2

2

Then each individual RBF network was trained to learn their related eigen vector, while the weights and all the parameters of Gaussian functions were updated until each network reaches global minimum error given by

$$\mathcal{E} = \frac{1}{N-1} \sum\_{\boldsymbol{\alpha}=1}^{N} \left[ \boldsymbol{\phi}\_{\boldsymbol{\alpha}}^{\boldsymbol{\prime}} (\boldsymbol{\chi}) - \boldsymbol{S}\_{\boldsymbol{\alpha}}^{\boldsymbol{\prime}} (\boldsymbol{\chi}) \right]^2. \tag{4.3.13}$$

where ( ) *<sup>i</sup> <sup>n</sup> S x* is an approximation to ( ) *<sup>i</sup> n x* , *i* th eigen function; and *N* is the total number of eigen values.

Figure 4.3 displays all the six eigen functions and Table 4.2 gives their functional forms. Figure 4.3 shows the eigen functions given by the KL theory (dots), obtained by solving the corresponding integral equation, overlaid with the outputs from the neural networks (lines), and the approximation functions are the same as the theoretically derived functions.

Figure 4.3. The approximated six eigen functions from BRF networks for equation (4.3.12) when 2 1 and *b*=0.1.

For the second example, we use the same covariance kernel with 2 *b* 0.2 and 1 , and decomposed the matrix for the domain [0, 1]. Figure 4.4 shows the covariance matrix calculated by the given covariance kernel. Based on the standard of choosing the number of significant eigen values, it can be seen that five significant eigen values together capture 99.9582% of the original variance. Thus, there are five RBF networks to be developed for the corresponding eigenvectors. Table 4.3 gives the value of each significant eigen value and the proportion of variance captured by the corresponding eigen value; Table 4.4 provides the eigen functions obtained from the neural networks; the graphical forms of the eigen functions are given by Figure 4.5. As before, the analytically derived eigen function as an very well approximated by the approximations from the neural networks. (The theoretical eigen functions are not displayed in Figure 4.5).

. (4.3.13)

*x* , *i* th eigen function; and *N* is the total number of

, and

2

128 in Porous Media - An Approach Based on Stochastic Calculus

Then each individual RBF network was trained to learn their related eigen vector, while the weights and all the parameters of Gaussian functions were updated until each network

> *i i n n*

*x Sx*

1 <sup>1</sup> () () <sup>1</sup>

Figure 4.3 displays all the six eigen functions and Table 4.2 gives their functional forms. Figure 4.3 shows the eigen functions given by the KL theory (dots), obtained by solving the corresponding integral equation, overlaid with the outputs from the neural networks (lines),

Figure 4.3. The approximated six eigen functions from BRF networks for equation (4.3.12)

decomposed the matrix for the domain [0, 1]. Figure 4.4 shows the covariance matrix calculated by the given covariance kernel. Based on the standard of choosing the number of significant eigen values, it can be seen that five significant eigen values together capture 99.9582% of the original variance. Thus, there are five RBF networks to be developed for the corresponding eigenvectors. Table 4.3 gives the value of each significant eigen value and the proportion of variance captured by the corresponding eigen value; Table 4.4 provides the eigen functions obtained from the neural networks; the graphical forms of the eigen functions are given by Figure 4.5. As before, the analytically derived eigen function as an very well approximated by the approximations from the neural networks. (The theoretical

For the second example, we use the same covariance kernel with 2 *b* 0.2 and 1

*n*

and the approximation functions are the same as the theoretically derived functions.

*N*

*i n* 

*N*

reaches global minimum error given by

*<sup>n</sup> S x* is an approximation to ( )

where ( ) *i*

eigen values.

when 2 

1 and *b*=0.1.

eigen functions are not displayed in Figure 4.5).


Table 4.2. The final formula from the developed RBF networks to approximate each significant eigen vector and their corresponding eigen values.

Figure 4.4. The covariance matrix of equation (4.3.12) when 2 1 and *b* = 0.2.

the development of the SSTM in Chapter 3. In Chapter 3 the covariance kernel given by equation (4.3.14) - we reproduce the equation here- constitutes an integral equation which

A Generalized Mathematical Model in One-Dimension 131

2

Figure 4.6. The covariance matrix calculated by the given equation (4.3.14) under the

Eigen values Values Contribution as a proportion

<sup>1</sup> 18.745 0.186

<sup>2</sup> 15.681 0.155

<sup>3</sup> 12.218 0.121

<sup>4</sup> 9.241 0.091

<sup>5</sup> 6.976 0.069

<sup>6</sup> 5.332 0.053

<sup>7</sup> 4.149 0.041

<sup>8</sup> 3.293 0.033

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

1 2 (,)

1 2

, (4.3.14)

*x x*

 

we solve analytically to obtain eigen values and eigen functions:

 and *b* have the same meanings as before. This covariance kernel is depicted graphically in Figure 4.6.

when <sup>2</sup> 

condition 2

1 .


Table 4.3. Relative amounts of variance in the data captured by each significant eigenvalue obtained from the KL expansion of equation (4.3.12) when <sup>2</sup> 1 and *b*=0.2.


Table 4.4. The eigen functions from the developed RBF networks to approximate each significant eigen vector for the kernel given by equation (4.3.12) when <sup>2</sup> 1 and *b*=0.2.

Figure 4.5. The approximated six eigen functions from BRF networks for equation (4.3.12) when 2 1 and 0.2 *b* .

From the previous two examples, we have seen that the covariance kernel given by equation (4.3.12) provides a relative small number of eigen functions and therefore one may say that the kernel given by equation (4.3.12) has fast convergence. This is quite a desirable property to have, especially in terms of computational efficiency of the algorithms. In the next example, we find the eigen values and the eigen functions of the covariance kernel we use in

*x e*

*x e*

3.70961 ( 0.536422) 3.5887 ( 0.459753)

14.0892 ( 0.212252) 14.9973 ( 0.226353)

*x x x e*

*x x x e*

2

2

*e e*

*e e*

1 and *b*=0.2.


2 2

2 2

*i*

4.39837 ( 0.820696)

*x*

*x*

2

2

1 and *b*=0.2.

3.81696 ( -0.559606)

*x*

11.2297 ( 0.826432)

*x*

( ) *x* for *x* 0,1

2

2

130 in Porous Media - An Approach Based on Stochastic Calculus

Table 4.3. Relative amounts of variance in the data captured by each significant eigenvalue

4.39837 ( 0.179304) 0.000620869 0.00124174 0.147148

13.4824 ( 0.374103) 0.201764 3.67217 10 0.197208

*x*

0.0357425 0.0604364 176.588 231.215 56.2021

0.195412 0.641926 0.834518 0.479846 0.649793

Table 4.4. The eigen functions from the developed RBF networks to approximate each

Figure 4.5. The approximated six eigen functions from BRF networks for equation (4.3.12)

From the previous two examples, we have seen that the covariance kernel given by equation (4.3.12) provides a relative small number of eigen functions and therefore one may say that the kernel given by equation (4.3.12) has fast convergence. This is quite a desirable property to have, especially in terms of computational efficiency of the algorithms. In the next example, we find the eigen values and the eigen functions of the covariance kernel we use in

*x*

 

Eigen values Values Contribution as a proportion

<sup>1</sup> 61.242 0.606

<sup>2</sup> 28.820 0.285

<sup>3</sup> 8.747 0.087

<sup>4</sup> 1.853 0.018

obtained from the KL expansion of equation (4.3.12) when <sup>2</sup>

<sup>5</sup> 0.29597 0.00293

Eigen values Values The analytical forms of eigen functions

0.147148

0.197207

for

<sup>1</sup> 61.242 <sup>2</sup> <sup>17</sup> 2.18805 ( 0.5) 0.0179462 7.11262 10 0.137587 *<sup>x</sup> x e*

*e*

*e*

 

significant eigen vector for the kernel given by equation (4.3.12) when <sup>2</sup>

when 2 

<sup>2</sup> 28.820

<sup>3</sup> 8.747

<sup>4</sup> 1.853

<sup>5</sup> 0.29597

1 and 0.2 *b* .

the development of the SSTM in Chapter 3. In Chapter 3 the covariance kernel given by equation (4.3.14) - we reproduce the equation here- constitutes an integral equation which we solve analytically to obtain eigen values and eigen functions:

$$\log(\mathbf{x}\_1, \mathbf{x}\_2) = \sigma^2 \ e^{-\frac{\left|\mathbf{x}\_1 - \mathbf{x}\_2\right|}{b}},\tag{4.3.14}$$

when <sup>2</sup> and *b* have the same meanings as before.

This covariance kernel is depicted graphically in Figure 4.6.

Figure 4.6. The covariance matrix calculated by the given equation (4.3.14) under the condition 2 1 .


In Table 4.5, we give the 32 eigen values which capture up to 94% of the only original variance; Table 4.6 gives only the first six eigen functions obtained from the networks for

A Generalized Mathematical Model in One-Dimension 133

2

*e e*

*e e*

*x e*

7.13097 ( 0.592419) 10.6256 ( 0.0883508)

*x x x e*

5.66061 ( 0.603547 ) 6.81302 ( 0.448057 )

2 2

*x e*

*x x x e*

7.12542 ( 0.5) 9.19718 ( 0.388233)

4.883 ( 0.215442 ) 1.85592 ( 0.0783536)

13.3089 ( 0.577832 ) 10.1574 ( 0.488994)

8.63616 ( 0.807765) 5.03775 ( 0.766933) 13.5694 ( 0.3499) 22.3578 ( 0.307475)

*x x x x*

*x x*

*x e*

 

*x e*

*x x*

*x x*

2

2

*e e e e*

Table 4.6. The final formula from the developed RBF networks to approximate the first 8

*e e*

*x e e e*

5.79605 ( 0.169048) 0.128953 0.257905 0.23274

*x*

0.297969 0.071349 1.99033 1.74606 0.339719

0.520716 0.848337 182.653 214.465 38.6999

1.10289 1.12328 10 55.2051

*e e*

99.6664 55.2051

1.30206 ( 0.0105042)

2.16846 0.148883 1451.69 2198.81 1013.85

13.1355 ( 0.350845)

*x*

31.8294 108.195 6255.59 7941.31 1838.53 282.847 52.6999

*x*

190.785 2539.25

822.58 496.01 6.74448

 

 

*e*

*e*

 

 

  *i*

2

2

2

2

2

2

5.79605 ( 0.830952 )

5.41338 ( 0.640904)

5.89655 ( 0.636135)

*x*

2

14.2086 ( 0.592557 )

*x*

9.37148 ( 0.813402 )

*x*

*x*

*x*

*x*

2 2

2 2

2 2

2 2

2 2 2 2

17.354 ( 0.590219)

*x*

10 9.19718 ( 0.611767 )

( ) *x* for *x* 0,1

brevity. Figure 4.7 shows the first eight eigen functions.

values Values The analytical forms of eigen functions

0.23274

3181.28

358.254

significant eigen vectors and their corresponding eigen values.

<sup>1</sup> 18.745 <sup>2</sup> <sup>16</sup> 1.19444 ( 0.5) 0.238941 7.92548 10 0.368191 *<sup>x</sup> x e*

*e*

Eigen

<sup>2</sup> 15.681

<sup>3</sup> 12.218

<sup>4</sup> 9.241

<sup>5</sup> 6.976

<sup>6</sup> 5.332

<sup>7</sup> 4.149

<sup>8</sup> 3.293


Table 4.5. The eigen values for the kernel given by equation (4.3.14) (32 significant eigen values which take up 94.29% of original variance)

132 in Porous Media - An Approach Based on Stochastic Calculus

<sup>9</sup> 2.661 0.026

<sup>10</sup> 2.188 0.022

<sup>11</sup> 1.826 0.018

<sup>12</sup> 1.545 0.015

<sup>13</sup> 1.323 0.013

<sup>14</sup> 1.145 0.011

<sup>15</sup> 0.9999 0.0099

<sup>16</sup> 0.881 0.0087

<sup>17</sup> 0.782 0.0077

<sup>18</sup> 0.699 0.0069

<sup>19</sup> 0.628 0.0062

<sup>20</sup> 0.568 0.0056

<sup>21</sup> 0.516 0.0051

<sup>22</sup> 0.471 0.0047

<sup>23</sup> 0.432 0.0043

<sup>24</sup> 0.398 0.0039

<sup>25</sup> 0.367 0.0036

<sup>26</sup> 0.341 0.0034

<sup>27</sup> 0.317 0.0031

<sup>28</sup> 0.295 0.0029

<sup>29</sup> 0.276 0.0027

<sup>30</sup> 0.259 0.0026

<sup>31</sup> 0.244 0.0024

<sup>32</sup> 0.229 0.0023

values which take up 94.29% of original variance)

Table 4.5. The eigen values for the kernel given by equation (4.3.14) (32 significant eigen


In Table 4.5, we give the 32 eigen values which capture up to 94% of the only original variance; Table 4.6 gives only the first six eigen functions obtained from the networks for brevity. Figure 4.7 shows the first eight eigen functions.

Table 4.6. The final formula from the developed RBF networks to approximate the first 8 significant eigen vectors and their corresponding eigen values.

1 2 1 2

*x x for x x*

Equation (4.3.15) is depicted in Figure 4.8, and Figure 4.9 shows the corresponding

Table 4.7 gives the most significant eigen values (the first nine values); Table 4.8 shows the functional forms of eigen functions and Figure 4.10 shows the graphical forms of eigen

Figure 4.9. The covariance matrix given by the empirical distribution (equation (4.3.15)).

 

A Generalized Mathematical Model in One-Dimension 135

1 0 0.5

*for x x*

0 0.8 1

1 2

, (4.3.15)

1 2 1 2 1 2

*Cov x x x x for x x*

0

3

covariance matrix.

Figure 4.8. An empirical distribution.

functions.

<sup>5</sup> ( , ) ( 0.8) 0.5 0.8

Figure 4.7. The approximated first eight eigen functions from RBF network for the equation (4.3.14) when 2 1 and *b* 0.1 .

We can also use an empirically derived covariance kernel. As an example, let us consider equation (4.3.15).

134 in Porous Media - An Approach Based on Stochastic Calculus

Figure 4.7. The approximated first eight eigen functions from RBF network for the equation

We can also use an empirically derived covariance kernel. As an example, let us consider

(4.3.14) when 2

equation (4.3.15).

1 and *b* 0.1 .

$$\text{Cov}(\mathbf{x}\_1, \mathbf{x}\_2) = \begin{cases} -\left| \mathbf{x}\_1 - \mathbf{x}\_2 \right| + 1 & \text{for} \quad 0 \le \left| \mathbf{x}\_1 - \mathbf{x}\_2 \right| \le 0.5\\ -\frac{5}{3} (\left| \mathbf{x}\_1 - \mathbf{x}\_2 \right| - 0.8) & \text{for} \quad 0.5 < \left| \mathbf{x}\_1 - \mathbf{x}\_2 \right| \le 0.8\\ 0 & \text{for} \quad 0.8 < \left| \mathbf{x}\_1 - \mathbf{x}\_2 \right| \le 1 \end{cases} \tag{4.3.15}$$

Equation (4.3.15) is depicted in Figure 4.8, and Figure 4.9 shows the corresponding covariance matrix.

Table 4.7 gives the most significant eigen values (the first nine values); Table 4.8 shows the functional forms of eigen functions and Figure 4.10 shows the graphical forms of eigen functions.

Figure 4.8. An empirical distribution.

Figure 4.9. The covariance matrix given by the empirical distribution (equation (4.3.15)).

<sup>7</sup> 0.395

<sup>8</sup> 0.390

<sup>9</sup> 0.201

given in equation (4.3.15).

2

2 2

8 4.2926 ( 0.649682 )

*x*

2

2 2

8.1404 ( 0.811728)

*x*

2

27.8868 ( 0.860872 )

*x*

2 2 2

*x e*

2

*e e*

8 3.99793 ( 0.630907 )

*x*

7 2.43237 ( 0.572872 )

7 5.91597 ( 0.47122 ) 6 8.42978 ( 0.395813) 6 8.86148 ( 0.38

*x x x*

2 8545)

Table 4.8. The final formulas from the developed RBF networks to approximate each

Figure 4.10. The approximated nine eigenfunctions from the BRF networks for the kernel

We have seen that some covariance kernels provide relatively small number of significant eigen values for the given domain [0,1] and whereas for others, obtaining the significant

151130 1.6007 10

 

*x x*

.396187 )

*x e e e*

8.57437 10 702189

 

6 5.69446 ( 0

<sup>2</sup> 0.100811)

4.26386 0.388207

386986 334630 1.13814 10

*e*

*e e*

<sup>2</sup>

<sup>6</sup> 1.15913 ( 0.0920272 ) 1.62219 10 *<sup>x</sup> e* 

6589.2 102.051 277640

*e e e*

<sup>2</sup>

3.56779 ( 0.288557 ) 25798.1 *<sup>x</sup> e* 

4.50902 ( 0.43218) 183.561 (

*x x*

*x e*

2

2

2

2 2

3.04348 ( 0.528018) 7 5.70678 ( 0.480254)

*x x*

*e e*

7 2.777 ( 0.591124 ) 7.32322 ( 0.57722 )

*x x*

7 25.9199 ( 0.765634) 7 26.0007 ( 0.764335) 7 26.0126 ( 0.764145)

 

*x x x*

A Generalized Mathematical Model in One-Dimension 137

2.02737 7.9356 5263.75

*e e e*

1.03168 10 7.83776 10 6.80653 10

1.52714 10

5.31436 10 6.4929 10

1.67931 10 2.8962 10 1.79879 10

 

significant eigenvector and their corresponding eigen values.

 


Table 4.7. Relative amount of variance in the data captured by each significant eigen value for the kernel in equation (4.3.15). (9 significant eigen values capture 98% of original variance)


*i*

2

2

8.33766 ( 0.951692 )

*x*

7 0.810094 ( 0.374112 )

*x*

2

2

( ) *x* for *x* 0,1

136 in Porous Media - An Approach Based on Stochastic Calculus

Table 4.7. Relative amount of variance in the data captured by each significant eigen value for the kernel in equation (4.3.15). (9 significant eigen values capture 98% of original variance)

0.0495323 3.28499 10 0.165441 *<sup>x</sup>*

2

2 2

3.15124 ( 0.0721231) 2.33563 ( 0.0182454)

4.31785 ( 0.0230752) 174.457 ( 0.074625 )

*x*

*x x*

*x x x e e e*

*x e e e*

> 7 6 2.78874 ( 1.40277 ) 1.17605 ( 0.688669) 0.360088 ( -0.448068) 0.206739 ( 0.405528) 2.08493 ( 0.139661)

*x x x x x*

 

*x e*

8.33761 ( 0.0483079) 0.0156028 0.0312062 0.127272

672.733 332.866 1.19376 10

 

7 0.807196 ( 0.375468) 6 0.798664 ( 0.379503 )

*x x*

*x*

 

> *e e*

93.5157 68.0722 813.179 1577.84 860.897

4.82583 4.6287 1.32608

 

9.87408 10 3.76629 10

*e e e e e*

 <sup>2</sup> <sup>7</sup> 2.10254 ( 0.336545) .06435 10 *<sup>x</sup> e* 

4.73704 0.196755

 

*x e*

2 2 3.71071 ( 0.0967842 )

2 2 14.504 ( -0.624954 )

*x*

*x*

*x e*

Eigen values Values Contribution as a proportion

<sup>1</sup> 66.022 0.6537

<sup>2</sup> 24.216 0.2398

<sup>3</sup> 3.163 0.0313

<sup>4</sup> 1.989 0.0197

<sup>5</sup> 0.447 0.0189

<sup>6</sup> 0.445 0.00443

<sup>7</sup> 0.395 0.0044

<sup>8</sup> 0.390 0.00391

<sup>9</sup> 0.201 0.00386

Eigen values Values The analytical forms of eigen functions

0.127273

1.60589 10 4.12218 10

1.39729 10 4.53545 10 3.41918 10 3.84429 10 1.13264 10

 

1

 

<sup>1</sup> 66.022 <sup>2</sup> <sup>16</sup> 1.38223 ( 0.5)

*e*

<sup>2</sup> 24.216

<sup>3</sup> 3.163

<sup>4</sup> 1.989

<sup>5</sup> 0.447

<sup>6</sup> 0.445


Table 4.8. The final formulas from the developed RBF networks to approximate each significant eigenvector and their corresponding eigen values.

Figure 4.10. The approximated nine eigenfunctions from the BRF networks for the kernel given in equation (4.3.15).

We have seen that some covariance kernels provide relatively small number of significant eigen values for the given domain [0,1] and whereas for others, obtaining the significant kernels

We investigate the effects of the kernels on the dispersivity values; we compute them using the stochastic inverse method (SIM) discussed in Chapter 3. Table 4.9 shows the results. For all practical purposes they are essentially the same. The mechanic of dispersion is more

A Generalized Mathematical Model in One-Dimension 139

and *b* . The mechanics of dispersion in general can also be assumed to be influenced by the mathematical form of the kernel. The both of these kernels are exponential decaying functions. Because of the case of computations, we continue to use kernel 2 in Table 4.9 in

Figure 4.11. The generalized SSTM 95% confidence intervals for the concentration

Figure 4.12. The generalized SSTM 95% confidence intervals for the concentration

= 0.1 and *b* =0.1.

 when <sup>2</sup> 

= 0.1 and *b* =0.1.

 when <sup>2</sup> 

1 2

 <sup>2</sup> 1 2

*x x b*

 

2

*e* *x x b*

 

2

*e* and *b* are allowed to vary, on both <sup>2</sup>

them

for a given *b* or if both <sup>2</sup>

influenced by <sup>2</sup>

realization for the kernel,

realization for the kernel,

the most of the work discussed in this book.

eigen functions could be tedious. The main point in this exercise is to show that any given covariance kernel can be used to obtain the eigen functions of the form given by equation (4.2.3). A corollary to this statement is that if we express the eigen functions in the form given by equation (4.2.3), then we can assume that these is an underlying covariance kernel responsible for these eigen functions. In deriving the SSTM in the form of equation (4.2.15), we assume that the form of equation (4.2.3) is given. But we see now that any covariancekernel driven SSTM can be represented by equation (4.2.15).

### **4.4 Effects of Different Kernels and** *<sup>x</sup> h*

We have seen in sections 4.2 and 4.3 that the SSTM developed in Chapter 2 and 3 can be recasted so that we could employ any given velocity kernel. In fact, we can even use an empirical set of data for the velocity covariance. We can anticipate that the generalized SSTM would behave quite similar to the one developed in Chapter 2 given that same covariance kernel is used. We compute the 95% confidence intervals for the concentration breakthrough curves (concentration realizations) at *x* = 0.5 m when the flow length is 1 m to compare the differences that occur in using different kernels in the generalized SSTM. The mean velocity is kept constant at 0.5 m/day, and the covariance kernel given by equation (4.3.14) is used to obtain Figure 4.11, and Figure 4.12 is obtained by employing the kernel, <sup>2</sup> 1 2 *x x* We

2 *b e* . First, the confidence intervals shown in Figure 4.11 are very similar to those ones could obtain by using the SSTM developed in Chapter 2. Comparing the effects of the kernels on the behaviours of the generalized SSTM, we see that the confidence interval bandwidth in Figure 4.12 is almost non-existent. The reason is that the kernel used has a faster convergence when decomposed in the eigen vector space. For smaller values of , the randomness in the concentration realization are minimal but as <sup>2</sup> is increased, we see increased randomness in the realizations. This also allows us to use the kernel used in Figure 4.12 for larger scale computations. We conclude that the choice of the velocity covariance kernel has a significant effect on the behaviour of the generalized SSTM increasing the flexibility of the SSTM.


Table 4.9. Comparison of the dispersivity values for the two kernels.

138 in Porous Media - An Approach Based on Stochastic Calculus

eigen functions could be tedious. The main point in this exercise is to show that any given covariance kernel can be used to obtain the eigen functions of the form given by equation (4.2.3). A corollary to this statement is that if we express the eigen functions in the form given by equation (4.2.3), then we can assume that these is an underlying covariance kernel responsible for these eigen functions. In deriving the SSTM in the form of equation (4.2.15), we assume that the form of equation (4.2.3) is given. But we see now that any covariance-

We have seen in sections 4.2 and 4.3 that the SSTM developed in Chapter 2 and 3 can be recasted so that we could employ any given velocity kernel. In fact, we can even use an empirical set of data for the velocity covariance. We can anticipate that the generalized SSTM would behave quite similar to the one developed in Chapter 2 given that same covariance kernel is used. We compute the 95% confidence intervals for the concentration breakthrough curves (concentration realizations) at *x* = 0.5 m when the flow length is 1 m to compare the differences that occur in using different kernels in the generalized SSTM. The mean velocity is kept constant at 0.5 m/day, and the covariance kernel given by equation (4.3.14) is used to obtain Figure 4.11, and Figure 4.12 is obtained by employing the kernel,

ones could obtain by using the SSTM developed in Chapter 2. Comparing the effects of the kernels on the behaviours of the generalized SSTM, we see that the confidence interval bandwidth in Figure 4.12 is almost non-existent. The reason is that the kernel used has a faster convergence when decomposed in the eigen vector space. For smaller values of

increased randomness in the realizations. This also allows us to use the kernel used in Figure 4.12 for larger scale computations. We conclude that the choice of the velocity covariance kernel has a significant effect on the behaviour of the generalized SSTM

*DL*

the randomness in the concentration realization are minimal but as <sup>2</sup>

1 2

0.0001 0.1 0.02305 0.04611 0.02460 0.04921 0.001 0.1 0.02699 0.05199 0.02513 0.05025 0.01 0.1 0.03059 0.06117 0.02782 0.05361 0.1 0.1 0.06852 0.13705 0.06407 0.12815

*x x*

2

Table 4.9. Comparison of the dispersivity values for the two kernels.

. First, the confidence intervals shown in Figure 4.11 are very similar to those

The Kernel 2

1 2 (,)

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

,

is increased, we see

<sup>2</sup> 1 2

*x x*

2

kernel driven SSTM can be represented by equation (4.2.15).

4.2

**4.4 Effects of Different Kernels and** *<sup>x</sup> h*

increasing the flexibility of the SSTM.

*DL*

The Kernel 1

1 2 (,)

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

 <sup>2</sup> 1 2

*x x b*

 

2

2 *b*

*e* We investigate the effects of the kernels on the dispersivity values; we compute them using the stochastic inverse method (SIM) discussed in Chapter 3. Table 4.9 shows the results. For all practical purposes they are essentially the same. The mechanic of dispersion is more influenced by <sup>2</sup> for a given *b* or if both <sup>2</sup> and *b* are allowed to vary, on both <sup>2</sup> and *b* . The mechanics of dispersion in general can also be assumed to be influenced by the mathematical form of the kernel. The both of these kernels are exponential decaying functions. Because of the case of computations, we continue to use kernel 2 in Table 4.9 in the most of the work discussed in this book.

Figure 4.11. The generalized SSTM 95% confidence intervals for the concentration realization for the kernel, 1 2 2 *x x b e* when <sup>2</sup> = 0.1 and *b* =0.1.

Figure 4.12. The generalized SSTM 95% confidence intervals for the concentration realization for the kernel, <sup>2</sup> 1 2 2 *x x b e* when <sup>2</sup> = 0.1 and *b* =0.1.

1 2

*db t db t*

A Generalized Mathematical Model in One-Dimension 141

.

*dB*

matrix

when,

Klebaner (1998)).

*a*

*m*

*K ij m*

*a* is a symmetric matrix.

*K*

1

for an infinitesimal time increment,

, ,

*x ik x jk*

*x ik*

*db t*

*<sup>m</sup> <sup>m</sup>*

,0 ,1 ,2

*F F F F* 

The drift and diffusion coefficients of the multi-dimensional SDEs are vector *Fx* and the

the

independent of *t* in many cases. Therefore, equation (4.2.4) is a linear multi-dimensional SDE. The associated with this SDE, the matrix *a* called the diffusion matrix can be defined,

> *<sup>T</sup> x x a GG*

when superscript *T* indicates the transposed matrix. Under the conditions such as the coefficients use locally Lipschitz, equation (4.5.4) has strong solutions (see Theorem 6.22 in

The diffusion matrix, *a* , is important to obtain the covariation of *xI* :

if

*G i j*

*GG i j*

if

,

0,1,2, *<sup>i</sup>* and *<sup>j</sup>* 0,1,2.

In obtaining equation (4.2.6), we employ the fact that independent Brownian motions have quadratic covariation. The most important use of quadratic covariation is that we can determine the movement of *I t <sup>x</sup>* with respect to time using the following well known results for Ito diffusions, which are also continuous Markov processes. It can be shown that,

 

. *x x x x*

*Gx* is independent of *t* , and the drift coefficient *Fx* can also be assumed to be

, , *i j i j ij d I I t dI t dI t a dt* (4.5.6)

, , , , *EI t I t F x i x i x i* (4.5.7)

, , , , , , *EI t I t I t I t I t a x i xi xj xj xi ij* (4.5.8)

(4.5.5)

. ,

1

### **4.5 Analysis of Ito Diffusions** *I xt <sup>i</sup>* ,

We have developed the Ito diffusions *I xt <sup>i</sup>* , in section 4.2 (see equation (4.2.20)), and we rewrite the diffusions in the differential form:

$$d\mathbf{l}\_{i}\left(\mathbf{x},t\right) = \mathbf{F}\_{i}\left(\mathbf{x},t\right)dt + \sigma\sum\_{j=1}^{n}\mathbf{G}\_{\left(j\right)}\left(\mathbf{x},t\right)db\_{j}\left(t\right), \text{ for } \left(i=0,1,2\right) \text{ and } \left(j=1,\ldots,m\right). \tag{4.5.1}$$

In equation (4.5.1), *F xt <sup>i</sup>* , are given by equations (4.2.21), (4.2.22) and (4.2.23), which can be considered as regular continuous and differentiable functions of *x* and *t* because we assumed the mean velocity *V xt* , to be a continuous, differentiable function with finite variation in the development of the SSTM.

For many situations, we can assume *V xt* , to be a function of *x* alone, and some regions of *x* it can be considered as a constant. *Gij* 's are continuous differentiable functions of *x* only, with finite variation with respect to *x* (see equations (4.2.24) to (4.2.26)). Therefore, for a fixed value of *x* , we can write equation (4.5.1) as,

$$d\mathcal{I}\_{\boldsymbol{x},i}\left(\mathbf{t}\right) = \mathcal{F}\_{\boldsymbol{x},i}\left(\mathbf{t}\right)dt + \sigma\sum\_{j=1}^{m}\mathcal{G}\_{\boldsymbol{x},j}db\_{j}\left(\mathbf{t}\right), \quad \left(\mathbf{i} = \mathbf{0}, \mathbf{1}, \mathbf{2}\right) \quad \text{and} \quad \left(j = 1, \ldots, m\right). \tag{4.5.2}$$

From the derivation of equation (4.2.19), we see that the *db t s <sup>j</sup>* are the same standard Wiener process increments for each *x i*, *I* and each *db t <sup>j</sup>* : equation (4.5.2) is a diffusion in *m* dimensions, and we can write *I t <sup>x</sup>* as a multidimensional stochastic differential equation (SDE) (Klebaner, 1998).

In coordinate form we can write,

$$d\boldsymbol{I}\_{\boldsymbol{x},i}\left(\boldsymbol{t}\right) = \boldsymbol{F}\_{\boldsymbol{x},i}\left(\boldsymbol{t}\right)dt + \sigma\left[\boldsymbol{\mathcal{G}}\_{\boldsymbol{x},i1} \cdot \boldsymbol{\mathcal{G}}\_{\boldsymbol{x},i2} \; \; \; \; \; \; \; \; \; \mathcal{G}\_{\boldsymbol{x},i\boldsymbol{w}}\right] \cdot \begin{bmatrix} db\_{1}\left(\boldsymbol{t}\right) \\ \boldsymbol{t}db\_{2}\left(\boldsymbol{t}\right) \\ \cdot \\ \cdot \\ \cdot \\ db\_{\boldsymbol{w}}\left(\boldsymbol{t}\right) \end{bmatrix}.\tag{4.5.3}$$

In matrix form, we can write,

$$d\underline{I}\_{\mathbf{u}}(t) = \underline{\mathbf{F}}\_{\mathbf{u}}dt + \sigma \underline{\mathbf{G}}\_{\mathbf{u}} dB(t) \tag{4.5.4}$$

$$\textbf{where}\prime\quad\underline{\mathbf{G}}\_{\underline{\mathbf{x}}}=\begin{bmatrix}\mathbf{G}\_{\underline{\mathbf{x}},01} & \mathbf{G}\_{\underline{\mathbf{x}},02} & \dots & \mathbf{G}\_{\underline{\mathbf{x}},0w} \\ \mathbf{G}\_{\underline{\mathbf{x}},11} & \mathbf{G}\_{\underline{\mathbf{x}},12} & \dots & \mathbf{G}\_{\underline{\mathbf{x}},1w} \\ \mathbf{G}\_{\underline{\mathbf{x}},21} & \mathbf{G}\_{\underline{\mathbf{x}},22} & \dots & \mathbf{G}\_{\underline{\mathbf{x}},2w} \end{bmatrix}\_{3\times w}$$

140 in Porous Media - An Approach Based on Stochastic Calculus

We have developed the Ito diffusions *I xt <sup>i</sup>* , in section 4.2 (see equation (4.2.20)), and we

In equation (4.5.1), *F xt <sup>i</sup>* , are given by equations (4.2.21), (4.2.22) and (4.2.23), which can be considered as regular continuous and differentiable functions of *x* and *t* because we assumed the mean velocity *V xt* , to be a continuous, differentiable function with finite

For many situations, we can assume *V xt* , to be a function of *x* alone, and some regions of *x* it can be considered as a constant. *Gij* 's are continuous differentiable functions of *x* only, with finite variation with respect to *x* (see equations (4.2.24) to

,

From the derivation of equation (4.2.19), we see that the *db t s <sup>j</sup>* are the same standard Wiener process increments for each *x i*, *I* and each *db t <sup>j</sup>* : equation (4.5.2) is a diffusion in *m* dimensions, and we can write *I t <sup>x</sup>* as a multidimensional stochastic differential

, , , 1 , 2 .... , . .

*dI t F dt G dB t x x* 

*xi xi x i x i x im*

*dI t F t dt G G G*

,

(4.2.26)). Therefore, for a fixed value of *x* , we can write equation (4.5.1) as,

1

*j*

*m*

, , ,

*xi xi x ij j*

,01 ,02 .... ,0 ,11 ,12 .... ,1 ,21 ,22 .... ,2 <sup>3</sup>

*xxx x m*

*GGG G*

*GG G*

*x x x m*

*x x x m <sup>m</sup>*

*GG G* 

*dI t F t dt G db t* 

for *<sup>i</sup>* 0,1,2 and *j m* 1,..., . (4.5.1)

*<sup>i</sup>* 0,1,2 and *j m* 1,..., . (4.5.2)

 

1 2

*db t db t*

.

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

(4.5.3)

*m*

*db t*

**4.5 Analysis of Ito Diffusions** *I xt <sup>i</sup>* ,

rewrite the diffusions in the differential form:

variation in the development of the SSTM.

equation (SDE) (Klebaner, 1998). In coordinate form we can write,

In matrix form, we can write,

when,

 1 , , , , *m*

*i i ij j j dI x t F x t dt G x t db t* 

$$\underline{dB} = \begin{bmatrix} db\_1(t) \\ db\_2(t) \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot b\_{\pi^\*}(t) \end{bmatrix}\_{\pi \times 1} \cdot \int\_{\pi \times 1} \cdot \int\_{\pi \times 1} \cdot \frac{1}{\cdot} \, d\theta\_{\pi^\*} \, d\theta\_{\pi^\*} $$
 
$$\underline{F}\_{\pi} = \begin{bmatrix} F\_{x,0} \\ F\_{x,1} \\ F\_{x,2} \end{bmatrix} .$$

The drift and diffusion coefficients of the multi-dimensional SDEs are vector *Fx* and the matrix *Gx* is independent of *t* , and the drift coefficient *Fx* can also be assumed to be independent of *t* in many cases. Therefore, equation (4.2.4) is a linear multi-dimensional SDE. The associated with this SDE, the matrix *a* called the diffusion matrix can be defined,

$$\underline{a} = \left(\sigma \underline{\underline{G}}\_{\times}\right) \left(\sigma \underline{\underline{G}}\_{\times}\right)^{\mathrm{r}} \tag{4.5.5}$$

when superscript *T* indicates the transposed matrix. Under the conditions such as the coefficients use locally Lipschitz, equation (4.5.4) has strong solutions (see Theorem 6.22 in Klebaner (1998)).

The diffusion matrix, *a* , is important to obtain the covariation of *xI* :

$$\mathrm{d}\left[\underset{}{\mathrm{I}}I\_{i\prime}\,\mathrm{I}\_{\,\,\,\,\,}\right](t) = \mathrm{d}I\_{i}\,\mathrm{d}\,\,\mathrm{I}\_{\,\,\,}\left(\mathrm{t}\right) = a\_{i\prime}dt,\tag{4.5.6}$$

when, 2 2 , 1 2 , , 1 if , if *m x ik K ij m x ik x jk K G i j a GG i j* 0,1,2, *<sup>i</sup>* and *<sup>j</sup>* 0,1,2.

### *a* is a symmetric matrix.

In obtaining equation (4.2.6), we employ the fact that independent Brownian motions have quadratic covariation. The most important use of quadratic covariation is that we can determine the movement of *I t <sup>x</sup>* with respect to time using the following well known results for Ito diffusions, which are also continuous Markov processes. It can be shown that, for an infinitesimal time increment,

$$E\left(I\_{x,i}\left(t+\Delta\right) - I\_{x,i}\left(t\right)\right) = F\_{x,i}\Delta + \mathcal{O}\left(\Delta\right),\tag{4.5.7}$$

$$E\left\langle \left( I\_{\mathbf{x},i}\left(\mathbf{t}+\boldsymbol{\Delta}\right) - I\_{\mathbf{x},i}\left(\mathbf{t}\right) \right) \left( I\_{\mathbf{x},j}\left(\mathbf{t}+\boldsymbol{\Delta}\right) - I\_{\mathbf{x},j}\left(\mathbf{t}\right) \right) \middle| I\_{\mathbf{x},i}\left(\mathbf{t}\right) \right\rangle = a\_{\boldsymbol{\eta}}\boldsymbol{\Delta} + \mathcal{O}\left(\boldsymbol{\Delta}\right),\tag{4.5.8}$$

Equation (4.5.11) provides realizations which have the statistical properties of the realizations depicted in Figure 4.13. The vector *B t* consists of independent Wiener

A Generalized Mathematical Model in One-Dimension 143

The statistical properties of these realizations are essentially the same to those of the

As we have mentioned earlier, *I t <sup>x</sup>* are only dependent on the velocity patterns in the medium, and it is important ask the question how the correlation length, *b* , affects the

diffusion form of the SDE given in equation (4.5.11). However, the correlation length *b* influence *I t <sup>x</sup>* nonlinearly through *Gx* , but this influence can always be captured by

the porous medium under study. We found that *b* =0.1 is suitable for our computational

Figure 4.14. Some realization of *I t <sup>x</sup>* for fixed *F s x i*, when 0.5 *x* based on equation (4.5.11).

*m*

*j I t F t G Bt* 

> ,1 ,1 ,1 1

*x x xjj*

*m*

*j I t F t G Bt* 

,

(4.5.12a)

(4.5.12b)

,

,0 ,0 ,0 1

*x x xjj*

(the square root of the variance of the kernel) acts as the multiplication factor to the

. Therefore, we can keep *b* at a constant value that is appropriate for

2 1 2

only). It is seen

2 *x x b*

 *e* 

processes; Figure 4.14 shows some realizations based on equation (4.5.11).

realization of *I t <sup>x</sup>* . (In this discussion we focus on the kernel

experiments in this chapter as well as in chapters 6, 7 and 8.

Equation (4.5.11) can be written in component forms,

realizations given in Figure 4.13.

that 

suitable changes in

and when *i j* , equation (4.6.8) becomes,

$$E\left|\left(I\_{\mathbf{x},i}\left(t+\Delta\right)-I\_{\mathbf{x},i}\left(t\right)\right)^{2}\Big|I\_{\mathbf{x},i}\left(t\right)\right|=a\_{i}\Delta+\mathcal{O}\left(\Delta\right),\tag{4.5.9}$$

As the solution to the SDE equation (4.5.4) exists, from equation (4.5.8) it is seen that *Fx* is the form of *xI* at time *t* , and *a* is the coefficients in the covariance of the infinitesimal displacement from *xI* . Using those results we can construct the realizations of *I t <sup>x</sup>* by using the fact that *I t <sup>x</sup>* are Gaussian processes. By dividing a given time interval into equidistant infinitesimal time interval, , we can generate normally distributed *<sup>x</sup> dI* increments for a given *x* , using the mean and variance obtained by equations (4.5.7) and (4.5.9). It should be noted that in generating the standard Wiener process increments, we use the zero-mean and -variance Gaussian increments.

We take 0 *xI t* when 0 *t* because

$$
\underline{\mathbf{U}}\_{\star}(t) = \int\_{0}^{t} \underline{\mathbf{F}}\_{\star}(t)dt + \sigma \int\_{0}^{t} \underline{\mathbf{G}}\_{\times}(t)d\underline{\mathbf{B}}(t). \tag{4.5.10}
$$

Figure 4.13 shows some realizations of *I t <sup>x</sup>* when 0.5 *x* .

Figure 4.13. Some realizations of *I t <sup>x</sup>* when 0.5 *x* .

The increments of *I t <sup>x</sup>* are Gaussian random variables having the mean and variance given by equations (4.5.7) and (4.5.8).

The SDE given in equation (4.5.4) is linear and strong solutions do exist. When *Fx i*, are not functions of *t* , then the solution of equation (4.5.4) is given by

$$
\underline{I}\_{\times}(t) = \underline{F}\_{\times}t + \sigma \underline{G}\_{\times} \underline{B}(t). \tag{4.5.11}
$$

(4.5.10)

(4.5.11)

, , , , *EI t I t I t a x i xi xi ii* (4.5.9)

142 in Porous Media - An Approach Based on Stochastic Calculus

<sup>2</sup>

As the solution to the SDE equation (4.5.4) exists, from equation (4.5.8) it is seen that *Fx* is the form of *xI* at time *t* , and *a* is the coefficients in the covariance of the infinitesimal displacement from *xI* . Using those results we can construct the realizations of *I t <sup>x</sup>* by using the fact that *I t <sup>x</sup>* are Gaussian processes. By dividing a given time interval into equidistant infinitesimal time interval, , we can generate normally distributed *<sup>x</sup> dI* increments for a given *x* , using the mean and variance obtained by equations (4.5.7) and (4.5.9). It should be noted that in generating the standard Wiener process increments, we use

exists,

 <sup>0</sup> <sup>0</sup> . *<sup>t</sup> <sup>t</sup> <sup>x</sup> <sup>x</sup> <sup>x</sup> I t F t dt G t dB t* 

The increments of *I t <sup>x</sup>* are Gaussian random variables having the mean and variance

The SDE given in equation (4.5.4) is linear and strong solutions do exist. When *Fx i*, are not

 . *x xx I t Ft G B t* 

and when *i j* , equation (4.6.8) becomes,

the zero-mean and -variance Gaussian increments.

Figure 4.13 shows some realizations of *I t <sup>x</sup>* when 0.5 *x* .

Figure 4.13. Some realizations of *I t <sup>x</sup>* when 0.5 *x* .

functions of *t* , then the solution of equation (4.5.4) is given by

given by equations (4.5.7) and (4.5.8).

We take 0 *xI t* when 0 *t* because

Equation (4.5.11) provides realizations which have the statistical properties of the realizations depicted in Figure 4.13. The vector *B t* consists of independent Wiener processes; Figure 4.14 shows some realizations based on equation (4.5.11).

The statistical properties of these realizations are essentially the same to those of the realizations given in Figure 4.13.

As we have mentioned earlier, *I t <sup>x</sup>* are only dependent on the velocity patterns in the medium, and it is important ask the question how the correlation length, *b* , affects the realization of *I t <sup>x</sup>* . (In this discussion we focus on the kernel 2 1 2 2 *x x b e* only). It is seen that (the square root of the variance of the kernel) acts as the multiplication factor to the diffusion form of the SDE given in equation (4.5.11). However, the correlation length *b* influence *I t <sup>x</sup>* nonlinearly through *Gx* , but this influence can always be captured by suitable changes in . Therefore, we can keep *b* at a constant value that is appropriate for the porous medium under study. We found that *b* =0.1 is suitable for our computational experiments in this chapter as well as in chapters 6, 7 and 8.

Figure 4.14. Some realization of *I t <sup>x</sup>* for fixed *F s x i*, when 0.5 *x* based on equation (4.5.11). Equation (4.5.11) can be written in component forms,

$$H\_{\mathbf{x},0}\left(\mathbf{t}\right) = F\_{\mathbf{x},0}\left(\mathbf{t}\right) + \sigma \sum\_{j=1}^{m} G\_{\mathbf{x},0} B\_{j}\left(\mathbf{t}\right),\tag{4.5.12a}$$

$$H\_{\mathbf{x},1}\left(\mathbf{t}\right) = F\_{\mathbf{x},1}\left(\mathbf{t}\right) + \sigma \sum\_{j=1}^{m} G\_{\mathbf{x},1/} B\_{j}\left(\mathbf{t}\right),\tag{4.5.12b}$$

Figure 4.16. The approximated *P xs ij* given by equations (4.2.10), (4.2.11) and (4.2.12)

A Generalized Mathematical Model in One-Dimension 145

When 0.1 *b* , the required number of significant eigen values is reduced to 6 and *Pij* functions are depicted in Figure 4.16. Similar observation as before, when 0.05 *b* , can be made. Figure 4.17 shows *Pij* functions when 0.2 *b* , and now the number of significant eigen values is 5 (i.e, 5 *m* ). The same observations can be made for *Pij* s when 0.2 *b* .

We produce the 3-dimensional graphs of *Pij* when 0,1 *i* and *j* 1,2,3,4,5 in Figure 4.18 and Figure 4.19; *b* is plotted as the y-axis. As one could expect, it is reasonable to assume that function surface of *Pij* is a smooth, continuous function of *b* . We can define

when 0.0 1.0 *x* for 0.1 *b* .

continuous functions of *x* and *b* to define *P s* <sup>2</sup> *<sup>j</sup>* .

$$\mathbf{I} \text{ and } \ I\_{\mathbf{x},2}\left(t\right) = F\_{\mathbf{x},2}\left(t\right) + \sigma \sum\_{j=1}^{m} G\_{\mathbf{x},2j} B\_{j}\left(t\right). \tag{4.5.12c}$$

In the above equations, *m* is the number of significant eigen functions and is dependent on *b* . From equation (4.2.24), (4.2.25) and (4.2.26), we see that *Gij* are related to *P x ij* which are given by equations (4.2.10), (4.2.11) and (4.2.12). For 0.05 *b* , 8 *m* , Figure 4.15 give *P xs ij* when 0.0 1.0 *x* . The following observations can be noted from Figure 4.15: (a) *P sij* , therefore *G sij* are sinesodial in nature; (b) amplitutes of *P sij* increase with *m* (eigen function number) but as eigen values decrease with *m* , *G sij* diminish with *m* (not shown); (c) *P s* <sup>2</sup> *<sup>j</sup>* are insignificant in comparision to *P s* <sup>0</sup> *<sup>j</sup>* and *P s* <sup>1</sup> *<sup>j</sup>* and therefore could be ignored; and (d) frequency of *Pij* functions increases as *m* increase. 

Figure 4.15. The approximated *P xs ij* given by equations (4.2.10), (4.2.11) and (4.2.12) when 0.0 1.0 *x* for 0.05 *b* .

.

(4.5.12c)

144 in Porous Media - An Approach Based on Stochastic Calculus

and ,2 ,2 ,2

*x x xjj*

In the above equations, *m* is the number of significant eigen functions and is dependent on *b* . From equation (4.2.24), (4.2.25) and (4.2.26), we see that *Gij* are related to *P x ij* which are given by equations (4.2.10), (4.2.11) and (4.2.12). For 0.05 *b* , 8 *m* , Figure 4.15 give *P xs ij* when 0.0 1.0 *x* . The following observations can be noted from Figure 4.15: (a) *P sij* , therefore *G sij* are sinesodial in nature; (b) amplitutes of *P sij* increase with *m* (eigen function number) but as eigen values decrease with *m* , *G sij* diminish with *m* (not shown); (c) *P s* <sup>2</sup> *<sup>j</sup>* are insignificant in comparision to *P s* <sup>0</sup> *<sup>j</sup>* and *P s* <sup>1</sup> *<sup>j</sup>* and therefore could be

(4.2.26),

Figure 4.15. The approximated *P xs ij* given by equations (4.2.10), (4.2.11) and (4.2.12)

when 0.0 1.0 *x* for 0.05 *b* .

ignored; and (d) frequency of *Pij* functions increases as *m* increase.

1

*j I t F t G Bt* 

*m*

Figure 4.16. The approximated *P xs ij* given by equations (4.2.10), (4.2.11) and (4.2.12) when 0.0 1.0 *x* for 0.1 *b* .

When 0.1 *b* , the required number of significant eigen values is reduced to 6 and *Pij* functions are depicted in Figure 4.16. Similar observation as before, when 0.05 *b* , can be made. Figure 4.17 shows *Pij* functions when 0.2 *b* , and now the number of significant eigen values is 5 (i.e, 5 *m* ). The same observations can be made for *Pij* s when 0.2 *b* .

We produce the 3-dimensional graphs of *Pij* when 0,1 *i* and *j* 1,2,3,4,5 in Figure 4.18 and Figure 4.19; *b* is plotted as the y-axis. As one could expect, it is reasonable to assume that function surface of *Pij* is a smooth, continuous function of *b* . We can define continuous functions of *x* and *b* to define *P s* <sup>2</sup> *<sup>j</sup>* .

Figure 4.17. The approximated *P xs ij* given by equations (4.2.10), (4.2.11) and (4.2.12) when 0.0 1.0 *x* for 0.2 *b* .

Figure 4.18. The 3-dimensional graphs of *Pij* when 0 *i* and *j* 1,2,3,4,5

A Generalized Mathematical Model in One-Dimension 147

146 in Porous Media - An Approach Based on Stochastic Calculus

Figure 4.17. The approximated *P xs ij* given by equations (4.2.10), (4.2.11) and (4.2.12)

when 0.0 1.0 *x* for 0.2 *b* .

Figure 4.18. The 3-dimensional graphs of *Pij* when 0 *i* and *j* 1,2,3,4,5 of

Figure 4.19. The 3-dimensional graphs of

1

where

when,

*n*

*dI t <sup>x</sup>* , and therefore can be written as,

*x x*

*x*

*x*

2 2 .

*x*

*x*

*C C F C*

 

*x*

Therefore, we can write,

*Pij* when 1 *i* and

*x*

As we recall that the SSTM can be expressed as a diffusion process with martingale

A Generalized Mathematical Model in One-Dimension 149

,0 ,1

*C <sup>C</sup> dC t C t dI dI dI <sup>x</sup> <sup>x</sup>*

We need to interpret this equation as follows: the solute concentration at a given point

consists of a combination of three diffusion processes which are solely based on velocity.

the combination prevailing

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

*C*

*x xn xn x n x n x n*

*x*

 

*C C dC t C dI FdI x x*

> 0*<sup>t</sup> C t FdI t <sup>x</sup> <sup>x</sup>* .

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

,

*n*

*t* denotes the discretized time and the spatial gradients act as the coefficients of

2

*x*

The spatial influence on the concentration is mediated through the prevailing concentration, and its spatial gradients. In keeping with the Ito definition of stochastic integration (Klebaner, 1998) we must use the concentration and its spatial gradients at a previous time.

*x*

*x*

properties in time dimensions when another set of diffusion processes, *I t <sup>x</sup>*

As a difference equation, we can write equation (4.5.13) as,

*x x x*

*j*

2 ,2 .

*x*

*x*

2 ,2

*n*

,

, (4.5.15)

,

*x*

(4.5.14)

2

*x t x t*

,0

*x*

*dI*

*dI*

*x*

*x*

,2

2 ,1

(4.5.13)

*x*

1,2,3,4,5

are used:

*x*

148 in Porous Media - An Approach Based on Stochastic Calculus

Figure 4.19. The 3-dimensional graphs of *Pij* when 1 *i* and *j* 1,2,3,4,5

As we recall that the SSTM can be expressed as a diffusion process with martingale properties in time dimensions when another set of diffusion processes, *I t <sup>x</sup>* are used:

$$d\mathbf{C}\_{\mathbf{x}}\left(t\right) = -\mathbf{C}\_{\mathbf{x}}\left(t\right)d\mathbf{I}\_{\mathbf{x},0} - \left(\frac{\partial \mathbf{C}\_{\mathbf{x}}}{\partial \mathbf{x}}\right)\_{\mathbf{x}}d\mathbf{I}\_{\mathbf{x},1} - \left(\frac{\partial^2 \mathbf{C}\_{\mathbf{x}}}{\partial \mathbf{x}^2}\right)\_{\mathbf{x}}d\mathbf{I}\_{\mathbf{x},2}.\tag{4.5.13}$$

We need to interpret this equation as follows: the solute concentration at a given point *x* consists of a combination of three diffusion processes which are solely based on velocity.

The spatial influence on the concentration is mediated through the prevailing concentration, and its spatial gradients. In keeping with the Ito definition of stochastic integration (Klebaner, 1998) we must use the concentration and its spatial gradients at a previous time. As a difference equation, we can write equation (4.5.13) as, follows: definition

$$\mathbf{C}\_{x}\left(t+1\right) - \mathbf{C}\_{x}\left(t\_{n}\right) = -\mathbf{C}\_{x}\left(t\_{n}\right)d\mathbf{I}\_{x,0}\left(t\_{n}\right) - \left(\frac{\partial \mathbf{C}}{\partial \mathbf{x}}\right)\_{\mathbf{x},t\_{n}}d\mathbf{I}\_{x,1}\left(t\_{n}\right) - \left(\frac{\partial^{2} \mathbf{C}}{\partial \mathbf{x}^{2}}\right)\_{\mathbf{x},t\_{n}}d\mathbf{I}\_{x,2}\left(t\_{n}\right),\tag{4.5.14}$$

where *nt* denotes the discretized time and the spatial gradients act as the coefficients of *dI t <sup>x</sup>* , and therefore can be written as,

$$d\mathbf{C}\_{x}(t) = \begin{bmatrix} -\mathbf{C}\_{x} & -\frac{\partial \mathbf{C}\_{x}}{\partial \mathbf{x}} & -\frac{\partial^{2} \mathbf{C}\_{x}}{\partial \mathbf{x}^{2}} \end{bmatrix} \begin{bmatrix} d\mathbf{I}\_{x,0} \\ d\mathbf{I}\_{x,1} \\ d\mathbf{I}\_{x,2} \end{bmatrix} = \underline{\mathbf{F}} d\underline{\mathbf{I}}\_{x}. \tag{4.5.15}$$
  $\underline{\mathbf{F}} = \begin{bmatrix} -\mathbf{C}\_{x} & -\frac{\partial \mathbf{C}\_{x}}{\partial \mathbf{x}} & -\frac{\partial^{2} \mathbf{C}\_{x}}{\partial \mathbf{x}^{2}} \end{bmatrix}.$ 

Therefore, we can write,

when,

$$C\_x
\left(t
\right) = 
\int\_o \underline{F} dI\_x
\left(t
\right).$$

As *P*<sup>2</sup> *<sup>j</sup>* functions are insignificant compared to the other two functions, we can approximate equation (4.5.15) as

$$d\mathbf{C}\_{x}\left(t\right) = \begin{bmatrix} -\mathbf{C}\_{x} & -\frac{\partial \mathbf{C}\_{x}}{\partial \mathbf{x}} \end{bmatrix} \begin{bmatrix} dI\_{x,0} \\ dI\_{x,1} \end{bmatrix}.\tag{4.5.16}$$

Similarly,

1

 

1 1

 

*k k*

*m m*

*k*

22 1

*ij a s* can also be evaluated using equation (4.5.24).

calculate the strong solutions for *I t <sup>x</sup>* .

The propagator of *I t <sup>x</sup>* can be defined as

Markov process; the propagator

To abbreviate the notation, we denote the propagator as

propagator function can be expressed as follows (Gillespie, 1992):

th propagator moment function.) ; and *P t*

**4.6 Propagator of** *I t <sup>x</sup>*

*i* 

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

moments of

*m*

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

*k k k k*

*x dx dx x*

2 .

2

(4.5.27)

(4.5.25)

2 .

*tI t t I t t I I t* (4.6.1)

(4.6.2)

and (4.6.3)

(4.6.4)

*i xi tI t t* ; , , is also a Markov process having a

is the probability density function of

*<sup>i</sup>t* . It can be shown that the

(4.5.26)

*k k*

 

2 2 11 2 2

*<sup>f</sup> df d f <sup>f</sup> <sup>a</sup> <sup>h</sup> <sup>h</sup>*

222 2 2

*df df a G f hf h dx dx*

33

*a*

*<sup>i</sup> t* has the analytical structure (Gillespie,1992),

 

*k x k x k*

*k kk xkk x k*

1

This menas that once we have eigen functions in the form given by equation (4.2.3), the diffusion matrix *a* can be evaluated and equation (4.5.7) and (4.5.9) can be used to

; , , , , , . *i xi x i xi xi*

 denotes the displacement of *I t x i*, during the infinitesimally small time interval *t* given the position (value) of *I t x i*, at time *t* . As *I t x i*, is an Ito diffusion, it is also a

Gaussian probability density function, which completely specifics the propagator variable.

, *<sup>n</sup> <sup>n</sup> E t d P t Bt t t <sup>i</sup> <sup>u</sup>*

when *Bn* s are well behaved functions of *t* for a given *I t x i*, ( *B t <sup>n</sup>* are also called the *n*

It can be shown that, for a given value of *I t x i*, at time *t* , the mean and variance of the

, *E t At t t <sup>i</sup> <sup>i</sup>*

, *Var t D t t t <sup>i</sup> <sup>i</sup>*

*k*

*m x j j*

 

2 2

*h*

. <sup>4</sup>

 *f*

A Generalized Mathematical Model in One-Dimension 151

This approximation enables us to reduce the computational time further in obtaining solutions. obtaining (4.2.3)

Equation (4.2.3) gives the general form of eigen functions, *f x <sup>j</sup>* s. By inspecting equation (4.2.10), (4.2.11) and (4.2.12), we can deduce the following relationships between *Pij* s and *<sup>j</sup> f* s:

$$P\_{0\_{/}}(\mathbf{x}) = \frac{df\_{/}(\mathbf{x})}{d\mathbf{x}} + h\_{\mathbf{x}} \frac{d^2 f\_{/}(\mathbf{x})}{d\mathbf{x}^2},\tag{4.5.17}$$

$$P\_{1/}\left(\mathbf{x}\right) = f\_i\left(\mathbf{x}\right) + h\_{\mathbf{x}} \frac{df\_i\left(\mathbf{x}\right)}{d\mathbf{x}},\tag{4.5.18}$$

$$\text{and}\quad P\_{2\nearrow}(\mathbf{x}) = \frac{h\_{\times}}{2} f\_{\rangle}(\mathbf{x}). \tag{4.5.19}$$

Therefore *Gij* are given by,

$$\mathbf{G}\_{0j} = \sqrt{\mathbf{\hat{\lambda}}\_{j}} \frac{d f\_{\rangle}}{d \mathbf{x}} + h\_{\mathbf{x}} \sqrt{\mathbf{\hat{\lambda}}\_{j}} \frac{d^2 f\_{\rangle}}{d \mathbf{x}^2} \,' \,. \tag{4.5.20}$$

$$\mathbf{G}\_{1j} = \sqrt{\mathcal{k}\_j} f\_j + h\_{\ge 1} \sqrt{\mathcal{k}\_j} \frac{df\_j}{d\mathbf{x}} \,\, ^\prime \tag{4.5.21}$$

$$\text{and}\quad G\_{2/} = \frac{h\_x}{2}\sqrt{\mathcal{k}\_{\parallel}}f\_{\cdot}.\tag{4.5.22}$$

The components of diffusion matrix *a* (see equation (4.5.5)) can be written as,

$$\mathfrak{a}\_{il} = \sigma^2 \sum\_{k=1}^{m} \mathbb{G}^2\_{(l-1)k'} \qquad\qquad \text{(i,1,2,3)}\,,\tag{4.5.23}$$

$$\mathbf{a} \text{ and } \ a\_{\boldsymbol{\psi}} = \sigma^2 \sum\_{k=1}^{m} \mathbf{G}\_{l-1k} \mathbf{G}\_{\boldsymbol{\phi}-1k\boldsymbol{\nu}} \qquad \qquad \left(\boldsymbol{j}, \mathbf{1}, \mathbf{2}, \mathbf{3}\right) \; . \tag{4.5.24}$$

For example,

$$\begin{aligned} a\_{11} &= \sigma^2 \sum\_{k=1}^m G\_{0k}^2 \\ &= \sigma^2 \sum\_{k=1}^m \left( \sqrt{\lambda\_k} \, \frac{df\_k}{d\boldsymbol{\chi}} + h\_\boldsymbol{x} \sqrt{\lambda\_k} \, \frac{d^2 f\_k}{d\boldsymbol{\chi}^2} \right)^2 \end{aligned}$$

$$\mathfrak{a}\_{11} = \sigma^2 \sum\_{k=1}^{m} \left( \lambda\_k \left( \frac{\partial f\_k}{\partial \mathbf{x}} \right)^2 + 2h\_\mathbf{x} \lambda\_k \left( \frac{df\_k}{d\mathbf{x}} \right) \frac{d^2 f\_k}{d\mathbf{x}^2} + h\_\mathbf{x}^2 \lambda\_k \left( \frac{\partial^2 f\_k}{\partial \mathbf{x}^2} \right)^2 \right). \tag{4.5.25}$$

Similarly,

Computational Modelling of Multi-Scale Non-Fickian Dispersion

(4.5.17)

(4.5.18)

*<sup>h</sup> P x fx* (4.5.19)

(4.5.20)

(4.5.21)

(4.5.22)

. (4.5.16)

,1

*x*

*x x*

*C dI*

*x dI*

150 in Porous Media - An Approach Based on Stochastic Calculus

As *P*<sup>2</sup> *<sup>j</sup>* functions are insignificant compared to the other two functions, we can

,0

This approximation enables us to reduce the computational time further in obtaining

Equation (4.2.3) gives the general form of eigen functions, *f x <sup>j</sup>* s. By inspecting equation (4.2.10), (4.2.11) and (4.2.12), we can deduce the following relationships between *Pij* s and *<sup>j</sup> f* s:

> <sup>2</sup> <sup>0</sup> <sup>2</sup> , *<sup>j</sup> <sup>j</sup> j x df x d f x P x <sup>h</sup>*

 <sup>1</sup> , *<sup>j</sup> j i x df x P x fx h dx*

<sup>0</sup> <sup>2</sup> , *<sup>j</sup> <sup>j</sup> j j x j df d f <sup>G</sup> <sup>h</sup>*

<sup>1</sup> ,

 and 2 . <sup>2</sup> *x j j j <sup>h</sup> G f* 

and <sup>2</sup> 1 1

*j*

*j jj x j df G fh dx* 

 

The components of diffusion matrix *a* (see equation (4.5.5)) can be written as,

 2 2 1 1

1

2 2 11 0 1

*k*

*m*

1

*k*

2

*k m*

*a G*

*m ij ikjk k a GG* 

*i*

*m ii i k k a G* 

*dx dx*

 

 

, ( ,1,2,3)

*dx dx*

*x j j*

2

*j*

, (4.5.23)

. (4.5.24)

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

,

2

*k k k x k*

 

*df d f <sup>h</sup> dx dx*

, ,1,2,3

*x x*

and 2 . <sup>2</sup>

*dC t C*

approximate equation (4.5.15) as

Therefore *Gij* are given by,

For example,

solutions.

$$\mathfrak{a}\_{22} = \sigma^2 \sum\_{k=1}^{m} \mathcal{G}\_{1k}^2 = \sigma^2 \sum\_{k=1}^{m} \left( \mathcal{A}\_k f\_k^2 + 2 \text{h}\_x \mathcal{A}\_k f\_k \frac{df\_k}{d\mathbf{x}} + \text{h}\_x^2 \mathcal{A}\_k \left(\frac{df\_k}{d\mathbf{x}}\right)^2 \right). \tag{4.5.26}$$

$$a\_{33} = \sigma^2 \sum\_{k=1}^{m} \frac{h\_x}{4} \mathcal{A}\_j f\_{\rangle}^2. \tag{4.5.27}$$

*ij a s* can also be evaluated using equation (4.5.24).

This menas that once we have eigen functions in the form given by equation (4.2.3), the diffusion matrix *a* can be evaluated and equation (4.5.7) and (4.5.9) can be used to calculate the strong solutions for *I t <sup>x</sup>* .

### **4.6 Propagator of** *I t <sup>x</sup>*

The propagator of *I t <sup>x</sup>* can be defined as

$$\left(\mathcal{L}\left(\Delta t; I\_{x,i}\left(t\right), t\right) \equiv \left(I\_{x,i}\left(t + \Delta t\right) - I\_{x,i}\right)\right|I\_{x,i}\left(t\right). \tag{4.6.1}$$

*i* denotes the displacement of *I t x i*, during the infinitesimally small time interval *t* given the position (value) of *I t x i*, at time *t* . As *I t x i*, is an Ito diffusion, it is also a Markov process; the propagator *i xi tI t t* ; , , is also a Markov process having a Gaussian probability density function, which completely specifics the propagator variable. To abbreviate the notation, we denote the propagator as *<sup>i</sup>t* . It can be shown that the moments of *<sup>i</sup> t* has the analytical structure (Gillespie,1992),

$$E\left(\boldsymbol{\varepsilon}\_{i}^{\boldsymbol{u}}\left(\Delta t\right)\right) \equiv \int\_{-\infty}^{\boldsymbol{u}} d\boldsymbol{\varepsilon} \boldsymbol{x}^{\boldsymbol{u}} P\left(\boldsymbol{\varepsilon} \,|\,\Delta t\right) = B\_{\boldsymbol{u}}\left(t\right)\Delta t + \mathcal{O}\left(\Delta t\right),\tag{4.6.2}$$

when *Bn* s are well behaved functions of *t* for a given *I t x i*, ( *B t <sup>n</sup>* are also called the *n* th propagator moment function.) ; and *P t* is the probability density function of *<sup>i</sup> t* .

It can be shown that, for a given value of *I t x i*, at time *t* , the mean and variance of the propagator function can be expressed as follows (Gillespie, 1992):

$$E\left(\varepsilon\_{i}\left(\Delta t\right)\right) = A\_{i}\left(t\right)\Delta t + O\left(\Delta t\right), \text{ and} \tag{4.6.3}$$

$$\operatorname{Var}\left(\varepsilon\_{\!\!\! }\left(\Delta t\right)\right) = D\_{\!\!\! }\left(\mathbf{t}\right)\Delta t + \operatorname{O}\left(\Delta t\right),\tag{4.6.4}$$

A Generalized Mathematical Model in One-Dimension 153

This is the first form of the Langevin equation for *I t x i*, .

*dw t* . Therefore, we can rewrite equation (4.7.3) as,

*d Var I t*

be simplified to equations (4.5.7) and (4.5.8).

The SSTM can be written as, for given *x* , (see equation (4.5.13)),

<sup>2</sup>

*<sup>C</sup> <sup>C</sup> dC t C t dI dI dI*

*x x x x x*

**4.8 The Evaluation of** *C t <sup>x</sup>* **Diffusions**

*dt*

Equation (4.7.5) is the second form of the Langevin equation.

 

It can be shown that, if *N Normal* 0,1 ,

*I t x i*, s more efficiently.

with initial conditions,

for derivations.)

2

<sup>2</sup>

 

Therefore, *dtN* is a normally distributed random variable having the density function *Normal dt* 0, . This is the density function of the standard Wiener process increments,

The advantage of this Langevin approximation for *I t x i*, over equation (4.5.2) is that it is not a multidimensional SDE but a one-dimensional one. This would allow us to compute

The moment evolution equations for *I t x i*, can be given as follows: (See Gillespie (1992)

, , *x i*

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

If *Fx i*, and *ii a* are deterministic functions of *x* then the moment evolution equations can

,0 ,1 2 ,2 *x x*

*x x*

*x x*

, (4.8.1)

*x i*

*E I tF E I t EF Ea*

*xi xi x i x i ii*

(4.7.7)

*E F*

,

*dEI t*

*dt*

, , . *x i x i ii dI t F dt a t dt N* (4.7.3)

*N Normal* , , (4.7.4)

, , . *xi xi ii i dI t F dt a dw t* (4.7.5)

(4.7.6)

*E I x i*, 0 0, and (4.7.8)

*Var I x i*, 0 0. (4.7.9)

when *Ai t* and *D t <sup>i</sup>* are independent of *t* . As the variance can not be negative, *D t <sup>i</sup>* should be a non-negative function. As mentioned before, *<sup>i</sup>t* is a Gaussian variable, i.e., it is normally distributed; therefore, we can write,

$$
\varepsilon\_i(\Delta t) = N \left( A\_i(t) \Delta t, D\_i(t) \Delta t \right), \tag{4.6.5}
$$

where " " indicates the random variable is generated from the normal density function having the mean, *At t <sup>i</sup>* and the variance *Dt t <sup>i</sup>* .

If we recall, equation (4.5.9) gives,

$$E\left|\left(I\_{\mathbf{x},i}\left(t+\Delta\right)-I\_{\mathbf{x},i}\left(t\right)\right)^{2}\Big|I\_{\mathbf{x},i}\left(t\right)\right|=a\_{i}\Delta t+O\Big(\Delta\big).\tag{4.6.6}$$

This can be written as,

$$E\left\{\varepsilon\_i^2\left(\Delta t\right)\middle|I\_{x,i}\left(t\right)\right\} = a\_i \Delta t. \tag{4.6.7}$$

when , *t t* i.e., 0, *t* therefore 0 *<sup>i</sup> t* .

$$E(\varepsilon\_i \mid \Delta t) = 0 \quad \text{at} \quad t = t\_\prime \dots$$

This leads to

$$\operatorname{Var}\left(\varepsilon\_{i}\left(\Delta t\right)\right) = a\_{i\cdot}\Delta + \operatorname{O}\left(\Delta\right). \tag{4.6.8}$$

By comparing equation (4.5.7) with equation (4.6.3), and equation (4.6.4) with equation (4.6.8), we deduce that,

$$A\_i(t) = F\_{x\_i \nu} \quad \text{and} \tag{4.6.9}$$

$$D\_i(t) = a\_{ii} \,\tag{4.6.10}$$

We can now write the Langevin equation for the Ito diffusion, *I t x i*, .

### **4.7 Langevin Equation for** *I t x i*,

From equation (4.6.1), the propagator can be written as,

$$dI\_{\mathbf{x},i} \left(\Delta t\right) \equiv I\_{\mathbf{x},i} \left(t + \Delta t\right) - I\_{\mathbf{x},i} \left(t\right), \text{ for a given } \ I\_{\mathbf{x},i} \left(t\right). \tag{4.7.1}$$

At the same time, we can write equation (4.6.5) as,

$$dI\_{\mathbf{x},i}(\mathbf{t}) = N \Big[\mathbf{a}\_{ii}(\mathbf{t}) \Delta \mathbf{t}\Big]^{\mathbf{y}^{\prime}} + F\_{\mathbf{x},i} \Delta t\_{\prime} \tag{4.7.2}$$

where *N* is a unit normal random variable generated from *N*0,1 density function. In the limit, *t* 0,

*<sup>i</sup> t N A t tD t t <sup>i</sup>* , , *<sup>i</sup>* (4.6.5)

(4.6.7)

(4.6.8)

, , *At F i xi* and (4.6.9)

, , , , *x i x i x i dI t I t t I t* for a given *I t x i*, . (4.7.1)

, , , *x i ii x i dI t N a t t F t* (4.7.2)

2

. *Dt a i ii* (4.6.10)

, , , . *EI t I t I t a t x i xi xi ii* (4.6.6)

*<sup>i</sup>t* is a Gaussian

152 in Porous Media - An Approach Based on Stochastic Calculus

when *Ai t* and *D t <sup>i</sup>* are independent of *t* . As the variance can not be negative,

where " " indicates the random variable is generated from the normal density function

<sup>2</sup>

 <sup>2</sup> , . *E tI t a t i xi ii*

*t* .

. *Var t a <sup>i</sup> ii*

<sup>1</sup>

where *N* is a unit normal random variable generated from *N*0,1 density function. In

By comparing equation (4.5.7) with equation (4.6.3), and equation (4.6.4) with equation

We can now write the Langevin equation for the Ito diffusion, *I t x i*, .

From equation (4.6.1), the propagator can be written as,

At the same time, we can write equation (4.6.5) as,

*D t <sup>i</sup>* should be a non-negative function. As mentioned before,

variable, i.e., it is normally distributed; therefore, we can write,

having the mean, *At t <sup>i</sup>* and the variance *Dt t <sup>i</sup>* .

when , *t t* i.e., 0, *t* therefore 0 *<sup>i</sup>*

i.e., therefore

If we recall, equation (4.5.9) gives,

This can be written as,

*<sup>i</sup>* / 0 at , *t t* .

(4.6.8), we deduce that,

the limit, *t* 0,

**4.7 Langevin Equation for** *I t x i*,

*E t* 

This leads to

$$dI\_{\mathbf{x},i}\left(\mathbf{t}\right) = \mathbf{F}\_{\mathbf{x},i}dt + \left[a\_{\boldsymbol{\mu}}\left(\mathbf{t}\right)dt\right]^{\mathbf{X}}\mathbf{N}.\tag{4.7.3}$$

This is the first form of the Langevin equation for *I t x i*, .

It can be shown that, if *N Normal* 0,1 ,

$$
\beta \text{N} + \mathcal{A} = \text{Normal}\{\mathcal{A}, \mathcal{B}^2\},
\tag{4.7.4}
$$

Therefore, *dtN* is a normally distributed random variable having the density function *Normal dt* 0, . This is the density function of the standard Wiener process increments, *dw t* . Therefore, we can rewrite equation (4.7.3) as,

$$dI\_{x,i}\left(t\right) = F\_{x,i}dt + \sqrt{a\_{\text{ul}}}dw\_i\left(t\right). \tag{4.7.5}$$

Equation (4.7.5) is the second form of the Langevin equation.

The advantage of this Langevin approximation for *I t x i*, over equation (4.5.2) is that it is not a multidimensional SDE but a one-dimensional one. This would allow us to compute *I t x i*, s more efficiently.

The moment evolution equations for *I t x i*, can be given as follows: (See Gillespie (1992) for derivations.)

$$\frac{d\left[E\left(I\_{x,i}\left(t\right)\right)\right]}{dt} = E\left[F\_{x,i}\right]\_{\prime} \tag{4.7.6}$$

$$\frac{d\left[\operatorname{Var}\left(I\_{x,i}\left(t\right)\right)\right]}{dt} = 2\left(E\left[I\_{x,i}\left(t\right)F\_{x,i}\right] - E\left[I\_{x,i}\left(t\right)\right]E\left[F\_{x,i}\right]\right) + E\left[a\_{ii}\right],\tag{4.7.7}$$

with initial conditions,

$$E\left[I\_{\mathbf{x},i}\left(\mathbf{0}\right)\right] = \mathbf{0}, \quad \text{and} \tag{4.7.8}$$

$$\operatorname{Var}\left[\boldsymbol{I}\_{\boldsymbol{x},i}\left(\mathbf{0}\right)\right] = \mathbf{0}.\tag{4.7.9}$$

If *Fx i*, and *ii a* are deterministic functions of *x* then the moment evolution equations can be simplified to equations (4.5.7) and (4.5.8).

### **4.8 The Evaluation of** *C t <sup>x</sup>* **Diffusions**

The SSTM can be written as, for given *x* , (see equation (4.5.13)),

$$d\mathbf{C}\_{\mathbf{x}}\left(\mathbf{t}\right) = -\mathbf{C}\_{\mathbf{x}}\left(\mathbf{t}\right)d\mathbf{I}\_{\mathbf{x},0} - \left(\frac{\partial \mathbf{C}\_{\mathbf{x}}}{\partial \mathbf{x}}\right)\_{\mathbf{x}}d\mathbf{I}\_{\mathbf{x},1} - \left(\frac{\partial^2 \mathbf{C}\_{\mathbf{x}}}{\partial \mathbf{x}^2}\right)\_{\mathbf{x}}d\mathbf{I}\_{\mathbf{x},2} \tag{4.8.1}$$

Following Klebaner (1998), the expectations of infinitesimal differences of *C t <sup>x</sup>* can be

A Generalized Mathematical Model in One-Dimension 155

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

 <sup>2</sup> <sup>222</sup> <sup>012</sup> . *EC t C t C t <sup>x</sup> <sup>x</sup> <sup>x</sup>*

Following the same arguments as in the case of deriving the Langevin equation for *I t x i*, ,

 <sup>222</sup> <sup>012</sup> , *<sup>x</sup> dC t dt*

> 222 <sup>012</sup> ,

 . *<sup>x</sup> <sup>x</sup> dC t dt dw t* 

This shows that the concentration at given *x* can be characterised by a Langevin type stochastic differential equation. This equation can be used to develop numerical solutions of

The time evolution of the probability density function of *C t <sup>x</sup>* , *P C tC t t xx x* , ,, <sup>0</sup> is described by Fokker-Plank equation (Klebaner, 1998; Gillespie, 1992). The Fokker-Plank

2 2

*P yt t y t y y* , , 0 00

Once we solve the Fokker-Plank equation (4.8.12) along with its initial condition, the time evolution of probability density function can be found. This is also a weak solution to equation (4.8.10). Equation (4.8.10) also has strong solutions which can be obtained by

, , <sup>1</sup> , <sup>2</sup> *P yty t <sup>x</sup> Px x x P t y y*

0 0

integrating the SDE (4.8.11) using Ito integration. The drift coefficient

where *y* denotes *C t <sup>x</sup>* and *Px* stands for *P yty t <sup>x</sup>* , , 0 0 .

 

where *dw t* is independent increments of the standard Wiener process, and if

 *<sup>x</sup>*  then

2

(4.8.12)

<sup>0</sup> , (4.8.13)

in equation

, (4.8.8)

*dw t* (4.8.10)

(4.8.11)

(4.8.9)

written as (see also equation (4.5.7) and (4.5.9)),

we can obtain the Langevin equation for *C t <sup>x</sup>* :

the concentration profiles.

equation for 0 0 , , *P yty t <sup>x</sup>* is,

where

Equation (4.8.12) has the initial condition,

is the Dirac-delta function.

(4.8.11) is a stochastic variable in *x* and *t* .

and

where subscript *x* refers to the first and second derivatives with respect to *x* . Note that the coefficients are functions of *t* for given *x* , and *x i*, *I* s are Ito diffusions based on the velocity structure. Equation (4.8.1) is a stochastic diffusion and a SDE which displays the interplay between the concentration profile and velocity structures in the medium. By substituting the equations of the form of equation (4.7.5) for *x i*, *dI* s we obtain,

$$\begin{split}d\mathbf{C}\_{\times}(t) &= -\mathbf{C}\_{\times}(t) \Big(F\_{\times,0}dt + \sqrt{a\_{\rm on}}dw\_{\rm o}(t)\Big) \\ &- \left(\frac{\partial \mathbf{C}}{\partial \mathbf{x}}\right)\_{\times} \Big(F\_{\times,1}dt + \sqrt{a\_{\rm in}}dw\_{\rm o}(t)\Big) \\ &- \left(\frac{\partial^{2} \mathbf{C}}{\partial \mathbf{x}^{2}}\right)\_{\times} \Big(F\_{\times,2}dt + \sqrt{a\_{\rm on}}dw\_{\rm o}(t)\Big). \end{split} \tag{4.8.2}$$

In equation (4.8.2), 0 1 *dw dw*, and 2 *dw* are independent standard Wiener increments. Equation (4.8.2) can be written as,

$$d\mathcal{C}\_x\left(t\right) = -\alpha \left(\mathcal{C}\_x\left(t\right), t\right) dt - \sum\_{k=0}^{2} \beta\_k \left(\mathcal{C}\_x\left(t\right), t\right) dw\_{k\prime} \tag{4.8.3}$$

where,

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

$$
\beta\_0 = \mathbb{C}\_{\times}(t) \sqrt{a\_{00}} \text{ .}\tag{4.8.5a}
$$

$$
\beta\_1 = \left(\frac{\partial \mathbb{C}}{\partial \mathbf{x}}\right)\_\mathbf{x} \sqrt{a\_{11}} \text{ and } \tag{4.8.5b}
$$

$$
\mathcal{J}\_2 = \left(\frac{\partial^2 \mathbf{C}}{\partial \mathbf{x}^2}\right)\_\mathbf{x} \sqrt{a\_{22}} \,. \tag{4.8.5c}
$$

Now the equation (4.8.3) can be written as

$$\begin{split} \, \mathrm{d}C\_{x} \begin{pmatrix} t \\ \end{pmatrix} &= -\alpha \mathrm{d}t - \sum\_{k=0}^{2} \beta\_{k} \mathrm{d}w\_{k} \, \mathrm{d}w\_{k} \, \mathrm{d}w \\\\ &= -\alpha \mathrm{d}t + \begin{bmatrix} -\beta\_{0} & -\beta\_{1} & -\beta\_{2} \end{bmatrix} \begin{bmatrix} \mathrm{d}w\_{0} \\ \mathrm{d}w\_{1} \\ \mathrm{d}w\_{2} \end{bmatrix} . \end{split} \tag{4.8.6}$$

Equation (4.8.6) gives a diffusion matrix,

$$\begin{aligned} f &= \begin{bmatrix} -\beta\_0 & -\beta\_1 & -\beta\_2 \end{bmatrix} \begin{bmatrix} -\beta\_0 & -\beta\_1 & -\beta\_2 \end{bmatrix}^\top, \\\ &= \beta\_0^2 + \beta\_1^2 + \beta\_2^2 \end{aligned} \tag{4.8.7}$$

Following Klebaner (1998), the expectations of infinitesimal differences of *C t <sup>x</sup>* can be written as (see also equation (4.5.7) and (4.5.9)),

$$E\left(\mathbf{C}\_x(t+\Delta) - \mathbf{C}\_x(t)\middle|\mathbf{C}\_x(t)\right) = -\mathbf{a}\Delta + \mathbf{O}\left(\Delta\right) \tag{4.8.8}$$

and

Computational Modelling of Multi-Scale Non-Fickian Dispersion

(4.8.2)

, (4.8.4)

(4.8.6)

(4.8.7)

154 in Porous Media - An Approach Based on Stochastic Calculus

where subscript *x* refers to the first and second derivatives with respect to *x* . Note that the coefficients are functions of *t* for given *x* , and *x i*, *I* s are Ito diffusions based on the velocity structure. Equation (4.8.1) is a stochastic diffusion and a SDE which displays the interplay between the concentration profile and velocity structures in the medium. By

 

*dC t C t F dt a dw t*

 

*x x*

<sup>2</sup>

*dC t C t t dt C t t dw*

 <sup>2</sup> ,0 ,1 <sup>2</sup> ,2 , *<sup>x</sup> x x <sup>x</sup> <sup>x</sup> C C C t t C tF F F*

1 11 *x*

> 2 2 2 22 *x*

*x*

*x k k k*

*dC t dt dw*

222 012

 

*C*

2 0

0 1 20 1 2

 

 

*x*

*a*

*a*

*dt dw*

0121

 

,

*C*

In equation (4.8.2), 0 1 *dw dw*, and 2 *dw* are independent standard Wiener increments.

*x*

,1 11 1

*<sup>C</sup> F dt a dw t*

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

0 , , , *<sup>x</sup> <sup>x</sup> k x <sup>k</sup> k*

*x x*

(4.8.3)

<sup>0</sup> *Ct a <sup>x</sup>* <sup>00</sup> , (4.8.5a)

, and (4.8.5b)

. (4.8.5c)

0

*dw*

*dw*

.

2

, *T*

 

*x x*

*<sup>C</sup> F dt a dw t*

,0 00 0

substituting the equations of the form of equation (4.7.5) for *x i*, *dI* s we obtain,

2

*x*

Equation (4.8.2) can be written as,

Now the equation (4.8.3) can be written as

Equation (4.8.6) gives a diffusion matrix,

*f*

where,

*x x x*

$$E\left(\left(\mathbf{C}\_x\left(t+\Delta\right)-\mathbf{C}\_x\left(t\right)\right)^2\Big|\mathcal{C}\_x\left(t\right)\right) = \left(\beta\_0^2 + \beta\_1^2 + \beta\_2^2\right)\Delta + \mathcal{O}\left(\Delta\right).\tag{4.8.9}$$

Following the same arguments as in the case of deriving the Langevin equation for *I t x i*, , we can obtain the Langevin equation for *C t <sup>x</sup>* :

$$d\mathcal{C}\_x\left(\mathbf{t}\right) = -\alpha dt + \sqrt{\beta\_0^2 + \beta\_1^2 + \beta\_2^2} \, dw\left(\mathbf{t}\right),\tag{4.8.10}$$

where *dw t* is independent increments of the standard Wiener process, and if

$$
\mathcal{J}\_x \equiv \sqrt{\mathcal{J}\_0^2 + \mathcal{J}\_1^2 + \mathcal{J}\_2^2} \quad \text{then}
$$

$$
d\mathcal{C}\_x(t) = -\alpha dt + \mathcal{J}\_x dw(t). \tag{4.8.11}
$$

This shows that the concentration at given *x* can be characterised by a Langevin type stochastic differential equation. This equation can be used to develop numerical solutions of the concentration profiles.

The time evolution of the probability density function of *C t <sup>x</sup>* , *P C tC t t xx x* , ,, <sup>0</sup> is described by Fokker-Plank equation (Klebaner, 1998; Gillespie, 1992). The Fokker-Plank equation for 0 0 , , *P yty t <sup>x</sup>* is,

$$\frac{\partial P\_{\text{x}}\left(y\_{\prime}t\middle|y\_{0\prime}t\_{0}\right)}{\partial t} = \frac{\partial \left(\alpha P\_{\text{x}}\right)}{\partial y} + \frac{1}{2} \frac{\partial^{2} \left(\beta\_{\text{x}}^{2}P\_{\text{x}}\right)}{\partial y^{2}},\tag{4.8.12}$$

where *y* denotes *C t <sup>x</sup>* and *Px* stands for *P yty t <sup>x</sup>* , , 0 0 .

Equation (4.8.12) has the initial condition,

$$P\{y\_\prime t = t\_0 | y\_{0\prime} t\_0\} = \mathcal{S}(y - y\_0)\_\prime \tag{4.8.13}$$

where is the Dirac-delta function.

Once we solve the Fokker-Plank equation (4.8.12) along with its initial condition, the time evolution of probability density function can be found. This is also a weak solution to equation (4.8.10). Equation (4.8.10) also has strong solutions which can be obtained by integrating the SDE (4.8.11) using Ito integration. The drift coefficient in equation (4.8.11) is a stochastic variable in *x* and *t* . be

verify the numerical solutions for the Fokker-Planck equation.

concertration *C t <sup>x</sup>* for a given *x* in the vicinity of x:

*x*

,1

is sufficiently small; the drift coefficient

*P Cx x* 0 1.0.

**4.9 Numerical Solutions** 

In equation (4.9.2),

continuous function.

*<sup>x</sup>* in equation (4.9.1) is given by,

and,

where,

0 1

A Generalized Mathematical Model in One-Dimension 157

when *E y* and *Var y* are given by equations (4.8.15) and (4.8.18). When time is zero,

Equation (4.8.19) should also be the solution to the Fokker-Plank equation (4.8.12). The Fokker-Plank equation (4.8.12) is a stochastic partial differential equation for which analytical equations can not easily be obtained. Therefore, we can make use of this fact to

We have seen in the previous section 4.8, the following SDE gives the time course of

 

where *dw t* is the standard Wiener increments with a zero mean and *dt* variance, if *dt*

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

,0 2

*x*

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

 

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

 , *<sup>x</sup> <sup>x</sup> dC t dt dw t* 

<sup>1</sup> , ,0 , exp . <sup>2</sup> <sup>2</sup> *<sup>x</sup> <sup>x</sup> y Ey P yty P yt Var y Var y*

2

 

is given by (see equation (4.8.4)),

(4.9.3)

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

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

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

(4.9.4)

*x x*

*x*

*x x*

<sup>2</sup>

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

where *V xt* , in the mean velocity which is assumed to be regular differentiable

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

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

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

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

2

(4.8.19)

(4.9.1)

(4.9.2)

Similar to equations (4.7.6) and (4.7.7), we can write time evolution of the moments *C t <sup>x</sup>* . The derivation of these equations are very similar to those given by Gillespie (1992).

The time evolution of the mean of *C t <sup>x</sup>* is given by,

$$\frac{d\left(E\left(\mathbb{C}\_{\times}(t)\right)\right)}{dt} = E\left[-\alpha\right],\tag{4.8.14}$$

therefore,

$$E\left(\mathbf{C}\_{x}\left(t\right)\right) = \bigcup\_{0}^{t} E\left[-\alpha\right]dt.\tag{4.8.15}$$

The mean of *C t <sup>x</sup>* at a given *x* is expressed as an integral of the expectation of . As can be seen from equation (4.8.4), is not only dependent on *C t <sup>x</sup>* and its first and second derivatives with respect to *x* , but also dependent on the mean velocity and its first and second derivatives with respect to *x* , according to equations (4.2.20), (4.2.21) and (4.2.22). The initial condition for equation (4.8.15) is *Cx* 0 , the value of the concentration at time is zero for a given *x* .

The evolution of the variance of *C t Var C t <sup>x</sup>* , *<sup>x</sup>* is given by,

$$\frac{d\left[\operatorname{Var}\left(\mathbb{C}\_{\boldsymbol{x}}\left(t\right)\right)\right]}{dt} = 2\left(E\left(-a\mathbb{C}\_{\boldsymbol{x}}\left(t\right)\right) + E\left(\mathbb{C}\_{\boldsymbol{x}}\left(t\right)\right)E\left(a\right)\right) + E\left[\mathcal{J}\_{\boldsymbol{x}}^{2}\right],\tag{4.8.16}$$

and the variance of *C t <sup>x</sup>* can be obtained by integrating equation (4.8.16) with respect to *t* ,

$$Var\left(\mathbb{C}\_{\boldsymbol{x}}\left(t\right)\right) = 2\oint\_{0} \left(-\alpha \mathbb{C}\_{\boldsymbol{x}}\left(t\right)\right)dt + 2\oint\_{0} \left(\mathbb{C}\_{\boldsymbol{x}}\left(t\right)\right)E\left(\alpha\right)dt + \int\_{0}^{t} \mathbb{E}\left[\left.\mathcal{J}\_{\boldsymbol{x}}^{2}\right]dt.\tag{4.8.17}$$

By Fubini's theorem (Klebaner,1998), for continuous stochastic variable such as , *C t <sup>x</sup>* and *<sup>x</sup>* , we can rewrite equation (4.8.17) as,

$$\operatorname{Var}\left(\mathbb{C}\_{\boldsymbol{x}}\left(\boldsymbol{t}\right)\right) = -2\operatorname{E}\left[\int\_{0}^{t} \alpha \,\mathbb{C}\_{\boldsymbol{x}}\left(\boldsymbol{t}\right)d\boldsymbol{t}\right] + 2\operatorname{E}\left[\operatorname{C}\_{\boldsymbol{x}}\left(\boldsymbol{t}\right)\right]\operatorname{E}\left(\alpha\right)d\boldsymbol{t} + \operatorname{E}\left[\int\_{0}^{t} \beta\_{\boldsymbol{x}}^{2}\right]d\boldsymbol{t}.\tag{4.8.18}$$

In equation (4.8.18), Fubini's theorem is only applied to the first and third term of equation (4.8.17).

Once we evaluate the mean and variance of *C t <sup>x</sup>* , we can obtain the probability density function, 0 0 , , *P yty t <sup>x</sup>* of *C t <sup>x</sup>* (y represents the value of *C t <sup>x</sup>* ). As Ito diffusions are Martingales with Markovian properties (Klebaner, 1998), and *C t <sup>x</sup>* is Gaussian,

$$P\_x\left(y, t \middle| y\_0, 0\right) = P\_x\left(y, t\right) = \frac{1}{\left[2\operatorname{Var}\left(y\right)\right]^{\frac{1}{2}}} \exp\left(\frac{-\left(y - E\left(y\right)\right)^2}{2\operatorname{Var}\left(y\right)}\right). \tag{4.8.19}$$

when *E y* and *Var y* are given by equations (4.8.15) and (4.8.18). When time is zero, *P Cx x* 0 1.0.

Equation (4.8.19) should also be the solution to the Fokker-Plank equation (4.8.12). The Fokker-Plank equation (4.8.12) is a stochastic partial differential equation for which analytical equations can not easily be obtained. Therefore, we can make use of this fact to verify the numerical solutions for the Fokker-Planck equation.

### **4.9 Numerical Solutions**

Computational Modelling of Multi-Scale Non-Fickian Dispersion

, (4.8.14)

. As

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

(4.8.15)

is not only dependent on *C t <sup>x</sup>* and its first and

 

> 

 *x*

156 in Porous Media - An Approach Based on Stochastic Calculus

Similar to equations (4.7.6) and (4.7.7), we can write time evolution of the moments *C t <sup>x</sup>* .

 *<sup>x</sup> dEC t E dt*

 0

The mean of *C t <sup>x</sup>* at a given *x* is expressed as an integral of the expectation of

The evolution of the variance of *C t Var C t <sup>x</sup>* , *<sup>x</sup>* is given by,

*t E C t E dt <sup>x</sup>*

second derivatives with respect to *x* , but also dependent on the mean velocity and its first and second derivatives with respect to *x* , according to equations (4.2.20), (4.2.21) and (4.2.22). The initial condition for equation (4.8.15) is *Cx* 0 , the value of the concentration at

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

and the variance of *C t <sup>x</sup>* can be obtained by integrating equation (4.8.16) with respect to *t* ,

 <sup>2</sup> 0 0 0

<sup>2</sup>

In equation (4.8.18), Fubini's theorem is only applied to the first and third term of equation

Once we evaluate the mean and variance of *C t <sup>x</sup>* , we can obtain the probability density function, 0 0 , , *P yty t <sup>x</sup>* of *C t <sup>x</sup>* (y represents the value of *C t <sup>x</sup>* ). As Ito diffusions are

*Var C t E C t dt E C t E dt E dt <sup>x</sup>*

*x x*

*x x*

By Fubini's theorem (Klebaner,1998), for continuous stochastic variable such as

Martingales with Markovian properties (Klebaner, 1998), and *C t <sup>x</sup>* is Gaussian,

*E C t EC t E E*

2 2 . *t t t Var C t E C t dt E C t E dt E dt <sup>x</sup>*

> 0 0 0 2 2 . *t t t*

.

*x x x*

, (4.8.16)

(4.8.18)

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

The derivation of these equations are very similar to those given by Gillespie (1992).

The time evolution of the mean of *C t <sup>x</sup>* is given by,

can be seen from equation (4.8.4),

*d Var C t*

*dt*

*<sup>x</sup>* , we can rewrite equation (4.8.17) as,

time is zero for a given *x* .

therefore,

and 

(4.8.17).

We have seen in the previous section 4.8, the following SDE gives the time course of concertration *C t <sup>x</sup>* for a given *x* in the vicinity of x:

$$d\mathbb{C}\_{\times}(t) = -adt + \mathcal{J}\_{\times} dw(t),\tag{4.9.1}$$

where *dw t* is the standard Wiener increments with a zero mean and *dt* variance, if *dt* is sufficiently small; the drift coefficient is given by (see equation (4.8.4)),

$$
\delta \mathbf{x} = \mathbf{C}\_x(t) 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{4.9.2}
$$

In equation (4.9.2),

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

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

and,

$$F\_{x,3} = \frac{h\_x}{2} \overline{V} \left( x, t \right);\tag{4.9.5}$$

where *V xt* , in the mean velocity which is assumed to be regular differentiable continuous function.

*<sup>x</sup>* in equation (4.9.1) is given by,

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

where,

$$
\beta\_0 = \mathbb{C}\_{\times}(t)\sqrt{a\_{00}}\,\tag{4.9.7}
$$

$$
\mathcal{O}\_1 = \left(\frac{\partial \mathcal{C}}{\partial \mathbf{x}}\right)\_\times \sqrt{a\_{11'}} \quad \text{and}$$

$$
\mathcal{O}\_2 = \left(\frac{\partial^2 \mathcal{C}}{\partial \mathbf{x}^2}\right)\_\times \sqrt{a\_{22'}} \tag{4.9.8}$$

By using a numerical scheme developed, we obtain realizations of *C t <sup>x</sup>* as strong solutions to equation (4.9.1). We can also obtain solutions to the Fokker-Plank equation

A Generalized Mathematical Model in One-Dimension 159

In this Chapter, we develop a generalized form of SSTM that can include any arbitrary velocity covariance kernel in principle. We have demonstrates that for a given kernel, a generalized analytical forms for eigen functions can be obtained by using the computational methods developed. We have also developed a Langevin form of the SSTM for a given *x* , and the time evolution of concentration, *C t <sup>x</sup>* , follows a stochastic differential equation

functions. In other words, if one monitors the concentration *C t <sup>x</sup>* at a given point in space, the data collected along with time would constitute a realization of the strong solution of the SDE. The solution is a function of the covariance kernel, i.e. a function of <sup>2</sup>

and *b* for an exponentially decaying kernel, and also a function of *C t <sup>x</sup>* itself and its first- and second derivatives with respect to *x* . This focus of SDE provides a very convenient and computationally efficient way to solve the stochastic partial differential

By deriving a Langevin form of the SSTM, we essentially prove that any time course of the concentration at a given point behaves according to the underlying SDE, which would characterize the nature of local porous medium and is a statement of mass concentration of

*<sup>x</sup>* which are again functions of *C t <sup>x</sup>* and eigen

(4.8.12) using the same finite differences.

and behaves

form along

**4.10 Remarks for the Chapter** 

having the coefficients

the solute.

equation associated with the SSTM.

The following equation gives the expressions of 00 11 *a a*, , and 22 *a* ,

 2 2 1 , ,0,1,2 , *m ii ij j aG i* (4.9.9)

where *m* is the number of effective eigen functions, and

$$G\_{\vec{\eta}} = \sqrt{\mathcal{k}\_{\vec{\eta}}} P\_{\vec{\eta}}, \quad \text{ (i, 0, 1, 2)}. \tag{4.9.10}$$

In the numerical solutions, we make use of the finite differences, for a given dependent variable, say *U* , based on the grid given in Figure 4.20

Figure 4.20. Space-time grid used in the numerical solutions with respect to y

$$\left(\frac{\partial \mathcal{L}I}{\partial y}\right)\_n^m = \frac{\left(\mathcal{U}\_n^m - \mathcal{U}\_{n-1}^m\right)}{\Delta y},\tag{4.9.11}$$

$$\left(\frac{\partial^2 \mathcal{U}}{\partial y^2}\right)\_u^m = \frac{\left(\mathcal{U}\_n^m - 2\mathcal{U}\_{n-1}^m + \mathcal{U}\_{n-2}^m\right)}{\Delta y^2},\tag{4.9.12}$$

and

$$\left(\frac{\partial \mathcal{U}}{\partial t}\right)\_{\boldsymbol{u}}^{\boldsymbol{m}} = \frac{\left(\boldsymbol{U}\_{\boldsymbol{u}}^{\boldsymbol{m}+1} - \boldsymbol{U}\_{\boldsymbol{u}}^{\boldsymbol{m}}\right)}{\Delta t}.\tag{4.9.13}$$

By using a numerical scheme developed, we obtain realizations of *C t <sup>x</sup>* as strong solutions to equation (4.9.1). We can also obtain solutions to the Fokker-Plank equation (4.8.12) using the same finite differences.

### **4.10 Remarks for the Chapter**

Computational Modelling of Multi-Scale Non-Fickian Dispersion

(4.9.8)

(4.9.10)

(4.9.9)

158 in Porous Media - An Approach Based on Stochastic Calculus

*a*

and

*a*

, ,0,1,2 ,

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

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

*x*

*C*

2 2

, ( ,0,1,2). *G Pi ij j ij* 

In the numerical solutions, we make use of the finite differences, for a given dependent

*C*

 

 

The following equation gives the expressions of 00 11 *a a*, , and 22 *a* ,

where *m* is the number of effective eigen functions, and

variable, say *U* , based on the grid given in Figure 4.20

1

Figure 4.20. Space-time grid used in the numerical solutions with respect to y

*n*

and

*n*

*n*

*U U U y y* 

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

*U UUU y y* 

*m mmm nnn*

 <sup>1</sup> , *m m m n n*

<sup>1</sup>

*m m m n n*

*U U U t t*

<sup>2</sup> ,

.

(4.9.11)

(4.9.12)

(4.9.13)

*m ii ij j aG i*

*x*

In this Chapter, we develop a generalized form of SSTM that can include any arbitrary velocity covariance kernel in principle. We have demonstrates that for a given kernel, a generalized analytical forms for eigen functions can be obtained by using the computational methods developed. We have also developed a Langevin form of the SSTM for a given *x* , and the time evolution of concentration, *C t <sup>x</sup>* , follows a stochastic differential equation having the coefficients and *<sup>x</sup>* which are again functions of *C t <sup>x</sup>* and eigen functions. In other words, if one monitors the concentration *C t <sup>x</sup>* at a given point in space, the data collected along with time would constitute a realization of the strong solution of the SDE. The solution is a function of the covariance kernel, i.e. a function of <sup>2</sup> and *b* for an exponentially decaying kernel, and also a function of *C t <sup>x</sup>* itself and its first- and second derivatives with respect to *x* . This focus of SDE provides a very convenient and computationally efficient way to solve the stochastic partial differential equation associated with the SSTM. collected

By deriving a Langevin form of the SSTM, we essentially prove that any time course of the concentration at a given point behaves according to the underlying SDE, which would characterize the nature of local porous medium and is a statement of mass concentration of the solute.

**Theories of Fluctuations and Dissipation** 

In the previous chapters, we see that the hydrodynamic dispersion is in fact a result of solute particles moving along a decreasing pressure gradient and encountering the solid surfaces of a porous medium. The pressure gradient provides the driving force which translates into kinetic energy, and the porous medium acts as the dissipater of the kinetic energy; any such energy dissipation associated with small molecules generates fluctuations among molecules. Looking at a molecular-level picture, the dissolved solute particles in water travelling through the porous medium slow down nearing a surface and then increase in velocity once the molecules get scattered after the impact with solid surface. Refining this picture a bit more, we see that the velocity boundary layers along the solid surfaces are helping this process. Not all the molecules hit solid surfaces either; some of these would be subjected to micro-level local pressure gradients and move away from the surfaces. A physical ensemble of these solute molecules would depict behaviours that are measurable using appropriate extensive variables. (Extensive variables depend on the extent of the system of molecules. i.e., the number of molecules, concentrations, kinetic energy etc., where as intensive variables do not change with size of the system, i.e., pressure, temperature, entropy etc.) These measurable quantities at macroscopic level have origins in microscopic level. Therefore, we can anticipate that molecular level description would justify the operational models that we develop at an ensemble level. Naturally one could expect that the statistical moments of the variables of an ensemble would lead to meaningful models of

In the development of the SSTM, we express the velocity of solute as the sum of the mean velocity and a fluctuating component around the mean. The mean velocity may then be evaluated by using the Darcy's law. We then express the fluctuating component in terms of the spectral expansion dependent on a covariance kernel. However, we need to understand that this type of picture in a more fundamental way should be based on the established theories. Towards that end, in this chapter, we review some of the fundamental theoretical frameworks associated with molecular fluctuation. We show the connectivity of thermodynamical, molecular and stochastic description of fluctuations and dissipations, and then we make use of Ito diffusions to obtain the models of statistical moments of relevant variables. While we do not cite the reference within this chapter -- as the works we refer to are well accepted knowledge in the disciplines such as thermodynamics, statistical mechanics and stochastic processes-- all the relevant works are given in the references list at the end of the book. However, we refer to Keizer's work (1987) primarily in this chapter.

works

**5.1 Introduction** 

the process we would like to observe.
