**Research article**
11 Nov 2021

**Research article** | 11 Nov 2021

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

Yaniv Edery Martin Stolar Giovanni Porta and Alberto Guadagnini

^{1},

^{1},

^{2},

^{2}

**Yaniv Edery et al.**Yaniv Edery Martin Stolar Giovanni Porta and Alberto Guadagnini

^{1},

^{1},

^{2},

^{2}

^{1}Faculty of Civil and Environmental Engineering, Technion, Haifa, Israel^{2}Department of Civil and Environmental Engineering, Politecnico di Milano, 20133, Milan, Italy

^{1}Faculty of Civil and Environmental Engineering, Technion, Haifa, Israel^{2}Department of Civil and Environmental Engineering, Politecnico di Milano, 20133, Milan, Italy

**Correspondence**: Yaniv Edery (yanivedery@technion.ac.il)

**Correspondence**: Yaniv Edery (yanivedery@technion.ac.il)

Received: 02 May 2021 – Discussion started: 17 May 2021 – Revised: 12 Aug 2021 – Accepted: 15 Sep 2021 – Published: 11 Nov 2021

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.

- Article
(1801 KB) -
Supplement
(149 KB) - BibTeX
- EndNote

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.

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 CO_{2} 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.1 Chemical model

We simulate a reactive transport scenario where calcite (CaCO_{3(s)}, subscript (s) denoting solid mineral) can dissolve or precipitate
locally in the presence of chemical equilibrium between dissolved carbonic acid (H_{2}CO_{3}) and pH. The amount of dissolved
H_{2}CO_{3} 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 ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ to produce H_{2}CO_{3} 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 Ca^{2+},
respectively. Thus, Ca^{2+} 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

where co denotes H_{2}CO_{3}, and *h* and *c* represent H^{+} and CaCO_{3(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 CaCO_{3(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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}=[\mathrm{1},\mathrm{3},\mathrm{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 $\mathrm{\Delta}/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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$, 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

where ** q**(

**) is the local Darcy flux, and vector**

*x***corresponding to spatial location. The local fluid velocity field is then obtained as $\mathit{v}=\mathbf{q}/\mathit{\theta}$, with a constant initial porosity**

*x**θ*= 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 10^{5} *h*
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 ($\widehat{v}$, as evaluated from Eq. 4 across the whole domain). A total number of particles evaluated as
$\text{Integer}({\mathrm{10}}^{\mathrm{5}}/L\cdot \mathrm{\Delta}t\cdot \widehat{v})$ 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 ${M}_{{\mathrm{H}}^{+}}$ = 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 10^{5} *h* particles).

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

where ** d** is particle displacement,

**(**

*x**t*

_{k}) is the vector identifying spatial coordinates of particle location at time

*t*

_{k},

**is fluid velocity at**

*v***(**

*x**t*

_{k}),

*δ*

*t*= $\mathit{\delta}s/v$ is the temporal displacement magnitude (

*v*is the norm of

*v*), and

*d*_{D}is the diffusive displacement. The latter is evaluated as

*d*

_{D}= $\mathit{N}[\mathrm{0},\mathrm{1}]\sqrt{\mathrm{2}{D}_{\mathrm{m}}\mathit{\delta}t}$, where

**[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**

*N**D*

_{m}= 10

^{−5}cm

^{2}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 ${M}_{{\mathrm{CaCO}}_{\mathrm{3}}}$ = 37 cm^{3} 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:

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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$), 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

where *t*_{w} is the waiting time of a particle within a given domain cell, and *t*_{1},*t*_{2}, 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).

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 *z*_{ij} evaluated in a given cell *i**j*. In the particle-tracking numerical simulations we quantify
*n*_{ij}(*t*) as the number of particles that have visited cell *i**j* along the simulation up to a given time *t*. Thus, we evaluate two relative
frequency (or empirical probability) distributions, i.e., *f*(*z*_{ij}) and *f*(*n**z*_{ij}), hereafter termed as unweighted and weighted distributions of
the variable *z*_{ij}, respectively. We define the weighted variable *n**z*_{ij}(*t*) = ${n}_{ij}\left(t\right){z}_{ij}/\stackrel{\mathrm{\u203e}}{n}\left(t\right)$, where $\stackrel{\mathrm{\u203e}}{n}$ is
the average value of *n*_{ij}. Note that the adopted weighting scheme corresponds to weighting *z*_{ij} by the solute mass distribution. Average
values of the weighted and unweighted distributions (hereafter denoted as $\stackrel{\mathrm{\u203e}}{z}$ and $\stackrel{\mathrm{\u203e}}{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 *n*_{R}, 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 *z*_{ij} is taken to correspond to either the cell log-conductivity *y*_{ij} or fluid velocity *v*_{ij}
in the results illustrated in Sect. 3.

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*(*y*_{ij}) (see, for example, Fig. 1a–c
(blue circles)). The shape of weighted conductivity distributions, *f*(*n**y*_{ij}), differs from the one of *f*(*y*_{ij}), 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, *n*_{ij}(*t*), associated with species *h* and co are
sampled across time.

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*(*y*_{ij}) (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*(*y*_{ij}) 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*(*n**y*_{ij})) and unweighted (*f*(*y*_{ij})) 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, ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$. When considering both distributions *f*(*y*_{ij}) and *f*(*n**y*_{ij}), associated with the case
${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}=\mathrm{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, ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}=\mathrm{3},\mathrm{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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ 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*(*y*_{ij}), 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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$, 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., ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ = 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 *ψ*(*t*_{w}) (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., ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ = 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 *ψ*(*t*_{w}) as compared to the case ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ = 1 (see also Fig. 2f and i for
*t* = 10 and 20 min, respectively).

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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$. 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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ analyzed. To analyze the influence of the preferential flow on the
velocity that is affecting particle transport, we consider the average value $\stackrel{\mathrm{\u203e}}{{n}_{\mathrm{R}}v}$, evaluated upon considering weighting by the
number of reactive particles, *n*_{R}, 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 $\stackrel{\mathrm{\u203e}}{{n}_{\mathrm{R}}v}$ are similar across the three heterogeneity levels examined, i.e., they are seemingly independent of
${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$. However, results in Fig. 4a also reveal that the average velocity $\stackrel{\mathrm{\u203e}}{{n}_{\mathrm{R}}v}$ displays distinct temporal histories
depending on ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$. In particular, the value of $\stackrel{\mathrm{\u203e}}{{n}_{\mathrm{R}}v}$ tends to attain an asymptotic value at time
*t* ≈ 7 min for ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}=\mathrm{5}$, while showing a sustained increasing trend for ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}=\mathrm{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(|Δ*k*_{ij}|) [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(|Δ*k*_{ij}|) is evaluated across all realizations for each of the
${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ 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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ (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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$. 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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$. 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 ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$, while a smaller initial heterogeneity
implies a longer transient period.

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.

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

The supplement related to this article is available online at: https://doi.org/10.5194/hess-25-5905-2021-supplement.

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.

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).

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).

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, https://doi.org/10.1007/BF00223627, 1989.

Aquino, T. and Bolster, D.: Localized point mixing rate potential in heterogeneous velocity fields, Transport Porous Med., 119, 391–402, https://doi.org/10.1007/s11242-017-0887-z, 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, https://doi.org/10.1016/j.jconhyd.2014.07.008, 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, https://doi.org/10.1029/2005RG000178, 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, https://doi.org/10.1002/2016RG000524, 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, https://doi.org/10.1029/2019GL084795, 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, https://doi.org/10.1029/2018WR022857, 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, https://doi.org/10.1029/2019WR026234, 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, https://doi.org/10.1016/j.jconhyd.2010.06.002, 2011.

Dreybrodt, W., Lauckner, J., Zaihua, L., Svensson, U., and Buhmann, B.: The kinetics of the reaction CO_{2} + H_{2}O → H^{+} + HCO${}_{\mathrm{3}}^{\mathrm{-}}$ as one of the rate limiting steps for the dissolution of calcite in the system H_{2}O-CO_{2}-CaCO_{3}, Geochim. Cosmochim. Ac., 60, 3375–3381, https://doi.org/10.1016/0016-7037(96)00181-0, 1996.

Edery, Y.: The Effect of varying correlation lengths on Anomalous Transport, Transport Porous Med., 137, 345–364, https://doi.org/10.1007/s11242-021-01563-9, 2021.

Edery, Y., Scher, H., and Berkowitz, B.: Dissolution and precipitation dynamics during dedolomitization, Water Resour. Res., 47, W08535, https://doi.org/10.1029/2011WR010551, 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, https://doi.org/10.1002/2013WR015111, 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, https://doi.org/10.1007/s11242-016-0684-0, 2016a.

Edery, Y., Geiger, S., and Berkowitz, B.: Structural controls on anomalous transport in fractured porous rock, Water Resour. Res., 52, 5634–5643, https://doi.org/10.1002/2016WR018942, 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, https://doi.org/10.1007/s11242-017-0923-z, 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, https://doi.org/10.1007/s11242-015-0456-2, 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, https://doi.org/10.1017/S0022112002007735, 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, https://doi.org/10.1007/978-94-011-1739-5_8, 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, https://doi.org/10.1103/PhysRevLett.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, https://doi.org/10.1016/0016-7037(88)90063-4, 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, https://doi.org/10.1038/s41598-017-18095-2, 2017.

Luquot, L. and Gouze, P.: Experimental determination of porosity and permeability changes induced by injection of CO_{2} into carbonate rocks, Chem. Geol., 265, 148–159, https://doi.org/10.1016/j.chemgeo.2009.03.028, 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, https://doi.org/10.1016/j.chemgeo.2018.01.016, 2018.

Morse, J. W. and Mackenzie, F. T.: Geochemical constraints on CaCO3 transport in subsurface sedimentary environments, Chem. Geol., 105, 181–196, https://doi.org/10.1016/0009-2541(93)90125-3, 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, https://doi.org/10.1016/j.jconhyd.2017.08.008, 2018.

Noiriel, C. and Daval, D.:
Pore-Scale Geochemical Reactivity Associated with CO_{2} Storage: New Frontiers at the Fluid–Solid Interface,
Accounts Chem. Res.,
50, 759–768, https://doi.org/10.1021/acs.accounts.7b00019, 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 CO_{2} storage, International Journal of Greenhouse Gas Control, 40, 292–311, https://doi.org/10.1016/j.ijggc.2015.06.014, 2015.

Rege, S. D. and Fogler, H. S.: Competition among flow, dissolution, and precipitation in porous media, AIChE J., 35, 1177–1185, https://doi.org/10.1002/aic.690350713, 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, https://doi.org/10.1016/j.advwatres.2009.07.003, 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, https://doi.org/10.1080/20016491089244, 2001.

Singurindy, O. and Berkowitz, B.: Dedolomitization and flow in fractures, Geophys. Res. Lett., 31, L24501, https://doi.org/10.1029/2004GL021594, 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, https://doi.org/10.1007/s00477-021-02006-z, 2021.

_{2}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.