A Lab on a Chip Experiment for Upscaling Diffusivity of Evolving Porous Media

Abstract

Reactive transport modelling is a powerful tool to assess subsurface evolution in various energy-related applications. Upscaling, i.e., accounting for pore scale heterogeneities into larger scale analyses, remains one of the biggest challenges of reactive transport modelling. Pore scale simulations capturing the evolutions of the porous media over a wide range of Peclet and Damköhler number in combination with machine learning are foreseen as an efficient methodology for upscaling. However, the accuracy of these pore scale models needs to be tested against experiments. In this work, we developed a lab on a chip experiment with a novel micromodel design combined with operando confocal Raman spectroscopy, to monitor the evolution of porous media undergoing coupled mineral dissolution and precipitation processes due to diffusive reactive fluxes. The 3D-imaging of the porous media combined with pore scale modelling enabled the derivation of upscaled transport parameters. The chemical reaction tested involved the replacement of celestine by strontianite, whereby a net porosity increase is expected because of the smaller molar volume of strontianite. However, under our experimental conditions, the accessible porosity and consequently diffusivity decreased. We propose a transferability of the concepts behind the Verma and Pruess relationship to be applied to also describe changes of diffusivity for evolving porous media. Our results highlight the importance of calibrating pore scale models with quantitative experiments prior to simulations over a wide range of Peclet and Damköhler numbers of which results can be further used for the derivation of upscaled parameters.

1. Introduction

Reactive transport modelling (RTM) has become a powerful and indispensable tool for predicting subsurface flow and transport processes in space and time for various energy-related applications including geothermal energy extraction [1], CO₂ sequestration [2,3], exploitation of natural gas reservoirs [4], or geological disposal of nuclear waste [5,6,7]. Although RTM can solve complex couplings between multi-phase fluid flow and mineralogical reactions [8], e.g., CO₂ gas injection into deep saline aquifers with feedbacks of mineralogical reactions on transport parameters [9], the empirical formulations to describe coupled processes are debated [10]. Indeed, RTM rationalizes the complex processes in porous media on a continuum-scale and considers a representative elementary volume [11] with the average chemical and transport properties of a rock, e.g., porosity and diffusivity, while the macroscopic behaviour of the system is controlled by pore scale features, e.g., pore connectivity and surface roughness of minerals. One of the big challenges of RTM is the development of constitutive equations for the effective transport properties of porous media and the reactivity of minerals that can be integrated into continuum analyses, i.e., “upscaling” [12,13,14].

The effects of mineral dissolution and precipitation on transport parameters in RTM are considered via empirical constitutive equations, e.g., Archie’s law [15] for diffusivity and the Kozeny Carman equation [16] for permeability. Experimental benchmarks [17,18,19,20,21] conducted to test these implementations showed that these relationships underestimate the impact of precipitation processes on transport parameters. Pore scale modelling approaches, with the versatility to define the pore architecture of a rock matrix down to the nanometre scale [22], the mineralogical distribution at the grain scale including surface roughness [23], were developed to investigate how the pore structure is affected by dissolution [24,25] and precipitation processes [26,27,28,29]. Pore scale modelling has been applied to derive transport parameters and constitutive equations that describe transport properties of evolving porous media. For example, Deng et al. [30] conducted a series of numerical experiments to investigate the impact of celestine nucleation mechanisms, surface precipitation (heterogenous nucleation) versus precipitation in pores (homogenous nucleation) on the transport properties of an inert porous medium. The study showed that classical Archie’s law cannot be applied for the studied system and proposed an extended law introducing the concept of critical porosity [31] and critical effective diffusivity, which refer to the inherent transport parameters of the altered porous medium after porosity clogging. To increase the efficiency of this approach, machine learning tools have been recently developed to parameterize constitutive equations [32,33]. The numerical datasets generated by pore scale modelling are used to train convolutional neural network models, which can then be further applied to derive constitutive equations. For example, Prasianakis et al. [34] used neural networks to investigate porosity–permeability relationships of a calcite containing porous medium undergoing dissolution. The numerical experiments covered a wide range of flow and reaction rate conditions, i.e., diverse Peclet and Damköhler numbers. Their model showed the inability of a single porosity–permeability relationship to be applied for the complete flow and reactivity conditions of the numerical experiments but instead proposed a fitted polynomial to characterize the permeability versus porosity evolution.

The derivation of constitutive transport equations is also possible by applying pore scale modelling on high quality images of the porous media. Steinwinder and Beckingham [35] applied pore scale modelling to X-ray tomography images of an Alberta Basin Sandstone that had undergone coupled mineral dissolution and precipitation to examine the impact of the variations of pore and pore throat size distributions on the evolution of porosity and permeability. The study highlighted the inappropriateness of the power law of the Kozeny Carman equation and proposed the Verma and Pruess relationship [36] or its extended version to reproduce the macroscopic permeability of the porous rock. The Verma and Pruess relationship introduces an additional parameter, the critical porosity, which in fact accounts for heterogeneities at the pore scale, e.g., localized clogging zones [37].

Various imaging methodologies have been developed for the visualization of reactive flows in porous media, including positron emission tomography [38], magnetic resonance imaging [37,38,39], and neutron tomography [40]. However, X-ray tomography is the most widely applied tool [3,41,42,43]. The principal reason is that beside conferring advantages such as non-invasive, non-destructive, and dynamic assessment of flow paths, it can also assess mineralogical changes. This identification is possible based on the attenuation coefficient of each material in the sampled material. While high atomic number containing minerals exhibit high attenuation and are therefore easily distinguishable, discerning lighter atoms containing minerals with similar densities can be challenging. For instance, a sandstone consisting of a mixture of calcite, gypsum, and quartz is complex to differentiate by X-ray tomography due to the very similar attenuation coefficients of the three minerals [44]. Recently the use of confocal Raman spectroscopy has been proposed as a non-destructive method to monitor the spatio-temporal evolution of porous media, e.g., assessing the effects of hydrological heterogeneities on nucleation mechanisms [45] or for identifying minerals during coupled dissolution and precipitation processes [46]. Given 4D capabilities and the high spatial resolution of Raman imaging (~500 nm), mineral identification can be used as a complementary analytical tool to X-ray tomography when the identification of mineral phases using the latter is demanding.

The aim of this work is therefore to develop an experiment which employs confocal Raman spectroscopy to monitor transport coupled mineralogical reactions, enabling the upscaling of diffusive transport parameters of evolving porous media. We focussed our study specifically on reactive media undergoing coupled mineral dissolution and precipitation as existing work so far have mostly tackled diffusion induced mineral precipitation [19,20,21,47]. Our experimental approaches follow the lab on a chip concept [45,48,49] whereby reactions occurring in microreactors are monitored by advanced analytical techniques.

2. Materials and Methods

2.1. Setup and Microreactor Design

The experimental setup consists of a microfluidic experiment monitored by high-resolution optical microscopy and confocal Raman spectroscopy. The micro-reactor was mounted on automated x-y-z stage and connected to a pump (NeMeSYS, Cetoni GmbH, Korbußen, Germany) and an effluent vessel as shown in Figure 1a. The 3D single use PDMS chip bonded to glass was designed and fabricated according to Poonoosamy et al. [47], consisting of four microreactors (Figure 1b) [48]. The supply channel of 150 µm width and a length of 12.7 mm connects the inert porous medium (Figure 1c) to a reservoir filled with a reactive porous medium (Figure 1d).

The inert porous medium was 1.92 mm long and 600 µm wide and consisted of an array of cylindrical pillars of 100 µm diameter distanced by 1 µm from each other. The reservoir had dimensions of 600 µm × 300 µm with 5 squared pillars integrated to maintain mechanical support. The depth of each microreactor was 10 µm. The reservoir was filled with naturally occurring celestine (SrSO₄) crystals with low porosity crushed to a size of 4–9 µm (see Poonoosamy et al. [17] for analysis of purity). A suspension of the crystals was injected via the crystal injection channel following the seeding procedure developed in Poonoosamy et al. [45]. The crystals were trapped in the reservoir by the pillars (distance 1 µm to each other), which acted as a sieve. The reactors were filled with celestine with initial porosities between 0.60 to 0.92 (see video in Supplement Video S1) of which the initial effective diffusivities were determined using pore scale modelling (cf. Section 2.4). The reactor with an initial porosity adjusted to 0.92 (average value for three sampled locations) was used for the reactive transport experiment (cf. Section 2.2).

2.2. Chemical Reactions

A solution of 1 M of sodium hydrogen carbonate (NaHCO₃) at pH ~ 9.0 was injected in the supply channel at a flow rate of 800 nL min⁻¹ for 24 h, leading to the following chemical reaction in the reactive porous medium:

$${\mathrm{CO}}_3{{}^{2-}}_{\mathrm{aq}}+{\mathrm{SrSO}}_{4\left(\mathrm{s}\right)}\to {\mathrm{SO}}_4{{}^{2-}}_{\mathrm{aq}}+{\mathrm{SrCO}}_{3\left(\mathrm{s}\right)}$$

(1)

The injection of the NaHCO₃ solution induces the dissolution of celestine and the precipitation of strontianite (SrCO₃) (Equation (1)), because of its much lower solubility product (cf.

Table 1. Properties of the different zones and composition of the injected solution in the 1D continuum scale reactive transport simulation.

Zone	Properties
Qa (inert porous medium)	length 1.9 × 10−3 m initial porosity, Φ, 0.21 volume fraction of inert matrix 0.79 De,initial 2.1 × 10−10 m2 s−1 exponent (Equation (6)), m = 1
Qb (intermediate zone)	length 5 × 10−5 m inital porosity, Φ, 0.19 volume fraction of celestine 0.1 volume fraction of inert matrix 0.71 De,initial 1.9 × 10−10 m2 s−1 exponent (Equation (6)), m = 1
Qc (reactive porous medium)	length 3 × 10−4 m initial porosity, Φ, 0.92 volume fraction of celestine 0.08 De,initial 8.3 × 10−10 m2 s−1 exponent (Equation (6)), m = 2.3
concentration NaHCO3	1 M with pH ~ 9
initial pH of porous medium	5.6
temperature	25 °C
solubility product, Ksp, celestine [50]	2.24 × 10−7 mol2 L−2
solubility product, Ksp, strontianite [50]	5.6 × 10−10 mol2 L−2

). This reaction is accompanied by an increase in “net porosity” (used in the following for porosity changes based solely on the molar volumes of the minerals involved) due to the smaller molar volume of strontianite (3.901 × 10⁻⁵ m³ mol⁻¹ [50]) compared to that of celestine (4.625 × 10⁻⁵ m³ mol⁻¹ [50]).

Prior to the injection of the reacting solution, the microreactor was purged with a saturated strontium sulphate solution for 24 h at a constant flow rate of 800 nL min⁻¹, during which the initial 3D representation of the reactive porous medium was collected using Raman spectroscopy (cf. Section 2.2). To prevent any reactions in the supply channel, pure water was injected at 800 nL min⁻¹ for 30 min before injecting the reacting solution. The experiment was conducted at room temperature (22 °C).

2.3. Imaging and Raman Spectroscopy

The microfluidic experiment was monitored using an automated inverted microscope (WITec alpha300 Ri Inverted Confocal Raman Microscope, WITec Wissenschaftliche Instrumente und Technologie GmbH, Ulm, Germany) with a 100× oil immersion objective with a numerical aperture (NA) of 1.25, working distance of 0.23 mm and cover glass correction for 0.17 mm from Nikon. The instrument is equipped with a Nd: YAG laser (λ = 532.236 nm), a thermoelectrically cooled CCD Camera, and an ultra-high-throughput spectrometer UHTS300. The laser power was set to 30 mW and a grating with 1800 grooves mm⁻¹ was chosen, yielding a spectral resolution of 2 cm⁻¹. The theoretical, diffraction-limited lateral and axial resolutions of the Raman measurements at the sampled surface were ~520 nm and ~2060 nm, respectively, considering Equations (3) and (4) in Everall [51], with a refraction index of the immersion medium of 1.51. The experimental setup is equipped with an enclosure to minimize interferences with room temperature fluctuations as well as cosmic rays.

Before the reaction, high resolution optical images of the compacted porous medium were taken, and the images were stitched for visualization of the structure of the reactive porous medium (Figure 1c,d). Three regions of interest (e.g., areas labelled 1 to 3 in Figure 1c) were selected and the mineralogical changes were monitored at regular time intervals using Raman spectroscopy. In situ hyperspectral Raman images were recorded from an area of 100 μm × 30 μm with a 500 nm step size (200 pixel × 60 pixel per image) in x and y direction. The individual Raman spectra were collected during continuous x–y stage movement. Raman intensities were recorded for 0.1 s in the wavenumber range from 300 to 1300 cm⁻¹. The total exposure time was 25 min for a single image. In addition, stacked Raman images of the same size were collected over a depth of 11 μm with a 1 µm step size. These images were used to reconstruct the 3D geometry of the skeleton of the porous medium. The voxel resolution of our measurement is 500 µm³. The heat introduced by the laser during the 3D Raman tomography has a negligeable effect on the dissolution of celestine as discussed in Poonoosamy et al. 2020 [45].

Raman image stacks were visualized with ImageJ 3D Viewer plugin.31. An image processing workflow was developed using ImageJ to measure the porosity (voxel without minerals) and volume fraction of each mineral from the stacks. For that, each mineral phase was converted to greyscale levels and segmented by default thresholding into two classes of pixels (“foreground” and “background”). In addition, the images were converted into binary text files using Matlab to simulate the flow through the porous media (Section 2.5).

2.4. Velocity Field in the Microfluidic Reactor

The 3D flow field in the microfluidic device was simulated with computational fluid dynamics using the software COMSOL Multiphysics 5.6 (COMSOL AB, Stockholm, Sweden, see Appendix A). This step enables the characterization of the transport regime that prevails in the micromodel.

2.5. Numerical Tracer Experiment Using Pore Scale Modelling

The effective diffusivity of the unreacted and reacted porous media were determined using pore scale modelling. Numerical inert tracer experiments were conducted on the experimental 3D reactive porous medium collected by Raman imaging (Figure 2a; the intermediate section with pillars was removed for this analysis, i.e., the simulations addressed only the pure reactive porous media). For the reacted porous medium, it is known that solute transport can also occur through the newly formed precipitates [52,53] and therefore diffusion through the precipitates is considered in our numerical investigations. The lattice Boltzmann method, LBM, which relies on the application of the Boltzmann transport equation to simulate fluid flow on discrete grids [54] was used for this purpose. Benefiting of the high efficiency of LBM for parallel computing, it has shown amazing performances to analyse fluid flow or heat and mass transfer in complex porous media during recent years [55,56].

We use LBM to solve the diffusion equation in the 3D porous media domains (from the Raman tomographs):

$$\frac{\partial C_p}{\partial t}=\nabla ·\left(D_i·\nabla C_p\right)$$

(2)

where $C_p$ denotes the local concentration of the inert tracer at pore scale. $D_i$ is the local diffusion coefficient given by:

$$D_{\mathrm{i}}=\{\begin{array}{cc}{D}_0^{},& in the eletrolyte in the pores\\ D_s^{},& in the newly formed precipitates\end{array},$$

(3)

where $D_0$ is the tracer’s diffusivity in free water and $D_s$ its effective diffusivity in the precipitates forming the existing solids precipitations (i.e., strontianite in this experiment). A constant tracer concentration of 2 mM was applied at the inlet. When the process of tracer diffusion reaches a steady state, the effective diffusion coefficient, D_e [m² s⁻¹], can be calculated using the total mass flux $J$ through the entire domain:

$$D_e=\frac{J·L}{S\left(C_{in}-C_{out}\right)},$$

(4)

with a length of the domain $L$ [m] and a total cross-sectional area $S$ [m²]. Equations (2) and (4) are numerically solved by our previously developed framework based on LBM, which has been validated in these references [22,57].

For the evaluation of the effective diffusivity (D_e) of the reactive porous medium that had undergone mineral dissolution and precipitation, two numerical case studies were considered, case 1 without solute transport through the newly formed strontianite (D_s = 0 m² s⁻¹) and case 2 with transport through the newly formed precipitate (D_s/D₀ = 0.01). This factor, rounded to a 2-digit number, is calculated from Deng et al. [30], who fitted an inherent diffusivity of the precipitates of 5 × 10⁻¹² m² s⁻¹ to match the experimental datasets of Chagneau et al. [19].

The correlation with the domain porosity, Φ, was also investigated for the derivation of a constitutive equation relating porosity and effective diffusivity.

2.6. Reactive Transport Model Setup at the Continuum Scale

In a final step, we aimed to compare the extended porosity-diffusivity model derived from 3D Raman images and pore scale modelling (Section 2.3 and Section 2.5) with the classical Archie’s law, which is applied to account for porosity changes due to mineralogical reactions on the solute diffusivity in the liquid phase. For this purpose, we conducted a 1D reactive transport simulation using the finite element package OpenGeoSys-6#Phreeqc [58,59], where an integration-point operator-splitting scheme [60] was used to solve the multi-species diffusion-reaction equations as follows [61,62]:

$$\frac{\partial \left(\Phi c_{T\alpha }\right)}{\partial t}=\nabla ·\left(D_{e,\alpha }\nabla c_{T\alpha }\right)+R_{\alpha }^{\min }\left(c_{T1},\dots ,c_{Tm}\right),$$

(5)

where $c_{T\alpha }$ ($\mathsf{\alpha }=1,\dots , m$) (mol m⁻³) is the total concentration of the ${\alpha }^{\mathrm{th}}$ primary species in the liquid; $D_{e,\alpha }$ (m² s⁻¹) is the effective diffusion coefficient of the primary species $\alpha$ (assumed to have the same value for all the species in our model); and $R_{\alpha }^{\min }\left(c_{T1},\dots ,c_{Tm}\right)$ is the heterogeneous reaction term accounting for the contributions of all mineral dissolution-precipitation reactions that involve the primary species $\alpha$.

The computational domain consists of three zones (Figure 2b). The leftmost zone ${\mathrm{Q}}_{\mathrm{a}}$ is made up of an inert pillared porous medium with a porosity of 0.21. The rightmost celestine-containing reservoir Q_c has a porosity of 0.92. The intermediate 50 µm zone, Q_b, onsisting of pillars with stacked crystals, has a porosity of 0.19. A 1 M NaHCO₃ solution (pH ~ 9) was allowed to diffuse through the microfluidic device into the reacting reservoir represented by Q_c. In Table 1, the properties of the porous media in the different zones and the prescribed boundary conditions are summarized. The entire domain was evenly discretized using linear line elements with a grid spacing of 50 µm. A fixed time step size of 0.2 s was used in the implicit backward Euler scheme.

Usually, the effective diffusion coefficient $D_e$ is parameterized as a function of porosity. In this simulation, we use the classical Archie’s law [15] to parameterize the effective diffusion coefficient [19]:

$$D_e=D_p{\Phi }^m,$$

(6)

where $D_p$ is the pore diffusion coefficient set to 1 × 10⁻⁹ m² s⁻¹, Φ [-] is the porosity, and m [-] is an empirical coefficient. In this study, m was set to 1 for the inert zone (Q_a) and the intermediate zone (Q_b) and to 2.3 (a fitted value for a sandy compacted porous medium [19]) for the reactive porous medium in zone Q_c.

In our continuum scale simulations, we used a simple reaction system where the dissolving celestine and precipitating strontianite were assumed to equilibrate simultaneously against the incoming solution (i.e., as soon as oversaturation is reached). The PSI/NAGRA 12/07 chemical thermodynamic database (summarized in Table 3 in [63]) was used for the speciation calculation.

3. Results

3.1. Evaluation of Velocity Field

The injection of the NaHCO₃ solution in the supply channels at constant rate yields a velocity field at steady state as shown in Figure 3a. The velocity is highest at ~1.2 × 10⁻³ m s⁻¹ in the solution supply channels and decreases to ~10⁻¹³ m s⁻¹ in the inert porous medium (i.e., left end of the Q_a region). The Peclet number (P_e, defined as the ratio of advective to diffusive transport rates) in the inert porous medium close to the reservoir was evaluated at 0.0001, with the consideration of a characteristic length of 1000 μm (i.e., a length over which the velocity is in the same order of magnitude). It can thus be inferred that diffusion dominates the transport of the reacting solutes into the reactive porous medium. Figure 3b shows the simulated concentration of sodium (Na, conservative solute transport) across the porous media with time. Thus, it takes three hours for a steady sodium concentration gradient to be established across the inert porous medium.

3.2. Mineralogical and Porosity Changes

The mineral phases in the reservoir were identified by the characteristic Raman-active modes. The spectra of strontianite and celestine are distinct. Figure 4a shows the temporal evolution of the distribution of celestine (identified in red) and newly formed strontianite (identified in yellow) in the three sampled areas labelled 1 to 3 (cf. Figure 1c). The dissolution of celestine followed by precipitation of strontianite is a fast process with only 14 ± 4% of the initial celestine remaining undissolved after 6 h of reaction (Figure 4b) and a complete replacement observed within 10 h. From the Raman images, precipitates were observed to initiate on the surface of larger crystals of celestine (e.g., in area i in Figure 4a) and as a very fine “fluffy” precipitate that recrystallizes into larger crystals (e.g., in area ii in Figure 4a). No precipitation was observed in the inert porous medium.

The accessible porosities, which we defined as the fraction of voxels unoccupied by minerals before and after reaction were calculated from the 3D Raman tomographic images of the porous media (Figure 5a) and are given in Figure 5b. The volume occupied by the newly formed precipitates is larger than the pristine celestine crystal (see video in Supplement Video S2). Thus, the average measured “accessible” porosity of the reacted porous medium decreases to ~0.54. In contrast, based on the difference in molar volume between celestine and strontianite, the final porosity or “net porosity” should increase from ~0.92 to ~0.93. Hence, the replacement of celestine by strontianite involves an increase in net porosity but a decrease in accessible porosity for transport.

3.3. Evaluation of the Effective Diffusivity

The tracer concentration in the reactive porous medium before and after reaction simulated using the lattice Boltzmann method is shown in Figure 6. The tracer concentration map across the domain before and after coupled dissolution and precipitation indicates a new pore architecture, whereby the final effective diffusivity is lower by a factor of ~0.5 (averaged for the three areas) as the initial effective diffusivity. When diffusion is allowed through the newly formed strontianite with D_s/D₀ = 0.01 (Figure 6c), this factor increases slightly to 0.52 (averaged for the three areas). In addition, sensitivity analyses, with D_s/D₀ equal to 0.04 and 0.1 were conducted. All simulated results indicated a decrease in the effective diffusivity of the reacted porous medium with the being less pronounced as D_s/D₀ increases (see tabulated datasets in Appendix C).

The correlation between the simulated effective diffusivities and the calculated porosities from the 3D images is shown in Figure 7. For the unreacted porous media of porosities between 0.60 to 0.92, the associated diffusivities are indicated on the graph by the red data points. Our analysis shows that Archie’s law parameterized for compacted sandy porous material reproduces well the simulated effective diffusion coefficient of the unreacted celestine medium at different porosities. In addition, we report the accessible porosities of the chemically altered media and the associated simulated effective diffusivities (grey data points in Figure 7). Here, the simulated effective diffusivities of the reacted porous media are rather invariant with respect to the assumption regarding the diffusivity of the newly formed strontianite precipitates.

In a continuum scale reactive transport model, the changes in diffusivities due to mineralogical reactions can only be accounted for based on the net porosity changes calculated from changes in molar volumes. Consequently, the net porosities (black data points in Figure 7) have to be considered instead of the accessible porosities in continuum scale reactive transport simulations. Here, Archie’s law fails to reproduce the experimental observations. Therefore, we tested an alternative approach for describing the porosity—diffusivity relationship. As mentioned earlier, Deng et al. [30] conducted a series of numerical experiments on diffusion induced precipitation of celestine and proposed the following porosity—diffusivity relationship:

$$D_e=D_{effc}+a{D}_p{\left(\Phi -{\Phi }_c\right)}^n,$$

(7)

where a and n are empirical coefficients, ${\Phi }_c$ and $D_{effc}$ are defined as the critical porosity and effective diffusivity post porosity clogging, respectively. In our experiment, where no porosity clogging is observed, we can assume that $D_{effc}$ = 0, and therefore the ratio of the initial (D_e,initial) and final (D_e,t) effective diffusivity can be simplified to:

$$\frac{D_{e, t}}{D_{e, initial}}={\left(\frac{{\Phi }_t-{\Phi }_c}{{\Phi }_0-{\Phi }_c}\right)}^n,$$

(8)

where $D_{e, initial}$ and ${\Phi }_0$ are the initial diffusivity and porosity, respectively. ${\Phi }_c$ and n are fitted to 0.65 and −6, respectively, to match the experimental data sets.

3.4. Reactive Transport Simulation at the Continuum Scale

Reactive transport modelling at the continuum scale was performed to evaluate the impact of coupled celestine dissolution and strontianite precipitation and consequences on the diffusivity evolution, using (i) the classical Archie’s law (Equation (6)) and (ii) using the extended diffusivity relationship (Equation (8)). The simulated mineralogical and associated porosity changes are the same for both studied cases and are presented in Figure 8a,b, respectively. The simulated mineralogical reactions with time are in good agreement with our experimental observations. The calculated Damköhler number, calculated using Equation (A2) in Appendix B, is around ~10⁵ indicating that the system in our experimental conditions is controlled by the reaction given by the very high concentration of NaHCO₃ injected.

Figure 8c,d show the simulated diffusive properties of the porous domain based on the two porosity-diffusivity models with a simulated diffusivity increase with Archie’s law and a diffusivity decrease with the extended diffusivity relationship (Equation (8)).

In the current study, the use of Archie’s law and the extended law had little impact of the simulated net mineralogical reactions within the considered domain and time scale. However, if the domain was extended to cm or m scale, differences between the mineralogical profiles and consequently porosity across the domain will develop.

4. Discussion

Archie’s law parameterized to describe diffusion in a compacted sandy porous medium reproduces well the diffusivity of the unreacted celestine. When precipitation occurs, the precipitates cement the porous media creating a new pore architecture and the parameterized Archie’s law is no more valid. The replacement of celestine by strontianite caused an unexpected decrease in accessible porosity which suggested that the newly formed precipitates were porous. The volume increase can be explained by the presence of fluid inclusions which are common for minerals precipitating at high saturation ratios [64,65], here calculated ~700 (cf. Appendix B for definition) and by the precipitation patterns. At high saturation ratios, both homogeneous (HON) and heterogeneous (HEN) nucleation can occur at competitive rates (9 × 10²² m³ s⁻¹ for HON versus 2 × 10¹⁹ m³ s⁻¹ for HEN using equations given in Appendix B) [66]. Heterogenous nucleation of celestine appears to be favoured on the larger celestine crystals. This process, termed interface coupled dissolution and precipitation [67,68], can lead to the formation of encrustation pseudomorphs [17] (cf. area i in Figure 4a), i.e., a crystal with a void inside occupying more pore space than the original celestine crystal. On the other hand, the dissolution of the smaller celestine crystals can be a source of Sr²⁺ to fuel homogeneous nucleation. HON resulted in the formation of the fine “fluffy” strontianite phase (cf. area ii in Figure 4a). The very fine strontianite seemed to undergo Oswald ripening forming larger cemented strontianite crystals, creating a new pore architecture. As such, the precipitation patterns can significantly reduce the experimentally measured accessible porosity and consequently the simulated diffusivity. Similar observations were reported for celestine dissolution followed by barite precipitation [64].

The newly formed precipitate will still enable diffusion of solutes to the surface of the initial minerals because of its inherent porosity [29,52], nevertheless, overall, the effective diffusivity decreases. The presented replacement study highlights the complexity of predicting coupled dissolution and precipitation processes and associated changes in transport properties of porous media. A similar behaviour was observed by Sabo and Beckingham [69] during which a coupled dissolution and precipitation with a net porosity increase was accompanied by a decrease in permeability. The authors proposed the Verma and Pruess relationship combined with a negative exponent to rationalize their observations. The Verma and Pruess relationship has already proven to reproduce permeability changes due to coupled mineral dissolution and precipitation processes [17,69,70] and should be considered in continuum scale reactive transport modelling. In fact, the proposed Equation (8) to characterize diffusivity is an adaptation of the Verma and Pruess relationship. Currently, there exist no tabulation for the range of values taken by the critical porosity or the exponent for diffusive transport. However, the fitted values to match the presented experimental datasets are within the range reported by Hommel et al. [70] for describing permeability changes as consequence of altered porosity.

In addition to the properties of the initial porous medium, the parameterization of equations relating porosity to transport parameters is also known to depend on Peclet and Damköhler numbers [71,72,73]. To develop process based predictive models, a compilation of datasets for a wide range of Peclet and Damköhler numbers is necessary [73]. However, this might be a tedious work if such parameterizations are based solely on experiments. A more optimized approach is to calibrate pore scale models with a few reactive transport experiments at the same scale, which can point out the unexpected system behaviour. Afterwards, these calibrated models can then be used to simulate various flow and reaction conditions, extending datasets over a wide range of Peclet and Damköhler numbers. In the absence of experimental datasets on evolving porous media, it is important to include the extended law as sensitivity case when modelling real-world larger scale problems, which should set a range of uncertainty of the coupled reactive transport modelling.

5. Conclusions

We developed a lab on a chip experiment to monitor and upscale transport parameters of evolving porous media. Our experiment, involving the replacement of celestine by strontianite, yielded a net porosity increase but a decrease in accessible porosity for transport and consequently diffusive transport. Conventional Archie’s law failed to reproduce the experimental results, but rather an extended law based on the concept of the Verma and Preuss relationships was applied to reproduce the experimental observations.

The presented micromodel design can be used to assess coupled dissolution precipitation processes in a diffusive transport regime. The replacement process of celestine by strontianite can be easily identified by 3D Raman tomographic imaging, overcoming the more challenging identification by X-ray tomography [74]. The evolution of the pore architecture was captured at a voxel resolution of 0.5 µm³ which could even be improved to 0.3 µm³ with an objective of 1.4 numerical aperture. Such an experiment can be carried out at elevated temperatures using a heating stage. Moreover, the presented chip design can be used to image the gas phase in variably saturated reactive porous media or for monitoring two phase flow systems [75]. Future work will focus on the development of microfluidic experiments combined with an automated image analysis and pipeline to a flow solver for the derivation of upscaled transport parameters using neural network.