Linking Self-Organized Heterogeneities to Solute Transport in Mixing-Induced Precipitated Porous Media

1. Introduction

Mineral precipitation is a fundamental geochemical process in porous media that can significantly alter the hydraulic and transport properties of aquifers [1,2,3]. This phenomenon typically occurs when waters with different chemical signatures mix, creating supersaturation conditions that promote the formation of new minerals in the porous matrix [4,5]. Over time, sustained precipitation can progressively clog the pores of the medium, reducing the effective porosity and permeability of the system, which, in turn, influences groundwater flow and solute transport [6,7,8]. Among the various types of mineral precipitation, carbonate precipitation has gained increasing attention due to its significant role in subsurface engineering applications [9,10,11], particularly in environments where calcium- and carbonate-rich waters interact [12]. Recent advancements have focused on controlled precipitation techniques, such as microbially induced calcium carbonate precipitation (MICP) [13,14,15] and enzymatic methods [16,17,18]. These techniques enable targeted mineral formation in specific subsurface environments, offering precise control over mineral deposition and its effects on permeability [19]. These precipitation techniques have found widespread applications across various disciplines [20]. For instance, in geological carbon sequestration (GCS), carbonate precipitation is intentionally induced to seal preferential flow paths in low-permeability cap rocks, reducing the risk of CO₂ leakage [21,22,23]. Similarly, in enhanced oil recovery (EOR), MICP has proven successful in modifying rock permeability by plugging high-permeability flow channels, thereby enhancing sweep efficiency and optimizing hydrocarbon extraction [24,25,26]. However, the efficiency of these techniques is strongly influenced by the spatial heterogeneity that mineral precipitation introduces within the porous medium [5,27,28]. This heterogeneity [29,30], significantly affects the spatial variability in hydraulic conductivity, which in turn impacts solute transport. Therefore, understanding the relationship between the heterogeneous structure of the porous medium formed through mixing-induced precipitation and its impact on solute transport is crucial for improving the use of precipitation-based techniques in environmental and engineering applications.

Prior research has identified two distinct categories of mineral precipitation mode, preferential precipitation, and uniform precipitation [31,32,33]. In the preferential precipitation mode, the solid phase tends to form scattered crystals on the rock-fluid interface [34,35], forming a heterogeneous precipitation texture that evolve across multiple spatial and temporal scales [5,27], whereas the uniform precipitation mode results in a uniform coating of precipitate [28]. Three predominant factors influence these mineral precipitation modes: solute mixing [5], mineral nucleation [28], and the relative timescale of reaction kinetics compared to ion transport [36].

Solute mixing plays a critical role in influencing the spatiotemporal evolution of reaction products [37,38,39]. This mixing affects the spatial distribution and concentration of dissolved ions, which, in turn, influences the precipitation process. Spatiotemporal patterns in the development of CaCO₃ deposits have been reported, including small isolated deposits on the surface of the porous medium grains, elongated structures following the main flow direction, and homogeneous precipitates along the flow front. Mineral nucleation, however, remains a complex phenomenon involving various physical processes such as fluid flow, species transport, interfacial processes, and geochemical reactions [33,40,41]. The nucleation rate plays a crucial role in determining the density and distribution of newly formed crystals [28]. A high nucleation rate leads to a more uniform distribution of the solid phase, while a lower nucleation rate results in fewer, larger crystals that grow rapidly, promoting mineral heterogeneity and causing a more pronounced reduction in hydraulic conductivity. Alongside nucleation, the relative timescale of reaction kinetics compared to ion transport is vital in shaping mineral precipitation [36,42,43]. When reaction rates are slow relative to transport, dissolved ion concentrations remain relatively uniform, fostering consistent mineral growth. However, when reaction kinetics surpass ion replenishment, concentration gradients develop, leading to localized precipitation.

The spatial variability of mineral precipitation, resulting from these factors, drives non-uniform changes in local reaction rates that modify the internal structure of the porous medium [29,30]. This, in turn, increases permeability variations at multiple scales, promoting the formation of preferential flow paths [44,45]. Experimental studies focused on sealing fractured rocks have demonstrated that mineral precipitation generates self-organizing flow channels within the fracture network [34,46]. As these areas fill, the flow is progressively redirected to regions where no precipitation occurs, forming new flow paths [34], with tortuous and braided pathways emerging within the fracture [46]. The heterogeneity caused by preferential precipitation in porous media has profound implications for solute transport. The irregular distribution of precipitates and channeling significantly modifies flow velocities and dispersive properties, leading to deviations from classical Fickian behavior. This non-Fickian transport manifests as anomalous breakthrough curves, which traditional models fail to accurately capture [47,48]. The challenge in modeling these systems lies in the difficulty of characterizing the spatial variability of permeability caused by freshly formed precipitates. Successfully addressing this challenge requires determining the appropriate level of complexity needed to accurately represent anomalous transport and identifying the key features that must be incorporated into numerical models.

This study investigates the relationship between precipitation-induced heterogeneities and solute transport by analyzing dominant precipitation patterns and their impact on breakthrough concentrations. Conservative tracer tests were conducted on an experimental aquifer before and after mixing-induced calcite precipitation to assess how heterogeneities at different scales evolve with time and impact solute transport behavior. These experiments build on previous experimental findings [48], which identified unusual breakthrough curves with dual peaks and extended tails, directly linked to precipitation-induced structural heterogeneity. To address these phenomena, this study presents the interpretation of the experiment through a solute transport model that simulates anomalous transport behavior resulting from mineral precipitation using a novel method to characterize the spatial variability of permeability determine from image analysis. In this context, this work aims to: (i) evaluate how to characterize the heterogeneity in order to reproduce the anomalous solute transport behavior after mineral precipitation, (ii) understand the level of complexity required to accurately reproduce solute transport in heterogeneous systems, and (iii) assess the key features of heterogeneity that must be represented in numerical models to correctly reproduce solute transport. To achieve these goals, this study introduces a numerical model that simulates solute transport in porous media affected by mineral precipitation. The model incorporates detailed representations of the evolving spatial heterogeneity resulting from the precipitation process. By leveraging laboratory-derived data on permeability fields influenced by calcite precipitation, this research offers a unique perspective on how varying levels of heterogeneity—driven by both the amount and distribution of precipitate affect solute transport. The methodology includes a comparative analysis between homogeneous and heterogeneous models, where the latter is defined by multiple permeability zones based on varying image intensities of the precipitation. This approach allows for a deeper understanding of the relationship between the distribution of the precipitated mineral and its influence on solute transport, while also highlighting the critical role of spatial variability in accurately predicting solute transport. The novelty of this approach lies in the ability to replicate the anomalous transport behavior observed in real-world systems that are affected by transient precipitation processes. Furthermore, this study provides valuable insights into the implications of incorporating heterogeneous permeability fields induced by precipitation into transport models, which can be applied to both field studies and laboratory experiments. In particular, the findings have direct relevance for applications in environmental engineering, such as carbon sequestration enhanced oil recovery (EOR), where mineral precipitation is often used to enhance solute retention and mobility.

The structure of this article is as follows. The Experimental Overview section summarizes previous experiments conducted on Mixing-Driven Precipitation (MDP) and tracer tests, including the setup, procedure, and main findings. The Methods section describes the methodology implemented to create the composite permeability fields, presents the flow and transport model, and explains the model calibration. The Results and Discussion section presents the main findings obtained from the dual- and triple-permeability composite field models. Finally, the Conclusions section summarizes the key insights, emphasizing the importance of accurately characterizing the heterogeneity caused by mineral precipitation to effectively represent the resulting anomalous solute transport.

2. Experimental Overview

A laboratory experiment was conducted by Guido González-Subiabre and Fernàndez-Garcia [48], to evaluate the impact of mineral precipitation induced by mixing on solute transport, consisting of three main stages (Figure 1a). First, an initial conservative tracer test (T₁) was performed to establish the systems initial transport conditions through a breakthrough curve (BTC). This was followed by the Mixing-Driven precipitation (MDP) experiment, during which synthetic solutions of $0.05\mathrm{mol}/\mathrm{kgw}$ CaCl₂ ($W_1$) and $0.1\mathrm{mol}/\mathrm{kgw}$ Na₂CO₃ ($W_2$) were simultaneously and uniformly injected in parallel, each through half of the inlet face of the tank, to form the calcite layer. This configuration induces a mixing zone that extends along the centerline of the tank as shown in Figure 1b. The solutions were selected based on previous experimental research [7,49,50]. Finally, a second conservative tracer test (T₂) was conducted to capture changes in the BTC resulting from the freshly formed calcite after the MDP stage. A detailed description of the experimental procedure, the setup and the results is provided by Guido González-Subiabre and Fernàndez-Garcia [48]. A two-dimensional horizontal tank made of plexiglass, designed to replicate a coarse sand granular confined aquifer, was used to perform the Mixing-Driven Calcite Precipitation (MDP) experiment and companion conservative tracer tests (Figure 1b). The tank was homogeneously wet-packed with spherical glass beads (d = 2 mm) and initially filled with deionized water (Milli-Q). The synthetic aquifer was designed to provide uniform flow conditions, featuring eight inlet and outlet ports spaced 2 cm apart. To monitor the evolution of the freshly precipitated calcite, image analysis was employed. The experimental setup was placed in a darkroom to minimize external light interference, while a 1550 Lm, 20 W Downlight LED positioned above the tank served as the light source for visualization. The intensity of the reflected light was then used to track the development of the precipitate over time. Images were captured every 30 seconds using a Nikon D7100 camera paired with a Tamron SP AF17-50mm F/2.8 XR Di II LD Aspherical (IF) Model A16 lens. For the tracer tests, a fluorescein solution with a uniform inlet concentration of c₀ = 3 mg/L was used. The concentrations at the outlet were measured with an Albillia FL24 fluorometer, which provided an integrated breakthrough curve for the entire system. Continuous monitoring with the fluorometer was conducted, taking measurements every second. Precipitation and tracer test experiments were conducted under advection-dominated transport conditions, supported by a grain Péclet number of Pe = 523. A prescribed flow boundary was established at the system entry with a total inflow rate of Q= 1.78 $\times {10}^{-7}{\mathrm{m}}^3/\mathrm{s}$, evenly distributed using 8 inlet ports. A constant head boundary condition was imposed (0.9 m above the tank elevation) at the tank outlet. The initial hydraulic conductivity ($K_0$) and porosity ($\varphi$) of the system were both experimentally estimated, with $K_0$ = 145.4 m/day and $\varphi$ = 0.34. For the experimental properties of the synthetic aquifer used in the experiment, consult Table 1. The results obtained in the experiment are shown in Figure 1(c). This finding indicated that after the MDP experiment, a symmetric, bell-shaped precipitation front was formed with a self-organized, heterogeneous pattern. This pattern was characterized by internal elongated carbonate structures that created preferential flow paths. Breakthrough curves revealed a shift in solute transport behavior from a classical Fickian-like distribution to a non-Fickian anomalous distribution, exhibiting a double-peak curve with pronounced tailing. Although the origin of this phenomenon is not fully understood, it is hypothesized to involve channeling and back-diffusion in low-permeability zones. Guido González-Subiabre and Fernàndez-Garcia [48] observed that precipitation-induced heterogeneity affects solute transport, yet the underlying mechanisms remain unclear. Simulating such heterogeneity in transport models is challenging due to the difficulty in capturing the evolving pore structure and permeability across multiple scales. This article provides insights to effectively represent mineral-induced heterogeneity in transport models. To understand the observed anomalies, we focus on modeling the heterogeneity caused by mineral precipitation. This is accomplished by processing the steady-state calcite precipitation image using a novel method to generate permeability fields that reflect the complexity of the induced heterogeneity. Different degrees of complexity are examined. The image processing technique, based on Schuszter et al. [51], assumes that the precipitate amount is proportional to the log of pixel intensity [51,52,53]. The generated permeability fields were designed to represent both the total area covered by the precipitate and regions with the highest concentration of precipitation, following a progressive approach from simple to complex representations. The breakthrough curve (BTC) obtained before mineral precipitation was used to determine the initial dispersivities (${\alpha }_L$ and ${\alpha }_T$), using a homogeneous permeability model characterized by a single hydraulic conductivity value ($K_0$ = 145.4 m/day), representing the initial hydraulic conductivity measured for the system. For the BTC obtained after mineral precipitation, two different permeability distribution scenarios were considered: a dual-zone and a triple-zone composite permeability field.

3. Materials and Methods

3.1. Construction of Composite Permeability Fields

The calcite image obtained after mineral precipitation was used to generate different permeability scenarios, which were then implemented in the solute transport model. The image were processed using the OpenCV library in Python [55,56]. Self-organized, heterogeneous porous media, characterized by elongated carbonate structures and preferential flow channels, were represented using two permeability distribution scenarios: a dual and a triple composite permeability field.. First, the dual composite permeability fields were created by processing the image according to the steps outlined in [51], which assumes that the precipitate amount is proportional to the logarithm of the pixel intensity: (a) converting the image from raw format (NEF) to a 16-bit image (TIFF); (b) crop the image to focus on the precipitated zone which resulted in a region of interest (ROI) with resolution of 3975 x 3000 pixels, each pixel with an area of $6.6\times {10}^{-3}$ mm²; (c) transforming red, green and blue pixel intensity ($R_{\mathrm{I}}$, $G_{\mathrm{I}}$, $B_{\mathrm{I}})$ into grayscale $\overline{I}(x,y,t)$ as

$$\overline{I}={\displaystyle \frac{R_{\mathrm{I}}+G_{\mathrm{I}}+B_{\mathrm{I}}}{3}},$$

(1)

which results in a single-band representation of the image with $\overline{I}\in [0,255]$; (d) normalization of the grayscale pixel intensity $\overline{I}$ to analyze image in the range [0,1] as

$$I_n(x,y,t)={\displaystyle \frac{\overline{I}(x,y,t)-{\overline{I}}_{\mathrm{back}}(x,y)}{{\overline{I}}_{\max }-{\overline{I}}_{\min }}},$$

(2)

where ${\overline{I}}_{\mathrm{back}}(x,y)$ is the grayscale light intensity distribution of a reference image representing a blank state obtained from the steady-state flow conditions prior to the MDP experiment, and ${\overline{I}}_{\max }$, ${\overline{I}}_{\min }$ correspond respectively to the maximum and minimum grayscale intensities found in the image; (e) apply the logarithm of the normalized intensity $I_n(x,y,t)$; f) create the histogram of the image (as shown in Figure 2a), to define the thresholds (Th) that will be used for the dual permeability fields; g) define four thresholds (Th) based on the percentile of the distribution of $log\left(In\right)$ values. The percentiles were chosen within the range $[-2,0]$, which represents the precipitated area. The first threshold, $T{h}_4$= -2, represents the $p_0$ and was chosen through visual inspection to define the total precipitate area. From this value, three additional thresholds were defined to representing the internal structure of the high-precipitation zones within the precipitated area. $T{h}_1$= -0.6, which represents the $p_{85}$ and covers 15% of the total precipitated area; $T{h}_2$= -0.75, which represents the $p_{70}$ and covers 30% of the total precipitated area; and $T{h}_3$= -0.9, which represents the $p_{55}$ and covers 45% of the total precipitated area. A threshold above $p_{85}$ represents a very small area, which can be disregarded by the model. While values below $p_{55}$ would represent a very large area, losing the ability to differentiate it from the total precipitate area defined by $T{h}_4$.; h) binarizing the image $B(x,y)$ based on each the selected thresholds as:

$$B(x,y)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}log\left(I_n(x,y,t)\right)\ge Th,\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

which results in four different images corresponding to the thresholds $T{h}_1$, $T{h}_2$, $T{h}_3$, and $T{h}_4$, where each pixel represents a value of 0 to represent the non-precipitated zone, or 1 to identify the precipitated zone; i) resize the binarized images (3000 x 3970) to the model size (200 x 265). Each pixel in the resized image $B_{\mathrm{m}}(i,j)$ represents the average of a 15x15 block from the original binarized images $B(x,y)$. If $B_m(i,j)\ge 0.5$, the pixel is set to 1; otherwise, if $B_m(i,j)<0.5$, it is set to 0; j) create maps of the four binary images representing the areas affected by precipitation in blue (1) and the unaffected areas in red (0) (as shown in Figure 2b, e).

The triple permeability fields were created by progressively combining the four resized (200x265) binarized dual permeability fields. Initially, Field 1 and Field 2 were combined, resulting in the Field 1-2. This process continued for all possible combinations, creating six different combinations (Field 1-2, Field 1-3, Field 2-3, Field 4-3, Field 4-2, and Field 4-1). In our approach, each pixel in the resulting composite permeability field was calculated by summing the corresponding pixel values from the two original images at the same (x, y) position. In the resulting triple composite permeability fields, 0 represents a non-precipitated zone, 1 represents a precipitated zone, and 2 indicates the region of highest precipitation. These fields were then represented as maps, where red represents areas unaffected by precipitation, green indicates areas affected by precipitation, and blue highlights areas with the highest precipitation (as shown in Figure 3). The Field 1-2, Field 1-3, and Field 2-3 do not include the total real area covered by the precipitation, so they only represent different scenarios of the internal structure of the precipitated layer. In contrast, the Fields 4-3, Field 4-2, and Field 4-1, by considering Field 4, represent on one hand the total area covered by the precipitation and, on the other hand, three different scenarios to discretize the internal structures (Field 1, Field 2, and Field 3).

Table 2. Pixel counts by field

Field	Red pixels	Blue pixels	Green pixels
1	51181	1819	-
2	48309	4691	-
3	44920	8080	-
4	34525	18475	-
2-1	48309	1819	2872
3-1	44920	1819	6261
3-2	44920	4691	3389
4-1	34525	1819	16656
4-2	34525	4691	13784
4-3	34525	8080	10395

Table 2. Pixel counts by field

Field	Red pixels	Blue pixels	Green pixels
1	51181	1819	-
2	48309	4691	-
3	44920	8080	-
4	34525	18475	-
2-1	48309	1819	2872
3-1	44920	1819	6261
3-2	44920	4691	3389
4-1	34525	1819	16656
4-2	34525	4691	13784
4-3	34525	8080	10395

3.2. Flow and Transport Model

The conservative breakthrough curves (BTCs) obtained before and after mineral precipitation were simulated through a flow and transport model. The flow model was developed under the assumption of steady-state conditions using MODFLOW-2005 [57], which solves the flow equation

$$\nabla ·\left(K\left(\mathbf{x}\right)\nabla h\right)=0,$$

(3)

where $K\left(\mathbf{x}\right)$ corresponds to the isotropic hydraulic conductivity, and h represents the piezometric head. The synthetic aquifer was simulated as a confined aquifer with a single layer, assuming steady-state condition. The dimensions of the flow model were 20 cm in width (W), 26.5 cm in length (L), and 1 cm in height (H). Flow and transport parameters implemented in the model are shown in Table 3. Cells of 1 mm in size were used, resulting in a grid with 200 rows and 265 columns. For the scenario representing the aquifer before mineral precipitation, a homogeneous model was applied using the initial hydraulic conductivity value measured in the experiment (${\mathit{K}}_{\mathit{0}}$ = 145.4 m/day). For the scenario representing the synthetic aquifer after mineral precipitation, a heterogeneous model was considered, utilizing either two or three hydraulic conductivities, depending on the model application scenario. For areas unaffected by mineral precipitation, the initial hydraulic conductivity value was used (${\mathit{K}}_{\mathit{0}}$). For the precipitated medium, the hydraulic conductivity values were determined through calibration. In the case of the dual-permeability model, only one hydraulic conductivity (${\mathit{K}}_{\mathit{1}}$) was calibrated, while for the triple-permeability model, two hydraulic conductivities representing the calcite layer were calibrated (${\mathit{K}}_{\mathit{1}}$ and ${\mathit{K}}_{\mathit{2}}$), where ${\mathit{K}}_{\mathit{2}}$ represents the areas of higher precipitation within the layer, and therefore the area of lower hydraulic conductivity. In terms of the flow boundary conditions, a prescribed flow boundary was established at the system entry, with a total constant inflow rate (Q= 0.0154 m³/day), distributed across 8 ports (injection wells) at the entry, each represented by 5 cells, and for the system outlet a constant head boundary condition was imposed (0.9 m above the tank elevation). The solute transport model was developed using MODPATH-RW [58], an extension of the original MODPATH code [59] that is capable to simulate the advection-dispersion equation (ADE), which assumes Fickian dispersion at the local scale, using the Random Walk Particle Tracking (RWPT) method,

$$\varphi {\displaystyle \frac{\partial c}{\partial t}}=-\mathbf{q}·\nabla c+\nabla ·(\varphi \mathbf{D}\nabla c),$$

(4)

where $\varphi$ is the medium porosity (dimensionless), c is the concentration of dissolved solute, $\mathbf{q}$ is the specific discharge and $\mathbf{D}$ is the hydrodynamic dispersion tensor, following the classic expression,

$$\varphi \mathbf{D}=({\alpha }_T\left|\mathbf{q}\right|+\varphi {\mathcal{D}}_m){\mathbf{I}}_d+\left({\alpha }_L-{\alpha }_T\right){\displaystyle \frac{\mathbf{q}{\mathbf{q}}^T}{\left|\mathbf{q}\right|}},$$

(5)

including both longitudinal (${\alpha }_L$) and transverse dispersivities (${\alpha }_T$), the identity matrix (${\mathbf{I}}_d$) and the apparent molecular diffusion ($D_m$), corrected by tortuosity. During the calibration of the homogeneous model (BTC obtained before mineral precipitation), both dispersivity coefficients (${\alpha }_L$ and ${\alpha }_T$) were adjusted. These values were then used to simulate the BTC obtained after calcite precipitation. The apparent molecular diffusion ($D_m$) was assigned a value of zero. The tracer test was simulated as a single specie, setting an injection time of 10 minutes in the eight injection wells at the entrance. The porosity ($\varphi$) was set at a value of 0.34, which corresponds to the value obtained experimentally.

3.3. Model Calibration

Parameter calibration of each individual permeability composite models are evaluated by using two statistical indicators to select the best result for each scenario. The first is the RMSE (Root Mean Squared Error), defined by equation RMSE, and the second is the coefficient of determination R², described in equation 7,

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2},$$

(6)

$${\mathrm{R}}^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{{\sum }_{i=1}^n{(y_i-\overline{y})}^2}},$$

(7)

where n is the total number of observations/predictions, $y_i$ are the observed values, ${\widehat{y}}_i$ are the values predicted by the model, and $\overline{y}$ is the mean of the observed values. The best result for each simulated scenario was selected considering on a high R² and a low RMSE.

We consider that changes in local dispersivities due to mineral precipitation are negligible compared to the associated changes in permeability. Therefore, local dispersivities were held constant for all models and calibrated based on the initial tracer test conducted prior to the Mixing-Driven Precipitation (MDP) experiment. The process was manual and guided by these statistical metrics to ensure an accurate fit of the arrival time, peak concentration, and overall shape of the BTC. As shown in Figure 4(a), the model provides a good representation of the experimental BTC, with a coincidence in the arrival time, a peak value similar to that of the experiment, and an excellent fit, as demonstrated by the statistical indicators (R²= 0.9939, RMSE= 0.01946). Regarding the tailing, the simulated BTC exhibits a slight delay compared to the experimental BTC, as seen in greater detail in Figure 4(b). This behavior may be attributed either to reduced sensor accuracy at low concentrations or to anomalous transport processes, as noted by Levy and Berkowitz [60], which can lead to tailing even in homogeneous media. The resulting values are shown in Table 4.

To select the best model among all the different permeability scenarios simulated, the Akaike Information Criterion (AIC) was used, calculated by equation 8.

$$\mathrm{AIC}=nln\left(S\right)+2M$$

(8)

where n is the total number of observations (562), M is the number of parameters adjusted in the model (2 for the model with 2 hydraulic conductivities and 3 for the model with 3 hydraulic conductivities), and S is the sum of squared residuals (SSR), which is calculated as

$$S=n\times {\mathrm{RMSE}}^2.$$

(9)

The model with the lowest AIC, considering all scenarios, including both dual-permeability and triple-permeability models, is considered the best, as it indicates the best balance between model fit and complexity.

4. Results

4.1. Dual Composite Permeability Models

The non-Fickian transport observed in the experimental BTC, following mineral precipitation, was first evaluated using four dual composite permeability fields. The best fit of simulated BTCs for each scenario are shown in Figure 5, with the corresponding model parameters listed in Table 5. The results are presented in order of the total precipitated area covered by the permeability field, with the precipitated area shown in blue and the non-precipitated area shown in red. BTCs are shown with concentration values normalized (C/$C_0$) and time expressed in pore volume ($P_V$). In addition, they are also shown on a logarithmic scale.

We begin describing the result obtained for the Field 1 (Figure 5a), which has the smallest area considered for mineral precipitation. The simulated BTC does not capture the anomalous transport behavior observed in the experimental BTC (Figure 5b). Instead, it exhibits Fickian behavior, with a delayed arrival and a significantly higher peak (C/$C_0$> 0.8). Regarding the tailing, the BTC does not account for the second peak observed in the experimental BTC and ends much earlier than the tail of the experimental BTC. It is interesting to highlight in this model that, although the BTC shows a Gaussian distribution, in Figure 5(c), it is evident that the low permeability zone results in an asymmetric distribution. To calibrate the BTC, the small area had to be assigned a very low hydraulic conductivity value (${\mathit{K}}_1$ = 8 m/day) for the precipitated area. resulting in the worst fit among the four fields ($R^2$= 0.7063, RMSE = 0.1022, AIC = 998.66). Field 2 (Figure 5d) covers a larger precipitated area. Although the simulated BTC still exhibits Fickian behavior, it begins to show slight characteristics of anomalous transport (Figure 5e). The arrival time is now closer to that of the experimental BTC, and the peak value is lower than the values reached in Field 1 ($C/C_0$< 0.8). The second peak, characteristic of the tailing, starts to appear subtly (Figure 5f). To calibrate the BTC, a slightly higher hydraulic conductivity (${\mathit{K}}_1$ = 17 m/day) was required, leading to a better fit compared to Field 1 (R²= 0.9180, RMSE= 0.0539, AIC= 279.52) . For Field 3 (Figure 5g), the simulated BTC no longer follows Fickian behavior, instead displaying characteristics of anomalous transport as reported in the experimental BTC (Figure 5h). Compared to Fields 1 and 2, the BTC better reproduces the arrival phase. However, the first peak remains poorly represented, with overestimated maximum values (C/$C_0$< 0.8), while the second peak in the tailing is now clearly observed, with values close to those obtained from the experimental BTC (Figure 5i). Calibration indicates that the BTC obtained for this scenario achieves a better fit than the BTCs for Field 1 and Field 2 (R² = 0.9595, RMSE= 0.0379, and AIC = -116.33). To calibrate the BTC to the considered area, a higher hydraulic conductivity value (${\mathit{K}}_1$ = 37 m/day) had to be assigned.

Finally, we present the simulated BTC for Field 4 (Figure 5j), which represents the entire area affected by precipitation. The simulated BTC exhibits non-Fickian behavior, more accurately capturing the anomalous transport observed in the experimental BTC (Figure 5k). Compared to the experimental BTC, it correctly represents the arrival, being the only one of the four fields to do so. Regarding the first peak, the model provides a good representation, slightly overestimating the values for the first peak ($C/C_0$∼0.7). Similarly to Field 3, the BTC captures the second peak but overestimates its values ($C/C_0$>0.1).

With respect to the tailing, as shown in Figure 5(k), the simulated BTC ends earlier than the experimental one. The BTC obtained for Field 4 (Figure 5j) provides the best fit among all fields (R²= 0.9838, RMSE= 0.0239, AIC= -634.57), requiring the highest hydraulic conductivity (${\mathit{K}}_1=60$ m/day). Although the results obtained using Field 4 provide excellent outcomes, the model has three key limitations: it fails to accurately represent the transition between peaks, underestimates the maximum values for both peaks, and does not fully capture the tailing behavior.

4.2. Triple Composite Permeability Models

We proceed with the analysis by evaluating the non-Fickian behavior, analyzing the results obtained from the triple composite permeability models. Here, we focus on improving the representation of the double peak transition, addressing the real values observed for both peaks, as well as correctly capturing the behavior of the tailing. The initial scenarios analyzed corresponded to the three fields that do not include the total real area covered by the precipitation, so they only represent different scenarios of the internal structure of the precipitated layer. The simulated BTCs for each scenario are shown in Figure 6, with the corresponding model parameters listed in Table 6. The colormap represents the hydraulic conductivity field, with blue indicating areas of higher mineral precipitation.

The first scenario corresponds to Field 1-2, which represents the smallest area affected by precipitation (refer to Figure 6a) from all triple composite permeability fields. In Figure 6(b), it can be seen that the simulated BTC does not accurately replicate the anomalous behavior. It exhibits a pattern similar to that observed for Field 2, with an arrival that does not match the experimental BTC, a higher peak value (C/$C_0$∼ 0.7), and a tailing that fails to reproduce the second peak, as shown in Figure 6(c). This BTC shows the worst fit among all the triple composite permeability fields (R²= 0.934, RMSE= 0.048, and AIC= -151.22). To calibrate the simulated breakthrough curve (BTC), the small area considered for precipitation had to be discretized into two regions with very low hydraulic conductivities (${\mathit{K}}_1$= 22 m/day and ${\mathit{K}}_2$= 13 m/day). Field 1-3 (refer to Figure 6d) covers a larger precipitated area compared to Field 1-2, and within the layer, a greater differentiation between hydraulic conductivities is observed, where the areas of higher precipitation represent only a small portion. The simulated BTC (Figure 6e) resembles the result obtained from Field 3 (binary field). It can reproduce the arrival time and the second peak in the tailing, but presents a slightly better fit for the first peak (C/$C_0$∼ 0.7). The simulated BTC for Field 1-3 (Figure 6d) achieves a better fit (R² = 0.964, RMSE= 0.035 and AIC= -203.8) than Field 1-2 by using higher hydraulic conductivity values ($K_1$ = 45 m/day and $K_2$ = 20 m/day) and considering a larger area affected by precipitation. Field 2-3 (refer to Figure 6g) covers a similar precipitated area compared to Field 1-3, but in this case, the areas of higher precipitation cover almost the entire layer, which prevents differentiation between the two hydraulic conductivities. The simulated BTC obtained by considering Field 2-3 (Figure 6h) presents a shape that visually resembles the experimental curve more closely than that of Field 1-3 but does not yield a better fit according to the calibration metrics (R² = 0.955, RMSE= 0.039, and AIC= -82.17). The arrival time resembles the shape obtained for the BTC for Field 3 and captures the second peak in the tailing, like Field 1-3, but it shows a weaker fit in reproducing the first peak (C/$C_0$∼ 0.7), as shown in Figure 6i. To calibrate the breakthrough curve (BTC), higher hydraulic conductivity values ($K_1$ = 65 m/day and $K_2$ = 30 m/day) were needed.

The evaluation proceeded with the next three proposed triple composite permeability fields, which included the total area covered by precipitation and three different approaches to represent the internal structure of the high-precipitation zones. The simulated BTCs for each scenario are shown in Figure 7, with the corresponding model parameters listed in Table 7.

Field 4-3 (see Figure 7a), compared to Field 4-1 and Field 4-2, covers a larger precipitated area to represent the internal precipitation structures formed within the calcite precipitate layer. The simulated breakthrough curve (BTC) yields better results than those obtained for the previous three cases with triple composite permeability fields (Figure 7b). In general, it exhibits behavior similar to that observed for Field 4. It successfully reproduces the arrival time and slightly overestimates the values obtained for the first peak (C/$C_0$∼ 0.7). Regarding the tailing, it represents the second peak but overestimates it (C/$C_0$∼ 0.1). As shown in Figure 7(c), the simulated BTC ends earlier than the experimental BTC. In terms of the fit, compared to Field 4, it presents a slightly better fit (R² = 0.988, RMSE = 0.019, AIC= -890.46). The best representation of the BTC (Figure 7b) was achieved using two different hydraulic conductivities, with the highest values considered for all triple composite permeability fields ($K_1$ = 88 m/day and $K_2$ = 48 m/day), which, on average, resemble the values used for Field 4 ($K_1$ = 60 m/day). Field 4-2 (see Figure 7d) covers a precipitated area that allows for better representation of preferential flow paths within the precipitated calcite layer, discretizing the internal precipitation structures more effectively. The simulated BTC (Figure 7e) shows behavior similar to that obtained for Field 4-3. It correctly reproduces the arrival time (Figure 7e) and maintains a slight overestimation of the first peak (C/$C_0$∼ 0.7), but it slightly improves the tail, better representing the second peak (Figure 7f). The transition between peaks improves significantly compared to Field 4, though it still does not represent the exact shape. The second peak tends to be slightly delayed relative to the experimental BTC, but this provides a better representation of the late-time behavior of the BTC, which improves significantly compared to Field 4-3. As observed in previous cases, the fit continues to improve as the analysis progresses, achieving the best results so far (R² = 0.985, RMSE= 0.022, and AIC= -725.68). For the BTC calibration, the hydraulic conductivities used were lower than those applied to calibrate the Field 4-1 curve ($K_1$ = 78 m/day and $K_2$ = 32 m/day).

Finally, we conclude with the results obtained for Field 4-1 (see Figure 7g), which represents the total real area covered by precipitation, along with the internal discretization of the calcite layer, considering only a small area of higher precipitation. The simulated BTC (see Figure 7h) represents the best result obtained across all field combinations, as it provides the best visual fit, along with the highest statistical values (R² = 0.991, RMSE= 0.017), and the lowest Akaike value (AIC= -1015.48). The obtained BTC successfully captures the anomalous behavior observed in the experimental BTC, accurately representing the arrival time and the first peak value (C/$C_0$∼ 0.67). Regarding the second peak, it is the BTC that best represents the real values, closely matching the experimental BTC (C/$C_0$< 0.1). As shown in Figure 7(i), the BTC provides an excellent representation of the tailing. To calibrate the BTC, the region with the highest mineral precipitation was assigned a low hydraulic conductivity value ($K_2$ = 17 m/day), while the rest of the layer, was assigned a higher value ($K_2$ = 73 m/day).

5. Discussion

The results obtained from the dual composite permeability models demonstrate that the extent of the precipitated area considered in the permeability field significantly influences the ability of the model to replicate the non-Fickian transport behavior observed in the experimental breakthrough curve (BTC). When using a permeability field that considers only a small area representing the highest precipitation zone, rather than the full extent of the calcite layer (Field 1), the model generates a BTC that exhibits Fickian behavior. This approach fails to accurately represent the arrival time, leads to an overestimation of the first peak, and does not capture the second peak that characterizes the anomalous behavior. When using a permeability field that represents an area smaller than the real precipitated area (Field 2 and Field 3), the BTC only partially represents the anomalous behavior. The arrival time improves slightly, but the value of the first peak is overestimated. There is significant difficulty in representing the second peak, with the tail consistently ending earlier than the experimental BTC. In contrast, when the total real extent of the affected area is considered in the permeability field (Field 4), the model is able to effectively capture non-Fickian behavior. This scenario provides an excellent fit for both the arrival time and the representation of the two peaks observed in the experimental curve. However, the model still struggles to fully reproduce the shape of the BTC during the transition between peaks and the maximum values of the first and second peaks, indicating limitations in capturing finer details of the observed behavior. The BTC obtained from Field 4 suggests that the first peak of the non-Fickian experimental BTC arises from bypass flow through non-precipitated zones, leading to an earlier arrival. Meanwhile, the second peak reflects flow concentration within the low-permeability calcite layer. This is the main reason why considering a permeability field with an area that does not correspond to the real extent of the precipitation leads to erroneous results. It would assume that flow is passing through non-precipitated areas, which results in an overestimation of the first peak.

On the other hand, the results obtained from the triple composite permeability models suggest that the best way to represent the non-Fickian BTC is by accurately defining the actual precipitated area and the internal structure of the high-precipitate zone. When a smaller area than the actual precipitated area is defined, even if the internal zone of the layer is discretized into two distinct hydraulic conductivities (Field 1-2, Field 1-3, and Field 2-3), the values of the first peak will always be overestimated. This follows the same reasoning as for the dual composite permeability fields: failing to account for the full precipitated area leads to less flow being concentrated in the considered layer and, consequently, an increased flow in the non-precipitated zone. An interesting observation arises when analyzing Field 1-2. Despite having a very small precipitated area, the internal discretization of the layer into two distinct hydraulic conductivities, unlike in Field 1, allows for the representation of a non-Fickian BTC. Although this does not provide a fully accurate representation, the analysis of Field 1 highlights the importance of internal discretization in capturing the anomalous BTC behavior, as it enables the reproduction of the two peaks observed in the experimental BTC. Further insights emerge when examining the BTC obtained for Field 1-3, which achieves a better fit for the tailing, as seen more clearly in the logarithmic scale. This suggests that tailing can be more accurately represented by discretizing only a small percentage of the layer, specifically targeting the areas with the highest precipitation. When analyzing the fields that consider the total precipitated area but differ in how they represent the preferential flow paths (Field 4-3, Field 4-2, and Field 4-1), better results are achieved, as all of them successfully capture the first peak. The best results are obtained when the internal zone of the layer represents a smaller area, meaning that focusing only on the regions with the highest calcite precipitation leads to improved accuracy. This reduces the overestimation of the first peak, enhances the transition between peaks, improves the representation of the second peak, and extends the tailing over a longer period. However, achieving these improvements requires a strong contrast in the assigned hydraulic conductivity values to clearly differentiate low-permeability zones.

6. Conclusions

Mixing-induced precipitation experiments conducted in an intermediate-scale Hele-Shaw cell resulted in the formation of a self-organized, heterogeneous porous medium, characterized by elongated carbonate structures and internal preferential flow channels aligned with the main flow direction. After precipitation, the obtained solute transport breakthrough curve (BTC) exhibited strong anomalous behavior, as indicated by earlier solute arrival, a distinct double peak, and pronounced tailing. In this paper, we investigated the relationship between the self-organized heterogeneous structure of the porous medium formed through mixing-induced precipitation and its impact on solute transport. To achieve this, we analyzed the spatial variability of hydraulic conductivity by implementing different permeability scenarios in a random walk particle tracking model. These scenarios, derived from image analysis of the precipitated structures, ranged from simple representations to increasingly complex configurations. Two different permeability distribution scenarios were considered: a dual composite permeability field and a triple composite permeability field. Finally, we conclude that:

1.

Two important features must be captured to accurately reproduce non-Fickian solute transport observed after mineral precipitation: a) identify the area of total extension of the precipitation; and b) represent the internal heterogeneous structure conducive to preferential channels.

2.

Dual composite permeability fields, which incorporate the entire area affected by precipitation, effectively capture non-Fickian behavior in the model. This approach provides an excellent fit for both the arrival time and the two peaks characteristic of dual-permeability media. However, the model still faces challenges in fully capturing the transition between the peaks and the maximum value of the first peak. The model suggests that the first peak is due to bypass flow through non-precipitated zones, leading to earlier arrival times. The second peak, on the other hand, reflects the concentration of flow within the low-permeability calcite layer.

3.

The best way to represent the non-Fickian behavior observed after mineral precipitation was achieved by considering triple composite permeability fields, taking into account the total actual area of the precipitate along with defining the high-precipitation area as a small region. In this way, the model reduces the overestimation of the first peak, enhances the transition between peaks, improves the representation of the second peak, and extends the tailing over a longer period.

This study highlights the critical role of precipitation induced heterogeneity in solute transport, emphasizing how varying levels and distributions of precipitate, driven by mineral precipitation processes, impact transport behavior. By integrating permeability fields that capture key aspects of the spatial variability of hydraulic conductivity caused by calcite precipitation into transport models, this research provides valuable insights applicable to both field studies and the improvement of multiple domain models. These findings are essential not only for understanding real-world systems affected by transient precipitation processes but also for practical applications in environmental engineering. Specifically, they can help enhance strategies for carbon sequestration and enhanced oil recovery (EOR), where controlling solute retention and mobility is crucial.