Articles | Volume 25, issue 11
Hydrol. Earth Syst. Sci., 25, 5905–5915, 2021
Hydrol. Earth Syst. Sci., 25, 5905–5915, 2021

Research article 11 Nov 2021

Research article | 11 Nov 2021

Feedback mechanisms between precipitation and dissolution reactions across randomly heterogeneous conductivity fields

Feedback mechanisms between precipitation and dissolution reactions across randomly heterogeneous conductivity fields
Yaniv Edery1, Martin Stolar1, Giovanni Porta2, and Alberto Guadagnini2 Yaniv Edery et al.
  • 1Faculty of Civil and Environmental Engineering, Technion, Haifa, Israel
  • 2Department of Civil and Environmental Engineering, Politecnico di Milano, 20133, Milan, Italy

Correspondence: Yaniv Edery (


Our study investigates interplays between dissolution, precipitation, and transport processes taking place across randomly heterogeneous conductivity domains and the ensuing spatial distribution of preferential pathways. We do so by relying on a collection of computational analyses of reactive transport performed in two-dimensional systems where the (natural) logarithm of conductivity is characterized by various degrees of spatial heterogeneity. Our results document that precipitation and dissolution jointly take place in the system, with the latter mainly occurring along preferential flow paths associated with the conductivity field and the former being observed at locations close to and clearly separated from these. High conductivity values associated with the preferential flow paths tend to further increase in time, giving rise to a self-sustained feedback between transport and reaction processes. The clear separation between regions where dissolution or precipitation takes place is imprinted onto the sample distributions of conductivity which tend to become visibly left skewed with time (with the appearance of a bimodal behavior at some times). The link between conductivity changes and reaction-driven processes promotes the emergence of non-Fickian effective transport features. The latter can be captured through a continuous-time random-walk model where solute travel times are approximated with a truncated power law probability distribution. The parameters of such a model shift towards values associated with increasingly high non-Fickian effective transport behavior as time progresses.


Regions of prevailing precipitation and dissolution are related to preferential flow patterns.

Large changes in non-Fickian transport parameters are observed while velocity variance displays modest variations.

Initial heterogeneity facilitates attaining an asymptotic average solute velocity value.

1 Introduction

Diagnosis and characterization of the feedback between geochemical precipitation/dissolution reactions and solute transport processes in heterogeneous subsurface systems is key to a variety of environmental and Earth science scenarios (Rege and Fogler, 1989; Berkowitz et al., 2016). A critical challenge is the emergence of complex dependencies between physical and chemical processes taking place across aquifer bodies (Saripalli et al., 2001). Heterogeneity of these systems promotes diverse patterns of precipitation and/or dissolution that may imprint a variety of dynamic system responses, including, for example, wormhole and oscillatory behaviors of system attributes such as porosity and permeability (Edery et al., 2011; Garing et al., 2015; Golfier et al., 2002). Examples of practical applications in this context include geologic CO2 storage (e.g., Pawar et al., 2015; Noiriel and Daval, 2017; Cabeza et al., 2020, and references therein), acid injection in production wells (e.g., Liu et al., 2017, and references therein), and reactive transport of contaminants (e.g., Ceriotti et al., 2018; Dalla Libera et al., 2020, and references therein).

Computational studies can assist the analysis of patterns of chemical transport across heterogeneous subsurface systems in the presence of precipitation/dissolution phenomena. While requiring an explicit description of the spatial heterogeneity of the system properties (Atchley et al., 2014), routine application of numerical simulations in practical settings is hampered by (i) our limited knowledge of the system attributes, resulting in uncertainty affecting the parameterization of the underlying physical and chemical processes and their variability, and (ii) the computational costs required to properly quantify such uncertainties and their propagation onto environmental quantities of interest. In this context, we rely on an effective approach to characterize the evolution of key features of solute transport in the presence of rock–fluid interactions across a porous medium whose spatially heterogeneous conductivity field is interpreted according to a commonly employed stochastic framework.

A critical element we tackle is related to the analysis of the dynamic feedback between reactive transport and spatially heterogeneous distributions of porous media attributes such as hydraulic conductivity. Following prior studies, we start by recognizing that, even under geochemical equilibrium conditions, the spatial heterogeneity of system attributes typically imprints an uneven spatial distribution of regions where chemical reactions take place, local fluctuations of conductivity being key to this element (Edery et al., 2016a). Further to this, our conceptualization of the setting is grounded on the observation that rendering of transport features in geological formations through effective formulations typically requires embedding non-Fickian features. To this end, we rely on an upscaled description of transport where solute travel/waiting times are approximated with a truncated power law probability density function (PDF), hereafter termed TPL (Berkowitz et al., 2006). This effective description is particularly relevant, because the emergence of non-Fickian transport features in heterogeneous formations has been observed at diverse scales of observation, including pore-, laboratory- and field-scale scenarios (e.g., Edery et al., 2011; Muljadi et al., 2018; Menke et al., 2018, and references therein).

In line with our objective, we rest on the framework of analysis developed in Edery et al. (2014, 2016a), where an effective depiction of transport processes is parameterized as a function of the statistics of solute residence times in randomly heterogeneous conductivity fields. A main element of this framework is that it yields a link between the TPL weighting times and the occurrence of preferential pathways that can be obtained from computational studies of transport in such conductivity fields (Edery et al., 2014). As a result, the methodology is conducive to an effective (or upscaled) representation of local features to identify signatures of non-Fickian transport (see also Dentz et al., 2011; Edery et al., 2014). To illustrate the main features associated with the scenario of interest, we consider a Darcy-scale formulation of a reactive transport setup, where precipitation and/or dissolution of minerals are driven by the injection of an acid compound establishing local equilibrium with the resident fluid and a solid matrix of the host porous medium which is considered to be composed of calcite mineral. While the geochemical processes we consider are somehow streamlined with respect to a field-scale scenario (see, for example, Lichtner, 1988; Dreybrodt et al., 1996), they embed the main elements characterizing the interplay between solute transport and rock–fluid interactions in Darcy-scale systems (e.g., Edery et al., 2011). Within this conceptual picture, our study aims at investigating (i) the interplay between the reactive process and the ensuing spatial distribution of preferential pathways associated with spatially heterogeneous conductivities and (ii) the link between locally occurring reaction-driven phenomena and emerging non-Fickian effective transport features, as captured by the TPL formulation of the PDF of particle travel/waiting times.

2 Methodology

2.1 Chemical model

We simulate a reactive transport scenario where calcite (CaCO3(s), subscript (s) denoting solid mineral) can dissolve or precipitate locally in the presence of chemical equilibrium between dissolved carbonic acid (H2CO3) and pH. The amount of dissolved H2CO3 as a function of pH (see Fig. S1 in the Supplement) is then governed by equilibrium conditions, which is tantamount to assuming a locally instantaneous reaction (i.e., as the reactive timescale approaches zero at equilibrium, this scenario corresponds to a local Damköhler number which tends to infinity). The formulation describing the chemical reactions can then be streamlined as


according to which two protons (H+) in Eq. (1b) react with CO32- to produce H2CO3 that in turn drives dissolution of the host calcium carbonate solid matrix. While simplified (see, for example, Lichtner, 1988; Dreybrodt et al., 1996), the chemical setup follows previous work by Edery et al. (2011), where there is an extensive assessment of the employed formulation (see also the Supplement). In this context and consistent with typical experimental practice, we consider the injected fluid and the porous medium to be associated with a source of H+ and an abundance of Ca2+, respectively. Thus, Ca2+ is not rate-limiting, and the spatial distribution of H+, as driven by transport and reaction, governs pH. The rate-limiting reaction is then Eq. (1b), that is controlled by the available H+ (or pH), similar to observations associated with other studies (Singurindy and Berkowitz, 2004; Edery et al., 2011). The chemical reaction system (Eqs. 1a and 1b) is here simplified (see, for example, Krauskopf and Bird, 1967) through

(2) co 2 h + c ,

where co denotes H2CO3, and h and c represent H+ and CaCO3(s), respectively.

2.2 Flow and transport modeling

Our computational setting is intended to mimic a laboratory-scale scenario where a 60 cm× 24 cm two-dimensional flow cell is filled with a porous system formed by a CaCO3(s) solid matrix. The system is initially fully saturated with water, and an injection of low-pH water takes place across the upstream side of the cell. To investigate the influence of the dissolution/precipitation reaction on solute transport, we consider a groundwater flow which is uniform in the mean, taking place within a two-dimensional domain where the (natural) logarithm of conductivity, y=ln (k), is considered a zero-mean, second-order stationary random field. The latter is further characterized by an isotropic, simple exponential, covariance function, with (normalized) correlation length l/L, with L being the length of the domain along the main flow direction. Various degrees of heterogeneity of the system are analyzed upon considering values of log-conductivity variance σ02=[1,3,5], with subscript 0 denoting that these values refer to the initially generated conductivity distributions (i.e., prior to the occurrence of reactions). The domain is discretized through 300 × 120 elements of uniform size Δ= 0.2 cm, yielding a field size of 60 cm× 24 cm. Each field is synthetically generated through the widely tested sequential Gaussian simulator GCOSIM3D (Gómez-Hernández and Journel, 1993) and is characterized by l/L= 0.016. This yields a value of Δ/l= 0.2, which is deemed adequate to capture the local features of the covariance of y and their impact on the main statistics of the velocity field and travel times (Ababou et al., 1989; Riva et al., 2009). We note that, while our study is representative of a laboratory-scale analysis, the dimensions of the domain have no particular implication, and they are only selected to ensure a meaningful description of the correlation structure that is included in the initially generated conductivity fields.

For each value of σ02, 20 random realizations of y are generated, each being then subject to a deterministic pressure drop (ΔH= 100 cm) between the inlet (left) and the outlet (right) sides. The local distribution of fluid velocity is computed through

(3) q ( x ) = 0 , q ( x ) = - k ( x ) h ( x ) ;

where q(x) is the local Darcy flux, and vector x corresponding to spatial location. The local fluid velocity field is then obtained as v=q/θ, with a constant initial porosity θ= 0.4 being here considered for the porous medium.

Solute transport is then simulated across each conductivity field by a particle-tracking approach (Le Borgne et al., 2008). A number of 105h particles (see Eq. 2), which is selected to represent a full pore volume (whose magnitude is evaluated through the initial condition, i.e, before porosity and permeability are altered by the reactive processes) at constant pH = 3.5, is divided by the domain length and multiplied by the mean velocity (v^, as evaluated from Eq. 4 across the whole domain). A total number of particles evaluated as Integer(105/LΔtv^) is then injected into the system at regular time intervals (Δt= 0.1 min). Particles are injected at the very beginning of each time step Δt and are flux-weighted according to the conductivity distribution at the inlet. The particles representing a full pore volume correspond to MH+= 10.79 moles of H+, the same amount being injected across the simulation course to obtain a constant pH = 3.5 in the injected fluid, while absence of h particles is taken to correspond to pH = 8. We then evaluate the pH value (or H+ molar mass) associated with each h particle by dividing the total number of H+ moles required to obtain a pH = 3.5 (i.e., 10.79 mol of H+) by the pore volume (as represented by 105h particles).

The upper and lower boundaries of the domain are reflective, while the outlet boundary is absorbing. Particle migration is simulated through

(4) d = v [ x ( t k ) ] δ t + d D ,

where d is particle displacement, x(tk) is the vector identifying spatial coordinates of particle location at time tk, v is fluid velocity at x(tk), δt=δs/v is the temporal displacement magnitude (v is the norm of v), and dD is the diffusive displacement. The latter is evaluated as dD=N[0,1]2Dmδt, where N[0,1] represents a two-dimensional vector of random variables, whose entries are mutually independent and sampled from a Gaussian distribution with zero mean and unit variance, with Dm= 10−5cm2 min−1 representing diffusion. The value of δs is selected to be an order of magnitude less than Δ, to accurately sample the velocity variability within a conductivity block. It is noted that while taking into account pore-scale processes within a continuum-scale model through a local-scale dispersion could be a modeling option, this would add an additional level of complexity to our numerical simulation without modifying the key elements of our work, which is focused on the interaction between flow patterns and reactive processes. In this context, our modeling choice is to represent macro-dispersive effects through averaging (in a multi-realization context) the effect of fluctuations of velocity arising between diverse realizations of the conductivity field. Consistent with this, we rely on a constant and isotropic diffusion coefficient that can be associated with an advection-dominated transport regime (as quantified in terms of a Péclet number, as seen in the following). This choice is also consistent with previous work (e.g., Aquino and Bolster, 2017; Wright et al., 2021).

Coupling between particle evolution and the geochemical setup illustrated in Sect. 2.1 is achieved in two steps. First, we advance all particles according to the displacement mechanism described above. Second, we satisfy the equilibrium condition (Eq. 2) by equilibrating both co and h within each cell, leading to precipitation or dissolution of a calcite mineral. The calcite volume to mole ratio is taken as MCaCO3= 37 cm3 mol−1 (Morse and Mackenzie, 1993), and the equilibrium between h and co particles (according to Eq. 2) leads to a local precipitation (or dissolution) of the solid. We update in time the spatial distribution of porosity assuming that it is characterized by a uniform change within each individual domain cell. We finally update conductivity through the Kozeny–Carman (KC) formulation:

(5) k ( ar ) i j = k ( br ) i j θ ( ar ) i j 3 1 - θ ( ar ) i j 2 1 - θ ( br ) i j 2 θ ( br ) i j 3 , θ ( ar ) i j = θ ( br ) i j ( 1 M CaCO 3 M H + ) ;

where k(ar)ij and θ(ar)ij are conductivity and porosity, respectively, after the reaction (ar) has taken place, while k(br)ij and θ(br)ij are their counterparts before the reaction is observed, with subscripts i and j being identifiers of a given cell. The process is repeated for each particle in each of the cells until an equilibrium between co and h is reached. We set an upper and a lower bound of 0.1 and 0.9, respectively, for porosity, to avoid the occurrence of unphysical porosity values. These constraints are set for consistency with our assumptions, i.e., to consider Darcy flow across the porous domain. A complete clogging (or opening) of a void space would require a different mathematical and conceptual treatment, which is beyond the scope of our study. Precipitation is treated numerically in a corresponding way.

We numerically calculate the updated local head and fluid velocity distributions from Eq. (4) at time intervals of 10 Δt, to reduce constraints associated with computational costs. The computational cost of each realization is between 1 to  3 d (depending on the value of σ02), upon relying on a 16-core Xeon 2.6 GHz processor with 64 GB RAM. With reference to the sensitivity of the results to the numerical parameters, we note that (a), when considering single realizations, the results showed only minute sensitivity to increasing the number of particles by a factor of 10 (less than 3 % difference in the results was observed); (b) results display a slightly larger sensitivity to the time step, otherwise resulting in a difference of less than about 5 % in the amount of reaction when decreasing the time interval by an order of magnitude; and (c) preliminary analyses aimed at assessing possible influences of the grid size on the key results of the study imbued us with confidence about the quality of the results obtained with the grid size employed in the study (details not shown), which we select as a good compromise between computational accuracy and execution time constraints, in light of our objectives.

The updated conductivity field is extracted and stored at the abovementioned regular intervals of 10 Δt. Transport of a non-reactive solute pulse is then simulated across each of these updated fields to capture the temporal evolution of the key parameters driving effective transport (see Sect. 2.3). While noting that natural porous media can exhibit complex relationships between permeability and porosity (Luquot and Gouze, 2009), which may not always be interpreted through the KC model (Eq. 6), we employ the latter formulation because it is considered a reference model in the literature and can serve as a proxy for alternative improved parameterizations (Erol et al., 2017). We also recall that we consider local equilibrium, this choice being consistent with the specific objective of this work which is related to the characterization of transport under dynamic evolution of continuum-scale quantities such as conductivity rather than being focused on a detailed characterization of the effects of the reactive process on the pore-scale structure (as reflected, for example, by local changes in specific surface area).

2.3 Quantities of interest

The workflow described in Sect. 2.2 enables one to extract computationally based quantities employed to characterize the analyzed reactive transport setup. As stated in Sect. 2.2, we simulate a tracer test within the original fields as well as within those modified by precipitation/dissolution. Particles are displaced through the action of advection and diffusion following a pulse (flux-weighted) injection at the inlet. These non-reactive transport simulations are performed to assess base values of parameters characterizing solute transport (a) prior to starting the reactive transport simulation as well as (b) at specific times after reaction changed the field. The empirical PDF of particle waiting times is assessed from the corresponding histogram starting by evaluating particle waiting times within a given domain cell through the inverse of the particle velocity computed at each time step multiplied by the cell length and weighted by the number of particles visiting the cell. This PDF is then used to estimate the parameters of the TPL model

(6) ψ ( t w ) = n t 1 exp ( - t / t 2 ) 1 + t t 1 1 + β ,

where tw is the waiting time of a particle within a given domain cell, and t1,t2, and β are model calibration parameters, which are estimated through a standard least square technique. Note that previous results have shown that the parameters obtained from Eq. (6) can be readily used to interpret breakthrough curves associated with non-reactive solutes (Edery et al., 2014).

Figure 1Relative frequency distributions f(yij) (blue circles) and f(nyij) (red circles) for a tracer test performed on the conductivity field prior to reaction and those associated with reactive simulation times of 10 and 20 min. Results correspond to σ02=1,3,5 (left, middle, and right columns, respectively) and to t= 0 (a–c), 10 (d–f), and 20 min (g–i). Mean and variance of these distributions are listed in Table 1.


The velocity fields are examined upon computing the evolution of the velocity and conductivity field statistics, as described in the following. Let us consider a discrete field of a generic quantity zij evaluated in a given cell ij. In the particle-tracking numerical simulations we quantify nij(t) as the number of particles that have visited cell ij along the simulation up to a given time t. Thus, we evaluate two relative frequency (or empirical probability) distributions, i.e., f(zij) and f(nzij), hereafter termed as unweighted and weighted distributions of the variable zij, respectively. We define the weighted variable nzij(t)=nij(t)zij/n(t), where n is the average value of nij. Note that the adopted weighting scheme corresponds to weighting zij by the solute mass distribution. Average values of the weighted and unweighted distributions (hereafter denoted as z and nz, respectively) can then be evaluated. In the following we perform particle weighting in the non-reactive as well as in the reactive transport scenarios. Distribution weighting by reactive particles is indicated by nR, meaning that weighting is performed based on the reactive transport simulations (i.e., considering h and co particles as explained above). The plain symbol n indicates weighting by non-reactive particles, employed to simulate conservative tracer tests as detailed above. The variable zij is taken to correspond to either the cell log-conductivity yij or fluid velocity vij in the results illustrated in Sect. 3.

Table 1Values of mean and variance of unweighted and weighted log-conductivity distributions and estimated parameters of the effected TPL model obtained through calibration of Eq. (6) against the computed distributions of particle waiting times in the domain cells. Results are listed for the three values of the initial log-conductivity variance (σ02) and are obtained from non-reactive transport simulations performed across conductivity fields resulting from reactive transport simulations at selected times.

Download Print Version | Download XLSX

Figure 2Sample and modeled probability density function ψ(tw) of particle waiting times for a tracer test performed on the conductivity field prior to reaction for those associated with reactive simulation times of 10 and 20 min. Results correspond to σ02=1,3,5 (left, middle, and right columns, respectively) and t=0 (a–c), 10 (d–f), and 20 min (g–i). Values of TPL parameters estimated by calibrating model Eq. (6) on the sample distributions are listed in Table 1.


3 Results

We start our analyses by simulating transport of a non-reactive tracer across the generated heterogeneous conductivity domains. As log-conductivity variance increases, the range of conductivity values naturally increases, this being reflected in the distribution f(yij) (see, for example, Fig. 1a–c (blue circles)). The shape of weighted conductivity distributions, f(nyij), differs from the one of f(yij), consistent with the observation that particles are chiefly channeled towards preferential flow pathways. The latter distributions tend to be shifted towards high conductivity values and are characterized by an enhanced mean conductivity value as compared against their generated (unweighted) counterparts (see conductivity mean and weighted mean values in Table 1 and the results corresponding to the blue and red circles depicted in Fig. 1a–c). This shift is imprinted onto the probability density function (PDF) of the waiting times and onto its associated TPL parameters (see Fig. 2a–c), consistent with prior studies (Edery et al., 2014, 2016b; Edery, 2021). We then simulate reactive transport across the collection of generated fields, allowing for precipitation (and/or dissolution) of calcite and assessing the evolution of the conductivity field according to the Kozeny–Carman formulation introduced in Sect. 2. Conductivity, head, and velocity fields, as well as particle visitations, nij(t), associated with species h and co are sampled across time.

Figure 3Heat map representing (a, b) log10(nRij), i.e., the number of h particles visiting each cell for σ02= 1 and 5, respectively, and (c–f) relative change in hydraulic conductivity at time t= 20 min (corresponding to the first pore volume) with respect to the initially generated values for σ02= 1 (c, e) and σ02= 5 (d, f). Panels (c) and (d) display positive changes in conductivity with respect to the initial field, while panels (e) and (f) display negative changes in conductivity, with both positive and negative changes being represented on a log scale. Results correspond to a selected realization of the log-conductivity fields. The highlighted box illustrates the separation between regions where precipitation or dissolution take place. Panels (g) and (h) display cells associated with a net decrease (green) and increase (red) of conductivity for cells outside (g) or within the PF (h) for σ02= 1.


After 200 Δt has elapsed (corresponding to a total simulation time of 20 min, i.e., a full pore volume), a set of h particles connecting the inlet to the outlet of the system is clearly visible (see Fig. 3a and b), these particles being non-uniformly distributed in space. Figure 3a and b depict a heat map of the h particles distribution at time t= 20 min (i.e., corresponding to the first pore volume), clearly evidencing the emergence of regions of preferential flow (PF). We also note that the number of h particles density (corresponding to concentration) tends to decrease with increasing distance from the inlet, these being replaced by co particles, consistent with the observation that they are consumed during the course of the reactive process which induces dissolution of the host solid matrix. The h and co particles attain equilibrium within cells away from the inlet. As such, reaction can only take place if a particle leaves (or enters) a cell under the action of advection and/or diffusion, leading to a new equilibrium state. When examining the alteration of conductivities due to the dissolution/precipitation reaction, we note that dissolution (corresponding to an increase in permeability values) is primarily tied to the preferential flow pathways. Otherwise, precipitation is seen to take place in regions close (on average) to these pathways. The highest strength of precipitation is observed in the proximity of the preferential pathways and decreases with distance from these.

Figure 3c and e depict the regions where conductivity has increased (due to dissolution) or decreased (due to precipitation), respectively. The h particles invading the domain closely follow the PFs displaying a fingering pattern, leading to a corresponding dissolution pattern associated with locally increased conductivities. Since conductivity values along the PFs are typically higher (on average), dissolution is increasing these conductivities even further, giving rise to a self-sustained enhancing mechanism. The concentration of h particles reaches a local (i.e., within a given cell) equilibrium with the produced co particles. Hence, dissolution will take place where transport induces shifts in concentration that need to be compensated by the dissolution/precipitation process to maintain local equilibrium. Such scenarios can be attained (i) by co particles exiting the preferential flow pathways due to the action of diffusion (i.e., they leave locations where concentration of h particles is large upon diffusing towards higher pH regions where they precipitate) or (ii) by h particles traveling through the fast preferential paths and advancing through these. Figure 3g and h display regions with dominating dissolution or precipitation for cells outside and within the PFs, respectively. Here cells associated with PFs are identified upon relying on particle visitations following (Edery et al., 2014). Dissolution dominates within the PFs (as indicated by the red cells in Fig. 3h), because h particles are injected through a flux-weighted boundary condition. On the other hand, the produced co particles do not precipitate at locations corresponding to the high h concentration residing in the PFs. These may precipitate away from these regions, where they experience low concentrations of h particles. Thus, we observe a reduction of conductivity taking place in regions adjacent to the PFs (Fig. 3b and g). In summary, our computational results document an increase in conductivity along the preferential pathways jointly with a conductivity reduction within regions close to these and along directions approximately normal to them.

Changes in conductivity values ensuing precipitation/dissolution are clearly visible by the broadening of the unweighted log-conductivity distribution f(yij) (see Fig. 1d–f and g–i (blue circles), evaluated at times t= 10 and t= 20 min, respectively). The reaction dynamics lead to a conductivity field characterized by a slightly increased average value, given that our computational analyses entail the injection of an acid fluid into the system (see Table 1 for details). Detailed inspection of Fig. 1d and f reveals that precipitation takes place across a slightly larger area than dissolution; that is, values of the frequency distribution f(yij) associated with low conductivities tend to increase at a larger rate rather than those corresponding to high conductivities (the left tail of the distributions becomes heavier than the right tail with the progress of reaction). The weighted (f(nyij)) and unweighted (f(yij)) distributions (red and blue circles in Fig. 1d–f and g–i at times t= 10 and t= 20 min, respectively) are visibly broadening, being associated with an average conductivity which is higher than one of the originally generated conductivity domains (see Table 1). As stated above, dissolution is focused along the preferential pathways, which comprise an area of limited extent with respect to the whole field.

The above-documented mechanism and its signature on the weighted and unweighted conductivity frequency distributions are sensitive to the initial log-conductivity variance, σ02. When considering both distributions f(yij) and f(nyij), associated with the case σ02=1, the distributions variance values are seen to increase in time, as compared to the values attained at the beginning of the simulation (i.e., prior to reaction; see Fig. 1a–d and g, as well as Table 1). Otherwise, as the initial heterogeneity increases (see, for example, σ02=3,5), mean and variance associated with the weighted and unweighted conductivity distributions display only minor changes (approximately 10 %) across the temporal window considered. The conductivity fields characterized by the lowest σ02 value are associated with preferential pathways that are not starkly recognizable when analyzed under non-reactive transport conditions. These channels become more clearly distinguishable as reactions induce an increase in the conductivities along the PFs. At the same time, precipitation causes a decrease in the conductivity outside the PF. This leads to an increased importance of the left tail of f(yij), corresponding to an increase in low conductivity values (see Fig. 1, left, middle, and bottom rows).

With reference to the highest conductivity variance analyzed, the reaction patterns for the precipitation and dissolution lead to a smaller relative change between conductivity frequency distributions evaluated prior and after the reaction. Relative changes between the unweighted and weighted conductivity frequency distributions (including the ensuing mean and variance listed in Table 1) evaluated before and after the reaction are less pronounced as the variance in the generated conductivity field increases (see Fig. 1 (middle and right columns) and Table 1). This is related to the observation that, as log-conductivity variance increases, preferential pathways in the originally generated field become markedly more distinct. Thus, relative differences between unweighted and weighted conductivity histograms are seen to diminish in time, because the flow field is already organized according to well-identified pathways and tends to preserve its initial pattern (Fig. 1d–f and g–i at times t= 10 and t= 20 min, respectively). Note that low-order statistics (i.e., mean and variance) of velocity and conductivity display only a minute evolution with the progress of reaction, in spite of the relevant changes exhibited by the tails of the frequency distributions (see Fig. 1) for all considered values of σ02, with the latter feature being relevant when addressing non-Fickian transport, as further discussed below.

As stated in Sect. 2.2, the conductivity fields altered through precipitation/dissolution and extracted at regular time intervals of 10 Δt are subject to non-reactive transport analyses, and the ensuing evolution of the parameters of the TPL model (Eq. 6) is analyzed. Key results of these analyses are listed in Table 1 with reference to the original (unaltered) conductivity fields and at the final simulation time (i.e., at time t= 200 Δt). Analysis of the results associated with transport across the log-conductivity field characterized by the smallest original variance (i.e., σ02= 1) and listed in Table 1 indicates that the changes in the sample log-conductivity PDF induced by the progress of the reaction are reflected by the parameters of the TPL model (Eq. 6). These transition from estimated values corresponding to an effective Fickian transport regime (corresponding to β=2; see also Fig. 2a) to values denoting a highly non-Fickian effective transport setting, manifested by the widening of the support of the waiting time PDF ψ(tw) (see also Fig. 2d and g for results obtained at t= 10 and 20 min, respectively). Effective transport in the domain with the highest variance (i.e., σ02= 5; see Table 1) is characterized by estimated TPL parameters corresponding to a non-Fickian signature also prior to the occurrence of precipitation/dissolution (Fig. 2c). Such a signature is then further enhanced after reaction has altered the conductivity field yet displaying a less marked evolution of the shape of the ψ(tw) as compared to the case σ02= 1 (see also Fig. 2f and i for t= 10 and 20 min, respectively).

Figure 4Temporal evolution of (a) the weighted mean velocity nRv (open symbols) and (b) the sum of all conductivity changes over 1 min. Results correspond to σ02= 1 (circles), 3 (squares), and 5 (diamonds).


The observed temporal changes in conductivity and the ensuing local dynamics of transport pattern yield global variations in the reaction rate. Consistent with prior studies and the imposed boundary conditions, the mean velocity associated with the originally generated conductivity domains increases with σ02. As the reaction progresses and the conductivity fields change, the increased area subject to dissolution leads to a slight increase in the mean velocity for all of the σ02 analyzed. To analyze the influence of the preferential flow on the velocity that is affecting particle transport, we consider the average value nRv, evaluated upon considering weighting by the number of reactive particles, nR, visiting each cell (where the term reactive particles denotes both h and co particles employed in the context of the reactive transport simulations). The weighted average velocity displays an initial increase over time due to the increase in conductivity within the preferential pathways. When considering the relative change across the whole simulation time, values of the temporal increase in nRv are similar across the three heterogeneity levels examined, i.e., they are seemingly independent of σ02. However, results in Fig. 4a also reveal that the average velocity nRv displays distinct temporal histories depending on σ02. In particular, the value of nRv tends to attain an asymptotic value at time t 7 min for σ02=5, while showing a sustained increasing trend for σ02=1. This result suggests that the feedback between reaction and flow patterns reaches an asymptotic condition faster in systems characterized by higher heterogeneity.

The temporal evolution of the velocity fields due to the precipitation/dissolution reaction and the resulting conductivity changes lead to a time-dependent reaction pattern and reaction rates. The Damköhler number is infinite on a local scale in our computational analyses, because the reaction is instantaneous. When considering the entire system, transport processes induce a net overall reaction rate that can be quantified as the sum of the total conductivity changes across time (Sum(|Δkij|) [cm s−1]). The latter incorporates both positive and negative changes in hydraulic conductivity and therefore quantifies the overall intensity of precipitation and dissolution processes in the domain. The quantity Sum(|Δkij|) is evaluated across all realizations for each of the σ02 values considered and is depicted in Fig. 4b as a function of time. These results indicate that the overall reaction rate increases in time with a similar rate for all considered values of σ02 (Fig. 4b) at early times. The observed increase is consistent with the initial advancement of the reaction front across the domain. We observe that the reactive process magnitude is proportional to σ02. For low initial levels of heterogeneity, conductivity values along the preferential pathways are closer to the average field conductivity than what can be observed for the highly heterogeneous domains. As such, the portion of the domain where precipitation or dissolution can take place increases at a rate proportional to σ02. As the reaction front reaches the domain outlet, the dissolving front found in the PF leaves the domain. Hence, the global variation in conductivity (which is proportional to the magnitude of reactive processes) tends towards an asymptotic value, corresponding to the diffusion-controlled solute exchange along a direction transverse to the preferential pathways (see Fig. 3). In agreement with results shown in Fig. 4a, this transition towards an asymptotic regime takes place earlier for larger values of σ02, while a smaller initial heterogeneity implies a longer transient period.

4 Conclusions

Our computational study tackles the quantitative characterization of the feedbacks between precipitation and dissolution reaction dynamics taking place in randomly heterogeneous conductivity fields associated with various degrees of spatial heterogeneity. Our work leads to the following key conclusions.

Joint occurrence of precipitation and dissolution is tightly coupled with the existence of preferential flow pathways. Conductivity increase due to the dissolution reaction along such paths leads to enhanced particle migration along these. The dominance of preexisting preferential flow regions on the (reactive) transport pattern across the field is therefore further reinforced and self-sustained across time. At the same time, diffusion promotes displacement of particles, leading to precipitation (and hence a progressive reduction over time of local conductivities) at locations in the proximity of these.

Reactive processes yield an increase over time of the range of conductivity values across the domain, eventually leading to a widening of the support of solute waiting times and conductivity distributions. The clear separation between regions where dissolution or precipitation takes place is reflected in sample distributions of conductivity which tend to become visibly left skewed with time, a feature which is associated with precipitation taking place in low-conductivity cells located in the proximity of existing preferential flow pathways.

We weight the conductivity and velocity distributions with the solute mass, to characterize the probability density function of particle travel/waiting times, which form the TPL model parameters, thus enabling us to capture effective (upscaled) non-Fickian transport behaviors. With the progress of precipitation/dissolution reactions, transport shifts towards an increasingly acute non-Fickian effective behavior (see Fig. 2 and ensuing parameter values listed in Table 1). The latter is then seen as a direct outcome of the documented feedbacks between transport and reactions taking place in heterogeneous porous media. The evolution of TPL model parameters towards a pronounced non-Fickian behavior is associated with only minor changes in the mean and variance of log-conductivity values. This result is consistent with the conceptual picture that the tails of flux and hydraulic conductivity distributions carry critical information to characterize transport while displaying only a minor effect on low-order statistics associated with these quantities. Our results suggest that this feature must be acknowledged to properly characterize transport in the presence of precipitation/dissolution. In this context, we recall that our study is not aimed at exploring the skill of any specific transport model to interpret non-reactive transport across the conductivity fields prior and/or after reaction takes place, with this particular analysis being deferred to a future study.

We observe the emergence of an asymptotic regime in highly heterogeneous systems, where the (averaged) solute velocity attains a constant value even in the presence of reaction. This suggests the occurrence of an equilibrium state between reactive processes and transport under the flow conditions analyzed. This regime is attained because the effects of locally occurring precipitation and dissolution balance each other at the overall scale of the system, so the ensuing (ensemble-averaged) solute velocity remains unaffected. The time required to attain such an asymptotic state increases with decreasing initial heterogeneity of the conductivity field, thus suggesting that pre-asymptotic behaviors may be more relevant in initially (i.e., prior to reaction taking place) homogeneous systems.

Our results are based on numerical simulations and may be used to inform upscaling approaches to capture the pre-asymptotic and asymptotic dynamics of reactive transport in heterogeneous systems through simplified models. Future computational studies might also include an assessment of the importance of initial/boundary conditions and/or of solute injection mode on the emergence of non-Fickian transport features, with special emphasis on locations close to the inlet boundary which might then have an impact on the overall solute residence times and related statistics.

Code and data availability

The code and data will be made available using a shared folder upon request.


The supplement related to this article is available online at:

Author contributions

YE performed numerical simulations, contributed to designing the research methodology, analyzed data, created figures, and wrote the first draft; MS analyzed data, performed numerical simulations, and created figures; GP contributed to designing the research methodology, discussed the results, and contributed to the writing; AG contributed to designing the research methodology, discussed the results, and contributed to the writing.

Competing interests

Some authors are members of the editorial board of Hydrology and Earth System Sciences. The peer-review process was guided by an independent editor, and the authors have also no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Yaniv Edery gives thanks to the support of the German-Israeli Foundation (grant no. I-2536-306.8) and the Israel Science Foundation (grant no. 801/20).

Financial support

This research has been supported by the German-Israeli Foundation for Scientific Research and Development (grant no. I-2536-306.8) and the Israel Science Foundation (grant no. 801/20).

Review statement

This paper was edited by Erwin Zehe and reviewed by Philippe Ackerer and one anonymous referee.


Ababou, R., McLaughlin, D., Gelhar, L. W., and Tompson, A. F. B.: Numerical simulation of three-dimensional saturated flow in randomly heterogeneous porous media, Transport Porous Med., 4, 549–565,, 1989. 

Aquino, T. and Bolster, D.: Localized point mixing rate potential in heterogeneous velocity fields, Transport Porous Med., 119, 391–402,, 2017. 

Atchley, A. L., Navarre-Sitchler, A. K., and Maxwell, R. M.: The effects of physical and geochemical heterogeneities on hydro-geochemical transport and effective reaction rates, J. Contam. Hydrol., 165, 53–64,, 2014. 

Berkowitz, B., Cortis, A., Dentz, M., and Scher, H.: Modeling non-Fickian transport in geological formations as a continuous time random walk, Rev. Geophys., 44, RG2003,, 2006. 

Berkowitz, B., Dror, I., Hansen, S. K., and Scher, H.: Measurements and models of reactive transport in geological media, Rev. Geophys., 54, 930–986,, 2016. 

Cabeza, Y., Hidalgo, J. J., and Carrera, J.: Competition Is the Underlying Mechanism Controlling Viscous Fingering and Wormhole Growth, Geophys. Res. Lett., 47, e2019GL084795,, 2020. 

Ceriotti, G., Guadagnini, L., Porta, G., and Guadagnini, A.: Local and Global Sensitivity Analysis of Cr (VI) Geogenic Leakage Under Uncertain Environmental Conditions, Water Resour. Res., 54, 5785–5802,, 2018. 

Dalla Libera, N., Pedretti, D., Tateo, F., Mason, L., Piccinini, L., and Fabbri, P.: Conceptual Model of Arsenic Mobility in the Shallow Alluvial Aquifers Near Venice (Italy) Elucidated Through Machine Learning and Geochemical Modeling, Water Resour. Res., 56, e2019WR026234,, 2020. 

Dentz, M., Gouze, P., and Carrera, J.: Effective non-local reaction kinetics for transport in physically and chemically heterogeneous media, J. Contam. Hydrol., 120–121, 222–236,, 2011. 

Dreybrodt, W., Lauckner, J., Zaihua, L., Svensson, U., and Buhmann, B.: The kinetics of the reaction CO2 + H2O → H+ + HCO3 as one of the rate limiting steps for the dissolution of calcite in the system H2O-CO2-CaCO3, Geochim. Cosmochim. Ac., 60, 3375–3381,, 1996. 

Edery, Y.: The Effect of varying correlation lengths on Anomalous Transport, Transport Porous Med., 137, 345–364,, 2021. 

Edery, Y., Scher, H., and Berkowitz, B.: Dissolution and precipitation dynamics during dedolomitization, Water Resour. Res., 47, W08535,, 2011. 

Edery, Y., Guadagnini, A., Scher, H., and Berkowitz, B.: Origins of anomalous transport in heterogeneous media: Structural and dynamic controls, Water Resour. Res., 50, 1490–1505,, 2014. 

Edery, Y., Porta, G. M., Guadagnini, A., Scher, H., and Berkowitz, B.: Characterization of Bimolecular Reactive Transport in Heterogeneous Porous Media, Transport Porous Med., 115, 291–310,, 2016a. 

Edery, Y., Geiger, S., and Berkowitz, B.: Structural controls on anomalous transport in fractured porous rock, Water Resour. Res., 52, 5634–5643,, 2016b. 

Erol, S., Fowler, S. J., Harcouët-Menou, V., and Laenen, B.: An Analytical Model of Porosity–Permeability for Porous and Fractured Media, Transport Porous Med., 120, 327–358,, 2017. 

Garing, C., Gouze, P., Kassab, M., Riva, M., and Guadagnini, A.: Anti-correlated Porosity–Permeability Changes During the Dissolution of Carbonate Rocks: Experimental Evidences and Modeling, Transport Porous Med., 107, 595–621,, 2015. 

Golfier, F., Zarcone, C., Bazin, B., Lenormand, R., Lasseux, D., and Quintard, M.: On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium, J. Fluid Mech., 457, 213–254,, 2002. 

Gómez-Hernández, J. J. and Journel, A. G.: Joint Sequential Simulation of MultiGaussian Fields, in: Geostatistics Tróia '92, vol. 5, edited by: Soares, A., Springer Netherlands, Dordrecht, 85–94,, 1993. 

Krauskopf, K. B. and Bird, D. K.: Introduction to geochemistry, Internat. 3 Rev. edn., McGraw Hill, New York, 640 pp., 1967. 

Le Borgne, T., Dentz, M., and Carrera, J.: Lagrangian Statistical Model for Transport in Highly Heterogeneous Velocity Fields, Phys. Rev. Lett., 101, 090601,, 2008. 

Lichtner, P.: The quasi-stationary state approximation to coupled mass transport and fluid-rock interaction in a porous medium, Geochim. Cosmochim. Ac., 52, 143–165,, 1988. 

Liu, P., Yao, J., Couples, G. D., Ma, J., and Iliev, O.: 3-D Modelling and Experimental Comparison of Reactive Flow in Carbonates under Radial Flow Conditions, Sci. Rep.-UK, 7, 17711,, 2017. 

Luquot, L. and Gouze, P.: Experimental determination of porosity and permeability changes induced by injection of CO2 into carbonate rocks, Chem. Geol., 265, 148–159,, 2009. 

Menke, H. P., Reynolds, C. A., Andrew, M. G., Pereira Nunes, J. P., Bijeljic, B., and Blunt, M. J.: 4D multi-scale imaging of reactive flow in carbonates: Assessing the impact of heterogeneity on dissolution regimes using streamlines at multiple length scales, Chem. Geol., 481, 27–37,, 2018. 

Morse, J. W. and Mackenzie, F. T.: Geochemical constraints on CaCO3 transport in subsurface sedimentary environments, Chem. Geol., 105, 181–196,, 1993. 

Muljadi, B. P., Bijeljic, B., Blunt, M. J., Colbourne, A., Sederman, A. J., Mantle, M. D., and Gladden, L. F.: Modelling and upscaling of transport in carbonates during dissolution: Validation and calibration with NMR experiments, J. Contam. Hydrol., 212, 85–95,, 2018. 

Noiriel, C. and Daval, D.: Pore-Scale Geochemical Reactivity Associated with CO2 Storage: New Frontiers at the Fluid–Solid Interface, Accounts Chem. Res., 50, 759–768,, 2017. 

Pawar, R. J., Bromhal, G. S., Carey, J. W., Foxall, W., Korre, A., Ringrose, P. S., Tucker, O., Watson, M. N., and White, J. A.: Recent advances in risk assessment and risk management of geologic CO2 storage, International Journal of Greenhouse Gas Control, 40, 292–311,, 2015. 

Rege, S. D. and Fogler, H. S.: Competition among flow, dissolution, and precipitation in porous media, AIChE J., 35, 1177–1185,, 1989. 

Riva, M., Guadagnini, A., Neuman, S. P., Janetti, E. B., and Malama, B.: Inverse analysis of stochastic moment equations for transient flow in randomly heterogeneous media, Adv. Water Resour., 32, 1495–1507,, 2009. 

Saripalli, K. P., Meyer, P. D., Bacon, D. H., and Freedman, V. L.: Changes in Hydrologic Properties of Aquifer Media Due to Chemical Reactions: A Review, Crit. Rev. Env. Sci. Tec., 31, 311–349,, 2001. 

Singurindy, O. and Berkowitz, B.: Dedolomitization and flow in fractures, Geophys. Res. Lett., 31, L24501,, 2004. 

Wright, E. E., Sund, N. L., Richter, D. H., Porta, G. M., and Bolster, D.: Upscaling bimolecular reactive transport in highly heterogeneous porous media with the LAgrangian Transport Eulerian Reaction Spatial (LATERS) Markov model, Stoch. Env. Res. Risk A., 35, 1529–1547,, 2021. 

Short summary
The interplay between dissolution, precipitation and transport is widely encountered in porous media, from CO2 storage to cave formation in carbonate rocks. We show that dissolution occurs along preferential flow paths with high hydraulic conductivity, while precipitation occurs at locations close to yet separated from these flow paths, thus further funneling the flow and changing the probability density function of the transport, as measured on the altered conductivity field at various times.