The scaling of mineral dissolution rates under complex flow conditions

Journal Pre-proofs The scaling of mineral dissolution rates under complex flow conditions Rong Li, Chen Yang, Dongfang Ke, Chongxuan Liu PII: DOI: Reference:

S0016-7037(20)30071-5 https://doi.org/10.1016/j.gca.2020.01.048 GCA 11628

To appear in:

Geochimica et Cosmochimica Acta

Received Date: Revised Date: Accepted Date:

9 September 2019 6 January 2020 23 January 2020

Please cite this article as: Li, R., Yang, C., Ke, D., Liu, C., The scaling of mineral dissolution rates under complex flow conditions, Geochimica et Cosmochimica Acta (2020), doi: https://doi.org/10.1016/j.gca. 2020.01.048

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier Ltd.

The scaling of mineral dissolution rates under complex flow conditions Rong Lia,b, Chen Yanga,b, Dongfang Kea,b, Chongxuan Liu*a,b a

State Environmental Protection Key Laboratory of Integrated Surface Water-

Groundwater Pollution Control, School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen, 518055, China b

Guangdong Provincial Key Laboratory of Soil and Groundwater Pollution Control,

School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen, 518055, China * Corresponding author: Chongxuan Liu (email: [email protected]; postal address: No 1088, Xueyuan Road, Xili, Nanshan District, Shenzhen, Guangdong, China)

Abstract: Mineral dissolution is an important process that provides materials and nutrients to ecosystems. However, a general rate scaling theory has not been developed that can be used to scale mineral dissolution rates from laboratory to field under variable water flow and solute transport conditions. In this study, a mathematical relationship has been derived that can be used to extrapolate the rate of mineral dissolution derived under well-mixed conditions to complex flow systems with spatial heterogeneity and temporal flow velocity transiency. The macroscopic mineral reaction rate in a flow 

system ( rFT ) was derived as rFT = 0 r ( ) f ( )d  adv , where f ( ) is the fluid travel time distribution (FTTD) function, r ( ) is the reaction rate as a function of reaction time in a well-mixed system, and  adv is the mean fluid travel time. Reactive transport simulations were performed to generate various scenarios of mineral dissolution under different flow conditions with permeability heterogeneity and flow velocity transiency using magnesite as an example. The macroscopic mineral dissolution rates calculated from the scaling relationship were compared with the rates averaged from the reactive transport simulations. The results indicate that the averaged magnesite dissolution rates decreased significantly in presence of permeability heterogeneity and under transient flow fluctuation conditions as compared to the rate determined from a homogeneous flow-through system with constant flow rate, showing a complex relationship between the averaged dissolution rate, local dissolution rate, and transient flow characteristics. Remarkably all the averaged rates that differ significantly under different flow conditions were converged to the scaling relationship, providing a consistent way to scale the rate of mineral dissolution from well-mixed reactors to complex flow systems. The results in this study were derived under a steady-state of porosity, permeability, and mineral surface reactivity, and further research is to incorporate their temporal variations.

1. Introduction Mineral dissolution is an important process in natural environment, affecting rock and soil chemical weathering (Anderson et al., 2002; Anderson, 2005; Navarre-Sitchler and Thyne, 2007), global CO2 cycling (Kump et al., 2000; Maher and Chamberlain, 2014, Colbourn et al. 2015), bio-utilization of mineral nutrients (Cole, 1995; Ehrenfeld et al., 2005; Sverdrup, 2009), and the release and transport of mineral-bound contaminants (Bearup et al., 2012; Putnis et al. 2013; Renard et al., 2018). The rates of mineral dissolution can differ significantly in different systems and change with observation scales. Up to 4-5 orders of magnitude difference can occur between the dissolution rates observed at the field and the laboratory scales (Malmstrom et al., 2000; White and Brantley, 2003; Zhu, 2005; Navarre-Sitchler & Brantley, 2007; Maher, 2010; Moore, et al., 2012). The discrepancy in mineral dissolution rate strongly limits the extrapolation of dissolution kinetic models and parameters characterized at the laboratory to natural systems. Extensive research has been performed to understand the rate difference at different scales. Various factors, including spatial and temporal variability of microscopic reactive mineral surface morphology (Levenson & Emmanuel, 2013; Fischer et al., 2014; Fischer & Luttge, 2017), mineral distribution heterogeneity (Li et al., 2008; Li et al., 2014; Molins et al., 2017), microscale solute concentration gradient (Meile & Tuncay 2006; Liu et al. 2015a, 2015b), can all affect the scale-dependent rates of mineral dissolution. In a reactive transport system, these factors are often in coupling with water flow and solute mass transport that collectively affect the mineral dissolution rates at different scales (Steefel, 2008; Battiato et al., 2009; Dentz et al., 2011; Liu et al., 2015a, 2015b; Wen & Li, 2017; Jung & Navarre-Sitchler, 2018; Deng et al., 2018). Advection, dispersion, and molecular-diffusion carry reactants and affect reactant concentration distributions that subsequently affect local reactions, while reactions consume reactants and release reaction products that in turn affect local reactant

concentrations and gradients and thus transport, leading to the strong coupling between transport processes and reactions. Such coupling occurs at the microscopic scale which is often below the resolution of the continuum-scale or macroscopic scale that is typically used to describe geochemical reactive transport in porous media. The macroscopic reaction rates calculated using the reactant concentrations at the macroscopic scale are therefore the coupled results of microscopic scale reactions and transport. The microscopic flow and transport conditions, and their coupling with reactions are therefore important to understand the difference of reaction rates between the microscopic and macroscopic scales, and to develop the relationship of dissolution rates between different scales. In previous studies, effective parameter approaches such as effective reaction rate constant (Luo et al., 2008; Liu et al., 2015b) or effective reactive surface area (Wen & Li, 2017) were often used to link the macroscopic and microscopic reaction rates. The effective parameters, however, have to be estimated from the microscopic distribution of reactants concentrations and reaction rates, which depend on the coupling of microscopic fluid flow, solute transport, and reactions (Liu et al., 2015a; Wen & Li, 2017, 2018; Yu & Hunt, 2017; Hunt & Sahimi, 2017; Deng et al., 2018). The applicability of such effective approaches is limited because the microscopic coupling is often difficult to obtain at the macroscopic scales, especially in heterogeneous porous media. A general approach is therefore needed that can be used to scale reaction rates from the microscopic to macroscopic scales under complex flow conditions. A theoretical scaling equation for the macroscopic mineral dissolution rate was derived mathematically in this study. The approach used a stream tube concept of solute transport and fluid travel time distribution (FTTD) (Simmons et al., 1995; Ginn et al., 1995; Ginn, 2001). To validate the FTTD-based scaling approach, the mineral dissolution processes under variable flow conditions and in physically heterogeneous porous media were first simulated using a reactive transport model at the microscopic scale, which were then used to calculate the macroscopic reaction rates. The calculated

macroscopic mineral dissolution rates were compared with the prediction from the FTTD-based scaling approach without microscopic information. The results indicated the FTTD-based scaling approach can describe macroscopic reaction rates. The approach is general and it can be applied to scale reaction rates to complex flow systems.

2. Derivation of FTTD-based scaling approach Average fluid travel time is often used to generate a Damkohler number Da that may be used to extrapolate a reaction rate to other systems (Maher, 2010; Salehikhoo et al., 2013; Wen et al., 2017; Yu & Hunt, 2017). However, under transient flow condition or in heterogeneous porous media, the average fluid travel time is not fully representative of flow velocity distribution. When a fluid flows through a physically heterogeneous domain, its travel time will depend on flow path that will be affected by local tortuosity and permeability along the flow path. Similarly, when the flow velocity is transient, the fluid travel time will change temporally even when the porous media is homogeneous. Under these conditions, a function of fluid travel time distribution (Simmons et al., 1995; Ginn et al., 1995; Cirpka & Kitanidis, 2000; Ginn, 2001) is therefore needed to fully describe the characteristics of a flow field. A mineral reaction rate depends on water-mineral contact time in batch reactors or fluid residence time (or travel time) in flow-through systems. With increasing contact time or fluid travel time, reactants will be gradually consumed and reaction products accumulated in aqueous phase, leading to the decrease of reaction rate as the reaction system proceeds toward chemical equilibrium. In a batch reactor, a reaction rate can be linked with water-mineral contact time τcontact using the following equation: r(τcontact)= Δc(τcontact) Vwater/(τcontactAm)

(1)

where r (mol/m2/s) is the mineral surface area-normalized averaged mineral reaction rate during the mineral-fluid contact time τcontact (d), Δc(τcontact) (mol/L) is the change of reactant or reaction product concentration during the τcontact, and Vwater (m3) and Am (m2) are the total aqueous phase volume and total reactive mineral surface area. In a batch

reactor, Δc and τcontact can be measured (Cama et al., 2002; Gislason and Oelkers, 2003; White and Brantley, 2003; Gudbrandsson et al., 2014; White et al., 2017), and thus an average reaction rate as a function of contact time can be obtained (Eq. 1). In a flow through system, the τcontact equals the fluid travel time τf in a stream tube if the reactive mineral properties are the same along the stream tube. The reactants/products concentration difference between outflow and inflow (ΔcFT), therefore, equals Δc(τf) as determined from the batch reactor. When there are many stream tubes from the inlet to outlet with each stream tube having temporal flow transiency or spatial heterogeneity in a flow system, the average reactant/product concentration change from the inlet to outlet,

cFT

, can be calculated using a stochastic

transport theory (Simmons et al., 1995; Ginn et al., 1995; Cirpka & Kitanidis, 2000; Ginn, 2001): 

cFT =  c( ) f ( )d 0





0

f ( )d

(2)

where f(τ) (1/d) is the FTTD of the fluid flow from the inlet to outlet of the flow system. The macroscopic average mineral reaction rate, rFT , can then be calculated using Eq. 3: rFT 

where

VFT (Ttotal )

VFT (Ttotal )cFT AmTtotal

(3)

(m3) is the total aqueous phase volume that flows through the system

during a time duration Ttotal (d). The variables related to total time duration, i.e., VFT (Ttotal ) Ttotal , in Eq. 3 can be replaced by average flow rate:

QFT =

VFT (Ttotal ) VPore = Ttotal  adv

where Vpore (m3) is the pore volume and the flow-through system. Defining





0

 adv

(4)

(d) is the average fluid advection time in

f ( )d =1 and by combining Eq. 2, 3 and 4, the

macroscopic mineral dissolution rate in the flow through system can be expressed as:



rFT 

VPore  c( ) f ( )d 0

Am adv

 =



0

r ( ) f ( )d

 adv

(5)

Eq. 5 indicates that the macroscopic (or average) mineral reaction rate in a flowthrough system rFT can be calculated using FTTD function ( f ( ) ) and microscopic mineral reaction rate r ( ) , which can be measured in a well-mixed batch reactor.

3. Simulation models Simulation models with microscopic reaction rates are used to generate various scenarios of mineral dissolution in flow-through systems under different flow conditions. The simulation results were used to calculate average reaction rates, which were then used to compare with the average reaction rates directly scaled from the microscopic reaction rate using Eq. 5 without microscopic reactive transport simulations.

3.1 Intrinsic (microscopic) reaction rate model Magnesite dissolution was used as an example to evaluate the applicability of Eq. 5 in scaling reaction rates. There are three reaction schemes involved in the magnesite dissolution (Table 1). An intrinsic reaction rate expression (local microscopic reaction rate) for the magnesite dissolution has the following expression (Plummer & Wigley, 1976; Chou et al., 1989; Salehikhoo et al., 2013): rMgn  (k1aH  k2 aH CO *  k3 )(1  2

3

IAP ) K eq

(6)

where rMgn (mol/m2/s) is the microscopic magnesite reaction rate (normalized to mineral surface area), k1, k2 and k3 (mol/m2/s) are the reaction rate constants and aH , 

aH CO * are the activities of H+ ion and H2CO3*(carbonic acid) in aqueous phase. Keq is 2

3

the equilibrium constant, and IAP is the ion activity product of the reactants/reaction products in aqueous solution. Table 1 provided the equilibrium constants and intrinsic

rate constants for magnesite dissolution. In modeling the dissolution reaction, H+, HCO3- , Mg2+, Na+ and Cl- were selected as the primary species, and the concentration of secondary species can be calculated from the equilibrium reactions provided in Table 1.

Table 1. Reactions and reaction parameters used in magnesite dissolution model. Mineral kinetic reactions*

Log Keq

k(mol/m2/s)**

MgCO 3 (s)  H   Mg 2+ +HCO 3 -

2.50

6.2  10 5

MgCO 3 (s)  H 2 CO 3*  Mg 2+ +2HCO 3 -

-3.85

5.25 106

MgCO3 (s)  Mg 2+ +CO32-

-7.83

1.0 1010

Aqueous speciation (considered equilibrium reactions)***

a(m2/g)**

1.87

as

H 2 O  H   OH 

-14.00

H 2 CO 3*  H   HCO 3 -

-6.35

HCO 3 -  H   CO 3 2-

-10.33

MgHCO3-  Mg 2+ +HCO3-

-1.04

MgCO3 (aq)  Mg 2+ +CO32-

-2.98

* Plummer & Wigley, 1976; Chou et al., 1989 ** Salehikhoo et al., 2013; Wen & Li, 2017, 2018 *** Wolery et al., 1990

3.2 Coupled model of transport and reactions The coupled model of transport and reactions was used to simulate magnesite dissolution in flow systems. The coupled model can be described using the following reactive transport equation:

ci c  2c  u i  D 2i  Ri t z z

(7)

where ci is the concentration of the i-th primary species, Ri (mol/L/s) is the corresponding microscopic reaction rate that can be calculated using Eq. 6 and parameters in Table 1, D (m2/s) is the dispersion coefficient, u (m/s) is the average pore water velocity that could be calculated using Darcy’s law with homogeneous or randomly generated heterogeneous permeability k. Three scenarios were considered in the simulations. In the first scenario, the magnesite dissolution was assumed to occur in a homogeneous, one dimensional (1D) porous medium under steady state flow conditions with relevant parameters provided in Table 2. 1D model was used here because it is widely used to effectively demonstrate the steady-state or transient coupling of reaction and transport (Liu et al., 2017; Heidari et al., 2017). The length (L), cross section area (A) and porosity of the porous media () are assumed to be L=0.22 m, A=0.01 m2,  =0.35, respectively. The volume fraction of the magnesite in the solid phase of the porous media ( f Mgn ) is assumed to be 10% and longitudinal dispersivity LD is assumed to be 0.001 m.

Table 2. Parameters of the 1D flow-through porous media Parameter

Value

L

0.22 m

A

0.01 m2



0.35

f Mgn

10 %

LD

0.001 m

The chemical species concentrations in the inlet solution are provided in Table 3 and the initial concentrations of the species in aqueous phase were assumed to be in equilibrium with the magnesite in solid phase. The effect of mineral dissolution on the changes in porosity, permeability and mineral surface area in the porous media were

ignored for simplicity. The flow velocity in the simulations was varied, ranging from 0.36 m/d to 36 m/d to evaluate the effects of reactant supply and reaction product removal rates on macroscopic mineral dissolution rate.

Table 3. Chemical species concentrations in the inlet solution Species

Inlet solution

pH

4.0

Dissolved Inorganic Carbon

1.09 105 mol/L

Mg2+

0 mol/L

Na+

8.96 104 mol/L

Cl-

1.0  103 mol/L

In the second scenario, the magnesite dissolution was assumed to occur in a physically heterogeneous 2D porous media (0.44 m × 0.22 m) with a random permeability field generated using sequential Gauss simulations. The 2D systems with random permeability distributions were used to demonstrate the effect of physical heterogeneity (Jung & Navarresitchler, 2018(a), 2018(b); Wen & Li, 2017) on reaction rate scaling and the effectiveness of the scaling equation (Eq. 5). In the sequential Gauss simulations, the permeability was assigned to follow a log-normal distribution and the variance of the log-normal distribution σ2(logK) was 3 in all simulated cases. Anisotropic exponential semivariogram was used to generate correlation structures. The vertical and horizontal correlation length λVertical and λHorizontal were assumed to be 0.16 m and 0.03 m in this study for a demonstration purpose. For comparison, the mineral composition, surface area, porosity and local dispersivity in this scenario were assumed to be identical to those in scenario 1 (Table 2). The only difference between 2D heterogeneous case (scenario 2) and 1D case (scenario 1) was the permeability heterogeneity. In the third scenario, the dynamic flow conditions were enforced in a

homogeneous flow-through system. The inlet water flux was set as periodic square waves with various maximum and minimum velocity and wave period as shown in Fig. 1 to mimic different dynamic flow scenarios. In Fig. 1 umax is the maximum flow velocity, umin is the minimum flow velocity, and uave is the average velocity, T is the time duration of one wave event, Tmax is the time duration when flow velocity is at the maximum, and Tmin is the time duration when flow velocity at the minimum. A new variable Rd was defined as Rd=Tmax/T, which represents the ratio of the time duration when flow is at the maximum vs. the total time duration of the one wave event. The average flow velocity can then be calculated as: uave  umax Rd  umin (1  Rd )

(8)

Eighty-four cases with different Rd, T, and umin were simulated in this study to systematically investigate the influence of dynamic water flow on mineral dissolution rates at the macroscopic scale. For comparison purpose, the average flow velocity in all the simulations were set to 3.6 m/d (see Table S1 in supporting information) by varying umin from 0 (corresponding to the stop/restart flow condition) to 2.7 m/d. With a known average velocity, the maximum velocity can be calculated using Eq. 8.

3.3 Macroscopic reaction rates The macroscopic reaction rates can be calculated using Eq. 5 without microscopic simulations. The macroscopic reaction rate RMgn (mol/d) in 1D column can also be calculated using the microscopic simulation results from the coupled model of transport and reaction using Eq. 9: RMgn 

nMgn Ttotal



Across {

Ttotal

0

L

(cout  cin )udt   [c(Ttotal , x)  c(0, x)]dx} 0

Ttotal

(9)

where nMgn is the accumulative mass of dissolved magnesite (mole) in the flowthrough system during total time duration, Ttotal , and c is the aqueous concentration of the dissolution product (i.e., Mg2+ for this case). When Ttotal is long enough, the second

term in the right hand side is negligible as compared to the first term. RMgn can thus be calculated using effluent Mg2+ concentration as a function of time: Across 

RMgn 

Ttotal

0

(cout  cin )udt Ttotal

(10)

When normalized to total magnesite surface area AMgn (m2) in the flow-through system, equation 10 becomes:

rMgn =

RMgn AMgn

 =

Ttotal

0

(cout  cin )udt

Lf Mgn  Mgn aMgnTtotal

(11)

where  Mgn is magnesite density (g/m3), aMgn is magnesite specific surface area (m2/g), other parameters are defined before. Under a steady state flow and dissolution condition in homogeneous porous media, Eq. 11 can be simplified to: rMgn =

(cout  cin )u Lf Mgn  Mgn aMgn

(12)

3.4 Fluid travel time distribution (FTTD) The scaling equation (Eq 5) requires fluid travel time distribution (FTTD) function, f ( ) . Under steady-state flow condition with or without physical heterogeneity, if a

single tracer pulse input Cin =M Tr (t ) was applied at inlet (  (t ) represents Dirac delta function), f ( ) can be calculated using the time dependent tracer concentration in outlet solution as Eq. 13 (Levenspiel 1999; Pawlowski et al. 2018): f ( ) 

QV CTr ( ) M Tr

(13)

where QV (m3/d) is the volumetric flow rate in the flow through system, CTr ( ) (mol/L) is the average tracer concentration in the outlet solution at time  and M Tr (mol) is the accumulative mass of tracer input. Under transient flow conditions, however, the fluid that enters the simulation

domain at different times may have different fluid travel times. With known temporal profiles of inlet flow velocity u(t) for a given 1D column, the FTTD function, f(t), can be calculated as described in appendix A. With the calculated FTTD, the macroscopic mineral reaction rate can be estimated using Eq. 5.

4. Results 4.1 Dissolution rate under steady state conditions Fig. 2(a) shows the macroscopic (average) magnesite dissolution rate ( rMgn in Eq. 12) calculated using the simulated effluent Mg2+ concentrations cout in the 1D flowthrough system under steady state flow conditions ( cin  0 in this case, Table 3). As water flow velocity increased from 0.9 m/d to 36 m/d, cout decreased from 2.86 × 10-4 mol/L to 1.0 × 10-4 mol/L and rMgn increased from 3.7 × 10-11 mol/m2/s to 5.2 × 10-10 mol/m2/s. The results reflected the effect of fluid travel time on the macroscopic dissolution rate. At a higher flow rate, travel time was shorter, reactant (H+) supply and removal of reaction production (Mg2+) were faster, so that the reactant (H+) concentration was higher and the reaction system was far away from the equilibrium state in the system, leading to a faster reaction rate, but a lower concentration of Mg2+ as it was quickly removed from the system. In a 2D system, 10 heterogeneous permeability fields that were generated randomly were used to simulate mineral dissolution under four different flow conditions with an average fluid velocity ranging from 1.2 m/d to 7.2 m/d (Fig. 3). The flow fields that under randomly generated permeability distribution contain preferential flow paths (Fig. 3a) where fluid travel time is shorter. Along these paths with shorter travel times, the pH changed less, and was close to that in the influent, and reaction product (Mg2+) concentration was lower. Consequently, magnesite dissolved faster in these preferential paths with the same reasons as described for the 1D system. In the

low permeability zone, the flow was nearly stagnant and the fluid travel time is long so that solution chemical composition was nearly in equilibrium with the magnesite. This resulted in a slow magnesite dissolution rate as it was constrained by the local chemical equilibrium. Because of the strong heterogeneity in flow, solute concentration, and local reaction rate distribution, the averaged macroscopic mineral dissolution rate in the simulation domain depended on permeability distribution and flow velocity (Fig. 3d). In general, the average magnesite dissolution rate in the heterogeneous porous media were slower than that in the homogeneous media (Fig. 9d). The difference can reach up to 55 % for the studied cases. Fig. 4(a) shows the FTTD functions calculated using Eq. 13 for the 2D steadystate flow-through cases. The FTTD function describes the probability of fluid travel times or water-mineral contact times along different flow paths. To use FTTD-based scaling approach (Eq. 5), microscopic mineral reaction rate r ( ) as a function of fluid travel time is also needed. Function r ( ) can either be measured from batch experiments or calculated from the kinetic simulations in a well-mixed reactor using intrinsic reaction parameters (Table 1). The calculated r ( ) was plotted in Fig. 4(b), which was then used to calculate the FTTD-based average reaction rate using Eq. 5. The FTTD-calculated reaction rates well matched with those calculated from the reactive transport simulations as shown in Fig. 4(c) with a relative error of <3%, indicating the effectiveness of the FTTD-based scaling approach. The FTTD-based scaling approach was also tested in a flow system with different size to evaluate the potential domain size effect. The total length of the simulation domain was 2.2 m that is comparable with the vertical size of groundwater flow or weathering front region in the field (Moore et al., 2012; Liu et al., 2017; Heidari et al, 2017). For this size of the domain, 5 heterogeneous permeability fields were generated randomly (Fig. 5a) and flow rate was 1.2 m/d and 12 m/d. rMgn and FTTD were calculated in the subdomains with different lengths ranged from 0.04 to 2.2 m (from the inlet). The results showed that the peaks in FTTD functions were delayed with

increasing subdomain length in general and rMgn decreased with subdomain size consequently. The calculated rate rMgn in subdomains again well matched with the rate calculated from FTTD-based model (Fig. 5c), further showing the effectiveness of the scaling equation. The potential effect of FTTD uncertainty on the effectiveness of the rate-scaling approach was also evaluated. The uncertainty of FTTD was generated by decreasing the number (Ns) of tracer samples in numerically sampling the tracer breakthrough curve simulated using the microscopic transport model. The baseline FTTD was generated using 1550 tracer samples (Fig. 5c and Fig. S3) and other FTTDs were generated with NS=25 to 388 (Fig. S4). Generally, the lower tracer sampling frequency resulted in higher FTTD uncertainty. The rFT,Mgn values calculated using Eq. 5 and these FTTDs indicated that when the sampling frequency decreased from NS=1550 to NS=97, the deviation between rFT,Mgn and

rMgn did not change significantly (Fig.

6b). However, when NS decreased from 97 to 25, the deviation increased significantly. This is because when the tracer sample number became too low, the estimated FTTD was no longer able to capture all the information of flow heterogeneity in the system, which subsequently led to the high deviation between the scaled reaction rate and that calculated using the microscopic reactive transport model.

4.2 Dissolution rate under transient flow conditions The effluent Mg2+ concentration changed dynamically under periodic flow velocity conditions in the 1D flow-through system (Fig. 7a). In general, Mg2+ concentration was higher during the lower velocity periods due to the long fluid travel time, and vice versa. The concentration change, however, lagged the change in flow velocity (Fig. 1), leading to the gradual change in effluent concentration because of the effect of dilution and local kinetics of mineral dissolution. With increasing time duration for each fluctuation wave (T), the lag effect normalized to T decreased. When

T was long enough, the lag effect nearly disappeared. Fig. 7(a) and (b) showed that the magnitude of the changes in effluent Mg2+ concentration increased with increasing difference between high and low flow velocity values u . This is because the change in fluid residence time became larger with increasing u and thus the effluent concentration difference increases with increasing changes in flow velocity. Under the same average flow velocity condition, the average rate of magnesite dissolution in the system, rMgn , calculated from the reactive transport simulation was lower under transient condition than that under steady-state condition (Fig. 8 vs. Fig. 2). In addition, the simulation results of rMgn decreased with increasing T (Fig. 8a) and with increasing difference between the maximum and minimum flow velocities within the T (Fig. 8b). The results showed that the average rate of dissolution rMgn varied significantly and was affected by many factors. The macroscopic mineral dissolution rate rMgn

under 1D steady-state flow

condition was a function of flow velocity uave only. (Fig. 2(a)). Under transient flow condition, however, the macroscopic mineral reaction rate is not only a function of average flow velocity, but also depends on the temporal fluctuation characteristics of flow. The average dissolution rate rMgn

increased with decreasing flow velocity

difference Δu. When Δu approached to 0 (steady-state flow condition), rMgn reaches the maximum (Fig. 8b). The result indicated that the average reaction rate determined under steady state flow condition is the upper limit for those under transient flow conditions with the same average flow velocity. The results above were calculated from the microscopic simulations. To test the applicability of the FTTD-based scaling approach for calculating the average dissolution rate, the FTTD function

f ( ) under transient flow conditions was

calculated using Eq. A-1 – A-4 in Appendix A (Fig. 9). The calculations involved

accumulative effluent QE(t) (Eq. A-2) and fluid travel time at different time  (t ) (Eq. A-3). These intermediate variables were also provided in Fig. 9(a) and (b). As shown in Fig. 9(b), with shorter T, the fluid travel time was almost unchanged during the flow fluctuation. Contrarily, with longer T，the fluid travel time showed strong fluctuation. The travel time distribution function f ( ) (Fig. 9(c)) for the specific cases studied here show multiple sharp peaks because of the repeated sharp change in flow velocity from umin to umax, and from umax to umin. The accumulative probability F ( ) , which can be calculated from the travel time distribution function (Eq. A-5) was also provided in Fig. 9(d). Using the calculated f ( ) (Fig. 9(c)) and microscopic reaction kinetics (Fig. 4(b)), the FTTD based average mineral reaction rate, rFT,Mgn , can be calculated (Fig. 10). The FTTD-based rates matched well with those calculated from microscopic reactive transport simulations rMgn with a maximum deviation of less than 8 %. The result indicated that the FTTD-based approach can be used to scale reaction rate under transient flow conditions.

5. Discussion This study demonstrated that the average or macroscopic rate of the mineral dissolution could be significantly higher in homogeneous than that in heterogeneous porous media with the same microscopic reaction rate under the same flow-through distance and average flow velocity (Fig. 3(d)). The decrease in macroscopic reaction rate in the heterogeneous porous media agreed qualitatively with previous simulation results of mineral dissolutions (Jung & Navarre-Sitchler, 2018a; Pandey & Rajaram, 2016; Liu et al., 2015a) and the results for other reactions in porous media (Ginn et al., 1995; Liu et al., 2013). However, a general method to scale reaction rate from a wellmixed system to heterogeneous porous media has not been fully developed. Two major

difficulties involved in developing such a method. One difficulty is to characterize the heterogeneity of a porous domain, especially at the microscopic scale in field that affects local flow velocity and reaction rates. The other difficulty is to quantitatively link microscopic flow and local reaction rate to determine the average or macroscopic reaction rate without direct reactive transport simulations at the microscopic level in heterogeneous porous media. The FTTD-based scaling approach developed here eliminates the need of performing reactive transport simulations at the microscopic level so that the microscopic characterization of porous media is not necessary. The tradeoff is that it requires the estimation of FTTD function in the porous media, which, as shown in this study, can be determined from the macroscopic behavior of tracer at effluent. Under field conditions, tracer concentrations and thus transport may be determined at the macroscopic level using various approaches (Uhlenbrook et al, 2002; Rinaldo et al., 2011; Mccallum et al., 2014). Another requirement of using the FTTDbased approach is to determine microscopic reaction rate as a function of fluid travel time, which can be readily determined from a well-mixed reactor or from kinetic simulations as shown in this study. Steady-state flow conditions were often used to explore the scaling effect on mineral dissolution rates (Simmons et al., 1995; Ginn et al., 1995; Cirpka & Kitanidis, 2000; Ginn, 2001). The results in this study showed that the flow transiency has a significant effect on the macroscopic mineral reaction rate. Under field conditions, the flow velocity fluctuations with variable fluctuating amplitudes and frequency are the common phenomenon as a result of rain water infiltration, groundwater-surface water exchanges, etc (Mayer, 1998; Li & Jiao, 2002; Ma et al., 2014; Schaefer et al., 2016; Yang et al., 2018; Huang et al., 2018). As demonstrated in this study, the average rate under transient flow condition was always lower than that under steady-state condition under the same average flow velocity conditions, and thus the dissolution rate under the steady state condition can be used as the upper limit of average reaction rates under transient flow conditions. The influence of flow transiency on the average reaction rate

was well predicted from its microscopic rate using the FTTD-based approach. In the field, the temporal variation of fluid flow may be more complicated as a result of the randomness of infiltration and water table fluctuation. Dynamic information of fluid flow, together with tracer measurements is required to characterize the FTTDs (Botter et al., 2011; Smith et al., 2018). Using FTTD based scaling approach, we could revisit the mechanism of reaction rate decrease due to the physical heterogeneity. As mentioned above, the local reaction rate increased in the fast flow paths (shorter travel time) and decreased in slower flow zones (longer travel time). The average reaction rate depends on their relative importance in the system. However, as shown in Fig. 4(a), the major peaks of f ( ) appeared in short travel time zone in most of FTTD curves. This indicated that major fraction of flowing water passed through the preferential paths. This part of flowing water can only contact with the reactive minerals in the preferential flow zones with insufficient contact time (decreased c ( ) value in Eq. 5 as shown in Fig. 4(b)). On the other hand, the mineral surface area in the low permeability zones could not effectively contribute on the overall reaction rate because of the slow water flow and mass transfer. The net result is that the physical heterogeneity resulted in the decrease of average reaction rate. The lower average reaction rate under transient flow condition could also be explained using Eq. 5. With velocity fluctuation, the travel time decreased under higher velocity stage, and vice versa. Although both the higher and lower velocity periods contributed to the overall average reaction rate, the decrease of c ( ) (or equivalently,

 r ( ) in Eq. 5) under shorter  , as shown in Fig. 4(b), was more than that can be compensated by the increase of c ( ) under longer  . The net result was the decrease of the average reaction rate. As mentioned in section 4.2, under a shorter fluctuation period T, the FTTD function was close to that under steady-state conditions. Thus the effect of flow transiency under short T has a less effect on the average rate. On the other hand, with a longer T, the flow becomes more transient and thus the average reaction

rate decreases more. A Dankoler number Da is often used to correlate the fluid travel time and average reaction rate (Maher, 2010; Salehikhoo et al., 2013). Such a correlation may be used to scale reaction rate if Da is known. Fig. 11 showed the average reaction rate estimated using the FTTD approach as a function of average Da (Fig. 11). The result indicated that the scaling approach can qualitatively match with observed change trend of mineral weathering rates from laboratory (Salehikhoo et al., 2013) to field scales (Maher, 2010; Wen & Li, 2018). More specifically, with a longer fluid residence time (Da>2 in this study), the FTTD-calculated reaction rates approached to a straight time, which is for transport-limit reaction cases ( rlim  csat uAcross AMgn  Da 1 , see Fig. 11) because under a long fluid travel time, the downstream of the flow-domain is near chemical equilibrium and transport controls reaction rates. With a shorter fluid travel time, however, the FTTD-calculated reaction rates are larger, showing that it increases with Da with a smaller slope ( r  Da  0.645 in this study). This is because both transport and reaction kinetics control the overall reaction rates (Salehikhoo et al., 2013; Wen et al., 2018). In general, transition from chemical equilibrium limitation to the kinetic limitation that was related to water-mineral contact time τ was quantitatively incorporated in term r(τ) or c(τ) in Eq. 5. Consequently, the FTTD-weighted average of r(τ) can capture the transition from transport limitation to the kinetic limitation. Due to the impact of physical heterogeneity and flow transiency discussed in this study, the average reaction rate could show significant divergence under identical Da number in the transport-kinetics mixed control region as shown in Fig. 11. In this case, the FTTD based scaling approaches could predict the average reaction rate better than the empirical correlations based on average Da number alone (this study; Salehikhoo et al., 2013; Wen et al., 2018). The accuracy of the macroscopic reaction rate calculated using the FTTD-based approach largely depends on the accuracy of the FTTD. In field, one of the major factors that may affect the accuracy of the estimated FTTD is the frequency of tracer sampling.

Our results indicated that the accuracy of the FTTD estimated from tracer measurements decreased with decreasing sampling frequency, which subsequently increased the uncertainty of the reaction rates calculated using the FTTD-based approach (Eq. 5). Therefore, sampling frequency and strategy needs to be carefully considered in tracer test for deriving FTTD. In this study, the FTTD-based approach was verified only in presence of spatial and temporal complexity of flow field. The geochemical reactions may lead to the changes in porosity, permeability and mineral surface area, and subsequently influence flow travel time distribution and macroscopic reaction rates (Molins et al., 2014; Szymczak & Ladd, 2014; Soulaine et al., 2017). The effects of the changes of physical properties (porosity and permeability change) on the macroscopic reaction rates could be captured using time-variable FTTDs that can be determined based on tracer concentrations measured at different times. The scaling approach (Eq. 5) can thus be used for such cases. The effects of the changes of chemical properties such as reactive mineral surface area were not considered in this study. Previous numerical studies (Li et al.; 2014; Wen & Li, 2017; Jung & Navarre-Sitchler, 2018b; Pandey & Rajaram, 2016) indicated that the spatial heterogeneity and temporal transiency of reactive mineral surface area can have a strong effect on average or macroscopic reaction rates. Thus, further studies are needed to incorporate the effects of chemical heterogeneity into the FTTD-based approach.

6. Conclusions In this study, a FTTD-based scaling approach was established for scaling reaction rates from the microscopic to macroscopic systems under complex flow conditions and in physically heterogeneous porous media. The average reaction rates calculated using FTTD-based approach were compared with those calculated from microscopic simulations. The result indicated that the FTTD-based approach can be used to scale

the reaction rate under conditions of physical heterogeneity and flow transiency. Although the magnesite dissolution was used as an example for demonstrative purpose, the approach is general and can be used for other reaction systems without need of performing microscopic reactive transport simulations. The trade-off is the requirement to determine FTTD function. This is effectively transform the reaction rate scaling to a problem of determining FTTD function. In this study, the FTTD function was calculated using simulated tracer effluent concentrations in flow-through systems. Under field conditions, the determination of FTTD function will require tracer tests or using natural tracer concentration information. The applicability of the FTTD-based approach therefore depends on the feasibility to find the FTTD function. The FTTD-based scaling approach as proposed in this study considered physical heterogeneity and flow transiency. Further research is needed to incorporate chemical heterogeneity into the FTTD-based approach for scaling reaction rates.

Supporting Information Table S1 Shows the input parameters for transient flow simulation cases. Fig. S1 shows the effluent Mg2+ concentrations as a result of magnesite dissolution in the column with variable flow velocity. Fig. S3 shows the supplementary results of tracer transport simulation and magnesite dissolution simulation for Fig. 5. Fig. S4 shows the FTTDs generated with different tracer sample number (NS). ACKNOWLEDGEMENT The funding support was provided by National Natural Science Foundation of China (project No. 41830861, 41907166). Additional support was from the Program for Guangdong Introducing Innovative and Entrepreneurial Teams (2017ZT07Z479) and Southern University of Science and Technology (Grant No. G01296001).

Appendix A: The FTTD calculation under transient flow condition Neglecting the dispersion, the accumulative effluent volume at time t in 1D

column can be calculated as: t

QE (t )  ACross  u (t ')dt ' 0

(A-1)

For the fluid entering the column at time t and leaving the column at t1, the fluid travel time is:

 (t )  t1  t

(A-2)

where the leaving time t1 can be calculated by solving the following equation: t1

QE (t1 )  QE (t )  ACross  u (t ')dt '  VPore  ACross L t

(A-3)

where Vpore is the pore volume of the entire column. The FTTD function f ( ) is the probability density of travel time . An important point is that the statistical counting procedure of f ( ) needs to be sampled over effluent water volume ΔQE instead of time Δt. The mathematical definition of f ( ) is f ( ) 

lim  0

QE 0

N (QE    '     ) N (QE    '     )  = lim  N Total QE,Total / QE  0

(A-4)

QE 0

where, N (QE    '     ) is the counting number of the samples (sample volume is QE ) that has fluid travel time between  and    , N Total is the total sample number and QE,Total is the total effluent volume of the samples. In addition with the probability density f ( ) , we can also define the accumulative probability function as Eq. A-5. 

F ( )   f ( )d  0

(A-5)

where F ( ) is the probability of fluid travel time less than a given value  , and it is also known as accumulative form of FTTD. Because f ( ) calculated in this study has sharp peaks at certain 

values (mathematically needs to be expressed as 

function), involving an accumulative function F ( ) may provide a clearer description of FTTDs. To apply the counting process of discretized  space with   103

f ( ) and F ( ) numerically, we

VPore as a sampling time step. AcrossVave

References Anderson, S. P., Dietrich, W. E., & Brimhall, G. H. (2002). Weathering profiles, massbalance analysis, and rates of solute loss: Linkages between weathering and erosion in a small, steep catchment. Geological Society of America Bulletin, 114(9), 1143-1158. Anderson, S. P. (2005). Glaciers show direct linkage between erosion rate and chemical weathering fluxes. Geomorphology, 67(1), 147-157. Battiato, I., Tartakovsky, D. M., Tartakovsky, A. M., & Scheibe, T. D. (2009). On breakdown of macroscopic models of mixing-controlled heterogeneous reactions in porous media. Advances in Water Resources, 32(11), 1664-1673. Bearup, L. A., Navarresitchler, A. K., Maxwell, R. M., & Mccray, J. E. (2012). Kinetic metal release from competing processes in aquifers. Environmental Science & Technology, 46(12), 6539-6547. Botter, G., Bertuzzo, E., & Rinaldo, A. (2011). Catchment residence and travel time distributions: The master equation. Geophysical Research Letters, 38(11). Cama, J., Metz, V., & Ganor, J. (2002). The effect of pH and temperature on kaolinite dissolution rate under acidic conditions. Geochimica et Cosmochimica Acta, 66(22), 3913-3926. Chou L., Garrels R. M. and Wollast R. (1989) Comparative study of the kinetics and mechanisms of dissolution of carbonate minerals. Chem. Geol. 78, 269–282. Cirpka, O. A., & Kitanidis, P. K. (2000). An advective ‐ dispersive stream tube approach for the transfer of conservative‐tracer data to reactive transport. Water Resources Research, 36(5), 1209-1220.Colbourn, G., Ridgwell, A., & Lenton, T. M. (2015). The time scale of the silicate weathering negative feedback on atmospheric CO2. Global Biogeochemical Cycles, 29(5), 583-596. Cole, D. W. (1995). Soil nutrient supply in natural and managed forests. Plant and Soil, 62, 43-53 Deng, H., Molins, S., Trebotich, D., Steefel, C. I., & Depaolo, D. J. (2018). Pore-scale numerical investigation of the impacts of surface roughness: Upscaling of reaction rates in rough fractures. Geochimica et Cosmochimica Acta, 239(15), 374-389.

Dentz, M., Borgne, T. L., Englert, A., & Bijeljic, B. (2011). Mixing, spreading and reaction in heterogeneous media: A brief review. Journal of Contaminant Hydrology, 120(1), 1-17. Ehrenfeld, J. G., Ravit, B., & Elgersma, K. J. (2005). FEEDBACK IN THE PLANTSOIL SYSTEM. Annual Review of Environment and Resources, 30(1), 75-115. Fischer, C., Kurganskaya, I., Schafer, T., & Luttge, A. (2014). Variability of crystal surface reactivity: What do we know? Applied Geochemistry, 43, 132-157. Fischer, C., & Luttge, A. (2017). Beyond the conventional understanding of water–rock reactivity. Earth and Planetary Science Letters, 457, 100-105. Ginn, T. R., Simmons, C. S., & Wood, B. D. (1995). Stochastic-Convective Transport with Nonlinear Reaction: Biodegradation with Microbial Growth. Water Resources Research, 31(11), 2689-2700. Ginn, T. R. (2001). Stochastic-convective transport with nonlinear reactions and mixing: finite streamtube ensemble formulation for multicomponent reaction systems with intra-streamtube dispersion. Journal of Contaminant Hydrology, 47(1), 1-28. Gislason, S. R., & Oelkers, E. H. (2003). Mechanism, rates, and consequences of basaltic glass dissolution: II. An experimental study of the dissolution rates of basaltic glass as a function of pH and temperature. Geochimica et Cosmochimica Acta, 67, 3817-3832. Gudbrandsson, S., Wolffboenisch, D., Gislason, S. R., & Oelkers, E. H. (2014). Experimental determination of plagioclase dissolution rates as a function of its composition and pH at 22°C. Geochimica et Cosmochimica Acta, 139, 154-172. Heidari, P., Li, L., Jin, L., Williams, J. Z., & Brantley, S. L. (2017). A reactive transport model for Marcellus shale weathering. Geochimica et Cosmochimica Acta, 217, 421-440. Huang, K., Liu, Y., Yang, C., Duan, Y., Yang, X., & Liu, C. (2018). Identification of Hydrobiogeochemical Processes Controlling Seasonal Variations in Arsenic Concentrations Within a Riverbank Aquifer at Jianghan Plain, China. Water Resources Research, 54(7), 4294-4308. Hunt, A. G., & Sahimi, M. (2017). Flow, Transport, and Reaction in Porous Media:

Percolation Scaling, Critical-Path Analysis, and Effective Medium Approximation. Reviews of Geophysics, 55(4), 993-1078. Jung, H., & Navarresitchler, A. K. (2018a). Physical heterogeneity control on effective mineral dissolution rates. Geochimica et Cosmochimica Acta, 227, 246-263. Jung, H., & Navarre-Sitchler, A. K. (2018b). Scale effect on the time dependence of mineral dissolution rates in physically heterogeneous porous media. Geochimica et Cosmochimica Acta, 234, 70-83. Kump, L. R., Brantley, S. L., & Arthur, M. A. (2000). Chemical Weathering, Atmospheric CO2, and Climate. Annual Review of Earth and Planetary Sciences, 28(1), 611-667. Levenson, Y., & Emmanuel, S. (2013). Pore-scale heterogeneous reaction rates on a dissolving limestone surface. Geochimica et Cosmochimica Acta, 119, 188-197. Levenspiel, O. Chemical reaction engineering, third ed., John Wiley & Sons, United States of America, 1999, pp. 257–320, http://dx.doi.org/10.1021/ie990488g. Li, H., & Jiao, J. J. (2002). Analytical solutions of tidal groundwater flow in coastal two-aquifer system. Advances in Water Resources, 25(4), 417-426. Li, L., Peters, C. A., & Celia, M. A. (2007). Effects of mineral spatial distribution on reaction rates in porous media. Water Resources Research, 43(1). Li, L., Steefel, C. I., & Yang, L. (2008). Scale dependence of mineral dissolution rates within single pores and fractures. Geochimica et Cosmochimica Acta, 72(2), 360377. Li, L., Salehikhoo, F., Brantley, S. L., & Heidari, P. (2014). Spatial Zonation Limits Magnesite Dissolution in Porous Media. Geochimica et Cosmochimica Acta, 126, 555-573. Liu, C., Shang, J., Kerisit, S. N., Zachara, J. M., & Zhu, W. (2013). Scale-dependent rates of uranyl surface complexation reaction in sediments. Geochimica et Cosmochimica Acta, 105, 326-341. Liu, Y., Liu, C., Zhang, C., Yang, X., & Zachara, J. M. (2015a). Pore and continuum scale study of the effect of subgrid transport heterogeneity on redox reaction rates. Geochimica et Cosmochimica Acta, 163, 140-155. Liu, C., Liu, Y., Kerisit, S. N., & Zachara, J. M. (2015b). Pore-Scale Process Coupling and Effective Surface Reaction Rates in Heterogeneous Subsurface Materials. Reviews in Mineralogy & Geochemistry, 80(1), 191-216.

Liu, Y., Xu, F., & Liu, C. (2017). Coupled Hydro-Biogeochemical Processes Controlling Cr Reductive Immobilization in Columbia River Hyporheic Zone. Environmental Science & Technology, 51(3). Luo, J., Dentz, M., Carrera, J., & Kitanidis, P. K. (2008). Effective reaction parameters for mixing controlled reactions in heterogeneous media. Water Resources Research, 44(2). Ma, R., Zheng, C., Liu, C., Greskowiak, J., Prommer, H., & Zachara, J. M. (2014). Assessment of controlling processes for field-scale uranium reactive transport under highly transient flow conditions. Water Resources Research, 50(2), 10061024. Maher, K. (2010). The dependence of chemical weathering rates on fluid residence time. Earth and Planetary Science Letters, 294, 101-110. Malmstrom, M., Destouni, G., Banwart, S. A., & Stromberg, B. (2000). Resolving the scale-dependence of mineral weathering rates. Environmental Science & Technology, 34(7), 1375-1378. Mccallum, J. L., Engdahl, N. B., Ginn, T. R., & Cook, P. G. (2014). Nonparametric estimation of groundwater residence time distributions: What can environmental tracer data tell us about groundwater residence time? Water Resources Research, 50(3), 2022-2038. Meile, C., & Tuncay, K. (2006). Scale dependence of reaction rates in porous media. Advances in Water Resources, 29(1), 62-71. Mitchell, R. J., & Mayer, A. S. (1998). A numerical model for transient-hysteretic flow and solute transport in unsaturated porous media. Journal of Contaminant Hydrology, 30 (3-4), 243-264. Molins, S., Trebotich, D., Yang, L., Ajofranklin, J. B., Ligocki, T. J., Shen, C., & Steefel, C. I. (2014). Pore-scale controls on calcite dissolution rates from flowthrough laboratory and numerical experiments. Environmental Science & Technology, 48(13), 7453-7460. Molins, S., Trebotich, D., Miller, G. H., & Steefel, C. I. (2017). Mineralogical and transport controls on the evolution of porous media texture using direct numerical simulation. Water Resources Research, 53(5), 3645-3661. Moore, J., Lichtner, P. C., White, A. F., & Brantley, S. L. (2012). Using a reactive transport model to elucidate differences between laboratory and field dissolution

rates in regolith. Geochimica et Cosmochimica Acta, 93, 235-261. Navarre-Sitchler, A. K., & Brantley, S. L. (2007). Basalt weathering across scales. Earth and Planetary Science Letters, 261(1-2), 321-334. Pandey, S., & Rajaram, H. (2016). Modeling the influence of preferential flow on the spatial variability and time-dependence of mineral weathering rates. Water Resources Research, 52(12), 9344-9366. Pawlowski, S., Nayak, N., Meireles, M., Portugal, C. A., Velizarov, S., & Crespo, J. G. (2018). CFD modelling of flow patterns, tortuosity and residence time distribution in monolithic porous columns reconstructed from X-ray tomography data. Chemical Engineering Journal, 350, 757-766. Plummer L. N. and Wigley T. M. L. (1976) The dissolution of calcite in CO2-saturated solutions at 25 oC and 1 atmosphere total pressure. Geochim. Cosmochim. Acta 40, 191–202. Putnis, C. V., Renard, F., King, H. E., Monteshernandez, G., & Ruizagudo, E. (2013). Sequestration of selenium on calcite surfaces revealed by nanoscale imaging. Environmental Science & Technology, 47(23), 13469-13476 Rinaldo, A., Beven, K., Bertuzzo, E., Nicotina, L., Davies, J., Fiori, A., Russo, D & Botter, G. (2011). Catchment travel time distributions and water flow in soils. Water Resources Research, 47(7). Salehikhoo, F., Li, L., & Brantley, S. L. (2013). Magnesite dissolution rates at different spatial scales: The role of mineral spatial distribution and flow velocity. Geochimica et Cosmochimica Acta, 108, 91-106. Schaefer, M. V., Ying, S. C., Benner, S. G., Duan, Y., Wang, Y., & Fendorf, S. (2016). Aquifer Arsenic Cycling Induced by Seasonal Hydrologic Changes within the Yangtze River Basin. Environmental Science & Technology, 50(7), 3521-3529. Simmons, C. S., Ginn, T. R., & Wood, B. D. (1995). Stochastic-convective transport with nonlinear reaction: Mathematical framework. Water Resources Research, 31(11), 2675-2688. Smith, A., Tetzlaff, D., & Soulsby, C. (2018). On the Use of StorAge Selection Functions to Assess Time ‐ Variant Travel Times in Lakes. Water Resources Research, 54(7), 5163-5185. Soulaine, C., Roman, S., Kovscek, A. R., & Tchelepi, H. A. (2017). Mineral dissolution and wormholing from a pore-scale perspective. Journal of Fluid Mechanics, 827,

457-483. Steefel. C.I. (2008) Geochemical kinetics and transport. In S. L. Brantley, J. D. Kubicki, A. F. White (Eds.), Kinetics of Water–Rock Interaction (pp. 539-584). New York: Springer Sverdrup, H. (2009). Chemical weathering of soil minerals and the role of biological processes.. Fungal Biology Reviews, 23(4), 94-100. Szymczak, P., & Ladd, A. J. (2014). Reactive-infiltration instabilities in rocks. Part 2. Dissolution of a porous matrix. Journal of Fluid Mechanics, 738, 591-630 Uhlenbrook, S., Frey, M. M., Leibundgut, C., & Maloszewski, P. (2002). Hydrograph separations in a mesoscale mountainous basin at event and seasonal timescales. Water Resources Research, 38(6), 31-1-31-14. Wen, H., & Li, L. (2017). An upscaled rate law for magnesite dissolution in heterogeneous porous media. Geochimica et Cosmochimica Acta, 289-305. Wen, H., & Li, L. (2018). An upscaled rate law for mineral dissolution in heterogeneous media: The role of time and length scales. Geochimica et Cosmochimica Acta, 235, 1-20. White, A. F., & Brantley, S. L. (2003). The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field?. Chemical Geology, 202, 479-506. Wolery T. J., Jackson K. J., Bourcier W. L., Bruton C. J., Viani B. E., Knauss K. G. and Delany J. M. (1990) Current status of the EQ3/6 software package for geochemical modeling. ACS Symp. Ser. 416, 104–116. Yang, C., Zhang, Y., Liu, Y., Yang, X., & Liu, C. (2018). Model‐Based Analysis of the Effects of Dam-Induced River Water and Groundwater Interactions on HydroBiogeochemical Transformation of Redox Sensitive Contaminants in a Hyporheic Zone. Water Resources Research, 54(9), 5973-5985. Yu, F., & Hunt, A. G. (2017). Damköhler number input to transport-limited chemical weathering and soil production calculations. Earth Space Chemistry, 1, 30-38.

Figures

Fig. 1. Periodic velocity fluctuation scenarios considered in the simulations

(a)

(b) Fig. 2. Magnesite dissolution rate and chemical concentration profile under steady state condition. (a) Effluent Mg2+ concentration and average dissolution rate in the column as a function of flow velocity under the steady-state condition. Solid squares denote the average reaction rate in the column and open circles denote Mg2+ concentration in effluent solution. As a comparison, experimental results of Mg2+ concentration from Salehikhoo et al. (2013) were also plotted (red open diamond). (b) Mg2+ concentration and pH distributions in the column under the steady-state condition. Symbol color indicates different flow velocity: 0.36 m/d (black), 3.6 m/d (green), and 36 m/d (red) flow velocities.

(a)

(b)

(c)

(d) Fig. 3. (a) Permeability distribution (log K) of two heterogeneous domains. (b) pH distribution in two heterogeneous domains under different flow velocities (uave=1.2 m/d, 3.6m/d). (c) local magnesite dissolution rate distribution in two heterogeneous domains under different flow velocities (uave=1.2 m/d, 3.6m/d). (d) Simulation results of macroscopic magnesite dissolution rate in 10 different heterogeneous domains under 4 different average flow velocities.

(a)

(b)

(c) Fig. 4. (a) Fluid travel time distributions f(τ) calculated from the effluent tracer concentration in the flow-through systems (see text). (b) Change of Mg2+ concentration and average dissolution rate in a batch reactor as a function of water-mineral contact time (the solid/solution ratio in the batch reactor was the same as in the column reactor). (c) Average reaction rates calculated from reactive transport simulations ( rMgn ) (line) and from the FTTD based model ( rFT,Mgn ) (symbols).

(a)

u=1.2 m/d

u=12 m/ (b)

(c)

Fig. 5 (a) Permeability distributions of five heterogeneous flow through domains with 2.2 m length and 0.44 m width. (b) Simulation results of FTTDs for different subdomains and flow velocities. The tracer break-through curves were observed for the subdomains with lengths 0.04, 0.1, 0.22, 0.44, 0.9, 1.5, 2.2 m (from the inlet) after a δ pulse tracer. (The FTTDs under hetro_1 case were presented as example. The rest FTTDs can be found in Fig. S3.) (d) Log-log plot for the average reaction rates calculated from reactive transport simulations ( rMgn ) approach ( rFT,Mgn ).

and those from the FTTD based

(a)

(b) Fig. 6. (a) FTTD based scaling result ( rFT,Mgn ) under different sampling frequencies (NS) (b) Maximum deviation between the average reaction rates calculated from reactive transport simulations ( rMgn ) and those from the FTTD based approach ( rFT,Mgn ) under different NS.

(a)

(b) Fig. 7. The effluent Mg2+ concentrations as a result of magnesite dissolution in the column with variable flow velocity. (a) umin =0 m/d, umax =18.0 m/d, Rd=0.2, (b) umin=1.8 m/d, umax =5.4 m/d, Rd=0.5. The time (x-axis) was normalized to total time duration of each fluctuation wave.

(a)

(b)

Fig. 8 (a) Change of average dissolution rate rMgn with the wave length T. (b) Change of rMgn with the velocity difference between the maximum and minimum velocity values Δu.

(a)

(b)

(c)

(d) Fig. 9 (a) Accumulative effluent volume QE, (b) the fluid travel time  of the fluid entered column system at time t, (c) the travel time distribution function (umin=0.9 m/d, umax=27.9 m/d, Rd=0.1), (d) the accumulative FTTD F ( ) .

f ( )

Fig. 10. The comparison between the average reaction rates obtained from microscopic reactive transport simulations ( rMgn ) and those from the FTTD-based approach ( rFT,Mgn ) under transient flow conditions.

Fig. 11. The FTTD-based reaction rates as a function of average Da   adv  react (dimensionless fluid travel time vs characteristic reaction time  react ). The FTTD functions for the cases denoted by black circle and red diamond symbols were obtained from the tracer pulse simulation in random heterogeneous domains (with length ranged from 0.04 to 2.2 m) with 1.2 m/d and 12 m/d flow velocities. The purple triangle and blue square cases were the results redrawn from Fig. 4 and Fig. 8. The black solid-line is the transport-limited rate that be calculated as rlim  csat uAcross AMgn  Da 1 , where csat is the equilibrium Mg2+ concentration. The green dash line is the empirical

correlation between average reaction rate and Da number under mixed-control condition

(0.01
r  5.0  10 11  Da 0.645 .

and

the

specific

correlation

equation

is

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: