On the localization of chemical reactions in multicontinuum media

Advances in Water Resources 167 (2022) 104286 Available online 28 July 2022 0309-1708/© 2022 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). On the localization of chemical reactions in multicontinuum media Jingjing Wang a , b , c , d , * , Jesus Carrera c , d , Maarten W. Saaltink b , d , Cristina Valhondo c , d , e a Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai 200092, China b Department of Civil and Environmental Engineering, Universitat Polit` ecnica de Catalunya (UPC), Jordi Girona 1-3, 08034 Barcelona, Spain c Geosciences Department, Institute of Environmental Assessment and Water Research (IDAEA), Severo Ochoa Excellence Center, Spanish Council for Scientific Research (CSIC), c/ Jordi Girona, 20, 08034 Barcelona, Spain d Associated Unit—Hydrogeology Group (UPC-CSIC), Jordi Girona, 08034 Barcelona, Spain e G´ eosciences Montpellier, Universit´ e de Montpellier, CNRS, 300 Avenue Emile Jeanbrau, CC MSE, 34095 Montpellier, France ARTICLE INFO Keywords: Chemical localization Multirate mass transfer Reactive transport Multicontinuum media Numerical modeling ABSTRACT Reactive transport (RT) through heterogeneous media, may cause chemical heterogeneity if water flux is slow through portions of the medium. In such cases, chemical localization (i.e., the occurrence of reactions that would not occur in well mixed media) may develop, which is especially relevant for biochemical reactions occurring in biofilms. The objective of this work is to study the conditions for chemical localization. We represent the impact of heterogeneity by means of the non-local multirate mass transfer (MRMT) model, which views the porous media as consisting of one mobile and many immobile zones. A dimensional analysis of the governing equations shows that the problem is characterized by reaction times and the distribution of residence times in immobile zones, relative to transport time. To analyze the interplay between them, we simulated simple RT problems in multicontinuum media. Results indicate that immobile zones with residence times much smaller than transport can be lumped together with the mobile zone by modifying the reaction rates, which reduces computations. More importantly, reactions driven by species that are not present in the inflowing water but are the result of previous reactions will take place preferentially in immobile zones, whose residence time is comparable to or larger than reaction times. In fact, daughter species may take a long time and distance to build up. That is, daughter species will not be largest near the inflow, where parent species display largest concentrations, but further downstream at isolated (long residence times) immobile zones. 1. Introduction Reactive transport (RT) in porous media is essential for understanding environmental processes, such as the fate of contaminants in aquifers, and for a proper design and management of remediation strategies, such as managed aquifer recharge (e.g., Valhondo et al., 2020) or CO 2 geological storage (e.g., Agrawal et al., 2021). Carrera et al. (2022) argue that the basic problem of RT is how to address the chemical and physical heterogeneity of natural media. Physical heterogeneity has received broad attention in the transport literature. It is well known that transport in heterogeneous porous media displays anomalous (non-Fickian) behavior (Kitanidis, 1988) both at field (Adams and Gelhar, 1992) and laboratory scales (Berkowitz and Scher, 2009; Levy and Berkowitz, 2003; Valocchi, 1985), displaying early arrivals, long tails and scale-dependent dispersion coefficients (Carrera, 1993). The problem is especially relevant for reactive transport because numerous reactions depend on local scale mixing. To overcome this problem, researchers have focused on small scale simulations in an attempt to simulate heterogeneity explicitly (e.g., Cirpka et al., 2015) or, even by pore scale simulations (e.g., Agrawal et al., 2021). Beyond the cost of such simulations, what emerges is that upscaling parameters does not suffice, but the transport equation may have to change and the advection-dispersion-reaction equation (ADRE) cannot be considered suitable for reactive transport (Carrera et al., 2022; Molins and Knabner, 2019; Soler-Sagarra et al., 2022). Many non-local methods have been developed to address the need for an alternative transport equation. These include continuous time random walks (CTRW) (Berkowitz et al., 2006; Berkowitz and Scher, 1998; Dentz et al., 2015, 2004), fractional advection-dispersion equations (fADE) (Benson et al., 2000; Benson and Meerschaert, 2009; Marseguerra and Zoia, 2008; Schumer et al., 2003a, 2003b), multirate mass transfer (MRMT) (Babey et al., 2017; Benson and Meerschaert, 2009; Dentz et al., 2011; Fernandez-Garcia and Sanchez-Vila, 2015; * Corresponding author. E-mail addresses: [email protected] (J. Wang), [email protected] (J. Carrera), [email protected] (M.W. Saaltink). Contents lists available at ScienceDirect Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres https://doi.org/10.1016/j.advwatres.2022.104286 Received 30 March 2022; Received in revised form 10 July 2022; Accepted 26 July 2022 Advances in Water Resources 167 (2022) 104286 2 Haggerty and Gorelick, 1995; Salamon et al., 2006; Wang et al., 2005), multiple interacting continua (MINC) (De Dreuzy et al., 2013; Pruess and Narasimhan, 1985), structured interacting continua (SINC) (Babey et al., 2015; Rapaport et al., 2017), memory functions (Carrera et al., 1998; Gouze et al., 2008; Haggerty et al., 2000; Willmann et al., 2008) and so forth. Although these methods use different approaches, they are essentially equivalent (Dentz and Berkowitz, 2003; Neuman and Tartakovsky, 2009; Silva et al., 2009). Unfortunately, their solution is often based on Laplace transform methods, which makes them poorly suited for multicomponent and/or non-linear RT. Therefore, most work on non-local methods deals with the behavior of breakthrough curves (BTCs) and conservative transport in heterogeneous porous media (Berkowitz and Scher, 2009; Dentz et al., 2015; Haggerty et al., 2000; Schumer et al., 2003b; Willmann et al., 2008). Still, MRMT, which can be considered the simplest paradigm of non-local transport methods, can be used for multicomponent RT. The medium is viewed as consisting of a mobile and several overlapping immobile zones. Each of these zones can be modeled numerically, which facilitates general RT simulation (Babaei and Islam, 2018; Massoudieh and Dentz, 2020; Soler-Sagarra et al., 2016; Willmann et al., 2010). Chemical heterogeneity has received a lot less attention than physical heterogeneity. Chemical heterogeneity may refer to large scale variability, usually known from geological understanding, or to pore scale variability. Portions of the porous media, grain aggregates or individual grains or portion of grains can have different mineral compositions, leading to different reaction rates and different concentrations at different parts of the pore space. Furthermore, pore scale variability, causes water fluxes to be spatially variable, which causes dispersion and affects solute transport. Areas with low water flux may develop chemical conditions different from those of high flux, which carry most of the water (Fig. 1). In fact, reacting surface tend to be larger in low flux than in high flux areas. That is, physical heterogeneity may induce chemical heterogeneity. Steefel et al. (2005) call microenvironments these parts of the domain. We term the phenomenon chemical localization to emphasize the local occurrence of reactions. Furthermore, heterogeneity can cause chemical localization, i.e., the occurrence of reactions that would not occur in well mixed media. The importance of chemical localization has been widely recognized in the literature. The discrepancy between measured mineral dissolution rates in laboratory and field is attributed to heterogeneity (Li et al., 2008). Modeling studies, sometimes combined with chemical and geophysical datasets, have shown the effect of heterogeneity on many particular systems, such as bioremediation of contaminated aquifers (Scheibe et al., 2006), reactive transport processes in fractures (Deng and Spycher, 2019) and mineral precipitation and dissolution in Nomenclature α Mass transfer rate, [T −1 ] f( α ) Probability density function of mass transfer rates τ α – – – 1/ α Residence time, [T] P( τ α ) Probability density function of residence times α 0 Scale parameter of f( α ) τα 0 Scale parameter of P( τ α ) β Shape parameter of f( α ) and P( τ α ) Γ(β) Gamma function c m Concentrations in the mobile zone, [ML −3 ] c im Concentrations in the immobile zones, [ML −3 ] ϕ m Porosity of mobile zone, [-] ϕ im Porosity of immobile zones, [-] q Darcy flux, [LT −1 ] D Hydrodynamic dispersion, [L 2 T −1 ] r m Sink-source term due to chemical reactions in the mobile zone, [ML −3 T −1 ] r im Sink-source term in the immobile zones, [ML −3 T −1 ] k m Reaction rate constant in the mobile zone, [T −1 ] k im ( τ α ) Local reaction rate constant in the immobile zone with τ α , [T −1 ] τ r – – – 1/k im ( τ α ) Reaction time in the immobile zone with τ α , [T] g(t) Memory function κ(t) Reaction rate kernel δ(t) Dirac delta function L c Characteristic length, [L] t c Characteristic time, [T] c c Characteristic concentration, [ML −3 ] q Modulus of q, [LT -1 ] D Longitudinal component of D, [L 2 T −1 ] v =q/ϕ m Mean fluid velocity in mobile zone, [LT −1 ] c mD Dimensionless concentrations in mobile zone (c mD =c m / c c ), [-] c imD Dimensionless concentrations in immobile zones (c imD =c im /c c ), [-] x D Dimensionless distance (x D =x/L c ), [-] t D Dimensionless time (t D =t/t c ), [-] η Ratio of total immobile to mobile porosity ( η =ϕ im /ϕ m ), [-] τ α D Dimensionless residence time ( τ α D = τ α /t c ), [-] τ rD Dimensionless reaction time ( τ rD = τ r /t c ), [-] c e Concentration of the inflowing water, [ML −3 ] cP m,1 Mobile concentration in the first node for the pulse injection, [ML −3 ] L Travel distance of solute, [L] t A Advective transport time, [T] ϕ =ϕ m +ϕ im Total porosity, [-] L{}(s)Laplace transform s Complex Laplace variable Fig. 1. Illustration of chemical localization: water flux tends to concentrate on a small fraction of the medium (blue arrows), often representing a small fraction of the reacting surfaces within a possible REV (Representative Elementary Volume). Proper mass balance requires model concentrations to be those of the flowing water, while some reactions may be localized in non-flowing portions of the pore space. The problem is especially relevant for biochemical reactions that occur in biofilms (green areas), possibly under redox conditions different from those of the mobile water. J. Wang et al. Advances in Water Resources 167 (2022) 104286 3 caprocks for CO2 storage (Tian et al., 2019). A clear example of chemical localization is presented by Luquot et al. (2016), who observed carbonate precipitation during injection of acidic CO2 brine into a shale core, which would not be possible if the acidic conditions associated to CO2 injection prevailed throughout the domain. Soler-Sagarra et al. (2016) showed that silicates weathering may cause sufficient alkalinity to facilitate carbonate precipitation in areas of low water flux. A similar situation may occur for biochemical reactions, which take place in biofilms. Biofilm growth generates immobile zones with limited access to electron acceptors, which facilitates the occurrence of reducing conditions (Bellin et al., 2014). The net result is that the microbial community is heterogeneous and reactions that would only occur under anaerobic conditions take place despite the overall aerobic water flux. In analyzing this problem, Dentz et al. (2011) concluded that transport under physical and chemical heterogeneity cannot be upscaled separately: upscaling reactions depends on physical heterogeneity and upscaling transport is affected by chemical heterogeneity. A similar conclusion was reached by Battiato et al. (2009). Furthermore, macroscale transport may be determined by microscale transport process in immobile zones (Gouze et al., 2008). The situation cannot be considered satisfactory. On the one hand, chemical localization is frequently observed and we conjecture that its importance will increase with the need to simulate biochemical RT (Carrera et al., 2022). Yet, it is not clear how to upscale the processes. What is clear is that a multicontinuum representation is needed to be able to host and represent different chemical environments. In fact, Lichtner and Kang (2007) adopted a lattice Boltzmann solution method with this goal. We contend that a MRMT representation can reach the same goal at a much lower cost. Still, the problem has never been analyzed formally, so that numerous questions need to be answered regarding the conditions for localization and when does it need to be modeled. The objective of this work is to study the conditions for localization of reactions in MRMT RT formulations. To do so, first, we establish the governing equations of reactive transport in multicontinuum media and deduce a dimensionless form of these equations. We then use the resulting model in two simple cases to analyze the conditions for localization. 2. Methodology We model reactive transport in heterogeneous porous medium using the MRMT approach (Haggerty and Gorelick, 1995). Every point (representative elementary volume) is viewed as consisting of a mobile zone and a distribution of immobile zones characterized by their residence times. Each immobile zone exchanges solute mass with the mobile zone proportionally to the mass transfer rate, α [T −1 ] (a complete description of variables employed in this paper is given in the Nomenclature). Therefore, it is natural to characterize this distribution by a probability density function (pdf) of mass transfer rates, f( α ) (see Fig. 2 left). To facilitate comparison to other non-local methods, we also characterize immobile zones by the distribution, P( τ α ), of residence times, τ α – – – 1/ α [T] (Haggerty et al., 2004) (see Fig. 2 right). P( τ α ) is given by (Dentz et al., 2015, 2011; Dentz and Berkowitz, 2003) P( τ α ) = τ −2 α f( α ) = τ β α 0 Γ(β) 1 τ β+1 α e− τα 0 τα (1) where β is the shape parameter, τα 0 is the scale parameter, and Γ(β) is the Gamma function. The residence time probability P( τ α )d τ α is the immobile zones frequency with residence times in the interval [ τ α , τ α +d τ α ]. Obviously, it satisfies the condition ∫∞ 0P( τα )d τα =1. 2.1. Governing equations The transport of any reactive species in a medium that consists of a mobile continuum and multiple immobile continua can be expressed as (Dentz et al., 2011; Donado et al., 2009; Willmann et al., 2010) ϕm ∂ cm ∂ t=Lt[cm] − ϕim ∫ ∞ 0 1 τ α P( τ α )[cm(x,t) − cim(x, τ α ,t)]d τ α +ϕmrm(2) where L t [c m ] = − [q ⋅ ∇c m ] + ∇ ⋅ [ϕ m D ⋅ ∇c m ] is the transport operator that describes advection and dispersion, q [LT −1 ] is Darcy flux, D [L 2 T −1 ] is the hydrodynamic dispersion tensor; ϕ m [-] and ϕ im [-] are the porosities of mobile and immobile zones, which denote the volume fraction of mobile pore and immobile pore water over the bulk volume, respectively; c m [ML −3 ] and c im [ML −3 ] are concentrations in the mobile and immobile zone, which are expressed as mass per unit volume of mobile and immobile water, respectively; r m [ML −3 T −1 ] is the sinksource term due to chemical reactions in the mobile zone and corresponds to the mass removed by reactions in the mobile zone per unit volume of mobile water per unit time. It should be noted that Eq. (2) remains valid for sorbing solutes, in which case ϕ m and ϕ im are the storage capacity (porosity times local retardation coefficient) of mobile and immobile zones, respectively. Mass balance in each immobile zone is given by ∂ cim(x, τ α ,t) ∂ t=1 τ α [cm(x,t) − cim(x, τ α ,t)] + rim(x, τ α ,t)(3) where r im [ML −3 T −1 ] is the reactions sink-source term in the immobile zone with residence time τ α , which corresponds to the mass removed by reactions in the immobile zone per unit volume of immobile water per unit time. Mass exchange between mobile and immobile zones is modeled by a first-order mass transfer term, represented by the continuous variable τ α Fig. 2. Pdf of mass exchange rates, α (Gamma distibution), and residence times, τ α , plotted in log-log scale for several values of τ α 0 and β =1/2. Note that the maximum of the latter is around τ α 0 (the mode is τ α 0 /(β +1) and the expected value τ α 0 /(1 −β)). Its slope at large τ α is β +1. J. Wang et al. Advances in Water Resources 167 (2022) 104286 4 and characterized by a distribution density function P( τ α ) for the immobile zones. Thus, the total mass exchange is the weighted sum over all immobile zones as expressed in the integral term of Eq. (2). Integrating Eq. (3) with weight P( τ α ) over τ α , multiplying by ϕ im , and substituting it into Eq. (2), yields the total solute mass balance, ϕm ∂ cm(x,t) ∂ t+ϕim ∫ ∞ 0 P( τ α ) ∂ cim(x, τ α ,t) ∂ td τ α =Lt[cm(x,t)] + r(x,t)(4) in which, r is the total reaction rate (now per unit volume of bulk porous medium) that integrates reactions in both mobile and immobile zones, i. e., r(x,t) = ϕmrm(x,t) + ϕim ∫ ∞ 0 P( τ α )rim(x, τ α ,t)d τ α (5) If the reaction follows first-order kinetics, the reaction rate in the mobile zone will be rm(x,t) = − kmcm(x,t)(6) where k m [T −1 ] is the reaction rate constant. Similarly, the local reaction rate in the τ α immobile zone becomes rim(x, τ α ,t) = − kim( τ α )cim(x, τ α ,t)(7) where k im ( τ α ) [T −1 ] is the local reaction rate constant. The reaction time, τ r – – – 1/k im ( τ α ) [T], is assumed to depend on the immobile zone. This dependence should be understood as indicating that, as suggested by Fig. 1, different reactions with different reaction rates may occur in different immobile zones. By inserting Eq. (7) into Eq. (3), then solving Eq. (3), we obtain the concentration in the immobile zone of τ α as a function of mobile concentration history cim(x, τ α ,t) = cim(x, τ α ,t=0)e −[1 τα +1 τ r]t +∫ t 0 φ(t−t′)cm(x,t′)dt′(8) where φ(t) is given by φ(t) = 1 τ α e −[1 τα +1 τ r]t (9) Assuming that the initial concentration in the immobile zones is zero (i. e., c im (x, τ α ,t =0) =0), then by substituting Eq. (8) into Eq. (2) and rearranging terms, we obtain the total solute mass balance, as a sole function of c m ϕm ∂ cm(x,t) ∂ t+ϕim ∂ ∂ t∫ t 0 g(t−t′)cm(x,t′)dt′ =Lt[cm(x,t)] − ∫ t 0 κ(t−t′)cm(x,t′)dt′(10) where the memory function (Carrera et al., 1998; Haggerty et al., 2000) and reaction rate kernel (Dentz et al., 2011) are given by g(t) = ∫ ∞ 0 P( τ α )φ(t)d τ α (11) and κ(t) = ϕm(x)1 τ r,m δ(t) + ϕim(x)∫ ∞ 0 1 τ r P( τ α )φ(t)d τ α (12) respectively, in which δ(t) is the Dirac delta function. The memory function can be viewed as the rate of change of concentration in the immobile zone that is caused by a unit change of concentration in the mobile zone at initial time t =0. In the presence of reactions, the memory function not only incorporates the distribution of local mass transfer rates but also the distribution of local reaction rates. Similarly, the reaction rate kernel represents the local reactions in immobile zones: it integrates the reactions in both mobile and immobile zones, and accounts for the mass transfer between mobile and immobile zones. Although the memory function is mainly controlled by the distribution of residence times in immobile zones, it is also affected by the local reaction times in immobile zones through φ(t) (Eq. (9)). Likewise, the reaction rate kernel is mainly controlled by the chemical reactions, but it is also affected by the distribution of residence times in immobile zones. Clearly, the governing equations of reactive transport under physical and chemical heterogeneity (10), (11) and (12) are controlled by the distribution of residence times and reaction times in immobile zones simultaneously. The physical and chemical heterogeneities interact with each other, and together they govern reactive transport (Dentz et al., 2011). 2.2. Dimensionless formulations We write the dimensionless form of Eqs. (10), (11)) and (12) using the definitions of Carrera et al. (2022), which allows us to avoid the need for a Peclet Number, so that all characteristic variables will be those describing the immobile zones distribution and reaction rates. Thus, we define the characteristic length Lc=ϕmD q(13) where q is the modulus of q and D is any norm of D (here its longitudinal component). And, the characteristic transport time tc=Lc v=ϕ2 mD q2(14) Note that in 1D, the characteristic length equals the longitudinal dispersivity α L if molecular diffusion is neglected in the mobile zone, because dispersion becomes D = α L v, where v =q/ϕ m is the mean fluid velocity in the mobile zone. Using these definitions to write all variables in dimensionless form (i.e., t D =t/t c , τ α D = τ α /t c , x D =x/L c , c mD =c m /c c , etc.) and defining η =ϕ im /ϕ m (see the nomenclature table at the end for the definition of all variables) and substituting them into Eq. (10), leads after some tedious but conventional algebra to the dimensionless form of the governing equations ∂ cmD ∂ tD + η∂ ∂ tD∫ tD 0 g(tD−t′ D)cmD(t′ D)dt′ D = − ∇DcmD + ∇2 DcmD −∫ tD 0 κ(tD−t′ D)cmD(t′ D)dt′ D (15) Note that in this dimensionless formulations, the physical and chemical heterogeneity are represented by the distribution of residence and reaction times (relative to the characteristic transport time t c ), and they are reflected simultaneously in the memory function and reaction rate kernel. Obviously, in the situation where the characteristic length L c and the dimensionless ratio of porosity of immobile zone to porosity of mobile zone η are invariable, the solution of Eq. (15) is governed by two characteristic dimensionless times, the dimensionless residence time τ α D (and its distribution) and the dimensionless reaction time τ rD in immobile zones. The solution of governing Eq. (15) in Laplace domain is expressed as J. Wang et al. Advances in Water Resources 167 (2022) 104286 5 L{cmD}(s)=exp[xD 2(1− 1+4[s[1+ η L{g}(s)] + L{κ}(s)] √)] (16) where L{⋅}(s)represents the Laplace transform of a function. The Laplace transform of the memory function and reaction rate kernel are defined as L{g}(s) = ∫ ∞ 0 P( τ α D)L{φ}(s)d τ α D(17) and L{κ}(s) = 1 τ r,mD + η ∫ ∞ 0 1 τ rD P( τ α D)L{φ}(s)d τ α D(18) in which, the Laplace transform of φ(t D ) is given by L{φ}(s) = 1 τ α D τ α D+ τ rD τ α D τ rD +s(19) If we substitute s for s[1+ η L{g}(s)] + L{κ}(s)in Eq. (16), then Eq. (16) becomes the solution of the transport equation in homogeneous media. Obviously, the retardation of localized physical heterogeneity on transport is reflected in the term s η L{g}(s), and the decay of localized chemical heterogeneity on reactive species is shown in the term L{κ}(s). 2.3. Steady state analytical solution A steady state solution is possible for the continuous injection (or cumulative concentration of the pulse injection). The relevant output for the chemical localization problem is the ratio between the species concentrations in immobile and mobile zone, given by ∫ ∞ 0 cimD(xD, τ α D,tD)dtD/∫ ∞ 0 cmD(xD,tD)dtD=L{φ}(s)|s=0(20) as can be obtained from the Laplace transform of Eqs. (8) and Eq. (9). The concentration in the immobile zone is the weighted integral of concentrations in all immobile zones, which is defined as cimD(xD,tD) = ∫ ∞ 0 P( τ α D)cim(xD, τ α D,tD)d τ α D(21) From Eq. (8), we can obtain that cimD(xD,tD) = ∫ tD 0 g(tD−t′ D)cm(xD,t′ D)dt′ D(22) Then, the ratio between the cumulative concentrations in the immobile and mobile zone becomes ∫ ∞ 0 cimD(xD,tD)dtD / ∫ ∞ 0 cmD(xD,tD)dtD=L{g}(s)|s=0(23) 3. Numerical models We have built several numerical models to study chemical localization in physical and chemical heterogeneous media. The models have been calculated by the code RT_MRMT_DSA (Wang et al., 2022), which is an object-oriented FORTRAN2003 code based on finite element method (FEM) and capable of simulating reactive transport in multicontinuum media with MRMT. Species transport and reactions are solved simultaneously by using the Direct Substitution Approach (DSA), based on the Newton-Raphson method (Saaltink et al., 2001; Yeh and Tripathi, 1989). Fig. 3. Mobile (solid lines) and average immobile concentrations (dash lines) versus time (left column) at x D =100L c in response to a pulse injection (left top row) and to a continuous injection (actually, cumulative concentration) (left bottom row) for transport of a conservative species in a medium with a gamma distribution of residence times in immobile zones with β =1/2, τ α 0D=102, 10 1 , 10 0 , and 10 −1 . The BTC peaks at the mobile advective time (t AD =100, red arrow) for slow exchange (i.e., large τ α 0D), but at the full porosity advective time (t AD =400, orange arrow) when exchange is fast, which also reduces the time it takes for immobile zones to equilibrate with inflow water. Cumulative concentration profile (colored, right column) versus dimensionless distance and dimensionless residence times (pdf with β =1/2 and τ α 0D=10−1) at t D =300. The white dash line indicates the immobile zone whose residence time equals the travel time t D . Note that immobile concentrations tend to mobile concentrations for τ α D <t D . J. Wang et al. Advances in Water Resources 167 (2022) 104286 6 All the models are one-dimensional with ϕ m =0.1 [-], ϕ im =0.3 [-], α L =1 [L], q =1 [LT −1 ]. Thus, we have L c =1 [L], t c =0.1 [T], η =3. The length of the domain L =100L c , and the simulated time T =10 4 t c . The mesh size Δx =1 (i.e., Δx D =L c ), and the time interval Δt =0.1 (i.e., Δt D =10t c ). At the inlet, we set a Cauchy boundary condition D∇cm(0,t)⋅n=q[ce−cm(0,t)],q⋅n<0 ce=δ(t)(24) where c e is the concentration of the inflowing water. The outlet is an open boundary condition D∇cm(L,t)⋅n=0,q⋅n≥0(25) wheren is a unit vector which is normal to the boundary, and ⋅ represents the inner product operator. Initial concentrations are set to zero in both mobile and immobile zones, that is c m (x,0) =c im (x, τ α ,0) =0 [ML −3 ], except for the first mobile node that has initial concentration c m,1 =1 [ML -3 ]. This simulates a pulse injection. In addition, we calculate the cumulative concentration ∫t 0cdt for both mobile and immobile zones, which is equivalent to a continuous injection, simulated through a Cauchy boundary condition at the inlet (Eq. (24)) with an inflow concentration, c e =c m,1 Δxϕ m /2q =0.05. We simulate three cases: a no reaction system (conservative transport), a single reaction system and a sequential reaction system. Unless specified otherwise, models use the same gamma distribution of residence times with β =1/2 and τα 0D=10. In the single reaction system, only one hypothetical reaction A → B occurs, following a first-order rate law with r k1 =k 1 c A . Thus, the degradation of species A follows r A = − k 1 c A , and the formation of species B, r B =k 1 c A . In the sequential reaction system, a second hypothetical reaction B → C is added, also with first-order kinetics, r k2 =k 2 c B . Thus, the degradation of species A follows r A = − k 1 c A , the formation and degradation of species B, r B =k 1 c A −k 2 c B , and the formation of species C, r C =k 2 c B . These are simple linear rate laws. In more realistic models rate laws may be more complex (Carrera et al., 2022). Also, the reaction rate coefficient can be different for each immobile (e.g., immobile zones with long residence times may tend to display more reducing conditions when simulating biochemical reactions). Eq. (7) expresses this dependence by writing k im ( τ α ). Nevertheless, this study focuses on the localization of chemical reactions. Therefore, to facilitate comprehension we used simple reactions with constant rate coefficients (k im ( τ α ) – – – k im ), although our code can handle more complicated reactions with variable rate coefficients. 4. Results and discussions We discuss the results for conservative transport (Section 4.1), which provides the basis for understanding chemical localization in single (Section 4.2), and sequential (Section 4.3) reaction systems. 4.1. No reaction system The impact of MRMT on conservative transport is well-known (e.g., Haggerty and Gorelick, 1995). We summarize here its main features to facilitate understanding the reactive case. Fig. 3 displays concentrations in both mobile and immobile zones in log scale for a pulse injection, and arithmetic scale for continuous injection. The log-log scale is usually preferred for tracer test analysis, because it emphasizes tailing. The late time slope of the pulse BTC is that of the memory function plus 1 (e.g., Haggerty et al., 2000). Therefore, by analyzing the BTC of conservative tracers, modelers can gain insight into the distribution of immobile zones. Still, solute remains in the immobile zones far longer, diffusing slowly back to the mobile zone. The effect of tailing is somewhat less dramatic for a continuous injection (cumulative concentration of a pulse). The breakthrough curve approaches asymptotically the input concentration (c e ) with the bulk of solute arising at t A . For slow exchange (large τα 0D), most of the solute mass pulse flows in the mobile region, so that the peak arrival time coincides with the advective time (red arrow in Fig. 3). Fast exchange (shorter τα 0D) allows equilibrating concentration in the mobile and immobile zones, so that the medium behaves as if the whole porosity was mobile. It is important to notice the above separation between fast and slow is relative to the advective transport time (t Am =Lϕ m /q, where L is the travel distance of solute). If t A ≫ τ α , the immobile region becomes accessible to the solute and the actual travel time reflects the total porosity (t AT =Lϕ/q, where ϕ =ϕ m +ϕ im is the total porosity). This effect is similar to retardation, and it is frequently stated that the impact of immobile zones is to retard the solute (e.g., Haggerty and Gorelick, 1995). In fact, the peak arrival for the τα 0D=10−1 BTC in Fig. 3 (orange arrow) is retarded by a factor of 4, which reflects that ϕ m +ϕ im equals 4ϕ m . Note, however, that (1) this kind of retardation has nothing to do with adsorption or any other partitioning process (in fact, ϕs in our equations can be viewed as retention capacities for sorbing solutes), and (2) the apparent retardation associated to diffusion in immobile zones will be distance and flux dependent. This kind of observations may explain why non-local models are less popular in unconsolidated, granular, aquifers than in fractured media. Fractured media with diffusion lengths of the order of meters will display residence times of the order t D =L 2 /D =(10m) 2 /10 −10 m 2 / s≅10 12 s≅3.10 5 years. That is, diffusion is never exhausted. Diffusion lengths in porous media are less than 1 cm, so that residence times in immobile zones will be of the order of tD= (10−2m)2/10−10m2/s= 106s∼4months. This time may be relevant for tracer tests, where travel time is of a few days, but is too short for natural groundwater flow. Support for this kind of observations is provided by Guimer` a and Carrera (2000). They observed that the “advective porosity” (i.e. the porosity derived from the peak arrival time) calculated from a broad collection of tracer tests in fractured rocks correlates with the peak arrival time. This implies that the apparent mobile porosity increases when the flow rate is reduced, which we take as indicative of fast immobile regions equilibrating with truly mobile zones. Increasing the travel time causes an increasing fraction of immobile regions to equilibrate with mobile zones, thus becoming effectively mobile. Further insight into the behavior of conservative solutes can be gained from Fig. 3 (right), which displays the immobile concentration (color scale) versus distance and immobile residence time for t D =300. The figure might be considered somewhat misleading in that the vertical axis refers to immobile porosity associated to small residence times (recall pdf of immobile zones for a gamma distribution, Fig. 2) and not distance to the mobile zone. However, as shown in Fig. 3 the concentration in immobile regions is independent of their volumetric fraction. What Fig. 3 (right) shows is that concentrations in the immobile zones with short residence times (from τ α D =10 −2 to 10 1 , much shorter than the dimensionless travel time of 300) are identical to concentrations in the mobile zone. On the other hand, concentration is negligible in the immobile zones whose residence time is much larger than travel time, for τ α D higher than 10 3 the concentration in the immobile zones is 0 even at the shorter distances. This might suggest that slow immobile regions might be neglected, whereas fast immobile regions might be lumped with the mobile domain. While this may be appropriate for inert tracers, it may not be suitable for reacting tracers because reactivity is usually higher in the immobile regions, certainly in the case of biofilms (Kone et al., 2014; Seifert and Engesgaard, 2007; Taylor and Jaff´ e, 1990), which motivates our work. 4.2. Single reaction system In this case we consider two solutes, A and B, A is the only solute present in the inflow water while B is produced by the degradation of A. J. Wang et al. Advances in Water Resources 167 (2022) 104286 7 Concentrations of the two species in mobile and immobile zones for a pulse and a continuous injection are shown in Fig. 4 for several reaction times higher than the mean residence time. The situation is now more complicated than in the conservative transport case because three sets of characteristic times are involved (residence times in immobile zones, reaction times, and advective transport time). The first immediate observation, is that maximum concentrations of species A are reduced when the reaction rate is increased (i.e., when the characteristic reaction time is reduced). Note that BTC concentrations of A become negligible when reaction time is much smaller than transport time (i.e., for case of τ rD =10 2 , when transport time t D >10 3 , that means τ rD ≪ t D , species A is almost reduced and cannot be detected). The behavior of the distribution of concentrations in the immobile regions may be less intuitive (see Fig. 5). Concentration of the parent species, A, is negligible in immobile zones with residence time much larger than the reaction time ( τ α D ≫ τ rD ), because the solute diffusing into these zones is degraded before a significant concentration can build up. It is also negligible in zones where residence time is much larger than the simulation time (t D ) because almost no solute has diffused into these regions yet. On the other extreme, concentration in fast immobile Fig. 4. Mobile (solid lines) and average immobile (dash lines) concentrations of species A (left column) and B (right column) versus time at x D =100L c in response to a pulse input (up row) and to a continuous injection (bottom row) for transport in a medium with a gamma distribution of residence times in immobile zones with β =1/2, and τ α 0D=10. Reaction times are uniformly distributed in immobile zones with τ rD =∞, 10 4 , 10 3 and 10 2 . Species B is absent when no reaction occurs ( τ rD =∞). Fig. 5. Cumulative concentration profiles versus dimensionless distance and dimensionless residence times for species A (left) and B (right) at t D =300 (top), and 10000 (bottom) after continuous injection of A in a medium with a gamma distribution of residence times in immobile zones with β =1/2, τ α 0D =10. Reaction A→B occurs with a characteristic time τ rD =10 2 . The white solid lines indicate immobile zones whose residence times equals the reaction time, τ α D = τ rD , and the white dash lines denote immobile zones whose residence times approximate the simulation time, τ α D =t D . J. Wang et al. Advances in Water Resources 167 (2022) 104286 8 regions is virtually equal to the mobile concentration, which suggests that these regions could be lumped with the mobile regions. More interesting is the behavior of daughter species B, which is washed away from the mobile zone and fast immobile zones before it becomes significant. As a result, B concentrations are dominant in immobile zones whose residence time lies between the reaction time and the simulation time; Solute A has arrived at immobile zones with residence time shorter than t D . At those immobile zones with residence time higher than reaction time, A has been transformed into B. This is the fundamental point to demonstrate that B has become localized in certain immobile zones, which motivates this work and is further analyzed for the case of sequential reactions. 4.3. Sequential reaction system The sequential reaction case becomes increasingly complicated by the multiplicity of times involved. Therefore, we start by considering the steady state solution to the sequential reaction system for a generic immobile region that exchanges mass with the mobile zone. The governing equations are 0= α (cmA −cA) − k1cA 0= α (cmB −cB) + k1cA−k2cB 0= α (cmC −cC) + k2cB (26) where we have dropped the dependence on x D ,t D and τ α D for simplicity and because it is self-evident (i.e., c A =c im,A (x D , τ α D ,t D ), c mA =c mA (x D ,t D ), etc.). From Eq. (26), we can obtain the analytical solution of each species at steady state, as a function of the mobile concentrations: cA= α cmA ( α +k1) cB= α cmB +k1cA ( α +k2) cC= α cmC +k2cB α (27) These expressions become clearer when the concentrations of daughter products in the mobile zone can be neglected (i.e., near the input, assuming they are absent in the inflow, and that the advective flux is large) Assuming c mB =c mC =0 in Eq. (27), we get Fig. 6. Steady state concentrations of species A (blue), B (orange), and C (yellow) as a function of residence time in immobile zones, τ α [T], assuming that mobile concentrations of B and C are zero. The first reaction time, τ r1 [T], for A→B transformation equals 10 (solid) and 1 (dash). The reaction time, τ r2 [T], for B→C reaction is 100 (left) and 1 (right). The corresponding reaction times are indicated by arrows. Fig. 7. Breakthrough curves of species A, B and C at L =100L c with one immobile (color lines) zone for τ r1D = τ r2D =10 3 . The black solid line represents results with a mobile zone with ϕ m =0.1 [-] and no reaction (in which case species B and C are never formed). The black dotted lines represent results with a mobile zone with ϕ m =0.4 [-] and τ r1,mD = τ r2,mD =1333 (i.e., k 1, m =k 2,m =0.0075). J. Wang et al. Advances in Water Resources 167 (2022) 104286 9 cA cmA = α ( α +k1) cB cmA = α k1 ( α +k1)( α +k2) cC cmA =k1k2 ( α +k1)( α +k2) (28) Fig. 6 displays the steady state concentrations of species A, B and C in immobile zones versus residence time ( τ α =1/ α ). Note that the sum of concentrations A, B and C equals one, which is the concentration of A in the inflow. This reflects the fact that this sum is conservative, because r A +r B +r C = − k 1 c A +k 1 c A −k 2 c B +k 2 c B =0 and the fact that at steady state the concentration of a conservative species in any immobile zone equals the mobile concentration. As can be observed, species A (blue), which is the only species present in the mobile zone, dominates immobile zones with residence times much smaller than the reaction time of the first reaction (say τ α < τ r1 /10). Therefore, immobile zones with residence times significantly lower than the first reaction characteristic time can effectively be considered as mobile. Species B (orange) dominates immobile zones with residence times between the two reaction times when τ r1 < τ r2 , but its concentration is negligible when τ r1 > τ r2 (i. e., when its decay is faster than its production). Species C (yellow) dominates those immobile zones with residence times larger than the reaction time of the second reaction ( τ α > τ r2 =1/k 2 ). Thus, when the second reaction is also fast (right graph of Fig. 6), the concentration of the intermediate species (B) becomes small, and the two reactions could probably be lumped and species B neglected (provided it does not affect other reactions). Results for the full transport case (A→B and B→C) are summarized in Fig. 7 and Fig. 8. Fig. 7 displays the BTCs for the transport conditions of Fig. 3 and Fig. 4, but only one immobile region to better isolate the effect of residence time. As in Fig. 7, when the exchange is slow, the BTC of species A (green and purple BTCs for τα 0D=104 and 10 3 , respectively) overlaps with that of a conservative solute in the mobile zone (black line), except at the tail. The concentrations of B and C are very small (note the vertical axis). The case of short residence times is also illustrative. The corresponding BTC (blue line) is virtually identical to the BTC (dotted black line) where all the porosity has been lumped as mobile. Obviously, the immobile zone reactions must be accounted for in this mobile (the effective reaction rate is derived in the Appendix A). That is, immobile regions with short residence time (compared to transport and reaction times, i.e., τ α D ≪ t D , and τ α D ≪ τ rD ) can be lumped and treated as mobile. The effective reaction rate in the case of Fig. 7 turns out to be, k am =k im ϕ im /(ϕ im +ϕ m ) (see Appendix A). The most significant observation from the point of view of chemical localization derives from Fig. 8, which displays immobile regions concentrations as a function of space and residence time. This figure makes it clear that C concentrations are only relevant in immobile zones with residence times comparable to reaction times, regardless of the transport time. Therefore, the three sets of times (transport time, reaction time, and residence times) are relevant when deciding the appropriate discretization of residence times. Note that this result is consistent with that of Fig. 6. The second reaction tends to become localized in zones where the residence time is greater than the first reaction time, and the final daughter species, C, tends to concentrate in zones where the residence Fig. 8. Cumulative concentration profiles versus dimensionless distance and dimensionless residence times for species A (left), B (middle) and C (right) for τ r1D = τ r2D =10 3 at t D =50, 200, 500, 2000 and 10000 from top to bottom. The medium is characterized by immobile zones of which the residence time follows a gamma distribution with β =1/2, τ α 0D=10. J. Wang et al.