Abstract. Flow and multicomponent reactive transport in saturated/unsaturated porous media are modeled by ensembles of computational particles moving on regular lattices according to specific random walk rules. The occupation number of the lattice sites is updated with a global random walk (GRW) procedure which spreads the particles from a lattice site with computational costs comparable to those for a single random walk step in sequential procedures. To cope with the nonlinearity and the degeneracy of the Richards equation the GRW flow solver uses linearization techniques similar to the \(L\)-scheme developed in finite element/volume approaches. Numerical schemes for reactive transport, coupled with the flow solver via numerical solutions for saturation and water flux, are implemented in splitting procedures. Diffusion-advection steps are solved by GRW algorithms using either biased or unbiased random walk probabilities. Since the number of particles in GRW simulations can be as large as the number of molecules involved in chemical reactions, one avoids the cumbersome problem of rescaling particle densities to approximate concentrations. Reaction steps are therefore formulated in terms of concentrations, as in deterministic approaches. The numerical convergence of the new schemes is demonstrated by comparisons with manufactured analytical solutions. Coupled flow and reactive transport problems of contaminant biodegradation described by the Monod model are further solved and the influence of flow nonlinearity/degeneracy and of the spatial heterogeneity of the medium is investigated numerically.

Keywords. Richards equation Coupled flow and transport Linearization Iterative schemes Global random walk Reactive flow Monod model

MSC. 76S05 35K57 65N12 86A05 65C35

1 Introduction

The development of accurate numerical methods for reactive transport in natural porous media faces challenges posed by the need to solve coupled systems of nonlinear equations in conditions of highly heterogeneous properties of the medium. Besides the inherent nonlinearity of the chemical reactions, their coupling to the advection and dispersion processes influences the accuracy of the predictions based on numerical simulations. For instance, an increase of the dispersion, due to the numerical diffusion of the scheme [ Bause and Knabner, 2004 , Brunner et al., 2012 ] or to inappropriate parameterizations of the dispersion coefficients [ Cirpka et al., 1999 ] , leads to an artificial mixing of the species and to an overprediction of the biodegradation rate. Further, reduced accuracy of the flow solutions, in conditions of parabolic/elliptic degeneracy of the Richards equation, affects the accuracy of the transport and reaction simulations [ Radu et al., 2009 ] . Even more challenges may appear in the case of simulations of coupled processes in heterogeneous media, when the accuracy of the flow solutions is not only conditioned by the performance of the numerical method but also by the parameters and correlation structure of the random-function representation of the heterogeneity [ Alecsa et al., 2020 , Suciu, 2020 ] .

Most of the numerical approaches in modeling reactive transport in porous media are based on finite element [ Knabner et al., 2003 , Brunner et al., 2012 , Radu et al., 2008 ] or finite volume [ Cirpka et al., 1999 , Lie, 2019 , Illiano et al., 2020 ] schemes. In these approaches, complex reaction problems with large systems of mobile and immobile species are solved by reduction schemes using linear combinations to construct non-reactive components [ Kreutle and Knabner, 2007 ] . The numerical diffusion in advection-dominated transport problems is controlled by using higher-order schemes and adaptive time steps [ Bause and Knabner, 2004 ] , mixed-finite element methods using Lagrange multipliers [ Brunner et al., 2012 ] , or streamtube approaches [ Cirpka et al., 1999 ] .

An attractive alternative to solve advection-dominated problems is provided by various particle methods, which by construction are free of numerical diffusion. Particle approaches to reactive transport were developed with sequential particle tracking (PT) [ Benson and Meerschaert, 2008 ] and continuous time random walk (CTRW) [ Edery et al., 2010 ] schemes, as well as with sequential biased random walk cellular automata (CA) [ Karapiperis and Blankleider, 1994 , Karapiperis, 1995 ] and global procedures [ Izsák and Lagzi, 2010 , Nagy and Izsák, 2011 ] based on global random walk (GRW) algorithms [ Vamoş et al., 2003 ] .

Within the sequential PT and CTRW approaches a reaction takes place when two particles representing different molecular species are close enough to interact. The decision is made by using the “interaction radius”, an empirical parameter related to the pore/grain size [ Edery et al., 2010 ] , or by using the colocation probability, i.e., the probability that the two particles will occupy the same position, constructed as a convolution of two Gaussian probability densities [ Benson and Meerschaert, 2008 ] . The latter approach was used to simulate biodegradation processes in a column experiment [ Ding and Benson, 2015 ] and in a real field application [ Ding et al., 2017 ] . In these simulations, the concentrations estimated by simple binning of particle numbers were further smoothed by averaging over PT runs. The efficiency of the collocation approach is enhanced if instead of stochastic birth/death processes the mass of the particle is adjusted upon reaction while the number of particles is preserved [ Bolster et al., 2016 ] or by using mass transfer models [ Engdahl et al., 2017 ] . A different paradigm in colocation-based PT algorithms for chemical reactions uses particles which carry a compositions of species concentration conjointly with mass-transfer techniques [ Benson and Bolster, 2016 , Engdahl et al., 2017 , Schmidt et al., 2019 ] . This approach, also useful in modeling small-scale influences on mixing-limited reactive transport [ Sole-Mari et al., 2020 ] , is similar to the method of “notional particles” intensively used in “probability density function” approaches to model complex reaction systems [ Pope, 1985 ] . The computational costs of the method utilizing mass-transfer and particles with composition increases at least quadratically with the number of particles [ Schmidt et al., 2018 ] . However, since rather moderate numbers of particles are required for a good accuracy, the method overperforms classical PT approaches with birth/death processes.

Other PT approaches estimate molecular species concentrations by smooth particle distributions provided by kernel density estimations and simulate the reaction step stochastically, with reaction probabilities depending on the distance between particles and the bandwidth of the spatial kernel [ Fernàndez-Garcia and Sanchez-Vila, 2011 , Sole-Mari et al., 2019 ] . In advection-dominated problems, the effect of the numerical diffusion induced by spatial redistribution of particles into concentrations in control volumes [ Cirpka et al., 1999 , Cui et al., 2014 , Sanchez-Vila and Fernández-Garcia, 2016 ] may be reduced by using locally adaptive particle support volumes oriented along the local velocities, as for instance in simulations of biodegradation reactions in heterogeneous aquifers [ Sole-Mari and Fernàndez-Garcia, 2018 ] . Kernel density estimations have also been used in PT modeling with reduction schemes [ Sole-Mari et al., 2021 ] similar to those introduced in finite element approaches for multicomponent reactive systems [ ]KreutleandKnabner2007.

Besides the need to map particle densities into concentrations, the coupling between transport and nonlinear reactions is another issue of concern in PT approaches [ Sanchez-Vila and Fernández-Garcia, 2016 ] . Accounting for dying and new created particles together with the advective-diffusive displacement of each particle can be a cumbersome task in solving complex problems for nonlinear processes. For instance, in modeling network reactions and multi-rate mass transfer in heterogeneous systems the PT method is limited to first-order reactions [ Henri and Fernandez-Garcia, 2015 ] . In modeling biodegradation processes described by multiple-Monod expressions and geochemical reactions, an operator splitting procedure has been successfully used as an alternative to PT models based on stochastic simulations of the reaction step [ Cui et al., 2014 ] . In this approach, rather than adding new particles to account for mass changes, variations of reactant concentrations are accounted for by changing the mass of the existing particles. The positions of the original system of species-particles tracked at all time steps is used to map particle densities into concentrations. Further, after each PT step the reactions are modeled deterministically, by assuming that all particles in a grid cell are in contact and react to each other. This approach has been validated by comparisons with experimental data, as well as with analytical and numerical solutions, for both saturated and unsaturated flow regimes (without transitions between them). However, the method still faces some technical challenges, as for instance, representing the background oxygen concentration which fills the entire domain by adding/removing particles accounting for differences from background concentration, or the need to use a “particle-splitting” procedure to preserve the accuracy in cells with small numbers of particles.

In approaches based on random walk on regular lattices, there is no need to search for nearby groups of particles which can interact. A reaction takes place when a sufficient amount of particles representing reacting species meet at a lattice site [ Izsák and Lagzi, 2010 , Karapiperis, 1995 ] . For instance, in the CA approach, assuming that the local chemical equilibrium is attained on a time scale much shorter than advection and diffusion scales, both the reaction probabilities \(P_r\) and the rate of reaction production are proportional to the physical rate constant \(k\), e.g. \(P_r\sim k\) and \(n_{c}\sim kn_{a}n_{b}^2\), for the reaction \(a+2b\rightarrow c\) [ Karapiperis, 1995 ] . The CA reaction rules use functions of occupation numbers \(n_{a}\) and \(n_{b}\) and reaction probabilities \(P_r\) to decide whether a reaction takes place if an \(a\)-particle and two \(b\)-particles meet at the same lattice site. Even though these rules are not unique, they approximate the law of mass action, which is recovered under certain limiting conditions [ Karapiperis and Blankleider, 1994 ] . The species concentration is retrieved by multiplying the particle densities at lattice sites by appropriate scaling factors.

In GRW-type procedures proposed by [ Izsák and Lagzi, 2010 ] and [ Nagy and Izsák, 2011 ] , the reaction step is performed deterministically. Comparing with the PT approach with deterministic reaction step of [ Cui et al., 2014 ] described above, we note that now the positions of all the particles are globally updated at every time step by advective and diffusive displacements. Moreover, GRW procedures use sufficiently large ensembles of particles to ensure the desired accuracy in representing concentrations by particle densities at lattice sites. As an illustration, in case of the three-component reaction considered above, the amount of product concentration in the time step \(\Delta t\), \(\Delta c_c=\Delta t kc_{a}c_{b}^2\), is directly converted into the amount of \(c\)-particles by representing the species concentrations as numbers of particles at lattice sites normalized by the number of particles used to represent a mole [ Nagy and Izsák, 2011 ] .

The GRW approaches to reactive transport presented in [ Izsák and Lagzi, 2010 , Nagy and Izsák, 2011 ] do not incorporate advection and the diffusion step is performed by generating binomial random variables. The latter represent numbers of particles undergoing unbiased random walk jumps from a lattice site to neighboring sites. The procedure corresponds to the complete and exact GRW algorithm for diffusion introduced in [ Vamoş et al., 2003 ] which becomes expensive in terms of computing time, with a slow linear increase for numbers of particles larger than \(10^{10}\) [ Vamoş et al., 2003 , Fig. 9 ] . Therefore, a limited number of particles was used in these diffusion-reaction simulations (e.g. \(1.5\cdot 10^5\) particles per lattice site for each molecular species in [ Nagy and Izsák, 2011 ] ) and a scaling factor was necessary to convert numbers of particles into moles. The GRW solvers for reactive transport introduced below in this article rely on the “reduced fluctuations algorithm” which can manage arbitrarily large numbers of particles [ Vamoş et al., 2003 ] , also used in our recent papers [ Suciu, 2020 , Suciu et al., 2020 , Suciu et al., 2021 ] . In this way, one achieves fast and highly accurate simulations, which cannot be done with the sequential PT and the biased random walk CA approaches described above, which are limited by the increase of the computational effort with the number of particles. Moreover, the number of particles representing a mole can be set to the actual value of the Avogadro number, so that no scaling factors are needed to relate particle densities and species concentrations. Also, for such large ensembles of particles the GRW simulations are self-averaging [ Vamoş et al., 2003 ] , i.e. there is no need to perform averages over runs [ Ding and Benson, 2015 , Ding et al., 2017 ] or kernel density estimations [ Fernàndez-Garcia and Sanchez-Vila, 2011 , Sole-Mari et al., 2019 ] to obtain smooth concentrations as in PT procedures.

In the present study we propose two-dimensional GRW solvers for nonlinear reactive transport, with a special emphasis on aerobic biodegradation. The solvers consist of explicit numerical schemes which solve the transport and the reaction steps by an operator splitting procedure. The flow velocity and the water content are provided by the GRW solver presented in [ Suciu et al., 2021 ] , where the treatment of the nonlinear terms in the Richards equation is based on the \(L\)-scheme linearization approach [ Pop et al., 2004 , List and Radu, 2016 ] . We consider here only the one-way coupling of Richards and transport equations in which the concentrations do not influence the flow.

With the unbiased GRW, the transport is modeled by updating the occupation numbers of the computational lattice with a global procedure which relocates all the particles from a lattice site in a single numerical procedure. First, all the particles undergo advective displacements, according to the local velocity evaluated at the lattice site, then they are distributed to neighboring sites by computing the number of unbiased random walk jumps in the spatial directions of the lattice. Considering the movement of a single particle, we note that it is similar to the PT algorithm, with the difference that the advective displacement is discrete and the diffusive one is modeled as unbiased random walk, rather than by realizations of the Wiener process (normally distributed random numbers) as in PT approaches. In fact, the GRW algorithm is a superposition of PT schemes (or weak solutions of Itô equation) on a regular lattice [ Suciu, 2019 , Sect. 3.3 ] .

Instead, if the jump probabilities are biased, the bias being proportional to the local velocity, the random walk models both the advective and diffusive displacements. This results in a biased global random walk (BGRW) algorithm which generalizes the CA approach of [ Karapiperis and Blankleider, 1994 ] as a global procedure. In the BGRW algorithm the particles are allowed to jump only to first-neighbor sites. In this way, one avoids jumps over several lattice sites with significant differences in local velocity and diffusion coefficients which produce overshooting errors. Thus, the BGRW algorithm is more accurate than the unbiased GRW, which instead, by using a coarser spatial discretization is faster.

In both unbiased GRW and BGRW solvers, the reaction step is solved deterministically by evaluating the reaction terms in the coupled system of reactive transport equations. In case of saturated flows, the nonlinear reaction terms are evaluated with an explicit linearization by computing them with concentration values from the previous time steps of the numerical scheme. Since numerical tests indicate that the explicit linearization is no longer accurate for variable water content, we propose an iterative approach for reactive transport in variably saturated flows similar to the GRW \(L\)-schemes introduced in [ Suciu et al., 2020 , Suciu et al., 2021 ] .

The paper is organized as follows. In Section 2 we introduce the governing equations for coupled flow and reactive transport with Monod reactions. Section 3 contains the presentation of the GRW solvers for multicomponent reactive transport. The solvers are further evaluated in Section 4 by comparisons with analytical solutions and numerical estimations of orders of convergence for saturated and unsaturated flow regimes, as well as in case of degenerate Richards equation describing the transition between unsaturated and saturated flows. In Section 5 the new solvers are tested on benchmark problems for biodegradation processes governed by Monod type models in saturated flows often used in studies published in the past. Further, the reactive transport solver will be coupled to Richards equation with randomly distributed parameters and the influence of the variable saturation and of the spatial heterogeneity of the soil/aquifer system will be investigated through numerical experiments in Section 6. Section 7 contains the conclusions of this study. The GRW/BGRW codes implemented in Matlab used in this study are stored in the Git repository MonodReactions [ Suciu and Radu, 2021 ] .

2 A simple model for reactive transport and biodegradation

The microbial biodegradation is a natural process which facilitates the in situ bioremediation of contaminated soils and groundwater systems, either as an intrinsic process or as part of an engineered remediation technology [ Wiedemeier et al., 1999 ] . The process occurs when both an electron acceptor (e.g. oxygen) and a microbian population are available. The contaminant to be degraded (e.g. benzene) acts as an electron donor and as a carbon source utilized by microorganisms [ Bause and Knabner, 2004 ] . Since the availability of the three elements depends on the water flow and the diffusion/dispersion processes, reliable models should consider the coupling of the reaction system to the flow and transport equations [ Radu, 2004 ] .

We consider a two-dimensional setting with an open domain \(\Omega \), with boundary \(\Gamma \), and a time interval \((0,T]\). Taking into account both the saturated and the vadose zone, the water flux is given by Darcy’s law

where \(z\) is the height against the gravitational direction. The pressure head \(\psi (x,z,t)\) solves the Richards equation

where \(\theta \) is the water content and \(K\) stands for the hydraulic conductivity of the medium. Equation (2) is closed by material laws expressing dependencies between \(K\), \(\theta \), and \(\psi \), based on experimental results. Further, we consider three molecular species, of which the first two are mobile and the third is immobile. The evolution of the corresponding concentrations \(c_\nu \), \(\nu =1,2,3\) is governed by the following system of coupled nonlinear equations

The constant \(k_d\) in Eq. (4) is the decay rate of the immobile component. The factor of the reaction rate \(R_3\) may be used to model the growth limitation of \(c_3\) (due for instance to limited pore space available), \(c_{3_{max}}\) is the maximum immobile concentration, and the coefficient \(\gamma \) is set to one in the model with growth limitation and to zero otherwise [ Radu, 2004 , Bause and Knabner, 2004 ] . In a basic scenario for aerobic biodegradation [ Bause and Knabner, 2004 , Brunner et al., 2012 , Cirpka et al., 1999 , Radu, 2004 ] , \(c_1\) represents the concentration of the electron donor, \(c_2\) corresponds to the electron acceptor (oxygen), \(c_3\) is the biomass concentration, and the growth of the biomass is assumed to follow a double Monod kinetics. The reaction term of Eq. (4) is expresses as \(R_3=Y\mu \), \(Y\) is the microbial yield coefficient per unit electron donor consumed, and the Monod term is defined by

where \(\mu _{max}\) is the maximum growth rate and \(M_1\), \(M_2\) are Monod constants. The reaction terms for the mobile species in Eqs. (3) are expressed as \(R_1=-\theta \alpha _1\mu \) and \(R_2=-\theta \alpha _2\mu \), where \(\alpha _1\), \(\alpha _2\) are stoichiometric constants.

As compared with more complex biodegradation models [ Bause and Knabner, 2004 , Cui et al., 2014 , Radu, 2004 ] , the model presented above is simplified in that only the degradation of a single contaminant is considered. Also, the inhibition terms which yield a slower microbian growth and the sorption terms are neglected. Instead, the coupling of the basic model process of aerobic degradation described by (3)-(5) with the flow model (1)-(2) allows investigations on the role of the variable saturation and of the heterogeneity of the soil-aquifer system.

3 Splitting algorithms for multicomponent reactive transport

In the biodegradation model presented above, the advection-dispersion-reaction problem is coupled in one way with the flow problem through the solution \((\theta , \mathbf{q})\) of the flow equations (1)-(2). The system of equations (1)-(5) can be solved with methods based on linearization procedures developed for fully coupled flow and transport equations [ Illiano et al., 2020 , Knabner et al., 2003 , Suciu et al., 2021 ] after replacing the general constitutive law for the nonlinear water content \(\theta (\psi ,c)\) by the simpler law \(\theta (\psi )\), corresponding to the one-way coupling. The resulting scheme has to be completed by the computation of the coupled nonlinear reaction terms at every time step and iteration in the transport solver. With these modifications, the scheme for fully coupled flow and transport from [ Suciu et al., 2021 ] becomes an iterative GRW scheme for multicomponent reactive transport. We note that the operator splitting PT method for reactive transport in variably saturated media proposed in [ Cui et al., 2014 ] also uses an iterative Newton method to solve the nonlinear system of ordinary differential equations of chemical kinetics. In the saturated flow regime, \(\theta \) no longer depends on the pressure \(\psi \) and, in the examples presented in this article, it will be taken as a constant parameter. For given flow velocity, the GRW solution can be constructed, similarly to the solution for reaction-diffusion problems from [ Nagy and Izsák, 2011 ] , as an explicit non-iterative scheme which evaluates nonlinear reaction terms from concentrations of mobile species obtained after the advection-dispersion step and immobile species concentrations from the previous time step (i.e., by an explicit linearization [ Haverkamp et al., 1977 ] ).

3.1 Non-iterative schemes for reactive transport in saturated media

Biased GRW algorithm for the transport step

For constant \(\theta \), the transport equation (3) for one mobile concentration \(c\), with a source/sink term \(f\) independent of \(c\) added in the right-hand side, can be equivalently written as

The new term \(f\) of the transport equation will be used in code verification tests presented in Section 4 to describe the result of applying the differential operators of the transport equation to a prescribed analytical solution.

The flow velocity \(\mathbf{q}\) is computed with the GRW algorithm for flows in saturated/unsaturated porous media proposed by [ Suciu et al., 2021 , Sect. 4.1 ] , with velocity on the boundaries approximated by forward finite differences, for the particular case of constant \(\theta \). The transport step, described by Eq. (6) without the reaction term, can be computed with explicit GRW algorithms.

The two-dimensional BGRW algorithm described below, which simulates the advection through biased jump probabilities [ Suciu, 2019 , Sect. 3.3.3 ] , is a particular case of the algorithm for fully coupled flow and transport presented in [ Suciu et al., 2021 , Sect. 4.2.1 ] . We assume a diagonal diffusion tensor with constant components \(D_1\) and \(D_2\) and by \(U\) and \(V\) we denote the components of the Darcy velocity \(\mathbf{q}\) 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. For constant \(\theta \leq 1\) the algorithm uses normalized diffusion and drift coefficients \(D_1/\theta \), \(D_2/\theta \), \(U/\theta \) and \(V/\theta \) [ Suciu et al., 2021 , Remark 5 ] .

The BGRW algorithm approximates the concentration by the distribution of an arbitrarily large number \(\mathcal{N}\) of particles undergoing biased random walk jumps on first-neighbor sites of the computational lattice. The number \(n_{i,j,k}\) of particles which at discrete time points \(k\Delta t\) are located at lattice sites \((i\Delta x,j\Delta z)\) approximates the concentration by \(c_{i,j,k}\approx n_{i,j,k}/\mathcal{N}\).

With the aid of the dimensionless diffusion coefficients and velocity components

the biased jump probabilities are defined as \(P(i\pm 1,j\; |\; i,j,k)=\frac{1}{2}(r_x\pm u_{i,j,k})\), for right/left jumps, \(P(i,j\pm 1\; |\; i,j,k)=\frac{1}{2}(r_y\pm v_{i,j,k})\), for up/down jumps, and \(P(i,j\; |\; i,j,k)=[1-(r_x+r_y)]\) is the no-jump probability. With these, the number of particles at lattice sites is updated by

where \(\lfloor \cdot \rfloor \) denotes the floor function.

The contributions to \(n_{i,j,k+1}\) in Eq. (8) are obtained with the BGRW algorithm which, at every time point \(k\Delta t\), redistributes the particles from sites \((l\Delta x,m\Delta z)\) according to

The quantities \(\delta n\) used in (9) to design the number of staying particles and the numbers of particles undergoing jumps on the four first-neighbor sites are binomial random variables with mean values given by

The realizations \(\delta n\) are approximated by a “reduced fluctuations algorithm” [ Vamoş et al., 2003 , Suciu, 2019 ] as follows: (i) initialize the real variable \(\sum _{\delta }=0\); (ii) compute \(\delta n\) with the relations (10) verified by the mean; (iii) split \(\delta n\) into \(\delta n=\lfloor \delta n\rfloor +\varepsilon _{\delta }\) at every lattice site; (iv) approximate \(\delta n\) by \(\lfloor \delta n\rfloor \); (v) sum up the remainders \(\varepsilon _{\delta }\) over all the sites; (vi) allocate one particle to the first lattice site where \(\sum _{\delta }=\sum \varepsilon _{\delta }\geq 1\) and replace \(\sum _{\delta }\) by \((\sum _{\delta }-1)\); (vii) save \(\sum _{\delta }\) and repeat the sequence (ii)-(vi) at every \(k\).

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

From the definition (7) of the BGRW parameters and the constraints (11) it follows that the allowable value of the local Péclet number is limited by \(\text{P\'{e}}\le 2\) [ Remark 4]Suciuetal2021.

Unbiased GRW algorithm for the transport step

The unbiased GRW algorithm moves groups of particles by advective displacements and diffusion jumps over one or more lattice sites. That is, it acts as a superposition of particle tracking schemes on lattices. The particle density at lattice points approximates the solution of the Fokker-Planck equation consisting of the homogeneous part of transport equation (6). Hence the unbiased GRW algorithm provides weak solutions to the associated Itô equation [ Suciu, 2019 , Sect. 3.3.2 ] .

The unbiased GRW redistributes the particles from lattice sites according to the relation

where \(d\) is a constant amplitude of diffusion jumps and the dimensionless variables \(r_x\), \(r_z\), \(u\) and \(v\) are defined by

The particle distribution is updated at every time step by

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

Comparing the relations (10) and (15), one can see that to the staying probability and the biased jump probabilities of the BGRW algorithm correspond the probability of the advective displacement and the unbiased probabilities \(r_x/2\) and \(r_y/2\) of jumps in the two spatial directions of the GRW algorithm.

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

Deterministic reaction step

Proceeding as in [ Nagy and Izsák, 2011 ] , we express the species concentration in moles per lattice site. The classical unit of moles per liter can simply be obtained by dividing the number of molecules by a fixed volume in liters and by the number of molecules in a mole, as done, for instance, in the simplified molecular dynamics approach to chemical kinetics proposed by [ Higham and Kloeden, 2021 ] . Since the (B)GRW implementations are not sensitive to the increase of the number of particles (Remark 1) we represent a mole of substance by the Avogadro number \(N_A\approx 6.022\cdot 10^{23}\), hereafter rounded to \(\mathcal{N}=10^{24}\). The concentrations of the mobile molecular species \(\nu =1, 2\) are computed as \(c_{\nu ,i,j,k+1/2}=n_{\nu ,i,j,k+1}/\mathcal{N}\), where \(n_{\nu ,i,j,k+1}\) is the distribution of the computational particles after the execution of the transport step obtained with either the BGRW algorithm (8) or the unbiased GRW algorithm (14). With this representation of the species concentrations, the reaction step in (6), for the particular case of the Monod model (3)-(5), will be performed deterministically (see also [ Nagy and Izsák, 2011 ] ) according to

where \(c_{3,i,j,k}\) is the concentration of the immobile species at the previous time step.

Following the reaction step, the numbers of particles representing mobile species are updated by \(n_{\nu ,i,j,k+1}=c_{\nu ,i,j,k+1}\mathcal{N}\), \(\nu =1,2\) and the splitting procedure is resumed by executing the new flow and transport steps.

3.2 Iterative scheme for reactive transport in variably saturated media

Since in case of variably saturated soils the unbiased GRW transport solver often requires extremely fine discretizations to obtain acceptable resolution of the discrete velocity components given by (13) [ Suciu et al., 2021 , Sect. 5.2.3 ] we shall consider in the following only the biased algorithm. The iterative BGRW \(L\)-scheme is obtained from a backward-time centered-space finite difference scheme for the transport equation for mobile species, Eq. (3) without reaction term and with a source term \(f\) added in the right-hand side, by adding the stabilization term \(L(c_{i,j,k}^{s+1}-c_{i,j,k}^{s})\), where \(L\) is a dimensionless constant parameter and \(s\), \(s=1,2,\ldots \), is the iteration index [ Sect. 4.2.1]Suciuetal2021. The number of particles at lattice sites is now updated with the following explicit relation with respect to the iteration index

where \(g^{s}_{i,j,k}=\mathcal{N}f_{i,j,k}\Delta t/L-\left[\theta (\psi _{i,j,k}^{s})n_{i,j,k}^{s} -\theta (\psi _{i,j,k-1})n_{i,j,k-1}\right]/L\), with pressure head \(\psi \) provided by the flow solver.

The parameters of the BGRW \(L\)-scheme are defined by

The BGRW \(L\)-scheme is completed by the same relations (9)-(11) as in case of non-iterative BGRW scheme from Section ??, which now have to be verified at every iteration \(s\). The iterations start with \(c_{i,j,k}^{1}=c_{i,j,k-1}\) and are stopped when the discrete \(L^2\) norms of the numerical solution \(c_{k}^{s}=(c_{i,j,k}^{s})\) verify

for some given tolerances \(\varepsilon _a\) and \(\varepsilon _r\).

This BGRW \(L\)-scheme is embedded into an operator splitting procedure alternating nonreactive transport for all mobile species \(c_\nu \) and reaction steps according to Eq. (16) at every iteration level \(s\).

4 Code verification tests

In the absence of exact solutions to boundary value problems for advective-dispersive transport coupled with Monod models, the codes can be verified by comparisons with analytical solutions to less complex problems, see e.g. [ Cui et al., 2014 ] , where one uses solutions for linear degradation and saturated flows with constant velocity. Using manufactured analytical solutions, instead, allows code verification tests for nontrivial situations such as flows governed by degenerate Richards equation and one-way coupling with nonlinear reactions [ Radu et al., 2009 ] or fully coupled one-dimensional flow and transport with transition from unsaturated to saturated conditions [ Suciu et al., 2020 ] . Here, we adopt this approach to test the new schemes for problems of increased complexity: constant flow and nonlinear reactions, saturated flow coupled with Monod biodegradation reactions, and space-time variable flow with transition from unsaturated to saturated regime coupled with Monod reactions.

4.1 Constant velocity and bimolecular nonlinear reactions

To verify the correctness of the transport and reaction steps in the splitting procedure we consider the system of coupled equations (3) with given parameters \(U=0\), \(V=-1\), \(D_1=D_2=D=0.1\), \(\theta =1\), and reaction terms given by \(R_1=-c_{1}c_{2}^{2}\) and \(R_2=-2c_{1}c_{2}^{2}\). The problem is solved in the domain \(\Omega =(0,2)\times (0,3)\) and the final simulation time is \(T=1\). The manufactured analytical solutions, previously used in [ Radu et al., 2008 , Brunner et al., 2012 ] ,

induce supplementary source terms in the right hand side of Eqs. (3) and define the initial and boundary conditions.

Numerical solutions obtained with the non-iterative BGRW algorithm presented in Section ?? are computed on regular lattices successively refined by halving five times the initial step size \(\Delta x = \Delta z = 0.2\). The accuracy of the numerical solutions at the final time \(t=T\) is quantified by discrete \(L^2\) norms and the estimated order of convergence (EOC) is given by the slope of decreasing error norms in logarithmic scale [ Alecsa et al., 2020 , Suciu et al., 2021 ] .

The results of the verification test presented in Table 1 show a very good accuracy of the BGRW solver for reactive transport and indicate the expected convergence of order 2. A similar test was done for the one-dimensional version of the BGRW solver, with manufactured solution similar to those for the two-dimensional case, where the factors containing the \(x\)-variable were removed. In this case, it is found that the numerical solutions converge with the order 2 as well, but the \(L^2\) errors are generally one order of magnitude larger (similarly to the one-dimensional solutions for the coupled flow and transport problem presented in [ Suciu et al., 2020 , Tables 10 and 11 ] ).

For the coarsest lattice used in this test (\(\Delta x=0.2\)) the local Péclet number \(\text{P\'{e}}=|V|\Delta x/D=2\) takes the maximum allowable value (see Section ??). Hence, as \(\Delta x\) decreases \(\text{P\'{e}}{\lt}2\), the relations (11) are verified, and the BGRW algorithm is free of numerical diffusion [ Suciu et al., 2021 , Table 9 ] . For smaller \(D\) the advection dominates and \(\Delta x\) has to be decreased to fulfil the condition \(\text{P\'{e}}\le 2\), which would slow down the BGRW computations. The alternative is to use the unbiased GRW solver presented Section ?? which is not constrained by a maximum \(\text{P\'{e}}\). By varying the dimensionless velocity \(v=\left\lfloor \Delta t |V|/\Delta z+0.5\right\rfloor \) and the amplitude of diffusion jumps \(d\) we have two degrees of freedom which can be used to accommodate different values of the velocity and diffusion coefficient. On the other side, while boundary conditions can be easily formulated for \(v=0\) and \(d=1\) [ Vamoş et al., 2003 ] , if \(v{\gt}1\) and \(d{\gt}1\) there are groups of particles jumping across the boundary from interior sites (e.g., \(\delta n_{i,j-v-d\; \mid \; i,j,k}\) for \(j-v-d{\lt}1\) in the test problem considered here), thus inevitably inducing numerical errors. In such cases, very fine lattices would be necessary to reduce the approximation errors. Note that this issue is avoided in simulations of transport in unbounded domains [ Suciu et al., 2021 , Sect. 6.3 ] , which renders the unbiased GRW very efficient in large scale simulations of transport in groundwater [ Suciu, 2019 ] .

The smallest errors of the GRW solutions are obviously obtained if \(v=1\) and \(d=1\). For the setup of the present test problem this choice is not possible because if \(v=1\) the amplitude of the diffusion jumps in (13) should be at least \(d=2\) to comply with the restriction \(r_x+r_z\le 1\), for \(\Delta x=0.2\), and even larger \(d\) is needed for smaller \(\Delta x\). Instead, for strongly advection-dominated problems with \(D=10^{-4}\), as those considered in the next section, the choice \(v=1,\; d=1\) is possible for all \(\Delta x\) considered in Table 1. The results of the unbiased GRW solutions, with the new dispersion coefficient and letting the other parameters unchanged, are presented in Table 2.

The larger error norms and the slower convergence of order 1 shown in the table can be attributed to the inherent errors of the GRW solver due to jumps over the lower boundary of the computational domain.

4.2 Saturated flow and Monod reactions

We consider in the following the same computational domain and final time as in Section 4.1 above. A more challenging test problem is formulated by considering an unsteady saturated flow with pressure head given analytically by

together with a constant water content \(\theta =0.3\), constant hydraulic conductivity \(K=0.05\), and constant dispersion coefficients \(D_1=D_2=0.025\). The nonlinear reaction is governed by the double Monod model (5) with parameters \(\mu _{max}=10^{-3}\), \(M_1=2\), \(M_2=0.2\), stoichiometric coefficients \(\alpha _1=1\), \(\alpha _2=3\), analytical solutions for the election donor and election acceptor concentrations verifying Eqs. (3)-(4) given by (20)-(21), and constant biomass concentration \(c_3=1\).

The coupled system of equations (2)-(3) with source terms in the right-hand side, computed analytically by differentiating the solutions (20)-(??), is solved with the flow GRW \(L\)-scheme from [ Suciu et al., 2021 , Sect. 4.1 ] coupled to the splitting algorithm consisting of the non-iterative BGRW scheme introduced in Section ?? and the deterministic reaction step described in Section  (??). The one-way coupling is realized by performing at every time step computations of the flow solution followed by computations of the species concentrations, using the velocity components provided by the flow scheme. The flow solutions are approximated, with parameter \(L=1\) and tolerances \(\varepsilon _a=\varepsilon _r=5\cdot 10^{-7}\), by performing a large number of iterations of the flow scheme, lying between 750 and 4500 for different time points and lattice steps \(\Delta x\). The \(L^2\) errors and the orders of convergence, estimated as in the previous section by halving the space step, are presented in Table 3.

For comparison, the same problem is solved with the iterative BGRW \(L\)-scheme from Section 3.2. For stabilization parameter \(L=1\) and tolerances \(\varepsilon _a=\varepsilon _r=10^{-6}\) in (??), the convergence of the liniarization procedure for different \(\Delta x\) values is already established after about 25 iterations. The corresponding convergence results are shown in Table 4. We note that there is a similar convergence behavior of the iterative and non-iterative solvers. The results shown in Tables 3 and 4 can be further improved, at the price of larger computation time, by increasing the stabilization parameter \(L\), which leads to smaller error norms and larger EOC [ Suciu et al., 2021 , Table 4 ] .

Further, we solve the Eqs. (3)-(4) for reactive transport with a given analytical solution of the flow problem. The flow velocity is obtained by differentiating the pressure head (??), according to (1). The results obtained with the non-iterative BGRW scheme, presented in Table 5, show a much improved order of convergence. Hence, the smaller EOC values of the concentration solutions shown Tables 3 and 4 can be attributed to the coupling and the limited accuracy of the flow solver. We mention that almost the same convergence results (not shown here) were obtained with the iterative scheme but the computation was several times slower (30 s vs. 4 s for the non-iterative solver). Hence, the non-iterative solver is more appropriate to solve Monod-type problems in case of prescribed velocity.

4.3 Saturated/unsaturated flows and Monod reactions

To model both unsaturated and saturated flow conditions, we consider the pressure head solution

As an example of nontrivial parameters which cause the degeneracy of the Richards equation, we consider the variable water content and hydraulic conductivity given by

The analytical solutions for the mobile species are again given by (20)-(21) and the reactions are governed by the same Monod model as in Section 4.2.

For the beginning, we investigate the convergence behavior of the BGRW \(L\)-scheme in conditions of strictly unsaturated flow by removing the last three terms from (??), with variable \(\theta \) given by (24) and a constant hydraulic conductivity \(K=0.05\). With the same parameters \(L\), \(\varepsilon _a\), and \(\varepsilon _r\) as for the saturated case, the iterative procedure already converges after 50 to 100 iterations for the flow solution and 23 to 28 iterations for the concentration solutions. The convergence behavior shown in Table 6 is similar to that in the saturated case shown in Table 4. The noticeable difference from the saturated case is given by the much shorter computing time of cca. 5 min, vs. cca 3 h in saturated case. Since, as we have seen above, the computation of the reactive transport alone takes only a few seconds, this difference is caused by the slower convergence of the flow solver in the saturated case.

In the following, we consider the challenging situation where the hydraulic conductivity \(K\) is given by the step function (27). Now, both \(\theta \) and \(K\) take constant values in saturated regions and Richards equation degenerates into an elliptic equation. The one-way coupling of the two-dimensional flow and reactive transport \(L\)-schemes is tested by comparisons with the exact concentration solutions (20)-(21) and the flow solution (??), which captures the transition from unsaturated to saturated regime. The setup of the reaction system is the same as that used in the previous sections to test particular cases of the general problem. The system of coupled equations (1)-(4), with source terms in Eqs. (2) and (3) derived from the analytical solution, is solved numerically with the flow and transport \(L\)-schemes. With the same \(\varepsilon _a\), and \(\varepsilon _r\) as used above and stabilization parameter set to \(L=10\) in the flow solver and to \(L=1\) in the transport solver, the desired accuracy of the iteration procedure is attained after \(10^3\) to \(2\cdot 10^4\) iterations of the flow scheme and 30 to 65 iterations of the transport scheme. The total computation time for the four discretization levels considered in this test is about 6 hours.

Convergence results for the flow and transport solvers are presented in Tables 7 and 8, respectively. One remarks a tendency towards second order of convergence of the pressure head and the concentration solutions. The \(L^2\) error norms of the velocity solutions are smaller by one order of magnitude than those for pressure and concentration. However, their convergence order is smaller, in the range from 1.4 to 1.7.

We note that the influence of the stabilization parameter \(L\) on accuracy and order of convergence estimations [ Sect. 5.1]Suciuetal2021 has not been investigated in the present study. Also, we did not tested the unbiased GRW non-iterative scheme presented in Section ?? by comparisons with analytical solutions of Monod problems. As seen in Section 4.1, the implementation of the boundary conditions in unbiased GRW schemes may distort the estimated order of convergence. The performance of the GRW scheme will be evaluated in Section 5 below through comparisons with BGRW solutions.

5 Biodegradation in saturated flows

In the following we consider a saturated regime with completely decoupled flow and reactive transport processes. The water content \(\theta \) equals the porosity and enters as a constant in the BGRW and GRW parameters given by relations (7) and (13), respectively. The two GRW algorithms will be used to build up numerical solutions to biodegradation problems in conditions of advection-dominated transport and heterogeneous structure of the aquifer system.

5.1 Biodegradation and advection-dominated transport

The artificial diffusion induced by numerical schemes is a serious challenge in simulations of biodegradation processes. We refer to [ Radu et al., 2011 ] for a comprehensive study on numerical diffusion for different discretization schemes. An artificial mixing of the species leads to an overprediction of the degradation or even to a complete degradation of the contaminant, contrary to what would be observed in practice. Mixed finite element methods and higher order schemes, together with stabilization techniques were proposed in the past to cope with the advection-dominated transport. Some drawbacks of these approaches have been also reported, such as the occurrence of negative concentrations (in mixed finite element approaches [ Brunner et al., 2012 ] ) or troubles with the construction of efficient iterative solvers (when using quadratic finite element schemes and streamline upwind stabilization [ Bause and Knabner, 2004 ] ).

As an alternative, we propose GRW solutions to the benchmark problem used in the papers cited above. Besides being free of numerical diffusion (provided that \(\text{P\'{e}}\le 2\) in case of BGRW algorithm), the GRW solvers working with indivisible particles do not produce negative concentrations, the conversion between particle densities and concentrations being straightforward, by using the same number of particles as molecules involved in reactions. Moreover, for saturated flows conditions, as in case of this benchmark problem, we have the advantage of using non-iterative BGRW/GRW schemes which, as seen in Section 4.2 above, are faster than the iterative schemes.

The problem is formulated in the domain \(\Omega =(0,0.5)\times (0,1)\). The reactive transport of the mobile species is governed by Eqs. (3) with parameters specified by \(\theta =0.2\), \(U=0\), \(V=-1\), \(D_1=D_2=D=10^{-4}\). The Monod reaction term \(\mu \) is modeled by Eq. (5) with \(M_1=M_2=0.1\), the biomass is kept constant, \(c_3=1\), and the reaction terms in Eqs. (3) are given by \(R_1=\theta \alpha _1\mu \) \(R_2=\theta \alpha _2\mu \), with \(\alpha _1=5\) and \(\alpha _2=0.5\). The final simulation time is \(T=5\), when one expects that the process reaches the steady state [ Bause and Knabner, 2004 ] . As initial conditions, one chooses \(c_1(x,z,0)=0\) and \(c_2(x,z,0)=0.1\) in \(\Omega \). The contaminant (electron donor) is injected at the middle part of the inflow boundary, where no electron acceptor (oxygen) is present, i.e., \(c_1(x,z,t)=1\) and \(c_2(x,z,t)=0\) for \((x,z,t)\in (0.225,0.275)\times \{ 1\} \times (0,T]\). Further, one sets \(c_1(x,z,t)=0\) and \(c_2(x,z,t)=0.1\) for \((x,z,t)\in (0,0.5)\backslash (0.225,0.275)\times \{ 1\} \times (0, T]\). Finally, the condition of vanishing concentration gradients is imposed on the left, right, and bottom boundaries.

The BGRW computations for the present setup require an extremely fine discretization (\(\Delta x=2\cdot 10^{-4}\) and \(\Delta t=10^{-5}\) for \(\text{P\'{e}}=2\)) which would render the computation until the final time \(T=5\) prohibitive. Therefore, we first restrict the final time to \(T=1\) and compare the BGRW solution to that provided by the unbiased GRW solver. The latter is obtained with \(\Delta x=2\cdot 10^{-3}\) (\(\text{P\'{e}}=20\)), \(\Delta t=4\cdot 10^{-4}\), and GRW parameters \(v=1\), \(d=1\). The comparison is done in terms of concentrations averaged over lattice sites, \(\langle c_{\nu }\rangle ,\; \nu =1,2\), by relative errors \(\varepsilon _{c_{\nu }}=(\langle c_{\nu }\rangle -\langle c_{\nu }^{b}\rangle )/\langle c_{\nu }^{b}\rangle \), where \(c_{\nu }^{b}\) is the BGRW solution. The errors obtained, \(\varepsilon _{c_1}=-0.0125\) and \(\varepsilon _{c_2}=-0.0016\), validate the GRW solver for the parameters of the benchmark problem described above and \(T=1\). Further, the problem is solved for \(T=5\) with the unbiased GRW solver. The electron donor and acceptor concentrations at the final time shown in Fig. 1 are similar to those obtained by quadratic finite element in [ Bause and Knabner, 2004 , Fig. 2 ] . One remarks that the contaminant is advected by the constant flow, with negligible diffusion and biodegradation.

Even though the GRW solver does not produce numerical diffusion, the effect of an artificial mixing which causes an overprediction of the biodegradation can be illustrated by increasing the diffusion coefficient. To simulate an enhancement of the mixing process, we repeat the computations after setting the diffusion coefficient to \(D=2.5\cdot 10^{-3}\). The results are presented in Fig. 2. The concentration profiles are now much more diffusive, the oxygen is largely consumed, as shown by the small values \(c_2\), and the more intense biodegradation in the lower part of the domain leads to the underestimation of the contaminant concentration \(c_1\) at the outflow boundary.

5.2 Influence of the heterogeneity on the biodegradation process

As seen in Section 5.1 above, in case of localized reaction zones the increase of the diffusion coefficient results in an increased biodegradation of the electron donor contaminant and consumption of the oxygen. The impact of the artificial increase of the diffusion in numerical schemes can be severe especially in the case of biodegradation controlled by transverse mixing in heterogeneous media [ Cirpka et al., 1999 ] . Classical approaches to cope with the numerical diffusion in advection-dominated problems are based on streamtube reactive transport models [ Cirpka et al., 1999 ] , local mesh refinement [ Klöfkorn et al., 2002 ] , or higher-order finite element used together with stabilization techniques [ Bause and Knabner, 2004 ] . In the following we propose GRW solutions to transport and biodegradation problems similar to those considered in these papers and investigate the influence of varying the parameters describing the heterogeneity.

We consider a problem similar to that formulated in [ Cirpka et al., 1999 ] , with a two-dimensional domain \(\Omega =(0,80)\times (0,20)\), a uniform initial distribution of biomass, a continuous injection of contaminant, constant concentration of oxygen in the water entering over the inflow boundary and in the domain, excepting the injection well. The water flow is driven by the gradient of the piezometric pressure and a variable hydraulic conductivity modeled as a realization of a random field. In the following, the distances are measured in meters, the time in days, and the concentration in moles.

To facilitate the computation of Monte Carlo ensembles of solutions, the flow velocity is approximated by a Kraichnan procedure [ Suciu, 2019 , Appendix C.3.2.2 ] as sum of 100 cosine random modes, with amplitudes, wavenumbers and phased determined by the parameters of the random hydraulic conductivity \(K\). The log-hydraulic conductivity field \(\ln K\) is a statistically homogeneous random function with isotropic Gaussian correlation. Its variance belongs to the set \(\sigma ^2\in \{ 0.1, 0.5,1\} \) and for its correlation length we consider \(\lambda =1\) and \(\lambda =2\). The mean conductivity is set to \(\langle K\rangle =25\). With a decrease of \(0.8\) of the pressure head between the inflow and outflow boundaries and a gradient \(J=-0.8/80\), the mean velocity becomes \(U=-J\langle K\rangle =0.25\). The dispersion is parameterized by a constant coefficient \(D=0.01\). Ensemble averages are estimated over 256 realizations of the \(\ln K\) field. These parameters provide fairly good approximations of the random Darcy velocity field and of the stochastic advection-dispersion process [ Schwarze et al, 2001 , Eberhard et al., 2007 ] . Since the porosity does not play a major role in saturated flows it is set to \(\theta =1\).

The setup of the chemical reactions closely follows that used by [ Cirpka et al., 1999 , Bause and Knabner, 2004 ] , with the only difference that, instead of mg/l, the concentrations are now given in moles. The parameters of the Monod model are given by \(M_1=2,\; M_2=0.2,\; Y=0.09,\; K_d=0.05\), \(\mu _{max}=5\), \(\alpha _1=1,\; \alpha _2=3\) [ Cirpka et al., 1999 ] . Additionally, following [ Bause and Knabner, 2004 ] we consider a growth limitation in the Monod model with \(c_{3_{max}}=1\) and \(\gamma =1\) in Eq. (4). The initial concentration of the contaminant is set to \(c_1=2\) at lattice sites inside a square of side equal to 2 centered at \(x=2.5\) and \(y=10\) and to \(c_1=0\) otherwise. The oxygen concentration is initially set to \(c_2=5\) in \(\Omega \), excepting the injection region where it is \(c_2=0\). The continuous injection is simulated, as in [ Bause and Knabner, 2004 ] , by resetting after every time step both \(c_1\) and \(c_2\) to their initial values in the injection region. The initial concentration of the biomass is set to \(c_3=0.001\) in \(\Omega \). In [ Bause and Knabner, 2004 ] , the boundary conditions for the mobile species are set to the initial values of \(c_1\) and \(c_2\) on the top, left and bottom boundaries. To avoid boundary effects in the unbiased GRW solver (see Section 4.1) which could lead to vanishing concentration \(c_2\) close to boundaries we also fix \(c_2\) to the initial value on a stripe of width 1.5 along the same boundaries. No boundary conditions are needed on the right boundary, because if no particles are introduced from from exterior of \(\Omega \) both the GRW and BGRW algorithms simulate the physical outflow situation by removing the particles crossing the boundary.

In Fig. 3 we present the concentration profiles of the three species at the final time \(T=300\) days computed with the non-iterative BGRW solver for a fixed realization of the velocity field with \(\lambda =2\) and \(\sigma ^2=0.5\). The results are obtained on a lattice with \(\Delta x=\Delta y=0.0333\) which ensures a maximum Péclet number Pé=1.75 verifying the relations (11). The corresponding time step, according to (7) is \(\delta t=0.0231\). The BGRW solution is further used as reference for comparison with that provided by the non-iterative unbiased GRW solver for the same realization of the velocity field. The latter is obtained by using a much coarser discretization, with \(\Delta x=\Delta y=0.1250\) (maximum Pé=6.57) and \(\delta t=1\) corresponding to \(d=2\) and \(u=U\delta t/\delta x=2\) in (13). With these, the GRW solver achieves a speedup of \(10^4\) over the BGRW solver.

The GRW solutions shown in Fig. 4 reproduce the general trend of the concentration profiles of the BGRW solution (Fig. 3) and the shape of the reaction front. However, as indicated by the middle panel of Fig. 4, there is an overshoot of the electron-donor particles in the GRW solution, which produces reactions beyond the reaction fronts predicted by the reference BGRW solution shown in Fig. 3. The more extended profiles of the biomass concentration \(c_3\) shown in the lower panel of Fig. 4 also can be explained by the numerical overshooting artifact. The effect of the overshooting can be quantitatively significant. In the present case, the relative errors of the lattice-averaged GRW solutions with respect to the BGRW solution are \(\varepsilon _{c_{1}}=-0.3109\), \(\varepsilon _{c_{2}}=0.0520\), and \(\varepsilon _{c_{3}}=0.8256\), that is, one order of magnitude larger than in case of the uniform velocity field analyzed in Section 5.1. The increased errors can be attributed to the spatial variation of the advection velocity, requiring larger parameters \(d\) and \(v\) which give rise to enhanced overshooting. The negative value \(\varepsilon _{c_{1}}\) quantifies the overprediction of the degradation caused by overshooting.

Keeping in mind that the GRW solutions overpredict the biodegradation, they can still be used to investigate the effect of heterogeneity through Monte Carlo simulations. The species concentrations averaged over 256 realizations of the \(\ln K\) field with \(\lambda =1\), \(\sigma ^2=0.1\) and the corresponding variances are shown in Figs. 5 and 6. The averaged concentrations shown in Fig. 5 are similar to the results obtained with uniform velocity. But the relative errors of the lattice-averaged mean concentrations from Fig. 5 with respect to the BGRW solutions for a constant longitudinal velocity set to \(U=0.25\) are close to \(\varepsilon _{c_{\nu }}\) in the comparison of the GRW/BGRW solutions for a single realization of the \(\ln K\) field presented above. The large variance \(\sigma ^{2}_{c_2}\) in Fig. 6, consistent with the inhomogeneity of the electron acceptor concentration \(c_2\) in a single realization shown in Fig. 4, is also partially due the overshooting in the GRW solver. Averages over lattice sites of the ensemble-mean species concentrations and variances for different correlation lengths \(\lambda \) and variances \(\sigma ^2\) of the \(\ln K\) field, computed at the final time \(T\), are presented in Table 9. The mean species concentrations show some variation, of about 22% (\(\langle c_1\rangle \)), 25% (\(\langle c_2\rangle \)), and 12% (\(\langle c_3\rangle \)), but there is no systematic dependence on \(\lambda \) and \(\sigma ^2\). The variances \(\sigma ^{2}_{c_{\nu }}\) slightly increase with \(\sigma ^2\) but do not show a systematic dependence on \(\lambda \).

The comparison of \(c_1\), \(c_2\), and \(c_3\) for the first realization in the Monte Carlo ensembles computed for the six combinations of \(\lambda \) and \(\sigma ^2\), presented in Fig. 7, shows that the spatial variability of the velocity field mainly influences the meandering of the reaction fronts and the amount of overshooting. As seen in both columns of panels in Fig. 7, for fixed \(\lambda \) the meandering is more pronounced as the variance \(\sigma ^2\) increases. Comparing the two columns, we remark an increase in meandering amplitude as the correlation length \(\lambda \) of the \(\ln K\) filed increases. We also remark that the amount of overshooting indicated by zones of low concentrations \(c_2\) of the electron acceptor outside the reaction fronts, increases with both \(\lambda \) and \(\sigma ^2\).

6 Biodegradation in variably saturated soils

In designing bioremediation strategies to clean up contaminated soils, the interplay between heterogeneity, variable saturation, and soil properties has to be taken into account in numerical models. As an illustration, we simulate in the following the infiltration and biodegradation of benzene in conditions of partial saturation for two types of soils.

We consider the same two-dimensional domain \(\Omega =(0,2)\times (0,3)\) used in verification tests presented in Section 4. The infiltration is driven by a variable pressure, given by \(\psi (x,z,t\leq t_1)=-0.2+2.7t/t_1\) and \(\psi (x,z,t{\gt} t_1)=0.2\), on the open part \(\Gamma _1=(0.5,1.5)\times \{ 3\} \) of the top boundary and by communication with the fresh water at hydrostatic equilibrium on the open parts of the domain \(\Gamma _{21}=\{ 0\} \times (0,0.5)\) and \(\Gamma _{22}=\{ 2\} \times (0,0.5)\). The zero flow condition is imposed on the closed boundary \(\Gamma \backslash (\Gamma _1\bigcup \Gamma _{21}\bigcup \Gamma _{22})\) and the initial pressure corresponds to hydrostatic equilibrium, \(\psi (x,z,0)=0.5-z\) [ Radu, 2004 ] . Hereafter, the time is measured in days and the distances in meters. The intermediate time for the variable pressure condition is set to \(t_1=1\) and the simulations are conducted for increasing times \(T\).

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

where \(\Theta = (\theta - \theta _{res})/(\theta _{sat} - \theta _{res})\) is the normalized water content, \(\theta _{res}\) and \(\theta _{sat}\) are the residual and the saturated water content, \(K_{sat}\) is the hydraulic conductivity of the saturated soil, and \(\alpha \), \(n\), \(m=1-1/n\) are model parameters depending on the soil type.

We consider here two sets of soil parameters, presented in Table 10. They correspond to a silt loam and a Beit Netofa clay, respectively, previously used in GRW simulations of coupled flow and nonreactive transport [ Suciu et al., 2021 , Sect. 5 ] .

The soil heterogeneity is modeled by a fixed realization of a random saturated hydraulic conductivity with mean \(K_{sat}\). The latter is computed by the Kraichnan method described in [ Suciu, 2019 , Appendix C3.1.2 ] , with an anisotropic Gaussian model of correlation lengths \(\lambda _x=0.1\) \(\lambda _z=0.01\), and variance \(\sigma ^2=4\). The dispersion is parameterized by the constant coefficient \(D=0.003\).

At the initial time, the concentration of benzene (electron donor) is set to \(c_1=2\) mole on the injection boundary \(\Gamma _1\) and to \(c_1=0\) otherwise, the concentration of oxygen (electron acceptor) is \(c_2=0\) on \(\Gamma _1\) and \(c_2=5\) mole otherwise, and the biomass concentration is \(c_3=0.001\) mole in \(\Omega \). At \(t\geq 0\), the concentrations \(c_1\) and \(c_2\) are kept at the initial values on the open boundaries and no flux conditions are imposed on \(\Gamma \backslash (\Gamma _1\bigcup \Gamma _{21}\bigcup \Gamma _{22})\). With these, one simulates a continuous injection of benzene through \(\Gamma _1\) and one assumes that the contaminated incoming water does not contain oxygen and the outgoing water is cleaned instantly at the contact with the fresh water on \(\Gamma _{21}\bigcup \Gamma _{22}\). The reactions are governed by the Monod model with the same parameters as those used in Section 5.2 for the simulations of reactive transport in aquifers.

The simulations are performed with the iterative BGRW \(L\)-scheme coupled with the GRW flow solver from [ Suciu et al., 2021 ] , by using the stabilization parameter \(L=1\) in both the flow and the transport solvers, \(\Delta x=\Delta z=0.05\), \(\Delta t=8.68\cdot 10^{-5}\), for loam soil parameters, and \(\Delta t=5.25\cdot 10^{-3}\) for clay soil parameters. The convergence of the \(L\)-schemes is achieved after tens to hundreds of iterations for \(T{\lt}3\), with more iterations in loam soil simulations, and decreases to only a few iterations for larger \(T\). As \(T\) increases from \(T=3\) to \(T=120\), the computing time ranges between 8 min and 2.5 h for the loam soil and between 50 s and 8 min for the clay soil.

The progress of the biodegradation and saturation processes for loam and clay soils is shown in Fig. 8. In the loam soil simulations the biodegradation process is quite fast, the oxygen is largely consumed, and the saturation is almost complete at \(T=120\). In case of clay soil simulations, even though the pressure profiles are similar, only the lower part of the domain \(\Omega \) becomes saturated and the biodegradation process is much slower. In addition, we observe that in case of loam soil the heterogeneity of the soil hydraulic properties influences not only the flow pattern but also the biodegradation process, as shown by the irregular contours of the benzene, oxygen, and biomass concentration. In clay soil simulations instead, even after 120 days the heterogeneity does not disturb significantly the regular concentration profiles, which are rather similar to those for homogeneous soils. These differences, consistent to those observed in similar simulations for nonreactive transport presented in [ Suciu et al., 2021 , Sect. 5.2.3 ] , can be explained not only by the lower hydraulic conductivity of the loam soil, but also by differences in \(\theta _{sat}\), \(\theta _{res}\), and parameters of the van Genuchten-Mualem model.

Based on the present simulations, different biodegradation scenarios can be constructed by slightly modifying the numerical setup, e.g., the injection conditions, the heterogeneity of the medium, as described by the random \(K_{sat}\), the dispersion coefficients, or the Monod parameters.

7 Conclusions

The GRW algorithms for reactive transport belong to a larger class of particle methods for solving partial differential equations. These methods use random walkers or trajectories of the Itô process to represent the solution and are, in principle, free of numerical diffusion. However, the GRW solvers differ in two important respects from other particle approaches. The first one is the global procedure which updates the positions of all the particles from a lattice site with computational costs that are practically the same as for a single random walk or PT step. The procedure allows us to use the same numbers of particles representing species concentrations as the number of molecules involved in chemical reactions. Consequently, there is no need to use spatial or ensemble averaging to obtain smooth concentration fields, the reaction steps are performed deterministically, and one avoids the cumbersome issue of representing the background oxygen concentration by particles filling the entire domain in PT simulations of biodegradation processes. The second advantageous feature is the robust coupling between the flow and reactive transport solvers, with flow solutions, provided by a GRW solver as well, computed at the sites of the lattice used in GRW transport simulations. One avoids thus the interpolation errors which occur when, for instance, the velocity field is computed with an external flow solver. This integrated approach is particularly useful in solving systems of coupled nonlinear flow and transport equations governing processes in saturated/unsaturated soils which require linearization procedures.

GRW schemes using unbiased random walk probabilities move large groups of particles from lattice sites by advective displacements and diffusive jumps, similarly to a PT scheme, with the difference that the particles undergo discrete movements between lattice sites, performed with a global procedure. BGRW schemes use random walk probabilities biased in the direction of the local velocity which simulate both the advection and the diffusion displacements through collective jumps on first-neighbor lattice sites. The two schemes can be used alternatively in solving specific transport problems. The unbiased GRW, which is significantly faster than the BGRW algorithm, is more efficient in solving large scale problems in unbounded domains but less efficient in solving boundary value problems. The BGRW instead is more accurate, allows straightforward implementations of the boundary conditions, and is more appropriate to solve coupled problems for soil systems with complex geometry.

Since the numbers of particles at lattice sites are readily converted into concentrations, we propose an operator splitting solution for the reaction step consisting of directly approximating the deterministic reaction terms and adding their contributions to the concentrations resulted from the advection-dispersion transport step. For saturated flow conditions the reactive transport solver can be implemented as a non-iterative scheme, with an explicit linearization of the reaction terms using concentration values provided by the transport step. Code verification tests for reactive transport with Monod reactions indicate the convergence towards manufactured analytical solutions of order one of the non-iterative GRW scheme and of order two of the BGRW scheme. For unsaturated flow conditions, the non-iterative schemes are no longer accurate and fail to converge (see Remark 2). Therefore, we designed iterative BGRW \(L\)-schemes coupled with GRW \(L\)-schemes for Richards equation, which can be easily adapted to bidirectional coupling between the flow and transport equations. The coupled \(L\)-schemes provide accurate solutions in case of variable saturation, as well as in case of degenerate Richards equation which describes transitions between unsaturated and saturated regimes. Code verification tests for transport and biodegradation indicate the approach to the convergence order 2 for pressure head and species concentrations solutions and a slower convergence of order 1.5 for the velocity components.

We applied the non-iterative GRW and BGRW schemes, which are faster than the iterative BGRW \(L\)-scheme, to solve a benchmark problem for biodegradation controlled by the transverse mixing in heterogeneous aquifers. The BGRW solution, for parameters describing a situation with moderate heterogeneity, served as reference to evaluate the faster but less accurate GRW solution. The comparison shows that the GRW solution reproduces the shape and the spatial distribution of electron donor, electron acceptor, and biomass concentrations but is affected by overshooting errors leading to the overestimation or the biodegradation process. Further investigations with the GRW solver show that the overshooting is enhanced by increasing the correlation lengths and the variance of the log-hydraulic conductivity. Since unbiased GRW and PT procedures move particles according to numerical schemes for Itô equation over spatial regions with varying properties (e.g., reaction fronts), the overshooting can be only hardly avoided at the cost of using very small time steps (see Remark 3). The occurrence of the overshooting errors in solving biodegradation problems remains a serious issue of concern for both unbiased GRW and PT approaches.

Finally, we solved the more challenging problem consisting of flow in soils, with transition from unsaturated to saturated conditions, and reactive transport with bidegradation governed by a double Monod model by using the flow GRW \(L\)-scheme coupled with the accurate BGRW \(L\)-scheme. The numerical solutions illustrate the progress of the saturation and biodegradation process in heterogeneous soils with different hydraulic properties and demonstrate the capability of the approach to simulate bioremediation scenarios.

Acknowledgements

The first author acknowledges the financial support of the under Grant SU 415/4-1 – 405338726 “Integrated global random walk model for reactive transport in groundwater adapted to measurement spatio-temporal scales”.

Alecsa et al., 2020

Alecsa, C.D., Boros, I., Frank, F., Knabner, P., Nechita, M., Prechtel, A., Rupp, A., Suciu, N., 2020. Numerical benchmark study for fow in heterogeneous aquifers. Adv. Water Resour., 138, 103558. https://doi.org/10.1016/j.advwatres.2020.103558

Bause and Knabner, 2004

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

Benson and Meerschaert, 2008

Benson, D.A., Meerschaert, M.M., 2008. Simulation of chemical reaction via particle tracking: Diffusion-limited versus thermodynamic rate-limited regimes. Water Resour. Res. 44(12):W12201. http://dx.doi.org/10.1029/2008WR007111

Benson and Bolster, 2016

Benson, D.A., Bolster, D., 2016. Arbitrarily complex chemical reactions on particles. Water Resour. Res., 52(11), 9190–9200. https://doi.org/10.1002/2016WR019368

Bolster et al., 2016

Bolster, D., Paster, A., Benson, D.A., 2016. A particle number conserving L agrangian method for mixing-driven reactive transport. Water Resour. Res., 52(2), 1518–1527. https://doi.org/10.1002/2015WR018310

Engdahl et al., 2017

Engdahl, N.B., Benson, D.A., Bolster, D., 2017. Lagrangian simulation of mixing and reactions in complex geochemical systems. Water Resour. Res., 53(4), pp.3513–3522. https://doi.org/10.1002/2017WR020362

Brunner et al., 2012

Brunner, F., Radu, F.A., Bause, M., Knabner, P., 2012. Optimal order convergence of a modified BDM1 mixed finite element scheme for reactive transport in porous media. Adv. Water Resour. 35, 163–171. http://dx.doi.org/10.1016/j.advwatres.2011.10.001

Cirpka et al., 1999

Cirpka, O.A., Frind, E.O., Helmig, R., 1999. Numerical simulation of biodegradation controlled by transverse mixing. J. Contam. Hydrol., 40(2), 159-182. https://doi.org/10.1016/S0169-7722(99)00044-3

Cui et al., 2014

Cui, Z., Welty, C., Maxwell, R.M., 2014. Modeling nitrogen transport and transformation in aquifers using a particle-tracking approach. Comput. Geosci., 70, 1–14. https://doi.org/10.1016/j.cageo.2014.05.005

Ding and Benson, 2015

Ding, D., Benson, D.A., 2015. Simulating biodegradation under mixing-limited conditions using Michaelis-Menten (Monod) kinetic expressions in a particle tracking model. Adv. Water Resour., 76, 109–119. https://doi.org/10.1016/j.advwatres.2014.12.007

Ding et al., 2017

Ding, D., Benson, D.A., Fernàndez-Garcia, D., Henri, C.V., Hyndman, D.W., Phanikumar, M.S., Bolster, D. (2017). Elimination of the reaction rate “scale effect”: Application of the Lagrangian reactive particle-tracking method to simulate mixing-limited, field-scale biodegradation at the Schoolcraft (MI, USA) Site. Water Resour. Res., 53(12), 10411–10432. https://doi.org/10.1002/2017WR021103

Eberhard et al., 2007

Eberhard, J.P., Suciu, N., Vamoş, C., 2007. On the self-averaging of dispersion for transport in quasi-periodic random media. J. Phys. Math. Theor., 40(4), 597–610. https://doi.org/10.1088/1751-8113/40/4/002

Edery et al., 2010

Edery Y., Scher, H., Berkowitz, B., 2010. Particle tracking model of bimolecular reactive transport in porous media. Water Resour. Res. 46(7):W07524. http://dx.doi.org/10.1029/2009WR009017

Fernàndez-Garcia and Sanchez-Vila, 2011

Fernàndez-Garcia, D., Sanchez-Vila, X., 2011. Optimal reconstruction of concentrations, gradients and reaction rates from particle distributions. J. Contam. Hydrol., 120, 99–114. https://doi.org/10.1016/j.jconhyd.2010.05.001

Haverkamp et al., 1977

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

Henri and Fernandez-Garcia, 2015

Henri, C.V., Fernàndez-Garcia, D., 2015. A random walk solution for modeling solute transport with network reactions and multi-rate mass transfer in heterogeneous systems: Impact of biofilms. Adv. Water Resour., 86, 119–132. https://doi.org/10.1016/j.advwatres.2015.09.028

Higham and Kloeden, 2021

Higham, D.J., Kloeden, P.E., 2021. An Introduction to the Numerical Simulation of Stochastic Differential Equations. SIAM. https://my.siam.org/Store/Product/viewproduct/?ProductId=32922725

Illiano et al., 2020

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

Izsák and Lagzi, 2010

Izsák, F., Lagzi, I., 2010. Models of Liesegang pattern formation. In: Precipitation Patterns in Reaction-Diffusion Systems. Kerala, India: Research Signpost, p. 207–217. http://bolyai.cs.elte.hu/ izsakf/english/proof_marc26.pdf

Karapiperis and Blankleider, 1994

Karapiperis, T., Blankleider, B., 1994. Cellular automaton model of reaction-transport processes. Physica D 78, 30–64. http://dx.doi.org/10.1016/0167- 2789(94)00093-X

Karapiperis, 1995

Karapiperis, T., 1995. Cellular automaton model of precipitation/dissolution coupled with solute transport. J. Stat. Phys. 81(1-2), 165–180. http://dx.doi.org/10.1007/BF02179974

Klöfkorn et al., 2002

Klöfkorn, R., Kröner, D., Ohlberger, M., 2002. Local adaptive methods for convection dominated problems. Int. J. Numer. Meth. Fluid., 40(1-2), 79–91. https://doi.org/10.1002/fld.268

Knabner et al., 2003

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

Kreutle and Knabner, 2007

Kräutle S, Knabner P., 2007. A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions. Water Resour. Res. 41(9), W09414. https://doi.org/10.1029/2005WR004465.

Lie, 2019

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

List and Radu, 2016

List, F. and 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

Nagy and Izsák, 2011

Nagy, N., Izsák, F., 2011. Stability of reaction fronts in random walk simulations. Applied Mathematics Research eXpress 2011:abr016. http://dx.doi.org/10.1093/amrx/abr016.

Natarajan et al., 2021

Natarajan, H., Popov, P.P., Jacobs, G.B., 2021. A high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangian drift-diffusion models coupled with Eulerian discontinuous spectral element method. Comput. Meth. Appl. Mech. Eng., 384, 114001. https://doi.org/10.1016/j.cma.2021.114001

Pop et al., 2004

Pop, I.S., Radu, F.A., Knabner, P., 2004. Mixed finite elements for the Richards’ equation: linearization procedure. Journal of computational and applied mathematics 168 (1-2), 365–373. https://doi.org/10.1016/j.cam.2003.04.008

Pope, 1985

Pope, S.B., 1985. PDF methods for turbulent reactive flows. Prog. Energy. Combust. Sci., 11(2), 119–192. https://doi.org/10.1016/0360-1285(85)90002-4

Radu, 2004

Radu, F.A., 2004. Mixed finite element discretization of Richards’ equation: error analysis and application to realistic infiltration problems. PhD Thesis, University of Erlangen-Nürnberg (2004)

Radu et al., 2008

Radu, F.A., Bause, M., Prechtel, A., Attinger, S., 2008. A mixed hybrid finite element discretization scheme for reactive transport in porous media. In Numerical mathematics and advanced applications (pp. 513-520). Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-540-69777-0_61

Radu et al., 2009

Radu, F.A., Pop, I.S., Attinger, S., 2010. Analysis of an Euler implicit-mixed finite element scheme for reactive solute transport in porous media. Numer. Meth. Part. Differ. Equat. 26(2), 320–344. https://doi.org/10.1002/num.20436

Radu et al., 2011

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 (1), 47–61. https://doi.org/10.1016/j.advwatres.2010.09.012

Sanchez-Vila and Fernández-Garcia, 2016

Sanchez-Vila, X., Fernández-Garcia, D., 2016. Debates-Stochastic subsurface hydrology from theory to practice: Why stochastic modeling has not yet permeated into practitioners?. Water Resour. Res., 52(12), 9246–9258. https://doi.org/10.1002/2016WR019302

Schmidt et al., 2018

Schmidt, M.J., Pankavich, S.D., Benson, D.A., 2018. On the accuracy of simulating mixing by random-walk particle-based mass-transfer algorithms. Adv. in Water Resour., 117, pp.115–119. https://doi.org/10.1016/j.advwatres.2018.05.003

Schmidt et al., 2019

Schmidt, M.J., Pankavich, S.D., Navarre-Sitchler, A., and Benson, D.A., 2019. A Lagrangian method for reactive transport with solid/aqueous chemical phase interaction. J. Comput. Phys., X, 2, 100021. https://doi.org/10.1016/j.jcpx.2019.100021

Schwarze et al, 2001

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

Sole-Mari and Fernàndez-Garcia, 2018

Sole-Mari, G., Fernàndez-Garcia, D., 2018, Lagrangian modeling of reactive transport in heterogeneous porous media with an automatic locally adaptive particle support volume, Water Resour. Res. 54(10), 8309–8331, http://dx.doi.org/10.1029/2018WR023033

Sole-Mari et al., 2019

Sole-Mari, G., Bolster, D., Fernàndez-Garcia, D., Sanchez-Vila, X., 2019. Particle density estimation with grid-projected and boundary-corrected adaptive kernels. Adv. Water Resour. 131, 103382. https://doi.org/10.1016/j.advwatres.2019.103382

Sole-Mari et al., 2020

Sole-Mari, G., Fernàndez-Garcia, D., Sanchez-Vila, X., Bolster, D., 2020. Lagrangian modeling of mixing-limited reactive transport in porous media: Multirate interaction by exchange with the mean. Water Resour. Res., 56(8), e2019WR026993. https://doi.org/10.1029/2019WR026993

Sole-Mari et al., 2021

Sole-Mari, G., Schmidt, M.J., Bolster, Fernàndez-Garcia, D., 2021. Random-walk modeling of reactive transport in porous media with a reduced-order chemical basis of conservative components. Water Resour. Res., 57, e2020WR028679. https://doi.org/10.1029/2020WR028679

Suciu, 2019

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

Suciu, 2020

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

Suciu et al., 2020

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

Suciu et al., 2021

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

Suciu and Radu, 2021

Suciu, N., Radu, F.A., 2021. MonodReactions. Git repository. https://doi.org/10.5281/zenodo.5104076

Vamoş et al., 2003

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.

Wiedemeier et al., 1999

Wiedemeier, T.H., Rifai, H.S., Newell, C.J., Wilson, J.T., 1999. Natural attenuation of fuels and chlorinated solvents in the subsurface. John Wiley & Sons.