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 [2, 4] or to inappropriate parameterizations of the dispersion coefficients , leads to an artificial mixing of the species and to overpredictions of the simulated biodegradation process. 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 . Even more challenges may appear in 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 [1, 37].
Most of the numerical approaches in modeling reactive transport in porous media are based on finite element [19, 4, 28] or finite volume [5, 22, 14] 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 [20, 21]. The numerical diffusion in advection-dominated transport problems is controlled by using higher-order schemes and adaptive time steps , mixed-finite element methods using Lagrange multipliers , or streamtube approaches .
An attractive alternative to solve advection-dominated problems is provided by various particles methods, which by construction are free of numerical diffusion. Particle approaches to reactive transport were developed with sequential particle tracking (PT)  and continuous time random walk (CTRW)  schemes, as well as with sequential biased random walk cellular automata (CA) [16, 17] and global procedures [15, 24] based on global random walk (GRW) algorithms .
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 , or by using the probability that the two particles will occupy the same position, constructed as a convolution of two Gaussian probability densities . The latter approach was used to simulate biodegradation processes in a column experiment  and in a real field application 
. In these simulations, the concentrations estimated by simple binning of particle numbers were further smoothed by averaging over PT runs.
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[11, 34]. In case of advection-dominated problems, the effect of the numerical diffusion induced by spatial redistribution of particles into concentrations in control volumes [5, 6, 31] 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 . Kernel density estimations have also been used in PT modeling with reduction schemes  similar to those introduced in finite element approaches for multicomponent reactive systems (e.g. ).
Besides the need to map particle densities to concentrations, the coupling between transport and nonlinear reactions is another issue of concern in PT approaches . 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 . 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 . 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 to 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 [15, 17]. 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 and the rate of reaction production are proportional to the physical rate constant , e.g. and , for the reaction . The CA reaction rules use functions of occupation numbers and and reaction probabilities to decide whether a reaction takes place if an -particle and two -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 . The species concentration is retrieved by multiplying the particle densities at lattice sites by appropriate scaling factors.
In GRW-type procedures presented in [15, 24], the reaction step is performed deterministically. Comparing with the PT approach with deterministic reaction step  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 , , is directly converted into the amount of -particles by representing the species concentrations as numbers of particles at lattice sites normalized by the initial number of moles of each species multiplied by the number of particles used to represent a mole .
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 which becomes expensive in terms of computing time, with a slow linear increase for numbers of particles larger than [41, Fig. 9]. Therefore, a limited number of particles was used in these diffusion-reaction simulations (e.g. for each molecular species in ) 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 , also used in our recent papers [37, 38]. In this way, one achieves fast and highly accurate simulations, which cannot be done with the sequential PT and biased random walk CA approaches described above. 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 , i.e. there is no need to perform averages over runs [7, 8] or kernel density estimations [11, 34] 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 , where the treatment of the nonlinear terms in the Richards equation is based on the -scheme linearization approach [26, 23]. 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 number 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 on 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[36, 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  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 coefficient 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, 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 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 -schemes introduced in [38, 39].
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 Mood 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 .
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 . 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 . 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 .
We consider a two-dimensional setting with an open domain , with boundary , and a time interval . Taking into account both the saturated and the vadose zone, the water flux is given by Darcy’s law
where is the height against the gravitational direction. The pressure head solves the Richards equation
where is the water content and stands for the hydraulic conductivity of the medium. Equation (2) is closed by material laws expressing dependencies between , , and , 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 , is governed by the following system of coupled nonlinear equations
The constant in Eq. (4) is the decay rate of the immobile component. The factor of the reaction rate may be used to model the growth limitation of (due for instance to limited pore space available), is the maximum immobile concentration, and the coefficient is set to one in the model with growth limitation and to zero otherwise [27, 2]. In a basic scenario for aerobic biodegradation [2, 4, 5, 27], represents the concentration of the electron donor, corresponds to the electron acceptor (oxygen), 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 , is the microbial yield coefficient per unit electron donor consumed, and the Monod term is defined by
where is the maximum growth rate and , are Monod constants. The reaction terms for the mobile species in Eqs. (3) are expressed as and , where , are stoichiometric constants.
As compared with more complex biodegradation models [2, 6, 27], 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 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 [14, 19, 39] after replacing the general constitutive law for the nonlinear water content by the simpler law , 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 modification in the scheme for fully coupled flow and transport from  one obtains an iterative GRW solver for multicomponent reactive transport. We note that the operator splitting PT method for reactive transport in variably saturated media proposed in 
also uses an iterative Newton method to solve the nonlinear system of ordinary differential equations of chemical kinetics. In the saturated flow regime,becomes a constant parameter and, for given flow velocity, the GRW solution can be constructed, similarly to the solution for reaction-diffusion problems from , 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 ).
3.1 Non-iterative schemes for reactive transport in saturated media
3.1.1 Biased GRW algorithm for the transport step
The transport equation (3) for one mobile concentration , with a source/sink term independent of added in the right-hand side, can be equivalently written as
The new term 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 is computed with the GRW algorithm for flows in saturated/unsaturated porous media presented in [39, Sect. 4.1], with velocity on the boundaries extended from near boundary interior sites or approximated by forward finite differences, for the particular case of constant . The transport step, described by Eq. (6) without the reaction term, can be computed with explicit GRW algorithms. In the following, a two-dimensional BGRW which simulates the advection through biased jump probabilities [16, 36]
is derived from the forward-time centered-space finite difference approximation of the transport equation. We assume a diagonal diffusion tensor with constant componentsand and by and we denote the components of the Darcy velocity along the horizontal axis and the vertical axis , by the time step, and by and the spatial steps. With these, the explicit scheme reads
As usual in GRW algorithms, we simulate the random walk of computational particles on the regular lattice to approximate the concentration by the density of the number of particles, . By using the finite difference scheme (3.1.1) and the dimensionless parameters
The contributions to in Eq. (3.1.1) are obtained with the BGRW algorithm
where, for consistency with the finite difference scheme (3.1.1), the quantities verify in the mean
The binomial random variables used in the BGRW algorithm are approximated by using the unaveraged relations (3.1.1), summing up the reminders of multiplication by and of the floor function , and allocating one particle to the lattice site where the sum reaches the unity (see also [36, 38]). By giving up the particle’s indivisibility and using the un-averaged relations (3.1.1) one obtains a deterministic BGRW algorithm. However, throughout the present study we shall use only the randomized implementation of the BGRW algorithm.
As follows from (3.1.1), the BGRW algorithm is subject to the following restrictions
3.1.2 Unbiased GRW algorithm for the transport step
The unbiased two-dimensional GRW algorithm for the transport step in Eq. (6) is obtained by globally moving groups of particles according to the rule
where is a constant amplitude of diffusion jumps and the dimensionless variables , , and are defined similarly to (8) by
The particles distribution is updated at every time step by
The averages over GRW runs of the terms from (13) are now related by
The binomial random variables 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 is chosen such that the dimensionless parameters and take integer values larger than unity which ensure the desired resolution of the velocity components [36, Sect. 220.127.116.11]. Further, the jumps’ amplitude is chosen such that the jump probabilities verify the constraint , imposed by the first relation (16).
3.1.3 Deterministic reaction step
The concentrations of the mobile molecular species are computed as , where is the distribution of the computational particles after the execution of the transport step obtained with either the BGRW algorithm (3.1.1) or the unbiased GRW algorithm (15). The total number of particles of the -species is chosen as precisely given by the initial number of -molecules in the reaction system, i.e., , where is the initial number of moles and is the Avogadro number. 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 ) according to
where 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 , 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 (14) [39, Sect. 5.2.3] we shall consider in the following only the biased algorithm. The iterative BGRW -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 added in the right-hand side, by adding the stabilization term , where is a dimensionless constant parameter and , , is the iteration index. Following the same procedure as in Section 3.1.1 one obtains
where , with pressure head provided by the flow solver (see also [39, Sect. 4.2.1]).
The parameters of the BGRW -scheme are defined by
The BGRW -scheme is completed by the same relations (10)-(12) as in case of non-iterative BGRW scheme from Section 3.1.1, which now have to be verified at every iteration . The iterations start with and are stopped when the discrete norms of the numerical solution verify
for some given tolerances and .
This BGRW -scheme is implemented into an operator splitting procedure alternating nonreactive transport for all mobile species and reaction steps according to Eq. (17) at every iteration level .
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. , 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  or fully coupled one-dimensional flow and transport with transition from unsaturated to saturated conditions . Here, we adopt this approach to test the new schemes for problems of increased complexity: constant flow and nonlinear reactions, saturated flows coupled with Monod model, and space-time variable flows with transition from unsaturated to saturated regime coupled with Monod biodegradation 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 , , , , and reaction terms given by and . The problem is solved in the domain and the final simulation time is . The manufactured analytical solutions
Numerical solutions obtained with the non-iterative BGRW algorithm presented in Section 3.1.1 are computed on regular lattices successively refined by halving five times the initial step sizes . The accuracy of the numerical solutions at the final time is quantified by discrete norms and the estimated order of convergence (EOC) is given by the slope of decreasing error norms in logarithmic scale [1, 39].
constant velocity and nonlinear reactive transport.
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 -variable were removed. In this case, it is found that the numerical solutions converge with the order 2 as well, but the errors are generally one order of magnitude larger (similarly to the one-dimensional solutions for the coupled flow and transport problem presented in [38, Tables 10 and 11]).
For the coarsest lattice used in this test () the local Péclet number takes the maximum allowable value (see Section 3.1.1). Hence, as decreases , the relations (12) are verified, and the BGRW algorithm is free of numerical diffusion [39, Table 9]. For smaller the advection dominates and has to be decreased to fulfil the condition , which would slow down the BGRW computations. The alternative is to use the unbiased GRW solver presented Section 3.1.2 which is not constrained by a maximum Pé. By varying the dimensionless velocity and the amplitude of diffusion jumps
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 forand , if and there are groups of particles jumping across the boundary from interior sites (e.g., for 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 [38, Sect. 6.3], which renders the unbiased GRW very efficient in large scale simulations of transport in groundwater .
The smallest errors of the GRW solutions are obviously obtained if and . For the setup of the present test problem this choice is not possible because if the amplitude of the diffusion jumps in (14) should be at least to comply with the restriction , for , and even larger is needed for smaller . Instead, for strongly advection-dominated problems with , as those considered in the next section, the choice is possible for all considered in Table 1. The results of the biased GRW solutions, with the new dispersion coefficient and letting the other parameters unchanged, are presented in Table 2.
for constant velocity and nonlinear reactive transport.
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 border 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 , constant hydraulic conductivity , and constant dispersion coefficients . The nonlinear reaction is governed by the double Monod model (5) with parameters , , , stoichiometric coefficients , , analytical solutions for the election donor and election acceptor concentrations verifying Eqs. (3)-(4) given by (21)-(22), and constant biomass concentration .
The coupled system of equation (2)-(3) with source terms in the right-hand side, computed analytically by differentiating the solutions (21)-(23), is solved with the flow GRW -scheme from [39, Sect. 4.1] coupled with the non-iterative BGRW scheme introduced in Section 3.1.1. The flow solutions are approximated, with parameter and tolerances , by performing a large number of iterations of the flow scheme, lying between 750 and 4500 for different time points and lattice steps . The errors and the orders of convergence, estimated as in the previous section by halving the space step, are presented in Table 3.
flow and Monod reactions.
For comparison, the same problem is solved with the iterative BGRW -scheme from Section 3.2. For stabilization parameter and tolerances , the convergence of the liniarization procedure for different values is already established after about 20 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.
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 (23), according to (1). The results obtained with the non-iterative BGRW scheme, presented in Table 5, show a much improved order of convergence. The smaller EOC values of the concentration solutions shown in case Table 3, can be certainly attributed to the moderate accuracy of the velocity provided by 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.
coupled saturated flow and Monod reactions.
analytical flow velocity and Monod reactions.
4.3 Saturated/unsaturated flows and Monod reactions
To model both unsaturated and saturated flow conditions, we consider the pressure head solution
the constant hydraulic conductivity , and a variable water content described by
For the beginning, we investigate the convergence behavior of the BGRW -scheme in conditions of strictly unsaturated flow by removing the terms from (24). With the same parameters , , and 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 orders of convergence shown in Table 6 are a bit larger than in saturated case (Table 4). The results can be in principle further improved by increasing the stabilization parameter [39, Sect. 5]. 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 tens of seconds, this difference is caused by the slower convergence of the flow solver in the saturated case.
The non-iterative BGRW scheme from Section 3.1 can be adapted for variable water content by adding to the right-hand side of Eq. (3.1.1) the source term . Performing a numerical test similar to that presented in Table 6 leads to relative errors of the numerical solutions for mobile concentrations and of about , even if the spatial step is further decreased from to . For comparison, the relative errors of the solutions provided by the non-iterative scheme for in the convergence test for saturated flow presented in Table 3 are about for and for . Therefore, while the non-iterative scheme is accurate and several times faster than the iterative BGRW scheme in case of saturated flows, it is not appropriate in problems for unsaturated flow conditions.
coupled unsaturated flow and Monod reactions.
Richards equation describing transitions from unsaturated to saturated flow regimes.
BGRW -scheme coupled with the flow solver from 
for saturated/unsaturated flow and Monod reactions.
The challenging situation where Richards equation degenerates into an elliptic equation when takes constant values in saturated regions has been considered in  for one-dimensional transport of a single molecular species. In the following, we test the one-way coupling of the two-dimensional flow and reactive transport -schemes by comparisons with the exact concentration solutions (21)-(22) and the flow solutions (24)-(25), which captures the transition from unsaturated and saturated regime, within the setup of Monod reactions 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 -schemes. With the same parameters , , and as used above, the desired accuracy of the iteration procedure is attained after 150 to 2200 iterations of the flow scheme and 25 to 35 iterations of the transport scheme. The computation of the coupled problem for the four discretization levels considered in this test takes about 1.17 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 concentration solutions. The error norms of the velocity solutions are smaller by one order of magnitude than those for pressure and concentration. However, their convergence order lies between 1.4 and 1.5. This could partially explain the decrease of the order of convergence of the concentration solutions shown in Table 8.
We note that the influence of the stabilization parameter on accuracy and order of convergence estimations (cf. [39, Table 4]) has not been investigated in the present study. Also, we did not tested the unbiased GRW non-iterative scheme presented in Section 3.1.2 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 equals the porosity and enters as a constant in the BGRW and GRW parameters given by relations (8) and (14), 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  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 ) or troubles with the construction of efficient iterative solvers (when using quadratic finite element schemes and streamline upwind stabilization ).
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 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 that the iterative schemes.
The problem is formulated in the domain . The reactive transport of the mobile species is governed by Eqs. (3) with parameters specified by , , , . The Monod reaction term is modeled by Eq. (5) with , the biomass is kept constant, , and the reaction terms in Eqs. (3) are given by , with and . The final simulation time is , when one expects that the process reaches the steady state . As initial conditions, one chooses and in . The contaminant (electron donor) is injected at the middle part of the inflow boundary, where no electron acceptor (oxygen) is present, i.e., and for . Further, one sets and for . 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 ( and for ) which would render the computation until the final time prohibitive. Therefore, we first restrict the final time to and compare the BGRW solution to that provided by the unbiased GRW solver. The latter is obtained with (), , and GRW parameters , . The comparison is done in terms of concentrations averaged over lattice sites, , by relative errors , where is the BGRW solution. The errors obtained, and , validate the GRW solver for the parameters of the benchmark problem described above and . Further, the problem is solved for with the unbiased GRW solver. The electron donor and acceptor concentrations at the final time shown in Fig. 2 are similar to those obtained by quadratic finite element in [2, 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 . 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 , and the more intense biodegradation in the lower part of the domain leads to the underestimation of the contaminant concentration 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 . Classical approaches to cope with the numerical diffusion in advection-dominated problems are based on streamtube reactive transport models , local mesh refinement , or higher-order finite element used together with stabilization techniques . 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 , with a two-dimensional domain , 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 [36, 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 . The log-hydraulic conductivity field
is a statistically homogeneous random function with isotropic Gaussian correlation. Its variance belongs to the setand for its correlation length we consider and . The mean conductivity is set to . With a decrease of of the pressure head between the inflow and outflow boundaries and a gradient , the mean velocity becomes . The dispersion is parameterized by a constant coefficient . Ensemble averages are estimated over 256 realizations of the field. These parameters provide fairly good approximations of the random Darcy velocity field and of the stochastic advection-dispersion process [32, 9]. Since the porosity does not play a major role in saturated flows it is set to .
The setup of the chemical reactions closely follows that used in [5, 2], 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 , , . Additionally, following  we consider a growth limitation in the Monod model with and in Eq. (4). The initial concentration of the contaminant is set to at lattice sites inside a square of side equal to 2 centered at and and to otherwise. The oxygen concentration is initially set to in , excepting the injection region where it is . The continuous injection is simulated, as in , by resetting after every time step both and to their initial values in the injection region. The initial concentration of the biomass is set to in . In , the boundary conditions for the mobile species are set to the initial values of and 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 close to boundaries we also fix 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 both the GRW and BGRW algorithms simulate the physical outflow situation by removing the particles crossing the boundary.
In Fig. 4 we present the concentration profiles of the three species at the final time days computed with the non-iterative BGRW solver for a fixed realization of the velocity field with and . The results are obtained on a lattice with which ensures a maximum Péclet number Pé=1.75 verifying the relations (12). The corresponding time step, according to (10) is . The BGRW solution is further used as reference for comparison with that provided by the non-iterative biased GRW solver for the same realization of the velocity field. The latter is obtained by using a much coarser discretization, with (maximum Pé=6.57) and corresponding to and in (14). Withe these, the GRW solver achieves a speedup of over the BGRW solver.
The GRW solutions shown in Fig. 4 reproduce the general trend of the concentration profiles of the BGRW solution (Fig. 4) 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. 4. The more extended profiles of biomass concentration 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 , , and , that is, one order of magnitude larger than in case of 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 and which give rise to enhanced overshooting. The negative value 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 field with , and the corresponding variances are shown in Figs. 6 and 6. The averaged concentrations shown in Fig. 6 are similar to the results obtained with uniform velocity. But the relative errors of the lattice-averaged mean concentrations from Fig. 6 with respect to the BGRW solutions for a constant longitudinal velocity set to are close to in the comparison of the GRW/BGRW solutions for a single realization of the field presented above. The large variance in Fig. 6, consistent with the inhomogeneity of the electron acceptor concentration 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 and variances of the field, computed at the final time , are presented in Table 9. The mean species concentrations show some variation, of about 22% (), 25% (), and 12% (), but there is no systematic dependence on and . The variances slightly increase with but do not show a systematic dependence on .
and variances averaged over lattice sites.
The comparison of , , and for the first realization in the Monte Carlo ensembles computed for the six combinations of and , presented in Fig. 9, 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. 9, for fixed the meandering is more pronounced as the variance increases. Comparing the two columns, we remark an increase in meandering amplitude as the correlation length of the filed increases. We also remark that the amount of overshooting indicated by zones of low concentrations of the electron acceptor outside the reaction fronts, increases with both and .
We note that the overshooting is an issue of concern for both unbiased GRW and particle tracking approaches. The GRW rule (13) moves particles in space, between the lattice sites, by summing discrete advective displacements and diffusive jumps, e.g., , in the same way as a numerical solution of the Itô equation, which is the mathematical model of the PT procedure. The only differences are, first, that and are discrete quantities which measure displacements in units of the spatial steps of the lattice, and second, instead of computing ensembles of PT trajectories, the GRW scheme counts the number of particles at lattice sites which approximate the solution of the partial differential equation to be solved. In both schemes, displacements of computational particles over space points where the problem data may have significant variations, e.g., diffusion fronts, are hardly avoidable. For instance, in the semi-Lagrangian method proposed by , to prevent a tracked particle from leaving the computational subdomain, the time step is restricted such that the sum of advective and dispersive displacements is smaller than the minimum grid spacing. Further, the polynomial solution of the transport problem associated to the particle has to be interpolated to the nodes of the Chebyshev quadrature where a coupled spectral element method provides the advection velocity. For comparison, in the BGRW method, which also restricts the time step, the smooth solution is represented by a sufficiently large number of particles which jump only to first neighbor lattice sites and the velocity is given at the same lattice sites (by another GRW solver in coupled flow and transport problems). One avoids in this way both the overshooting and the interpolation procedure.
are discrete quantities which measure displacements in units of the spatial steps of the lattice, and second, instead of computing ensembles of PT trajectories, the GRW scheme counts the number of particles at lattice sites which approximate the solution of the partial differential equation to be solved. In both schemes, displacements of computational particles over space points where the problem data may have significant variations, e.g., diffusion fronts, are hardly avoidable. For instance, in the semi-Lagrangian method proposed by
, to prevent a tracked particle from leaving the computational subdomain, the time step is restricted such that the sum of advective and dispersive displacements is smaller than the minimum grid spacing. Further, the polynomial solution of the transport problem associated to the particle has to be interpolated to the nodes of the Chebyshev quadrature where a coupled spectral element method provides the advection velocity. For comparison, in the BGRW method, which also restricts the time step, the smooth solution is represented by a sufficiently large number of particles which jump only to first neighbor lattice sites and the velocity is given at the same lattice sites (by another GRW solver in coupled flow and transport problems). One avoids in this way both the overshooting and the interpolation procedure.
6 Biodegradation in variably saturates 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 used in verification tests presented in Section 4. The infiltration is driven by a variable pressure, given by and , on the open part of the top boundary and by communication with the fresh water at hydrostatic equilibrium on the open parts of the domain and . The zero flow condition is imposed on the closed boundary and the initial pressure corresponds to hydrostatic equilibrium, . Hereafter, the time is measured in days and the distances in meters. The intermediate time for the variable pressure condition is set to and the simulations are conducted for increasing times .
The relationships defining the water content and the hydraulic conductivity are given by the van Genuchten-Mualem model
where is the normalized water content, and are the residual and the saturated water content, is the hydraulic conductivity of the saturated soil, and , , 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 [39, Sect. 5].
|Silt loam||Beit Netofa clay|
The soil heterogeneity is modeled by a fixed realization of a random saturated hydraulic conductivity with mean . The latter is computed by the Kraichnan method described in Section 5.2, with an anisotropic Gaussian model of correlation lengths , and variance . The dispersion is parameterized by the constant coefficient .
At the initial time, the concentration of benzene (electron donor) is set to mole on the injection boundary and to otherwise, the concentration of oxygen (electron acceptor) is on and mole otherwise, and the biomass concentration is mole in . At , the concentrations