Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media

4 years ago

(original), Convergence order, global random walk, ISI/JCR, Numerical Analysis, Numerical Modeling, paper

Abstract

In this article, we present new random walk methods to solve flow and transport problems in saturated/unsaturated porous media, including coupled flow and transport processes in soils, heterogeneous systems modeled through random hydraulic conductivity and recharge fields, processes at the field and regional scales. The numerical schemes are based on global random walk algorithms (GRW) which approximate the solution by moving large numbers of computational particles on regular lattices according to specific random walk rules. To cope with the nonlinearity and the degeneracy of the Richards equation and of the coupled system, we implemented the GRW algorithms by employing linearization techniques similar to the L $L$ -scheme developed in finite element/volume approaches. The resulting GRW L-schemes converge with the number of iterations and provide numerical solutions that are first-order accurate in time and second-order in space. A remarkable property of the flow and transport GRW solutions is that they are practically free of numerical diffusion. The GRW solvers are validated by comparisons with mixed finite element and finite volume solvers in one- and two-dimensional benchmark problems. They include Richards’ equation fully coupled with the advection-diffusion-reaction equation and capture the transition from unsaturated to saturated flow regimes.

Authors

Nicolae Suciu
Mathematics Department, Friedrich-Alexander University of Erlangen-Nurnberg, Erlangen, Germany
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj-Napoca, Romania

Davide Illiano
Department of Mathematics, University of Bergen, Norway

Alexander Prechtel
Department of Mathematics, University of Bergen, Norway

Florin A. Radu
Department of Mathematics, University of Bergen, Norway

Keywords

Richards equation; Coupled flow and transport; Linearization; Iterative schemes; Global random walk

References

see the expanding block below

PDF

Extended version published in Arxiv:2011.12889v3

Cite this paper as:

N. Suciu, D. Illiano, A. Prechtel, F.A.Radu, Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media, Advances in Water Resources, 152 (2021), 103935, https://doi.org/10.1016/j.advwatres.2021.103935

About this paper

Journal

Elsevier

Publisher Name

Advances in Water Resources

DOI

https://doi.org/10.1016/j.advwatres.2021.103935

Print ISSN

Not available yet.

Online ISSN

ISSN: 0309-1708

Not available yet.

Google Scholar Profile

Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media

Nicolae Suciu*,** suciu@math.fau.de , Davide Illiano***, Alexander Prechtel*, Florin A. Radu***
*Mathematics Department, Friedrich-Alexander University of Erlangen-Nürnberg, Cauerstraße. 11, 91058 Erlangen, Germany
**Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Fântanele 57, 400320 Cluj-Napoca, Romania
***Department of Mathematics, University of Bergen, Allégaten 41, 5007 Bergen, Norway

Abstract

In this article, we present new random walk methods to solve flow and transport problems in saturated/unsaturated porous media, including coupled flow and transport processes in soils, heterogeneous systems modeled through random hydraulic conductivity and recharge fields, processes at the field and regional scales. The numerical schemes are based on global random walk algorithms (GRW) which approximate the solution by moving large numbers of computational particles on regular lattices according to specific random walk rules. To cope with the nonlinearity and the degeneracy of the Richards equation and of the coupled system, we implemented the GRW algorithms by employing linearization techniques similar to the $L$ -scheme developed in finite element/volume approaches. The resulting GRW $L$ -schemes converge with the number of iterations and provide numerical solutions that are first-order accurate in time and second-order in space. A remarkable property of the flow and transport GRW solutions is that they are practically free of numerical diffusion. The GRW solvers are validated by comparisons with mixed finite element and finite volume solvers in one- and two-dimensional benchmark problems. They include Richards’ equation fully coupled with the advection-diffusion-reaction equation and capture the transition from unsaturated to saturated flow regimes.

Keywords:

Richards equation , Coupled flow and transport , Linearization , Iterative schemes , Global random walk MSC: 76S05 , 65N12 , 86A05 , 65C35 , 76R50

^†^†journal: Advances in Water Resources

1 Introduction

The accuracy and the robustness of the numerical schemes is the primary requirement for reliable and meaningful results of the current efforts to improve the understanding of the complexity and interdependence of the flow and transport processes in subsurface hydrology through numerical investigations. Numerical solvers for partial differential equations modeling individual or coupled processes are often used as basic elements in the formulation of the more complex problems of practical interest, such as parameter identification [17], hydraulic tomography [7], Monte Carlo approaches for systems with randomly distributed parameters [25], or upscaling for mutiphase flows in heterogeneous subsurface formations [15], among others.

A central issue in subsurface hydrology is the need of robust and computationally efficient numerical models for partially saturated soil-groundwater systems. The transition between unsaturated and saturated zones is particularly challenging. In unsaturated flows the water content $\theta$ and the hydraulic conductivity $K$ depend nonlinearly on the pressure head $\psi$ through material laws based on experiments, as far as $\psi<0$ . The evolution of $\psi$ is governed by the parabolic Richards’ equation which degenerates to a (generally) linear elliptic equation (i.e. the equation for steady-state flow in aquifers) if $\psi\geq 0$ [4]. Since the regions where degeneracy takes place depend on the evolution of the pressure $\psi$ in time and space, they are not known a priori. To cope with the nonlinearity and degeneracy of the Richards’ equation, different linearization methods are needed, such as the Newton scheme [34, 15, 19], which is second-order convergent but converges only locally (requires a starting point close enough to the solution) or the more robust but only first-order convergent Picard, modified Picard, or $L$ schemes [37, 28, 24, 31].

Explicit and implicit schemes proposed for nonlinear flows in unsaturated regime provide solutions with comparable accuracy but are generally ambiguous to compare in terms of computing time. Since they do not need to solve systems of linear algebraic equations at every time step, explicit schemes are in principle faster [23] but their speed may be seriously affected by the need to use very small time steps [16, 11, 2]. The time step in explicit schemes is constraint by stability conditions [39, 23, 11] and has to be significantly reduced to ensure small local Péclet number (Pé), defined with respect to the space step. Large (global) Pé characterizes advection-dominated transport problems [5, 21]. In such cases, reducing the local Pé is a remedy to avoid the numerical diffusion and the oscillatory behavior of the solution [29]. The criterion of small local Pé is also recommended for numerical schemes solving the pressure equation in saturated flows [14] and, since Richards’ equation has the structure of the advection-diffusion equation, the recommendation holds for the unsaturated flows as well.

Well known approaches to avoid the numerical diffusion are the particle tracking in continuous space and the discrete random walk on lattices [40]. The accuracy of these schemes is determined by the number of computational particles undergoing random jumps in continuous space or on discrete lattices. In random walk schemes, the increase of the computation time with the number of particles is simply avoided by randomly distributing the particles along the spatial directions with a global procedure, according to appropriate jump probabilities. In this way, one obtains a global random walk (GRW) which performs the spreading of all the particles from a given site with computational costs that are practically the same as for generating the jump of a single random walker in sequential procedures [44]. In particular cases (e.g., when using biased jump probabilities to account for variable coefficients or for advective displacements) the GRW algorithms are equivalent to explicit finite difference schemes with time step size constrained by stability requirements. In unbiased GRW schemes for transport problems with variable coefficients, which still satisfy stability conditions, no restrictions on the time step are needed to reduce the local Pé number, which renders the approach particularly efficient in large scale simulations of transport in groundwater (see [40] for details and examples).

The elliptic and parabolic equations governing the pressure head for flows in unsaturated/saturated porous media are essentially diffusion equations with second order operator in Stratonovich form. They can be recast as Fokker-Planck equations, with drift augmented by the row derivative of the coefficient tensor, and further solved by random walk approaches [40]. An alternative approach starts with a staggered finite difference scheme, further used to derive biased random walk rules governing the movement on a regular lattice of a system of computational particles. The particle density at lattice sites provides a numerical approximation of the pressure head solution. This approach has been already illustrated for flows in saturated porous media with heterogeneous hydraulic conductivity [2, 41].

In this article, we present new GRW schemes for nonlinear and non-steady flows in soils which model the transition from unsaturated to saturated regime in a way consistent with the continuity of the constitutive relationships $\theta(\psi)$ and $K(\psi)$ . Following [24, 31], the nonlinearity of the Richards equation is solved with an iterative procedure similar to the $L$ -scheme used in finite element/volume approaches. Numerical tests demonstrate the convergence of the $L$ -scheme for unsaturated/saturated flows. For fully saturated flow regime with constant water content $\theta$ and time independent boundary conditions the GRW $L$ -scheme is equivalent to a transient finite difference scheme.

Coupled flow and reactive transport problems for partially saturated soils rise new stability and consistency issues and demand augmented computational resources. Our GRW approach in this case consists of coupling the flow solver described above with existing GRW transport solvers [40] adapted for nonlinear problems, which are implemented as $L$ -schemes as well. The flow and transport solvers are coupled via an alternating splitting procedure [18] which successively iterates the corresponding $L$ -schemes until the convergence of the pressure head and concentration solutions is reached, within the same tolerance, at every time step. Code verification tests using analytical manufactured solutions are employed to verify the convergence of the iterations and the accuracy of the splitting scheme.

The GRW scheme for one-dimensional solutions of the Richards equation, which captures the transition from unsaturated to saturated flow regimes is validated by comparisons with solutions provided by Richy software, based on the mixed finite element method (MFEM), with backward Euler discretization in time and Newton linearization, developed at the Mathematics Department of the Friedrich-Alexander University of Erlangen-Nürnberg [33, 34]. For the particular case of unsaturated flows, the one-dimensional GRW solutions are also tested by comparisons with experimental data [48, 49] and exact solutions published in the literature [46, 47]. The two-dimensional GRW solutions are compared on benchmark problems with two-point flux approximation (TPFA) finite volume solvers using backward Euler discretization in time and $L$ -scheme linearization [18]. The TPFA codes are implemented in MRST, the MATLAB Reservoir Simulation Toolbox [22].

The paper is organized as follows. Section 2 presents the GRW algorithm and the linearization approach for one-dimensional flow problems. The one-dimensional solver is further validated through comparisons with MFEM solutions, experimental data, and exact analytical solutions in Section 3. Two dimensional GRW algorithms for fully coupled and decoupled flow and transport problems are introduced in Section 4. Code verification tests and comparisons with TPFA solutions for benchmark problems are presented in Section 5. Some examples of flow and transport solutions for groundwater problems at the field and regional scale are presented in Section 6. The main conclusions of this work are finally presented in Section 7. GRW codes implemented in Matlab for model problems considered in this article are stored in the Git repository RichardsEquation [43].

2 One-dimensional GRW algorithm for unsaturated/saturated flow in soils

We consider the water flow in unsaturated/saturated porous media described by the one-dimensional Richards equation [16, 33, 19] in the space-time domain $[0,L_{z}]\times[0,T]$ ,

\frac{\partial\theta(\psi)}{\partial t}-\frac{\partial}{\partial z}\left[K(% \theta(\psi))\frac{\partial}{\partial z}(\psi+z)\right]=0,

(1)

where $\psi(z,t)$ is the pressure head expressed in length units, $\theta$ is the volumetric water content, $K$ stands for the hydraulic conductivity of the medium, and $z$ is the height oriented positively upward. According to (1), the water flux given by Darcy’s law is $q=-K(\theta(\psi))\frac{\partial}{\partial z}(\psi+z)$ .

To design a GRW algorithm, we start with the staggered finite difference scheme with backward discretization in time which approximates the solution of Eq. (1) at positions $z=i\Delta z$ , $i=1,\ldots,I$ , $I=L_{z}/\Delta z$ , and time points $t=k\Delta t$ , $k=1,\ldots,T/\Delta t$ , according to

\theta(\psi_{i,k})-\theta(\psi_{i,k-1})=\frac{\Delta t}{{\Delta z}^{2}}\left\{% \left[K(\psi_{i+1/2,k})(\psi_{i+1,k}-\psi_{i,k})-K(\psi_{i-1/2,k})(\psi_{i,k}-% \psi_{i-1,k})\right]+\left(K(\psi_{i+1/2,k})-K(\psi_{i-1/2,k})\right)\Delta z% \right\}.

(2)

To cope with the double nonlinearity due to the dependencies $K(\theta)$ and $\theta(\psi)$ we propose an explicit scheme similar to the linearization approach known as “ $L$ -scheme”, originally developed for implicit methods [e.g. 28, 24, 31]. The approach consists of the addition of a stabilization term $L(\psi_{i,k}^{s+1}-\psi_{i,k}^{s})$ , $L=const$ , in the left-hand side of (2) and of performing successive iterations $s=1,2,\ldots$ of the modified scheme until the discrete $L^{2}$ norm of the solution $\psi_{k}^{s}=(\psi_{i,k}^{s},\ldots,\psi_{I,k}^{s})$ verifies

\|\psi_{k}^{s}-\psi_{k}^{s-1}\|\leq\varepsilon_{a}+\varepsilon_{r}\|\psi_{k}^{% s}\|

(3)

for some given tolerances $\varepsilon_{a}$ and $\varepsilon_{r}$ . The adapted $L$ -scheme reads

\psi_{i,k}^{s+1}=\left[1-(r_{i+1/2,k}^{s}+r_{i-1/2,k}^{s})\right]\psi_{i,k}^{s% }+r_{i+1/2,k}^{s}\psi_{i+1,k}^{s}+r_{i-1/2,k}^{s}\psi_{i-1,k}^{s}+\left(r_{i+1% /2,k}^{s}-r_{i-1/2,k}^{s}\right)\Delta z-\left(\theta(\psi_{i,k}^{s})-\theta(% \psi_{i,k-1})\right)/L,

(4)

where

r_{i\pm 1/2,k}^{s}=K(\psi_{i\pm 1/2,k}^{s})\Delta t/(L{\Delta z}^{2}).

(5)

For fixed time step $k$ , the iterations start with the solution after the last iteration at the previous time $k-1$ , $\psi_{i,k}^{1}=\psi_{i,k-1}$ , $i=1,\ldots,I$ . Note that, unlike implicit $L$ -schemes (e.g., [37, 28, 24]), the explicit scheme (4) uses forward increments of $\psi$ . In this way, the solution $\psi_{i,k}^{s+1}$ is obtained from values of $\psi$ and $r$ evaluated at the previous iteration, without solving systems of algebraic equations.

The solution $\psi_{i,,k}^{s}$ is further represented by the distribution of $\mathcal{N}$ computational particles at the sites of the one-dimensional lattice, $\psi_{i,k}^{s}\approx n_{i,k}^{s}a/\mathcal{N}$ , with $a$ being a constant equal to a unit length, and the $L$ -scheme (4) becomes

n_{i,k}^{s+1}=\left[1-\left(r_{i+1/2,k}^{s}+r_{i-1/2,k}^{s}\right)\right]n_{i,% k}^{s}+r_{i+1/2,k}^{s}n_{i+1,k}^{s}+r_{i-1/2,k}^{s}n_{i-1,k}^{s}+\left\lfloor% \mathcal{N}f^{s}\right\rfloor,

(6)

where the source term is defined as $f^{s}=\left(r_{i+1/2,k}^{s}-r_{i-1/2,k}^{s}\right)\Delta z-\left[\theta(n_{i,k% }^{s})-\theta(n_{i,k-1}\right]/L$ and $\left\lfloor\cdot\right\rfloor$ denotes the floor function.

The physical dimension of the parameter $L$ of the scheme is that of an inverse length unit to ensure that $r_{i\pm 1/2,k}^{s}$ defined by (5) are non-dimensional parameters, as needed in random walk approaches [44, 40]. By imposing the constraint $r_{i\pm 1/2,k}^{s}\leq 1/2$ , the parameters $r_{i\pm 1/2,k}^{s}$ can be thought of as biased jump probabilities. Hence, the contributions to $n_{i,k}^{s+1}$ from neighboring sites $i\pm 1$ summed up in (6) can be obtained with the GRW algorithm which moves particles from sites $j$ to neighboring sites $i=j\mp 1$ according to the rule

\displaystyle n_{j,k}^{s}=\delta n_{j,j,k}^{s}+\delta n_{j-1,j,k}^{s}+\delta n% _{j+1,j,k}^{s}.

(7)

For consistency with (6), the quantities $\delta n^{s}$ in (7) have to satisfy in the mean [40, Sect. 3.3.4.1],

\overline{\delta n_{j,j,k}^{s}}=\left[1-\left(r_{j-1/2,k}^{s}+r_{j+1/2,k}^{s}% \right)\right]\overline{n_{j,k}^{s}},\;\;\overline{\delta n_{j\mp 1/2,j,k}^{s}% }=r_{j\mp 1/2,k}^{s}\overline{n_{j,k}^{s}}.

(8)

The quantities $\delta n^{s}$ are binomial random variables approximated by using the unaveraged relations (8) for the mean, summing up the reminders of multiplication by $r$ and of the floor function $\lfloor\mathcal{N}f^{s}\Delta t\rfloor$ , and allocating one particle to the lattice site where the sum reaches the unity.

Remark 1

The finite difference $L$ -scheme (4) and the derived GRW relation (6) are explicit iterative schemes for Richards equation in mixed form (1). The essential difference of the $L$ -schemes from explicit schemes in mixed formulation designed for unsaturated regime [16, 11, 23] is that they apply to both unsaturated and saturated flow conditions.

Remark 2

Consider the saturated regime, $\theta=const$ , with space-variable hydraulic conductivity $K$ and a given source term $f$ . With the parameter $L$ set to $L=1/a$ , after disregarding the time index $k$ the scheme (6) solves the following equation for the hydraulic head $h=\psi+z$ ,

\frac{1}{a}\frac{\partial h}{\partial s}-\frac{\partial}{\partial z}\left[K% \frac{\partial h}{\partial z}\right]=f.

(9)

For boundary conditions independent of $s$ , the solution of Eq. (9) approaches a steady-state regime corresponding to the saturated flow (see also [2, 41]). The modified GRW scheme (6) is equivalent to a convergent finite difference scheme first order accurate in time and second order in space [42, Remark 1].

3 Validation of the one-dimensional GRW flow algorithm

3.1 Transition from unsaturated to saturated flow regime

The one-dimensional algorithm for flow in unsaturated/saturated soils is validated in the following by comparisons with MFEM solutions obtained with the Richy software [33, 34]. For this purpose, we solve one-dimensional model-problems for the vertical infiltration of the water through both homogeneous and non-homogeneous soil columns [38], previously used in [34] to assess the accuracy and the convergence of the MFEM solutions.

We consider the domain $z\in[0,2]$ and the boundary conditions specified by a constant pressure $\psi(0,t)=\psi_{0}$ at the bottom of the soil column and a constant water flux $q_{0}$ at the top. Together, these constant conditions determine the initial pressure distribution $\psi(z,0)$ as solution of the steady-state flow problem. For $t>0$ , the pressure $\psi_{0}$ is kept constant, at the bottom, and the water flux at the top of the column is increased linearly from $q_{0}$ to $q_{1}$ until $t\leq t_{1}$ and is kept constant for $t>t_{1}$ .

Refer to caption — Figure 1: Time steps for Scenario (1) and Scenario (2).

For the unsaturated regions ( $\psi<0$ ) we consider the constitutive relationships given by the simple exponential model [13]

\theta(\psi)=\theta_{res}+(\theta_{sat}-\theta_{res})e^{\alpha\psi},

(10)

K(\theta(\psi))=K_{sat}\frac{\theta(\psi)-\theta_{res}}{\theta_{sat}-\theta_{% res}},

(11)

where $\theta=\theta_{sat}$ and $K=K_{sat}$ denote the constant water content respectively the constant hydraulic conductivity in the saturated regions ( $\psi\geq 0$ ) and $\theta_{res}$ is the residual water content. The more complex and physically sounded van Genuchten-Mualem parameterization model will be used for two-dimensional problems in the following sections.

The flow problem for Eq. (1) with the parameterization (10-11) is solved in two Scenarios: (1) homogeneous soil, with $K_{sat}=2.77\cdot 10^{-6}$ , $\theta_{res}=0.06$ , $\theta_{sat}=0.36$ , $\alpha=10$ , $q_{0}=2.77\cdot 10^{-7}$ , $q_{1}=2.50\cdot 10^{-6}$ , which are representative for a sandy soil, and (2) non-homogeneous soil, with the same parameters as in Scenario (1), except the saturated hydraulic conductivity, which takes two constant values, $K_{sat}=2.77\cdot 10^{-6}$ for $z<1$ and $500K_{sat}$ for $z\geq 1$ (modeling, for instance, a column filled with sand and gravel). To capture the transition from unsaturated to saturated regime, the pressure at the bottom boundary is fixed at $\psi_{0}=0.5$ . For the parameters of the one-dimensional flow problems solved in this section we consider meters as length units and seconds as time units. The simulations are conducted up to $T=10^{4}$ (about 2.78 hours) and the intermediate time is taken as $t_{1}=T/10^{2}$ .

We consider a uniform GRW lattice with $\Delta z=10^{-2}$ , equal to the length of the linear elements in the MFEM solver. The GRW computations are initialized by multiplying the initial condition by $\mathcal{N}=10^{24}$ particles. Since, as shown by (11), the hydraulic conductivity varies in time, the length of the time step determined by (5) for the maximum of $K$ at every time iteration and by specifying a maximum $r_{max}=0.8$ of the parameter $r_{i\pm 1/2,k}$ may vary in time (see Fig. 2). The parameter of the regularization term in the $L$ -scheme is set to $L=1$ for the computation of the initial condition (solution of the stationary problem, i.e. for $\partial\theta/\partial t=0$ in (1)) and to $L=2$ for the solution of the non-stationary problem. In both cases, the convergence criterion (3) is verified by choosing $\varepsilon_{a}=0$ and a relative tolerance $\varepsilon_{r}=10^{-9}$ .

The comparison with the MFEM solutions presented in Figs. 2-4 shows a quite good accuracy of the GRW solutions for pressure, water content, and water flux. The relative errors, computed with the aid of the $L^{2}$ norms by $\varepsilon_{\psi}=\|\psi^{GRW}-\psi^{MFEM}\|/\|\psi^{MFEM}\|$ , and similarly for $\theta$ and $q$ , are presented in Table 1.

Table 1: Error norms of the GRW solutions.

	$\varepsilon_{\psi}$	$\varepsilon_{\theta}$	$\varepsilon_{q}$
Scenario (1)	1.81e-02	2.20e-02	3.50e-02
Scenario (2)	5.20e-03	2.35e-02	2.07e-02

The $L$ -scheme converges with speeds depending on the problem. To solve the problem for the initial condition, one needs $3.5\cdot 10^{4}$ iterations in Scenario (1) and $6.5\cdot 10^{6}$ iterations in Scenario (2). Instead, to solve the non-stationary problem for a final time $T=10^{4}$ , one needs about 70 iterations in Scenario (1) and about 700 iterations in Scenario (2) (see Fig. 6 and Fig. 6). The convergence of the iterative GRW $L$ -scheme can be further investigated through assessments of the computational order of convergence of the sequence of successive correction norms $\|\psi_{k}^{s}-\psi_{k}^{s-1}\|$ [9, 10]. Estimations provided in [42, Appendix A] indicate a linear convergence for Scenario (1) but only a power law convergence $\sim s^{-1}$ , which is slower than the linear convergence [10], for Scenario (2).

Supplementary tests done in Scenario (1) indicate the existence of a lower bound of the constant $L$ which ensures the convergence [42, Sect. 3.1]. It is found that increasing $L$ above the value which ensures the convergence of the GRW $L$ -scheme with a desired accuracy only results in increasing number of iterations and more computing time. The parameter $L$ has to be established experimentally by checking the convergence and, as highlighted by the examples presented in Section 5 below, it depends on the complexity of the problem to be solved.

3.2 Comparison with experiments and exact solutions for unsaturated flows

An experiment consisting of free drainage in a 600 cm deep lysimeter filled with a material with silty sand texture conducted at the Los Alamos National Laboratory [1] is often used to validate one-dimensional schemes for unsaturated flows (see e.g., [48, 49, 11]). This example is provided with the Hydrus 1D software [36], which is also used for validation purposes in the papers cited above.

The relationships defining the water content $\theta(\psi)$ and the hydraulic conductivity $K(\theta(\psi))$ are given by the van Genuchten-Mualem model

\Theta(\psi)=\begin{cases}\left(1+(-\alpha\psi)^{n}\right)^{-m},&\psi<0\\ 1,&\psi\geq 0,\end{cases}

(12)

K(\Theta(\psi))=\begin{cases}K_{sat}\Theta(\psi)^{\frac{1}{2}}\left[1-\left(1-% \Theta(\psi)^{\frac{1}{m}}\right)^{m}\right]^{2},&\psi<0\\ K_{sat},&\psi\geq 0,\end{cases}

(13)

where $\theta_{res}$ , $\theta_{sat}$ , and $K_{sat}$ represent the same parameters as for the exponential model considered in Section 3.1, $\Theta=(\theta-\theta_{res})/(\theta_{sat}-\theta_{res})$ is the normalized water content, and $\alpha$ , $n$ and $m=1-1/n$ are model parameters depending on the soil type.

With the parameters given in the Hydrus 1D example, $\theta_{res}=0.0$ , $\theta_{sat}=0.331$ , $K_{sat}=25$ cm/d, $\alpha=0.0143$ cm^-1, $n=1.5$ , for initial and boundary conditions for free drainage given by $\psi(z,0)=0$ cm and $q(0,t)=0$ cm/d [48], the solutions provided by the GRW $L$ -scheme (6-8) for simulation times from 1 d to 100 d are obtained with stabilization parameter $L=0.5$ after a number of 9 to 35 iterations (tolerance specified by $\varepsilon_{a}=\varepsilon_{r}=5\cdot 10^{-6}$ in (3)). The spatial resolution is set to $\Delta z=10$ cm, while the time step varies slightly between $10^{-2}$ d and $3.16\cdot 10^{-2}$ d, according to (5). The results are compared in Figs. 8 and 8 with Hydrus 1D results and experimental data. The water content profiles (Fig. 8) are quite close to measurements and similar to those presented in [49, 11]. The pressure profiles (Fig. 8) deviate from experiment, mainly for $T=1$ d and $T=100$ d, with approximately the same amount as in [11, Fig. 12]. An improved prediction of the pressure profiles is obtained in [48] with slightly modified parameters of the van Genuchten-Mualem model, but with the price of larger deviations for the water content.

The $\theta$ -based form of Richards equation has shown significantly improved performance in numerical schemes for unsaturated flows in spatially homogeneous soils (e.g., constant $K_{sat}$ ), especially in modeling infiltration into dry media [48], and is well suited to analytical approaches [26, 46, 32]. Philip [26] derived an exact solution for infiltration problems expressed in the form $z(\theta,t)$ , that is, the depth where the water content takes specified values at given time points $t$ . Philip’s solution has been used in [46] to construct a table of coefficients which allow the computation of $z(\theta,t)$ for three different $\theta$ and arbitrary $t$ . The solution verifies the dimensionless form the $\theta$ -based Richards equation

\frac{\partial\Theta}{\partial t}=\frac{\partial}{\partial z}\left[D(\Theta)% \frac{\partial}{\partial z}(\Theta)\right]-\frac{dK(\Theta)}{d\Theta}\frac{% \partial\Theta}{\partial z},

(14)

where $z$ is positive downward, $D(\Theta)=K(\Theta)d\psi/d\Theta$ , and $K(\Theta)$ is given according to the van Genuchten-Mualem model by the upper branch of (13). Such analytical solutions have been used in [27, 11] to verify various one-dimensional numerical schemes based on finite volume and finite element approaches. In order to test the GRW $L$ -scheme (6-8), we solve the same infiltration problem (soil column 100 cm deep, constant unsaturated initial water content $\theta_{i}$ , and infiltration imposed by $\psi=0$ on the upper boundary). We use a van Genuchten-Mualem parameter $n=1.5$ together with the parameters of the hypothetical loam soil used in [46]: $K_{sat}=6\cdot 10^{-4}$ cm/s, $\theta_{sat}=0.45$ , $\theta_{res}=0.1$ , $\theta_{i}=0.17$ , $\alpha=0.01$ cm^-1. The pressure corresponding to the initial water content is obtained by (12), $\psi(\theta_{i})=-24.87$ cm. The computations are carried out with $\Delta z=1$ cm, $\Delta t$ between $9.26\cdot 10^{-4}$ h and $5.23\cdot 10^{-4}$ h, $L=0.2$ , and the convergence is achieved after a number of 15 to 160 iterations ( $\varepsilon_{a}=\varepsilon_{r}=5\cdot 10^{-6}$ ). The analytical solutions $\tilde{z}(\theta,t)$ for $\theta=$ 0.24, 0.31, and 0.38 at successive times between 0.5 h and 2 h are obtained with the coefficients for $n=1.5$ and $\Theta(\theta_{i})=0.2$ given in [46, Table 3]. The GRW results $z(\theta,t)$ for the same $\theta$ and $t$ are obtained by linear interpolation of the numerical solution $\theta(z,t)$ . Relative errors $(z-\tilde{z})/\tilde{z}$ of the numerical solution $z(\theta,t)$ with respect to the analytical solution $\tilde{z}(\theta,t)$ are shown in Table 2.

Table 2: Relative errors of the GRW solution

z(\theta,t)

with respect to the exact solution from [46].

$t$ (h)	$\theta=0.24$	$\theta=0.31$	$\theta=0.38$
0.5	5.31e-02	5.31e-02	5.69e-02
1.0	-2.46e-03	2.175e-02	5.10e-02
1.5	-5.70e-02	-1.41e-02	4.55e-02
2.0	-9.69e-02	-3.88e-02	4.91e-02

An exact solution for constant flux infiltration with dry initial condition $\Theta(z,0)=0$ has been derived in [32] and further used to verify the numerical solution provided by a pressure formulation of the Richards equation [47]. The solution solves Eq.(14) with coefficient given by Fujita’s model [12],

D(\Theta)=D_{0}/(1-v\Theta)^{2},

where $D_{0}$ and $v$ are positive constants. Since $\Theta(z,0)=0$ implies $\psi(\Theta(z,0))=\infty$ as initial condition for the numerical scheme in pressure formulation, the singularity was avoided in [47] by considering $\Theta(z,0)=3.4483\cdot 10^{-6}$ as a numerical simulation parameter. As for the GRW scheme (6-8), we would have $K(\Theta(z,0))=0$ and, according to (5), the condition $r^{s}_{i\pm 1/2,k}\leq 1/2$ implies $\Delta t=\infty$ , for finite $\Delta z$ . Using the same initial $\psi$ as in [47] requires a very fine discretization which would slow down considerably the computation. Therefore, we opt for the direct approach of solving (14) as a diffusion equation with drift coefficient defined by $V(\Theta)=dK(\Theta)/d\Theta$ . The latter will be computed analytically from the parameterization $K(\Theta)$ used in [47].

Proceeding as in Section 2, we start with a forward-time centered-space finite difference discretization of Eq. (14),

	$\displaystyle\Theta_{i,k+1}-\Theta_{i,k1}=$	$\displaystyle+\frac{\Delta t}{{\Delta z}^{2}}\left[D(\Theta_{i+1/2,k})(\Theta_% {i+1,k}-\Theta_{i,k})-D(\Theta_{i-1/2,k})(\Theta_{i,k}-\Theta_{i-1,k})\right]$
		$\displaystyle-\frac{\Delta t}{{2\Delta z}}V_{i,k}(\Theta_{i+1,k}-\Theta_{i-1,k% }),$

we approximate the solution by a distribution of $\mathcal{N}$ particles on a regular lattice, $\Theta_{i,k}\approx n_{i,k}/\mathcal{N}$ , and end up with

	$\displaystyle n_{i,k+1}=$	$\displaystyle[1-(r_{i+1/2,k}+r_{i-1/2,k})]n_{i,k}$
		$\displaystyle+\frac{1}{2}(r_{i+1/2,k}-v_{i,k})n_{i+1,k}+\frac{1}{2}(r_{i-1/2,k% }+v_{i,k})n_{i-1,k}.$		(15)

The dimensionless parameters in Eq. (3.2) are given by

r_{i\pm 1/2,k}=\frac{2\Delta t}{{\Delta z}^{2}}D_{i\pm 1/2,k},\;\;v_{i,k}=% \frac{\Delta t}{\Delta z},\;\;r_{i\pm 1/2,k}\leq 1,\;\;|v_{i,k}|\leq r_{i\pm 1% /2,k}.

Equation (3.2) sums up contributions of random walkers jumping on the lattice according to the rule

n_{j,k}=\delta n_{j\mid j,k}+\delta n_{j-1\mid j,k}+\delta n_{j+1\mid j,k},

(16)

which defines a biased global random walk algorithm (BGRW) [40, Sect. 3.3.3]. The numbers of particles $\delta n$ in (16) are binomial random variables determined by the same procedure as in Section 2 and their ensemble averages verify

\overline{\delta n_{j\mid j,k}}=[1-(r_{i+1/2,k}+r_{i-1/2,k})]\overline{n_{i,k}% },\;\;\overline{\delta n_{j\pm 1\mid j,k}}=\frac{1}{2}(r_{i\pm 1,k}\mp v_{i,k}% )\overline{n_{i,k}}.

Following [47], we set on the top boundary the constant flux condition $Q=q/(\theta_{sat}-\theta_{res})=0.2759$ cm/min, with $\theta_{sat}=0.35$ , $\theta_{res}=0.06$ , and consider the constant parameters $D_{0}=2.75862$ cm²/min and $v=0.85$ of the Fujita’s model. The BGRW results for the final time $T=0.3625$ min, obtained with $\Delta z=10^{-2}$ cm and $\Delta t$ between $1.51\cdot 10^{-5}$ min and $1.23\cdot 10^{-5}$ min, are compared in Table 3 with the analytical solution presented in [47, Table 1].

Table 3: GRW solution

\theta(z,t)

compared to
the analytical solution from [32].

$z$ (cm)	$\tilde{\theta}(z,t)$	$\theta(z,t)$	$(\theta-\tilde{\theta})/\tilde{\theta}$
0	0.0907	0.0929	2.53e-02
-0.2	0.0861	0.0884	2.68e-02
-0.4	0.0819	0.0842	2.70e-02
-0.6	0.0782	0.0802	2.60e-02
-0.8	0.0748	0.0766	2.40e-02
-1.0	0.0719	0.0734	2.12e-02
-2.0	0.0631	0.0635	6.40e-03

The tests for unsaturated one-dimensional flows presented above are completed in [42, Sect. 5.2.4] by convergence investigations and estimations of convergence order of the GRW algorithms for fully coupled nonlinear flow and transport problems for saturated/unsaturated porous systems.

4 Two-dimensional GRW solutions

4.1 Two-dimensional GRW algorithm for flow in soils and aquifers

In two spatial dimensions the pressure head $\psi(x,z,t)$ satisfies the equation

\frac{\partial}{\partial t}\theta({\psi})-\nabla\cdot\left[K(\theta(\psi)% \nabla(\psi+z)\right]=0.

(17)

The two-dimensional GRW algorithm on regular staggered grids ( $\Delta x=\Delta z$ ) which approximates the solution of (17) by computational particles, $\psi\approx na/\mathcal{N}$ , is constructed similarly to (6-8). The solution at iteration $s+1$ is obtained by gathering particles from neighboring sites according to

$\displaystyle n_{i,j,k}^{s+1}$	$\displaystyle=$	$\displaystyle\left[1-\left(r_{i+1/2,j,k}^{s}+r_{i-1/2,j,k}^{s}+r_{i,j+1/2,k}^{% s}+r_{i,j-1/2,k}^{s}\right)\right]n_{i,j,k}^{s}$	(18)
		$\displaystyle+r_{i+1/2,j,k}^{s}n_{i+1,j,k}+r_{i-1/2,j,k}^{s}n_{i-1,j,k}$
		$\displaystyle+r_{i,j+1/2,k}^{s}n_{i,j+1,k}+r_{i,j-1/2,k}^{s}n_{i,j-1,k}+\left% \lfloor\mathcal{N}f^{s}\right\rfloor,$

where the source term is defined as $f^{s}=\left(r_{i,j+1/2,k}^{s}-r_{i,j-1/2,k}^{s}\right)\Delta z-\left[\theta(n_% {i,j,k}^{s})-\theta(n_{i,j,k-1})\right]/L$ . The two-dimensional GRW rule which at time $k$ moves particles from sites $(l,m)$ to neighboring sites $(l\mp 1,m\mp 1)$ reads as follows,

n_{l,m,k}^{s}=\delta n_{l,m|l,m,k}^{s}+\delta n_{l-1,m|l,m,k}^{s}+\delta n_{l+% 1,m|l,m,k}^{s}+\delta n_{l,m-1|l,m,k}^{s}+\delta n_{l,m+1|l,m,k}^{s}.

(19)

For consistency with (18), the numbers of particles $\delta n^{s}$ verify in the mean

	$\displaystyle\overline{\delta n_{l,m\|l,m,k}^{s}}=\left[1-\left(r_{l-1/2,m,k}^{% s}+r_{l+1/2,m,k}^{s}+r_{l,m-1/2,k}^{s}+r_{l,m+1/2,k}^{s}\right)\right]% \overline{n_{l,m,k}^{s}}$
	$\displaystyle\overline{\delta n_{l\mp 1,m\|l,m,k}^{s}}=r_{l\mp 1/2,m,k}^{s}% \overline{n_{l,m,k}^{s}}$
	$\displaystyle\overline{\delta n_{l,m\mp 1\|l,m,k}^{s}}=r_{l,m\mp 1/2,k}^{s}% \overline{n_{l,m,k}^{s}}.$		(20)

The parameters $r_{l\mp 1/2,m,k}^{s}$ and $r_{l,m\mp 1/2,k}^{s}$ , defined by relations similar to (5), are dimensionless positive real numbers. They represent biased jump probabilities on the four allowed spatial directions of the GRW lattice and are constraint by the first relation (4.1) such that their sum be less or equal to one. A sufficient condition would be that each of them verifies $r\leq 1/4$ .

The binomial random variables variables $\delta n$ are approximated in the same way as in the one-dimensional case. By giving up the particle indivisibility, one obtains deterministic GRW algorithms which represent the solution $n$ by real numbers and use the unaveraged relations (4.1) for the computation of the $\delta n$ terms. In the following we use this deterministic implementation of the GRW algorithm to compute flow solutions for unsaturated/saturated porous media.

Remark 3

After disregarding the index $k$ and letting $L=1/a$ , $\theta=const$ , the algorithm (18-4.1) becomes a transient scheme to solve the equation governing flows in saturated porous media [2, 41] (see also Remark 2).

4.2 GRW algorithms for two-dimensional fully coupled flow and surfactant transport

Let the pressure $\psi(x,z,t)$ and the concentration $c(x,z,t)$ solve the equations of the following model of fully coupled flow and surfactant transport in unsaturated/saturated porous media [20, 18],

\frac{\partial}{\partial t}\theta(\psi,c)-\nabla\cdot\left[K(\theta(\psi,c)% \nabla(\psi+z)\right]=0,

(21)

\frac{\partial}{\partial t}\left[\theta(\psi,c)c\right]-\nabla\cdot\left[D% \nabla c-\mathbf{q}c\right]=R(c),

(22)

where $\mathbf{q}=-K(\theta(\psi,c)\nabla(\psi+z)$ is the water flux (Darcy velocity) and $R(c)$ is a nonlinear reaction term. Equations (21-22) are coupled in both directions through the nonlinear functions $\theta(\psi,c)$ and $\theta(\psi,c)c$ . The pressure equation (21) is solved with the GRW $L$ -scheme described in the previous subsection, with a slight modification due to the dependence of $\theta$ on both $\psi$ and $c$ . New algorithms are needed instead to solve the coupled, nonlinear transport equation (22).

4.2.1 Biased GRW algorithm for transport problems

To derive a GRW algorithm for the transport equation, we start with a backward-time central-space finite difference scheme for Eq. (22). Considering a diagonal diffusion tensor with constant components $D_{1}$ and $D_{2}$ , and denoting by $U$ and $V$ the components of the Darcy velocity along the horizontal axis $x$ and the vertical axis $z$ , by $\Delta t$ the time step, and by $\Delta x$ and $\Delta z$ the spatial steps, the scheme reads as

	$\displaystyle\theta(\psi_{i,j,k},c_{i,j,k})c_{i,j,k}-\theta(\psi_{i,j,k-1},c_{% i,j,k-1})c_{i,j,k-1}=$
	$\displaystyle-\frac{\Delta t}{2\Delta x}\left(U_{i+1,j,k}c_{i+1,j,k}-U_{i-1,j,% k}c_{i-1,j,k}\right)-\frac{\Delta t}{2\Delta z}\left(V_{i,j+1,k}c_{i,j+1,k}-V_% {i,j-1,k}c_{i,j-1,k}\right)$
	$\displaystyle+\frac{D_{1}\Delta t}{\Delta x^{2}}\left(c_{i+1,j,k}-2c_{i,j,k}+c% _{i-1,j,k}\right)+\frac{D_{2}\Delta t}{\Delta z^{2}}\left(c_{i,j+1,k}-2c_{i,j,% k}+c_{i,j-1,k}\right)+R(c_{i,j,k})\Delta t=$
	$\displaystyle-\left(\frac{2D_{1}\Delta t}{\Delta x^{2}}+\frac{2D_{2}\Delta t}{% \Delta z^{2}}\right)c_{i,j,k}$
	$\displaystyle+\left(\frac{D_{1}\Delta t}{\Delta x^{2}}-\frac{\Delta t}{2\Delta x% }U_{i+1,j,k}\right)c_{i+1,j,k}+\left(\frac{D_{1}\Delta t}{\Delta x^{2}}+\frac{% \Delta t}{2\Delta x}U_{i-1,j,k}\right)c_{i-1,j,k}$
	$\displaystyle+\left(\frac{D_{2}\Delta t}{\Delta z^{2}}-\frac{\Delta t}{2\Delta z% }V_{i,j+1,k}\right)c_{i,j+1,k}+\left(\frac{D_{2}\Delta t}{\Delta z^{2}}+\frac{% \Delta t}{2\Delta z}V_{i,j-1,k}\right)c_{i,j-1,k}+R(c_{i,j,k})\Delta t.$		(23)

Next, similarly to the scheme for the flow equation, we add a regularization term $L(c^{s+1}_{i,j,k}-c^{s}_{i,j,k})$ in Eq. (4.2.1), define the dimensional parameters

r_{x}=\frac{2D_{1}\Delta t}{L\Delta x^{2}},\;\;r_{z}=\frac{2D_{2}\Delta t}{L% \Delta z^{2}},\;\;u_{i\pm 1,j,k}^{s}=\frac{\Delta t}{L\Delta x}U_{i\pm 1,j,k}^% {s},\;\;v_{i,j\pm 1,k}^{s}=\frac{\Delta t}{L\Delta z}V_{i,j\pm 1,k}^{s},

(24)

approximate the concentration by the density of the number of computational particles, $c^{s}_{i,j,k}\approx n^{s}_{i,j,k}/\mathcal{N}$ , and finally we obtain

$\displaystyle n_{i,j,k}^{s+1}=$	$\displaystyle\left[1-\left(r_{x}+r_{z}\right)\right]n_{i,j,k}^{s}$
	$\displaystyle+\frac{1}{2}\left(r_{x}-u_{i+1,j,k}^{s}\right)n_{i+1,j,k}^{s}+% \frac{1}{2}\left(r_{x}+u_{i-1,j,k}^{s}\right)n_{i-1,j,k}^{s}$
	$\displaystyle+\frac{1}{2}\left(r_{z}-v_{i,j+1,k}^{s}\right)n_{i,j+1,k}^{s}+% \frac{1}{2}\left(r_{z}+v_{i,j-1,k}^{s}\right)n_{i,j-1,k}^{s}+\left\lfloor% \mathcal{N}g^{s}\right\rfloor,$	(25)

where $g^{s}=R(n_{i,j,k}^{s})\Delta t/L-\left[\theta(\psi_{i,j,k}^{s},n_{i,j,k}^{s})n% _{i,j,k}^{s}-\theta(\psi_{i,j,k-1},n_{i,j,k-1})n_{i,j,k-1}\right]/L$ , with $\psi$ approximated by the distribution of particles in the flow solver for Eq. (21). Note that the definition of the dimensionless numbers (24) implies that the parameter $L$ has to be a dimensionless number as well.

The contributions to $n_{i,j,k}^{s+1}$ in Eq. (4.2.1) are obtained with the BGRW algorithm

n_{l,m,k}^{s}=\delta n_{l,m\mid l,m,k}^{s}+\delta n_{l-1,m\mid l,m,k}^{s}+% \delta n_{l+1,m\mid l,m,k}^{s}+\delta n_{l,m-1\mid l,m,k}^{s}+\delta n_{l,m+1% \mid l,m,k}^{s},

(26)

where, for consistency with the finite difference scheme (4.2.1), the quantities $\delta n$ verify in the mean

	$\displaystyle\overline{\delta n_{l,m\mid l,m,k}^{s}}=\left[1-\left(r_{x}+r_{z}% \right)\right]\mbox{ }\overline{n_{i,j,k}^{s}},\;\;\;$	$\displaystyle\overline{\delta n_{l\pm 1,m\mid l,m,k}^{s}}=\frac{1}{2}(r_{x}\mp u% _{l,m,k}^{s})\overline{n_{l,m,k}^{s}},$
		$\displaystyle\overline{\delta n_{l,m\pm 1\mid l,m,k}^{s}}=\frac{1}{2}(r_{z}\mp v% _{l,m,k}^{s})\overline{n_{l,m,k}^{s}}.$		(27)

The binomial random variables $\delta n$ used in the BGRW algorithm are approximated similarly to the algorithms described in the previous sections, by summing up to unity reminders of multiplication and floor operations. A deterministic BGRW algorithm can be obtained, similarly to the flow solver presented in Section 4.1 above, by giving up the particle’s indivisibility and using the un-averaged relations (4.2.1). However, for the computations presented in the next section, we use a randomized implementation of the BGRW algorithm.

As follows from (4.2.1), the BGRW algorithm is subject to the following restrictions

r_{x}+r_{z}\leq 1,\mbox{ }\left|u_{l,m,k}^{s}\right|\leq r_{x},\mbox{ }\left|v% _{l,m,k}^{s}\right|\leq r_{z}.

(28)

Remark 4

The constraints (28) impose a limitation on the maximum allowable value of the local Péclet number. Assume a constant flow velocity $-V$ and a constant diffusion coefficient $D$ . Then, according to (28) and (24), the condition $v\leq r$ implies $\text{P\'{e}}=V\Delta z/D\leq 2$ .

Remark 5

Taking into account that the iterations start with $n_{i,j,k-1}$ , setting $L=1$ , $\theta=1$ , and dropping the superscripts $s$ , the relation (4.2.1) becomes

$\displaystyle n_{i,j,k}=$	$\displaystyle\left[1-(r_{x}+r_{z})\right]n_{i,j,k-1}$
	$\displaystyle+\frac{1}{2}\left(r_{x}-u_{i+1,j,k-1}\right)n_{i+1,j,k-1}+\frac{1% }{2}\left(r_{x}+u_{i-1,j,k-1}\right)n_{i-1,j,k-1}$
	$\displaystyle+\frac{1}{2}\left(r_{z}-v_{i,j+1,k-1}\right)n_{i,j+1,k-1}+\frac{1% }{2}\left(r_{z}+v_{i,j-1,k-1}\right)n_{i,j-1,k-1}+\left\lfloor\mathcal{N}R(n_{% i,j,k-1})\Delta t\right\rfloor.$	(29)

Relation (5), together with (26-28), define a BGRW algorithm for (decoupled) reactive transport described by Eq. (22) with $\theta(\psi,c)=1$ .

4.2.2 Unbiased GRW algorithm for transport problems

The unbiased GRW algorithm is obtained by globally moving groups of particles according to the rule

$\displaystyle n_{i,j,k}^{s}=$	$\displaystyle\;\delta n_{i+u_{i,j,k}^{s},j+v_{i,j,k}^{s}\mid i,j,k}^{s}$	(30)
	$\displaystyle+\delta n_{i+u_{i,j,k}^{s}+d,j+v_{i,j,k}^{s}\mid i,j,k}^{s}+% \delta n_{i+u_{i,j,k}^{s}-d,j+v_{i,j,k}^{s}\mid i,j,k}^{s}$
	$\displaystyle+\delta n_{i+u_{i,j,k}^{s},j+v_{i,j,k}^{s}+d\mid i,j,k}^{s}+% \delta n_{i+u_{i,j,k}^{s},j+v_{i,j,k}^{s}-d\mid i,j,k}^{s},$

where $d$ is a constant amplitude of diffusion jumps and the dimensionless variables $r_{x}$ , $r_{z}$ , $u$ and $v$ are defined similarly to (24) by

r_{x}=\frac{2D_{1}\Delta t}{L(d\Delta x)^{2}},\;\;r_{z}=\frac{2D_{2}\Delta t}{% L(d\Delta z)^{2}},\;\;u_{i,j,k}^{s}=\left\lfloor\frac{\Delta t}{L\Delta x}U_{i% ,j,k}^{s}+0.5\right\rfloor,\;\;v_{i,j,k}^{s}=\left\lfloor\frac{\Delta t}{L% \Delta z}V_{i,j,k}^{s}+0.5\right\rfloor.

(31)

The particles distribution is updated at every time step by

n_{l,m,k}^{s+1}=\delta n_{l,m,k}^{s}+\sum_{i\neq l,j\neq m}\delta n_{l,m\mid i% ,j,k}^{s}+\left\lfloor\mathcal{N}g^{s}\right\rfloor.

(32)

The averages over GRW runs of the terms from (30) are now related by

	$\displaystyle\overline{\delta n_{i+u_{i,j,k}^{s},j+v_{i,j,k}^{s}\mid i,j,k}^{s% }}=\left[1-\left(r_{x}+r_{z}\right)\right]\mbox{ }\overline{n_{i,j,k}^{s}},$
	$\displaystyle\overline{\delta n_{i+u_{i,j,k}^{s}\pm d,j+v_{i,j,k}^{s}\mid i,j,% k}^{s}}=\frac{r_{x}}{2}\hskip 2.84544pt\overline{n_{i,j,k}^{s}},$
	$\displaystyle\overline{\delta n_{i+u_{i,j,k}^{s},j+v_{i,j,k}^{s}\pm d\mid i,j,% k}^{s}}=\frac{r_{z}}{2}\hskip 2.84544pt\overline{n_{i,j,k}^{s}}.$		(33)

Comparing with the BGRW relations (4.2.1), we remark that (31) defines unbiased jump probabilities $r_{x}/2$ and $r_{y}/2$ on the two spatial directions.

The unbiased GRW algorithm for decoupled transport is obtained by letting $L=1$ and dropping the superscripts $s$ (see also Remark 5).

The binomial random variables $\delta n$ used in the unbiased GRW algorithm are approximated by the procedure used for the flow solver and for the BGRW algorithm presented in the previous subsection. For fixed space steps, the time step is chosen such that the dimensionless parameters $u_{i,j,k}^{s}$ and $v_{i,j,k}^{s}$ take integer values larger than unity which ensure the desired resolution of the velocity components [40, Sect. 3.3.2.1]. Further, the jumps’ amplitude $d$ is chosen such that the jump probabilities verify the constraint $r_{x}+r_{z}\leq 1$ , imposed by the first relation (33).

The unbiased GRW, as well as the BGRW algorithm introduced in Section 4.2.1 above, have been tailored to solve problems with constant diffusion coefficients, as those considered in Sections 5.2 and 6.3 below. In case of diagonal diffusion tensors with space-time variable coefficients $D_{1}$ and $D_{2}$ , the algorithms for the transport problem are straightforwardly obtained by assigning to $r_{x}$ and $r_{z}$ superscripts $s$ and appropriate subscripts $i, j, k$ .

5 Validation of the two-dimensional GRW algorithms

5.1 GRW flow solutions

For the beginning, we conduct verification tests of the GRW flow code by comparisons with an analytical solution and compute numerical estimates of the order of convergence. The results are further compared with those obtained by a TPFA code implemented in the MRST software [22, 18]. The two codes are tested by solving a problem with manufactured solution previously considered in [30]. The domain is the unit square $[0,1]\times[0,1]$ and the final time is $T=1$ . The manufactured solution for the pressure head $\psi_{m}$ is given by

\psi_{m}(x,z,t)=-t\ x\ (x-1)\ z\ (z-1)\ -\ 1.

(34)

The water content $\theta$ and the conductivity $K$ are expressed as

\theta(\psi)=\frac{1}{1-\psi},\quad K(\theta(\psi))=\psi^{2}.

(35)

The analytical solution (34) defines the boundary and initial conditions and induces a source term $f$ , computed analytically from Eq. (17) with parameters given by the expressions (35).

We start the computations on a uniform mesh with $\Delta x=\Delta z=0.1$ and halve the mesh size step three times successively. The accuracy of the numerical solutions, at the final time $t=T$ , is quantified by the $L^{2}$ norm $\varepsilon_{l}=\|\psi^{(l)}-\psi_{m}\|$ , $l=1,\ldots,4$ , where $l=1$ corresponds to the original mesh. The estimated order of convergence (EOC) that describes the decrease of the error in logarithmic scale is computed according to

EOC=\log\left(\frac{\varepsilon_{l}}{\varepsilon_{l+1}}\right)/\log(2),\quad l% =1,\ldots,3.

(36)

The computations with the TPFA code start with a time step $\Delta t=0.1$ which is also halved at each refinement of the mesh. The parameters of the convergence indicator (3) are set to $\varepsilon_{a}=10^{-6}$ and $\varepsilon_{r}=0$ . Finally, the linearization parameter $L$ is set equal to $1/2$ and the convergence of the $L$ -scheme is achieved after circa 100 iterations per time step, independently of the mesh size.

In the GRW computations we use the same spatial refinement of the grid and tolerances $\varepsilon_{a}$ and $\varepsilon_{r}$ as above but, according to (5), we have to use adaptive time steps $\Delta t=\mathcal{O}(\Delta z^{1/2})$ (see discussion in Section 3.1). The convergence criterion (3) is already fulfilled by the GRW $L$ -scheme with parameter $L=1$ for numbers of iterations increasing from $s=2$ to $s=5$ as the space step decreases. The accuracy $\varepsilon_{l}$ instead is strongly influenced by $L$ . For $L<800$ the $\varepsilon_{l}$ values may increase with the refinement of the mesh, leading to negative EOC, that is, the GRW solution does not converge to the exact solution $\psi_{m}$ . However, it is found that the increase of $\varepsilon_{l}$ is prevented by using a sufficiently large parameter $L$ .

The results presented in Table 4 indicate the convergence of order 1 in space for TPFA and of order 2 for the GRW solutions. The higher order of convergence also leads to much smaller errors of the GRW code after the first refinement of the mesh.

Table 4: Estimated order of convergence of the TPFA and GRW flow solvers.

	$\varepsilon_{1}$	EOC	$\varepsilon_{2}$	EOC	$\varepsilon_{3}$	EOC	$\varepsilon_{4}$
TPFA	8.45e-03	0.94	4.40e-03	0.97	2.25e-03	0.97	1.15e-03
GRW (L=800)	7.20e-03	2.24	1.52e-03	3.21	1.65e-04	0.50	1.17e-04
GRW (L=1000)	9.24e-03	2.22	1.99e-03	2.83	2.80e-04	1.66	8.84e-05
GRW (L=1200)	8.89e-03	2.23	1.90e-03	2.80	2.72e-04	2.14	6.16e-05

Further, we solve the benchmark problem from [24, Sect. 4.2], which describes the recharge of a groundwater reservoir from a drainage trench in a two-dimensional geometry. The groundwater table is fixed by a Dirichlet boundary condition on the right hand side. The drainage process is driven by a Dirichlet boundary condition changing in time on the upper boundary of $\Omega$ .

The precise structure of the domain is defined by

\begin{split}\Omega&=(0,2)\times(0,3),\\ \Gamma_{D_{1}}&=\{(x,z)\in\partial\Omega\ |\ x\in[0,1]\wedge z=3\},\\ \Gamma_{D_{2}}&=\{(x,z)\in\partial\Omega\ |\ x=2\wedge z\in[0,1]\},\\ \Gamma_{D}&=\Gamma_{D_{1}}\cup\Gamma_{D_{2}},\\ \Gamma_{N}&=\partial\Omega\setminus\Gamma_{D}.\end{split}

The Dirichlet and Neumann boundary conditions on $\Gamma_{D}$ and $\Gamma_{N}$ , respectively, as well as the initial condition consisting of hydrostatic equilibrium are specified as follows:

\begin{split}&\psi(x,z,t)=\begin{cases}-2+2.2t/\Delta t_{D},&\quad\ \text{on}% \ \Gamma_{D_{1}},T\leq\Delta t_{D},\\ 0.2,&\quad\ \text{on}\ \Gamma_{D_{1}},T>\Delta t_{D},\\ 1-z,&\quad\ \text{on}\ \Gamma_{D_{2}},\end{cases}\\ &-K(\theta(\psi(x,z,t))\nabla(\psi(x,z,t)+z)\cdot\mathbf{n}\ =\ 0,\quad\text{% on}\ \Gamma_{N},\\ &\psi(x,z,0)\ =\ 1-z,\quad\text{on}\ \Omega,\\ \end{split}

where $\mathbf{n}$ represents the outward pointing normal vector.

We consider here two sets of soil parameters, presented in Table 5, which correspond to a silt loam and a Beit Netofa clay, respectively.

Table 5: Simulation parameters

	Silt loam	Beit Netofa clay
Vam Genuchten parameters:
$\theta_{sat}$	0.396	0.446
$\theta_{res}$	0.131	0
$\alpha$	0.423	0.152
$n$	2.06	1.17
$K_{sat}$	$4.96\cdot 10^{-2}$	$8.2\cdot 10^{-4}$
Time parameters:
$\Delta t_{D}$	1/16	1
$\Delta t$	1/48	1/3
$T$	3/16	3

The time unit is 1 day and spatial dimensions are given in meters. Furthermore, we consider a regular mesh consisting of 651 nodes (i.e., $\Delta x=\Delta z=0.1$ ).

By setting the stabilization parameters to $L=0.5$ for loam and for $L=0.12$ for clay, the convergence criterion (3) with $\varepsilon_{a}=\varepsilon_{r}=5\cdot 10^{-6}$ is fulfilled after about 120 iterations of the GRW $L$ -scheme, for both soil models (Figs. 10 and 10). The results shown in Figs. 12 and 12 are as expected for this benchmark problem (see [34, 24]): the drainage process in the clay soil is much slower, so that the pressure distribution after three days is similar to that established in the loam soil after 4.5 hours.

The results obtained with the TPFA $L$ -scheme, with $L=1$ for both soil models, are used as reference to compute the relative errors $\varepsilon_{\psi}$ , $\varepsilon_{\theta}$ , $\varepsilon_{q_{x}}$ , and $\varepsilon_{q_{z}}$ shown in Table 6. One remarks that $\varepsilon_{\psi}$ and $\varepsilon_{\theta}$ are close to the corresponding errors for the one-dimensional case presented in Table 1, but $\varepsilon_{q_{x}}$ and $\varepsilon_{q_{z}}$ are one order of magnitude larger than $\varepsilon_{q}$ in shown in Table 1. A possible explanation could be the occurrence of the numerical diffusion in the flow TPFA code (see discussion at the end of Section 5.2.3 below). The computational times of the GRW code are 1 second and 1.6 seconds for loam and clay cases, respectively. The times of the TPFA runs, on the same computer, are one order of magnitude larger, i.e., 25 seconds and 38 seconds, respectively.

Table 6: Comparison of GRW and TPFA solutions of the
flow benchmark problem.

	$\varepsilon_{\psi}$	$\varepsilon_{\theta}$	$\varepsilon_{q_{x}}$	$\varepsilon_{q_{z}}$
loam	5.73e-02	4.00e-03	2.30e-01	1.04e-01
clay	5.48e-02	6.71e-04	4.73e-01	1.14e-01

5.2 GRW/BGRW solutions for fully coupled flow and transport problems

5.2.1 Code verification tests

The code verification tests for coupled flow and transport problems are conducted similarly to those for the flow solver presented in the previous subsection, by considering, along with the exact flow solution (34), the exact solution for the concentration field given by

c_{m}(x,z,t)=t\ x\ (x-1)\ z\ (z-1)\ +\ 1.

(37)

After setting $R=0$ and $D=1$ , the coupled system of equations (21-22) is solved in the unit square for a total time $T=1$ , with source terms, initial conditions, and boundary conditions resulted from the exact solutions (34) and (37) with a new parameterization given by

\theta(\psi,c)=\frac{1}{1-\psi-c/10},\quad K(\theta(\psi))=\psi^{2}.

(38)

The GRW flow-algorithm (18-4.1), with $\theta$ and $K$ given by (38), is coupled with the BGRW transport-algorithm (4.2.1-28) initialized with $\mathcal{N}=10^{24}$ particles into an alternating splitting scheme [18]. The approach alternates iterations of flow and transport solvers until the convergence criterion (3) with $\varepsilon_{a}=10^{-6}$ and $\varepsilon_{r}=0$ is fulfilled by the numerical solutions for both $\psi$ and $c$ . In order to highlight the approach to the convergence order 2, the stabilization parameters of the flow and the transport solvers are set to $L_{p}=L_{c}=100$ . The GRW results presented in Tables 7 and 8 are compared with results obtained with a TPFA solver applying the same alternating linearized splitting procedure with parameters $L_{p}=L_{c}=1$ which ensure the convergence of order 1.

The GRW flow solver approximates the Darcy velocity by centered differences only in the interior $\Omega$ of the computational domain. Therefore, the velocity $\mathbf{q}|_{\partial\Omega}$ , needed to compute the number of biased jumps from the boundary $\partial\Omega$ in the BGRW relation (4.2.1) has to be provided in some way. The straightforward approach is to compute the velocity by using an approximate forward finite difference discretization of Darcy’s law. Another option is to extend on the boundary the velocity from the first neighboring interior site. Thanks to the manufactured solution (34) on which the code verification test is based, we also have the exact velocity computed analytically. The latter allows accuracy assessments for the above approximations. We note that the GRW results for the pressure solver obtained with analytical, approximate, and extend $\mathbf{q}|_{\partial\Omega}$ are identical in the precision of three significant digits (Table 7). For the concentration solutions (Table 8), we note the remarkably good performance of approximate and extended $\mathbf{q}|_{\partial\Omega}$ .

Table 7: Estimated order of convergence of the TPFA and GRW solvers: pressure solutions.

	$\varepsilon_{1}$	EOC	$\varepsilon_{2}$	EOC	$\varepsilon_{3}$	EOC	$\varepsilon_{4}$
TPFA	8.14e-03	0.93	4.27e-03	0.95	2.20e-03	0.97	1.12e-03
GRW	3.71e-03	2.02	9.18e-04	1.94	2.40e-04	1.45	8.78e-05

Table 8: Estimated order of convergence of the TPFA and GRW solvers: concentration solutions.

	$\varepsilon_{1}$	EOC	$\varepsilon_{2}$	EOC	$\varepsilon_{3}$	EOC	$\varepsilon_{4}$
TPFA	6.26e-03	0.83	3.52e-03	0.89	1.90e-03	0.91	1.01e-03
GRW (analytical $\mathbf{q}\|_{\partial\Omega}$ )	3.92e-03	2.00	9.78e-04	1.83	2.74e-04	1.05	1.32e-04
GRW (approximate $\mathbf{q}\|_{\partial\Omega}$ )	4.72e-03	1.99	1.19e-03	1.85	3.29e-04	1.17	1.46e-04
GRW ( $\mathbf{q}\|_{\partial\Omega}$ from $\mbox{int}(\Omega)$ )	5.26e-03	2.00	1.31e-03	1.87	3.59e-04	1.23	1.53e-04

5.2.2 Estimates of numerical diffusion

The small errors shown in Table 8 indicate that the numerical diffusion in solving the transport step of the coupled problem does not play a significant role. This is somewhat expected for the small Péclet numbers of order Pé= $10^{-2}$ encountered in these computations. But for the numerical setup of the benchmark problem presented in Section 5.1 and realistic transport parameters Pé can be significantly larger than unity. Therefore we proceed to estimate the numerical diffusion of the codes compared here by following the procedure used in [29].

Table 9: Estimation of numerical diffusion for
BGRW, GRW and TPFA codes.

	$\Delta x$	$T/\Delta t$	Pé	$\varepsilon_{D_{x}}$	$\varepsilon_{D_{z}}$
BGRW	0.1	2	3.31	7.55e-02	2.60e-01
	0.05	9	1.65	1.90e-16	1.48e-15
	0.01	239	0.33	4.16e-16	1.02e-15
	0.005	960	0.17	2.93e-15	3.63e-15
GRW	0.1	4	3.31	1.94e-16	6.14e-16
	0.05	4	1.65	6.60e-17	8.05e-16
	0.01	19	0.33	1.94e-16	4.79e-16
	0.005	39	0.17	2.10e-15	8.92e-16
TPFA	0.1	5	3.31	9.16e-03	1.99e-01
	0.05	10	1.65	4.69e-03	9.94e-02
	0.01	50	0.33	9.58e-04	1.99e-02
	0.005	100	0.17	5.38e-04	9.89e-03

We consider the analytical Gaussian solution $c(x,z,t)$ of Eq. (22) with $\theta=1$ , $R=0$ , and constant coefficients $D=0.001$ and $V=-0.0331$ , corresponding to the Cauchy problem with a Dirac initial concentration pulse located at the coordinates (1,2.1). The constant velocity $V$ , oriented downwards along the $z$ -axis, is the steady-state solution of the benchmark flow problem from Section 5.1 with $K=K_{sat}$ corresponding to the loam soil, initial condition $\psi(x,z,0)=1-z/3$ , Dirichlet boundary conditions $\psi(x,0,t)=1$ , $\psi(x,3,t)=0$ , and no-flow Neumann conditions on the vertical boundaries. The initial condition $c(x,z,0)$ is the same Gaussian function evaluated at $t=1$ and the final time is $T=3$ . For decreasing mesh sizes $\Delta x$ and $\text{P\'{e}}=V\Delta x/D$ , the number of time steps was restricted by the requirement that the support of the numerical solution does not extend beyond the boundaries $\partial\Omega$ (to mimic diffusion in unbounded domains). The effective diffusion coefficients $D_{x}$ and $D_{z}$ are computed from the spatial moments along the $x$ - and $z$ -directions of the numerical solution (see [29, Eqs. (38-41)]). The numerical diffusion is estimated by relative errors $\varepsilon_{D_{x}}=|D_{x}-D|/D$ and $\varepsilon_{D_{z}}=|D_{z}-D|/D$ averaged over the time interval $[0,T]$ . Table 9 shows that while the TPFA results are strongly influenced by the mesh size, similarly to the finite-volume results from [29], the unbiased GRW algorithm is practically unconditionally-free of numerical diffusion. The BGRW algorithm is also free of numerical diffusion provided that $\text{P\'{e}}\leq 2$ (see also Remark 4). We also note that $\Delta x=0.05$ defines the coarsest grid acceptable for solving the benchmark problem for coupled flow and transport with BGRW and TPFA codes.

5.2.3 Fully coupled water flow and surfactant transport

In the following we solve the coupled flow and transport problem (21-22) by using the setup of the benchmark flow problem problem from Section 5.1 completed by parameters and initial/boundary conditions modeling a situation of coupled water flow and surfactant transport. The surfactant concentration in the domain $\Omega$ has a stratified distribution described by the plane $c(x,z,0)=z/1.2$ . Further, the concentration is set to $c=1$ on the Dirichlet boundary $\Gamma_{D_{1}}$ and to $c=0$ on $\Gamma_{D_{2}}$ , and no-flow Neumann conditions are imposed on the vertical boundaries.

The flow and transport are coupled in both directions through the van Genuchten-Mualem parameterization (12-13) with $\theta(\psi,c)=\theta(\gamma(c)\psi)$ , where $\gamma(c)=1/[1-b\ln(c/a+1)]$ models the concentration-dependent surface tension between water and air [20]. The constant parameters of $\gamma(c)$ are set to $a=0.44$ and $b=0.0046$ [18]. To describe a more realistic heterogeneous soil, the saturated conductivity $K_{sat}$ is modeled as a log-normal space random function with a small variance $\sigma^{2}=0.5$ and Gaussian correlation of correlation lengths $\lambda_{x}=0.1$ m and $\lambda_{z}=0.01$ m in horizontal and vertical directions, respectively. The $\ln K$ field is generated by summing up 100 random periodic modes with the Kraichnan algorithm presented in [40, Appendix C.3.1.2]. The diffusion coefficient is set to a constant value, $D=10^{-3}$ m/day, which is representative for soils and aquifers [29, 34, 40]. Following [18], the nonlinear reaction term is specified as $R(c)=10^{-3}c/(1+c)$ . Instead of using a fixed number of time steps, as in the flow benchmark presented in Section 5.1, now we fix the total time to $T=3$ days, set the intermediate time controlling the drainage process to $\Delta t_{D}=T/3$ , and keep the original time steps $\Delta t$ which ensure the appropriate resolution for contrasting fast and slow processes in loam and clay soils, respectively.

Preliminary tests showed that, in order to obtain an acceptable resolution of the velocity components in the benchmark setup, the unbiased GRW requires extremely fine discretizations with $\Delta x=\mathcal{O}(10^{-5})$ . Therefore the transport step is solved with the BGRW algorithm for the mesh size $\Delta x=0.05$ suggested by the above investigations on numerical diffusion. The velocity $\mathbf{q}|_{\partial\Omega}$ on boundaries is approximated by forward finite differences.

The convergence of the flow and transport $L$ -schemes using GRW algorithms requires relatively large linearization parameters, $L_{p}=L_{c}=20$ , for loam soil, and $L_{p}=L_{c}=100$ for clay soil models. These are two order of magnitude larger than for the decoupled-flow benchmark presented in Section 5.1, probably due to the increased complexity of the coupled problem. By setting the tolerances of the convergence criterion (3) to $\varepsilon_{a}=\varepsilon_{r}=5\cdot 10^{-6}$ the convergence is achieved after about 2000 iterations for the loam soil and about 14000 iterations for the clay soil (see Figs. 14 - 16).

The results obtained by coupling the GRW-flow and BGRW-transport solvers are presented in Figs. 18-26. The randomness of $K_{sat}$ is especially felt by the pressure distribution in the more permeable loam soil (Fig. 18), while in the clay soil the pressure remains almost stratified (Fig. 18). The same contrast is shown by the water content, with almost saturated loam soil (Fig. 20) and partially stratified saturation in the clay soil (Fig. 20). Since the Darcy velocity is proportional to the gradient of the random pressure, the heterogeneity of the advective component of the transport process is mainly manifest in the final distribution of the concentration in the loam and clay soils (compare Figs. 22 and Fig. 22). Significant differences between the loan and clay soils are also illustrated by the spatial distribution of the velocity components (Figs. 24 - 26).

The results obtained with the GRW/BGRW flow and transport solvers are compared with those provided by a TPFA code using $L_{p}=L_{c}=1$ , for both soils, and $L_{c}=2L_{p}$ . The convergence is achieved in reasonable computing times of 263 seconds (loam) and 177 seconds (clay) only when using the Anderson acceleration procedure [3, 45, 8]. Note that the GRW times on the same computer are of the same order of magnitude (526 and 178, respectively), without appealing to the acceleration procedure.

The errors for pressure, water content and velocity components shown in Table 10 are more or less similar to those for the flow benchmark problem given in Table 6. The difference of one order of magnitude between the $\varepsilon_{c}$ values for the two soils can be traced back to the amount of numerical diffusion of the TPFA transport solver (see Table 9). The estimated mean Péclet number for the loam soil, $\text{P\'{e}}\approx 1.3$ , is much larger than the value $\text{P\'{e}}\approx 4\cdot 10^{-3}$ estimated for the clay soil and can partially explain the larger $\varepsilon_{c}$ value in the first case. Since the pressure equation is essentially an advection-diffusion equation with velocity given by the derivatives of the coefficient $K$ [e.g., 14, 41], the errors $\varepsilon_{q_{x}}$ and $e_{q_{z}}$ , of order $~{}10^{-1}$ also could be produced by numerical diffusion, in the flow solver. In the setup of the benchmark problems, for both coupled flow and transport and decoupled flow, we estimate a mean Péclet number $\text{P\'{e}}\approx 0.9$ for both loam and clay soil models (for comparison, in the one dimensional case with smaller $\varepsilon_{q}$ , Pé was about 0.03 in Scenario 1 and 0.3 in Scenario 2). Since the flow and transport solvers implemented in MRST basically use the same TPFA finite volume method, we may expect that the flow solver produces a numerical diffusion comparable to that of the transport solver shown in Table 9.

Table 10: Comparison of GRW and TPFA solutions of the
coupled flow-transport benchmark problem.

	$\varepsilon_{\psi}$	$\varepsilon_{c}$	$\varepsilon_{\theta}$	$\varepsilon_{q_{x}}$	$\varepsilon_{q_{z}}$
loam	2.89e-02	4.79e-01	7.25e-05	3.15e-01	2.18e-01
clay	5.95e-02	3.77e-02	7.61e-04	3.66e-01	5.36e-01

A one-dimensional version of the benchmark problem for flow and surfactant transport can be readily obtained and solved with one-dimensional GRW algorithms [42, Sect. 5.2.4]. Even though the lateral heterogeneity of the two-dimensional benchmark is ignored, the main features are also revealed by the one-dimensional drainage model: the discrepancy between fast-loam and slow-clay flow and transport processes, the same intervals of variation of the solutions, and similar behavior on the vertical direction.

6 Two-dimensional GRW solutions for groundwater flow and transport at regional and field scales

For saturated aquifers ( $\theta=const$ ) Eq. (17) reduces to a linear equation solved by the steady state hydraulic head solution in $h(x,y)$ , under time independent boundary conditions. As noted in Remark 3, the GRW $L$ -scheme (18-4.1) becomes, in this case, a transient scheme for the linear flow equation. In the following examples, we consider flow problems formulated in two-dimensional domains, $(x,y)\in[0,L_{x}]\times[0,L_{y}]$ , with Dirichlet boundary conditions $h(0,y)=H_{1}$ and $h(L_{x},y)=H_{2}$ and no-flow Neumann conditions on top and bottom boundaries. In the saturated flow regime, the transport Eq. (22) is also linear and decoupled from the linear flow equation. Decoupled transport problems can be solved by either biased- or unbiased-GRW algorithms (see Remark 5 and Section 4.2.2) on the same lattice as that used to compute the flow velocity.

6.1 Flow in heterogeneous aquifers at regional scale

For the beginning, we follow the setup for regional scale used in [17] to compare approaches for inverse modeling of groundwater flow. The domain and the boundary conditions are specified by $L_{x}=4900$ m, $L_{y}=5000$ m, $H_{1}=0$ m, $H_{2}=5$ m. The hydraulic conductivity $K$ is a log-normally distributed random field defined by the mean $\langle K\rangle=12\cdot 10^{-4}$ m/s, the correlation length $\lambda=500$ m, and the variance $\sigma^{2}=1$ of the $\ln K$ -field. The $K$ -field is generated, as in Section 5.2.3 above, by summing 100 random periodic modes with the Kraichnan algorithm. Besides the exponential correlation considered in [17], we also investigate the behavior of the flow solution for Gaussian correlation of the $\ln K$ field with the same correlation length, as well as in case of the smaller variance $\sigma^{2}=0.1$ , for both correlation models.

The two correlation models of the $\ln K$ -field are of the form $C(r)=\sigma^{2}\exp[-(r/\lambda)^{\alpha}]$ , where $r=(r_{x}^{2}+r_{y}^{2})^{1/2}$ is the spatial lag, the exponent $\alpha=1$ corresponds to the exponential model, and $\alpha=2$ to the Gaussian one. Since the correlation functions depend on spatial variables through $r/\lambda$ , the computation can be done for spatial dimensions scaled by $\lambda$ , that is, fields of dimensionless correlation length $\lambda^{*}=1$ and a domain $[0,L_{x}/\lambda]\times[0,L_{y}/\lambda]$ . The results on the original grid are finally obtained after the multiplication by $\lambda$ of the solution $h(x,y)$ and of the spatial coordinates.

The solutions $h(x,y)$ of the stationary equation (17) corresponding to $\theta=const$ , for given realizations of the $K$ -field with $\sigma^{2}=0.1$ , are obtained under the initial condition $h_{0}(x,y)$ , which is the plane defined by the Dirichlet boundary conditions $h(0,y)=0$ and $h(L_{x}/\lambda,y)=H_{2}/\lambda$ . With space steps set to $\Delta x=\Delta y=0.2$ m, the steady state is reached after about $4\cdot 10^{5}$ iterations of the GRW solver. The relative errors of the solution $h$ obtained with the scaled geometry with respect to the solution of the unscaled problem are of the order $10^{-14}$ , that is, close to the machine precision [42, Sect. 6.1].

To estimate the order of convergence of the GRW scheme for this particular flow problem, we use manufactured analytical solutions provided in the Git repository https://github.com/PMFlow/FlowBenchmark and, similarly to estimations performed in Section 5.1, we compute the EOC according to (36) by successively halving the space steps from $\Delta x=\Delta y=2\cdot 10^{-1}$ up to $\Delta x=\Delta y=2.5\cdot 10^{-2}$ .

Table 11: Computational order of convergence of the GRW scheme estimated according to (36).

Correlation model	$\sigma^{2}$	$\varepsilon_{1}$	EOC	$\varepsilon_{2}$	EOC	$\varepsilon_{3}$	EOC	$\varepsilon_{4}$
Exponential	$0.1$	1.35e+01	3.67	1.06e+00	1.86	2.92e-01	0.66	1.85e-01
Exponential	$1$	1.80e+02	3.24	1.90e+01	2.09	4.47e+00	1.96	1.15e+00
Gaussian	$0.1$	7.37e-02	1.98	1.87e-02	1.63	6.03e-03	1.14	2.73e-03
Gaussian	$1$	1.31e-01	1.59	4.35e-02	1.51	1.53e-02	1.47	5.51e-03

Table 12: Computational order of convergence of the TPFA solver estimated according to (36).

Correlation model	$\sigma^{2}$	$\varepsilon_{1}$	EOC	$\varepsilon_{2}$	EOC	$\varepsilon_{3}$	EOC	$\varepsilon_{4}$
Exponential	$0.1$	4.67e+00	1.71	1.43e+00	1.95	3.70e-01	0.48	2.65e-01
Exponential	$1$	1.01e+02	2.23	2.14e+01	3.11	2.48e+00	0.41	1.86e+00
Gaussian	$0.1$	9.22e-02	2.00	2.30e-02	2.00	5.75e-03	2.00	1.44e-03
Gaussian	$1$	1.84e-01	2.00	4.61e-02	2.00	1.16e-02	2.00	2.89e-03

We note that the EOC approach presented here differers somewhat from that used in [2, 41]. The reference solution is now the manufactured solution, instead of the solution on the finest grid, and the error norm is no longer computed after the first iteration but after large numbers of iterations (from $10^{5}$ to more than $10^{7}$ ), when the GRW solution approaches the stationarity. Due to the limited number of iterations, the solutions are not yet strictly stationary and the order of convergence may be not accurately estimated in some cases. Therefore we also use a TPFA flow solver to compute EOC values for the same scenarios.

The results presented in Tables 11 and 12 show significant differences between the two correlation models. For Gaussian correlation the errors obtained with the two approaches are relatively small in all cases. Instead, for exponential correlation, despite the strong EOC obtained after the first two refinements, the errors are extremely large for $\sigma^{2}=1$ and become smaller than one only for $\sigma^{2}=0.1$ , after the second refinement of the grid. These results are consistent with those presented in [2], where similar benchmark problems were solved for a larger range of parameters of the $\ln K$ field.

6.2 Flow in conditions of random recharge

We consider in the following a flow problem formulated for the same geometry and boundary conditions as in the previous subsection, which has been used in [25] to design a new Monte Carlo approach for flow driven by spatially distributed stochastic sources. Now the hydraulic conductivity is constant, $K=12\cdot 10^{-4}$ m/s, and the groundwater recharge is described by a source term $f$ in Eq. (17), modeled as a random space function of mean $\langle f\rangle=362.912$ mm/year, log-normally distributed with exponential correlation specified by different correlation lengths and variances of the $\ln f$ field. Among different scenarios presented in[25], we consider for comparison with the present computations only the case $\lambda=500$ m and the variance $\sigma^{2}=1$ .

As in the previous subsection, we use the setup for the problem’s geometry scaled by $\lambda$ , for which the random recharge problem with $\sigma^{2}=1$ is solved with relative errors of the order $10^{-15}$ [42, Sect. 6.2].

In a first validation test, we compare the GRW and TPFA solutions of the random recharge problem on the computational domain scaled by $\lambda=500$ m, for single-realizations of the random recharge with both exponential and Gaussian correlation of the $\ln f$ field and two variances, $\sigma^{2}=0.1$ and $\sigma^{2}=1$ . The absolute and relative differences, $\varepsilon_{a}=\|h^{GRW}-h^{TPFA}\|$ and $\varepsilon_{r}=\|h^{GRW}-h^{TPFA}\|/\|h^{TPFA}\|$ , presented in Table 13 indicate a good agreement between the two approaches.

Table 13: Comparison of GRW and TPFA solutions of the
random recharge problem.

Correlation model	$\sigma^{2}$	$\varepsilon_{a}$	$\varepsilon_{r}$
Exponential	0.1	63.44	5.97e-2
Exponential	1	101.71	9.82e-2
Gaussian	0.1	84.12	8.72e-2
Gaussian	1	137.09	1.62e-2

Further, we perform statistical inferences of the mean and variance obtained from an ensemble of 100 Monte Carlo simulations within the setup of [25] for random recharge term with exponential correlation and variance $\sigma^{2}=1$ . The mean and the variance of the hydraulic head $h$ are computed as averages over realizations of the $\ln f$ field followed by spatial averages, with standard deviation estimated by spatial averaging. The results presented in Table 14 show, again, that the GRW and TPFA results are in good statistical agreement.

Table 14: Statistical moments of the hydraulic head
(Monte Carlo and spatial averages).

	mean	variance
GRW	21.51 $\pm$ 9.17	41.11 $\pm$ 27.82
TPFA	19.74 $\pm$ 7.84	32.09 $\pm$ 21.33

Finally, we compare the mean and the variance estimated at the center of the computational domain by GRW and TPFA simulations with the results presented in [25]. As seen in Table 15, the mean values compare quite well but both the GRW and TPFA approaches overestimate the variance computed for the same parameters in [25, Fig. 6]. This discrepancy can be attributed either to the large errors expected for exponential correlation model (see Tables 11 and 12) or to the statistical inhomogeneity of the Monte Carlo ensemble of 100 realizations indicated by the large standard deviations shown in Table 14.

Table 15: Statistical moments of the hydraulic head
(MC averages at the center of the domain).

	mean	variance
GRW	31.67	65.14
TPFA	28.31	53.39
(Passeto et al., 2011)	31.05	40.08

6.3 Flow and advection-dispersion transport in aquifers

In the following we consider an incompressible flow in the domain $[0,20]\times[0,10]$ , driven by Dirichlet boundary conditions $h(0,y)=1$ and $h(20,y)=0$ and zero Neumann conditions on top and bottom boundaries. The hydraulic conductivity is a random space function with mean $\langle K\rangle=15$ m/day, with Gaussian correlation of the $\ln K$ field, correlation length $\lambda=1$ m, and variance $\sigma^{2}=0.1$ , generated by summing 10 random modes with the Kraichnan algorithm. An ensemble of velocity fields corresponding to 100 realizations of the $K$ field is obtained with the flow solver used in Section 6.1, for the resolution of the GRW lattice defined by space steps $\Delta x=\Delta y=0.1$ .

Further, Monte Carlo simulations of advection-diffusion are carried out using the velocity realizations and the isotropic local dispersion coefficient $D=0.01$ m²/day. The linear transport equation obtained by setting $\theta=const$ in Eq. (22) is solved with the unbiased GRW algorithm described in Section 4.2.2 by using $\mathcal{N}=10^{24}$ particles to represent the concentration. The final time $T=10$ days is chosen such that the support of the concentration does not reach the boundaries during the simulation. Hence, the Monte Carlo inferences can be compared with results of linear theory which provides first-order approximations of dispersion coefficients for small variances $\sigma^{2}$ [6]. In turn, such linear approximations are accurately retrieved by averaging over ensembles of particle tracking simulations of diffusion in realizations of velocity fields approximated to the first-order in $\sigma^{2}$ by a Kraichnan procedure [35]. Following this approach, to infer dispersion coefficients in linear approximation, we use an ensemble of $10^{4}$ realizations of Krainchan velocity fields, computed with 100 random modes by the algorithm described in [40, Appendix C.3.2.2], and the unbiased GRW solver, with $\mathcal{N}=10^{24}$ particles in each realization. Longitudinal and transverse “ensemble” dispersion coefficients, $D_{x}$ and $D_{y}$ , are computed as half the slope of the ensemble average of the second spatial moments of the concentration distribution, centered at the ensemble average center of mass [6, 29, 35]. The results presented in Fig. 28 show a that, in spite of relatively small ensemble of velocity realizations, the ensemble dispersion coefficients obtained with the 100 GRW solutions of the full flow problem are quite close to the reference linear results.

The computation of the velocity realizations with the transient GRW flow solver requires $10^{4}$ to $10^{5}$ iterations to fulfill the convergence criterion (3) with tolerances $\varepsilon_{a}=\varepsilon_{r}=5\cdot 10^{-7}$ and about 160 seconds per realization. For the chosen discretization, $\Delta x=\Delta y=0.1$ , the unbiased GRW transport solver requires, according to (31), a relatively rough time discretization of $\Delta t=0.5$ . This leads to a total computation time of about 1.4 seconds for the estimation of the dispersion coefficients by averaging over the 100 realizations of the statistical ensemble. By comparison, the TPFA codes needs about 3.8 seconds to compute a velocity realization and about 13 seconds for a single transport realization, by using the same spatial resolution and a time step $\Delta t=0.05$ . But the TPFA estimates of the dispersion coefficients deviate by more than one order of magnitude from the linear reference solution. Since reducing the spatial steps and the local Pé to reduce the numerical diffusion dramatically increases the computational burden for the TPFA codes, we solved a rescaled problem. So, to preserve the mean and the spatial variability of the velocity field, we chose a smaller domain $[0,2]\times[0.1]$ , correlation length of the $\ln K$ field $\lambda=0.1$ , and a new Dirichlet condition, $h(0,y)=0.1$ . Now, the TPFA codes require about 60 seconds to compute one flow realization and about 3 hours for a transport realization, with $\Delta x=\Delta y=0.001$ and $\Delta t=0.0005$ . The computation times for the GRW codes to solve the rescaled problem by using $\Delta x=\Delta y=0.01$ and $\Delta t=0.07$ are practically unchanged. Figure 28 shows that the GRW estimations of the dispersion coefficients are again close to the linear approximation. Instead the TPFA results overestimate the linear approximation by 10% to 20%. The deviations of the TPFA coefficients shown in Fig. 28 are comparable with the numerical diffusion (estimated for constant velocity) in case of the longitudinal coefficient $D_{x}$ but two orders of magnitude larger for the transverse coefficient $D_{y}$ [42, Table 17].

7 Conclusions

The GRW schemes for simulating flow in either unsaturated or saturated porous media are equivalent to finite-difference schemes, in their deterministic implementation, or for sufficiently large numbers of particles in randomized implementations. The same, in case of BGRW solver for transport problems. Instead, the unbiased GRW is a superposition of Euler schemes for Itô equation [40], which is no longer equivalent with a finite difference scheme, unless the coefficients of the transport equation are constant. In simulations of reactive transport, GRW algorithms can use huge numbers of computational particles, even as large as the number of molecules involved in reactions, allowing simple and intuitive representations of the process.

While unbiased GRW algorithms are mainly efficient in obtaining fast solutions for large-scale transport in aquifers, BGRW solvers are appropriate for computing solutions of fully coupled flow and transport problems in soil systems with fine variation of the parameters. The algorithms are implemented as iterative $L$ -schemes which linearize the Richards equation and describe the transition from unsaturated to saturated regime. The GRW/BGRW solutions are first-order accurate in time and second-order accurate in space. For saturated regimes, the flow solver becomes a transient scheme solving steady-state flows in aquifers.

Since the GRW algorithms are explicit schemes which do not need to solve systems of algebraic equations, they are simpler and, in some cases, faster than finite element/volume schemes. The GRW $L$ -schemes for non-steady coupled problems for flow and transport in soils, as well as for transport simulations in saturated aquifers, are indeed much faster than the TPFA codes used as reference in this study. However, the flow solutions for saturated porous media in large domains (e.g. field or regional scale) require much larger computing time than classical numerical schemes, due to the large number of iterations needed to achieve the convergence of the transitory scheme used to compute steady-state solutions (see also a detailed analysis in [2]).

The obvious advantage of the GRW schemes is that they are practically free of numerical diffusion. This is demonstrated by the results for decoupled transport presented in Table 9. But, as shown by the discussion at the end of Section 5.2.3, the flow solvers also can be affected by numerical diffusion, which is difficult to isolate from other errors occurring in coupled flow and transport problems. Such errors are avoided by GRW algorithms, which prevent the occurrence of the numerical diffusion by using consistent definitions of the jump probabilities as functions of the coefficients of the flow and transport equations.

Acknowledgements

The authors are grateful to Dr. Emil Cătinaş for fruitful discussions on convergent sequences and successive approximation approaches. Nicolae Suciu acknowledges the financial support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Grant SU 415/4-1 – 405338726 “Integrated global random walk model for reactive transport in groundwater adapted to measurement spatio-temporal scales”. VISTA, a collaboration between the Norwegian Academy of Science and Letters and Equinor, funded the research of Davide Illiano, project number 6367, project name: “Adaptive model and solver simulation of enhanced oil recovery”.

References

[1] Abeele, W., 1984. Hydraulic Testing of Crushed Bandelier Tuff. Technical Report LA 10037-MS. Los Alamos National Laboratory, Los Alamos, New Mexico. https://www.osti.gov/biblio/60528-hydraulic-testing-crushed-bandelier-tuff
[2] Alecsa, C.D., Boros, I., Frank, Knabner, P., Nechita, M., Prechtel, A., Rupp, A., Suciu, N., 2019. Numerical benchmark study for fow in heterogeneous aquifers. Adv. Water Resour., 138, 103558. https://doi.org/10.1016/j.advwatres.2020.103558
[3] Anderson, D.G., 1965. Iterative procedures for nonlinear integral equations. J. ACM, 12(4), 547–560. https://doi.org/10.1145/321296.321305
[4] Alt, W., Luckhaus, H., 1983. Quasilinear elliptic-parabolic differential equations, Math. Z., 183(3), 311–341. https://doi.org/10.1007/BF01176474
[5] Bause, M., Knabner, P., 2004. Numerical simulation of contaminant biodegradation by higher order methods and adaptive time stepping. Comput. Visual. Sci., 7(2), 61–78. https://doi.org/10.1007/s00791-004-0139-y
[6] Bellin, A., Salandin, P., Rinaldo, A., 1992. Simulation of dispersion in heterogeneous porous formations: Statistics, first-order theories, convergence of computations. Water Resour. Res. 28(9), 2211–2227. https://doi.org/10.1029/92WR00578
[7] Bellin, A., Fiori, A., Dagan, G., 2020. Equivalent and effective conductivities of heterogeneous aquifers for steady source flow, with illustration for hydraulic tomography. Adv. Water Resour., 142, 103632. https://doi.org/10.1016/j.advwatres.2020.103632
[8] Both, J.W., Kumar, K., Nordbotten, J.M., Radu, F.A., 2019. Anderson accelerated fixed-stress splitting schemes for consolidation of unsaturated porous media. Comput. Math. Appl., 77(6), 1479–1502. https://doi.org/10.1016/j.camwa.2018.07.033
[9] Cătinaş, E., 2019. A survey on the high convergence orders and computational convergence orders of sequences. Appl. Math. Comput., 343, 1–20. https://doi.org/10.1016/j.amc.2018.08.006
[10] Catinas, E., 2020. How many steps still left to x*?, SIAM Rev., to appear.
[11] Caviedes-Voullième, D., Garci, P., Murillo, J., 2013. Verification, conservation, stability and efficiency of a finite volume method for the 1D Richards equation. J. Hydrol., 480, 69–84. https://doi.org/10.1016/j.jhydrol.2012.12.008
[12] Fujita, H., 1952. The exact pattern of a concentration-dependent diffusion in a semi-infinite medium, part II. Textil Res. J., 22(12), 823–827. https://doi.org/10.1177/004051755202201209
[13] Gardner, W.R., 1958. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci., 85(4), 228–232. https://journals.lww.com/soilsci/toc/1958/04000
[14] Gotovac, H., Cvetković, V., Andričevicć, R., 2009. Adaptive Fup multi-resolution approach to flow and advective transport in highly heterogeneous porous media: Methodology, accuracy and convergence. Adv. Water Resour. 32(6), 885–905. https://doi.org/10.1016/j.advwatres.2009.02.013
[15] Hajibeygi, H., Olivares, M.B., HosseiniMehr, M., Pop, S., Wheeler, M., 2020. A benchmark study of the multiscale and homogenization methods for fully implicit multiphase flow simulations. Adv. Water Resour., 143, 103674. https://doi.org/10.1016/j.advwatres.2020.103674
[16] Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P. J., Vachaud, G., 1977. A comparison of numerical simulation models for one-dimensional infiltration 1. Soil. Sci. Soc. Am. J., 41(2), 285–294. https://doi.org/10.2136/sssaj1977.03615995004100020024x
[17] Hendricks Franssen, H.J., Alcolea, A., Riva, M., Bakr, M., Van der Wiel, N., Stauffer, F., Guadagnini, A., 2009. A comparison of seven methods for the inverse modelling of groundwater flow. Application to the characterisation of well catchments. Adv. Water Resour., 32(6), 851–872. https://doi.org/10.1016/j.advwatres.2009.02.011
[18] Illiano, D., Pop, I.S., Radu, F.A., 2020. Iterative schemes for surfactant transport in porous media. Comput. Geosci. https://doi.org/10.1007/s10596-020-09949-2
[19] Knabner, P., Angermann, L., 2003. Numerical Methods for Elliptic and Parabolic Partial Differential Equations. Springer, New York.
[20] Knabner, P., Bitterlich, S., Iza Teran, R., Prechtel, A., Schneid, E., 2003. Influence of Surfactants on Spreading of Contaminants and Soil Remediation. In: Jäger W., Krebs H.J. (eds) Mathematics – Key Technology for the Future. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-55753-8_12
[21] Kuzmin, D., 2009. Explicit and implicit FEM-FCT algorithms with flux linearization. J. Comput. Phys., 228(7):2517–2534. http://dx.doi.org/doi:10.1016/j.jcp.2008.12.011.
[22] Lie, K.-A., 2019. An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the MATLAB Reservoir Simulation Toolbox (MRST). Cambridge University Press. https://doi.org/10.1017/9781108591416
[23] Liu, F., Fukumoto, Y., Zhao, X., 2020. Stability Analysis of the Explicit Difference Scheme for Richards Equation. Entropy, 22(3), 352. https://doi.org/10.3390/e22030352
[24] List, F., Radu, F.A., 2016. A study on iterative methods for solving Richards’ equation. Comput. Geosci., 20(2), 341–353. https://doi.org/10.1007/s10596-016-9566-3
[25] Pasetto, D., Guadagnini, A., Putti, M., 2011. POD-based Monte Carlo approach for the solution of regional scale groundwater flow driven by randomly distributed recharge. Adv. Water Resour., 34(11), 1450–1463. https://doi.org/10.1016/j.advwatres.2011.07.003
[26] Philip, J.R., 1969. Theory of infiltration. Advances in hydroscience, 5, 215–296. https://doi.org/10.1016/B978-1-4831-9936-8.50010-6
[27] Phoon, K.K., Tan, T.S., Chong, P.C., 2007. Numerical simulation of Richards equation in partially saturated porous media: under-relaxation and mass balance. Geotech. Geol. Eng., 25(5), 525-541. https://doi.org/10.1007/s10706-007-9126-7
[28] Pop, I.S., Radu, F.A., Knabner, P., 2004. Mixed finite elements for the Richards’ equation: linearization procedure, J. Comput. Appl. Math., 168(1), 365–373. https://doi.org/10.1016/j.cam.2003.04.008
[29] Radu, F.A., Suciu, N., Hoffmann, J., Vogel, A., Kolditz, O., Park, C.-H., Attinger, S., 2011. Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study. Adv. Water Resour., 34, 47–61. http://dx.doi.org/10.1016/j.advwatres.2010.09.012.
[30] Radu, F.A., Wang, W., 2014. Convergence analysis for a mixed finite element scheme for flow in strictly unsaturated porous media. Nonlinear Anal. R. World Appl., 15, 266–275. https://doi.org/10.1016/j.nonrwa.2011.05.003
[31] Radu, F.A., Kumar, K. , Nordbotten, J.M., Pop, I.S., 2018. A robust, mass conservative scheme for two-phase flow in porous media including hölder continuous nonlinearities. IMA J. Numer. Anal., 38(2), 884–920. https://doi.org/10.1093/imanum/drx032
[32] Sander, G.C., Parlange, J.Y., Kühnel, V., Hogarth, W.L., Lockington, D., O’kane, J.P.J., 1988. Exact nonlinear solution for constant flux infiltration. J. Hydrol., 97(3–4), 341–346. https://doi.org/10.1016/0022-1694(88)90123-0
[33] Schneid, E., Prechtel, A., Knabner, P., 2000. A comprehensive tool for the simulation of complex reactive transport and flow in soils. Land Contam. Reclamat., 8, 357–365. https://doi.org/10.2462/09670513.570
[34] Schneid, E., 2000. Hybrid-gemischte finite-elemente-diskretisierung der Richards-Gleichung. Doctoral dissertation, Naturwissenschaftliche Fakultät der Friedrich-Alexander-Universität Erlangen-Nürnberg.
[35] Schwarze, H., Jaekel, U., Vereecken, H., 2001. Estimation of macrodispersion by different approximation methods for flow and transport in randomly heterogeneous media. Transport Porous Media, 43(2), 265-287. https://doi.org/10.1023/A:1010771123844
[36] Simunek, J., Sejna, M., Saito, H., Sakai, M., van Genuchten, M., 2008. The Hydrus-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media, Version 4.0. Department of Environmental Sciences, University of California Riverside. https://www.pc-progress.com/Downloads/Pgm_hydrus1D/HYDRUS1D-4.08.pdf
[37] Slodicka, M., 2002. A robust and efficient linearization scheme for doubly non-linear and degenerate parabolic problems arising in flow in porous media. SIAM J. Numer. Anal., 23 (5), 1593–1614. https://doi.org/10.1137/S1064827500381860
[38] Srivastava, R., Yeh, T.C.J., 1991. Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils. Water Resour. Res., 27(5), 753–762. https://doi.org/10.1029/90WR02772
[39] Strikwerda, J.C., 2004. Finite Difference Schemes and Partial Differential Equations. SIAM. https://doi.org/10.1137/1.9780898717938
[40] Suciu, N., 2019. Diffusion in Random Fields. Applications to Transport in Groundwater. Birkhäuser, Cham. https://doi.org/10.1007/978-3-030-15081-5
[41] Suciu, N., 2020. Global Random Walk Solutions for Flow and Transport in Porous Media, in F.J. Vermolen, C. Vuik (eds.), Numerical Mathematics and Advanced Applications ENUMATH 2019, Lecture Notes in Computational Science and Engineering 139. Springer Nature, Switzerland. https://doi.org/10.1007/978-3-030-55874-1_93
[42] Suciu, N., Illiano, D., Prechtel, A., Radu, F.A., 2020. Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media (extended version). arXiv:2011.12889
[43] Suciu, N., Illiano, D., Prechtel, A., Radu, F.A., 2021. RichardsEquation, Git repository. https://doi.org/10.5281/zenodo.4709693
[44] Vamoş, C., Suciu, N., Vereecken, H., 2003. Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J. Comput. Phys., 186, 527–544. https://doi.org/10.1016/S0021-9991(03)00073-1.
[45] Walker, H.F., Ni, P., (2011). Anderson acceleration for fixed-point iterations. SIAM J. Numer. Anal., 49(4), 1715–1735. https://doi.org/10.1137/10078356X
[46] Warrick, A.W., Lomen, D.O., Yates, S.R., 1985. A generalized solution to infiltration. Soil. Sci. Soc. Am. J., 49(1), 34–38. https://doi.org/10.2136/sssaj1985.03615995004900010006x
[47] Watson, K.K., Sardana, V.A., Sander, G.C., 1995. Comparison of analytical and numerical results for constant flux infiltration. J. Hydrol., 165(1–4), 101–112. https://doi.org/10.1016/0022-1694(94)02580-5
[48] Zadeh, K.S., 2011. A mass-conservative switching algorithm for modeling fluid flow in variably saturated porous media. J. Comput. Phys., 230(3), 664-679. https://doi.org/10.1016/j.jcp.2010.10.011
[49] Zambra, C.E., Dumbser, M., Toro, E.F., Moraga, N.O., 2012. A novel numerical method of high-order accuracy for flow in unsaturated porous media. Int. J. Numer. Meth. Engng., 89(2), 227–240. https://doi.org/10.1002/nme.3241

[1] W. Abeele, Hydraulic Testing of Crushed Bandelier Tuff, Technical Report LA 10037-MS, Los Alamos National Laboratory, Los Alamos, New Mexico (1984) Google Scholar

[2] C.D. Alecsa, I. Boros, P. Knabner Frank, M. Nechita, A. Prechtel, A. Rupp, N. Suciu, Numerical benchmark study for fow in heterogeneous aquifers, Adv. Water Resour., 138 (2019), p. 103558, 10.1016/j.advwatres.2020.103558

[3] W. Alt, H. Luckhaus, Quasilinear elliptic-parabolic differential equations, Math. Z., 183 (3) (1983), pp. 311-341, 10.1007/BF01176474

[4] D.G. Anderson, Iterative procedures for nonlinear integral equations, J. ACM, 12 (4) (1965), pp. 547-560, 10.1145/321296.321305

[5] M. Bause, P. Knabner, Numerical simulation of contaminant biodegradation by higher order methods and adaptive time stepping, Comput. Visual. Sci, 7 (2) (2004), pp. 61-78, 10.1007/s00791-004-0139-y

[6] A. Bellin, A. Fiori, G. Dagan, Equivalent and effective conductivities of heterogeneous aquifers for steady source flow, with illustration for hydraulic tomography, Adv. Water Resour., 142 (2020), p. 103632, 10.1016/j.advwatres.2020.103632

[7] A. Bellin, P. Salandin, A. Rinaldo, Simulation of dispersion in heterogeneous porous formations: statistics, first-order theories, convergence of computations, Water Resour. Res., 28 (9) (1992), pp. 2211-2227, 10.1029/92WR00578

[8] J.W. Both, K. Kumar, J.M. Nordbotten, F.A. Radu, Anderson accelerated fixed-stress splitting schemes for consolidation of unsaturated porous media, Comput. Math. Appl., 77 (6) (2019), pp. 1479-1502, 10.1016/j.camwa.2018.07.033

[9] D. Caviedes-Voullième, P. Garci, J. Murillo, Verification, conservation, stability and efficiency of a finite volume method for the 1D Richards equation, J. Hydrol., 480 (2013), pp. 69-84, 10.1016/j.jhydrol.2012.12.008

[10] E. Cătinaş, A survey on the high convergence orders and computational convergence orders of sequences, Appl. Math. Comput., 343 (2019), pp. 1-20, 10.1016/j.amc.2018.08.006

[11] E. Cătinaş, How many steps still left to x*? SIAM Rev (2020)

[12] H.J. Hendricks Franssen, A. Alcolea, M. Riva, M. Bakr, N. Van der.Wiel, F. Stauffer, A. Guadagnini, A comparison of seven methods for the inverse modelling of groundwater flow. Application to the characterisation of well catchments, Adv. Water Resour., 32 (6) (2009), pp. 851-872, 10.1016/j.advwatres.2009.02.011

[13] H. Fujita, The exact pattern of a concentration-dependent diffusion in a semi-infinite medium, Part II, Textil Res. J., 22 (12) (1952), pp. 823-827, 10.1177/004051755202201209

[14] W.R. Gardner, Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table, Soil Sci., 85 (4) (1958), pp. 228-232, CrossRef View Record in Scopus Google Scholar, https://journals.lww.com/soilsci/toc/1958/04000

[15] H. Gotovac, V. Cvetković, R. Andričevicć, Adaptive Fup multi-resolution approach to flow and advective transport in highly heterogeneous porous media: methodology, accuracy and convergence, Adv. Water Resour., 32 (6) (2009), pp. 885-905, 10.1016/j.advwatres.2009.02.013, Article Download PDF View Record in Scopus Google Scholar

[16] H. Hajibeygi, M.B. Olivares, M. Hosseini Mehr, S. Pop, M. Wheeler, A benchmark study of the multiscale and homogenization methods for fully implicit multiphase flow simulations, Adv. Water Resour., 143 (2020), p. 103674, 10.1016/j.advwatres.2020.103674, Article Download PDF View Record in Scopus Google Scholar

[17] Haverkamp, Vauclin, Touma, Wierenga, Vachaud, 1977, R. Haverkamp, M. Vauclin, J. Touma, P.J. Wierenga, G. Vachaud, A comparison of numerical simulation models for one-dimensional infiltration 1, Soil. Sci. Soc. Am. J., 41 (2) (1977), pp. 285-294, 10.2136/sssaj1977.03615995004100020024x, CrossRef View Record in Scopus Google Scholar

[18] D. Illiano, I.S. Pop, F.A. Radu, Iterative schemes for surfactant transport in porous media, Comput. Geosci. (2020), 10.1007/s10596-020-09949-2, Google Scholar

[18] P. Knabner, L. Angermann, Numerical Methods for Elliptic and Parabolic Partial Differential Equations, Springer, New York (2003), Google Scholar

[19] P. Knabner, S. Bitterlich, R. Iza Teran, A. Prechtel, E. Schneid, Influence of surfactants on spreading of contaminants and soil remediation Jäger W., Krebs H.J. (Eds.), Mathematics – Key Technology for the Future, Springer, Berlin, Heidelberg (2003), 10.1007/978-3-642-55753-8_12, Google Scholar

[20] D. Kuzmin, Explicit and implicit FEM-FCT algorithms with flux linearization, J. Comput. Phys., 228 (7) (2009), pp. 2517-2534, 10.1016/j.jcp.2008.12.011, Article Download PDF View Record in Scopus Google Scholar

[21] K.A. Lie, An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the MATLAB Reservoir Simulation Toolbox (MRST), Cambridge University Press (2019), 10.1017/9781108591416, Google Scholar

[21] F. List, F.A. Radu, A study on iterative methods for solving Richards’ equation, Comput. Geosci., 20 (2) (2016), pp. 341-353, 10.1007/s10596-016-9566-3, CrossRef View Record in Scopus Google Scholar

[22] F. Liu, Y. Fukumoto, X. Zhao, Stability analysis of the explicit difference scheme for Richards equation, Entropy, 22 (3) (2020), p. 352, 10.3390/e22030352, CrossRef View Record in Scopus Google Scholar

[23] D. Pasetto, A. Guadagnini, M. Putti, POD-based monte carlo approach for the solution of regional scale groundwater flow driven by randomly distributed recharge, Adv. Water Resour., 34 (11) (2011), pp. 1450-1463, 10.1016/j.advwatres.2011.07.003, Article Download PDF View Record in Scopus Google Scholar

[24] J.R. Philip, Theory of infiltration, Adv. Hydrosci., 5 (1969), pp. 215-296, 10.1016/B978-1-4831-9936-8.50010-6, Article Download PDF View Record in Scopus Google Scholar

[25] K.K. Phoon, T.S. Tan, P.C. Chong, Numerical simulation of richards equation in partially saturated porous media: under-relaxation and mass balance, Geotech. Geol. Eng., 25 (5) (2007), pp. 525-541, 10.1007/s10706-007-9126-7, CrossRef View Record in Scopus Google Scholar

[26] I.S. Pop, F.A. Radu, P. Knabner, Mixed finite elements for the Richards’ equation: linearization procedure, J. Comput. Appl. Math., 168 (1) (2004), pp. 365-373, 10.1016/j.cam.2003.04.008, Article Download PDF View Record in Scopus Google Scholar

[27] F.A. Radu, K. Kumar, J.M. Nordbotten, I.S. Pop, A robust, mass conservative scheme for two-phase flow in porous media including hölder continuous nonlinearities, IMA J. Numer. Anal., 38 (2) (2018), pp. 884-920, 10.1093/imanum/drx032, CrossRef View Record in Scopus Google Scholar

[28] F.A. Radu, N. Suciu, J. Hoffmann, A. Vogel, O. Kolditz, C.-H. Park, S. Attinger, Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study, Adv. Water Resour., 34 (2011), pp. 47-61, 10.1016/j.advwatres.2010.09.012, Article Download PDF View Record in Scopus Google Scholar

[29] F.A. Radu, W. Wang, Convergence analysis for a mixed finite element scheme for flow in strictly unsaturated porous media, Nonlinear Anal. R. World Appl., 15 (2014), pp. 266-275, 10.1016/j.nonrwa.2011.05.003, Article Download PDF View Record in Scopus Google Scholar

[30] G.C. Sander, J.Y. Parlange, V. Kühnel, W.L. Hogarth, D. Lockington, J.P.J. O’kane, Exact nonlinear solution for constant flux infiltration,J. Hydrol., 97 (3–4) (1988), pp. 341-346, 10.1016/0022-1694(88)90123-0, Article Download PDF View Record in Scopus Google Scholar

[31] E. Schneid, Hybrid-gemischte finite-elemente-diskretisierung der richards-gleichung, Naturwissenschaftliche Fakultät der Friedrich-Alexander-Universität Erlangen-Nürnberg (2000), Doctoral dissertation, Google Scholar, Schneid

[32] E. Schneid, A. Prechtel, P. Knabner, A comprehensive tool for the simulation of complex reactive transport and flow in soils, Land Contam. Reclamat., 8 (2000), pp. 357-365, 10.2462/09670513.570, View Record in Scopus Google Scholar

[33] H. Schwarze, U. Jaekel, H. Vereecken, Estimation of macrodispersion by different approximation methods for flow and transport in randomly heterogeneous media,Transp. Porous Media, 43 (2) (2001), pp. 265-287, 10.1023/A:1010771123844, View Record in Scopus Google Scholar

[34] J. Simunek, M. Sejna, H. Saito, M. Sakai, M. van Genuchten, The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media, Version 4.0, (2008), Google Scholar

[35] M. Slodicka, A robust and efficient linearization scheme for doubly non-linear and degenerate parabolic problems arising in flow in porous media, SIAM J. Numer. Anal., 23 (5) (2002), pp. 1593-1614, 10.1137/S1064827500381860, View Record in Scopus Google Scholar

[36] R. Srivastava, T.C.J. Yeh, Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils, Water Resour. Res., 27 (5) (1991), pp. 753-762, 10.1029/90WR02772, View Record in Scopus Google Scholar

[37] J.C. Strikwerda, Finite difference schemes and partial differential equations, SIAM (2004), 10.1137/1.9780898717938, Google Scholar

[38] N. Suciu, Diffusion in Random Fields. Applications to Transport in Groundwater, Birkhäuser, Cham (2019), 10.1007/978-3-030-15081-5, Google Scholar

[39] N. Suciu, Global Random Walk Solutions for Flow and Transport in Porous Media, Springer Nature, Switzerland (2020), 10.1007/978-3-030-55874-1_93, Google Scholar, Numerical Mathematics and Advanced Applications ENUMATH 2019, Lecture Notes in Computational Science and Engineering 139

[40] Suciu, N., Illiano, D., Prechtel, A., Radu, F. A., 2020. Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media (extended version). arXivpreprint: 2011.12889., Google Scholar

[41] Suciu, N., Illiano, D., Prechtel, A., Radu, F. A.2021. https://github.com/PMFlow/RichardsEquation Git repository. doi:10.5281/zenodo.4709693., Google Scholar

[42] C. Vamoş, N. Suciu, H. Vereecken, Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comput. Phys., 186 (2003), pp. 527-544, 10.1016/S0021-9991(03)00073-1, Article Download PDF View Record in Scopus Google Scholar

[43 H.F. Walker, P. Ni, Anderson acceleration for fixed-point iterations, SIAM J. Numer. Anal., 49 (4) (2011), pp. 1715-1735, 10.1137/10078356X, CrossRef View Record in Scopus Google Scholar

[44] A.W. Warrick, D.O. Lomen, S.R. Yates, A generalized solution to infiltration, Soil. Sci. Soc. Am. J., 49 (1) (1985), pp. 34-38, 10.2136/sssaj1985.03615995004900010006x, CrossRef View Record in Scopus Google Scholar

[45] K.K. Watson, V.A. Sardana, G.C. Sander, Comparison of analytical and numerical results for constant flux infiltration, J. Hydrol., 165 (1–4) (1995), pp. 101-112, 10.1016/0022-1694(94)02580-5, Google Scholar

[46] K.S. Zadeh, A mass-conservative switching algorithm for modeling fluid flow in variably saturated porous media, J. Comput. Phys., 230 (3) (2011), pp. 664-679, 10.1016/j.jcp.2010.10.011, Google Scholar

[47] C.E. Zambra, M. Dumbser, E.F. Toro, N.O. Moraga, A novel numerical method of high-order accuracy for flow in unsaturated porous media, Int. J. Numer. Meth. Eng., 89 (2) (2012), pp. 227-240, Google Scholar, 10.1002/nme.3241

2021

Global Random Walk Solutions for Flow and Transport in Porous Media

October 22, 2021

Approximation properties of a family of integral type operators

May 14, 2021

Accurate spectral collocation solutions to 2nd-order Sturm–Liouville problems

March 2, 2021

	$\displaystyle\overline{\delta n_{l,m\|l,m,k}^{s}}=\left[1-\left(r_{l-1/2,m,k}^{% s}+r_{l+1/2,m,k}^{s}+r_{l,m-1/2,k}^{s}+r_{l,m+1/2,k}^{s}\right)\right]% \overline{n_{l,m,k}^{s}}$
	$\displaystyle\overline{\delta n_{l\mp 1,m\|l,m,k}^{s}}=r_{l\mp 1/2,m,k}^{s}% \overline{n_{l,m,k}^{s}}$
	$\displaystyle\overline{\delta n_{l,m\mp 1\|l,m,k}^{s}}=r_{l,m\mp 1/2,k}^{s}% \overline{n_{l,m,k}^{s}}.$		(20)

Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media

Abstract

Authors

Keywords

References

PDF

Cite this paper as:

About this paper

Journal

Publisher Name

DOI

Print ISSN

Online ISSN

Google Scholar Profile

Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media

Abstract

Keywords:

1 Introduction

2 One-dimensional GRW algorithm for unsaturated/saturated flow in soils

Remark 1

Remark 2

3 Validation of the one-dimensional GRW flow algorithm

3.1 Transition from unsaturated to saturated flow regime

3.2 Comparison with experiments and exact solutions for unsaturated flows

4 Two-dimensional GRW solutions

4.1 Two-dimensional GRW algorithm for flow in soils and aquifers

Remark 3

4.2 GRW algorithms for two-dimensional fully coupled flow and surfactant transport

4.2.1 Biased GRW algorithm for transport problems

Remark 4

Remark 5

4.2.2 Unbiased GRW algorithm for transport problems

5 Validation of the two-dimensional GRW algorithms

5.1 GRW flow solutions

5.2 GRW/BGRW solutions for fully coupled flow and transport problems

5.2.1 Code verification tests

5.2.2 Estimates of numerical diffusion

5.2.3 Fully coupled water flow and surfactant transport

6 Two-dimensional GRW solutions for groundwater flow and transport at regional and field scales

6.1 Flow in heterogeneous aquifers at regional scale

6.2 Flow in conditions of random recharge

6.3 Flow and advection-dispersion transport in aquifers

7 Conclusions

Acknowledgements

References

References

Related Posts