## Abstract

This article presents a new approach to solve the equations of flow in heterogeneous porous media by using random walks on regular lattices. The hydraulic head is represented by computational particles which are spread globally from the lattice sites according to random walk rules, with jump probabilities determined by the hydraulic conductivity. The latter is modeled as a realization of a random function generated as a superposition of periodic random modes. One- and two-dimensional numerical solutions are validated by comparisons with analytical manufactured solutions. Further, an ensemble of divergence-free velocity fields computed with the new approach is used to conduct Monte Carlo simulations of diffusion in random fields. The transport equation is solved by a global random walk algorithm which moves computational particles representing the concentration of the solute on the same lattice as that used to solve the flow equations. The integrated flow and transport solution is validated by a good agreement between the statistical estimations of the first two spatial moments of the solute plume and the predictions of the stochastic theory of transport in groundwater.

## Authors

**Nicolae Suciu, **Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

## Keywords

global random walk; flow in heterogeneous porous media; Monte Carlo simulations.

## References

*Flow and Transport in Porous Formations*. Springer, Berlin (1989) https://doi.org/10.1007/978-3-642-75015-1

*Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study*. Adv. Water Resour.

**34**, 47–61 (2011) http://dx.doi.org/10.1016/j.advwatres.2010.09.012CrossRefGoogle Scholar

*Code verification by the method of manufactured solutions*. J. Fluids Eng.

**124**(1):4–10 (2002) http://dx.doi.org/10.1115/1.1436090CrossRefGoogle Scholar

*Diffusion in random velocity fields with applications to contaminant transport in groundwater*. Adv. Water Resour.

**69**(2014) 114-133. http://dx.doi.org/10.1016/j.advwatres.2014.04.002CrossRefGoogle Scholar

*Diffusion in Random Fields. Applications to Transport in Groundwater.*Birkhäuser, Cham (2019) https://doi.org/10.1007/978-3-030-15081-5

*Evaluation of overshooting errors in particle methods for diffusion by biased global random walk.*Rev. Anal. Num. Th. Approx. (Romanian Academy)

**35**, 119–126 (2006) https://ictp.acad.ro/jnaat/journal/article/view/2006-vol35-no1-art17

*Generalized random walk algorithm for the numerical modeling of complex diffusion processes.*J. Comput. Phys.,

**186**, 527–544 (2003) https://doi.org/10.1016/S0021-9991(03)00073-1.MathSciNetCrossRefGoogle Scholar

##### Cite this paper as:

N. Suciu, *Global Random Walk Solutions for Flow and Transport in Porous Media*, in: Numerical Mathematics and Advanced Applications ENUMATH 2019, pp. 939-947, doi: https://doi.org/10.1007/978-3-030-55874-1_93

## About this paper

##### Journal

Springer

##### Publisher Name

##### Print ISSN

Not available yet.

##### Online ISSN

##### Google Scholar Profile

^{1}

^{1}institutetext: Nicolae Suciu, Department of Mathematics, Friedrich-Alexander University Erlangen-Nuremberg,

Cauerstraße. 11, 91058 Erlangen, Germany,

^{1}

^{1}email: suciu@math.fau.de

Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy,

Fantanele 57, 400320 Cluj-Napoca, Romania

# Global Random Walk Solutions for Flow and

Transport in Porous Media

###### Abstract

This article presents a new approach to solve the equations of flow in heterogeneous porous media by using random walks on regular lattices. The hydraulic head is represented by computational particles which are spread globally from the lattice sites according to random walk rules, with jump probabilities determined by the hydraulic conductivity. The latter is modeled as a realization of a random function generated as a superposition of periodic random modes. One- and two-dimensional numerical solutions are validated by comparisons with analytical manufactured solutions. Further, an ensemble of divergence-free velocity fields computed with the new approach is used to conduct Monte Carlo simulations of diffusion in random fields. The transport equation is solved by a global random walk algorithm which moves computational particles representing the concentration of the solute on the same lattice as that used to solve the flow equations. The integrated flow and transport solution are validated by a good agreement between the statistical estimations of the first two spatial moments of the solute plume and the predictions of the stochastic theory of transport in groundwater.

## 1 Introduction

Global random walk (GRW) algorithms, unconditionally stable and free of numerical diffusion, solve parabolic partial differential equations by moving computational particles on regular lattices according to random walk rules Vamosetal2003 . The GRW approach, intensively used in Monte Carlo simulations of diffusion in random velocity fields Suciu2014 , is based on the relationship between the ensemble of trajectories of the diffusion process governed by the Itô equation and the probability density of the process which verifies a Fokker-Planck equation. The strength of the approach is that instead of sequentially generating random walk trajectories one counts the number of particles at lattice sites, updated according to the binomial distributions which describe the number of jumps to neighboring lattice sites. Since the total number of particles is arbitrarily large, one obtains in this way highly accurate numerical approximations of the probability density, which, in case of transport simulations, is the normalized concentration of the solute that is transported by advection and diffusion.

The second order differential operator of the Fokker-Planck equation coincides with the diffusion operator in the mass balance equation with closure given by the Fick’s law only if the diffusion coefficient is constant. Otherwise, in order to rewrite the diffusion equation as a Fokker-Planck equation, the drift coefficients have to be augmented with the spatial derivatives of the diffusion coefficients (e.g., (Suciu2019, , Sect. 4.2.1)). The same procedure can be used to write the pressure equation for flow in porous media based on Darcy’s law (e.g., Raduetal2011 ) as a Fokker Planck equation. Further, numerical solutions of the flow problem can be obtained by solving the problem formulated for the equivalent Fokker-planck equation by GRW methods.

In case of highly variable coefficients, the particles may jump over several lattice sites, the variability of the coefficients is smoothed, and the GRW solutions are affected by overshooting errors SuciuandVamos2006 . Such errors are avoided by using biased-GRW algorithms, where the particles are only allowed to jump to neighboring sites and the advective displacements are accounted for by biased jump probabilities that are larger in the direction of the drift (Suciu2019, , Sect. 3.3.3). A different remedy, appropriate for GRW solutions of the flow problem, consists of approximating directly the second order operator in diffusion form, with flux terms estimated at the middle of the lattice intervals by using staggered grids (Suciu2019, , Sect. 3.3.4.1).

We consider a two-dimensional domain, $(x,y)\in \mathrm{\Omega}=[0,{L}_{x}]\times [0,{L}_{y}]$ and the pressure equation for flow in saturated porous media with constant porosity,

$$\frac{\partial h}{\partial t}-\left[\frac{\partial}{\partial x}\left(K\frac{\partial h}{\partial x}\right)+\frac{\partial}{\partial y}\left(K\frac{\partial h}{\partial y}\right)\right]=f,$$ | (1) |

where $h(t,x,y)$ is the hydraulic head, $K(x,y)$ is an isotropic hydraulic conductivity, and $f(x,y)$ is a source/sink term. For an arbitrary initial condition (IC) and time independent boundary conditions (BC), the non-steady solution of Eq. (1) approaches a stationary hydraulic head $h(x,y)$. With $f=0$ and BC given by

$h(0,y)$ | $=H,h({L}_{x},y)=0,\forall y\in [0,{L}_{y}]$ | (2) | ||

$\frac{\partial h}{\partial y}}(x,0)$ | $={\displaystyle \frac{\partial h}{\partial y}}(x,{L}_{y})=0,\forall x\in [0,{L}_{x}],$ | (3) |

the gradient of the hydraulic head $h(x,y)$ determines the components of the divergence-free velocity field according to Darcy’s law

$${V}_{x}=-K\frac{\partial h}{\partial x},{V}_{y}=-K\frac{\partial h}{\partial y}.$$ | (4) |

The GRW solutions presented in the following are tested against analytical manufactured solutions Roache2002 obtained from data files and codes given in the Git repository https://github.com/PMFlow/FlowBenchmark. The isotropic hydraulic conductivity is a space random function with mean $\u27e8K\u27e9=15$ m/day. The random fields $\mathrm{ln}K$ are normally distributed, with Gaussian correlations of fixed correlation length $\lambda =1$ m and increasing variances ${\sigma}^{2}=0.1,\mathrm{\hspace{0.33em}1},\mathrm{\hspace{0.33em}2}$, with realizations generated as sums of a finite number $N=100$ of cosine random modes (Suciu2019, , Appendix C.3.1.2).

## 2 One-dimensional GRW algorithms

The three different GRW approaches discussed in Section 1 are tested in the following by solving the one-dimensional (1D) version of Eq. (1) written as an advection-diffusion equation,

$$\frac{\partial h}{\partial t}+\frac{\partial}{\partial x}\left(\frac{\partial K}{\partial x}h\right)=\frac{{\partial}^{2}}{\partial {x}^{2}}(Kh)+f,$$ | (5) |

in the interval $(0,L)$, $L=200\lambda $. The source $f$ and the BCs $h(0)$ and $h(L)$ are determined by the manufactured solution $\stackrel{~}{h}(x)=3+\mathrm{sin}(x)$ (see Git repository, /FlowBenchmark/Manufactured_Solutions/Matlab/Gauss 1D/). The precision of the numerical solution $h(x)$ is quantified by the discrete ${L}^{2}$ norm $\epsilon ={\Vert h-\stackrel{~}{h}\Vert}_{{L}^{2}}$.

GRW algorithms approximate the hydraulic head by the distribution $n(i,k)$ at lattice sites $i$ and times $k$ of a system of $\mathcal{N}$ random walkers, $h(i\mathrm{\Delta}x,k\mathrm{\Delta}t)\approx n(i,k)/\mathcal{N}$. The unbiased GRW moves groups of particles on the lattice according to

$$n(j,k)=\delta n(j+{v}_{j},j,k)+\delta n(j+{v}_{j}-d,j,k)+\delta n(j+{v}_{j}+d,j,k),$$ | (6) |

$$n(i,k+1)=\delta n(i,i,k)+\sum _{j\ne i}\delta n(i,j,k)+\lfloor \mathcal{N}f\mathrm{\Delta}t\rfloor ,$$ | (7) |

where $\lfloor \cdot \rfloor $ is the floor function, ${v}_{j}=\lfloor {V}_{j}\mathrm{\Delta}t/\mathrm{\Delta}x+0.5\rfloor $, ${V}_{j}=\frac{\partial K}{\partial x}(j\mathrm{\Delta}x)$, are discrete advective displacements, and $d$ is the amplitude of the diffusive jumps. The time step $\mathrm{\Delta}t$ and the space step $\mathrm{\Delta}x$ are related by

$$K(j\mathrm{\Delta}x)={r}_{j}\frac{{(d\mathrm{\Delta}x)}^{2}}{2\mathrm{\Delta}t},$$ | (8) |

where ${r}_{j}$, $0\le {r}_{j}\le 1$, is a variable jump probability. The number of particles undergoing diffusion jumps, $\delta n(j+{v}_{j}\mp d,j,k)$, and the number of particles waiting at $j+{v}_{j}$ over the $k$ time step, $\delta n(j+{v}_{j},j,k)$, are binomial random variables with mean values given by (Suciu2019, , Sect. 3.2.2)

$$\overline{\delta n(j,j+{v}_{j},k)}=(1-{r}_{j})\overline{n(j,k)},\overline{\delta n(j+{v}_{j}\mp d,j,k)}=\frac{1}{2}{r}_{j}\overline{n(j,k)}.$$ |

In the biased GRW algorithm $d=1$ and Eq. (6) is replaced by

$$n(j,k)=\delta n(j,j,k)+\delta n(j-1,j,k)+\delta n(j+1,j,k).$$ |

The quantities $\delta n$ verify in the mean

$$\overline{\delta n(j,j,k)}=(1-{r}_{j})\overline{n(j,k)},\overline{\delta n(j\mp 1,j,k)}=\frac{1}{2}({r}_{j}\mp {v}_{j})\overline{n(j,k)},$$ |

where ${v}_{j}={V}_{j}\delta t/\delta x$. The biased jump probabilities are given by $({r}_{j}\mp {v}_{j})$ with ${r}_{j}$ defined by (8) and, in addition to ${r}_{j}\le 1$, one imposes the constraint $\mid {v}_{j}\mid \le {r}_{j}$ (Suciu2019, , Sect. 3.3.3).

Alternatively, a GRW solution of the 1D version of the pressure equation (1),

$$\frac{\partial h}{\partial t}-\frac{\partial}{\partial x}\left(\frac{K\partial h}{\partial x}\right)=f,$$ |

can be obtained as steady-state limit of the staggered finite difference scheme

$\frac{1}{\mathrm{\Delta}t}}[h(i,k+1)-h(i,k)]=$ | ||

$\frac{1}{\mathrm{\Delta}{x}^{2}}}\{[K(i+1/2)(h(i+1,k)-h(i,k))]-[K(i-1/2)(h(i,k)-h(i-1,k))]\}-f.$ |

With jump probabilities defined by $r(i\mp 1/2)=K(i\mp 1/2)\mathrm{\Delta}t/\mathrm{\Delta}{x}^{2}$, $r\le 1/2$, the staggered scheme becomes

$n(i,k+1)$ | $=$ | $\left[1-r(i-1/2)-r(i+1/2)\right]n(i,k)$ | (9) | ||

$+r(i-1/2)n(i-1,k)+r(i+1/2)n(i+1,k)-\lfloor \mathcal{N}f\mathrm{\Delta}t\rfloor .$ |

The contributions to lattice sites $i$ from neighboring sites summed up in (9) are obtained with the GRW algorithm which moves particles from sites $j$ to neighboring sites $i=j\mp 1$ according to

$n(j,k)=\delta n(j,j,k)+\delta n(j-1,j,k)+\delta n(j+1,j,k).$ | (10) |

For consistency with the staggered scheme (9), the quantities $\delta n$ in (11) have to satisfy in the mean (Suciu2019, , Sect. 3.3.4.1),

$$\overline{\delta n(j,j,k)}=\left[1-r(j-1/2)-r(j+1/2)\right]\overline{n(j,k)},\overline{\delta n(j\mp 1,j,k)}=r(j\mp 1/2)\overline{n(j,k)}.$$ | (11) |

The three GRW algorithms from above approximate the binomial random variables $\delta n$ by rounding of the unaveraged relations for the mean, e.g., (11), summing up the reminders of multiplication by $r$ and of the floor function $\lfloor \mathcal{N}f\mathrm{\Delta}t\rfloor $, and allocating a particle to the lattice site where the sum reaches the unity. By giving up the particle indivisibility, one obtains deterministic GRW algorithms which represent the solution $n$ by real numbers and use the anaveraged relations for $\delta n$.

Since the unbiased GRW algorithm is prone to large overshooting errors, due to the high variability of the coefficients $\frac{\partial K}{\partial x}$ and $K$, it is not appropriate to solve the flow problem. For instance, to solve the test problem for ${\sigma}^{2}=0.1$ with errors $\epsilon \sim {10}^{-1}$, an extremely small step $\mathrm{\Delta}x={10}^{-5}$ is needed and one time iteration takes about 8 minutes. Efficient solutions can be obtained with the biased GRW algorithm (Method 1) and the GRW on staggered grids in both the random (Method 2) and the deterministic implementation (Method 3). The stationary regime of the GRW simulations, indicated by constant total number of particles and ${L}^{2}$ error (see Fig 1), is reached after $T=2\cdot {10}^{7}$ time iterations. The comparison shown in Table 1 indicates that the deterministic implementation of the GRW on staggered grids is the most efficient approach to solve flow problems for heterogeneous porous media.

Method 1 | Method 2 | Method 3 | |

$\epsilon $ | 4.31e-02 | 1.60e-2 | 1.60e-02 |

CPU (min) | 37.12 | 14.03 | 6.48 |

## 3 Two-dimensional GRW solutions of the flow problem

The 2D GRW algorithm on staggered grids is defined similarly to (9-11) as follows:

$n(i,j,k+1)$ | $=$ | $\left[1-r(i-1/2,j)-r(i+1/2,j)-r(i,j-1/2)-r(i,j+1/2)\right]n(i,j,k)$ | (12) | ||

$+r(i-1/2,j)n(i-1,j,k)+r(i+1/2,j)n(i+1,j,k)$ | |||||

$+r(i,j-1/2)n(i,j-1,k)+r(i,j+1/2)n(i,j+1,k)-\lfloor \mathcal{N}f\mathrm{\Delta}t\rfloor ,$ |

$n(l,m,k)$ | $=$ | $\delta n(l,m|l,m,k)$ | (13) | ||

$+\delta n(l-1,m|l,m,k)+\delta n(l+1,m|l,m,k)$ | |||||

$+\delta n(l,m-1|l,m,k)+\delta n(l,m+1|l,m,k),$ |

$\overline{\delta n(l,m,k)}=\left[1-r(l-1/2,m)-r(l+1/2,m)-r(l,m-1/2)-r(l,m+1/2)\right]\overline{n(l,m,k)}$ | (14) | ||||

$\overline{\delta n(l\mp 1,m|l,m,k)}=r(l\mp 1/2,m)\overline{n(l,m,k)}$ | |||||

$\overline{\delta n(l,m\mp 1|l,m,k)}=r(l,m\mp 1/2)\overline{n(l,m,k)}.$ |

The problem (1-3), with ${L}_{x}=20\lambda $, ${L}_{y}=10\lambda $, $H=1$, $\mathrm{\Delta}x=\mathrm{\Delta}y={10}^{-1}$, and IC given by the manufactured solution, ${h}_{0}(x,y)=\stackrel{~}{h}(x,y)=1+\mathrm{sin}(2x+y)$, is solved by the deterministic GRW on staggered grids. The ${L}^{2}$ errors with respect to the manufactured solution are shown in Table 2.

${\sigma}^{2}=0.1$ | ${\sigma}^{2}=1$ | ${\sigma}^{2}=2$ | |
---|---|---|---|

$\epsilon $ | 3.12e-02 | 4.08e-02 | 1.35e-01 |

$T$ | 2e06 | 1e07 | 2e07 |

CPU (min) | 11.26 | 55.49 | 111.54 |

The convergence of the 2D GRW solutions is investigated for the homogeneous problem ($f=0$, and IC given by the constant slope ${h}_{0}(0,y)=1,{h}_{0}(L,y)=0$), for increasing ${\sigma}^{2}$, by successively halving the step size five times from $\mathrm{\Delta}x=\mathrm{\Delta}y={10}^{-1}$ to $\mathrm{\Delta}x=\mathrm{\Delta}y=3.125\cdot {10}^{-3}$. Using the errors with respect to the solution on the finest grid, ${\epsilon}_{k}={\Vert {h}^{(k)}-{h}^{(6)}\Vert}_{{L}^{2}}$, the estimated order of convergence (EOC) is computed according to $EOC=\mathrm{log}\left({\epsilon}_{k}/{\epsilon}_{k+1}\right)/\mathrm{log}(2),k=1,\mathrm{\dots},4$. The results presented in Table 3 indicate the convergence of order 2 for all the three values of ${\sigma}^{2}$.

${\sigma}^{2}$ | ${\epsilon}_{1}$ | EOC | ${\epsilon}_{2}$ | EOC | ${\epsilon}_{3}$ | EOC | ${\epsilon}_{4}$ | EOC | ${\epsilon}_{5}$ |
---|---|---|---|---|---|---|---|---|---|

$0.1$ | 2.78e-03 | 1.9979 | 6.96e-04 | 2.0167 | 1.72e-04 | 2.0687 | 4.10e-05 | 2.3219 | 8.20e-06 |

$1$ | 2.63e-03 | 1.9902 | 6.62e-04 | 2.0131 | 1.64e-04 | 2.0685 | 3.91e-05 | 2.3219 | 7.82e-06 |

$2$ | 2.67e-03 | 1.9839 | 6.75e-04 | 2.0150 | 1.67e-04 | 2.0654 | 3.99e-05 | 2.3219 | 7.98e-06 |

## 4 Statistical inferences

A Monte Carlo ensemble of solutions is obtained by solving the 2D homogeneous problem ($f=0$) for 100 realizations of the random $K$-field with fixed ${\sigma}^{2}=0.1$. A realization of the hydraulic head is shown in Fig 2. The components of the velocity field computed according to Darcy’s law (4) are shown in Fig. 3.

The Monte Carlo estimates presented in Table 4 are close to the first-order theoretical predictions of the stochastic theories of flow and transport in groundwater Dagan1989 , e.g. variances of the order ${\sigma}_{h}^{2}\sim {(\sigma \lambda H/{L}_{x})}^{2}\sim {10}^{-4}$ for the hydraulic head, ${\sigma}_{{V}_{x}}^{2}=\frac{3}{8}{\sigma}^{2}=3.75\cdot {10}^{-2}$, and ${\sigma}_{{V}_{y}}^{2}=\frac{1}{8}{\sigma}^{2}=1.25\cdot {10}^{-2}$, for longitudinal and transverse velocity components, respectively.

$h$ | ${V}_{x}$ | $Vy$ | |
---|---|---|---|

mean | 5.00e-01$\pm $ 1.74e-01 | 9.98e-01$\pm $ 1.91e-02 | -7.32e-04$\pm $ 1.24e-02 |

variance | 3.27e-04$\pm $ 6.27e-05 | 3.68e-02$\pm $ 5.50e-03 | 1.30e-02$\pm $ 1.90e-03 |

## 5 Validation of the flow and transport GRW solutions

The numerical setup from Section 4 corresponds to a possible scenario of contaminant transport in groundwater systems with low heterogeneity Suciu2014 . With a typical value of the local dispersion coefficient ${D}_{0}=0.01$ m^{2}/day, Monte Carlo simulations of advection-diffusion are carried out using the ensemble of 100 velocity fields, on the same lattice as that used in Section 4 to solve the flow problem. The transport problem is solved with the 2D generalization of the unbiased GRW algorithm (6-7) (Suciu2019, , Appendix A.3.2). The Monte Carlo estimates of the velocity of the plume’s center of mass and of the effective dispersion coefficients presented in Fig. 4 are in good agreement with the first-order results obtained from simulations using a linear approximation of the velocity field (Suciu2019, , Appendix C.3.2.2).

## References

- (1) Dagan, G.: Flow and Transport in Porous Formations. Springer, Berlin (1989) https://doi.org/10.1007/978-3-642-75015-1
- (2) Radu, F.A., Suciu, N., Hoffmann, J., Vogel, A., Kolditz, O., Park, C-H., Attinger, S.: Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study. Adv. Water Resour. 34, 47–61 (2011) http://dx.doi.org/10.1016/j.advwatres.2010.09.012
- (3) Roache, P.J.: Code verification by the method of manufactured solutions. J. Fluids Eng. 124 (1):4–10 (2002) http://dx.doi.org/10.1115/1.1436090
- (4) Suciu, N.: Diffusion in random velocity fields with applications to contaminant transport in groundwater. Adv. Water Resour. 69 (2014) 114-133. http://dx.doi.org/10.1016/j.advwatres.2014.04.002
- (5) Suciu, N.: Diffusion in Random Fields. Applications to Transport in Groundwater. Birkhäuser, Basel (2019) https://doi.org/10.1007/978-3-030-15081-5
- (6) Suciu, N., Vamoş, C.: Evaluation of overshooting errors in particle methods for diffusion by biased global random walk. Rev. Anal. Num. Th. Approx. (Romanian Academy) 35, 119–126 (2006) https://ictp.acad.ro/jnaat/journal/article/view/2006-vol35-no1-art17
- (7) Vamoş, C., Suciu, N., Vereecken, H.: Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J. Comput. Phys., 186, 527-544 (2003) https://doi.org/10.1016/S0021-9991(03)00073-1.