资源环境科技发展态势分析平台

Anomalous Transport in Dissolving Porous Media: Transitions Between Fickian and Non-Fickian Regimes

https://doi.org/10.1029/2025GL115940
2025-07-30
查看原文
Geophysical Research Letters . Volume 52 , issue 15
Jingxuan Deng, Rishabh P. Sharma, Piotr Szymczak, Peter K. Kang

Abstract

Mineral dissolution is a fundamental geologic process that alters pore structures, significantly impacting fluid flow and solute transport in porous media. Depending on the interplay between advection, diffusion, and reaction rates, mineral dissolution produces distinct dissolution patterns, such as wormholing and uniform dissolution. These structural changes directly influence the flow field, which in turn controls solute transport behavior. In this study, we conducted pore network modeling to investigate the effects of initial pore network heterogeneity and dissolution regimes on solute transport dynamics. Wormholing enhances network heterogeneity by creating preferential flow paths and stagnation zones, resulting in a transition from Fickian to non-Fickian transport. In contrast, uniform dissolution homogenizes the pore network and the flow field, driving a transition from non-Fickian to Fickian transport, even in networks with high initial heterogeneity. These transitions are governed by initial network heterogeneity and the Damköhler number.

Key Points

  • Pore network modeling shows the distinctive impacts of dissolution regimes on solute transport behavior

  • Dissolution in the wormholing regime leads to non-Fickian transport, while dissolution in the uniform regime results in Fickian transport

  • The transition between transport regimes is governed by initial network heterogeneity and Damköhler number

Plain Language Summary

Mineral dissolution reshapes the internal structure of porous materials, creating distinct patterns such as wormholes. This process drives key geological phenomena, including the formation of caves, sinkholes, and subsurface drainage systems characteristic of karst terrains. The resulting structural changes significantly impact fluid flow and solute transport through porous materials. In this study, we investigated how the initial pore structure and two specific dissolution regimes—wormholing and uniform dissolution—affect solute movement. Our findings reveal that wormholing, through the development of preferential flow paths, gives rise to complex transport patterns (non-Fickian behavior). In contrast, uniform dissolution reduces structural heterogeneity, leading to more predictable solute movement (Fickian behavior), even in initially highly heterogeneous systems. We show that the transition between these transport behaviors is controlled by the interplay between initial pore structure and the prevailing dissolution regime. These findings improve our understanding of solute transport in subsurface systems, with important implications for managing groundwater resources, storing carbon dioxide underground, and mitigating environmental contamination.

1 Introduction

Mineral dissolution is a geologic process that occurs across diverse environments, and often governs biogeochemical and hydrogeological processes (Schuler et al., 2024; Seiffert et al., 2014; Wild et al., 2022). The dissolution of host rock alters structural heterogeneity and flow patterns. Dissolution can proceed in a spatially uniform or non-uniform manner, depending on the interplay between advection, diffusion, and reaction rates. The dynamics of dissolution are commonly characterized by the Damköhler number (Da) or the reactant penetration length (lp). When the reactant penetration length is smaller than the system length (large Da), dissolution becomes concentrated along the most conductive flow paths, forming wormholes that extend deep into the porous medium. These wormholes serve as preferential pathways for reactants, enabling deeper penetration into the system. In contrast, when the penetration length is much longer than the system length (small Da), reactants are able to penetrate the entire system before being consumed. This results in uniform dissolution, which homogenizes the porous medium (Golfier et al., 2002; Hoefner & Fogler, 1988; Roded et al., 2020; Szymczak & Ladd, 2009).

Structural heterogeneity in rocks, such as variations in pore size distribution, plays an important role in controlling solute transport behavior. In homogeneous media, where the flow field is generally uniform, overall solute transport can be effectively modeled using an average advective velocity combined with Gaussian spreading. This process is known as Fickian transport and can be well described by the classical advection-dispersion equation. However, in heterogeneous media, variations in permeability induce velocity heterogeneity, leading to preferential flow pathways and stagnant zones. Solutes preferentially move through high-permeability regions, while some become trapped for long periods in low-permeability zones. This results in anomalous, or non-Fickian transport behavior, which cannot be adequately captured by traditional advection-dispersion models (Berkowitz et al., 2006; Dentz et al., 2004; P. K. Kang, Dentz, et al., 2015). The solute breakthrough curve (BTC), which represents the concentration of solute arriving at a downstream location over time, is a key measure for distinguishing between Fickian and non-Fickian transport regimes (Berkowitz et al., 2000, 2006; Edery et al., 2014; Margolin & Berkowitz, 2000). Non-Fickian transport is often manifested by early arrivals and power-law late-time tailing, while Fickian transport typically produces an inverse Gaussian distribution without power-law tailing (Dentz et al., 2004; Margolin & Berkowitz, 2000; Scher & Montroll, 1975). Understanding the mechanisms underlying Fickian and non-Fickian transport is vital for numerous applications that rely on accurate prediction of solute transport. These include geological carbon sequestration (Q. Kang et al., 2009; Matter & Kelemen, 2009), groundwater remediation (Jeen et al., 2011; Kim et al., 2024; Mayer et al., 2001), and in situ mining of critical minerals (Ranjith et al., 2017; Seredkin et al., 2016). However, most studies on anomalous transport have primarily focused on cases where the structural heterogeneity of the rock remains constant.

Biogeochemical processes can alter pore structures, making it important to understand solute transport behavior in evolving rocks (Hidalgo et al., 2015; Perez et al., 2022; Seigneur et al., 2019; Yang et al., 2024; Yoon et al., 2015; Zareei et al., 2022). Regarding dissolution-induced pore structure evolution, most prior research has primarily focused on how mineral dissolution alters pore structure, permeability, and flow patterns, with limited emphasis on solute transport (Hyman et al., 2024; Menke et al., 2016, 2023; Roded et al., 2020; Szawełło et al., 2024). For example, Roded et al. (2020) found that dissolution in the uniform regime tends to homogenize the medium structure, leading to a more uniform velocity distribution. Menke et al. (2016) showed that initial structure and velocity field heterogeneity significantly influence dissolution dynamics and identified a new dissolution regime termed the channeling regime. Szawełło et al. (2024) quantified flow focusing across different dissolution regimes and provided a quantitative measure for identifying wormholing, channeling, and uniform regimes. There are also many important studies that have examined how medium heterogeneity affects the dissolution process (Cheung & Rajaram, 2002; Hanna & Rajaram, 1998; Hyman et al., 2022; Maheshwari et al., 2013; Upadhyay et al., 2015; Z. Wang et al., 2022).

Despite these efforts, few studies have investigated how changes in pore structure due to mineral dissolution influence solute transport behavior. Recent studies have combined experiments and simplified models to explore how dissolution and/or precipitation alter transport properties (Edery et al., 2021; Muljadi et al., 2018; Saeibehrouzi et al., 2025; Shavelzon & Edery, 2024; Trabucchi et al., 2023). While these studies offer valuable insights, a comprehensive predictive framework for how dissolution governs solute transport is still lacking. In this study, we simulate mineral dissolution and solute transport in pore networks under various scenarios, elucidating how the interplay between medium heterogeneity and dissolution regimes—specifically uniform and wormholing—shapes solute transport behavior.

2 Methods

The 2D pore network represents the porous media as a series of interconnected cylindrical capillaries with a constant length l 0 ${l}_{0}$ , as illustrated in Figure 1a. The points where pores intersect are denoted as nodes, which are assumed to be volumeless and therefore non-reactive. In this study, the system consists of 100 × ${\times} $ 100 pores where the nodes on the left and right of the system serve as inlet and outlet nodes, respectively. During the dissolution process, the pore diameter increases proportionally to local reactant consumption (Budek & Szymczak, 2012; Sharma et al., 2023). Heterogeneity is introduced into the network by randomly assigning initial pore diameters from a lognormal distribution with mean diameter d 0 ${d}_{0}$ and standard deviation σ d ${\sigma }_{\mathrm{d}}$ . Here, we explored three levels of initial pore diameter heterogeneity ( σ $\sigma $ d = 0.001, 0.1, 1), whose underlying velocity fields lead to transport behaviors ranging from Fickian to non-Fickian. The diameter probability density function (PDF) and corresponding velocity distributions for the three initial networks are shown in Figure 1b.

Details are in the caption following the image
Figure 1
Open in figure viewerPowerPoint

(a) A schematic representation of a 2D pore network model with boundary conditions. The network consists of 100 × 100 interconnected cylindrical tubes, whose diameters increase proportionally to local reactant consumption. These tubes represent pores, and they intersect each other at nodes (black dots), which are assumed to be volumeless so that all reactions occur exclusively within the pores. Flow and reactant are injected from the inlet nodes (red dots) with a constant total flux Q and a constant inlet concentration cin. Green nodes indicate the side boundary nodes subject to periodic boundary conditions. (b) Heterogeneity is introduced into the network by randomly assigning initial pore diameters from a lognormal distribution. The main plot shows the Probability Density Function (PDF) of pore diameters for three different initial heterogeneity conditions with σ d ${\sigma }_{\mathrm{d}}$  = 0.001, 0.1, and 1. The inset presents the velocity PDFs that correspond to the initial pore diameter distributions. The error bars show the variability in values across 10 different realizations.

2.1 Dissolution Model

The pore-network model used in this work is based on the implementation of Budek and Szymczak (2012). The dissolution process is modeled by sequentially solving for flow, reactant concentration, and geometry evolution at each time step. These processes can be solved sequentially because, in typical geological systems, the time needed for dissolution to induce geometrical changes that affect the flow field is much longer than the flow and concentration relaxation times. For the boundary conditions, the reactive fluid with constant concentration, c = cin, is injected into the system through the inlet nodes, where a pressure p in ( t ) ${p}_{\text{in}}(t)$ is applied. The inlet pressure p in ( t ) ${p}_{\text{in}}(t)$ is adjusted at each time step to maintain a constant total volumetric flow rate Q through the system. The outlet has a constant pressure p out ( t ) = 0 ${p}_{\text{out}}(t)=0$ . Periodic boundary conditions are applied along the side boundaries of system (green nodes shown in Figure 1a). The fluid is assumed to be incompressible, and the flow field is assumed to be laminar and stationary at each dissolution time step. The flow field is then computed using the Hagen-Poiseuille equation in each pore and the continuity equation at each intersection:
q i j = − π d i j 4 128 μ l 0 p j − p i $\begin{array}{c}{q}_{ij}=-\frac{\pi {{d}_{ij}}^{4}}{128\mu {l}_{0}}\left({p}_{j}-{p}_{i}\right)\end{array}$ (1)
∑ i q i j = 0 $\begin{array}{c}\sum\limits _{i}{q}_{ij}=0\end{array}$ (2)
where q i j ${q}_{ij}$ is the volumetric flux of pore i j $ij$ , connecting to nodes i $i$ and j $j$ , p j − p i $\left({p}_{j}-{p}_{i}\right)$ is the pressure difference between nodes i $i$ and j $j$ , d i j ${d}_{ij}$ and l 0 ${l}_{0}$ are the diameter and length of the channel respectively, and μ $\mu $ is the dynamic viscosity of the fluid. Once the flow field is calculated, the concentration of the reactant along a pore i j $ij$ can be obtained from the mass balance equation by neglecting the axial diffusion:
q i j d c d x = − π d i j R $\begin{array}{c}{q}_{ij}\frac{dc}{dx}=-\pi {d}_{ij}R\end{array}$ (3)
where x is the coordinate along the axis of the pore. The reactant is consumed within the pores, with a rate R = k eff c $R={k}_{\text{eff}}c$ , assumed to be linear in concentration. Here, k eff ${k}_{\text{eff}}$ is the effective reaction rate constant, which takes into account the effect of transverse diffusion of reactant from bulk to the mineral surface. It is defined as:
k eff = k 1 + k d i j / D · Sh $\begin{array}{c}{k}_{\text{eff}}=\frac{k}{1+k{d}_{ij}/D\mathit{\cdot }\text{Sh}}\end{array}$ (4)

Here, k $k$ denotes the surface reaction rate constant, and D $D$ is the molecular diffusion coefficient of the reactant. Sh $\text{Sh}$ stands for the Sherwood number, which is a non-dimensional mass transfer coefficient (Bejan, 2013; Gupta & Balakotaiah, 2001). In principle, due to the buildup of the boundary layer, the Sherwood number can depend on the axial coordinate, particularly near the entrance of the pore (Ebadian & Dong, 1998). However, for the range of flow and reaction rates considered in this study, the entrance length is small relative to the pore length, allowing us to neglect entrance effects. Additionally, the value of the Sherwood number depends on the reaction rate itself, but the variation of Sh with k d / D $kd/D$ is relatively small (Gupta & Balakotaiah, 2001; Hayes & Kolaczkowski, 1994), bounded by two asymptotic limits: 5.39 for high reaction rates (transport limit) and 4.86 for low reaction rates (reaction limit) (Ebadian & Dong, 1998). This effect has negligible influence on the overall dissolution behavior observed in our simulations. In this work, we approximated Sh $\text{Sh}$ with a constant value of 5.

The effective reaction rate constant, k eff ${k}_{\text{eff}}$ , depends on the relative importance of reactive and diffusive effects. Correspondingly, in sufficiently narrow pores, diffusion is rapid enough to maintain an almost uniform concentration field across the pore diameter, making dissolution reaction-limited, with k eff ≈ k ${k}_{\text{eff}}\mathit{\approx }k$ . Conversely, in wider pores or in cases of fast reactions, the reaction rate becomes limited by the diffusive transport of the reactant across the aperture, resulting in k eff ≈ D Sh d ${k}_{\text{eff}}\mathit{\approx }\frac{D\text{Sh}}{d}$ . Assuming each pore has a uniform diameter, the concentration profile along a pore can be described by integrating Equation 3:
c ( x ) = c i j in e − π d i j k eff x q i j $\begin{array}{c}c(x)={c}_{ij}^{\text{in}}{e}^{-\tfrac{\pi {{d}_{ij}k}_{\text{eff}}x}{{q}_{ij}}}\end{array}$ (5)
Here, c i j in ${c}_{ij}^{\text{in}}$ denotes the concentration of reactant at the inlet of a pore. The concentration of reactant at the outlet of the pore can then be described in terms of the inlet concentration as c i j out = c i j in e − π d i j k eff l 0 q i j ${c}_{ij}^{\text{out}}={c}_{ij}^{\text{in}}{e}^{-\tfrac{\pi {{d}_{ij}k}_{\text{eff}}{l}_{0}}{{q}_{ij}}}$ . To calculate the concentration values at the nodes, we implement the complete mixing rule, which assumes that incoming reactant is instantaneously and homogeneously mixed before being distributed into the outgoing pores. The outgoing concentration at a given node i $i$ , which also serves as the inlet concentration for all downstream pores i j $ij$ connected to it, is determined by equation:
c i j in = ∑ k q k i c k i out ∑ k q k i , $\begin{array}{c}{c}_{ij}^{\text{in}}=\frac{\sum\limits _{k}{q}_{ki}{c}_{ki}^{\text{out}}}{\sum\limits _{k}{q}_{ki}},\end{array}$ (6)
where the sum is over all the incoming pores of node i $i$ and c k i out ${c}_{ki}^{\text{out}}$ is the outlet concentration of pore k i $ki$ . Based on the concentration profile, the erosion rate of a pore can then be calculated by:
∂ t d i j / 2 = k eff ν c sol c , $\begin{array}{c}{\mathit{\partial }}_{t}\left({d}_{ij}/2\right)=\frac{{k}_{\text{eff}}}{\nu {c}_{\text{sol}}}c,\end{array}$ (7)
where ν $\nu $ is the stoichiometric coefficient and c sol ${c}_{\text{sol}}$ is the molar concentration of soluble material. The total volume of dissolved mineral ∆ V $\mathit{{\increment}}V$ for a pore over time ∆ t $\mathit{{\increment}}t$ can therefore be calculated by integrating Equation 7 along the pore length. To keep the model tractable, we assume that the pores are uniformly broadening along their length. Given that the shortest penetration length in our model corresponds to about 10 pore lengths, the assumption of uniform pore broadening remains well-justified at the scale considered. Therefore, the enlargement of the pore diameter is given by:
∆ d = ∆ V π d i j l 0 $\begin{array}{c}\mathit{{\increment}}d=\frac{\mathit{{\increment}}V}{\pi {d}_{ij}{l}_{0}}\end{array}$ (8)
We rescale our variables based on d 0 ${d}_{0}$ , c 0 ${c}_{0}$ , and l 0 ${l}_{0}$ . Accordingly, the dimensionless time is defined as:
t ˆ = t τ d , $\begin{array}{c}\widehat{t}=\frac{t}{{\tau }_{d}},\end{array}$ (9)
where τ d ${\tau }_{d}$ is the characteristic dissolution time and is defined as:
τ d = 2 k γ d 0 , $\begin{array}{c}{\tau }_{d}=\frac{2k\gamma }{{d}_{0}},\end{array}$ (10)
where γ = c in / ν c sol $\gamma ={c}_{\text{in}}/\nu {c}_{\text{sol}}$ is the acid capacity number, describing the volume of solid dissolved by a unit volume of reactant. Here, τ d ${\tau }_{d}$ corresponds to the time needed for the pores at the inlet to increase their diameter by a factor of 2. We terminate dissolution simulations once the system's permeability increases by a factor of 1,000. Detail of the model implementation can be found in Budek & Szymczak, 2012; Sharma et al., 2023.

2.2 Passive Solute Transport Through Dissolving Pore Network

To examine the impact of dissolution on the evolution of transport behavior, we conduct solute transport simulations at multiple stages: at the onset (before dissolution), at an intermediate stage, and at the end (after the system's permeability has increased by a factor of 1,000). During these simulations, the pore network structure remains unchanged, assuming the transport timescale is significantly shorter than the dissolution timescale. Using the flow field derived from Equations 1 and 2, we conduct passive tracer transport simulations using particle tracking (P. K. Kang, Dentz, et al., 2015; P. K. Kang et al., 2017). For the transport simulation, 106 solute particles are introduced uniformly throughout the inlet nodes. We verified that transport statistics, such as breakthrough curves and dispersion coefficients, showed no significant change beyond 105 particles, confirming that our choice of 106 particles ensures statistical convergence. Each inlet node receives the same number of particles, and particles move from node to node with an average velocity in a given pore. The particle movement in space and time are described by the following equations:
x n + 1 = x n + l 0 u i j | u i j | $\begin{array}{c}{\boldsymbol{x}}_{n+1}={\boldsymbol{x}}_{n}+{l}_{0}\frac{{u}_{ij}}{\vert {u}_{ij}\vert }\end{array}$ (11)
τ n + 1 = τ n + l 0 | u i j | $\begin{array}{c}{\tau }_{n+1}={\tau }_{n}+\frac{{l}_{0}}{\vert {u}_{ij}\vert }\end{array}$ (12)
where x n ${\boldsymbol{x}}_{n}$ is the particle position after n transitions, and τ n ${\tau }_{n}$ is the total travel time after n transitions. u i j ${u}_{ij}$ is the average velocity through a pore ij that the particle is going through at the n + 1 $n+1$ step. The sign of u i j ${u}_{ij}$ is positive when the particle moves from node i $i$ to node j $j$ and negative when moving in the opposite direction, from j $j$ to i $i$ . The velocity u i j ${u}_{ij}$ is determined by:
u i j = 4 q i j π d i j 2 $\begin{array}{c}{u}_{ij}=\frac{4{q}_{ij}}{\pi {d}_{ij}^{2}}\end{array}$ (13)
At each node, particle chooses an outgoing pore following the full mixing rule. The transition probability of particle visiting pore ij is described using the following:
p i j = | q i j | ∑ k | q i k | $\begin{array}{c}{p}_{ij}=\frac{\vert {q}_{ij}\vert }{{\sum }_{k}\vert {q}_{ik}\vert }\end{array}$ (14)
where ∑ k | q i k | ${\sum }_{k}\vert {q}_{ik}\vert $ is the summation of flux over the outgoing pores of node i $i$ .

2.3 Dimensionless Numbers Characterizing the Dissolution Process

2.3.1 Effective Damköhler Number Da eff ${\text{Da}}_{\text{eff}}$ and Parameter G

The effective Damköhler number characterizes the relative importance of reaction kinetics to mass transport processes in a reactive system. It serves as a dimensionless number for analyzing and predicting dissolution patterns. In this work, the effective Damköhler number is defined as:
Da eff = π d 0 k L / q 0 1 + G $\begin{array}{c}{\text{Da}}_{\text{eff}}=\frac{\pi {d}_{0}kL/{q}_{0}}{1+\mathrm{G}}\end{array}$ (15)
where d 0 ${d}_{0}$ and L $\mathrm{L}$ are the initial mean diameter of the pore and the length of the system, respectively. Next, q 0 ${q}_{0}$ is a characteristic initial volumetric flow rate in a pore, q 0 = Q / 2 N inlet ${q}_{0}=Q/2{N}_{\text{inlet}}$ , where N inlet ${N}_{\text{inlet}}$ is the number of inlet nodes. The parameter G is defined as:
G = k d 0 D Sh $\begin{array}{c}\mathrm{G}=\frac{k{d}_{0}}{D\hspace*{.5em}\text{Sh}}\end{array}$ (16)
and quantifies the relative magnitude of reaction rate and the transverse diffusion rate across the pore, characterizing the hindering effect of the diffusive transport on the dissolution rate.

2.3.2 Penetration Length l p ${l}_{\mathrm{p}}$

Penetration length refers to the distance over which a reactive fluid significantly dissolves the mineral within a porous medium before the reaction rate decreases substantially. It is a measure of how far the reactant can penetrate into the medium. It is linked to the effective Damköhler number by
l p = L Da eff $\begin{array}{c}{l}_{\mathrm{p}}=\frac{L}{{\text{Da}}_{\text{eff}}}\end{array}$ (17)

A small Daeff corresponds to slow reaction relative to advection, resulting in a penetration length that can exceed the system size, and thus leads to uniform dissolution across the domain. Conversely, a large Daeff leads to shorter penetration lengths, where dissolution becomes concentrated along the most conductive flow paths and forms wormholes (Golfier et al., 2002; Hoefner & Fogler, 1988; Roded et al., 2020; Szymczak & Ladd, 2009).

3 Results and Discussions

To highlight the influence of dissolution regimes on the evolution of solute transport behavior, we first focus on two scenarios: (a) dissolution in a wormholing regime at Da eff = 10 and l p = 0.1 L ${\text{Da}}_{\text{eff}}=10\,\text{and}\,{l}_{\mathrm{p}}=0.1L$ , with an initial network diameter heterogeneity of σ d = 0.001 ${\sigma }_{\mathrm{d}}=0.001$ , corresponding to a conductivity heterogeneity of σ ln K = 0.04 ${\sigma }_{\ln\,K}=0.04$ ; (b) dissolution in a uniform regime at Da eff = 0.01 and l p = 100 L ${\text{Da}}_{\text{eff}}=0.01\,\text{and}\,{l}_{\mathrm{p}}=100L$ with an initial diameter heterogeneity of σ d = 0.1 ${\sigma }_{\mathrm{d}}=0.1$ , corresponding to a conductivity heterogeneity of σ ln K = 3.3 ${\sigma }_{\ln\,K}=3.3$ . In other words, we start with a homogeneous network for the wormholing regime and a heterogeneous network for the uniform regime. After analyzing these two cases, we expand the study to explore a broader parameter space.

3.1 Dissolution Induces the Transition of Transport Regimes

Figure 2 presents the probability density functions (PDFs) of solute breakthrough times and spatial maps of particle locations obtained by the particle tracking for the two scenarios representing wormholing and uniform dissolution regimes. Because the shape and temporal evolution are identical to breakthrough curves (BTCs), we use the terms BTC and PDF interchangeably in the following text. BTCs are obtained by ensemble-averaging over 10 independent realizations. Comparison of results from 20 realizations confirmed that 10 realizations are sufficient. The tracer BTCs from initial, intermediate, and final stages of the dissolution simulation show significant alteration in transport regimes due to dissolution. In the wormholing regime, dissolution causes a transition from Fickian to non-Fickian transport. As seen in Figure 2a, the BTC changes from a distribution with limited spreading (indicative of Fickian transport) to a distribution featuring early arrival and late-time tailing. The final BTC shows a unique feature of a plateau followed by a power-law slope of approximately 1.7. Note that non-Fickian transport is characterized by a power-law tail with a slope less than 3 (Dentz et al., 2004; Margolin & Berkowitz, 2000). The inset of Figure 2a shows an increase in the variance (second moment) of solute breakthrough times as a function of dissolution time, indicating the substantial enhancement of solute spreading due to dissolution. This finding is further supported by additional analyses of the BTC tail slope, skewness, and log-variance of particle arrival times (Figures S1 and S2 in Supporting Information S1). Figures 2b and 2c are snapshots of particle locations at three different times, illustrating particle transport dynamics through the network at the initial and final stages of dissolution, respectively, with the diameter field shown in the background. In the initial, undissolved network (Figure 2b), particle movement occurs with limited longitudinal spreading due to the lack of heterogeneity. However, at the final stage of dissolution (Figure 2c), longitudinal spreading is significantly enhanced. Particles primarily travel through a major wormhole (green), leading to early breakthrough, followed by transport in secondary wormholes (yellow), while particles in unreacted regions (blue) get stagnant and move slowly (Movie S1).

Details are in the caption following the image
Figure 2
Open in figure viewerPowerPoint

Transitions between Fickian and non-Fickian transport during mineral dissolution in the wormholing and uniform regimes. (a, d) Probability density functions (PDFs) of solute breakthrough times—equivalent to breakthrough curves (BTCs)—at different stages of dissolution. (b, c, e, f) Spatial maps of particle transport through the pore network at the initial and final dissolution stages. Particle locations at three representative breakthrough times (τ) are shown in green, yellow, and blue. (a–c) Wormholing regime: The BTC evolves from a narrow distribution (black) to one with early arrivals and a power-law tail (slope ≈ 1.7, red), indicating a transition from Fickian to non-Fickian transport. Initially (b), particles exhibit limited longitudinal spreading due to low heterogeneity, consistent with Fickian behavior. At the final stage (c), pronounced spreading emerges as particles primarily travel through a dominant wormhole (green), with later arrival via secondary channels (yellow) and stagnation in unreacted zones (blue), characteristic of non-Fickian transport. (d–f) Uniform dissolution regime: The BTC transitions from a broad distribution with power law tailing (black, slope ≈ 3) to a narrower one without significant tailing (red), indicating a shift from non-Fickian to Fickian transport. Initially (e), strong longitudinal spreading arises from structural heterogeneity. By the final stage (f), uniform dissolution homogenizes the network, reducing spreading and restoring Fickian behavior. Insets: Evolution of the variance of normalized breakthrough times ( σ τ / τ peak 2 ${\sigma }_{\tau /{\tau }_{\text{peak}}}^{2}$ ) during dissolution. Colored circles (black, green, red) indicate the initial, intermediate, and final stages. All BTCs are normalized by the peak arrival time, and breakthrough times in spatial maps are likewise normalized. In (c), the three selected times are marked on the final BTC in (a).

The plateau observed in the BTC is a strong indication of anomalous transport and it corresponds to particles that enter the major wormhole from neighboring locations as well as those moving through the secondary wormholes (refer to Movie S1 and Figure S3 in Supporting Information S1). The arrival of particles through these secondary wormholes leads to a secondary peak in the BTC, appearing as a plateau after ensemble averaging of 10 realizations. This secondary peak is more evident when examining a single realization (see the final BTC in Movie S1). The double-peaked pattern in solute BTCs has been widely reported in field and Darcy-scale numerical studies (Dewaide et al., 2017; Edery et al., 2021; Guihéneuf et al., 2017; Mathias et al., 2007). Various mechanisms have been proposed to explain this phenomenon. In fractured systems, contributing factors include the recirculation zones, which trap solutes and prolong residence time (L. Wang et al., 2020) and the presence of multiple flow pathways (Mathias et al., 2007). In karst systems, additional factors such as underground lakes, waterfalls, and solutes being diverted into auxiliary conduits due to blockages in the main channel have been identified as causes of double- or multi-peaked BTCs (Dewaide et al., 2017; Field & Leij, 2012). Our study suggests that the secondary dissolution channels can also contribute to the double-peaked BTC behavior.

In the uniform dissolution regime, dissolution drives the transition from non-Fickian to Fickian transport. Figure 2d illustrates a dramatic transition in BTC behavior: from a distribution exhibiting both early arrivals and late-time tailing, characterized by a power-law tail with a slope of approximately 3, to a distribution with limited spreading and without significant tailing. The inset of Figure 2d shows a decrease in the variance of solute breakthrough time during dissolution (see also additional analyses of particle arrival times in Figures S1 and S2 in Supporting Information S1). Particle spreading in Figures 2e and 2f show a significant reduction in longitudinal spreading due to the homogenization of the network through uniform dissolution (For further visualization, please refer to Movie S2).

By comparing Figures 2c and 2e, along with the shapes of the final BTC in Figure 2a and the initial BTC in Figure 2d, we observe that the non-Fickian behavior induced by dissolution is fundamentally different from that in a heterogeneous network with random heterogeneity. Specifically, the plateau regime (or the secondary peak) in the final BTC of the wormholing dissolution is not observed in the BTC of the initially heterogeneous network with randomly assigned pore sizes. This suggests that different upscaling strategies may be necessary to effectively capture anomalous transport in heterogeneous networks depending on the origin of heterogeneity. The non-Fickian behavior observed in the initially heterogeneous system in this study is similar to transport in heterogeneous porous rocks, soils, or fracture networks without significant secondary alteration. Anomalous transport in such systems is often captured by traditional stochastic models, which assume independent and identically distributed probability density functions (Dentz et al., 2004). In contrast, the non-Fickian behavior caused by dissolution may resemble transport in systems with spatially correlated structures (e.g., wormholes), often leading to double-peaked BTCs, such as in karst systems. Models considering velocity correlation will be an appropriate model for such cases (de Anna et al., 2013; Davy et al., 2024; P. K. Kang et al., 2017; P. K. Kang, Le Borgne et al., 2015; Morales et al., 2017). This observation indicates that the nature of subsurface heterogeneity—whether induced by dissolution or pre-existing—can potentially be inferred from BTC characteristics.

3.2 Mechanisms Leading to the Transition of Transport Regimes by Dissolution

To elucidate the mechanisms driving the transition between Fickian and non-Fickian transport during dissolution process, we now analyze the flow field and dissolution patterns in detail. Dissolution in both uniform and wormholing regimes leads to distinct dissolution patterns and heterogeneity alterations. Figures 3 and 4 illustrate the alterations in pore diameter and flow fields caused by the dissolution in wormholing and uniform regimes, respectively.

Details are in the caption following the image
Figure 3
Open in figure viewerPowerPoint

Dissolution patterns and heterogeneity alteration caused by dissolution in wormholing regime. (a) Diameter field at the initial, intermediate, and final stage of dissolution. (b) The corresponding flow rate fields. The thickness of the lines in the diameter and flow rate maps is proportional to the magnitude of diameter and flow rate, respectively. (c) Velocity distribution (PDF) at the initial and final stages of dissolution, showing the increase in the velocity heterogeneity due to the wormholing. Inset: evolution of velocity variance during the dissolution process. The three circles with black, green, and red color correspond to the initial, intermediate, and final stages of dissolution.

Details are in the caption following the image
Figure 4
Open in figure viewerPowerPoint

Alterations in dissolution patterns and heterogeneity caused by dissolution in uniform regime. (a) Diameter field at the initial and final stage of dissolution. (b) The corresponding flow rate field showing flow field homogenization. (c) The velocity distribution showing the decrease in the velocity heterogeneity due to the uniform dissolution. The inset figure is the diameter distribution (PDF) of initial network and network after dissolution with G = 0 and G = 100, demonstrating the difference between two homogenization mechanisms. At G = 0, without the effect of transverse diffusion, the network is homogenized because dissolution decreases the pore diameter ratio. The diameter PDF after dissolution is shifted without changing its shape, with the variance remaining mostly unchanged from 1.05e-2 to 1.08e-2. At G = 100, the network is further homogenized because the transverse diffusion limits the growth of wider pores relative to the thinner ones. The shape of diameter PDF after dissolution becomes narrower, with variance reduced to 1.09e-3.

The wormholes form due to reactive-infiltration instability, based on a positive feedback loop between flow, transport, and dissolution: a small perturbation to the dissolution front results in a locally increased flow rate, which, in turn, leads to increased dissolution, amplifying the perturbation. As a result, the flow and reactant transport become focused in highly conductive channels that advance quickly into the porous matrix. As reactive fluid preferentially flows through these wormholes, other areas receive significantly less reactant, resulting in limited dissolution outside the wormholes (Hoefner & Fogler, 1988; Szymczak & Ladd, 2009). These undissolved regions eventually become stagnant zones trapping particles. Figure 3a shows that dissolution changes the network pattern from a homogeneous to a highly heterogeneous network with one leading wormhole and multiple secondary wormholes near the inlet of the system. The resulting flow field, as shown in Figure 3b, also evolves from being uniform to preferential flows through wormholes. The velocity PDFs in Figure 3c confirm that wormholing increases flow field heterogeneity, leading to the observed transition from Fickian to non-Fickian transport. The inset in Figure 3c shows that the variance of the velocity field increases dramatically as dissolution progresses.

The formation of wormholes significantly alters the flow field, leading to distinct transport behaviors at different dissolution stages. At the intermediate stage, wormholes extend halfway through the system, leaving a large undissolved region ahead of their tips. As a result, flow is concentrated within the wormholes before redistributing uniformly beyond the tips (see Figure 3b and Movie S1 at the intermediate stage). This flow focusing leads to distinct particle arrival times within the wormholes, producing multiple peaks in the BTC. The similar findings have been reported using entropy-based analyses (Shavelzon & Edery, 2024; Zehe et al., 2021). After ensemble averaging, these peaks appear as plateau, which resembles later part of the BTC in the final stage (see the intermediate and final BTCs in Figure 2a and Movie S1). Unreacted regions act as trapping zones, governing the late-time transport behavior.

At the final stage of dissolution, when a leading wormhole breaks through the system, the flow field is divided into three distinct transport regions (as shown in Figure 2c and Movie S1): (a) a highly permeable preferential flow channel that carries the majority of the flow, (b) secondary flow channels with relatively lower permeability, and (c) low-velocity regions with minimal dissolution. This division corresponds to the three different tracer breakthrough behaviors observed in Figures 2a and 2c, characterized by early arrival followed by a plateau region, and then power-law tailing. Early arrival occurs due to particles transported through the major wormhole, while the plateau results from particles near the major wormhole that preferentially enter it, as well as those transported through secondary wormholes. Particles that enter the low-velocity regions take the longest time to break through, leading to power-law tailing. Thus, the shape of BTC shown in Figure 2a can be understood as the combined contribution of these three distinct transport zones (Figure S3 in Supporting Information S1 further confirms this relationship). Additionally, the change in BTC shape from the intermediate to the final stage is mainly driven by the breakthrough of the leading wormhole, which induces strong flow focusing and leads to the emergence of early arrival and an extended plateau in the BTC. As shown in Figures S3 and S4 in Supporting Information S1, the transport behaviors in Figures S3b and S3c in Supporting Information S1 are governed by the leading wormhole, whereas those in Figures S3d and S3e in Supporting Information S1 correspond to transport patterns that also exist in the intermediate stage (Figures S4b and S4c in Supporting Information S1).

In the uniform regime, dissolution homogenizes the network, in contrast to flow focusing observed in the wormholing regime. At low Daeff (large lp/L), the flow rate is significantly higher than the reaction rate, allowing the reactant to penetrate deep into the network before noticeable dissolution occurs. Since both large and small pores experience similar reactant concentrations, all pores dissolve at comparable rates and grow by a similar magnitude. Consequently, as dissolution progresses, the ratio of pore diameters decreases, leading to a reduction in the log variance of the pore size distribution (see Figure S5 in Supporting Information S1, G = 0 case, in the Supporting Information S1). Additionally, the hindering effect of diffusion further accelerates the homogenization of the pore space. In the case shown in Figures 2d–2f, the G parameter is 1, indicating that transverse diffusion and surface reaction rate have comparable influences. For larger pores, reactant diffusion from the bulk to the pore walls becomes less efficient. This effect is captured in Equation 4, where a larger pore diameter d $d$ corresponds to a lower effective reaction rate. As a result, smaller pores dissolve faster than larger ones, which enhances the homogenization of the pore structure over time (Roded et al., 2020).

Figures 4a and 4b show that, although initial diameter heterogeneity creates multiple preferential flow paths, both the pore diameters and flow rates become more uniform by the final stage of dissolution. Accordingly, the velocity PDF becomes significantly narrower as dissolution progresses (Figures 4c and in Supporting Information S1), indicating that uniform dissolution reduces velocity heterogeneity and drives the transition from non-Fickian to Fickian transport behavior (refer to Figure 2d for BTCs). To distinguish between the two homogenization mechanisms, we conducted additional dissolution simulations under two conditions: without transverse diffusion effects (G = 0) and with strong hindering effects from transverse diffusion (G = 100). Both the G = 0 and G = 100 cases exhibit homogenization effects, with G = 100 displaying significantly stronger homogenization due to the hindering effect. As shown in the inset of Figure 4c, for the G = 0 case, where all pores dissolve at a similar rate, the diameter distribution shifts toward larger diameters without changing its shape. In contrast, for the G = 100 case, the effect of transverse diffusion hindrance leads to a significant narrowing of the diameter distribution. Figures S5 and S7 in Supporting Information S1 further demonstrate the homogenization effect. The homogenization is more evident in a log scale (Figure S5 in Supporting Information S1), with the log variance decreasing from 7.02e-1 to 7.97e-3 for G = 0 and from 7.02e-1 to 7.48e-4 for G = 100.

3.3 Impact of Initial Heterogeneity and Damköhler Number

We extend our findings by conducting simulations with varying degrees of initial heterogeneity at the two different Da values. Based on the results, we constructed a phase diagram that characterizes solute transport regime transitions (Figure 5) by comparing BTC behavior of the initial network and the network at the final stage of dissolution. Figure 5 shows that, regardless of the initial heterogeneity, wormholing leads to non-Fickian transport, while uniform dissolution results in Fickian transport. Depending on the combination of initial heterogeneity and Damköhler number, four types of transport behavior are observed: persistent Fickian, non-Fickian to Fickian, Fickian to non-Fickian, and persistent non-Fickian. As demonstrated in the previous section, dissolution in the wormholing regime increases heterogeneity by forming dissolution channels, leading to non-Fickian transport. For an initially homogeneous network, where transport is Fickian, wormholing dissolution shifts the behavior from Fickian to non-Fickian, as shown in Figure 5c. For an initially heterogeneous network with non-Fickian transport, the behavior remains non-Fickian after dissolution (Figures 5a and 5b), but the BTC shape changes and the power-law tailing intensifies. In the uniform regime, dissolution reduces flow heterogeneity, inducing Fickian transport. For an initially heterogeneous network with non-Fickian transport, dissolution transitions the behavior to Fickian (Figures 5d and 5e). Meanwhile, for an initially homogeneous network where transport is Fickian, the behavior remains Fickian after dissolution (Figure 5f) in a uniform dissolution regime.

Details are in the caption following the image
Figure 5
Open in figure viewerPowerPoint

Phase diagram of BTCs showing the impact of initial heterogeneity and Damköhler number on transitions in solute transport regimes. (a–c) In the wormholing regime (Daeff = 10), dissolution leads to non-Fickian transport. Final BTCs (red) for all three levels of initial heterogeneity (σd = 1, 0.1, and 0.001) exhibit early arrival and late-time tailing. However, the extent of the plateau region and the power-law slope of the tail vary with initial heterogeneity. (d–f) In the uniform regime (Daeff = 0.01), dissolution leads to Fickian transport. Final BTCs (red) become narrow and lack power-law tails, though greater initial heterogeneity hinders complete homogenization. Transitions in transport behavior are governed by the interplay between initial network heterogeneity and Damköhler number and there are four types of transport behavior based on the combination of initial heterogeneity and Damköhler number: (a, b) persistent non-Fickian, (c) Fickian to non-Fickian, (d, e) non-Fickian to Fickian, and (f) persistent Fickian.

The impact of initial heterogeneity on the dissolution pattern and shape of BTC are different between uniform and wormholing regime. For uniform dissolution, the BTCs at the final stage for all three levels of initial heterogeneity become narrow and lack power-law tailing, as demonstrated in Figures 5d–5f. Additionally, the diameter heterogeneity decreases throughout the dissolution, but achieving full homogenization becomes more challenging with increasing initial heterogeneity. For the σ d = 1 ${\sigma }_{\mathrm{d}}=1$ cases, even after a 1,000-fold increase in permeability, the level of heterogeneity remains significantly higher than in the σ d = 0.1 ${\sigma }_{\mathrm{d}}=0.1$ and 0.001 $0.001$ cases. In the wormholing regime, initial heterogeneity affects the dissolution pattern and, consequently, controls the shape of BTC. As shown in Figures S8 and S9 in Supporting Information S1, dissolution leads to the formation of a dominant wormhole across all three levels of initial heterogeneity. However, the shape of the leading wormhole, as well as the structure and number of secondary wormholes, vary depending on initial heterogeneity. As previously discussed, the secondary wormholes and underlying velocity heterogeneity are the primary contributors to the plateau regimes and late-time tailing of the BTC. As shown in Figures 5a–5c, the extent of the plateau regime and the power-law slope of late-time tailing are sensitive to initial heterogeneity, underscoring its influence on transport behavior.

4 Conclusions

In this study, we explored the impact of dissolution on solute transport behavior by examining various combinations of initial heterogeneity, Damköhler number, and G parameter. Through simulations covering a wide range of scenarios, we observed significant changes in flow patterns and transport behavior due to dissolution. These changes were strongly influenced by the initial heterogeneity and the dissolution regime. In the wormholing regime, where the penetration length is shorter than the system length, dissolution led to the formation of preferential flow paths, strongly enhancing network heterogeneity. This increase in heterogeneity caused a transition from Fickian to non-Fickian transport behavior, characterized by early solute arrival and long-tailed BTCs. Unique to wormholing regime, we observed distinctive non-Fickian transport behaviors, such as plateaus and double peaks in BTCs, which have not been reported in networks with random heterogeneity. These phenomena were explained by the interplay of main wormholes, secondary wormholes, and stagnation zones between the wormholes.

Conversely, in the uniform dissolution regime, where the penetration length is longer than the system length, dissolution homogenized the network by equalizing pore sizes over time. The homogenization was driven by two mechanisms. First, under reaction-limited conditions, pores grow at similar rates, causing the ratio of pore sizes to decrease with time. Second, under transport-limited conditions, transverse diffusion hindrance limits the growth of larger pores, while smaller pores grow more rapidly, further decreasing pore size variability. Such homogenization mechanisms reduce network heterogeneity, leading to a more uniform flow field, and induce a transition from non-Fickian to Fickian transport behavior, as evidenced by the narrowing of the BTCs. A phase diagram created by varying initial heterogeneity and Damköhler number revealed four possible types of transport behavior: persistent Fickian, non-Fickian to Fickian, Fickian to non-Fickian, and persistent non-Fickian. These results demonstrate that dissolution not only alters the pore structure but also significantly impacts solute transport behavior, highlighting the critical role of initial conditions and the reaction-flow interplay in determining transport regimes. The findings suggest that the evolution of transport dynamics could be predicted based on the initial heterogeneity and Damköhler number.

This study provides new insights into transport behavior in evolving porous media, where dissolution is prevalent. The non-Fickian transport behavior caused by wormholing highlights how wormholes create highly localized transport pathways, leading to distinctive BTC patterns. These pathways fundamentally alter solute migration, with significant implications for predicting contaminant dispersion and managing resources in such systems. For example, such unique transport behavior suggests a practical application: monitoring tracer transport in karst systems or fractured aquifers could serve as an early warning mechanism for wormholing risks, particularly beneath hydraulic structures such as dams in gypsum terrains. Furthermore, the distinctive BTC shapes identified in this study indicate that key subsurface processes and heterogeneity can manifest in BTCs. Specifically, the extent of dissolution and the prevailing dissolution regime could be inferred by analyzing BTC patterns. Although this study is based on a 2D pore network, both wormholing and uniform dissolution have commonly observed in 3D systems (Fredd & Fogler, 1998; Hoefner & Fogler, 1988). However, dimensionality can influence flow focusing and wormhole morphology (Cohen et al., 2008), and the full implications of 3D transport effects remain to be explored.

Acknowledgments

Research supported as part of the Center on Geo-process in Mineral Carbon Storage, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences at the University of Minnesota under award #DE-SC0023429. Peter K. Kang acknowledges the National Science Foundation under Grant EAR-2046015 and the STEAM Program (RS-2024-00461440) from the National Research Foundation of Korea for partial support of this research. Jingxuan Deng acknowledges partial support for this research provided by a Data Science Graduate Assistantship with funding made available by the MnDRIVE initiative through the University of Minnesota Data Science Initiative. Piotr Szymczak acknowledges the support of the National Science Centre (Poland) under research Grant 2022/47/B/ST3/03395. He also acknowledges the support through the Fulbright STEM Award program for funding his visit at the University of Minnesota.

    Data Availability Statement

    All simulation data files are available through the Data Repository for University of Minnesota (DRUM) (Deng et al., 2025), https://conservancy.umn.edu/items/a78964c3-e6c3-418c-a308-e85371aabb4d.

      相关资源