Relating Darcy-Scale Chemical Reaction Order to Pore-Scale Spatial Heterogeneity.

Due to spatial scaling effects, there is a discrepancy in mineral dissolution rates measured at different spatial scales. Many reasons for this spatial scaling effect can be given. We investigate one such reason, i.e., how pore-scale spatial heterogeneity in porous media affects overall mineral dissolution rates. Using the bundle-of-tubes model as an analogy for porous media, we show that the Darcy-scale reaction order increases as the statistical similarity between the pore sizes and the effective-surface-area ratio of the porous sample decreases. The analytical results quantify mineral spatial heterogeneity using the Darcy-scale reaction order and give a mechanistic explanation to the usage of reaction order in Darcy-scale modeling. The relation is used as a constitutive relation of reactive transport at the Darcy scale. We test the constitutive relation by simulating flow-through experiments. The proposed constitutive relation is able to model the solute breakthrough curve of the simulations. Our results imply that we can infer mineral spatial heterogeneity of a porous media using measured solute concentration over time in a flow-through dissolution experiment.

Introduction

Geochemical reactions such as mineral dissolution play an essential role in determining water chemistry, soil formation, biogeochemical cycling, and global climate (Wen and Li 2017). Mineral reactions can also occur by injecting CO in subsurface reservoirs during geothermal energy extraction (Randolph and Saar 2011; Adams et al. 2021; Ezekiel et al. 2022) or calcite mineralization (Pogge von Strandmann et al. 1983). One of the most significant obstacles to understanding the geochemical reactivity of natural subsurface environments stems from the multitude of spatial scales that have to be considered (Noiriel et al. 2012). Due to spatial scaling effects, mineral dissolution rates are known to be 3–6 orders of magnitude lower in the field than when measured in the laboratory (White and Brantley 2003; Navarre-Sitchler and Brantley 2007; Maher 2010; Moore et al. 2012). The deviation in mineral dissolution rates strongly limits the extrapolation of kinetic dissolution models and parameters characterized at the laboratory to natural systems (Li et al. 2020).

Many factors are responsible for the spatial scaling effects of mineral dissolution rates. This work focuses on how pore-size distribution—which contributes to hydraulic heterogeneity—and spatial mineral distribution causes such spatial scaling effects. Li et al. (2007) performed simulations with various settings of spatial mineral distribution and concluded that spatial mineral distribution has a significant scaling effect when the reactive minerals are of small but typical proportions. Experimental studies using a column packed with quartz and magnesite have confirmed the significant role of spatial heterogeneities in subsurface reactive transport and can be used to quantify the effect of spatial mineral distribution on dissolution rates (Salehikhoo et al. 2013; Li et al. 2014; Li and Salehikhoo 2015). X-ray micro-tomography provides observations of the impact of physical and chemical heterogeneity on reaction rates in multimineral porous media (Tutolo et al. 2015; Luhmann et al. 2017; Al-Khulaifi et al. 2017, 2018, 2019; Menke et al. 2016, 2018). Fischer et al. (2014) and Fischer and Luttge (2017) studied how mineral surface roughness at the nanometer scale affects surface reaction rates and proposed to upscale the mineral reaction rate using Monte Carlo simulations. Ma et al. (2019, 2021) quantified the accessible surface area of minerals in a sandstone using scanning electron microscopy (SEM) images and Brunauer–Emmett–Teller (BET) surface area measurements.

Regarding the influence of hydraulic heterogeneity, Wen and Li (2017) and Jung and Navarre-Sitchler (2018a) performed reactive transport simulations on stochastically generated permeability fields and studied how hydrologic heterogeneity affects mineral dissolution rates. Using Monte Carlo simulations, Jung and Navarre-Sitchler (2018b) further studied the time dependency of mineral dissolution rates, and Wen and Li (2018) developed an upscaled rate law for mineral dissolution in heterogeneous media under variable residence-time and length-scale conditions. Li et al. (2020) upscaled mineral dissolution rates in a porous medium with a random permeability field using the fluid travel-time distribution function. The works mentioned above used the transition state theory (Lasaga 1998) to model mineral dissolution rates, with a macroscale reaction order of unity . This is reasonable since there are no physical explanations why this macroscale reaction order should not be one (Lasaga 1998; Brantley and Conrad 2008).

However, there are rate models with a macroscale reaction order of 2 in kinetics of crystal growth (Nancollas 1968; Reddy 1975, 1977). Such second-order kinetics are used for modeling surface spiral growth (Nielsen 1984). Considering calcite as our mineral of particular interest, fitting experimental data using a reaction order larger than 1 is common, especially when the saturation is close to equilibrium (Plummer and Wigley 1976; Plummer et al. 1978; Palmer 1991; Svensson and Dreybrodt 1992). A higher reaction order is also observed in modeling calcite dissolution in seawater (e.g., Subhas et al. 2015; Naviaux et al. 2019).

In this work, we use analytical techniques to develop a constitutive relation of mineral dissolution kinetics in porous media based on models with a reaction order . We characterize hydraulic heterogeneity and mineral spatial heterogeneity by the longitudinal dispersivity and the reaction order. We use the bundle-of-tubes analogy to show how the reaction order relates to both hydraulic and spatial mineral heterogeneity in porous media. Furthermore, we simulate experimental scenarios involving advective and dispersive transport using such a constitutive relation. Our results show how concentration breakthrough curves of the reactive species reveal information of both the hydraulic and chemical heterogeneity of porous media.

Materials and Methods

This section first introduces reactive transport models at the pore scale and the Darcy scale. Then, we lay out statistical distributions of pore sizes and effective-surface-area ratios. Such distributions can thus define the volume-averaged concentration. Constitutive relations based on the Darcy-scale reaction order is established using Taylor series expansions. We check the applicability of the constitutive relations by comparing the modeled concentration and the volume-averaged concentration using a goodness-of-fit measure, the Jensen–Shannon divergence. Finally, we explain how one can apply the proposed constitutive relation using a flow-through experiment.

Reactive Transport at the Pore Scale

At the pore scale, we model the transport of a mineral-forming solute bywhere is the solute concentration in the fluid (mol m), is the fluid velocity (m s), D is the molecular diffusivity (m s). In Eq. (1), we consider as the unknown variable, and the fluid velocity, , is given. The heterogeneous reaction is defined aswhere is the unit normal vector pointing outwards from the fluid to the solid surface, and is the rate of the heterogeneous reaction per surface area (mol m s) as a function of . Since we consider a bundle of capillary tubes as a model of a porous medium, we introduce Eq. (1) described in cylindrical coordinates:and the boundary condition, Eq. (2),where r is the radius of the capillary tube (m). We assume fully developed axisymmetric fluid flow, , constant molecular diffusivity, and axisymmetric solute concentration. Furthermore, we assume the heterogeneous reaction causes negligible change to the tube radius and the mineral surface area. Regarding dissolution and precipitation reactions leading to changes in the pore geometry of porous media and fractures, we refer the reader to the works of Sallés et al. (1993), Békri et al. (1995, 1997) for further discussion. Following Arce et al. (2005), the area-averaged concentration of the capillary tube is

To develop an area-averaged or upscaled equation, we perform area-averaging of Eq. (3) to obtain

Evaluate the integral and rearrange:

Following Paine et al. (1983), we further expand the convective term using Gray (1975)’s representation:where and are the spatial deviation terms of the fluid velocity and solute concentration, respectively. Considering only transport, Paine et al. (1983) showed that the dispersive transport term can be represented by the Taylor–Aris theory of dispersion:with the constraints and . Ananthakrishnan et al. (1965) delineated dispersion effects for and introduced a time-dependent dispersion coefficient (m s) to model dispersion in the regime not covered by the Taylor–Aris dispersion. We utilize the time-dependent dispersion coefficient obtained by the method of moments (Barton 1983) and asymptotic techniques (Vrentas and Vrentas 1988)where is the nth root of Bessel function of the first kind of order 1 (Meng and Yang 2017). The time-dependent term in serves as a correction term for Taylor–Aris dispersion. Such a definition of dispersive transport is suitable for simple initial conditions, e.g., injecting a pulse of solute. Considering general initial conditions under steady flow conditions, one has to include source terms to address the complexity of dispersive transport in a cylindrical tube (Taghizadeh et al. 2020). Since the averaged convection term can be described by averaged quantities, we divert our focus to the source term on the right-hand side of Eq. (7):

The heterogeneous reaction term depends on the solute concentration at the solid–fluid boundary in the radial direction. For reaction–diffusion systems in a cylindrical tube with a linear irreversible heterogeneous reaction, , Arce et al. (2005) performed an order-of-magnitude estimate and claimed that when , the solute concentration on the boundary can be approximated bywhere k is a rate constant (m s). The constraint gives us a good range of k, since the pore radius is usually some tens of micrometers. Combining Eqs. (7)–(12), the area-averaged solute transport equation in a cylindrical pore reads

This is an ad hoc approach of developing the averaged model. Though not perfect, such a one-dimensional (1D) expression does represent reactive transport in a pore throat in many pore-network models (Algive et al. 2010; Raoof et al. 2012, 2013; Varloteaux et al. 2013a, b; Qin and Hassanizadeh 2015; Bekri et al. 2015; Gostick et al. 2016; Xiong et al. 2016; Esteves et al. 2020). Rigorous upscaling of reactive flow in thin geometries has been performed with the following considerations: general mass action kinetics (van Duijn and Pop 2004), dominant Péclet and Damköhler numbers (Mikelić et al. 2006; van Duijn et al. 2008), changes in pore-scale geometry (van Noorden 2009b; Kumar et al. 2011), changes in pore-scale geometry with non-isothermal effects (Bringedal et al. 2015, 2016), changes in pore-scale geometry with two-phase flow (von Wolff and Pop 2021), and coexisting homogeneous reactions (Boso and Battiato 2013). For perforated porous media, rigorous upscaling of reactive flow involving dissolution or precipitation processes has been performed (Kumar et al. 2016), with considerations of changing pore-scale geometry (van Noorden 2009a) as well as multiphase reaction–diffusion systems (Redeker et al. 2016).

The dimensionless number, kr/D, is often identified as the pore-scale Thiele modulus squared or the Damköhler number Da. Balakotaiah et al. (1995) used the invariant manifold expansion for advection and heterogeneous reaction in a cylindrical tube and showed that when , dispersion effects can be modeled by the Taylor–Aris theory. For reaction–diffusion systems in porous media, Valdés-Parada et al. (2017) showed that the constraint can be loosened to by modifying the effective diffusivity and the effective reaction rate constant at the macro scale, and Bourbatache et al. (2020) recovered classical homogenized diffusive equations for small values of Damköhler numbers (defined by the length scale of the representative elementary volume).

Second-Order Dissolution Kinetics

Throughout this work, the heterogeneous reaction of interest is mineral dissolution. A common mineral dissolution model is of second order:where is the dissolution rate constant (mol m s) and is the precipitation rate constant of second-order kinetics (mol m s). We assume the solution is dilute, such that the activity coefficient of the solute is unity. We relate the dissolution model with a more prevalent formulation involving the solubility product constant, , and the ion activity product, IAP:

One can switch between the two forms by stating and . Using calcite as an example, when pH, the dominant dissolving species at the diffusion boundary layer are Ca and CO (Sjöberg and Rickard 1984). Thus, at intermediate pH values, . The assumption of electroneutrality close to the mineral surface yields , which results in second-order dissolution kinetics (Ebigbo et al. 2012; Levenson and Emmanuel 2013). For magnesite dissolution at neutral to alkaline pH regimes, such second-order kinetics is also suitable (Salehikhoo et al. 2013; Wen and Li 2017). Second-order dissolution/precipitation kinetics also appears in upscaling reactive transport processes in porous media with attention to moving solid–fluid interfaces (Ray et al. 2019; Bringedal et al. 2020; Gärttner et al. 2020). Combining Eqs. (13) and (14) yields

We relate the dimensional and nondimensional quantities bywhere t is the nondimensional time, z is the nondimensional space, C is the nondimensional solute concentration, and is the length of the cylindrical pore (m). Then, we nondimensionalize Eq. (16) by the following scaling of time and concentrationwhere [t] is the diffusive time scale (s). The variables in square brackets remove the physical dimension of the starred variables and refer to characteristic quantities. Thus we have a nondimensional equation of reactive transportwhere is the aspect ratio, Pe is the Péclet number, and Da is the Damköhler number. In Eq. (20), we further relate Da to and the pore-scale Thiele modulus (of second-order reactions), denoted as , to show that for such a definition of the Damköhler number, does not necessarily break the assumption of a small Thiele modulus. For a slender cylindrical pore, , and a controlled flow rate, , we neglect the effect of dispersion. The term in the dispersion correction term also indicates such a correction term vanishes rapidly over time. Hence, the nondimensional reactive transport equation reads

First-Order Dissolution Kinetics

First-order-kinetic models are uncommon, since chemical reactions often involve two reagents (Cussler 2009). If one wanted to make use of first-order kinetics, one would have to assume the concentration of a mineral-forming ion is in excess or constant (Meile and Tuncay 2006) or limit the usage of first-order kinetics to low solute concentrations (Kaufmann and Dreybrodt 2007). Nonetheless, we introduce the model of first-order mineral dissolutionwhere is the precipitation-rate constant of first-order kinetics (m s). First-order dissolution kinetics is utilized to model evolution of karst aquifers (Gabrovšek and Dreybrodt 2001). They apply to dissolution rates of various minerals such as gypsum, rocksalt, calcium carbonate, and quartz (Jeschke and Dreybrodt 2002). Considering a general surface reaction, first-order kinetic models resemble adsorption and desorption kinetics studied by, e.g., Zhang et al. (2017), as well as the Noyes–Whitney model of drug dissolution (Dokoumetzidis and Macheras 2006). First-order kinetics is also applied to heterogeneous reactions between living cells and extracellular fluids (Santos-Sánchez et al. 2016). Combining Eqs. (13) and (22) and ignoring dispersion effects, we scale time and concentration bysuch that the Péclet and the Damköhler numbers are

Notice the pore-scale Thiele modulus squared, , appears in the form of the precipitation-rate constant, . Since is proportional to , the assumption of is fulfilled, similar to the relation between the Damköhler number and the pore-scale Thiele modulus of second-order kinetics. The nondimensional reactive transport equation of first-order kinetics is therefore

As for other possible models of heterogeneous kinetics, Qiu et al. (2017) provided a comprehensive review of upscaling reactive transport processes in porous media, e.g., Michaelis–Menten-type kinetics (Wood et al. 2007; Dadvar and Sahimi 2007), Monod-type kinetics (Heße et al. 2009), nonlinear kinetics that reduce to first-order kinetics when the reaction order is 1 (Guo et al. 2015).

We define the average velocity in a cylindrical pore using the Hagen–Poiseuille equationwhere is the dynamic viscosity of the fluid (Pa s) and is the pressure difference between the inlet and the outlet (Pa).

Reactive Transport at the Darcy Scale

We use the bundle-of-tubes analogy to model reactive transport at the Darcy scale (Kozeny 1927). We define the specific mineral surface area (m) of the porous medium aswhere A is the mineral surface area, V is the bulk volume of the porous media (m), is the ratio between mineral surface area and the total surface area of a pore, and N is the total number of pores. We assume that the pores have the same length, , as the porous medium, such that the tortuosity is 1. The porosity of the porous medium iswhere is the fluid volume. Dividing Eq. (27) by Eq. (28), we obtain

Recall the Damköhler number for first-order kinetics, Eq. (24), the Damköhler number at the Darcy scale is thereforewhere the characteristic length L is chosen to be the length of the porous medium, . The Darcy-scale Damköhler number for second-order kinetics is

Since we consider the porous medium as a bundle of tubes, the seepage velocity of the porous medium can be defined using a volume-averaged velocity,

Hence, the Péclet number at the Darcy scale is

Dispersion effects arise when the pore sizes are not uniform (Carbonell 1979; Arriaza and Ghezzehei 2013; Meng and Yang 2017). Therefore, we introduce a longitudinal dispersion coefficient (m s),where is the longitudinal dispersivity with the unit of length (m). We derive the longitudinal dispersion coefficient using the spatial moments of the averaged solute concentration of the bundle of tubes model

The shape factor f(r) abbreviates the expression of pore sizes in Eq. (35). Detailed derivations are presented in “Appendix 1”. The longitudinal dispersivity iswhere is the center of mass of an instantaneous source injected in the porous domain. The spatial moment analysis showed that the longitudinal dispersivity increases linearly as the distance traveled of an instantaneous source, which is reasonable in an infinite domain. In recent studies based on pore-network modeling, Mahmoodlu et al. (2020) observed the longitudinal dispersivity increases as travel distance increases.

By scaling the longitudinal dispersion coefficient by the molecular diffusivity yields the nondimensional reactive transport model at the Darcy scale,where is the reaction term that is yet to be defined. Such a Darcy-scale model resembles the upscaled formulation of heterogeneous reaction systems including convective transport in porous media, which incorporates the effect of pore geometry in the dispersivity tensor and the effective reaction-rate constant (Wood et al. 2007; Battiato et al. 2009; Battiato and Tartakovsky 2011; Valdés-Parada et al. 2011; Guo et al. 2015; Qiu et al. 2017). Our proposed model assumes the reaction term does not contribute to the overall dispersion coefficient. The flow-through experiment introduced in Sect. 2.7 circumvents such effects by eliminating spatial variation of solute concentration. Note that when all pores have the same radius and effective-surface-area ratio, Eq. (37) reduces to the equation for a single pore. In the next section, we discuss statistical distributions of pore sizes and effective-surface-area ratios.

Statistical Distributions of Pore Sizes and Effective-Surface-Area Ratios

The pore sizes of porous media usually follow a log-normal distribution (Shi et al. 1991; Hefny et al. 2020). The probability density function of log-normally distributed pore sizes, r, iswhere and are the mean and standard deviation of the variable’s natural logarithm, respectively. The sum of pore sizes raised to the power of m can be described by a raw statistical moment,

If the pore-size distribution is known, the intrinsic permeability, , of the bundle can be calculated as

Using the moments of the lognormal distribution, we can now clarify the shape factor,

Since we are interested in not only the pore-size distribution but also the effective surface area of the reactive minerals, we utilize the effective-surface-area ratio, , which has a value between zero and one. A convenient choice for modeling the distribution of a variable bounded by zero and one is the beta distribution. The probability density function of the beta distribution iswhere is the gamma function. Variables a and b shape the beta distribution.

The specific surface area over porosity, , which is a part of the Darcy-scale Damköhler number, Eqs. (30) and (31), can be defined as

Such a definition of depends on the pore geometry. Inferring from Hussaini and Dvorkin (2021)’s compilation of specific surface area versus porosity using digital images of natural rocks, can range from 0.15 (Fountainebleau sandstone) to 0.7 (Kentucky sandstone) when E.

A Constitutive Relation for Reaction at the Darcy Scale

A constitutive relation is an additional equation that specifies properties of a material. For example, the longitudinal dispersivity describes the dispersion of solute introduced by variations in fluid velocities in a porous material. In this section, we develop a constitutive relation that models the average solute concentration of a porous medium with varying reaction rates in pores.

First-Order Kinetics

Before we approach the full reactive transport problem—Eq. (37)—we start with a single pore involving only mineral dissolution,where is the initial concentration of the mineral-forming solute in the pore. We consider an aspect ratio, , small enough such that a volume-averaged concentration is representative. Using Eq. (24), the Damköhler number in a cylindrical pore iswhere defines the pore-scale geometry and defines the physics of the problem. In this work, we focus on analyzing how pore-scale geometry affects solute concentration and reaction rate over time, assuming constant . We abbreviate to . The solution for the single-pore reaction problem is

We define the volume-averaged solution as

The bounds of the integral over pore size, r, should correspond to the bounds of the prescribed pore-size distribution. Consider and expand the volume-averaged solution using a Taylor series around :

Let us now describe the Darcy-scale reaction in the same manner as we did for a single pore, Eq. (44). In this case, the definition of the Darcy-scale Damköhler number,leads to a first-order approximation of the volume-averaged concentration. To better approximate , we propose a nonlinear reaction-rate model as a constitutive relation,where n is the reaction order (Lasaga 1998). Many researchers have attempted to explain values of the reaction order in terms of dissolution or precipitation processes (Blum and Lasaga 1987; Teng et al. 2000). However, attributions of a process on the basis of this exponent are generally not defensible without further observations (Brantley 2008). For further discussion, see Brantley (2003, 2008).

The solution to the nonlinear reaction-rate model, Eq. (54), isIts Taylor series expansion around while isWe observe that the Darcy-scale Damköhler number, Eq. (53), still matches the first-order term of the volume-averaged solution, Eq. (52). If we definethen the nonlinear reaction-rate model approximates the volume-averaged concentration to the second order with respect to time. By Cauchy–Schwarz inequality, we infer , which agrees with experimental observations. The inverse square root of this definition of the reaction order, n, is also known as Tucker’s congruence coefficient, which assesses similarity between two variables (Lorenzo-Seva and ten Berge 2006).

Figure 1 shows scatter plots of pore sizes and effective-surface-area ratios. Each point represents an observation of the pore size and the effective-surface-area ratio in a porous sample. Tucker’s congruence coefficient, , measures the similarity between pore size and effective-surface-area ratio. From the leftmost figure to the rightmost figure, the congruence coefficient decreases as the observations become less similar, or more heterogeneous. Since the reaction order, n, is the squared inverse of the congruence coefficient, the reaction order increases as the heterogeneity increases. Such a definition of the reaction order is a function of the geometric variables and r. Thus we can use the reaction order to infer pore-scale spatial heterogeneity of minerals.

Fig. 1

This figure shows scatter plots of pore sizes and effective-surface-area ratio. The title of each plot shows the reaction order, n, and Tucker’s congruence coefficient,

Second-Order Kinetics

We model the single-pore problem with second-order kinetics bywhere the solution is

Using Eq. (20), the Damköhler number in a cylindrical pore iswhere , which we abbreviate to , defines the physics of the problem for second-order kinetics. The volume-averaged concentration is defined using Eq. (47), and we apply Taylor series expansion to the volume-averaged concentration,where is the Bernoulli number (Oldham et al. 2009). The Taylor series expansion of the hyperbolic tangent function, Eq. (60), converges for , which is not of concern since we utilize only the derivatives of .

We propose the following constitutive relation that describes the solute concentration at the Darcy scale,We do not attempt a solution for . However, we can still expand around with Taylor series:

See “Appendix 2” for detailed derivations. Comparing the third-order term of Eq. (64) with that of Eq. (66) yieldsOne can utilize the quadratic formula to explicitly determine n,where we consider only the larger value of n as a solution. Note that the reaction orders for first- and second-order kinetics are non-dimensional and, most importantly, independent of the length scale and the reaction rate constants.

Goodness of Fit Between the Pore-Scale and the Darcy-Scale Concentrations

There exists a variety of goodness-of-fit measures between models and experimental observations. For example, the coefficient of determination, , is often used to determine the kinetic rate law when applying the integral method (Brantley and Conrad 2008; Zhao and Skelton 2014). The mean-squared error is also a goodness-of-fit measure, and the least-squares approach tends to minimize such a metric. We use another goodness-of-fit measure, the Jensen–Shannon divergence, which is based on the Kullback–Leibler divergence (Kullback and Leibler 1951). The Kullback–Leibler divergence between some unknown distribution, p(x), and an approximating distribution, q(x), is:

The Kullback–Leibler divergence satisfies with equality if, and only if, (Bishop 2006). Although , it may diverge to infinity depending on the underlying densities (Nielsen 2020). Thus we use the Jensen–Shannon divergence,which is bounded between 0 and 1 when using base-2 logarithms (Lin 1991). Throughout this work, we use the Jensen–Shannon distance, which is defined as the square root of the Jensen–Shannon divergence.

Such a metric measures the distance between probability distributions (Endres and Schindelin 2003; Österreicher and Vajda 2003; Levene and Kononovicius 2019). The following describes how we apply this measure to solute concentration over time, and . Suppose we regard solute concentration over time as cumulative distribution functions. In that case, we measure the Jensen–Shannon divergence of their derivatives, which can be seen as the probability density functions or the reaction rates over time.

When the observed solute concentration is not monotonically increasing over time like a cumulative distribution function, we simply use the root mean square error (RMSE) as a quality measure of the constitutive relation,

Determination of the Darcy-Scale Damköhler Number, , and the Reaction Order, n, using power series

We test the constitutive relation using the volume-averaged concentration , which is an analogy of solute concentration measurement from a dissolution experiment. Assume the solute concentration can be described by the constitutive relation within a certain error. Then we can use the Taylor series expansions, Eqs. (57) and (66), to obtain the Darcy-scale Damköhler number,which is the initial rate of reaction. When the kinetics is of first order, the reaction order is obtained by differentiating Eq. (57) twice,

Differentiating Eq. (66) thrice and rearranging yields the reaction order for second-order kinetics,

This method of determining and n utilizes power-series expansion and requires only the derivatives of concentration at , given .

We consider three sets of log-normally distributed pore sizes, , , and , which have but different variances. The pore sizes are chosen such that they range from 10 to 80 m (Gong et al. 2020). Likewise, we assume the effective surface area follows the beta distribution, where considers most pores are fully reactive, assumes a larger variance of mineral surface area in the pores, and implies that the reactive mineral constitutes a small portion of the porous sample. Figure 2 shows the details of the aforementioned probability distributions.

Fig. 2

This figure shows the probability density functions of the distribution of pore sizes and effective-surface-area ratio. The legends state the essential parameters for generating the probability density functions. The definitions of the probability density functions are defined in Eqs. (38) and (42)

The products of the random variables R and form nine scenarios of the bundle-of-tubes model, which can be used as benchmarks for our power-series approach to obtain the Darcy-scale Damköhler number and the reaction order. We compare this method with a goodness-of-fit minimization using both and n as unknowns, similar to the ideas of nonlinear least-squares model fitting (Fogler 2016). Initially, the pores are filled with dissolving fluid with no solute concentration, . Then the mineral starts to dissolve into the fluid, such that the solute concentration increases. We assume we can observe the average concentration, , without transport effects. To capture the full reaction behavior, the simulation ends when concentration is larger than 0.99. The physics-related parameters, , is set as 50, such that the Darcy-scale Damköhler numbers of the scenarios are at a similar scale.

Flow-Through Experiment

In the previous section, we test the constitutive relation considering only mineral reaction. To measure the solute concentration of the fluid in a porous sample, one has to push the fluid out of the porous sample. We now discuss the full reactive transport problem. Consider a flow-through experiment, i.e., injecting fluid that dissolves the mineral in a porous sample. We collect the fluid from the outlet and measure the solute concentration over time. Before the experiment, the porous sample should be saturated by the dissolving fluid, which will result in an initial concentration of the solute in the fluid. Then we start the flow-through experiment by injecting fluid without solute under high Péclet number. This process aims to flush out all residual concentrations, such that the initial condition of the concentration corresponds to the upscaling approach considering only reaction. We take fluid samples from the outlet and measure the solute concentration. Since we are injecting fluid with zero solute concentration, we expect the measured solute concentration at the outlet to decrease over time. When the outlet concentration is close to zero, we reduce the Péclet number to 0.1–1% of the original Péclet number in order to observe an increase of the outlet solute concentration. The flow-through experiment creates a V-shaped curve of outlet concentration over time.

We use FEniCS (Alnaes et al. 2015) to solve the transient initial boundary value problem of each pore and apply Eq. (47) to obtain the average concentration. We discuss two types of porous samples, namely the combinations of the pore-size distribution and the effective-surface-area ratio, and . We test the porous samples for the low and high Damköhler-number cases, described in the following subsections. In all cases, we consider the molecular diffusivity, D, as  m s.

The Low Damköhler-Number Case,

The Darcy-scale Damköhler number, Eq. (53), is proportional to the length scale squared. This case is suitable when the mineral has low reaction rates or when the porous domain is short (small length scale), e.g., a 5  cm rock sample in a laboratory. We consider this as our “small length-scale scenario” for which the outlet solute concentration can be measured. Since the Damköhler number is low, the solute concentration during injection of the dissolving fluid should be far from chemical equilibrium. Therefore, we assume the reaction is of first order.

Following the procedure of the flow-through experiment, we set the initial Péclet number as 10, and the Darcy-scale Damköhler number is 0.1. We reduce the Péclet number to 0.01 at 0.15 nondimensional time. The simulation ends at 1.5 nondimensional time, which for the 5 cm rock sample mentioned above would correspond to about 8.7 days.

The solute concentration over time at the outlet of the porous sample is collected from the simulation data. Knowing the Péclet number and that the mineral reaction is of first order, we fit the observed concentration over time with the reactive transport model,

The shape factor, f(r), is defined using Eq. (35). Utilizing the optimization procedures in SciPy (Virtanen et al. 2020), we find and n by minimizing the RMSE between the observation and the model.

Another method of fitting and n is to utilize a part of , where the diffusion effects are dominant enough () such that we can treat the concentration as constant over space. Owing to the divergence theorem and the boundary condition results in

Then, we can perform the least-squares fitting on the left-hand side to determine and n. Figure 3 shows the outlet concentration of the flow-through experiment and highlights the region in which we consider diffusion to be dominant. Such an approach relies heavily on the strong-diffusion assumption and is therefore not suitable for the high Damköhler-number case, discussed in the next section.

Fig. 3

The outlet concentration over time of the low Damköhler-number problem. The orange line shows the part where we apply the least-squares fitting technique Eq. (77). The vertical dashed line indicates the time, , when the injection rate is reduced

The High Damköhler-Number Case,

Consider a field experiment in an aquifer, where the length scale is on the order of meters. Here, a fluid is injected in one well and produced at another well. In this case, it may not be possible to observe a concentration breakthrough at the outlet. Thus, we perform simulations of a flow-through experiment and measure the concentration at the inlet over time. The Darcy-scale Damköhler number is . We assume a second-order kinetics model for the mineral reaction, and our reactive transport model is

We consider three cases of initial Péclet numbers, , , and . The injection stops at dimensionless time, and the simulation ends at dimensionless time, which is roughly 231.5 days considering a 100 m simulation domain.

Though it is practically not possible to observe the inlet concentration during the injection phase, we perform the fitting of , n, and f(r) using all information of at the inlet. In the latter phase when injection stops (rising limb of the curve in Fig. 3), the concentration at the inlet can be determined by sampling the fluid in the injection well.

Results

In this section, we show the benchmarks of the power-series approach and the simulations of flow-through experiments. Then, we discuss the results in Sect. 4.

Benchmark of the Power-Series Approach

We benchmark the power-series approach that obtains the reaction order, n, and the Darcy-scale Damköhler number, , using a bundle of tubes characterized by the distributions of pore sizes and effective-surface-area ratio described in Sect. 2.6.

Figures 4 and 5 show the benchmark of the first-order kinetics and second-order kinetics, respectively. In the upper part of the figures, we plot the contour lines of log-scaled Jensen–Shannon distance between and . The red points indicate the approximation of and n using the power-series approach. In an ideal case, the red points should be in the minimum of the Jensen–Shannon distance. In the lower part of the figures, we plot the concentration over time of all scenarios.

Fig. 4

The upper part of the figure shows the contours of the log-scaled Jensen–Shannon distance between the observed concentration and the modeled concentration of first-order reactions using the constitutive relation. The red points are the Darcy-scale Damköhler number and the reaction order approximated by the power-series approach (Eqs. (73) and (74)). The legend shows the value of the log-scaled Jensen–Shannon distance using a grayscale, corresponding to the brightness of the colored contour lines. The lower part of the figure shows the concentrations over time and of all pore sizes and effective-surface-area ratios scenarios

Fig. 5

The upper part of the figure shows the contours of the log-scaled Jensen–Shannon distance between the observed concentration and the modeled concentration of second-order kinetics using the constitutive relation. The red points are the Darcy-scale Damköhler number and the reaction order approximated by the power series approach [Eqs. (73) and (75)]. The legend shows the value of the log-scaled Jensen–Shannon distance using a grayscale, corresponding to the brightness of the colored contour lines. The lower part of the figure shows the concentrations over time and of all pore sizes and effective-surface-area ratios scenarios

Figure 6 shows the observed concentration, , and the modeled concentration, . Table 1 shows the Darcy-scale Damköhler number, the reaction order, and shape factor obtained by direct calculation [Eqs. (53), (58), and (68)], and by RMSE minimization of the low and high Damköhler-number cases ( and , respectively) corresponding to the and the distributions. In the low Damköhler-number scenarios, we use the least-squares curve fitting method, Eq. (77), to obtain and n. The results of the least-squares curve fitting are , for the scenario, and , for the scenario.

Fig. 6

The figures show the concentration (solid lines) and the modeled concentration (dashed lines) of the low Damköhler-number cases on the top panel. The high Damköhler-number cases are summarized in the middle and the bottom panel. The left and right panels show the and the cases, respectively

Table 1

This table summarizes the Darcy-scale Damköhler number, , the reaction order, n, and the shape factor, f(r), obtained by fitting the concentration-over-time curve using the constitutive relation with the RMSE metric

R3Ω3
Initial Pe_d	10	4000	8000	40,000
Direct calculation
Da_d	0.1		40,000
n	1.40		1.77
f(r)		blank	4.081×10^-2
Minimum RMSE
Da_d	0.1	52,528	45,154	38,464
n	1.57	1.57	1.66	1.30
Minimum RMSE with shape factor fitting
Da_d	0.1	50,887	42,049	38,950
n	1.68	1.55	1.59	1.53
f(r)	2.39×10^-2	1.402×10^-1	1.222×10^-1	3.168×10^-2
R_3Ω_3
Direct calculation
Da_d	0.1		40,000
n	1.74		2.60
f(r)			1.718
Minimum RMSE
Da_d	0.13	25,376	24,691	10^-6*
n	4.60	1.0	1.43	4.89*
Minimum RMSE with shape-factor fitting
Da_d	0.12	35,791	25,025	29,111
n	3.87	1.56	1.0	1.0
f(r)	16.8	2.091	1.433	0.463

The values of the direct calculation are the result of prescribing , the distribution, and the distribution. In the low case, we assumed first-order kinetics. Hence we use Eq. (58) to calculate the reaction order. We assumed second-order kinetics for the high case, and Eq. (68) is used for calculating the reaction order for second-order kinetics. The direct calculation of f(r) uses Eq. (41). The results of the least-squares curve fitting are , for the scenario, and , for the scenario

*The fitting in this case does not yield reasonable results

Discussion

In Fig. 4, we observed that the power-series approach obtains and n close to the minimum Jensen–Shannon distance. The connections between nonlinear fitting of the parameters, and n, and the geometric information of the porous medium are established since the power-series approach is exact for retrieving expected values of the pore-size and the effective-surface-area distributions.

In Fig. 5, we observed general agreement of and n obtained by the power-series approach to those at the minimum Jensen–Shannon distance. As the variance of the pore-size increases, the obtained and n (red points) stray from the minimum Jensen–Shannon distance. In the concentration plots of the , , and cases, the modeled concentration fits well when . Some discrepancy between and is present when . Comparing to the results of first-order kinetics, where fits well throughout all concentrations, our averaged model of second-order kinetics can only fit reactions far from equilibrium (e.g., ).

The power-series approach of obtaining and n suffers from the fact that: Therefore, we require other nonlinear fitting methods by minimizing the divergence between models and observations. In the next section, we discuss the intricacies of nonlinear fitting of solute concentration during flow-through experiments.

reaction rates at zero concentration can be hard to obtain, and

numerical differentiation of higher-order derivatives can yield spurious results.

The Scenario

The top-left panel of Fig. 6 shows the low Damköhler-number, first-order-kinetics case, where both the RMSE minimization and the least-squares fitting method fit the outlet concentration . Both methods of obtaining and n are accurate within 15% relative error as confirmed by direct calculation from pore-size and effective-surface-area distributions. The good agreement can be attributed to the fact that the solute can be mostly flushed out from the porous domain, due to the low Damköhler number. Hence, the solute concentration in each tube goes down to almost zero, and spatial gradients of the solute concentration in the z direction are negligible. This creates a situation similar to the problem considering only reaction effects, which is not the case for the scenarios of high Damköhler-number, second-order kinetics.

Focusing on the high Damköhler-number scenarios, we observe general agreement of fitted parameters. For the case of initial , the error in can be attributed to the fact that the volume injected is not enough, such that the dissolution in smaller pores is not observed in . The case of initial approximated closer to the prescribed value than the previous case. The third case of initial , performed the best in retrieving . This may be caused by, similar to the low Damköhler-number case, that the initial solute concentration of all pores is pushed out from the porous domain, such that is close to zero. The purpose of comparing these three cases is to emphasize that the inversion of parameters is influenced by how we perform the injection test, namely, by the selection of the initial Péclet number.

Furthermore, we performed a fitting in all cases considering the shape factor as unknown. For the low Damköhler number case, we have a good fit of f(r) within an order of magnitude. However, the obtained value for n exhibits a bigger error. This signals the vagueness of the reaction order and the shape factor in minimizing RMSE, where increasing or decreasing either one of the variables leads to similar RMSE. For the high Damköhler number cases, the obtained improved slightly comparing to the cases with a prescribed shape factor. All fitted shape factors are within an order of magnitude compared to the theoretical calculations. In particular, for the case of initial , we observe a good fit of the retrieved shape factor.

This scenario tests how effective the pore structure can be inferred in a more dispersive setup. Table 1 summarizes the obtained and n for both the low and high Damköhler-number cases. For the low Damköhler-number case, the minimum RMSE methods do not yield decent approximations of , n, and f(r). The top right panel of Fig. 6 shows that all methods failed to fit the V-shaped region (), where dispersive effects matter the most. The pore-size distribution mainly consists of smaller pores, as compared to , which means that the solute to be harder to flush out. This dispersion effect is also demonstrated in the right panel of Fig. 8, where the dispersion model fails to capture the long-tailed concentration profile. However, the least-squares fitting method yields exact , and reasonable n within 10% relative error. The least-squares fitting method performs well for the low Damköhler number in both the and cases.

Fig. 8

Comparison of the volume-averaged solute concentration and the modeled concentration. The left panel shows the case of pore size distribution , and the right panel shows the case of pore size distribution

For the high Damköhler-number scenarios, fitting the averaged concentration using the model, Eq. (78), does not recover pore-scale information in general. In Fig. 6, the concentration profile of the initial with a prescribed shape factor does not fit the averaged concentration. Compared to the results with shape-factor fitting, the prescribed shape factor yields too much dispersion for the given Péclet number. Recall that the longitudinal dispersion coefficient is defined as .

The case of initial with shape-factor fitting gives us the best result in retrieving the model parameters. As shown in both Figs. 6 and 7, the initial Péclet number is not enough to flush out the solute in all pores. However, the dispersion effect is not as pronounced as the initial case, due to the lower initial Péclet number. The key to better knowledge of the model parameters and pore-scale information is to control the initial Péclet number for less dispersion. We demonstrate that the rate of injection influences the parameter fitting of the averaged model, due to dispersion caused by the variance of pore sizes. How to properly choose the initial injection rate for inference of pore-scale information could be the subject of future studies.

Fig. 7

This figure shows the concentration of the pores, C, the averaged concentration, , and the modeled concentration, , of the selected scenarios with shape factor fitting. The left panel shows the solute concentration in the porous domain during fluid injection. The middle panel shows the solute concentration after fluid injection has stopped. The right panel shows the last time step of the simulation, as indicated by the increasing inlet concentration

Figure 7 shows the concentration in pores, the averaged concentration, and the modeled concentration of the high Damköhler number scenarios. When the injection stops, the increase of the inlet solute concentration is not only due to the reaction, but also due to the diffusion of solute from the reservoir to the inlet. We attribute the underestimated n to our reactive transport model, Eq. (78), not being able to capture the averaged diffusion effects in each pore, which results in a lower reaction order. Certainly, this situation is not as ideal as the low Damköhler-number situation, which decouples transport and reaction by flushing out almost all of the residual concentration. The modeled concentration, , still represents the average behavior of the pore concentrations.

To summarize, although imperfect, our method of parameter estimation using a solute concentration breakthrough curve is useful for modeling the average behavior of reactive transport in porous media. The results suggest it is possible to infer pore-scale information using the inversion of averaged parameters.

In all of the flow-through experiments, we considered only one observation point, either the fluid inlet or the fluid outlet. This work serves as a demonstration of the base case with only one observation. To improve the fitting of the reaction order of the high Damköhler-number cases (especially the ones with distribution), one can incorporate more observation points, spatially distributed within the domain.

The Applicability of the Constitutive Relations

By adding an exponent to the pore-scale reaction model, the Darcy-scale reaction model effectively describes the averaged behavior of reactions taking place independently in different pores of the porous medium under a specific flow-through experiment. We examine the applicability of the simple approach using the averaged concentration of first-order kinetics as an example. The averaged concentration, Eq. (48), can be considered as a continuous mixture of exponential distributions,where is a parameter that characterizes the exponential distributions. The finite mixture is known as the hyperexponential distribution, which is utilized for fitting long-tail distributions (Feldmann and Whitt 1998; Okada et al. 2020). If is a gamma distribution, then is a Pareto distribution (Balakrishna and Lai 2009). The concentration of our proposed model, Eq. (55), has the following form when ,which is the cumulative distribution function of the Pareto distribution, also called the Lomax distribution. There exists a particular ratio distribution of the effective surface area and the pore sizes, , that satisfies . Such an existence contributes to the effectiveness of the constitutive relation for the first-order kinetics.

In contrast, the theoretical basis of constitutive relations for second-order kinetics is less distinct. We simply followed the derivations of first-order kinetics and exploited the oddity of and to obtain a second-order approximation. The solution for second-order kinetics in a single pore, Eq. (60), can be recast intowhich is a cumulative distribution function of a scaled logistic distribution. If we consider time a semi-infinite domain , such a distribution is also known as the half logistic distribution (Balakrishnan 1985). Though we did not find or derive the relationships of the logistic distribution mixture, we denote the possibility of approximating such a mixture using the Pareto distribution, Eq. (80), since the hyperbolic tangent function can be represented by a Laurent serieswhich is a mixture of exponential distributions. Such expansion techniques would avoid relying on derivatives around as is done in this study, which may be advantageous but requires more research.

Limitations and Outlook

The main limitation of this study is that the reaction model we considered is simple and may not be adequate to describe complex geochemical processes, e.g., a rock sample which consists of multiple dissolving minerals. We also assumed the amount of dissolved mineral is minimal such that the geometry of the pores is not affected. Though we considered two common models of first- and second-order kinetics, we treat the rate constant, and , as constant, and indeed, and may change as temperature, pH, or ionic strength changes.

Moreover, the assumption that chemical reactions in each tube occur independently of the other tubes is idealized. A better approach is utilizing pore-network models. In pore-network modeling, the porous medium is discretized as a network of pore bodies and pore throats, and the coordination number is defined as number of connections to each pore body. Experimental studies have reported that the average coordination number of a sandstone is (Ioannidis and Chatzis 2000; Øren and Bakke 2003; Hefny et al. 2020). Our bundle-of-tubes model has an average coordination number of 0 (or 1, if one considers the pore bodies reside at the fluid inlet and outlet boundaries), which is the base case for pore-network modeling. For this base case, the reaction order does not depend on or . The dependence of the reaction order on or for larger coordination numbers requires further studies.

Nonetheless, our simple model reveals a possible mechanistic explanation to the usage of the Darcy-scale reaction order larger than one, and how it can reveal geometric information of the porous medium using the solute breakthrough curve. We suggest considerations of the aforementioned limitations as potential topics for future research.

We propose another possible application of this work in view of energy conservation in a porous sample, where there is only one definition of temperature of the fluid, and the porous solid, . A special case of energy conservation without pressure work and viscous heating is analogous to Eq. (13):where is the thermal diffusivity of the fluid, and is a heat source introduced by the porous solid. Usually, such an energy conservation model of heat tracer tests assumes thermal equilibrium between the fluid and the porous solid, , (Shook 2001; Anderson 2005; Saar 2011). However, studies and modeling on thermal disequilibrium between fluid and solid phases have gained interests lately (Karani and Huber 2017; Koch et al. 2021). If we consider a heat tracer test where we create a breakthrough curve like those in Fig. 3, we can model the behavior bywhere is a certain average of thermal diffusivity of the porous media, and is a characteristic length that defines the heat flux between the solid–fluid boundary. The formulation is similar to our study of reactive transport with first-order kinetics. Therefore, it is possible to apply the same techniques described in this work to obtain the “reaction order” and infer pore-scale information.

Conclusion

Mineral reaction kinetics defined at the pore scale are not necessarily valid at the Darcy scale. We utilize a bundle-of-tubes model to study the modeling of dissolution kinetics in porous media at the Darcy scale. By adding an exponent, n (i.e., the Darcy-scale reaction order), to first- and second-order kinetics, the resulting constitutive relation approximates the average dissolution rate of the bundle-of-tubes model under specific experimental conditions. Using the pore-size and the effective-surface-area ratio distributions to characterize the tube bundles, we expand the solute concentration of dissolving species with Taylor series and thus relate the Darcy-scale Damköhler number, , and reaction order with the distribution moments. The Taylor series expansions show that the Darcy-scale reaction order of first-order kinetics is the inverse square root of Tucker’s congruence coefficient (also known as the cosine similarity) between the pore sizes and effective-surface-area ratios. Therefore, an increase of reaction order indicates an increase of pore-scale heterogeneity. Such a relation gives a mechanistic meaning to the reaction order.

Furthermore, we simulate flow-through experiments of dissolving porous media at the laboratory as well as the field scale and discuss how one can utilize the constitutive relation by fitting a solute concentration breakthrough curve with and n as unknowns. As an additional benefit, we discuss cases of the flow-through experiments where the shape factor of longitudinal dispersivity is also considered as a fitting parameter. The inversion is successful, and the fitted parameters are close to the prescribed parameters calculated by the moments of pore-size and effective-surface-area ratio distributions. We infer that: The relations we derived provide us a quantitative approach to measure the spatial heterogeneity of a porous domain using the Darcy-scale reaction order and reveal a mechanistic explanation for .

detailed pore-scale information (characterized by functions of moments) can be inferred using averaged Darcy-scale quantities (such as solute concentration), and

by analyzing the solute concentration of dissolving minerals over time using flow-through experiments, we can acquire the Darcy-scale reaction order and the dispersion coefficient, which represent heterogeneity at the pore scale.