Posts by Nicolae Suciu

Abstract

Nonlinearities of coupled flow and transport problems for partially saturated porous media are solved with explicit iterative L-schemes. Their behavior is analyzed with the aid of the computational orders of convergence. This approach allows highlighting the influence of the truncation errors in the numerical schemes on the convergence of the iterations. Further, by using manufactured exact solutions, error-based orders of convergence of the iterative schemes are assessed and the convergence of the numerical solutions is demonstrated numerically through grid-convergence tests.

Authors

Nicolae Suciu
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj-Napoca, Romania

Florin A. Radu
Center for Modeling of Coupled Subsurface Dynamics, University of Bergen, Norway

Emil Cătinaş
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj-Napoca, Romania

Keywords

Richards’ equation; coupled flow and transport; finite differences; global random walk; iterative schemes; convergence order

Paper coordinates

N. Suciu, F.A. Radu, E. Cătinaş, Iterative schemes for coupled flow and transport in porous media – Convergence and truncation errors, J. Numer. Anal. Approx. Theory, 53 (2024) no. 1, pp. 158-183, https://doi.org/10.33993/jnaat531-1429

PDF

https://ictp.acad.ro/jnaat/journal/article/view/1429/1496

About this paper

Journal

Journal of Numerical Analysis and Approximation Theory

Publisher Name

Romanian Academy Publishing House

Editura Academiei Romane

DOI

https://doi.org/10.33993/jnaat531-1429

Print ISSN

2457-6794

Online ISSN

2501-059X

google scholar link

]1] C.D. Alecsa, I. Boros, F. Frank, P. Knabner, M. Nechita, A. Prechtel, A.RuppandN. Suciu,Numerical benchmark study for fow in heterogeneous aquifers,Adv. Water Resour.,138(2020), 103558.https://doi.org/10.1016/j.advwatres.2020.103558
[2] E. Catinas, A survey on the high convergence orders and computational convergenceorders of sequences, Appl. Math. Comput.,343(2019), pp. 1–20.

https://doi.org/10.1016/j.amc.2018.08.006
[3] E. Catinas, How many steps stil l left to x*?, SIAM Rev.,63(2021) no. 3, pp. 585–624.https://doi.org/10.1137/19M1244858

[4] D. Illiano, I.S. Pop and F.A. Radu,Iterative schemes for surfactant transport inporous media, Comput. Geosci.,25(2021), pp. 805–822.https://doi.org/10.1007/s10596-020-09949-2

[5] P. Knabner, S. Bitterlich, R. Iza Teran, R., A. Prechtel and E. Schneid, Influence of surfactants on spreading of contaminants and soil remediation, pp.152–161,in Jager, W., Krebs, H.J. (Eds.),Mathematics–Key Technology for the Future, Springer,Berlin, Heidelberg, 2003.https://doi.org/10.1007/978-3-642-55753-8_12

[6] P. Knabne rand L. Angermann, Numerical Methods for El liptic and Parabolic PartialDifferential Equations, Springer, New York, 2003.https://doi.org/10.1007/b97419

[7] F. List and F.A. Radu, A study on iterative methods for solving Richards’ equa-tion, Comput. Geosci.,20(2016) no. 2, pp. 341–353.https://doi.org/10.1007/s10596-016-9566-3

[8] I.S. Pop, F.A. Radu and P. Knabner,Mixed finite elements for the Richards’ equation:linearization procedure, J. Comput. Appl. Math.,168(2004) no. 1, pp. 365–373.https://doi.org/10.1016/j.cam.2003.04.008
[9] F.A. Radu, K. Kumar, J.M. Nordbotten and I.S. Pop, A robust, mass con-servative scheme for two-phase flow in porous media including Holder continuous nonlinearities, IMA Journal Numer. Anal.,38(2018) no. 2, pp. 884–920. url-https://doi.org/10.1093/imanum/drx032
[10] P.J. Roache,Code verification by the method of manufactured solutions, J. Fluids Eng.,124(2002) no. 1, pp. 4–10.http://dx.doi.org/10.1115/1.1436090
[11] C.J. Roy,Review of code and solution verification procedures for computational simu-lation, J. Comput. Phys.,205(2005), pp. 131–156.
https://doi.org/10.1016/j.jcp.2004.10.036
[12] J.S. Stokke, K. Mitra, E. Storvik, J.W. BothandF.A. Radu, An adaptivesolution strategy for Richards’ equation, Comput. Math. Appl.,152(2023), pp. 155–167.https://doi.org/10.1016/j.camwa.2023.10.020
[13] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, SIAM,2004.https://doi.org/10.1137/1.9780898717938
[14] N. Suciu,Diffusion in Random Fields. Applications to Transport in Groundwater,Birkhauser, Cham, 2019.https://doi.org/10.1007/978-3-030-15081-5
[15] N. Suciu, Global Random Walk Solutions for Flow and Transport in Porous Media,inF.J. VermolenandC. Vuik(eds.),Numerical Mathematics and Advanced Applica-tions ENUMATH 2019, Lecture Notes in Computational Science and Engineering,139,Springer Nature, Switzerland, 2020.https://doi.org/10.1007/978-3-030-55874-1_93
[16] N. Suciu, D. Illiano, A. Prechtel and F.A. Radu, Global random walk solversfor ful ly coupled flow and transport in saturated/unsaturated porous media, Adv. WaterResour.,152(2021), 103935.https://doi.org/10.1016/j.advwatres.2021.103935
[17] N. Suciu and F.A. Radu, Global random walk solvers for reactive transport andbiodegradation processes in heterogeneous porous media, Adv. Water Res.,166(2022),104268.https://doi.org/10.1016/j.advwatres.2022.104268
[18] N. Suciu, F.A. Radu and I.S. Pop, Space-time upscaling of reactive transport inporous media, Adv. Water Resour.,176(2023), 104443.http://dx.doi.org/10.1016/j.advwatres.2023.104443
[19] N. Suciu, F.A. Radu, J.S. Stokke, E. Catinas and A. Malina, Computa-tional orders of convergence of iterative methods for Richards’ equation, arXiv preprintarXiv:2402.00194 (2024),https://doi.org/10.48550/arXiv.2402.00194
[20] C. Vamos, N. Suciu and M Peculea, Numerical model ling of the one-dimensional diffusion by random walkers, Rev. Anal. Numer. Theor. Approx.,26(1997) nos. 1–2, pp. 237–247.https://ictp.acad.ro/jnaat/journal/article/view/1997-vol26-nos1-2-art32/
[21] C. Vamos, N. Suciu and H. Vereecken, Generalized random walk algorithm for thenumerical modeling of complex diffusion processes, J. Comput. Phys.,186(2003), pp.527–544.https://doi.org/10.1016/S0021-9991(03)00073-1Received by the editors: May 17, 2024; accepted: June 26, 2024; published online: July11, 2024.

Paper in HTML form

2024_1_Suciu_Radu_Catinas+(3)

ITERATIVE SCHEMES FOR COUPLED FLOW AND TRANSPORT IN POROUS MEDIA - CONVERGENCE AND TRUNCATION ERRORS

NICOLAE SUCIU $^{1}$ , FLORIN A. RADU $^{2}$ , and EMIL CĂTINAŞ $^{1}$

Abstract

Nonlinearities of coupled flow and transport problems for partially saturated porous media are solved with explicit iterative $L$ -schemes. Their behavior is analyzed with the aid of the computational orders of convergence. This approach allows highlighting the influence of the truncation errors in the numerical schemes on the convergence of the iterations. Further, by using manufactured exact solutions, error-based orders of convergence of the iterative schemes are assessed and the convergence of the numerical solutions is demonstrated numerically through grid-convergence tests.

MSC. 35G31, 76S05, 76R50, 65M06, 65C35.
Keywords. Richards' equation, coupled flow and transport, finite differences, global random walk, iterative schemes, convergence order.

1. INTRODUCTION

Flow and reactive transport in porous media are modeled by systems of coupled nonlinear partial differential equations. The seepage velocity through the bulk of porous medium is proportional to the hydraulic conductivity and to the gradient of the pressure head, according to Darcy's law, which combined with the conservation equation results in the Richards' equation. In its general, mixed formulation, Richards' equation contains two unknown functions, the volumetric water content and the pressure head. Therefore, it needs to be closed by constitutive relationships for the dependence of the water content on the pressure head and the dependence of the hydraulic conductivity on the water content. Such relationships are highly nonlinear, e.g., the exponential model or the van Genuchten-Mualem model, with discontinuous derivatives when the negative pressure, in unsaturated flow regime, changes into the positive pressure characterizing the saturated regime with constant water content and hydraulic conductivity. Consequently, the resulting Richards' equation for the pressure head is non-linear and parabolic-elliptic degenerate.

The concentration of the chemical compounds transported in soils and aquifers is governed by parabolic advection-diffusion-reaction equations, with flow velocity provided by the solution of the Richards' equation. In the presence of capillarity effects, for instance in problems of surfactant transport, the water content depends on pressure as well as on concentration and the flow and transport equations are coupled in a nonlinear way in both directions. Solving these coupled nonlinear problems requires iterative schemes that use different linearization techniques, such as Newton's method, (modified) Picard method,

L

-scheme or their combinations. While Newton's method is second order (under usual circumstances) but only locally convergent, which means that it converges only if the initial choice is not too far from the solution, the

L

-scheme is a robust method unconditionally first order convergent, the robustness of the Picard method being somewhere in between [7]. A significant improvement of the convergence can be achieved by starting with

L

-scheme iterations and, once an appropriate initial solution is found, switching to the faster and higher order convergent Newton's method [12].

The

L

-scheme can be thought as a quasi-Newton method, with derivative of the water content with respect of the pressure head replaced by a positive constant

L

. This makes it robust and independent of the initial choice, at the price of being only first order convergent. Given its simplicity, the

L

-scheme is particularly suitable in solving intricate and fully coupled nonlinear problems. Most of the implementations of the

L

-scheme are implicit linearizations based on finite element [8,7,12] or finite volume [9,4] discretizations. Explicit

L

-schemes based on finite difference (FD) and random walk methods, introduced in [16], have been successfully used to solve coupled problems of flow and biodegradation in partially saturated soils and saturated heterogeneous aquifers [17]. The reason of using explicit

L

-schemes is that numerical solutions of the flow equations are obtained on the same regular grid on which reactive transport is solved with random walk algorithms free of numerical diffusion [16].

Theoretical proofs of convergence of the linearization schemes essentially show that the errors of the iterative solutions with respect to the exact solution of the nonlinear problem can be made arbitrarily small [7]. However, even if the numerical analysis of the scheme provides estimates of the convergence rates, in general it does not produce a thorough and comprehensive characterization of the type of convergence behavior, by predicting the upper/lower limits of the rates as in the theory of the convergent sequences

[2, 3]

. The analysis can be completed either by assessing computational orders of convergence, with sequences of successive corrections, or by assessing error-based orders, when exact solutions of the problem are available [19]. This approach could be mainly useful when the constraints required by the theoretical proofs are hard to meet, or when there are no convergence proofs at all, but, in practice, the convergence of the scheme is demonstrated numerically.

We investigate in this framework the convergence of the explicit

L

-scheme through computational and error-based orders of convergence on three numerical examples: solutions of Richards' equation, fully coupled flow and surfactant transport, and coupled flow and transport with van Genuchten-Mualem parametrization.

2. COUPLED FLOW AND TRANSPORT IN POROUS MEDIA

We consider in the following the fully coupled one-dimensional problem of flow governed by Richards' equation and transport of a single chemical species governed by an advection-diffusion-reaction equation,

\begin{matrix} (1) & \frac{\partial}{\partial t} θ (ψ, c) - \frac{\partial}{\partial z} [K (θ (ψ, c)) \frac{\partial}{\partial z} (ψ + z)] = 0 \\ (2) & \frac{\partial}{\partial t} [θ (ψ, c) c] - \frac{\partial}{\partial z} [D \frac{\partial}{\partial z} c - q c] = R (c) \end{matrix}

where

ψ (z, t)

is the pressure head, with

z

being the height against the gravitational direction,

θ

is the water content,

q = - K (θ (ψ, c)) \frac{\partial}{\partial z} (ψ + z)

is the water flux (also called Darcy velocity),

K

stands for the hydraulic conductivity of the porous medium,

c

is the concentration,

D

is the diffusion/dispersion coefficient, and

R

is the reaction term. The system (1)-(2) is closed via constitutive relationships defining the dependencies

θ (ψ, c)

and

K (θ (ψ, c))

Equations (1)-(2) are fully coupled through the terms

θ (ψ, c)

and

[θ (ψ, c) c]

. Richards' equation (1) is parabolic-elliptic degenerate, with variable

θ (ψ, c)

ψ < 0

and

θ =

const if

ψ \geq 0

. Since the system is strongly nonlinear, through the functional dependencies

θ (ψ, c)

and

K (θ (ψ, c))

, linearization methods have to be used to compute numerical solutions.
2.1. Explicit

L

-schemes for one-dimensional Richards' equation. Following [16], we start with the staggered FD scheme with backward discretization in time of the Richards' equation (1) decoupled

^{1}

from the transport equation (2),

\begin{aligned} θ (ψ_{i, k}) - θ (ψ_{i, k - 1}) = \\ = \frac{Δ t}{Δ z^{2}} {[K (θ (ψ_{i + 1 / 2, k})) (ψ_{i + 1, k} - ψ_{i, k}) - K (θ (ψ_{i - 1 / 2, k})) (ψ_{i, k} - ψ_{i - 1, k})] \\ (3) & + (K (θ (ψ_{i + 1 / 2, k})) - K (θ (ψ_{i - 1 / 2, k}))) Δ z} \end{aligned}

To obtain an iterative method for solving Richard's equation, we denote by

ψ_{i, k}^{s}

the approximate solution after

s

iterations and add to the l.h.s. a stabilization term

L (ψ_{i, k}^{s + 1} - ψ_{i, k}^{s})

, where

L

is a constant which has the dimension of an inverse length. Whenever the iterative method converges, this stabilization term vanishes and the limit is the solution of the nonlinear scheme

(3). Defining the dimensionless numbers

r_{i \pm 1 / 2, k}^{s} = K (ψ_{i \pm 1 / 2, k}^{s}) Δ t / (L Δ z^{2})

, one obtains

\begin{matrix} ψ_{i, k}^{s + 1} = [1 - (r_{i + 1 / 2, k}^{s} + r_{i - 1 / 2, k}^{s})] ψ_{i, k}^{s} + r_{i + 1 / 2, k}^{s} ψ_{i + 1, k}^{s} + r_{i - 1 / 2, k}^{s} ψ_{i - 1, k}^{s} \\ (4) & + (r_{i + 1 / 2, k}^{s} - r_{i - 1 / 2, k}^{s}) Δ z - (θ (ψ_{i, k}^{s}) - θ (ψ_{i, k - 1})) / L \end{matrix}

In practice, we run successive iterations

s = 1, 2, \dots

until a threshold condition is fulfilled, e.g.,

‖ ψ_{k}^{s} - ψ_{k}^{s - 1} ‖ \leq ε_{a} + ε_{r} ‖ ψ_{k}^{s} ‖

, where

ψ_{k}^{s}

denotes the solution vector of components

ψ_{i, k}^{s}, i = 1, \dots, I, ‖ \cdot ‖

is the

l^{2}

norm, and

ε_{a} \geq 0

and

ε_{r} \geq 0

are the absolute and relative tolerances (see e.g., [7]). The relation (4) is an explicit

L

-scheme that allows the computation of the approximate solution recursively and avoids the need to solve systems of linear equations, as in implicit finite element [7] or finite volume [9]

L

-schemes. Consequently, the explicit scheme (4) vas found to be one order of magnitude faster than the implicit finite volume

L

-scheme when solving typical problems for partially saturated flows [16].

To transform the FD scheme (4) into a random walk scheme we represent the solution

ψ_{i,, k}^{s}

N

particles distributed over the lattice sites,

ψ_{i, k}^{s} \approx n_{i, k}^{s} a / N

, where

a

is a unit length that will be disregarded in the following. This results in the following relation summing contributions from neighboring sites to the number of particles at the site

i

and time step

k

n_{i, k}^{s + 1} = [1 - (r_{i + 1 / 2, k}^{s} + r_{i - 1 / 2, k}^{s})] 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} + ⌊ N f^{s} ⌋

, where

f^{s} = (r_{i + 1 / 2, k}^{s} - r_{i - 1 / 2, k}^{s}) Δ z - [θ (n_{i, k}^{s}) - θ (n_{i, k - 1})] / L

and

⌊ \cdot ⌋

is the floor function.

In order for (5) to be a random walk scheme, the coefficients multiplying numbers of particles at lattice sites

[1 - (r_{i + 1 / 2, k}^{s} + r_{i - 1 / 2, k}^{s})], r_{i + 1 / 2, k}

, and

r_{i - 1 / 2, k}

should be normalized probabilities. This implies the following restriction

[1 - (r_{i + 1 / 2, k}^{s} + r_{i - 1 / 2, k}^{s})] \leq 1

. A sufficient condition for that is

\begin{matrix} (6) & r_{i \pm 1 / 2, k}^{s} \leq 1 / 2 \end{matrix}

The first three terms in r.h.s. of (5) represent groups of particles jumping on the site

i

form the right and from the left, and the ratio of the initial number of particles at the site

i

which do not undergo jumps,

\begin{matrix} (7) & δ n_{i ∣ i \pm 1, k} = r_{i \pm 1 / 2, k}^{s} n_{i \pm 1, k}^{s}, δ n_{i ∣ i, k} = [1 - (r_{i + 1 / 2, k}^{s} + r_{i - 1 / 2, k}^{s})] n_{i, k}^{s} \end{matrix}

With (7), (5) becomes

\begin{matrix} (8) & n_{i, k}^{s + 1} = δ n_{i ∣ i, k} + δ n_{i ∣ i + 1, k} + δ n_{i ∣ i - 1, k} + ⌊ N f^{s} ⌋ \end{matrix}

The numbers

δ n

of random walkers are binomially distributed random variables. Unlike in classical random walk approaches, consisting of generating random walk trajectories and counting particles at lattice sites, we directly evaluate the binomial variables

δ n

and construct a global random walk (GRW)
[21]. Since there are different probabilities

r_{j \pm 1 / 2, k}

for the right/left jumps, one obtains a biased global random walk (BGRW) [14]. The resulting BGRW algorithm moves the particles from a given lattice site

j

according to the rule

\begin{matrix} (9) & n_{j, k}^{s} = δ n_{j ∣ j, k}^{s} + δ n_{j + 1 ∣ j, k}^{s} + δ n_{j - 1 ∣ j, k}^{s}, \end{matrix}

where

\begin{matrix} (10) & δ n_{j \pm 1 ∣ j, k} = r_{j \pm 1 / 2, k}^{s} n_{j, k}^{s}, δ n_{i ∣ i, k} = [1 - (r_{j + 1 / 2, k}^{s} + r_{j - 1 / 2, k}^{s})] n_{i, k}^{s} \end{matrix}

The random variable

δ n_{j + 1 ∣ j, k} = r_{j + 1 / 2, k}^{s} n_{j, k}^{s}

representing the number of right jumps in (9) is binomially distributed, with parameters

n_{j, k}^{s}, p = r_{j + 1 / 2, k}^{s} / r_{i, k}^{s}

, and

q = 1 - p = r_{j - 1 / 2, k}^{s} / r_{i, k}^{s}

, where

r_{i, k}^{s} = r_{j + 1 / 2, k}^{s} + r_{j - 1 / 2, k}^{s}

. The number of left jumps is the difference between the total number of jumps and the number of right jumps

δ n_{j - 1 ∣ j, k} = r_{i, k}^{s} n_{j, k}^{s} - δ n_{j + 1 ∣ j, k}

. Since the computation of the exact binomial distribution could be computationally prohibitive for very large

n_{j, k}^{s}

, one uses its approximation via a "reduced-fluctuations algorithm" similar to that proposed for the unbiased GRW algorithm in [21] as follows:

algorithm 1. Reduced fluctuations algorithm

(i) initialize a real variable rest

= 0

;
(ii) compute

δ n

with the relations (10);
(iii) split

δ n

into

δ n = ⌊ δ n ⌋ + {δ}

at every lattice site;
(iv) approximate

δ n

⌊ δ n ⌋

;
(v) sum up the remainders into the variable rest

= \sum {δ}

over all the sites;
(vi) allocate one particle to the first lattice site where rest

\geq 1

;
(vii) save rest

=

rest -1 and repeat the sequence (ii)-(vi) at every

k

Specifically, the algorithm is implemented in the Matlab code as follows:

Code 1. Flow-BGRW with reduced fluctuations

restr=0; rest1=0; rest2=0;
D=K(tht); D=D(1:I-1); % conductivity function of water content
r=dt*D/dx^2/L;
rloc=[1-2*r(1),1-(r(1:I-2)+r(2:I-1)),1-2*r(I-1)];
restr=rloc.*n+restr; nn=floor(restr); restr=restr-nn;
rest1=r(2:I-1).*n(3:I)+rest1;
njumpleft=floor(rest1); rest1=rest1-njumpleft; % left jumps
nn(2:I-1)=nn(2:I-1)+njumpleft;
rest2=r(1:I-2).*n(1:I-2)+rest2;
njumpright=floor(rest2); rest2=rest2-njumpright; % right jumps
nn(2:I-1)=nn(2:I-1)+njumpright;

Code 1 computes the solution in the interior of the domain. The Darcy velocity on the boundaries can be computed in several ways: using an approximate forward finite difference discretization of Darcy's law, extending on the boundary the velocity from the first neighboring interior site, computing
the exact velocity analytically when a manufactured solution for the pressure head is available [16, Sect. 5.2.1].

Remark 1. Since the BGRW is an explicit scheme, the source term can be evaluated first, then the first tree terms of (8) which model the transport will be computed according to Code 1. The reduced-fluctuations algorithm preserves the indivisibility of the particles by summing up to unity the fractional parts

{\cdot}

of all the terms defined by (7). Therefore, the floor function in (8) is no longer necessary and the source term becomes

N f^{s}

. Dividing (8) by

N

, using (7), and taking the average, one retrieves the

F D L

-scheme (4) for the unknown

ψ_{i, k}^{s} = \overset{―}{n_{i, k}^{s}} / N

. A condition of kind

K Δ t / Δ z^{2} \leq 1 / 2

ensures the stability of the forward-time central-space scheme for the heat equation with constant coefficient

K

and is assumed (without carrying out the stability analysis) to hold in case of variable coefficients as well [13]. In the present context, the condition (6) guaranties that particles cannot be created or destroyed otherwise than through boundary conditions or in the presence of an internal source/sink, that is, it ensures the stability of the BGRW L-scheme (5). Since the FD Lscheme (4) is the ensemble average of the BGRW scheme, its stability is also ensured by the constraint (6).

The FD scheme is also retrieved when all floor operators are removed in Code 1. The resulting code sequence can be recast as the FD scheme presented in Code 2 below, where

p = n / N

and

p p = n n / N

Code 2. Flow-FD scheme

D = K (

); D = D (1 : I - 1); %

conductivity function of water content

r = dt * D / {dx}^{\land} 2 / L

;

r l o c = [1 - 2 * r (1), 1 - (r (1 : I - 2) + r (2 : I - 1)), 1 - 2 * r (I - 1)]

;
pp=rloc.*p;

pp (2 : I - 1) = pp (2 : I - 1) + r (2 : I - 1) \cdot * p (3 : I) + r (1 : I - 2) \cdot * p (1 : I - 2);

Remark 2. In the particular case of the saturated regime,

θ =

const, with space-variable hydraulic conductivity

K

and a given source term

f

, after setting

L = 1

and disregarding the time index

k

in (4) one obtains a transient scheme to solve the equation for the hydraulic head

h = ψ + z

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

This scheme has been used to solve coupled flow and transport in saturated aquifers with moderate heterogeneity [15]. However, because the transient scheme may be very slow for highly heterogeneous coefficients

K

[1], a more efficient strategy to solve coupled problems for saturated porous media is the coupling between an implicit FD flow scheme [1] and a GRW transport solver, see https://github.com/PMFlow/Coupled_FDM_GRW.

Matlab codes for one- and two-dimensional solutions of Richards' equation, verification tests, and benchmark problems presented in [16] are posted at https://github.com/PMFlow/RichardsEquation.
2.2. Explicit BGRW scheme for reactive transport. Similarly to the derivation of the flow scheme in Section 2.1, we start with a FD scheme with backward time discretization of the equation (2) with a constant diffusion coefficient

D

θ (ψ_{i, k}, c_{i, k}) c_{i, k} - θ (ψ_{i, k - 1}, c_{i, k - 1}) c_{i, k - 1} = -

\begin{matrix} (11) & = - t frac Δ t 2 Δ z (q_{i + 1, k} c_{i + 1, k} - q_{i - 1, k} c_{i - 1, k}) + \frac{D Δ t}{Δ z^{2}} (c_{i + 1, k} - 2 c_{i, k} + c_{i - 1, k}) + R (c_{i, k}) = \end{matrix}

- \frac{2 D Δ t}{Δ z^{2}} c_{i, k} + (\frac{D Δ t}{Δ z^{2}} - \frac{Δ t}{2 Δ z} q_{i + 1, k}) c_{i + 1, k} + (\frac{D Δ t}{Δ z^{2}} + \frac{Δ t}{2 Δ z} q_{i - 1, k}) c_{i - 1, k} + Δ t R (c_{i, k})

,
where

ψ_{i, k}

is the solution of the coupled flow solver. Further, we denote by

c_{i, k}^{s}

the approximate solution after

s

iterations, add to the l.h.s. of (11) a stabilization factor

L (c_{i, k}^{s + 1} - c_{i, k}^{s})

, where

L

is now a dimensionless constant, and define the dimensionless parameters

\begin{matrix} (12) & r = \frac{2 D Δ t}{L Δ z^{2}}, v_{i \pm 1, k}^{s} = \frac{Δ t}{L Δ z} q_{i \pm 1, k}^{s} \end{matrix}

Representing the concentration by the distribution of

N

particles on the lattice,

c_{i, k}^{s} \approx n_{i, k}^{s} / N

, and using the parameters (12), the stabilized version of the scheme (11) becomes an explicit

L

-scheme,

\begin{matrix} (13) & n_{i, k}^{s + 1} = (1 - r) n_{i, k}^{s} + \frac{1}{2} (r - v_{i + 1, k}^{s}) n_{i + 1, k}^{s} + \frac{1}{2} (r + v_{i - 1, k}^{s}) n_{i - 1, k}^{s} + ⌊ g^{s} ⌋ \end{matrix}

where

g^{s} = N {Δ t R (n_{i, k}^{s}) / L - [θ (ψ_{i, k}^{s}, n_{i, k}^{s}) n_{i, k}^{s} - θ (ψ_{i, k - 1}, n_{i, k - 1}) n_{i, k - 1}] / L}

.
With the constrains

\begin{matrix} (14) & r \leq 1, | v_{i, k}^{s} | \leq r \end{matrix}

which ensure the normalization of the jump probabilities, the explicit

L

-scheme (13) for reactive transport is defined as a BGRW algorithm. The gathering and the spreading of groups of particles at/from a lattice point are described by relations similar to (8) and (9). The binomial random variables

δ n

are defined by

\begin{matrix} (15) & δ n_{j \pm 1 ∣ j, k}^{s} = \frac{1}{2} (r \pm v_{j, k}^{s}) n_{j, k}^{s}, δ n_{j ∣ j, k}^{s} = (1 - r) n_{j, k}^{s} \end{matrix}

Remark 3. The biased jump probabilities

\frac{1}{2} (r \pm v_{j, k}^{s})

in (15) simulate the advective displacement by moving larger numbers of particles in the positive flow direction. With an unbiased GRW, the advection is simulated by moving groups of particles over

v_{j, k}^{s}

lattice sites, i.e.,

δ n_{j + v ∣ j, k}^{s} = (1 - r) n_{j, k}^{s}

, where

v = v_{j, k}^{s} [16, 17]

. This procedure is faster but yields overshooting errors when there are important changes in velocity over a time step. The BGRW algorithm which allows only first-neighbor jumps is accurate but can be much slower than the unbiased GRW when a fine discretization is required to describe the spatial variability of the solution. Moreover, the second constraint (14) and the definitions of the parameters (12) imply an upper bound for the local Péclet number,

P e ́ = q_{i, k}^{s} Δ z / D \leq 2

The binomial random variables

δ n

defined by (15) are generated with the Algorithm 1. The Matlab implementation of the advection-diffusion steps of the transport solver is presented in Code 3 below.

Code 3. Transport-BGRW with reduced fluctuations

u=q*dt/Lc/dx; % computation of the BGRW parameters (12)
ru=2*D*dt/Lc/dx^2*ones(1,I);
restr=0; restjump=0;
for x=1:I
if n(x) > 0
    r=ru(x);
    restr=n(x)*(1-r)+restr; nsta=floor(restr);
    restr=restr-nsta; njump=n(x)-nsta;
    nn(x)=nn(x)+nsta;
    if(njump)>0
        restjump=njump*0.5*(1-u(x)/r)+restjump;
        nj(1)=floor(restjump); restjump=njump-nj(1);
        nj(2)=floor(restjump); restjump=restjump-nj(2);
        if x==1
            nn(2)=nn(2)+nj(2);
        elseif x==I
            nn(I-1)=nn(I-1)+nj(1);
        else
            for i=1:2
                xd=x+(2*i-3);
                nn(xd)=nn(xd)+nj(i);
            end
        end
    end
end
end

Remark 4. The ensemble average of the BGRW L-scheme for reactive transport also leads to an equivalent

F D

scheme. However, we do not pursue this approach which gives up the particle indivisibility. The numerical schemes developed in [16, 17] aim at providing a "microscopic" kinematic description of the reactive transport with piecewise-analytic time functions associated to each particle. Using a space-time averaging over this microscopic description we are then able to construct a macroscopic description through almost everywhere continuous fields that can be used to model experimental measurements or to achieve a spatio-temporal upscaling [18].

Matlab codes for coupled flow and transport and test problems are posted in the repository https://github.com/PMFlow/RichardsEquation. BGRW/GRW solvers for decoupled reactive transport with applications to biodegradation
processes in soils and aquifers, presented in [17], may be found in the repository https://github.com/PMFlow/MonodReactions.
2.3. Discussion on convergence. We look for convergence conditions by investigating the evolution of the difference

e_{ψ}^{s} = ψ^{s} - ψ

of the solutions of (4) and (3) and of the difference

e_{c}^{s} = c^{s} - c

of the solutions of (13) and (11) at fixed position

i

and time

k

Let

L_{ψ}

and

L_{c}

be the stabilization constants flow the flow and transport schemes, respectively. After denoting by

\frac{Δ t}{Δ z^{2}} D_{ψ} (ψ)

the r.h.s. of (3) and subtracting (3) from (4), we have

\begin{matrix} (16) & L_{ψ} (e_{ψ}^{s + 1} - e_{ψ}^{s}) + θ (ψ^{s}, c^{s}) - θ (ψ, c) = \frac{Δ t}{Δ z^{2}} (D_{ψ} (ψ^{s}) - D_{ψ} (ψ)) \end{matrix}

Similarly, with

\frac{Δ t}{2 Δ z} D_{q} (c)

and

\frac{Δ t}{Δ z^{2}} D_{D} (c)

denoting the advection and diffusion operators in the r.h.s. of (11), converting the numbers of particles to concentrations,

c_{i, k}^{s} = n_{i, k}^{s} / N

, and subtracting (11) from (13) we have

L_{c} (e_{c}^{s + 1} - e_{c}^{s}) + θ (ψ^{s}, c^{s}) c^{s} - θ (ψ, c) c = \frac{Δ t}{2 Δ z} (D_{q} (c^{s}) - D_{q} (c)) + \frac{Δ t}{Δ z^{2}} (D_{D} (c^{s}) - D_{D} (c))

.
Assuming

\begin{aligned} (A1) \partial θ / \partial ψ > 0, \partial θ / \partial c > 0, \partial (θ c) / \partial ψ > 0, \partial (θ c) / \partial c > 0 \\ (A2) | D_{ψ} (ψ^{s}) - D_{ψ} (ψ) | \leq L_{D_{ψ}} | e_{ψ}^{s} |, D_{D} (c^{s}) - D_{D} (c) | \leq L_{D_{D}} | e_{c}^{s} ∣ \end{aligned}

(16) and (17) give

\begin{matrix} (18) & {(e_{ψ}^{s + 1})}^{2} \leq {(e_{ψ}^{s})}^{2} {(1 - \frac{1}{L_{ψ}} \frac{\partial θ (ψ, c)}{\partial ψ} + \frac{L_{D_{ψ}} Δ t}{L_{ψ} Δ z^{2}})}^{2} + {(e_{c}^{s})}^{2} {(\frac{1}{L_{ψ}} \frac{\partial θ (ψ, c)}{\partial c})}^{2} \end{matrix}

\begin{matrix} (19) & {(e_{c}^{s + 1})}^{2} \leq {(e_{c}^{s})}^{2} {(1 - \frac{1}{L_{c}} \frac{\partial (θ (ψ, c) c)}{\partial c} + \frac{L_{D_{q}} Δ t}{2 L_{c} Δ z} + \frac{L_{D_{D}} Δ t}{L_{c} Δ z^{2}})}^{2} + {(e_{ψ}^{s})}^{2} {(\frac{1}{L_{c}} \frac{\partial (θ (ψ, c) c)}{\partial ψ})}^{2} . \end{matrix}

Let

F

denote one of the factors of

{(e_{ψ}^{s})}^{2}

and

{(e_{c}^{s})}^{2}

in (18) and (19). Assuming for all factors the condition

F \leq \frac{1 - ϵ}{2}

, where

ϵ \in (0, 1]

, one obtains

\begin{aligned} inf_{ψ} {\partial θ (ψ, c) / \partial ψ} & \geq L_{D_{ψ}} \frac{Δ t}{Δ z^{2}} + 1 - \sqrt{\frac{1 - ϵ}{2}} \\ sup_{ψ} {\partial θ (ψ, c) / \partial ψ} & \leq L_{D_{ψ}} \frac{Δ t}{Δ z^{2}} + 1 + \sqrt{\frac{1 - ϵ}{2}} \\ (20) & \frac{1}{L_{ψ}} \frac{\partial θ (ψ, c)}{\partial c} & \leq \sqrt{\frac{1 - ϵ}{2}} \end{aligned}

\begin{aligned} inf_{c} {\partial (θ (ψ, c) c) / \partial c} & \geq L_{D_{q}} \frac{Δ t}{2 Δ z} + L_{D_{D}} \frac{Δ t}{Δ z^{2}} + 1 - \sqrt{\frac{1 - ϵ}{2}} \\ sup_{c} {\partial (θ (ψ, c) c) / \partial c} & \leq L_{D_{q}} \frac{Δ t}{2 Δ z} + L_{D_{D}} \frac{Δ t}{Δ z^{2}} + 1 + \sqrt{\frac{1 - ϵ}{2}} \\ (21) & \frac{1}{L_{c}} \frac{\partial (θ (ψ, c) c)}{\partial ψ} & \leq \sqrt{\frac{1 - ϵ}{2}} \end{aligned}

Adding the inequalities (18) and (19) and using (20) and (21) one gets

\begin{matrix} (22) & {(e_{ψ}^{s + 1})}^{2} + {(e_{c}^{s + 1})}^{2} \leq (1 - ϵ) [{(e_{ψ}^{s})}^{2} + {(e_{c}^{s})}^{2}] \end{matrix}

With

{‖ e^{s} ‖}^{2} = \sum_{i = 1}^{I} [{(ψ_{i, k}^{s} - ψ_{i, k})}^{2} + {(c_{i, k}^{s} - c_{i, k})}^{2}]

being the squared error norm of the solution vector

x^{s} = (ψ_{k}^{s}, c_{k}^{s}) \in R^{2 I}

, the inequality (22) implies

\begin{matrix} (23) & {‖ e^{s + 1} ‖}^{2} \leq (1 - ϵ) {‖ e^{s} ‖}^{2} \end{matrix}

and thus a contractive property of the coupled

L

-schemes. Further, (23) can be used to show that

{x^{s}}

is a convergent Cauchy sequence (see e.g. [6, Sect. 8.1]) and its limit is unique, because

R^{2 I}

is a complete metric space.

For given

Δ z

, the time step

Δ t

is constrained by the stability conditions of the flow and transport schemes. With the positive constant

α

defined as

α = max (2 r_{i \pm 1 / 2, k}^{s})

, the stability condition (6) for the flow scheme becomes

α \leq 1

. Further, from

α = 2 max (K (ψ_{k}^{s}) Δ t_{ψ} / (L_{ψ} Δ z^{2})

, we determine the corresponding time step

Δ t_{ψ}

\begin{matrix} (24) & Δ t_{ψ} = \frac{α L_{ψ} Δ z^{2}}{2 max (K (ψ_{k}^{s}))} \end{matrix}

Similarly, with

β = max (r)

, the stability condition (14) becomes

β \leq 1

β = 2 D Δ t_{c} / L_{c} Δ z^{2}

, and the time step

Δ t_{c}

is given by

\begin{matrix} (25) & Δ t_{c} = \frac{β L_{c} Δ z^{2}}{2 D} \end{matrix}

The schemes (4) and (13) are coupled by alternating flow and transport steps at every time iteration [4, 16], which require a unique time step, defined as

\begin{matrix} (26) & Δ t = min (Δ t_{ψ}, Δ t_{c}) \end{matrix}

Using (24)-(26) we introduce the factor

γ

\begin{matrix} (27) & γ = \frac{Δ t}{Δ z^{2}} = min (\frac{α L_{ψ}}{2 max (K (ψ_{k}^{s}))}, \frac{β L_{c}}{2 D}) \end{matrix}

Inserting (27) and

γ Δ z

in (20) and (21), the convergence conditions explicitly depend on the parameters

L_{ψ}

and

L_{c}

and become more robust, with only one term (coming from the advection term of (2)) depending on discretization through

Δ z

. By relating the parameters

L_{ψ}

and

L_{c}

to the physical parameters of the problem, the conditions (20) and (21) also provide a general frame for an adaptive

L

-scheme.

We note that, even though the conditions (20) and (21) show that the coupled

L

-schemes should converge, they are of limited use in practice: actually
they only can be verified if the solution of the problem is known. Therefore, in the following we investigate the convergence by using sequences of iterates and the definitions of the computational and error-based orders of convergence.

3. NUMERICAL EXAMPLES

We illustrate the explicit iterative methods presented in Section 2 by three numerical examples. First, we solve a problem with transition between unsaturated and saturated flow regimes for the degenerate Richards' equation with both the FD and the BGRW scheme presented in Section 2.1. Then, we consider a fully coupled flow and surfactant transport problem [5] with an academic example of constitutive laws

θ (ψ)

and

K θ (ψ)

, and a one way coupled flow and transport problem with the realistic van Genuchten-Mualem parametrization of

θ

and

K

. The last two problems are solved by the coupled FD-flow and BGRW-transport schemes presented in Section 2.1 and Section 2.2.

The convergence of the iterative schemes is investigated numerically with tools from the theory of abstract convergent sequences [2]. The convergence of an arbitrary sequence of real numbers

x_{s} \to x^{*} \in R

is characterized by the behavior of the successive errors

e_{s} = | x^{*} - x_{s} |

. The sequence

{x_{s}}

converges with the (classical)

C

-order

p \geq 1

\begin{matrix} (28) & (C) lim_{s \to \infty} \frac{e_{s + 1}}{{(e_{s})}^{p}} = Q_{p} \in {\begin{cases} (0, \infty), & p > 1 \\ (0, 1), & p = 1, \end{cases} \end{matrix}

which implies the (weaker)

{order}^{2} Q

(see [3]),

\begin{matrix} (29) & (Q) lim_{s \to \infty} \frac{\ln e_{s + 1}}{\ln e_{s}} = p \end{matrix}

While Eq. (29) is a common way to identify the

Q

-order of convergence, Eq. (28) can subsequently be used to check whether a classical convergence order exists.

If the limit

x^{*}

is not known, the errors

e_{s}

are replaced by the corrections

| x_{s + 1} - x_{s} |

and (28)-(29) define "computational orders of convergence", denoted by a prime,

C^{'}

and

Q^{'}

. If

p > 1

, computational and error-based orders of convergence are equivalent, denoting this by using curly braces, and they are related by

{C, C^{'}} \underset{⇍}{\Rightarrow} {Q, Q^{'}} .

However, when

p = 1

, the above equivalences do not hold in general [3] and the existence of the error-based orders has to be checked on a case by case basis.

Results from the theory of convergent sequences in

R^{n}

[2] cannot be directly extrapolated to flow and transport iterative schemes for porous media, where

the convergence is investigated in function spaces with norms that often depend on the particular problem formulation (e.g., [7, 9, 12]). Therefore we rely on the sequences of positive real numbers provided by iterative methods, e.g.,

‖ ψ^{s} - ψ^{s - 1} ‖

and

‖ c^{s} - c^{s - 1} ‖

. If the method converges, which in a restricted sense means that the successive corrections become arbitrarily small (e.g., [9]), the corrections form sequences

{x_{s}}

converging to zero that can be analyzed with the definitions of the error-based orders of convergence (28) and (29) [19]. Further, based on numerical evaluations of the

Q

C

orders of convergence of the sequences

{x_{s}}

we may assume the

Q^{'}

C^{'}

convergence of the solution (

ψ^{s}, c^{s}

) of the coupled problem.

Remark 5. Let

s^{*}

be the iteration count for which (28) written for the sequence of successive corrections

‖ ψ^{s + 1} - ψ^{s} ‖

defines a constant quotient

Q_{1} (s^{*}) = Q

within a given absolute tolerance

ε_{a}

. Then, from (28), we have

\begin{aligned} ‖ ψ^{s + 1} - ψ^{s} ‖ & \leq Q ‖ ψ^{s} - ψ^{s - 1} ‖ \leq Q^{2} ‖ ψ^{s - 1} - ψ^{s - 2} ‖ \leq \dots \\ (30) & \leq Q^{s - s^{*}} ‖ ψ^{s^{*} + 1} - ψ^{s^{*}} ‖ . \end{aligned}

Using (30) we also have, for any

s > s^{*}, m \in N

\begin{aligned} ‖ ψ^{s + m} - ψ^{s} ‖ = \\ = ‖ ψ^{s + 1} - ψ^{s} + ψ^{s + 2} - ψ^{s + 1} + \dots + ψ^{s + m} - ψ^{s + m - 1} ‖ \\ \leq (Q^{s - s^{*}} + Q^{s - s^{*} + 1} + \dots + Q^{s - s^{*} + m - 1}) ‖ ψ^{s^{*} + 1} - ψ^{s^{*}} ‖ \\ = Q^{s - s^{*}} (1 + Q + \dots + Q^{m - 1}) ‖ ψ^{s^{*} + 1} - ψ^{s^{*}} ‖ \\ \leq Q^{s - s^{*}} (\sum_{m = 0}^{\infty} Q^{m}) ‖ ψ^{s^{*} + 1} - ψ^{s^{*}} ‖ \\ (31) & = Q^{s - s^{*}} lim_{m \to \infty} \frac{1 - Q^{m}}{1 - Q} ‖ ψ^{s^{*} + 1} - ψ^{s^{*}} ‖ = \frac{Q^{s - s^{*}}}{1 - Q} ‖ ψ^{s^{*} + 1} - ψ^{s^{*}} ‖ . \end{aligned}

According to (31),

lim_{s \to \infty} ‖ ψ^{s + m} - ψ^{s} ‖ = 0

, hence

{ψ^{s}}

is a convergent Cauchy sequence and has a unique limit in the complete space

R^{I}

. The convergence of the iterates of the

L

-scheme implies the cancelation of the stabilization terms, e.g.,

L (ψ^{s + 1} - ψ^{s})

, so that the limit of the sequence

{ψ^{s}}

coincides with the solution of the nonlinear scheme (compare (3) and (4)). Thus, the

C^{'}

-linear convergence of the

L

-scheme is equivalent to the convergence of the successive approximations

{ψ^{s}}

to the unique solution of the nonlinear scheme.

REMARK 6. From a theoretical perspective the

C

-linear convergence of the

L

-scheme assumes the convergence of the

L

-scheme to the exact solution. However, the numerical solution is only an approximation of the exact solution depending on the fineness of the discretization. It follows that a limit quotient

Q_{1} < 1

in (28), which defines the

C

-linear convergence, cannot be practically reached in numerical simulations which, necessarily use finite discretizations. For instance, the errors with respect to the exact solution of the second order accurate schemes used in the examples below are of the order

O (Δ z^{2})

. As the
errors

e_{s}

in (28) reach this threshold, they cannot be decreased by further iterations. One obtains thus estimations of the quotient

Q_{1} \approx 1

, corresponding to the sublinear convergence of the iterations [2]. The approach to the exact solution is instead quantified by the decay of the errors with respect to the discretization. The common approach

[11, 1, 16, 19]

is to assess estimated orders of convergence (EOC) by the slope of the decay of the errors in the logarithmic scale after successively halving the discretization parameter

Δ z

\begin{matrix} (32) & EOC = \log (\frac{e_{l}}{e_{l + 1}}) / \log (2), l = 1, 2, \dots \end{matrix}

3.1. First example: Richards' equation. For the beginning, let us consider the one-dimensional Richards' equation (1) closed by one of the simplest models, with exponential dependence of the water content on the pressure head and linear dependence of the hydraulic conductivity on the water content [16],

\begin{matrix} (33) & Θ (ψ) = {\begin{cases} e^{α ψ}, & ψ < 0 \\ 1, & ψ \geq 0, \end{cases} \\ (34) & K (θ (ψ)) = {\begin{cases} K_{s a t} Θ (ψ), & ψ < 0 \\ K_{s a t}, & ψ \geq 0, \end{cases} \end{matrix}

where

Θ = (θ - θ_{res}) / (θ_{sat} - θ_{res})

is the normalized water content,

θ = θ_{sat}

and

K = K_{sat}

denote the constant water content respectively the constant hydraulic conductivity in the saturated regions (

ψ \geq 0

), and

θ_{res}

is the residual water content.

We solve a typical problem modeling the experimental study of the flow in variably saturated soil columns, formulated in the domain

[0, 2]

. The boundary conditions are specified by a constant pressure

ψ (0, t) = ψ_{0}

at the bottom of the soil column and a constant water flux

q_{0}

on the top. Together, these constant conditions determine the initial pressure distribution

ψ (z, 0)

as solution of the steady-state flow problem. For

t > 0

, the pressure

ψ_{0}

is kept constant at the bottom of the column and on the top the water flux is increased linearly from

q_{0}

q_{1}

until

t \leq t_{1}

and is kept constant for

t > t_{1}

. The texture of an homogeneous soil, representative for a column filled with sand, is parameterized with

K_{sat} = 2.77 \cdot 10^{- 6}, θ_{res} = 0.06, θ_{sat} = 0.36, α = 10

. The water fluxes on the top boundary are given by

q_{0} = 2.77 \cdot 10^{- 7}

and

q_{1} = 2.50 \cdot 10^{- 6}

. To capture the transition from unsaturated to saturated regime, the pressure at the bottom boundary is fixed at

ψ_{0} = 0.5

. For the parameters specified above, 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 use a regular lattice with

Δ z = 10^{- 2}

. The BGRW computations are initialized by multiplying the initial condition by the total number of particles

N

. Since the hydraulic conductivity varies in time according to (34), the time step is determined according to (24) for the maximum of

K

at every time
iteration and by specifying a maximum

α = 0.8

of the parameters

r_{i \pm 1 / 2, k}

. 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 θ / \partial t = 0

in (1)) and to

L = 2

for the solution of the non-stationary problem. In both cases, the convergence criterion is verified by choosing

ε_{r} = 0

and a fixed absolute tolerance

ε_{a} = 10^{- 18}

To illustrate the consistency of the representation of the solution by ensembles of random walkers and the equivalence of the BGRW and FD schemes, we compare their solutions for increasing

N

between

10^{3}

and

10^{24}

at the final time

T

(Fig. 1). The errors

| ψ^{B G R W} (z) - ψ^{F D} (z) |

of the BGRW solution with respect to the FD solution shown in Fig. 2 decrease as

N^{- 1}

and approach the machine precision

\sim 10^{- 16}

for

N \geq 10^{18}

. The

l^{2}

norm

‖ ψ^{B G R W} - ψ^{F D} ‖

of the difference of the two solution vectors with

I = 2 / Δ z + 1 = 101

components decreases at a rate of

\sim 10 N^{- 1}

and reaches a plateau at about

10^{- 14}

(see Table 1). Since the FD solution is the ensemble average of the BGRW solution (see Remark 1), this behavior also illustrates the self-averaging property of the GRW algorithms [21].

$N$	1 e 3	1 e 6	1 e 10	1 e 18	1 e 24
$\| ψ^{B G R W} (z) - ψ^{F D} (z) \|$	$1.70 e - 03$	$4.24 e - 06$	$2.52 e - 10$	$7.13 e - 16$	$9.21 e - 16$
$‖ ψ^{B G R W} - ψ^{F D} ‖$	$2.27 e - 02$	$6.52 e - 05$	$3.32 e - 09$	$1.13 e - 14$	$1.40 e - 14$

Table 1. Mean errors and error norms for Section 3.1, Fig. 2.

Let

{ψ^{s}}

be the sequence of iterates of the BGRW

L

-scheme and

ψ^{B G R W} = lim_{s \to \infty} ψ^{s}

. By the triangle inequality

‖ ψ^{s} - ψ^{s - 1} ‖ \leq ‖ ψ^{s} - ψ^{F D} ‖ + ‖ ψ^{F D} - ψ^{s - 1} ‖

we have

lim_{s \to \infty} ‖ ψ^{s} - ψ^{s - 1} ‖ \leq 2 ‖ ψ^{B G R W} - ψ^{F D} ‖

. Hence, the order

\sim 20 N^{- 1}

of convergence of the BGRW solution to the FD solution is an upper limit for corrections. Fig. 3 and Fig. 4 show the decrease of the correction norms for

N = 10^{10}

and

N = 10^{18}

, respectively. In the first case the corrections reach a plateau of about

10^{- 9}

, consistent with the upper bound

20 N^{- 1}

. In the second case, the corrections show a plateau at the machine precision, which, as indicated by Fig. 2, is the same as the plateau for the FD scheme.

Remark 7. The convergence behavior illustrated above can be explained with the reduced fluctuations Algorithm 1. The sum of fractional parts rest

= \sum_{i = 1}^{I} {δ_{i}}

can be used to evaluate the truncation errors of the algorithm. At a given lattice site

i, {δ_{i}} \in [0, 1)

with a mean

\overset{―}{{δ_{i}}} \approx 0.5

. The sum over the lattice sites is thus rest

\approx 0.5 I

. The total numbers of particles on the lattice is ntot

= \sum_{i = 1}^{I} n_{i, k}^{s} = \sum_{i = 1}^{I} | ψ_{i, k}^{s} | N = \overset{―}{| ψ |} N I

. With these, the truncation error can be evaluated as rest/ntot

\approx 0.5 / \overset{―}{| ψ |} N^{- 1}

, which is independent of the
discretization

^{3}

(i.e., the number of lattice sites I). For verification purposes, we remove all the floor functions in the Code 1, add the truncation sequence rest=nn+rest; nn=floor(rest); rest=rest-nn;
at the end, and perform the computation with

N = 10^{10}

. The relative error of the solution

ψ_{T / Δ t}

obtained in this way with respect to the solution of the original Code 1, evaluated with the

l^{2}

norm, is of about

4.25 \cdot 10^{- 4}

. At the middle of the time interval

[0, T]

we evaluate

\overset{―}{| ψ |} = 6.45 \cdot 10^{- 2}

and obtain the truncation error

0.5 / \overset{―}{| ψ |} N^{- 1} = 7.75 \cdot 10^{- 10}

, which corresponds to the average value of the plateau shown by the corrections

‖ ψ^{s} - ψ^{s - 1} ‖

in this case. Since the original Code 1 contains three truncation sequences, we take rest

/ I \approx 1.5

and obtain the truncation error

1.5 / \overset{―}{| ψ |} N^{- 1} = 2.30 \cdot 10^{- 9}

. This error is almost the same as the plateau of the correction from Fig. 3. The plateau arises when the correction norms reach the threshold given by the truncation errors and cannot be further decreased by the iterative method. We note that this threshold is also close to the relative errors of the BGRW solution with respect to the FD solution at the final time

T

shown in Fig. 2 and Table 1. Thus, the accuracy of the BGRW L-scheme is limited by the truncation errors produced by the reduced fluctuations Algorithm 1.

The plateau of the order

10 N^{- 1}

also arises if one adds to the three terms depending on the parameter

r

in the FD Code 2 a noisy term

randn / N

, where randn is a random variable drawn from the standard normal distribution. The results (not shown here) are similar to those from Fig. 3 and Fig. 6. One expects therefore that the accuracy of the random walk methods based on random number generators is limited as well by truncation errors.

Fig. 5 and Fig. 6 show that the length of the sequence of corrections for

N = 10^{10}

not affected by oscillations is too short to allow the estimation of the order of convergence

p

with (29) and the evaluation of the quotient

Q_{1}

according to (28). For

N = 10^{18}

, the sequence not affected by oscillations is three times longer and the behavior of the estimate

p (s) \approx 1

and of the quotient

Q_{1} (s) < 1

suggest the

C^{'}

-linear convergence (see Fig. 7 and Fig. 8). This is, according to Remark 5, a numerical argument for the convergence of the iterates to the unique solution which solves the nonlinear scheme (3).

In the absence of an exact solution for this example, we investigate the existence of the "error-based" order

C

on sequences of errors with respect to the FD solution taken as a reference. We compute BGRW solutions for

N = 10^{16}

and form the sequence

{‖ ψ^{s} - ψ^{F D} ‖}

. The errors presented in Fig. 9 maintain a linear trend in semi-logarithmic representation over 130 iterations

s

and allow the evaluation, according to (28), of a subunit quotient

Q_{1}

in the same interval, shown in Fig. 10. We can therefore assume that the solution of the BGRW scheme converges with the linear

C

-order to the reference FD solution.

Fig. 1. Section 3.1: comparison of the BGRW (5) and FD (4) solutions.

Fig. 3. Section 3.1: corrections

(N = 10^{10})

Fig. 5. Section 3.1: estimated order of convergence (

N = 10^{10}

Fig. 2. Section 3.1: errors of the BGRW solutions.

Fig. 4. Section 3.1: corrections

(N = 10^{18})

Fig. 6. Section 3.1: estimation of the quotient

Q_{1} (N = 10^{10})

3.2. Second example: fully coupled flow and surfactant transport. The next example is the fully coupled problem of variably saturated flow, with transition from saturated to unsaturated regime, and surfactant transport

Fig. 7. Section 3.1: estimated order of convergence (

N = 10^{18}

Fig. 9. Section 3.1: errors of the BGRW solution (

N = 10^{16}

Fig. 8. Section 3.1: estimation of the quotient

Q_{1} (N = 10^{18})

Fig. 10. Section 3.1: error-based quotient

Q_{1}

of the BGRW solution.

governed by Eqs. (1) and (2) with the exact manufactured solutions

ψ_{m} (z, t) = - t z (z - 1) + \frac{x}{4}, c_{m} (z, t) = t z (z - 1) + 1

and the water contend parameterized by

θ (ψ, c) = {\begin{cases} \frac{1}{1 - ψ - c / 10}, & ψ < 0 \\ 0.3, & ψ \geq 0 . \end{cases}

The hydraulic conductivity is constant

K (θ (ψ, c)) = 1

and the diffusion coefficient in Eq. (2) is set to

D = 1

(in arbitrary units). The manufactured solutions verify Eqs. (1) ans (2) with supplementary source terms

[10, 1, 16, 17]

. These source terms, generated by inserting the manufactured solutions into the equations, are calculated analytically with the Matlab symbolic toolbox.

We first solve the problem in the domain

[0, 1]

, on a regular lattice with

δ z = 2.5 \cdot 10^{- 2}

, with the coupled FD scheme (4) and BGRW transport scheme (13). The time step

Δ t = 1.56 \cdot 10^{- 3}

is determined according to (26) with

L_{ψ} = L_{c} = 100, α = 1

, and

β = 0.5

. The corrections are evaluated in the

l^{2}

norm with the absolute tolerance

ε_{a} = 10^{- 18}

and the convergence is assessed with the condition

max (‖ ψ^{s} - ψ^{s - 1} ‖, ‖ c^{s} - c^{s - 1} ‖) \leq ε_{a}

The results obtained for

N = 10^{10}

particles used in the transport scheme presented in Figs. 11 and 12 show that the truncation errors generated by the transport Code 3 yield a plateau

\sim 10^{- 10}

of the concentration corrections and also, due to the coupling, a plateau higher than the machine precision of about

10^{- 13}

of the pressure corrections. The oscillations arising when the corrections reach the plateau do not allow reliable estimations of the order of convergence and of the trend of the quotient

Q_{1}

, for both pressure (Figs. 13 and 14) and concentration (Figs. 15 and 16). Increasing the number of particles in the transport scheme to

N = 10^{24}

lowers the plateau to the machine precision for both pressure and concentration corrections (Figs. 17 and 18), renders more reliable the estimates of the order of convergence

p = 1

(Figs. 19 and 21), and indicates limits

Q_{1} < 1

, consistent with the

C^{'}

linear convergence (Figs. 20 and 22). According to Remark 5, this result is an indication of the convergence to the solution of the nonlinear schemes (3) and (11).

Tests (not presented here) with a similar two dimensional problem [16, Sec. 5.2.1] also show the plateau

\sim 10^{- 13}

of the pressure corrections produced by the coupling of flow and transport for

N = 10^{10}

particles in the transport code and a plateau at the machine precision for

N = 10^{24}

particles.

Fig. 11. Section 3.2: corrections of the pressure solution

(N = 10^{10})

Fig. 12. Section 3.2: corrections of the concentration solution (

N = 10^{10}

The convergence of the iterations does not yet ensure the convergence of the numerical solution to the manufactured solution. This can be verified by a grid convergence test performed by successively halving the space step

Δ z

. To do that, we fix the number of particles in the transport code to

N = 10^{24}

, which, as seen above, ensures that the truncation errors do not affect the convergence of the corrections, and perform iterations of the linearization schemes (4) and (13) until the corrections fall below the threshold

ε_{a} = 10^{- 18}

. The solutions for successive

Δ z

are further used to evaluate the EOC by the rate of logarithmic decay of the deviations from the manufactured solutions measured with the

Fig. 13. Section 3.2: estimation of the convergence order

p

, pressure (

N = 10^{10}

Fig. 15. Section 3.2: estimation of the convergence order

p

, concentration (

N = 10^{10}

Fig. 17. Section 3.2: corrections of the pressure solution

(N = 10^{24})

Fig. 14. Section 3.2: estimation of the quotient

Q_{1}

, pressure (

N = 10^{10}

Fig. 16. Section 3.2: estimation of the quotient

Q_{1}

, concentration

(N = 10^{10})

Fig. 18. Section 3.2: corrections of the concentration solution

(N = 10^{24})

discrete

L^{2}

norm

\sqrt{Δ z} ‖ \cdot ‖_{l^{2}}

. The results presented in Table 2 indicate the second order convergence of the solution

{ψ, q, c}

Fig. 19. Section 3.2: estimation of the convergence order

p

, pressure (

N = 10^{24}

Fig. 21. Section 3.2: estimation of the convergence order

p

, concentration (

N = 10^{24}

Fig. 20. Section 3.2: estimation of the quotient

Q_{1}

, pressure (

N = 10^{24}

Fig. 22. Section 3.2: estimation of the quotient

Q_{1}

, concentration

(N = 10^{24})

$Δ z$	$‖ ψ - ψ_{m} ‖$	EOC	$‖ q - m_{m} ‖$	EOC	$‖ c - c_{m} ‖$	EOC
$1.00 e - 1$	$4.59 e - 02$	-	$1.23 e - 01$	-	$3.43 e - 02$	-
$5.00 e - 1$	$1.14 e - 02$	2.01	$3.34 e - 02$	1.89	$1.02 e - 02$	1.75
$2.50 e - 2$	$2.85 e - 03$	2.00	$8.68 e - 03$	1.94	$2.76 e - 03$	1.88
$1.25 e - 2$	$7.11 e - 04$	2.00	$2.21 e - 03$	1.98	$7.21 e - 04$	1.94

Table 2. Orders of convergence of the numerical solution to Section 3.2.

The knowledge of the exact solution allows us to compute error-based orders of convergence as well. Using a finer discretization,

Δ z = 1.25 \cdot 10^{- 2}

, and a larger tolerance

ε_{a} = 10^{- 3}

we form sequences of discrete

L^{2}

norms of errors and corrections for pressure (Figs. 23 and 24) and concentration (Figs. 25 and 26). We note that as the convergence criterion for errors is fulfilled, the corrections are four to five orders of convergence smaller. Estimations (not shown here) of the order of convergence

p

with (29) were found to be consistent with the definition of the

Q

linear convergence, for both errors and the corresponding corrections. The quotient

Q_{1}

computed with (28) from the
sequences of successive errors approaches the unity, indicating the

C

sublinear convergence of the linearization schemes for both flow and transport (Figs. 27 and 29 ). We have thus an illustration of the observation made in Remark 6 that

C

-linear convergence cannot be obtained in numerical simulations. For the corresponding corrections,

Q_{1} < 1

indicates the

C^{'}

-linear convergence (Figs. 28 and 30), which, according to Remark 5, indicates the convergence of the iterative schemes to the solutions of the coupled nonlinear schemes (3) and (11).

Fig. 23. Section 3.2: errors of the pressure solution.

Fig. 25. Section 3.2: errors of the concentration solution.

Fig. 24. Section 3.2: corrections of the pressure solution.

Fig. 26. Section 3.2: corrections of the concentration solution.

3.3. Third example: One-way coupled flow and transport problem with van Genuchten-Mualem parameterization. In the following, we consider the one-way coupled system (1)-(2), i.e.,

θ

independent of

c

, with an experimentally based parameterization model. The relationships defining the water content

θ (ψ)

and the hydraulic conductivity

K (θ (ψ))

are given by the

Fig. 27. Section 3.2: estimation of the quotient

Q_{1}

, pressure errors.

Fig. 28. Section 3.2: estimation of the quotient

Q_{1}

, pressure corrections.

Fig. 29. Section 3.2: estimation of the quotient

Q_{1}

, concentration errors.

Fig. 30. Section 3.2: estimation of the quotient

Q_{1}

, concentration corrections.

van Genuchten-Mualem model

\begin{matrix} (35) & Θ (ψ) = {\begin{cases} {(1 + (- α ψ)^{n})}^{- m}, & ψ < 0 \\ 1, & ψ \geq 0, \end{cases} \\ (36) & K (Θ (ψ)) = {\begin{cases} K_{s a t} Θ (ψ)^{\frac{1}{2}} {[1 - {(1 - Θ (ψ)^{\frac{1}{m}})}^{m}]}^{2}, & ψ < 0 \\ K_{s a t}, & ψ \geq 0, \end{cases} \end{matrix}

where

Θ = (θ - θ_{res}) / (θ_{sat} - θ_{res})

is the normalized water content defined in the same way as for the exponential model (33)-(34),

α, n, m = 1 - 1 / n

are model parameters depending on the soil type, and

K_{sat}

is the constant hydraulic conductivity in the saturated regime. In this example, we consider the parameters of a silt loam used in previously published studies [7, 16, 17],

θ_{sat} = 0.396, θ_{res} = 0.131, α = 0.423, n = 2.06, K_{sat} = 4.96 \cdot 10^{- 2}

We chose the manufactured solutions

ψ_{m} (z, t) = - t z (z - 1) - \frac{x}{4} + \frac{x}{2}, c_{m} (z, t) = t z (z - 1) + 1

which illustrate again the degeneracy of Richards' equation, with transitions from saturated to unsaturated flow regime. The diffusion coefficient in the transport equation (2) is set to

D = 0.01

. The source terms generated by the manufactured solutions are calculated as in Section 3.2 above with the Matlab symbolic toolbox.

Since the flow is now independent of transport it is not influenced the truncation produced by the Code 3. Therefore, we focus here on the error-based orders of convergence of the linearization schemes. The coupled problems are solved in the domain

[0, 1]

, with

Δ z = 1.25 \cdot 10^{- 2}

and

Δ t = 3.91 \cdot 10^{- 3}

. The convergence is assessed with discrete

L^{2}

norms of pressure and concentration errors with respect to the manufactured solutions for the tolerance

ε_{a} = 10^{- 3}

. Preliminary tests with increasing numbers of particles show that the transport scheme (13) fails to converge if

N < 10^{8}

. To completely eliminate the effect to the truncation errors, the computations are carried out with

N = 10^{18}

particles.

As in the previous example, we found that the correction norms reach values that are systematically smaller by several orders of magnitude than the error norms, for both pressure and concentration and the estimates obtained with (29) indicate the

Q

order of convergence

p = 1

. The quotient

Q_{1}

computed with (28) approaches the unity for the errors of both the pressure and the concentration (Figs. 31 and 33) indicating the sublinear convergence. For the corresponding corrections (Figs. 32 and 34)

Q_{1} < 1

, hence the

C^{'}

linear convergence of both iterative schemes may be assumed. These results illustrate again the observations made in Remarks 5 and 6.

Fig. 31. Section 3.3: estimation of the quotient

Q_{1}

, pressure errors.

Fig. 32. Section 3.3: estimation of the quotient

Q_{1}

, pressure corrections.

The grid convergence test, conducted in the same way as in Section 3.2, presented in Table 3 indicates the quadratic convergence of the solution

{ψ, q, c}

towards the exact manufactured solution.

Fig. 33. Section 3.3: estimation of the quotient

Q_{1}

, concentration errors.

Fig. 34. Section 3.3: estimation of the quotient

Q_{1}

, concentration corrections.

$Δ z$	$‖ ψ - ψ_{m} ‖$	EOC	$‖ q - m_{m} ‖$	EOC	$‖ c - c_{m} ‖$	EOC
$1.00 e - 1$	$5.85 e - 02$	-	$2.28 e - 03$	-	$1.23 e - 02$	-
$5.00 e - 1$	$1.42 e - 02$	2.05	$5.69 e - 04$	2.00	$3.54 e - 03$	1.80
$2.50 e - 2$	$3.27 e - 03$	2.11	$1.48 e - 04$	1.94	$9.19 e - 04$	1.95
$1.25 e - 2$	$8.09 e - 04$	2.02	$4.50 e - 05$	1.72	$2.38 e - 04$	1.95

Table 3. Orders of convergence of the numerical solution to Section 3.3.

4. CONCLUSIONS

In this study, we demonstrated through numerical examples the convergence of the explicit

L

-schemes for nonlinear and degenerate problems of coupled flow and transport processes in porous media. The convergence issue is two-fold. First it is the convergence of the iterative method which solves the nonlinearity of the problem and provides a solution of the nonlinear numerical scheme, for instance the

L

-scheme for Richards' equation (4) which approximates the solution of the nonlinear scheme (3). The second aspect is the convergence of the numerical solution to the exact solution of the problem, provided that there exists a unique solution.

The convergence of the iterative schemes has been investigated through estimations of computational and error-based orders of convergence. The three examples presented in Section 3 provide numerical arguments for the

C^{'}

-linear convergence. According to Remark 5, the

C^{'}

-linear convergence indicates the uniqueness of the iterative solutions and their convergence to the solutions of the nonlinear schemes. The sequences of errors with respect to exact manufactured solutions of the coupled problem presented in Section 3.2 and Section 3.3 are only sublinearly

C

convergent and illustrate the impossibility of achieving the

C

-liner convergence via numerical simulations, pointed out in Remark 6. However, with the FD solution as reference, we have shown that the random walk BGRW

L

-scheme for Richards' equation could be

C

-linear convergent (Fig. 9). Using grid convergence tests and manufactured solutions
we get strong numerical evidence of second order convergence (

\sim (Δ z)^{2}

) of the

L

-schemes for coupled flow and transport problems (Table 2 and Table 3).

The random walk method for Richards' equation is rather of academic interest. It has been used, for instance, to prove the stability of the equivalent FD scheme (Remark 1). Instead, the random walk approach is essential in modeling the transport because it facilitates a discrete description at the molecular level of the chemical reactions [17] and provides a basis for the space-time upscaling of the reactive transport in porous media [18]. An issue of concern of this approach are the truncation errors produced by the Algorithm 1 which preserves the indivisibility of the particles. We have shown that the truncation errors induce a plateau in the decrease of the successive corrections of the linearization method (Remark 7). The plateau is inversely proportional to the number of particles,

\sim N^{- 1}

, and may hinder the evaluation of the convergence orders if

N

is not large enough. In modeling reactive transport, where the amount of molecules is of the order of Avogadro's number, the plateau reaches the machine precision. However, the truncation-limited precision becomes significant in applications for small systems of particles (e.g., [20]), for instance in modeling the isotopic separation of heavy water where the number of deuterium isotopes can be of the order

N \sim 10^{6}

or even smaller.

Acknowledgements. F.A. Radu acknowledges the support of the VISTA program, The Norwegian Academy of Science and Letters, and Equinor.

REFERENCES

[1] C.D. Alecsa, I. Boros, F. Frank, P. Knabner, M. Nechita, A. Prechtel, A. Rupp and N. Suciu, Numerical benchmark study for fow in heterogeneous aquifers, Adv. Water Resour., 138 (2020), 103558. https://doi.org/10.1016/j.advwatres. 2020.103558 중
[2] E. Cătinaş, A survey on the high convergence orders and computational convergence orders of sequences, Appl. Math. Comput., 343 (2019), pp. 1-20. https://doi.org/10. 1016/j.amc.2018.08.006
[3] E. Cătinaş, How many steps still left to

x^{*}

?, SIAM Rev., 63 (2021) no. 3, pp. 585-624. https://doi.org/10.1137/19M1244858
[4] D. Illiano, I.S. Pop and F.A. Radu, Iterative schemes for surfactant transport in porous media, Comput. Geosci., 25 (2021), pp. 805-822. https://doi.org/10.1007/ s10596-020-09949-2 준
[5] P. Knabner, S. Bitterlich, R. Iza Teran, R., A. Prechtel and E. Schneid, Influence of surfactants on spreading of contaminants and soil remediation, pp.152-161, in Jäger, W., Krebs, H.J. (Eds.), Mathematics-Key Technology for the Future, Springer, Berlin, Heidelberg, 2003. https://doi.org/10.1007/978-3-642-55753-8_12 [x
[6] P. Knabner and L. Angermann, Numerical Methods for Elliptic and Parabolic Partial Differential Equations, Springer, New York, 2003. https://doi.org/10.1007/b97419 중
[7] F. List and F.A. Radu, A study on iterative methods for solving Richards' equation, Comput. Geosci., 20 (2016) no. 2, pp. 341-353. https://doi.org/10.1007/ s10596-016-9566-3 중
[8] I.S. Pop, F.A. Radu and P. Knabner, Mixed finite elements for the Richards' equation: linearization procedure, J. Comput. Appl. Math., 168 (2004) no. 1, pp. 365-373. https: //doi.org/10.1016/j.cam.2003.04.008 중
[9] F.A. Radu, K. Kumar, J.M. Nordbotten and I.S. Pop, A robust, mass conservative scheme for two-phase flow in porous media including Hölder continuous nonlinearities, IMA Journal Numer. Anal., 38 (2018) no. 2, pp. 884-920. urlhttps://doi.org/10.1093/imanum/drx032
[10] P.J. Roache, Code verification by the method of manufactured solutions, J. Fluids Eng., 124 (2002) no. 1, pp. 4-10. http://dx.doi.org/10.1115/1.1436090 중
[11] C.J. Roy, Review of code and solution verification procedures for computational simulation, J. Comput. Phys., 205 (2005), pp. 131-156. https://doi.org/10.1016/j.jcp. 2004.10.036 중
[12] J.S. Stokke, K. Mitra, E. Storvik, J.W. Both and F.A. Radu, An adaptive solution strategy for Richards' equation, Comput. Math. Appl., 152 (2023), pp. 155-167. https://doi.org/10.1016/j.camwa.2023.10.020 [X]
[13] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, SIAM, 2004. https://doi.org/10.1137/1.9780898717938
[14] N. Suciu, Diffusion in Random Fields. Applications to Transport in Groundwater, Birkhäuser, Cham, 2019. https://doi.org/10.1007/978-3-030-15081-5 중
[15] N. Suciu, Global Random Walk Solutions for Flow and Transport in Porous Media, in F.J. Vermolen and C. Vuik (eds.), Numerical Mathematics and Advanced Applications ENUMATH 2019, Lecture Notes in Computational Science and Engineering, 139, Springer Nature, Switzerland, 2020. https://doi.org/10.1007/978-3-030-55874-1_93 준
[16] N. Suciu, D. Illiano, A. Prechtel and F.A. Radu, Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media, Adv. Water Resour., 152 (2021), 103935. https://doi.org/10.1016/j.advwatres.2021.103935 [
[17] N. Suciu and F.A. Radu, Global random walk solvers for reactive transport and biodegradation processes in heterogeneous porous media, Adv. Water Res., 166 (2022), 104268. https://doi.org/10.1016/j.advwatres.2022.104268 [
[18] N. Suciu, F.A. Radu and I.S. Pop, Space-time upscaling of reactive transport in porous media, Adv. Water Resour., 176 (2023), 104443. http://dx.doi.org/10.1016/ j.advwatres.2023.104443 [
[19] N. Suciu, F.A. Radu, J.S. Stokke, E. Cătinaş and A. Malina, Computational orders of convergence of iterative methods for Richards' equation, arXiv preprint arXiv:2402.00194 (2024), https://doi.org/10.48550/arXiv.2402.00194
[20] C. Vamoş, N. Suciu and M Peculea, Numerical modelling of the onedimensional diffusion by random walkers, Rev. Anal. Numér. Théor. Approx., 26 (1997) nos. 1-2, pp. 237-247. https://ictp.acad.ro/jnaat/journal/article/view/ 1997-vol26-nos1-2-art32/
[21] C. Vamoş, N. Suciu and H. Vereecken, Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comput. Phys., 186 (2003), pp. 527-544. https://doi.org/10.1016/S0021-9991(03)00073-1 중

Received by the editors: May 17, 2024; accepted: June 26, 2024; published online: July 11, 2024.

$^{1}$ Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Fântânele str. no. 57, Cluj-Napoca, Romania, e-mail: nsuciu@ictp.acad.ro, ecatinas@ictp.acad.ro.
$^{2}$ Center for Modeling of Coupled Subsurface Dynamics, University of Bergen, Allégaten 41, Bergen, Norway, e-mail: florin.radu@uib.no.
$^{1}$ In the fully coupled problem $θ (ψ_{i, k})$ is replaced by $θ (ψ_{i, k}, c_{i, k})$ , where $c_{i, k}$ is the solution of (2) and the coupled scheme consists of alternating flow and transport steps.
$^{2}$ This is usually denoted by $Q_{L}$ (from logarithm) to distinguish from other types of $Q$ orders. Since only definition (29) will be used in the following, we disregard the subscript $L$ .
$^{3}$ This can be easily verified by doubling the number of lattice sites and taking the appropriate parameter $L = 100$ in the Matlab code.

Iterative schemes for coupled flow and transport in porous media -Convergence and truncation errors

10 months ago

Abstract Nonlinearities of coupled flow and transport problems for partially saturated porous media are solved with explicit iterative L-schemes. Their…

Global random walk for one-dimensional one-phase Stefan-type moving-boundary problems: Simulation results

12 months ago

Abstract This work presents global random walk approximations of solutions to one-dimensional Stefan-type moving-boundary problems. We are particularly interested in…

Comment on “Modeling groundwater flow with random hydraulic conductivity using radial basis function partition of unity method” by Shile et al. (2024)

1 year ago

Abstract The goal of this article is to explain some striking discrepancies between the Monte Carlo inferences for flow in…

Posts by Nicolae Suciu

Abstract

Authors

Keywords

Paper coordinates

PDF

About this paper

Journal

Publisher Name

DOI

Print ISSN

Online ISSN

References

Paper in HTML form

ITERATIVE SCHEMES FOR COUPLED FLOW AND TRANSPORT IN POROUS MEDIA - CONVERGENCE AND TRUNCATION ERRORS

Abstract

1. INTRODUCTION

2. COUPLED FLOW AND TRANSPORT IN POROUS MEDIA

algorithm 1. Reduced fluctuations algorithm

3. NUMERICAL EXAMPLES

4. CONCLUSIONS

REFERENCES

Related Posts

Iterative schemes for coupled flow and transport in porous media -Convergence and truncation errors

Global random walk for one-dimensional one-phase Stefan-type moving-boundary problems: Simulation results

Comment on “Modeling groundwater flow with random hydraulic conductivity using radial basis function partition of unity method” by Shile et al. (2024)

Controlling numerical diffusion in solving advection-dominated transport problems

Computational orders of convergence of iterative methods for Richards’ equation

Space-time upscaling of reactive transport in porous media

Global Random Walk Solutions for Flow and Transport in Porous Media

Global random walk solvers for reactive transport and biodegradation processes in heterogeneous porous media

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

Some relations between the microscopic and the macroscopic modelling of thermodynamical processes