Lagrangian simulation of multi-step and rate-limited chemical reactions in multi-dimensional porous media

Bing-qing Lu,Yong Zhng,*,Hong-gung Sun,Chun-mio Zheng

aDepartment of Geological Sciences,University of Alabama,Tuscaloosa,AL 35487,USA

bCollege of Mechanics and Materials,Hohai University,Nanjing 210098,China

cSchool of Environmental Science and Engineering,Southern University of Science and Technology,Shenzhen 518055,China

Abstract Management of groundwater resources and remediation of groundwater pollution require reliable quantification of contaminant dynamics in natural aquifers,which can involve complex chemical dynamics and challenge traditional modeling approaches.The kinetics of chemical reactions in groundwater are well known to be controlled by medium heterogeneity and reactant mixing,motivating the development of particle-based Lagrangian approaches.Previous Lagrangian solvers have been limited to fundamental bimolecular reactions in typically one-dimensional porous media.In contrast to other existing studies,this study developed a fully Lagrangian framework,which was used to simulate diffusion-controlled,multi-step reactions in one-,two-,and three-dimensional porous media.The interaction radius of a reactant molecule,which controls the probability of reaction,was derived by the agent-based approach for both irreversible and reversible reactions.A flexible particle tracking scheme was then developed to build trajectories for particles undergoing mixing-limited,multi-step reactions.The simulated particle dynamics were checked against the kinetics for diffusion-controlled reactions and thermodynamic wellmixed reactions in one-and two-dimensional domains.Applicability of the novel simulator was further tested by(1)simulating precipitation of calcium carbonate minerals in a two-dimensional medium,and(2)quantifying multi-step chemical reactions observed in the laboratory.The flexibility of the Lagrangian simulator allows further re finement to capture complex transport affecting chemical mixing and hence reactions.

Keywords:Lagrangian framework;Chemical reaction;Diffusion-limited process;Multi-step reactions;Interaction radius

1.Introduction

Accurate simulation of reactive transport in natural media remains one of the greatest challenges in hydrology,mainly due to its multi-scale process including molecular-scale reaction and the representative elementary volume(REV)-scale transport(Dentz et al.,2011).Recent efforts have focused on the fundamental,second-order bimolecular reaction A+B⇌C(Willingham et al.,2008;Luo et al.,2008;Benson and Meerschaert,2008;Bolster et al.,2017).To predict the concentration evolution of compounds undergoing bimolecular reaction in laboratory-and regional-scale porous media,the following advection-dispersion-reaction(ADR)continuum model combining transport and reaction at the same REV scale has been widely used(e.g.,Ham et al.,2004):

where v[LT-1]is the mean flow velocity;D[L2T-1]is the macrodispersion coefficient;Kf[M-1L3T-1]and Kr[T-1]are the forward and reverse kinetic coefficients of reaction,respectively;and CA（x，t）,CB（x，t）,and CC（x，t）[ML-3]denote the concentrations of species A,B,and C,respectively,at time t and position x.

The assumption underlying the ADR model(Eq.(1))is the perfect mixing of reactants,which has been found to be invalid for groundwater,since the saturated porous medium is usually inhomogeneous(Tartakovsky et al.,2009;Barnard,2017).The incomplete mixing of reactant particles also explains the wellknown fact that the reaction rate measured with batch experiments(in solutions)is much higher than that observed in sand columns,and the reaction rate measured in the laboratory(under the condition of complete mixing)is typically much higher than that observed in the field(Kapoor et al.,1997;Raje and Kapoor,2000;Cirpka,2002;Gramling et al.,2002),challenging the application of the continuum model(Eq.(1))in quantifying real-world reactive kinetics in aquifers(Cirpka and Valocchi,2016).

Full particle tracking-based Lagrangian approaches have been developed by various researchers for decades to quantify the imperfect mixing of reactants undergoing bimolecular reactions,or the so-called diffusion-controlled bimolecular reaction kinetics in porous media. For example,Smoluchowski(1918)assumed that two reactant molecule particles can react whenever they are close to each other.Gillespie(1977)developed perhaps the first random walkbased stochastic approach to simulate coupled chemical reactions with inherent fluctuations and correlations in spatially homogeneous systems.Toussaint and Wilczek(1983)quantified kinetics of diffusive particles undergoing irreversible annihilation.The pioneering work of Smoluchowski(1918)motivated various stochastic simulation algorithms in the chemical physics and biology communities(Erban and Chapman,2009),including,for example,the compartment models(where molecules within the same compartment can react)(Hattne et al.,2005;Isaacson and Peskin,2006)and the grid-free methods for the simulation of motion of individual molecules(Andrews and Bray,2004;Tournier et al.,2006;Gillespie,2009).Recently,in the hydrological community,Benson and Meerschaert(2008)have developed a probabilitybased scheme by accounting for the overlapped effective volume of two reactant molecules.Edery et al.(2009,2010)introduced a memory function into stochastic simulation algorithms,to capture the non-Fickian transport missed by previous simulations.Ding et al.(2013)developed a particle tracking scheme using two probabilities dictated by the physics of transport and energetics of reaction.Paster et al.(2013)proved theoretically that the particle method proposed by Benson and Meerschaert(2008)does match the governing equations for chemical transport.Benson et al.(2017)compared Lagrangian schemes and Eulerian nonlinear reactive transport in two-dimensional(2D)porous media.Engdahl et al.(2017)then modeled multicomponent reactions in one-dimensional(1D)media.Most of the previous Lagrangian approaches have been limited to bimolecular reactions and/or numerical analysis.To the best of our knowledge,none of the previous Lagrangian solvers has been checked against real world,complex(such as multi-step or higher than second-order)chemical reactions,and most successful applications are limited to 1D media.These limitations motivated this study.

The rest of this paper is organized as follows:Section 2 describes how the Lagrangian framework was developed to approximate reversible,multi-step chemical reactions,including derivation of the spatial range in any dimension where reaction may occur.A fully Lagrangian scheme is then described in section 3,as used to model the coupled advection,dispersion,and reaction at two different scales(i.e.,the molecular-scale reaction and REV-scale transport).In section 4,we check reactive kinetics against classical rate equations in both one and two dimensions,and the application ofthe Lagrangian approximation ofdiffusioncontrolled chemical reactions to simulate a two-step reaction as observed in the literature(Oates and Harvey,2006)and a 2D mineral precipitation process.In section 5,we discuss the ways reaction kinetics change with medium dimensions and possible model extensions.Conclusions are drawn in section 6.

2.Methodology development:Derivation of interaction radius

Fig.1.Interaction volume for 1D,2D,and 3D porous media(arrows show the jump direction of B molecules around an A molecule,and the black circle in(c)denotes the A molecule).

This study conceptualized chemical reaction as a random process,for the following reasons.According to collision theory(Trautz,1916),reactant particles can collide if their nuclei get closer than a certain distance,which is called the interaction radius,denoted as R(Fig.1).Chemical reactions will occur if the colliding particles have enough activation energy and the right orientation to break existing bonds and finish the transformation.External properties,such as the salinity,temperature,and pH,can also affect the reaction rate by affecting the fraction of collisions.Hence,the chemical reaction might be approximated by a probability-controlled,random process.The interaction radius therefore needs to be defined,in order to capture the probabilities of the encounter and combination of reactant particles,which control the reaction rate of incompletely-mixed reactions.

2.1.Interaction radius for complex reversible reactions under arbitrary dimensions

The agent-based approach developed by Pogson et al.(2006)is extended here to define the interaction radius for reversible,multi-step,chemical reactions.The interaction radius also represents the maximum distance for two molecules to collide.Any two molecules separated at a distance larger than R cannot interact during a single time step.In other words,R should be related to the time step.

For illustration purposes,we use the following chain reactions measured by Oates and Harvey(2006):

where Ti denotes the molecular formula of Tiron(1,2-dihydroxybenzene-3,5-disulfonic acid),Mo is the molecular formula of molybdate,MoTi is the 1:1 complex,and MoTi2is the second chelate formed.Eq.(2a)describes the first complexation reaction,whose product(MoTi)then undergoes an additional reaction with Tiron to produce the complex MoTi2expressed by Eq.(2b).Note that reactions in Eq.(2)are reversible.The two-step reaction in Eq.(2)is considered because the reaction kinetics observed by Oates and Harvey(2006)can be used to check the feasibility ofthe Lagrangianschemedevelopedbelow.Themethodology proposed in this paper,including the derivation of the interaction radius,can be conveniently extended for other complex chemical reactions(see discussion in section 4).

For simplicity of description,here we use symbols A,B,C,and D to represent chemicals Ti,Mo,MoTi,and MoTi2,respectively.Eq.(2)can then be rewritten as

It is noteworthy that the rate equation(Eq.(4))is valid for chemical transport in media with any dimension.During a typically small time step Δt[T],the change in the concentration of A is

The proportion of A molecules interacting with B and C at this time step therefore is

If the volume of the medium(with any dimension)is U[Ld](where the exponent d denotes the dimensionality),we obtain the correspondence between the following two dimensionless ratios(Pogson et al.,2006;Zhang et al.,2013,2014):

where UAdenotes the proportion of volume in U where A interacts with B and C.

In a typical,fully Lagrangian solver,the chemical concentration is interchangeable with the corresponding particle number density(LaBolle et al.,1996,2000).Hence,the concentration of A at time t can be calculated by the following equation:

where NA（t）[dimensionless]denotes the number of A particles at time t,and mA[M]is the(molar)mass carried by each A particle.Inserting Eqs.(6)and(8)into(7),one obtains

An A(i.e.,Ti)particle may react with a B(i.e.,Mo)or C(i.e.,MoTi)particle if and only if this B or C particle is located in an interaction volume U*surrounding this A particle of the following size:

Inserting Eq.(9)into Eq.(10)provides the final solution for the interaction volume for reactant A particles at time t,which leads to the interaction radius R varying with the domain dimensionality.For example,for reactions occurring in a 1D sand column(with a total length of L),we have U*=2R(i.e.,the length of a section)(Fig.1(a)),and the rate-limited interaction radius R for reactant A particles is

For a 2D medium with volume U(Fig.1(b)),the interaction radius R is

For a 3D medium with volume U(Fig.1(c)),the interaction radius R is

Eqs.(11)through(13)show three interesting results.First,the reaction kinetics change with time,since the interaction radius R is now a function of time.For a system initially filled with reactants A and B only,the interaction radius R reaches its maximum at the beginning(i.e.,it is dominated by the forward, first complexation reaction),and then gradually decreases in time due to the increased concentration of MoTi and MoTi2(representing an increasing probability for reverse reaction).The value of R eventually approaches zero,representing equilibrium between forward and reverse reactions.Hence,Lagrangian schemes for reactions with the interaction radius changing as in Eqs.(11)through(13)can quantify dynamics of complex reactions.Second,there are four terms on the right hand side(RHS)of Eqs.(11)through(13),representing four procedures that can affect the mass of A.In particular,the first forward reaction A+B→C(i.e.,Ti+Mo→MoTi,captured by the first term on the RHS of Eqs.(11)through(13))and the second forward reaction C+A→D(i.e.,MoTi+Ti→MoTi2,captured by the third term on the RHS of Eqs.(11)through(13))consume A particles,resulting in two positive components in R(that can enhance the reaction rate).Meanwhile,the first reverse reaction A+B←C(see the second term on the RHS of Eqs.(11)through(13))and the second reverse reaction C+A←D(see the fourth term on the RHS of Eqs.(11)through(13))generate A particles,resulting in two negative components in R(that can decrease the reaction rate).If each of the four competing procedures can be incorporated into the particle tracking scheme(instead of using the overall interaction radius R to control reactions for A),it may provide an alternative Lagrangian solver(details of these two options will be discussed in the next section).Third,the interaction radius R does not change linearly with dimensionality(note the different power-law relationships between R and the relevant parameters in Eqs.(11)through(13)),implying the potential scaling behavior when increasing from one dimension to higher dimensions.The dimensionality-dependent scaling behavior is actually well known in the physics community;see,for example,the pioneering work in Kang and Redner(1984,1985),which showed that the reactant concentration decays as(where C0Ais the initial concentration of A,and d＜4 denotes the spatial dimension)for reaction A+B→inert,when the initial densities are equal and when the diffusion is isotropic.In the following sections,we will check the feasibility of R.Due to the lack of reliable observation of 3D reactive kinetics,we leave the validation of Eq.(13)for a future study.

2.2.Simplification for one-dimensional bimolecular reactions

For comparison purposes,the derivation described above is reduced to the well-studied,fundamental bimolecular reactions in 1D media,as the particle-based simulation schemes have been described in the literature(Benson and Meerschaert,2008;Edery et al.,2009,2010).The rate equation for a 1D irreversible reaction A+B→C can be obtained from Eq.(1)after neglecting advection and dispersion:

Following the procedure in section 2.1,the resultant 1D interaction radius for reactant A particles is

where N0A[dimensionless]is the initial number of A particles.The interaction radius in Eq.(16)is analogous to the empirical probability of forward reaction(denoted as P)proposed by Benson and Meerschaert(2008):

where Ω[L]is the domain size,N0is the initial number of particles(note that here C0A=C0B),and μ（s）[L-1]is associated with the density of two reactant particles separated by a distance s.

For irreversible reactions with unequal initial concentrations between A and B,the interaction radius R expressed by Eq.(15)is no longer a constant,but is self-adjusted according to the relative concentration between A and B.In other words,the ratio of initial concentrations between reactants can affect the interaction radius R and therefore the reaction kinetics.

It is also noteworthy that the interaction radius for B particles,denoted as RB,can be different from that for A particles for unequal initial reactant concentrations.Following the argument proposed above,the equation for RBis

The interaction radius RBin Eq.(18)for B particles(where the B particle is located in the center of the interaction volume)is equal to the interaction radius R in Eq.(15)for A particles(where the A particle is located in the center of the interaction volume)only for equal initial reactant(molar)concentrations.For unequal initial concentrations,the interaction radius between A and B changes,due to the different numbers of A and B particles.The resultant reaction rate and kinetics,however,remain the same,no matter whether we put the A or the B particle in the center of the interaction volume(shown in Fig.1).

3.Development of a fully Lagrangian simulator for reactions

The interaction radius in Eqs.(11)through(13)can be incorporated into fully particle tracking-based Lagrangian schemes,which are developed to quantify reactive kinetics for pollutants in aquifers involving both molecular-scale reaction and REV-scale transport.The operator-splitting technique is used to describe reactive transport,where the displacement of a contaminant due to advection and/or diffusion/dispersion is separated from the mixing-controlled reaction and can be tracked using the standard Markov process(LaBolle et al.,1996,2000).Here we assume that chemical reactions do not change the flow velocity or dispersion,and a coupled reaction and transport model can be built by refreshing the transport strength at the beginning of each time step.

A flowchart for the Lagrangian framework is shown in Fig.2.There are two options to handle the reverse component in reversible reactions.The first option is to split the reverse and forward reactions,where the reverse reaction must be accounted for before particle displacement and forward reaction(Fig.2).In this option,the interaction radii for the first forward reaction A+B→C and the second forward reaction C+A→D are the first term and the third term on the RHS of Eqs.(11)through(13),respectively.Meanwhile,the interaction radii for the first reverse reaction A+B←C and the second reverse reaction C+A←D are the second term and the fourth term on the RHS of Eqs.(11)through(13),respectively.Note that the sequence of reverse and forward reactions is critical.If the reverse reaction is simulated after the forward one,the reverse reaction will be overestimated because the newly generated reverse products do not have the right probability to collide in the same time step.In contrast,the reverse reaction tends to be underestimated if the reverse reaction occurs before the forward reaction.This is because the reverse products,such as the reactants A and B for the first reverse reaction,are forced to undergo forward reaction immediately,without experiencing any motion process.Numerical experiments(not shown here)also revealed that splitting the reverse reactions into different portions cannot obviate the problem.

While simulating reversible reactions,the probability of reactions of two adjacent reactants relates to the interaction radius R defined by the second and fourth terms on the RHS of Eqs.(11)through(13),where a random number will be generated to control randomly whether the reverse reaction will proceed.Note that the reaction does not necessarily occur even if two molecules are closer than R,since other properties,such as the orientation and bimolecular attraction,also affect the reaction.

Fig.2.Flowchart of Lagrangian method for simulating reversible bimolecular reactions(for irreversible reactions,options 1 and 2 will be deleted).

The second option is to apply the complete expression of Eqs.(11)through(13)to calculate the forward reaction(Fig.2).The reverse reaction is not needed in this method,since Eqs.(11)through(13)account for the overall, final probabilities of reactive particles.For example,when the reversible reaction reaches equilibrium,R in Eqs.(11)through(13)decreases to zero,implying that the concentration of reactants reaches a constant.Both methods will be checked in section 4.

4.Validation and application

The Lagrangian framework described above has been extensively tested,with some examples shown below.The modeled reaction rate was also compared to the classical thermodynamic reaction rate,to evaluate the impact of incomplete mixing of reactants on chemical reactions.

4.1.Reactive kinetics in 1D porous media:Lagrangian framework versus classical rate equations

In numerical experiments,the Lagrangian solver captures the particle dynamics,which can be checked against the known kinetics of diffusion-controlled,simple reactions such as the 1D bimolecular reaction.We first check the irreversible reaction A+B→C when A and B have the same initial density(Fig.3).Dimensionless parameters were used for numerical experiments in this study,for description simplicity.Over a relatively short time,the Lagrangian solution matches the analytical solution of the rate equation(which is CA（t）=C0A/（1+C0AKft）)(Fig.3),since there are enough A and B particles remaining in the system that can mix well.A relatively small Kfwas selected here,resulting in a relatively small interaction radius.A larger Kfwas tested with the same conclusion,except that the discrepancy between the Lagrangian solution(for the diffusion-controlled reaction)and the analytical solution(assuming perfect mixing)appeared at a later time.At later times(when the particles cannot maintain the same high degree of mixing),the Lagrangian solution deviates significantly from the analytical solution given above(see Fig.3(b)),confirming the well-known fact that a mixingcontrolled reaction(in real-world heterogeneous media)exhibits a reaction rate apparently lower than the reaction with perfect mixing as assumed by the classical thermodynamic rate law.An increase in the dispersion coefficient(leading to higher mobility)and the number of particles(leading to more mixing)will enhance the contact of reactants,and hence accelerate the reaction,as shown by Fig.4.The separate impacts of particle number and the dispersion coefficient on reaction kinetics are described in Zhang et al.(2013),and therefore the similar results are not shown here.

Fig.3.Lagrangian solutions(symbols)versus analytical solutions(solid lines)for irreversible reaction A+B→C with equal initial reactant concentrations(parameters are as follows:Kf=0.02,D=0.01,Δt=1,N0A=5000,and the initial reactant concentration is 2).

Fig.4.Lagrangian solutions(symbols)versus analytical solutions(solid lines)for irreversible reaction A+B→C with equal initial reactant concentrations and different numbers of particles(the dashed lines denote the asymptotic concentration,and the domain size L=64).

We then modeled the irreversible reaction A+B→C with unequal initial reactant concentrations(Fig.5).The modeled concentration for reactant A either decreases exponentially,when the initial concentration of A is smaller than the initial concentration of B:C0A＜C0B(Fig.5(a),where C0A=0.005 and C0B=0.01),or it reaches the limit C0A-C0Bif C0A＞C0B(Fig.5(b),where C0A=0.005 and C0B=0.0046).The analytical solution for the rate equation in both Fig.5(a)and(b)can be written as CA（t）= α/［（1+ α/C0A）exp（αKft）-1］,where the factor α=C0A-C0B(Zhang et al.,2013).

For reversible reactions,we compare the two simulation options discussed in section 3.When both the reverse and forward reactions are modeled(i.e.,option 1),the resultant concentration at later times contains more noise than that using the fully time-dependent R(t)(i.e.,option 2),as shown in Fig.6(a).This discrepancy is due to the additional step in option 1,which accounts for the random,reverse reaction.Apparently,option 2 is more computationally efficient,although this does not imply that the resulting smoother concentration is superior to that generated by option 1.Fig.6(b)shows the evolution of the interaction radius.The analytical solution for the rate equation is as follows(Kang and Redner,1985;Zhang et al.,2013):

4.2.Two-dimensional bimolecular reaction

We extended the 1D Lagrangian solver validated above to multi-dimensional reactions with orthogonal Fickian spreading rates.The corresponding interaction radius(Fig.1)can be transformed conveniently among dimensions using Eqs.(11)through(13).

We first checked the evolution of reactant concentration in a 2D domain.Some numerical examples are shown in Fig.7.Reactants A and B particles are released randomly in a square,and the motion of reactants is driven by isotropic diffusion.The Lagrangian solutions capture the full range of reaction kinetics for all degrees of mixing:fast reaction when reactant particles are well mixed,and the gradual transition to the t-1/2scaling law with incomplete mixing.

We then checked the pattern of reactant particles in forward reaction A+B→C in a 2D domain using the Lagrangian solver(Fig.8).Segregation is observed at later times for reactant particles that were randomly distributed at the beginning(e.g.,Fig.8(c)),recovering the well-known segregation effect for 2D reaction first identified by the chemical physics community(Kopelman,1988).While the classical papers for particle segregation were limited to simple annihilation reactions(Toussaint and Wilczek,1983;Kopelman,1988),the fully Lagrangian solver developed in this study can be extended to multi-step,complex reactions.

Fig.5.Lagrangian solutions(symbols)versus analytical solutions(solid lines)for irreversible reaction A+B→C with unequal initial reactant concentrations.

Fig.6.In fluence of interaction radius R on chemical concentration for reversible reaction A+B⇌C for constant R and time-dependent R(t)(CA（t）and CC（t）were solved by the analytical approach(lines)and the Lagrangian methods(symbols),respectively,with parameters as follows:Kf=0.004,Kr=0.015,D=0.1,Δt=1,N0A=2000,and C0A=C0B=5).

Finally,we applied the Lagrangian solver to simulation of precipitation of calcium carbonate minerals observed in a 2D flow cell by Scheibe et al.(2007).In their laboratory experiment,two solutes,Na2CO3and CaCl2,were injected into a flow cell filled with 0.5 mm quartz sand.The bimolecular irreversible reaction caused precipitation of CaCO3within a narrow zone at the center of the experimental cell.The simulated evolutions of reactants and products are shown in Fig.9.The apparent segregation of reactants and products is consistent with that observed in the laboratory(see Fig.1 in Scheibe et al.(2007)).The laboratory experiments conducted by Scheibe et al.(2007)showed that the interface between Na2CO3and CaCl2is not the regular rectangle(with straight/smooth sides)predicted by the standard continuum model(i.e.,the 2D extension of Eq.(1))or a standard Eulerian approach,but a curved interface filled with random noise,just like our numerical simulation shown in Fig.9(a)and(b).We emphasize here that,to the best of our knowledge,the Lagrangian approach is the only viable tool so far for capturing the irregular interface observed in the experiments conducted by Scheibe et al.(2007).Fig.10(a)and(b)show the uneven distribution of reactants and product particles in the horizontal and vertical directions,respectively,due to the mixing-controlled reaction.Note that no actual data were provided by Scheibe et al.(2007);here,the numerical solutions can only be compared visually with their results.

4.3.Complex chemical reactions

The Lagrangian solver described in section 3 was used to quantify the multi-step reactions(Eq.(2))described by Oates and Harvey(2006),where the colorimetric reaction of Tiron and molybdate is suitable for optical quantification of chemical reaction during fluid- fluid mixing(Fig.11).The transport model used in the Lagrangian method is the classical 2ndorder advection-dispersion equation.A 0.05 mol/L Tiron solution was injected into a 1D,saturated,glass-bead column filled with a 0.025 mol/L molybdate solution at pH 6.1.Snapshots of chemicals at four times(8,18,32,and 46 min),including both the two colorless reactants(i.e.,Ti and Mo)and the two colored,soluble products(MoTi and MoTi2),were digitally imaged and quantified using light absorbance(see symbols in Fig.12).

Fig.7.Lagrangian solutions(symbols)for irreversible reaction A+B→C with equal initial reactant concentrations in a 2D domain(model parameters are as follows:Kf=0.1,C0A=C0B=1,and the domain size L=1×1;and the dashed lines denote the asymptotic concentration).

Fig.8.Simulated particle position for irreversible reaction A+B→C in a 2D domain at different times(model parameters are as follows:Kf=50,D=0.01,Δt=0.1,N0A=N0B=4000,and C0A=C0B=1;and reactant particles are located randomly in the 2D domain at time t=0).

Our numerical approach is described below.Due to the limited information about chemical transport,we assumed that all chemicals moved with water at the same velocity v,and dispersed in space with the same dispersion coefficient D.The velocity v=0.67 cm/min and the dispersion coefficient D=0.50 cm2/min estimated by Oates and Harvey(2006)were adopted here.The best- fitforward and reverse reaction rates are as follows:=9.8 L/(mmol·min)and=0.02 min-1for the firstcomplexation reaction(Eq.(2a)),and=1.9 L/(mmol·min)and=0.01 min-1for the subsequent reaction(Eq.(2b)).These reaction rates were fitted using the first snapshot for product CC+CDat time t=8 min,and then used to predict the snapshots at later times(t=18,32,and 46 min)(Fig.12).The first forward rate constant is about 4 times larger than the second forward rate constant.The two best- fit forward rate constants are on the same order as those found by Oates and Harvey(2006),and the two best- fit reverse kinetic coefficients are relatively small since the products are relatively stable.It is also noteworthy that the classical continuum model,which assumes perfect mixing for chemical particles,significantly overestimates the reaction rate and the product concentration(Oates and Harvey,2006)(see the dotted lines in Fig.12).Our numerical results are acceptable since(1)they significantly improve the standard model in capturing the observed trend of product snapshots(affected by incomplete mixing),and(2)the Lagrangian parameters fitted by the first snapshot predict well all of the following data.

Fig.9.Snapshots of bimolecular reaction in a 2D flow cell(with a domain size of 1×2)at t=10 and 100(model parameters are as follows:Dx=Dy=0.01(with Dxand Dyrepresenting the dispersion coefficient along the x and y directions,respectively),Kf=50,and the drift velocity is zero,i.e.,Fickian diffusion with re flective boundaries).

Fig.10.Latitudinal and longitudinal distributions of particles in a 2D flow cell with a domain size of 1×2(model parameters are as follows:Dx=Dy=0.01,Kf=50,and the drift velocity is zero).

5.Discussion

5.1.Dimension-dependent reaction kinetics

Our numerical experiments show that the diffusioncontrolled chemical reaction kinetics are sensitive to the medium dimensions.For bimolecular reactions in onedimensional media,the Lagrangian solution of reactant concentration at later times declines much slower than at earlier times,following a power-law function CA（t）～ t-1/4shown in Fig.4.This trend is consistent with the kinetics of mixingcontrolled reactions(Sung and Yethiraj,2005,among many others).This well-known 1/4 scaling law of the reactant concentration decay rate is likely due to the speed of the interface of reactant particle clouds in 1D media(where reactant particles are distributed as islands),which expands in space at the rate of approximately t-1/4.

Fig.11.Lagrangian simulation of snapshots for all four chemicals at t=32 min(the dashed line shows the numerical solution of MoTi2).

For bimolecular reactions in 2D media,reactant concentration at later times declines as slow as a power-law function CA（t）～ t-1/2.At later times,reactions can only occur around the interface of different reactant segregations(e.g.,Fig.8(c)),implying a pattern of self-organization and self-protection of species.The expansion of the 2D plume interface may explain the 1/2 scaling law of the reactant concentration decay rate(Kang and Redner,1984,1985),as demonstrated in Fig.7.

5.2.Model extension:Higher-order reaction and non-Fickian transport

Fig.12.Best- fit product concentration using Lagrangian method(solid lines)versus laboratory measurements(symbols,Oates and Harvey,2006)at time t=8,18,32,and 46 min(the dashed lines represent the solutions of the classical continuum model assuming perfect mixing in Oates and Harvey(2006)).

The Lagrangian scheme developed in this study can be extended to other complex reactions.For example,a forward third-order reaction A+B+D→E might be handled as a sequence of A+B→C and C+D→E,as suggested by Benson and Meerschaert(2008).The multi-step Lagrangian approach developed in this study can be conveniently extended to capture this higher-order reaction.We will check the feasibility of this extension in a future study.

The flexibility of the Lagrangian simulator allows further re finement to capture complex transport affecting chemical mixing and hence reactions within and across scales.For example,the transport of tracers in natural porous and fractured media is typically non-Fickian(see the extensive review byBerkowitzetal.(2006)andNeumanand Tartakovsky(2009)),which can be captured by space and/or time nonlocal transport models.Non-Fickian transport can be added to the calculation of displacement for particles,using the renewal-reward approach proposed by Zhang et al.(2015).In addition,the local-scale simulation of chemical reactions can be combined with the regional scale,nonlocal particle tracking scheme in natural aquifers with any size and boundary conditions(Zhang et al.,2016),which may provide a multi-scale simulation of subsurface biogeochemical processes.These extensions will be pursued in a future study.

5.3.Physical meaning of interaction radius R

The interaction radius R affects the reaction rate.We explain the impact of R on the reaction rate using the irreversible reaction A+B→C in a 1D domain with size L,which is modeled with two steps in the Lagrangian simulation.We start from one A particle(A is put at the center of the reaction volume,as shown by Fig.1)and look for the surrounding B particle.The distance between each pair of reactants is evaluated by Zhang et al.(2013):

where RAis the interaction radius for A,and ds denotes the distance between A and B particles.Second,a uniform[0,1]random number ω is generated and compared with P*:if ω≤P*,A and B particles combine to produce a C particle;otherwise,no chemical reaction can occur.As demonstrated by Zhang et al.(2013),this two-step scheme(with independent probabilities)leads to the average reaction probabilityfor all reactants:

where RAcan also be regarded as the maximum range to search for a B particle,and L/NA（t）denotes the average space occupied by each A particle at time t in the domain.

In the Lagrangian simulation,one A particle interacts with one B particle to form one C particle,leading to the following mass balance law:

where ΔCA（t）and ΔCB（t）are,respectively,the magnitudes of change in concentrations of A and B in time step Δt.The relative change of reactant concentrations for A and B in Δt can be defined as follows(Zhang et al.,2013):

Inserting Eq.(24)into Eq.(23),one obtains

which can be re-arranged as follows:

The derivation above(from Eqs.(23)through(26))shows that,due to mass balance,different reactants should have different interaction radii,if they have unequal initial concentrations.The higher the initial concentration of one reactant is,the smaller its interaction radius.If there are more A particles than B,then,for an A particle to collide with an adjacent B particle,the search radius needs to be relatively large(because the B particles are sparse).However,if we put B at the center of the reaction volume and search the surrounding A particles for reaction,the search radius is shorter due to the denser distribution of A particles.Eq.(25)shows that it is not the equivalence between the interaction radius(RAand RB),but the equivalence between the product RACA（t）NA（t）and RBCB（t）NB（t）that is required to maintain mass balance.It is also noteworthy that Eq.(26)can be directly derived by dividing Eq.(15)by Eq.(18),which also validates the interaction radius derived by this study.

In addition,we reiterate that the interaction radius can change with time.For example,for the 1D reversible reaction A+B⇌C,the corresponding R is defined by Eq.(11):

At early times,the first term on the RHS of Eq.(27)denoting the forward reaction is positive and generates the largest interaction radius R(t).With the evolution of time,the second term on the RHS of Eq.(27)(which is negative and represents the reverse reaction)strengthens(since there are more C particles),decreasing R(t).When the forward and reverse reactions reach equilibrium,R(t)approaches its minimum(zero)(Fig.6(b)),so that the concentrations of all reactants and products remain stable after this time.Numerical tests show that the decreasing R(t)generates Lagrangian solutions matching the analytical solutions,as shown,for example,by Fig.6(a).

6.Conclusions

This study developed a fully Lagrangian framework to simulate multi-step,multi-scaling,mixing-limited chemical reactions in multi-dimensional saturated porous media,where the random motion of chemical particles is assumed to follow Brownian motion.The particle tracking results were then checked against known reaction kinetics,and used for highdimensional or multi-step chemical reactions documented in the literature.Three major conclusions are obtained from this study.

(1)The core of the Lagrangian approach is to define the interaction radius R,which can be estimated using the agentbased approach for media with any dimension and reactions with multiple steps.The derived R is a function of time for most types of reactions.For example,for the case of irreversible bimolecular reactions,R is a constant related to the initial reactant density,while R decreases nonlinearly with time for all the other reactions considered in this study,including(1)irreversible bimolecular reactions with unequal initial reactant concentrations,(2)reversible bimolecular reactions,and(3)multi-step bimolecular reactions.The interaction radius R is also dimension-dependent,scaling as a power-law function with an increasing dimension.

(2)Numerical experiments show that the novel Lagrangian solver can capture major reaction kinetics.In particular,tests show that the diffusion strength and particle number density significantly affect the reaction rate,and the power-law scaling rate for reactant mass changes from 1/4 in 1D media to 1/2 in 2D media,due likely to the expansion rate of plumes varying depending on the number of dimensions.Particle dynamics simulated by the Lagrangian solver also reveal the discrepancy between diffusion-limited reactions and thermodynamic ratelimited reactions.Such a discrepancy is consistent with that shown in the literature,validating the numerical solver.

(3)Applications show that the Lagrangian solver can ef ficiently simulate precipitation of calcium carbonate minerals in a 2D medium,and quantify multi-step chemical reactions observed in the laboratory.The flexibility of the Lagrangian simulator allows further re finement to capture complex transport affecting chemical reactions,such as higher-order chemical reactions or non-Fickian transport observed frequently in groundwater.