Surface‐Volume Scaling Controlled by Dissolution Regimes in a Multiphase Flow Environment

Fluid‐rock dissolution occurs ubiquitously in geological systems. Surface‐volume scaling is central to predicting overall dissolution rate R involved in modeling dissolution processes. Previous works focused on single‐phase environments but overlooked the multiphase‐flow effect. Here, through limestone‐based microfluidics experiments, we establish a fundamental link between dissolution regimes and scaling laws. In regime I (uniform), the scaling is consistent with classic law, and a satisfactory prediction of R can be obtained. However, the scaling for regime II (localized) deviates significantly from classic law. The underlying mechanism is that the reaction‐induced gas phase forms a layer, acting as a barrier that hinders contact between the acid and rock. Consequently, the error between measurement and prediction continuously amplifies as dissolution proceeds; the predictability is poor. We propose a theoretical model that describes the regime transition, exhibiting excellent agreement with experimental results. This work offers guidance on the usage of scaling law in multiphase flow environments.

due to reaction.Continuum-scale transport models demonstrated the potential in handling dissolution problems at large-scale application scenarios (Esteves et al., 2022;Steefel & Maher, 2009;Steefel et al., 2015;Wen & Li, 2017).Despite this, there is a significant disparity of several orders of magnitude between laboratory measurements and predictions made by reactive transport models regarding reaction rate (Le Traon et al., 2021;Liu et al., 2022).Previous investigations have recognized that the key limitation comes from the parameters (e.g., reaction rate constant, specific surface area, etc.) used in the continuum-scale model (Beckingham et al., 2016;Deng, Molins, et al., 2018;L. Li et al., 2008;Qin & Beckingham, 2021).These parameters require the knowledge of the pore-scale reactive transport process and are difficult to be determined by continuum-scale experiments solely.Thus, a comprehensive understanding of the pore-scale reactive transport mechanism is crucial in providing a more accurate estimation of the parameters used in the large-scale model and further enhancing the predictivity of the continuum-scale model (Seigneur et al., 2019;Steefel et al., 2015;Yoon et al., 2012).
Accurate estimation of the reaction rates is a prerequisite for the successful simulation of continuum-scale reactive transport.Reactive transport code, such as CrunchFlow ®, estimates the reaction rate R with where A is the reactive surface area, k is the reaction rate constant, and f e is a function of Gibbs free energy or the thermodynamic driving force, which is related to the concentration of reactants (Beckingham et al., 2016;Lasaga, 1984;Ling et al., 2022).Since the reaction rate constant and thermodynamic properties of mineral dissolution reactions are generally well understood, the principal uncertainties in modeling mineral reaction rates come from the estimation of mineral reactive surface area A (Beckingham et al., 2017;Qin & Beckingham, 2021).This reactive surface area A is usually estimated using a classic relationship: where Φ is the porosity and Φ m is the individual mineral volume fraction (Ling et al., 2022;Qin & Beckingham, 2021;Seigneur et al., 2019;Steefel & Lichtner, 1998).The superscript "i" refers to the corresponding quantity evaluated at t = 0.However, this equation is derived under the assumption that mineral grain is dissolved uniformly and the predictability of the scaling law under complex systems needs to be further studied (Noiriel & Daval, 2017;Seigneur et al., 2019;Steefel & Lichtner, 1998).
Reactive transport models mainly deal with the dissolution process under the single-phase flow condition (Deng et al., 2018;Molins et al., 2021;Noiriel & Soulaine, 2021;Steefel et al., 2015).Previous works have shown significant effects of multiphase flow on dissolution dynamics (P.Li et al., 2022;Soulaine et al., 2018;Thompson & Gdanskl, 1993;J. Xu & Balhoff, 2022).Soulaine et al. (2018) develop a pore-scale model to reproduce the generation of CO 2 bubbles during the dissolution process of calcite, revealing that bubbles prevent the emergence of the wormhole and limit the overall dissolution rate in porous media.P. Li et al. (2022) also simulated the dissolution processes under multiphase flow conditions, showing that the suppression of the overall dissolution rate is due to the occupation of the reactive surface area by non-reactive bubbles and the limited transport by the recirculation zones in troughs.Jiménez-Martínez et al. (2020) showed that multiphase flow significantly modifies the competition between flow and reaction, and the local dissolution rate tends to be uniform in the flow domain.These works highlight the key role of multiphase flow in dissolving confined geometry.However, how the pore-scale multiphase flow controls the interface-volume relationship (Equation 2) used in the continuum-scale reactive transport modeling is not well understood.
Here, we aim to elucidate how multiphase flow affects reactive transport via interface-volume relationship quantification.By conducting experimental observations and quantification, we correlate the dissolution regimes with surface-volume scaling laws and propose a theoretical model for the regime transition.Our findings reveal that the classic scaling's suitability is dependent on the dissolution regimes.Our research bridges the gap between pore-scale dissolution dynamics and the overall dissolution rate prediction in a multiphase flow system.

Materials and Methods
We perform multiphase flow-through dissolution experiments in the microfluidic-based visualization system (Figure S1 in Supporting Information S1).We fabricate soluble microfluidics using rectangular limestone slices (the thickness is 200 μm), as shown in Figures S2 and S3 in Supporting Information S1.To conduct a dissolution experiment, we saturate the microfluidics with deionized (DI) water.Then the HCl solution (pH ranging from 1 to 4, Codow Standard titration solution) is injected into the flow cell by a syringe pump (PhD Ultra 70-3007, Harvard Apparatus) with constant flow rates.To improve the differentiation between the liquid and gas phases (Figure S1d in Supporting Information S1), a fluorescent dye is added to the acid solution (0.1 mg/L, FLUORES-CENCE Super Fluor 488).As the acidic solution is in contact with the carbonate rock slices, dissolution occurs, and the CO 2 gas phase forms expressed as The multiphase flow-dissolution process is recorded by the fluorescent microscopy visualization techniques (ZEISS AXIO Observer) coupled with a high-speed CMOS camera (Axiocam 305 color) (Figure S1d in Supporting Information S1).Through image processing (Figure S4 in Supporting Information S1), we calculate the solid-fluid interfaces A(t), the fracture volume Φ(t), and the overall dissolution rate during dissolution.
We employ the classic Péclet number Pe and Damkohler number Da II to quantify the interplay among advection, diffusion, and reaction, namely Pe = vd 0 /D m , Da II = kd 0 /D m , where v is the average velocity, v = Q/[(2d 1 + d 2 )h 0 ]; d 0 = 200 μm is the width; D m is the molecular diffusion coefficient, D m = 2 × 10 −9 m 2 /s (Panga et al., 2005); k is the reaction rate constant of limestone depending on pH.The values of reaction rate constant are k = 1 × 10 −6 m/s, 1 × 10 −7 m/s, 1 × 10 −8 m/s and 1 × 10 −9 m/s for pH = 1, pH = 2, pH = 3 and pH = 4, respectively (Alarji et al., 2022).We conducted a total of 27 experiments, varying the flow rate Q from 0.24 to 24 μL/min, the pH from 1 to 4, and the geometry parameter λ (λ = d 2 /d 1 ) from 0.26 to 2.6.The experimental conditions and the related parameters are listed in Table S1 in Supporting Information S1.

Observation of Dissolution Regimes
From Figure 1, we observe two distinct dissolution regimes forming in the multiphase flow environment, namely regime I (Figure 1c) and regime II (Figure 1b).In regime I, the limestone is mainly dissolved in the perpendicular (vertical) channels (the initial aperture is d 1 ), and the aperture of these two channels increases significantly.
On the contrary, in regime II, the dissolution mainly occurs in the longitudinal (horizontal) channel (the initial aperture is d 2 ), leading to a localized dissolution channel in the flow direction.The typical dissolution processes of rock surfaces with injected time are shown in Figures 1d and 1e.The emergence and transformation of the two different dissolution regimes are influenced by the interplay between reaction and multiphase flow.
To understand the role of multiphase flow in dissolving rough fractures, we also present the fluid phase distribution, namely the HCl solution (fluid) and CO 2 (gas), in Figure 1a.For the lowest flow rate considered (Q = 0.24 μL/min, Figure 1a 9 -1a 12 ), regime I forms for all λ.In general, due to the lower flow resistance, the HCl solution is anticipated to flow mainly through the channels with a higher permeability.In the case of large λ (λ > 1) that the longitudinal channel has the larger permeability, HCl solution primarily flows through the longitudinal channel, leading to preferential dissolution.However, our experiments (Figures 1a 9 -1a 12 ) indicate that the dissolution primarily takes place in the perpendicular channels, regardless of the λ involved.This inconsistency highlights the substantial impact of multiphase flow on shaping the flow channels.In a multiphase environment, two factors dictate the infiltration of HCl solution into the channels: pressure drop and flow resistance from CO 2 bubbles.The CO 2 bubbles, resulting from the HCl and limestone reaction, create a non-wetting phase that occupies the longitudinal channels, leading to additional flow resistance and hindering HCl solution infiltration.
The above discussion indicates that in the case of the lowest flow rate (Q = 0.24 μL/min, Figure 1a 9 -1a 12 ), the flow resistance created by the CO 2 bubbles exceeds the pressure drop generated by the injected flow.This indicates that the multiphase flow resistance dominates the flow-reaction system, outweighing the advection from the flow rate.As a result, the longitudinal channel cannot be invaded by the HCl solution, while the perpendicular channels are able to be invaded, leading to the formation of regime I.As Q and λ values increase, the pressure As Q further increases up to Q = 24 μL/min (Figures 1a 1 -1a 4 ), the critical λ corresponding to regime transition decreases (a point between Figures 1a 2 and 1a 3 ).As the pressure drop rises, a greater flow resistance (smaller λ) is required to counteract this drop.In a word, our experimental findings and analysis demonstrate that multiphase flow plays a crucial role in shaping the formation and transformation of dissolution regimes by inducing extra flow resistance.

Scaling Laws Controlled by Dissolution Regimes
As indicated in the Introduction, surface-volume scaling is a central relationship to estimating the overall dissolution rate involved in a continuum-scale transport model.In this subsection, we will study how dissolution regimes (quantified in Figure 1) affect this scaling relationship.We calculate the evolution of A(t) and fracture volume Φ(t) and present the results in Figure 2. We also present the classic law obtained in single-phase flow system for comparison.Interestingly, for the first time, we find a fundamental link between dissolution regimes (Figure 1) and surface-volume scaling laws.In regime I, all of the experimental data is consistent with the classic scaling law with the same exponent 1/3: However, the surface-volume scaling in regime II deviates significantly from the classic law, with the exponent being 1/5: The detailed derivation of the classic scaling law (Equation 3) is given in Text S1 in Supporting Information S1.It is shown in Figure 1 that multiphase flow significantly affects the formation and transition of the dissolution patterns through the additional flow resistance induced by the CO 2 phase.In addition to the experimental observations and analyses, there are still unanswered mechanistic questions, such as why the scaling law in regime I in multiphase-flow systems is consistent with the classic law in single-phase flow system, and why the scaling law in regime II deviates significantly from the classic law.
We note that the underlying assumption of the classic scaling law under a single-phase flow environment is that rock samples have a high carbonate content and the mineral grains are uniformly dissolved (Ling et al., 2022;Noiriel & Daval, 2017;Seigneur et al., 2019;Steefel & Lichtner, 1998;Steefel et al., 2015), as shown in Figure 2.
In realistic application scenarios, the presence of complex multiphase hydrodynamics, mineralogy, and pore geometry weakens the uniformity assumption and leads the classic scaling law to deviation (Ling et al., 2022;Noiriel & Daval, 2017).Our study reveals that despite the significant influence of multiphase flow on the dissolution processes, the dissolution dynamics of regime I exhibit a similar morphology to that of the single-phase flow environment.The materials used in this work have a high carbonate content (92.87%, Figure S2 in Supporting Information S1), driving the experimental data to be consistent with the classic law.We acknowledge that bubbles are expected to appear in both perpendicular and horizontal channels.The mobility of the generated bubbles differs due to flow heterogeneity, that is, whether the channel is parallel or perpendicular to the incoming flow, The scaling relationship between A(t) and Φ(t) depends on dissolution regimes.The scaling relationship of regime I is consistent with classical law, namely −dA/dΦ ∼ Φ 1/3 (Text S1 in Supporting Information S1).The classic scaling law that assumes uniform dissolution of mineral grains is also shown for comparison purposes.In regime II, the scaling relationship deviates from the classical law, with the power index decreasing from 1/3 (classic law or regime I) to 1/5 (regime II).The slope of the A-ϕ scaling curves remains relatively constant within the range of porosity (0.1 ≤ ϕ ≤ 0.6) for both two regimes, with −0.98 ± 0.09 for regime I, and −0.53 ± 0.10 for regime II (Figure S11 in Supporting Information S1).Thus, the power-law scaling curves exhibit straight lines in linear scale plots.
and geometrical anisotropy, that is, differences in channel aperture.This disparity in mobility causes the bubbles to either be flushed away or merge, ultimately resulting in a deviation of the dissolution pattern toward the two observed regimes.Specifically, regime I is characterized by limestone dissolution primarily occurring in the perpendicular (vertical) channels, resulting in a significant increase in aperture size in these channels.In this case, the reaction-induced CO 2 bubbles would be flushed away perpendicular to the incoming flow.The effective area, that is, the acid-rock interface area, fluctuates as the dissolution proceeds, as evidenced by Figure 1e.However, regime II is characterized by dissolution predominantly occurring in the longitudinal (horizontal) channel.It causes the merged bubbles to shrink by strong incoming flow, eventually forming a CO 2 film as a barrier (Figure S7 in Supporting Information S1).Such an effect leads to the fact that the acid-rock interface area gradually decreases as dissolution progresses, as shown in Figure 1g.Therefore, the dissolution dynamics in regime II break the uniform dissolution assumption of the classic scaling law, leading to a significant deviation.

Model Predictability Controlled by Dissolution Regimes
As previously indicated, the dissolution regime significantly affects the scaling relationship of the reactive surface area and thus has a direct impact on the overall reaction rate.To investigate how the dissolution regimes control the applicability and predictability of the scaling laws used in the reactive transport model, we compare the predicted reaction rate, R t,pre , with the measured reaction rate R t,exp , for the two dissolution regimes.The measured rate, R t,exp , can be obtained via post-processing of the fluorescent images (Figure S4 in Supporting Information S1): where ΔV t is the dissolved volume of the limestone within the time interval Δt (Figure S4 in Supporting Information S1).
The predicted rate, R t,pre , can be obtained by substituting Equation 2into Equation 1: where k = 1 × 10 −6 m/s for pH = 1 (Alarji et al., 2022); f e can be considered as the inflow concentration C in under acidic pH, f e = C in (Noiriel et al., 2007;Panga et al., 2005); A(t) is the using the measured reactive surface area A in Figure 2.
To quantify the applicability and predictability of the scaling laws controlled by dissolution regimes, we introduce a correction factor C f and its trend of the variation k avg : where k avg is the average slope of C f curve to represent the trend of the variation of C f ; ΔC f,i is the variation of C f during the time interval Δt i ; n is the total number of time intervals.A larger value of k avg denotes the increases of C f with a higher growth rate; a lower value of k avg shows that C f fluctuates.
Figure 3a presents the variations of C f during dissolution, showing that there is one order of magnitude difference between the predicted and measured reaction rate.This discrepancy is expected because some factors cannot be considered using Equation 1, including the transport limitation depending on the scale (L.Li et al., 2008), precipitation-induced passivation (Noiriel et al., 2007), and surface area accessibility (Beckingham et al., 2017).
Here, we focus on the trend of the variation of C f affected by dissolution regimes during multiphase reactive flow.Qualitatively, the arrows stand for the trend of C f are shown in Figure 3.We observe that in dissolution regime I, the correction factor C f fluctuates around a stable value (Figures 3a 1 -3a 9 ), indicating that a satisfactory prediction for the reaction rate can be approached by shifting a constant.In other words, the classic model can be applied by introducing a constant correction factor C f .However, in dissolution regime II, such correction factor C f continuously increases as dissolution proceeds (Figures 3b 1 -3b 3 ), indicating that the difference between R t,pre and R t,exp amplifies.It thus implies that the classic model fails to estimate the dissolution rate in the regime II.
The dependence of the applicability and predictability of the scaling laws on the regimes can be further indicated by the contour plot of k avg (Figure 3b).From Figure 3b, we observe that as the dissolution regime shifts from I to II, the average slope increases, implying that C f increases with a higher growth rate in regime II.Consistent with Figure 3a, it evidences that the contour plot of k avg overlaps with the phase diagram of dissolution regimes in Figure 1.As dissolution regime I occurs, the dissolution rate can be predicted by modifying the classic reaction rate model with a constant C f due to its fluctuation with the normalized effective area (Figure 1f and Figure S8a in Supporting Information S1).However, in regime II, C f cannot reach a stable plateau, and the error of reaction rate predicted by the model (Equation 1) amplifies as the dissolution processes because the dissolution rate is suppressed by the barrier effect (Figure 1g and Figure S8b in Supporting Information S1).Although this shifting method using the correction factor C f is practiced in some previous studies (Beckingham et al., 2016;Hyman et al., 2022), our work quantifies its predictability controlled by dissolution regimes in a multiphase flow environment for the first time.

Theoretical Model for the Regime Transition
As previously indicated in Section 3.1, regime I is distinguished by the dissolution of limestone in channels perpendicular to the flow direction, whereas in regime II, transport and dissolution take place primarily in the longitudinal channel.We propose a theoretical model for the regime transition to generalize our experimental findings under various pH (or Da II ) conditions.The regime to which the dissolution dynamic belongs is determined by whether the bubble in the longitudinal channel moves or not.Therefore, we consider the critical state when a bubble forms in the longitudinal channel (Figure S9 in Supporting Information S1).Under such a state, the fluid flow in this channel is controlled by the interplay between driving force and resistance.For the multiphase flow in wetting media (imbibition), the driving force is from two components (Hu et al., 2018(Hu et al., , 2019)), including the viscous pressure P vis along the characteristic length δx and capillary pressure P cap on the water-CO 2 interface.The resistance P gas is from CO 2 bubbles generated during a unit time period Δt.The fluid flow behav ior in the longitudinal channel can be described as: where Q in and A in are the flow rate and area of the inlet; K is the effective permeability of the channel; μ is viscosity and L is the characteristic length scale in the longitudinal direction.
Da II = kd 0 /D m and λ = d 2 /d 1 into Equation 8, we establish a theoretical model for the regime transition (Text S2 in Supporting Information S1): where  CO 2 is the molar mass of CO 2 ;  CO 2 is the density of CO 2 ; β is the capillary coefficient; σ is the interfacial tension; R is the ideal gas constant; T denotes the temperature; α is an attenuation factor that accounts for the wall roughness, friction or other flow rate loss effects.This attenuation factor α is the only parameter that needs to be calibrated in the theoretical model.Here, α is first calibrated as 0.004 using the experimental data at pH = 1 (Figure 4a) and then is used to predict the dissolution regimes for pH = 2, 3, and 4. As shown in Figures 4b-4d, the model prediction agrees very well with the experiments, indicating that the calibrated value of α is reasonable.
The theoretical model (Equation 9a and 9b) provides a generality extension of the regime boundary as a function of geometry (λ), capillarity (σ), and reactivity (Ψ or k).We validate this theoretical model using experimental patterns under various pH conditions (Figure 4 and Figure S10 in Supporting Information S1).All the parameters are kept the same as those used for the experiments are listed in the Table S2 in Supporting Information S1.As shown in Figure 4, the regime boundary predicted by the model agrees very well with the experimental results across different pH experiments.The model accurately depicts the shift in regimes from high reactivity to low reactivity.Furthermore, as the pH increases (from high to low reactivity), the classic dissolution regime (regime I) expands in area.This expansion correlates with observable physical phenomena and suggests that as reactivity diminishes, the regime toward a uniform pattern, thereby reinforcing the applicability of the classic surface-volume law.In summary, Equation 9a and 9b provides a generalized boundary for the dissolution regimes across various dissolution scenarios.It also provides a theoretical method to identify the predictability of the classic surface-volume law in a multiphase flow environment.

Conclusions
We establish a fundamental connection between the regimes and surface-volume scaling laws.If the dissolution exhibits a uniform pattern (regime I), a long-term dissolution rate can be predicted.However, if the dissolution exhibits the localized regime (regime II), the scaling deviates significantly from the classic law.We identify a pore-scale barrier mechanism that suppresses the overall dissolution rate and leads to this deviation.We propose a theoretical model that provides a generalized boundary for the dissolution regimes, and it further provide a theoretical method to identify the predictability of the classic surface-volume law.This work improves our understanding of how dissolution regimes control the scaling laws and their applicability and predictability in a multiphase-flow environment.Our theoretical model is a significant contribution to using the surface-volume scaling laws in reactive transport modeling.It has practical implications for geologic carbon sequestration and acid stimulation in carbonate reservoirs.In such applications, it is crucial to accurately estimate the overall dissolution rate under multiphase flow conditions.The multiphase flow barrier effect in regime II (localized) would suppress the increase of permeability and thus reduce the CO 2 leakage risk and the gas/oil recovery.This barrier effect should be considered to improve the precision of the reactive flow modeling at the field scales.S2 in Supporting Information S1.

Figure 1 .
Figure 1.Experimental phase diagram of dissolution regimes for pH = 1.(a) Dissolution morphologies at the time of Φ = 2Φ i under three flow rates (Q) and four geometry parameters (λ = d 2 /d 1 ) conditions.The Pe number (or the flow rate Q) increases from the bottom to the top.The initial solid-fluid interfaces are denoted as black dashed lines.Hydrochloric acid (HCl) solution, carbonate slice, and CO 2 phase (gas) are shown as green, dark, and light brown, respectively.Two distinct dissolution regimes are observed, denoted as regime I (c) and regime II (b).In regime I (uniform), dissolution mainly occurs in the two perpendicular (vertical) channels (c), whereas the dissolution of limestone occurs along the longitudinal (horizontal) channel in regime II (b).(d, e) Typical dissolution processes with injected time for the two regimes (Movies S1 and S2).The evolution of the rock surfaces with the injected volume for all experimental conditions is shown in Figure S5 in Supporting Information S1.The variations of the normalized total surface area (fluid-rock interface) and the normalized effective area (acid-rock interface) are presented in (f) and (g) for regime I (d) and regime II (e).The reproducibility of the experimental results is shown in Figure S6 in Supporting Information S1.The experimental phase diagrams of dissolution regimes for pH = 2, pH = 3, and pH = 4 are shown in Figure S10 in Supporting Information S1.

Figure 2 .
Figure2.The fundamental link between dissolution regimes and scaling laws.The scaling relationship between A(t) and Φ(t) depends on dissolution regimes.The scaling relationship of regime I is consistent with classical law, namely −dA/dΦ ∼ Φ 1/3 (Text S1 in Supporting Information S1).The classic scaling law that assumes uniform dissolution of mineral grains is also shown for comparison purposes.In regime II, the scaling relationship deviates from the classical law, with the power index decreasing from 1/3 (classic law or regime I) to 1/5 (regime II).The slope of the A-ϕ scaling curves remains relatively constant within the range of porosity (0.1 ≤ ϕ ≤ 0.6) for both two regimes, with −0.98 ± 0.09 for regime I, and −0.53 ± 0.10 for regime II (FigureS11in Supporting Information S1).Thus, the power-law scaling curves exhibit straight lines in linear scale plots.

Figure 3 .
Figure 3. Model predictability controlled by dissolution regimes.(a) Variations of correction factor C f for all experimental conditions.C f is ratio of the predicted to the measured reaction rates.In regime I, the correction factor C f fluctuates around a stable value.In regime II, C f continuously increases, and the difference between the predicted and measured reaction rates amplifies.The classic model fails to estimate the dissolution rate in this regime.(b) Contour plot of the average slope of C f (k avg ).k avg is used to quantify the trend of C f affected by dissolution regimes.As the dissolution regime shifts from I to II, C f increases, indicating that the predictability decreases.

Figure 4 .
Figure 4.The theoretical model for the transition of the dissolution regimes in a multiphase flow environment.From the comparison between the predicted regime boundaries (solid line) and the experimental results (squares and triangles) with pH = 1 (a), pH = 2 (b), pH = 3 (c), and pH = 4 (d), the theoretical model (Equation 9a and 9b) can well capture the transition of the regime transition affected by the geometry (λ) and reactivity (Ψ or k).The parameters used in the theoretical prediction are listed in TableS2in Supporting Information S1.