1. Introduction
The increasing emission of carbon dioxide (CO2) into the atmosphere has been identified as a primary contributor to the greenhouse effect, driving global climate change (Celia et al. Reference Celia, Bachu, Nordbotten and Bandilla2015; De Silva, Ranjith & Perera Reference De Silva, Ranjith and Perera2015). To mitigate this issue, several strategies have been implemented to reduce carbon emissions, with carbon capture and storage through geological CO2 sequestration emerging as an efficient solution (He et al. Reference He, Xu, Ji and Jiang2022). Conventional geological sequestration methods include storing CO2 in deep saline aquifers, depleted oil and gas reservoirs, and unmineable coal seams (Akindipe, Saraji & Piri Reference Akindipe, Saraji and Piri2022). Of these, saline aquifers are considered the most suitable due to their wide distribution, large storage capacity, and proximity to CO2 emission sources (Cui et al. Reference Cui, Hu, Ning, Jiang and Wang2023). During storage in saline aquifers, supercritical CO2 is injected into the reservoir, displacing the pre-existing brine. Following injection, CO2 undergoes various trapping mechanisms, such as physical trapping (structural and capillary) and chemical trapping (dissolution and mineralisation), which facilitate its secure storage within the aquifer (Izadpanahi et al. Reference Izadpanahi, Blunt, Kumar, Ali, Gaeta Tassinari and Sampaio2024). The success of the injection process is crucial to ensure that the subsequent trapping mechanisms are effective, making efficient CO2 injection a critical component of the sequestration process. Injectivity, defined as the rate at which CO2 can be injected, is a key metric for evaluating a wellbore’s capability for CO2 sequestration (Hajiabadi et al. Reference Hajiabadi, Bedrikovetsky, Borazjani and Mahani2021). High injectivity is essential for maintaining sequestration efficiency and is influenced by reservoir properties such as porosity and permeability (He et al. Reference He, Xu, Ji and Jiang2022). However, salt precipitation often occurs in the near-wellbore regions during CO2 injection, significantly reducing reservoir injectivity (Baumann, Henninges & De Lucia Reference Baumann, Henninges and De Lucia2014; Grude, Landrø & Dvorkin Reference Grude, Landrø and Dvorkin2014). Addressing this challenge requires identifying the key factors driving salt precipitation and developing effective strategies to mitigate its impact.
Salt precipitation during CO2 injection is a complex multiphysics phenomenon in porous media (Li, Steefel & Jun Reference Li, Steefel and Jun2017; Nooraiepour et al. Reference Nooraiepour, Fazeli, Miri and Hellevang2018). As CO2 is injected into saline aquifers, it displaces the resident brine and induces its evaporation due to the dry nature of the injected gas. This evaporation increases the salt concentration in the brine. Once the concentration exceeds the saturation threshold, salt begins to precipitate. The formation of salt can obstruct the porous structure of the reservoir, hindering both CO2 injection and brine migration. This process involves intricate mechanisms, including gas–brine multiphase flow, vapour mass transfer, phase changes, salt concentration evolution, and mineral crystallisation (Yang et al. Reference Yang, Lei, Wang, Xu, Chen and Luo2023).
To better understand the effects and mechanisms of salt precipitation, extensive studies have been conducted at both laboratory scales (Oh et al. Reference Oh, Kim, Han, Kim, Kim and Park2013; Tang et al. Reference Tang, Yang, Du and Zeng2014; Zhang et al. Reference Zhang, Kang, Selvadurai and You2019) and pore scales (Kim, Sell & Sinton Reference Kim, Sell and Sinton2013; Norouzi Rad, Shokri & Sahimi Reference Norouzi Rad, Shokri and Sahimi2013). Among these, pore-scale studies have proven particularly effective in exploring the multiphysics mechanisms underlying salt precipitation. Their strength lies in their ability to probe the intricate depths of porous media, enabling detailed observation of the morphological evolution of salt precipitation and gas–brine phase transitions. With advancements in observational techniques, a variety of pore-scale experiments have been performed to characterise salt precipitation. These include two-dimensional microfluidics (Kim et al. Reference Kim, Sell and Sinton2013; Miri et al. Reference Miri, van Noort, Aagaard and Hellevang2015; Nooraiepour et al. Reference Nooraiepour, Fazeli, Miri and Hellevang2018; He et al. Reference He, Jiang and Xu2019, Reference He, Xu, Ji and Jiang2022; Ho & Tsai Reference Ho and Tsai2020) and three-dimensional X-ray micro-computed tomography (micro-CT) (Roels, Ott & Zitha Reference Roels, Ott and Zitha2014; Ott, Roels & de Kloe Reference Ott, Roels and de Kloe2015; Rad et al. Reference Norouzi Rad, Shokri, Keshmiri and Withers2015; Akindipe et al. Reference Akindipe, Saraji and Piri2021, Reference Akindipe, Saraji and Piri2022). These investigations have provided valuable insights into the mechanisms of salt precipitation in porous media. For instance, Akindipe, Saraji & Piri (Reference Akindipe, Saraji and Piri2021) used micro-CT to observe salt precipitation during CO2 injection into saline aquifers. Their findings identified three stages in the precipitation process: advection-dominated flow, transition, and diffusion-dominated flow. During the advection-dominated stage, CO2 displaced brine. In the transition stage, salt concentration increased at the evaporation front. Finally, in the diffusion-dominated stage, salt solutes accumulated at the gas–brine interface, leading to precipitation. Similarly, Ho & Tsai (Reference Ho and Tsai2020) observed a drainage–evaporation–precipitation sequence in microfluidic experiments. These studies shed light on the interplay of gas–brine flow, evaporation and crystallisation processes, offering theoretical insights to mitigate salt precipitation. Pore-scale investigations have also provided detailed information on the morphology and location of salt precipitation. Miri et al. (Reference Miri, van Noort, Aagaard and Hellevang2015) identified two primary forms of crystallisation: large single crystals in the bulk brine, and aggregate growth within the gas phase, driven by brine absorption, commonly referred to as salt self-enhancement. He, Jiang & Xu (Reference He, Jiang and Xu2019) found that aggregated precipitation occurred in hydrophilic porous media due to brine film backflow. Using micro-CT imaging, Rad et al. (Reference Norouzi Rad, Shokri, Keshmiri and Withers2015) demonstrated that salt precipitation patterns are influenced by the grain size of the porous media, highlighting the importance of capillary-driven brine flow and preferential evaporation sites. These comprehensive studies have significantly advanced understanding of salt precipitation patterns and the factors influencing them, including pore structure and brine distribution.
Most pore-scale studies have focused on porous media matrices. However, in reservoirs used for CO2 sequestration, fractures are commonly formed through either natural evolution or engineering drilling. The disparity in size between fractures and the pore structures in the matrix creates notable differences in gas–water flow and evaporation characteristics, which significantly impact salt precipitation. He et al. (Reference He, Xu, Ji and Jiang2022) investigated salt precipitation in fractured porous media. Their study revealed that capillary forces predominantly direct CO2 flow through fractures, displacing brine, while brine within the matrix undergoes evaporation rather than displacement. In hydrophilic porous media, brine backflow leads to salt accumulation within fractures, obstructing flow channels. Similar observations were reported by Miri et al. (Reference Miri, van Noort, Aagaard and Hellevang2015). These findings provide a preliminary understanding of salt precipitation behaviour in fractured porous structures. However, salt precipitation induced by CO2 injection is governed by complex multiphysics mechanisms that differ between fractures and pores due to their size variations. Traditional experimental methods are insufficient to capture these complexities, highlighting the need for numerical simulations to gain a deeper understanding of the underlying processes.
While pore-scale experiments have provided valuable insights into the mechanisms of salt precipitation, their scope is limited to observable phenomena. Critical aspects, such as the evolution of salt concentration and its effects on precipitation kinetics, remain insufficiently understood. To achieve a more comprehensive and quantitative understanding of these phenomena, pore-scale numerical simulations are essential. These simulations enable detailed investigation of reaction and transport mechanisms within porous media. With advancements in computational techniques, various numerical models have been developed to study salt precipitation processes, including direct numerical simulations (Deng et al. Reference Deng, Tournassat, Molins, Claret and Steefel2021), pore network models (Dashtian, Shokri & Sahimi Reference Dashtian, Shokri and Sahimi2018), and the lattice Boltzmann (LB) method (Prasianakis et al. Reference Prasianakis, Curti, Kosakowski, Poonoosamy and Churakov2017). Among these, the LB method stands out for pore-scale numerical simulation (Xu et al. Reference Xu, Dai, Yang, Liu and Shi2022; Lei et al. Reference Lei2023; Fei, Derome & Carmeliet Reference Fei, Derome and Carmeliet2024) for its superior parallel efficiency (Latt et al. Reference Latt2021) and its ability to integrate multiphysics models (Yang et al. Reference Yang, Xu, Kou, Wang, Lei, Wang, Li and Luo2024). Consequently, it has been widely adopted for simulating salt precipitation processes (Yoon et al. Reference Yoon, Valocchi, Werth and Dewers2012; Chen et al. Reference Chen, Kang, Carey and Tao2014; Prasianakis et al. Reference Prasianakis, Curti, Kosakowski, Poonoosamy and Churakov2017; Fazeli et al. Reference Fazeli, Masoudi, Patel, Aagaard and Hellevang2020; Masoudi et al. Reference Masoudi, Fazeli, Miri and Hellevang2021; Nooraiepour, Masoudi & Hellevang Reference Nooraiepour, Masoudi and Hellevang2021). For instance, Fazeli et al. (Reference Fazeli, Masoudi, Patel, Aagaard and Hellevang2020) utilised the LB method to introduce a probabilistic model based on classical nucleation theory (CNT). This model explored the dynamics of salt precipitation under varying supersaturation conditions, growth kinetics and brine flow rates. Masoudi et al. (Reference Masoudi, Fazeli, Miri and Hellevang2021) extended this approach to three-dimensional micro-CT images, enabling detailed analysis of how grain size influences salt precipitation characteristics. However, most existing studies focus on single-phase brine flow and do not adequately account for the complexities arising from CO2-induced brine evaporation and the associated changes in salt concentration.
To address this gap, our recent work (Yang et al. Reference Yang, Lei, Wang, Xu, Chen and Luo2023) introduced an LB numerical model to simulate salt precipitation during brine evaporation. This model integrates key processes, including gas–brine multiphase flow, brine evaporation, salt concentration evolution, and the formation of salt precipitation. We extended the application of this model to simulate salt precipitation in a porous matrix, where the numerical results closely aligned with experimental observations. As highlighted earlier, the dynamics of a matrix containing fractures exhibits distinct characteristics compared to that of a homogeneous matrix without fractures. However, to date, pore-scale numerical studies have not adequately addressed the interaction of physical processes during CO2 injection into fractured matrices. This represents a critical gap that requires further investigation to enhance understanding and optimise CO2 sequestration strategies.
In this study, pore-scale numerical investigations were conducted to examine salt precipitation in a porous matrix with fractures during CO2 injection, building on previous microfluidic experiments (He et al. Reference He, Xu, Ji and Jiang2022). The study comprehensively considered multiple physical processes, including gas–brine multiphase flow, brine evaporation, salt concentration evolution, and salt precipitation growth. The numerical results were rigorously compared with experimental data to elucidate the mechanisms underlying salt precipitation. Additionally, simulations of the CO2 injection process under various reservoir conditions were performed to evaluate the impact of salt precipitation on injection performance.
2. Governing equations and numerical models
2.1. Description of physical problem
The physical processes investigated in this study are illustrated schematically in figure 1. As shown in figure 1(a), the numerical simulation replicates the microfluidic experimental process conducted byHe et al. (Reference He, Xu, Ji and Jiang2022). The computational domain comprises a porous media matrix and an underlying fracture. The total domain size is
$5\,\text{mm}\times10\,\text{mm}$
, with fracture width
$L=400$
μm, selected to align with the geometric parameters of the microfluidic chip in He et al. (Reference He, Xu, Ji and Jiang2022), where the fracture width was specified as 500 µm. The porous matrix consists of particles measuring 300 μm in diameter, with overall porosity 0.48. The narrowest pore-throat width,
$l=85\,\unicode{x03BC} \text{m}$
, corresponds to the flow channel dimensions of the microfluidic chip. Initially, the porous matrix was saturated with NaCl brine. Dry CO2 gas was injected through the fracture inlet to induce brine evaporation and subsequent salt precipitation in the matrix. Figure 1(b) illustrates schematically the mechanisms of brine evaporation and salt precipitation. When dry CO2 is injected, the mass fraction (or partial pressure) of water vapour in the gas phase decreases, triggering brine evaporation. This evaporation leads to an increase in salt concentration near the gas–brine interface as salt species undergo advection and diffusion. This concentration increase initiates the nucleation and growth of salt crystals. The resulting salt precipitation alters the pore structure within the matrix, thereby influencing subsequent physical processes.
(a) Schematic description of the physical problem. The computational domain is designed based on microfluidic experiments. (He et al. Reference He, Xu, Ji and Jiang2022). Black, blue and red colours represent matrix grain, gas and brine, respectively. (b) Schematic representation of mechanisms of brine evaporation and salt precipitation during CO2 injection.
A pore-scale numerical simulation was performed using LB models to elucidate the physical mechanisms underlying the experimentally observed phenomena. The multiphysics mechanisms depicted in figure 1(b) were incorporated comprehensively, including gas–brine multiphase flow, brine evaporation with vapour mass transfer, salt species mass transfer in brine, and the nucleation and growth of salt precipitation. Individual numerical submodels were developed to describe each of these processes, as detailed in the next subsection. These submodels were fully coupled to simulate the entire CO2 injection process, ensuring an integrated analysis of all relevant interactions. Mechanisms such as the evolution of salt concentration and evaporation rates were analysed in detail. Additionally, a regime analysis of salt precipitation patterns was conducted under conditions both replicating and extending beyond the experimental set-up.
2.2. Numerical assumptions and governing equations
To mathematically articulate the physical processes, a set of assumptions was established to construct the governing equations. (i) Temperature changes resulting from CO2 injection were considered negligible, and heat transfer was excluded. The temperature was maintained constant at the experimental benchmark of
$50\,^\circ\text{C}$
(He et al. Reference He, Xu, Ji and Jiang2022). (ii) Brine was treated as a diluted solution, with Fick’s law employed to describe the mass transfer of salt species (Yang et al. Reference Yang, Lei, Wang, Xu, Chen and Luo2023). (iii) The mass transfer of salt within the brine was assumed to have a negligible effect on fluid properties, and the influence of gravity on fluid dynamics was considered minimal. (iv) The gas phase was conceptualised as an ideal gas, with the mass transfer of brine vapour in the gas phase contributing to brine evaporation (Fei et al. Reference Fei, Qin, Zhao, Derome and Carmeliet2022). (v) The formation of salt precipitates was governed by CNT, while the movement of these precipitates within the fluid flow was not considered (Deng, Poonoosamy & Molins Reference Deng, Poonoosamy and Molins2022). (vi) The dissolution of salt precipitates in the brine was not accounted for. Based on these assumptions, the governing equations for gas–brine multiphase flow, brine vapour mass transfer in the gas phase, and salt species mass transfer in the brine phase are presented as follows:
\begin{align} &\qquad\qquad\qquad\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\right)=0, \end{align}
\begin{align} \frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}+\nabla \cdot \left(\rho \mathbf{uu}\right)&=-\nabla p+\nabla \cdot \left(\rho \upsilon \left(\nabla \mathbf{u}+\left(\nabla \mathbf{u}\right)^{T}\right)\right)+\nabla \left(\rho \left(\xi -\frac{2}{3}\upsilon \right)\left(\nabla \cdot \mathbf{u}\right)\right), \end{align}
\begin{align} &\frac{\partial \rho Y^{{b}}}{\partial t}+\nabla \cdot \left(\rho Y^{{b}}\mathbf{u}\right)=-\nabla \left(\rho D^{{b}}\,\nabla Y^{{b}}\right)\quad\text{in }\Omega _{{g}}, \end{align}
\begin{align} &\qquad\frac{\partial C}{\partial t}+\nabla \cdot \left(C\mathbf{u}\right)=\nabla \cdot \left(D_{{s}}\,\nabla C\right)\quad\text{in }\Omega _{{b}}, \end{align}
where
$\rho$
denotes the fluid density,
$t$
denotes the time,
$\mathbf{u}$
represents the fluid velocity, and
$p$
indicates the pressure involving the capillary force determined by surface tension
$\kappa$
. Also,
$\upsilon$
and
$\xi$
are kinematic and bulk viscosity. Here,
$Y^{{b}}$
denotes the mass fraction of brine vapour in the gas-phase domain
$\Omega _{{g}}$
. It is defined as
$Y^{{b}}=\rho ^{{b}}/\rho$
, where
$\rho ^{{b}}$
indicates the density of brine vapour component. Also,
$D^{{b}}$
denotes the diffusivity of brine vapour in the gas phase,
$C$
represents the salt concentration in the brine phase domain
$\Omega_{{b}}$
, and
$D_{{s}}$
is the diffusivity of salt. At the gas–brine interface,
$Y^{{b}}$
should be a constant according to the saturation vapour pressure of brine, and the evaporation mass flux
$m'$
at the interface is calculated as (Sugimoto et al. Reference Sugimoto, Sawada, Kaneda and Suga2021)
\begin{equation} m'=\frac{\rho D^{{b}}}{1-Y^{{b}}}\nabla Y^{{b}}\cdot \mathbf{n}\quad\text{in }\Omega_{{g}}, \end{equation}
where
$\mathbf{n}$
denotes the vector normal to the interface. Because salt exists in the brine phase only at the gas–brine interface, the salt mass flux should be zero (Maes & Soulaine Reference Maes and Soulaine2018):
\begin{equation} C\left(\mathbf{u}-\mathbf{w}\right)-D_{{s}}\,\nabla C=\mathbf{0}\quad\text{in }\Omega _{{b}}, \end{equation}
where
$\mathbf{w}$
denotes the velocity of the gas–brine interface.
When the salt concentration exceeds the solubility limit, nucleation and growth of salt precipitates occur. The nucleation rate
$J_{{n}}$
can be calculated using CNT as (Li et al. Reference Li, Fernandez-Martinez, Lee, Waychunas and Jun2014)
\begin{equation} J_{{n}}=J_{0}\exp \left(-\frac{\Delta G_{{n}}}{k_{{B}}T_{{n}}}\right)=J_{0}\exp \left(-\frac{\beta \left(V_{{m}}/N_{{A}}\right)^{2}\gamma ^{3}}{k_{{B}}^{3}T_{{n}}^{3}\left(\ln\Phi \right)^{2}}\right), \end{equation}
where
$J_{0}$
denotes the pre-exponential factor, related to the frequency of collision,
$\Delta G_{{n}}$
is the free-energy barrier nucleus formation,
$k_{{B}}=1.38\times 10^{-23}\,\textrm{J}\,\textrm{K}^{-1}$
is the Boltzmann constant,
$T_{{n}}$
is the temperature,
$\beta$
denotes a geometric factor based on the shape of the nucleus,
$V_{{m}}$
is the molar volume of the salt,
$N_{A}=6.02\times 10^{23}\text{ mol}^{-1}$
is the Avogadro’s number,
$\gamma$
represents the interfacial free energy, and
$\ln\phi$
is the supersaturation, defined as
\begin{equation} \ln\Phi =\ln \left(C^{2}/K_{{sp}}\right), \end{equation}
where
$K_{{sp}}$
is the solubility product of salt. After nucleation, salt precipitation begins to grow at the existing nucleation sites at a rate
$r$
given by (Starchenko Reference Starchenko2022)
\begin{equation} r=A_{{s}}k\left(\frac{C^{2}}{K_{{sp}}}-1\right)=-A_{{s}}D_{{s}}\frac{\partial C}{\partial \mathbf{n}}, \end{equation}
where
$A_{{s}}$
represents the surface area of salt precipitation, and
$k$
is the kinetic coefficient. Equations (2.1)–(2.9) describe the multiphysical processes during CO2 injection, including multiphase flow with phase changes, salt concentration evolution and salt precipitation. By introducing the characteristic quantities – length
$L$
, characteristic velocity
$U$
, characteristic density
$\rho _{{ch}}$
, characteristic mass fraction
$Y_{{ch}}^{{b}}$
, and characteristic concentration
$C_{{ch}}$
– dimensionless parameters, marked with asterisks, can be derived by dimensionless analysis as
\begin{equation} \begin{array}{l} x^{*}=\displaystyle\frac{x}{L},\quad y^{*}=\frac{y}{L},\quad z^{*}=\frac{z}{L},\quad t^{*}=\frac{tU}{L},\quad \mathbf{u}^{*}=\frac{\mathbf{u}}{U}, \\ \displaystyle\rho ^{*}=\frac{\rho }{\rho _{{ch}}},\quad C^{*}=\frac{C}{C_{{ch}}},\quad p^{*}=\frac{pL}{\kappa },\quad Y^{{b*}}=\frac{Y^{{b}}-Y_{{eq}}^{{b}}}{Y_{{ch}}^{{b}}-Y_{{eq}}^{{b}}},\\ \displaystyle Ca=\frac{\rho _{{ch}}\upsilon U}{\kappa },\quad Pe^{{b}}=\frac{UL}{D^{{b}}},\quad Da_{{s}}=\frac{kC_{{ch}}L}{K_{{sp}}D_{{s}}}, \end{array} \end{equation}
where
$Y_{{eq}}^{{b}}$
denotes the brine vapour mass fraction at the brine–gas interface, corresponding to the saturated vapour pressure. Here,
$Y_{{ch}}^{{b}}$
was set as the brine vapour mass fraction at the inlet
$Y_{{in}}^{{b}}$
, whereas the characteristic salt concentration
$C_{{ch}}$
was set as
$C_{{ch}}=C_{0}^{2}/K_{{sp}}^{0.5}$
to represent its role in determining the precipitation rate. The dimensionless version of the governing equations can be found in (S1)–(S4) of the supplementary material. The key characteristic numbers include the capillary number
$Ca$
, Péclet number
$Pe$
, and Damköhler number
$Da$
. The Péclet number in the gas phase
$Pe^{{b}}$
is related to the transport of brine vapour. This characterises the competition between advection and diffusion of brine vapour, directly influencing the evaporation rate. Additionally, the Damköhler number in the brine phase
$Da_{{s}}$
reflects the competition between salt precipitation and salt species mass transfer. These characteristics provide a framework for analysing the relationships between the key physical processes involved in salt precipitation.
2.3. Numerical models
Pore-scale numerical models based on the LB framework were developed to solve the governing equations outlined above. To simulate gas–brine multiphase flow and brine evaporation described by (2.1)–(2.3) and (2.5), a multiphase pseudo-potential LB model was employed (Li et al. Reference Li, Luo, Kang, He, Chen and Liu2016). This model enables gas–brine multiphase flow with a high-density ratio of up to
$O(10^2)$
, and is capable of calculating the advection–diffusion process of vapour in the gas phase. The mass transfer of salt species within the brine was simulated using the continuum species transport (CST)-LB model (Yang et al. Reference Yang, Dai, Xu, Liu, Shi and Long2021) to solve (2.4) and (2.6). This model effectively handles complex interfacial mass transfer in multiphase systems. The nucleation of salt precipitates was calculated using CNT (Starchenko Reference Starchenko2022), and salt precipitation growth was achieved by dynamically updating the salt volume in each computational grid. The numerical models and physical parameters are briefly introduced in this subsection, with further details available in our previous study (Yang et al. Reference Yang, Lei, Wang, Xu, Chen and Luo2023). To the best of our knowledge, this numerical model is the first to comprehensively simulate complex salt precipitation in porous media, integrating processes such as brine evaporation, salt concentration increases, and the nucleation and growth of salt precipitates.
2.3.1. Pseudopotential model for gas–brine multiphase flow and brine evaporation
The multiphase LB model was utilised to simulate gas–brine multiphase flow and brine evaporation. Over the years, various multiphase LB models have been developed, including the colour-gradient model (Gunstensen et al. Reference Gunstensen, Rothman, Zaleski and Zanetti1991), pseudo-potential model (Shan & Chen Reference Shan and Chen1993; Li et al. Reference Li, Luo, Kang, He, Chen and Liu2016), free-energy model (Swift et al. Reference Swift, Orlandini, Osborn and Yeomans1996) and phase-field model (He, Chen & Zhang Reference He, Chen and Zhang1999). Among these, the pseudo-potential model is particularly notable for its simplicity in implementation and broad applicability. It effectively handles multicomponent multiphase flows with high-density ratios and phase change dynamics (Mukherjee, Berghout & Van den Akker Reference Mukherjee, Berghout and Van den Akker2019; Fei et al. Reference Fei, Qin, Wang, Huang, Wen, Zhao, Luo, Derome and Carmeliet2023; Qin et al. Reference Qin, Fei, Zhao, Kang, Derome and Carmeliet2023). As a result, a multicomponent pseudo-potential model was employed in this study. The LB equation for the
$\sigma$
th component (where
$\sigma=g$
for gas,
$\sigma=b$
for brine), using the multiple-relaxation-time (MRT) collision operator, is expressed as (Guo & Zheng 2008)
\begin{align} \mathbf{f}^{\sigma \ast }\left(\mathbf{x},t\right)=\mathbf{f}^{\sigma }\left(\mathbf{x},t\right)&-\mathbf{M}^{-1}\mathbf{\Lambda }^{\sigma }\left[\mathbf{m}^{\sigma }\left(\mathbf{x},t\right)-\mathbf{m}^{\sigma, \textrm{eq}}\left(\mathbf{x},t\right)\right]+\Delta t\cdot \mathbf{M}^{-1}\left(\mathbf{I}\mathbf{-}\frac{\mathbf{\Lambda }^{\sigma }}{2}\right)\mathbf{S}_{F}^{\sigma }\notag\\&\quad f_{\alpha }^{\sigma }\left(\mathbf{x}+\mathbf{e}_{\alpha }\Delta t,t+\Delta t\right)=f_{\alpha }^{\sigma \ast }\left(\mathbf{x},t\right) \end{align}
where
$\mathbf{f}^{\sigma }=[f_{0}^{\sigma },\ldots, f_{\alpha }^{\sigma },\ldots, f_{q-1}^{\sigma }]^{T}$
contains the density distribution function of the
$\sigma$
th component
$f_{\alpha }^{\sigma }$
at discrete velocity direction
$\alpha$
, position
$\mathbf{x}$
, and time
$t$
. In this work, the D3Q19 discrete velocity model is employed with
$q=19$
,
$\mathbf{m}^{\sigma }=\mathbf{M}\mathbf{f}^{\sigma }$
is the density distribution function moment obtained by transformation matrix
$\mathbf{M}$
,
$\mathbf{m}^{\sigma, \textrm{eq}}$
is the equilibrium moment,
$\mathbf{S}_{F}^{\sigma }$
is the source term to introduce the force
$\mathbf{F}^{\sigma }$
on the fluid,
$\mathbf{\Lambda }^{\sigma }$
is the diagonal relaxation matrix containing parameters related to the kinematic viscosity
$\upsilon ^{\sigma }$
and bulk viscosity
$\xi ^{\sigma }$
,
$\mathbf{e}_{\alpha }$
denotes the discrete velocity vector, and
$\Delta x$
and
$\Delta t$
are lattice space and time steps. The density and velocity of each component can be calculated as
\begin{equation} \rho ^{\sigma }=\sum f_{\alpha }^{\sigma },\quad \rho ^{\sigma }\mathbf{u}^{\sigma }=\sum f_{\alpha }^{\sigma }\mathbf{e}_{\alpha }+\frac{\Delta t}{2}\mathbf{F}^{\sigma }. \end{equation}
The total density and mixture velocity are then obtained as (Chai & Zhao Reference Chai and Zhao2012)
\begin{equation} \rho =\sum \rho ^{\sigma },\quad \mathbf{u}=\sum \rho ^{\sigma }\mathbf{u}^{\sigma }\Big/\sum \rho ^{\sigma }. \end{equation}
For the pseudo-potential model, phase separation and transitions were achieved by introducing interaction forces between the density distribution functions. In this study, both the intracomponent force
$\mathbf{F}_{\sigma \sigma }$
and the intercomponent force
$\mathbf{F}_{\sigma \overline{\sigma }}$
were considered as (Li et al. Reference Li, Luo, Kang, He, Chen and Liu2016)
\begin{align} \mathbf{F}_{\sigma \sigma }=-G_{\sigma \sigma }\,\psi ^{\sigma }\left(\mathbf{x}\right)\sum _{\alpha }w\left(\left| \mathbf{e}_{\alpha }\right| ^{2}\right)\psi ^{\sigma }\left(\mathbf{x}+\mathbf{e}_{\alpha }\right)\mathbf{e}_{\alpha }, \end{align}
\begin{align} \mathbf{F}_{\sigma \overline{\sigma }}=-G_{\sigma \overline{\sigma }}\,\varphi ^{\sigma }\left(\mathbf{x}\right)\sum _{\alpha }w\left(\left| \mathbf{e}_{\alpha }\right| ^{2}\right)\varphi ^{\overline{\sigma }}\left(\mathbf{x}+\mathbf{e}_{\alpha }\right)\mathbf{e}_{\alpha }, \end{align}
where
$w$
denotes the weight factor. The pseudo-potential
$\psi ^{\sigma }$
and
$\varphi ^{\sigma }$
were computed using the fluid density of each component as follows:
\begin{equation} \psi ^{\sigma }\left(\rho ^{\sigma }\right)=\sqrt{\frac{2}{G_{\sigma \sigma }}\left(p_{{EOS}}^{\sigma }-c_{s}^{2}\rho ^{\sigma }\right)},\quad \varphi ^{\sigma }\left(\rho ^{\sigma }\right)=1-\exp \left(-\rho ^{\sigma }/\rho _{0}^{\sigma }\right), \end{equation}
where
$G_{\sigma \sigma }=\text{sign}(p_{{EOS}}^{\sigma }-c_{s}^{2}\rho ^{\sigma })$
and
$c_{s}=\Delta x/\sqrt{3}\,\Delta t$
. Here,
$p_{{EOS}}^{\sigma }$
is calculated according to the equation of state (Yuan & Schaefer Reference Yuan and Schaefer2006), and
$\rho _{0}^{\sigma }$
is the numerical parameter selected according to the fluid density of each component. Employing the Chapman–Enskog analysis, it is feasible to derive the Navier–Stokes equations (2.1) and (2.2) from the LB equations, incorporating the total pressure (Li, Luo & Li Reference Li, Luo and Li2013):
\begin{equation} p=\sum _{\sigma }\left(p_{EOS}^{\sigma }+\frac{1}{2}G_{\sigma \overline{\sigma }}\varphi _{\sigma }\varphi _{\,\overline{\sigma }}\right). \end{equation}
The total pressure accounts for both the thermodynamic pressure and capillary forces. Surface tension
$\kappa$
is controlled by adjusting the parameter
$G_{\sigma \overline{\sigma }}$
, and its value can be validated through simulations based on the Laplace law. Similarly, the advection–diffusion equation (2.3) for brine vapour mass transfer can be derived from the LB equation (2.11) using Chapman–Enskog analysis (Guo et al. Reference Guo, Zhang, Lei and Galindo-Torres2022). This derivation inherently incorporates the phase transition process described in (2.5) within its framework (Yang et al. Reference Yang, Lei, Wang, Xu, Chen and Luo2023).
2.3.2. The CST-LB model for salt species mass transfer in brine
The mass transfer of salt species within the brine represents a typical advection–diffusion process in a multiphase system, as described by (2.4). To address this complexity, the CST–LB model was employed (Yang et al. Reference Yang, Dai, Xu, Liu, Shi and Long2021). This model has been validated rigorously for its effectiveness in simulating species mass transfer in multiphase systems. The CST-LB equation with the MRT collision operator is given by
\begin{equation} \begin{array}{c} \mathbf{g}^{\ast }\left(\mathbf{x},t\right)=\mathbf{g}\left(\mathbf{x},t\right)-\mathbf{N}^{-1}\left(\boldsymbol{\Lambda }_{{s}}\left(\mathbf{n}\left(\mathbf{x},t\right)-\mathbf{n}^{{eq}}\left(\mathbf{x},t\right)\right)+\Delta t\left(\mathbf{I }-\dfrac{\mathbf{\Lambda }_{{s}}}{2}\right)\mathbf{S}_{{CST}}\right),\\ g_{\alpha }\left(\mathbf{x}+\mathbf{e}_{\alpha }\,\Delta t,t+\Delta t\right)=g_{\alpha }^{\ast }\left(\mathbf{x},t\right), \end{array} \end{equation}
where
$\mathbf{g}=[g_{0},\ldots, g_{\alpha },\ldots, g_{q-1}]^{T}$
contains the concentration distribution function
$g_{\alpha }$
. In this study, the D3Q7 model was adopted (
$q=7$
) to solve the mass-transfer issue without a significant accuracy loss (Chen et al. Reference Chen, Zhang, Kang and Tao2020). Here,
$\mathbf{N}$
denotes a transformation matrix with
$\mathbf{n}=\mathbf{Ng}$
, and
$\mathbf{n}^{{eq}}$
denotes the equilibrium moment. The diagonal relaxation matrix
$\mathbf{\Lambda }_{{s}}$
incorporates the parameters that are intricately linked to the diffusion coefficient
$D_{{s}}$
. A critical aspect of this model is the additional CST term
$\mathbf{S}_{{CST}}$
, which facilitates the accurate representation of the salt concentration jump at the gas–brine interface in accordance with Henry’s law expressed as
$C^{{b}}=HC^{{g}}$
. Here,
$H$
signifies Henry’s coefficient. For the D3Q7 model, the additional CST term can be written as (Yang et al. Reference Yang, Lei, Wang, Xu, Chen and Luo2023)
\begin{equation} \mathbf{S}_{{CST}}=\frac{1}{4}C\frac{H-1}{Hx_{{b}}+\left(1-x_{{b}}\right)}\left[0,\frac{\partial x_{{b}}}{\partial x},\frac{\partial x_{{b}}}{\partial y},\frac{\partial x_{{b}}}{\partial z},0,0,0\right]^{T}, \end{equation}
where
$x_{{b}}$
is the phase fraction of the brine, which can be calculated using the pseudo-potential
$\varphi$
as
\begin{equation} x_{{b}}=\frac{\varphi ^{{b}}}{\varphi ^{{b}}+\varphi ^{{g}}}. \end{equation}
The salt concentration is calculated as
\begin{equation} C=\sum g_{\alpha }. \end{equation}
Using Chapman–Enskog analysis, the above LB equation can be converted into the advection–diffusion equation as (Yang et al. Reference Yang, Dai, Xu, Liu, Shi and Long2021)
\begin{equation} \frac{\partial C}{\partial t}+\nabla \cdot \left(C\mathbf{u}\right)=\nabla \left(D_{{s}}\,\nabla C-C\frac{D_{{s}}\left(H-1\right)}{Hx_{{b}}+\left(1-x_{{b}}\right)}\,\nabla x_{{b}}\right). \end{equation}
This equation is equivalent to the conventional CST model (Haroun, Legendre & Raynal Reference Haroun, Legendre and Raynal2010) and accounts for the concentration discontinuity at the interface. In this model, the Henry coefficient is assigned an exceedingly high value (
$H=1000$
) to ensure that the salt solute is exclusively present in the brine phase, with
$C^{{g}}\approx 0$
. This strategic approach allows for the simultaneous resolution of (2.4) and (2.6).
2.3.3. Salt precipitate nucleation and growth model
In this work, both heterogeneous (on a solid surface) and homogeneous (in bulk brine) nucleation of salt precipitates were accounted for by calculating the nucleation rate
$J_{{n}}$
based on (2.7). The nucleation process at each lattice site was modelled as a stochastic event, governed by a probability distribution function (Starchenko Reference Starchenko2022)
\begin{equation} {\textrm{d}}P_{{n}}=J_{{n}}\exp \left(-J_{{n}}t\right)\,{\textrm{d}}t. \end{equation}
The nucleation probability
$P_{{n}}$
is calculated by
\begin{equation} P_{{n}}\left(t+\Delta t\right)=P_{{n}}\left(t\right)+\frac{{\textrm{d}}P_{{n}}}{{\textrm{d}}t}\cdot \Delta t. \end{equation}
During each time interval
$N\,\Delta t$
,
$P_{{n}}$
is accumulated using (2.24) for each fluid node. Following every
$N$
time steps,
$P_{{n}}$
is compared with the probability threshold
$P_{{cr}}$
in each node. If
$P_{{n}}\geq P_{{cr}}$
, then nucleation occurs at this node with an initial critical nucleus diameter
$l_{0}$
. If
$P_{{n}}\lt P_{{cr}}$
, then
$P_n$
is reset to zero, and the calculation based on (2.24) is repeated in the next
$N$
time steps. In this study,
$P_{{cr}}$
is assigned a random value within a specific range
$[0.5,1.5]$
at each fluid node.
After nucleation, salt precipitation occurs within the nucleus at a reaction rate described by (2.9). To align with the lattice grid structure, the nuclei are modelled as cubes at each node. The surface area of the nucleus can be calculated as
$A_{{s}}=\chi l^{2}=\chi V_{{s}}^{2/3}$
, where
$V_{{s}}$
is the salt volume at each node. Here,
$\chi =6$
for nuclei in the bulk fluid node (homogeneous nucleation), and
$\chi =5$
for nuclei in the fluid node adjacent to the matrix solid node (heterogeneous nucleation). The salt volume in each fluid node is updated as
\begin{equation} V_{{s}}\left(t+\Delta t\right)=V_{{s}}\left(t\right)+rV_{{m}}\,\Delta t, \end{equation}
where
$V_{{m}}$
denotes the molar volume of the salt. When
$V_{{s}}$
reaches the lattice volume
$\Delta x^{3}$
, the fluid node transitions into a salt node. Adjacent to these new salt nodes, salt growth continues according to the reaction surface area
$A_{{s}}=\Delta x^{2}$
adopted in (2.9).
In bulk fluid nodes, containing a growing nucleus, the salt concentration decreases because the growth reaction is treated as a mass-source term in (2.18). Concurrently, at the surface of the salt nodes or solid nodes with heterogeneous nuclei, the salt mass flux follows the reaction flux as
\begin{equation} -D_{{s}}\,\nabla C=-r\mathbf{n}, \end{equation}
where
$\mathbf{n}$
denotes the surface normal vector of the solid or salt nodes. This boundary condition facilitates the calculation of the unknown concentration distribution function at the solid or salt surface, as demonstrated in our previous work (Yang et al. Reference Yang, Dai, Xu, Liu, Shi and Long2021).
2.3.4. Numerical implementation and validation
The numerical simulation process is illustrated in figure 2. Each time step began with the simulation of gas–brine multiphase flow and phase change, governed by (2.11). The resulting brine phase fraction
$x_{{b}}$
and fluid velocity
$\mathbf{u}$
as inputs for (2.18) to calculate the evolution of mass concentration in the brine, incorporating the salt growth reaction described in (2.26). Based on the calculated salt concentration, the probability of salt nucleation was determined for each
$N$
using (2.24). When the total time step reached an integer multiple of
$N$
, the nucleation probability
$P_{{n}}$
was compared with the threshold
$P_{{cr}}$
to assess whether nucleation occurred. The salt volume
$V_{{s}}$
at each fluid node was updated using (2.25). If
$V_{{s}}$
reached the lattice volume, then the fluid node was converted into a salt node. The resulting solid structure was then integrated into the computational domain, influencing subsequent simulations.
Flowchart of numerical implementation.
At the inlet and outlet of the fracture, the Zou–He boundary condition (Zou & He Reference Zou and He1997) was applied to specify the densities of CO2 and brine vapour. This approach enabled control over the brine vapour mass fraction at the inlet, as well as the injection rate. At the inlet, the density of brine vapour was set to
$\rho _{{in}}^{{b}}=0.117\,\text{kg}\,\textrm{m}^{-3}$
, and the density of CO2 was set to
$\rho _{{in}}^{{g}}=117\,\textrm{kg}\,\textrm{m}^{-3}+\Delta \rho$
, ensuring that the brine vapour mass fraction was
$Y_{{in}}^{{b}}=0.001$
. At the outlet, the brine vapour density was set equal to that in the neighbouring node (outflow boundary condition), while the CO2 density was fixed at
$\rho _{{out}}^{{g}}=117\,\text{kg}\,\textrm{m}^{-3}$
. Different injection rates were achieved by adjusting the pressure difference between the inlet and outlet
$\Delta \rho$
. Since the salt concentration in the gas phase was near zero in the numerical set-up, the inlet salt concentration was set to
$C_{{in}}=0\,\text{mol}\,\textrm{l}^{-1}$
, and the outlet boundary condition was also set to
$\left.\nabla C\cdot \mathbf{n}\right| _{{out}}=0$
.
The numerical procedures outlined in this study were implemented in C++ and enhanced with parallel programming using a message-passing interface to improve computational efficiency. The physical properties used in the simulations are summarised in table 1. Each submodel was rigorously validated for precision by simulating the Stefan problem and brine droplet evaporation (Yang et al. Reference Yang, Lei, Wang, Xu, Chen and Luo2023). Furthermore, the numerical model successfully replicated the microfluidic experiment on salt precipitation, capturing the observed trends in salt precipitation dynamics reported by Ho & Tsai (Reference Ho and Tsai2020) in our previous work (Yang et al. Reference Yang, Lei, Wang, Xu, Chen and Luo2023). A detailed comparison between the numerical and experimental results was conducted, as detailed in the next section. The strong agreement between the numerical predictions and experimental observations reinforces the reliability and accuracy of the developed numerical models.
Physical properties in the simulation.
3. Results and discussions
3.1. Salt precipitation mechanisms
3.1.1. Salt precipitation under experimental conditions – ‘displacement pattern’
To investigate the control mechanisms of salt precipitation during CO2 injection, we first simulated the process at injection rate
$U=8\times10^{-3}\,\text{m}\,\text{s}^{-1}$
, equivalent to the experimental condition (He et al. Reference He, Xu, Ji and Jiang2022) 10 μl min−1 reported by He et al. (Reference He, Xu, Ji and Jiang2022). This injection rate corresponds to a Péclet number
$Pe^{{b}}=0.45$
. The initial salt concentration in the brine was set to 4.0 mol l−1, corresponding to a salt mass fraction of 20 %. The Damköhler number is
$Da_{{s}}=0.85$
.
The numerical results are presented in figure 3(a), with the dynamic evolution process shown in supplementary movie 1. At
$t=50\,\text{s}$
, following the injection of dry CO2, the brine in the matrix adjacent to the inlet begins to evaporate, leading to an elevated salt concentration near the gas–brine interface. By
$t=200\,\text{s}$
, with continued CO2 injection, further brine evaporation occurs, resulting in an additional increase in salt concentration at the gas–brine interface. This elevated concentration induces salt precipitation within the matrix near the fracture. Due to the higher gas flow rate within the fracture, brine in close proximity to the fracture experiences enhanced evaporation compared to the brine deeper within the matrix. Consequently, the regions near the fracture are more prone to salt precipitation. To quantitatively analyse the distribution of salt precipitation, the computational domain was segmented into two regions: a near-fracture region (1 mm wide) and an inside-matrix region, as depicted in figure 3(b). Salt precipitate saturation
$S_{{s}}$
was monitored across these regions. Here,
$S_{{s}}$
is defined as the ratio of the salt precipitation volume in each region to the total pore volume of the entire matrix, providing a measure of salt precipitation. During the initial 200 s of CO2 injection, the near-fracture region exhibited a significantly higher increase in salt precipitation compared to the inside-matrix region. This trend indicates an initial preferential deposition of salt near the fracture, followed by a gradual inward spread of salt precipitation.
(a) Salt concentration evolution during CO2 injection at high injection rate
$U = 8\times10^{-3}\,\text{m}\,\text{s}^{-1}$
and initial concentration 4.0 mol l −1, corresponding to
$Pe^b = 0.45$
,
$Da_s = 0.85$
. White represents salt precipitation. Dark blue corresponding to
$C=0$
mol l−1 represents the gas phase. (b) Salt precipitation saturation variation curves in the near-fracture and inside-matrix regions.
At 200 s, significant salt accumulation was observed in the areas downstream of the fracture, leading to matrix clogging. This clogging hindered downstream brine evaporation and flow, causing evaporation to shift predominantly upstream after 500 s, with the evaporation front gradually advancing inwards. Salt precipitation within the matrix increased more slowly than that for the near-fracture region. Despite partial obstruction near the downstream fracture, most flow channels within the matrix remained open. The differing rates of salt precipitation between the regions were influenced by the competitive interplay of dewetting processes (including evaporation and displacement) and precipitation mechanisms. In the near-fracture zone, higher salt concentrations accelerated precipitation relative to dewetting. Conversely, in the matrix, lower salt concentrations promoted faster dewetting, which impeded stable salt precipitation and nucleation, thereby slowing growth. Additional details on these processes are provided in § 3.1.4.
This configuration enables gradual brine evaporation within the matrix, creating spaces suitable for CO2 sequestration. This mode of salt precipitation is identified as the ‘displacement pattern’, where CO2 effectively displaces brine, allowing it to be sequestered within the reservoir. The numerical results align closely with the experimental findings reported by He et al. (Reference He, Xu, Ji and Jiang2022). Comparison of the experimental and simulation results is shown in figure 4. To match the experimental time scale, a reduced brine vapour diffusion coefficient
$D^{{b}}=0.87\times 10^{-6}\,\textrm{m}^{2}\,\textrm{s}^{-1}$
was employed. Additionally, the injection rate and initial salt concentration were adjusted to ensure that
$Pe^{{b}}$
and
$Da_{{s}}$
remained consistent with the values shown in figure 3 (
$Pe^{{b}}=0.45, \ Da_{{s}}=0.85$
). The numerical results in figures 3 and 4 exhibit similar displacement patterns, further validating the reliability of the parameter adjustments. For comparison,
$t_{0}$
was designated as the moment when the evaporation front first appeared in the experiment. At high
$Pe^{{b}}$
, both the experimental observations and numerical simulations revealed that evaporation initially occurred in the brine near the inlet, followed by a gradual expansion of the evaporation front into the porous matrix. Salt precipitation began near the fracture and subsequently spread into the matrix. The close agreement between the numerical simulations and experimental observations substantiates the credibility of the proposed numerical models for further analysis of salt precipitation mechanisms during CO2 injection.
Comparison of experimental (He et al. Reference He, Xu, Ji and Jiang2022) and numerical results. Experimental conditions are U = 50 μl min−1, C 0 = 4 mol l−1, corresponding to numerical settings Pe b = 0.45, Da s = 0.85. The red line represents the evaporation front in experiment and simulation. Experimental and numerical results agree qualitatively on characteristics of brine evaporation and salt precipitation. Initially, evaporation occurs near the inlet, followed by expansion of the evaporation region into the interior of the matrix, accompanied by formation of salt precipitation. As evaporation continues, salt precipitation propagates downstream along the fracture. Notably, salt precipitation does not completely block flow channels, displaying displacement patterns in both experimental and numerical results.
(a) Salt concentration evolution during CO2 injection at low injection rate U = 5 × 10−4 m s−1 and initial concentration 4.0 mol l−1, corresponding to Pe b = 0.03, Da s = 0.85. (b) Distribution of total salt volume in each cross-section of the near-fracture region at low injection rate U = 5 × 10−4 m s−1 (sealing pattern). (c) Distribution of total salt volume in each cross-section of the near-fracture region at high injection rate U = 8 × 10−3 m s−1 (displacement pattern).
3.1.2. Salt precipitation under lower injection rates – ‘sealing pattern’
In the experimental study, all operational gas injection rates exceeded 10 μm min−1, which is significantly higher than the rates typically employed in near-well injection scenarios for practical engineering applications (He et al. Reference He, Xu, Ji and Jiang2022). However, in real-world CO2 injection processes, certain regions within reservoirs – such as low-permeability zones – often exhibit lower flow rates (Cui et al. Reference Cui, Hu, Ning, Jiang and Wang2023). Presently, experimental studies investigating the mechanisms of salt precipitation under reduced gas injection rates are limited. This knowledge gap underscores the pressing need for further research in this area. In response, we employed numerical simulations to model the CO2 injection process at lower flow rates, aiming to provide valuable insights into the behaviour of salt precipitation under these conditions.
Figure 5(a) presents the numerical results for a lower injection rate U = 5 × 10−4 m s−1, compared with the higher injection rate case discussed in § 3.1.1. At this reduced rate, the Péclet number is reduced to
$Pe^{{b}}=0.03$
. The dynamic process of salt precipitation under these conditions is shown in supplementary movie 2. The mechanism of salt precipitation differed significantly from that observed under high injection rates. Due to the slower pace of gas injection, the introduction of dry CO2 had a limited impact on downstream regions, resulting in minimal brine evaporation near the inlet. At t = 50 s, a noticeable increase in the salt concentration in the brine near the inlet was observed. With continued CO2 injection, salt precipitation began in the matrix near the inlet, leading to clogging by t = 100 s. The accumulation of salt precipitates impeded further brine evaporation upstream. By t = 200 s, brine evaporation near the inlet ceased entirely due to salt clogging. This obstruction shifted the evaporation process downstream to regions without salt precipitation.
To visually analyse the movement of the evaporation leading edge, the total volume of salt precipitates in the near-fracture region was measured across various cross-sections, as shown in figure 5(b). Initially, salt precipitation occurred upstream, and as CO2 injection continued, the salt-blocked area gradually advanced downstream. Consequently, the evaporation front, characterised by high salt concentrations, also moved downstream. This progression is highlighted by the yellow box in figure 5(a), which depicts the advancement of the evaporation front from 500 s to 1000 s. The downstream movement of the evaporation front led to localised salt precipitation near the fracture, inhibiting further brine evaporation and expulsion into the interior matrix regions. Ultimately, by 2000 s, salt precipitates had accumulated in the vicinity of the fracture, resulting in matrix clogging, as shown in figure 5(b). This mode of salt precipitation differs markedly from the displacement pattern, in which upstream flow channels remain open, allowing brine evaporation to expand further, with most of the salt accumulating downstream. This is illustrated in figure 5(c). In contrast, at low flow rates, salt completely blocks the flow channel from upstream to downstream, trapping the majority of the brine within the matrix, and severely hindering CO2 injection. This precipitation behaviour is referred to as the ‘sealing pattern’.
3.1.3. Salt precipitation under higher salt concentrations – ‘breakthrough pattern’
In addition to the gas injection rate, the initial salt concentration in the brine is a critical factor influencing salt precipitation. As the initial salt concentration decreases, the dewetting process gradually becomes comparable to, and eventually dominates, salt precipitation. As illustrated in figure 6(a), lower initial salt concentrations (
$Da_{{s}}=0.05$
) in the brine significantly reduced the likelihood of salt precipitation, even under low-flow-rate conditions. This mitigation prevented matrix clogging. By t = 1000 s, no substantial salt precipitation was observed (the detailed evolution of salt concentration at different moments is provided in figure S1 of the supplementary material). Under such conditions, the overall precipitation process transitions from the sealing pattern to the displacement pattern, allowing better brine displacement and maintaining open flow channels for CO2 injection.
(a) Salt concentration at t = 1000 s during CO2 injection at low injection rate U = 5 × 10−4 m s−1 and lower initial concentration 1.0 mol l−1, corresponding to Pe b = 0.03, Da s = 0.05. Salt precipitation shifts from sealing to displacement patterns. (b) Salt concentration at t = 1000 s during CO2 injection at injection rate 1.7 × 10−3 m s−1 and higher initial concentration 5.75 mol l−1, corresponding to Pe b = 0.10, Da s = 1.76. Salt precipitation shows a breakthrough pattern.
Conversely, higher initial salt concentrations in the brine increased the likelihood of salt precipitation during evaporation. Figure 6(b) illustrates this phenomenon at initial salt concentration 5.75 mol l−1 under high-injection-rate conditions with
$Pe^{{b}}=0.10,\ Da_{{s}}=1.76$
. The evaporation of brine resulted in intensified salt precipitation near the fracture, ultimately leading to matrix clogging. Before significant brine displacement could occur, the matrix became fully obstructed. This obstruction allowed only a small fraction of the gas to penetrate the salt precipitation blockage near the inlet and enter the matrix (further details are provided in supplementary figure S2 and movie 3). As a result, the evaporation process was largely restricted to areas adjacent to the protruding gas. The additional salt deposits further impeded flow pathways. Eventually, brine displacement by the gas phase was confined to a narrow region near the advancing gas front, as shown in figure 6(b).
In contrast to the displacement pattern observed in figure 3 with an initially low salt concentration, the brine evaporation region in figure 6(b) was confined to a localised flow channel, leaving the matrix only partially permeable. While this condition allowed for continued CO2 injection, the dewetting process was less efficient compared to the displacement pattern. This specific dynamic of salt precipitation, characterised by localised flow and partial permeability, is referred to as the ‘breakthrough pattern’.
3.1.4. Mechanism analysis and comparison of different salt precipitation patterns
The evolution of salt precipitation during CO2 injection, along with the brine vapour mass fraction, is illustrated in figure 7(a). At a low injection rate
$U=5\times 10^{-4}\,\textrm{m}\,\textrm{s}^{-1}$
, the corresponding Péclet number is reduced to
$Pe^{{b}}=0.03$
, indicating a diffusion-controlled regime for the distribution of brine vapour. As illustrated in figure 7(a), a significant mass-fraction gradient of brine vapour exists within the fracture. Near the inlet, the brine vapour mass fraction is comparatively low, promoting rapid evaporation. Due to this rapid evaporation upstream, the salt concentration near the inlet is higher, increasing the likelihood of salt precipitation. In this case,
$Da_{{s}}=0.85$
is relatively high, reflecting that the salt precipitation rate dominates the mass transfer rate, making precipitation more likely. Because the evaporation rate is low due to the minimal CO2 injection rate, salt deposits stabilise near the fracture, gradually obstructing the upstream matrix. This obstruction causes dry CO2 to diffuse further downstream, shifting the evaporation and salt precipitation zones downstream within the matrix. By 2000 s, the matrix becomes progressively fully obstructed by salt precipitation, culminating in a sealing precipitation pattern. To further analyse the competition between dewetting and salt precipitation, we calculated their respective rates,
$J_{{d}}$
and
$J_{{s}}$
, under three distinct salt precipitation patterns at 500 s, as depicted in figure 7(e). Concurrently, we recorded the amount of salt precipitation in the near-fracture and inside-matrix regions at 1000 s to represent the salt distribution characteristics, as depicted in figure 7( f). As indicated in figure 7(e), under the sealing pattern, the dewetting process occurs at a significantly slower rate, being dominated by salt precipitation. Consequently, as shown in figure 7
(f), salt primarily accumulates in the near-fracture region, impeding CO2 injection and further reducing the dewetting rate. The inability of the gas to penetrate the matrix interior results in zero salt precipitation in the inside-matrix region, highlighting the localised nature of the sealing pattern.
Evolution of brine vapour mass fraction for (a) low injection rate U = 5 × 10−4 m s−1 (Pe b = 0.03), and (b) high injection rate U = 8 × 10−3 m s−1(Pe b = 0.45). (c) Streamlines at t = 500 s for low injection rate U = 5 × 10−4 m s−1. (d) Streamlines at t = 200 s for high injection rate U = 8 × 10−3m s−1. (e) Dewetting and salt precipitation rates at t = 500 s for three typical salt precipitation patterns. ( f) Salt saturation in near-fracture and inside-matrix regions at t = 1000 s for three typical salt precipitation patterns.
In contrast, at higher flow rates, as shown in figure 7(b), with Péclet number
$Pe^{{b}}=0.45$
, the mass transfer of brine vapour transitions to advection control. Rapid CO2 injection maintains consistently low brine vapour levels throughout the fracture, resulting in relatively high evaporation rates in both upstream and downstream areas. While mass transfer in the fracture is advection-controlled, the slower flow rate within the matrix leads to a diffusion-controlled mechanism for brine vapour distribution. This creates a pronounced vapour mass-fraction gradient in the matrix, as depicted in figure 7(b). Consequently, the evaporation rate of brine within the matrix is slower than that in the fracture. As illustrated by the streamline diagram in figure 7(d), upstream brine is displaced downstream by the injected gas alongside brine evaporation, as the capillary number reaches
$Ca=2.86$
. In contrast, at low flow rates, the displacement effect is negligible, as shown in figure 7(c). At high injection rates, a significant reduction in brine volume is attributed to its displacement by the gas phase. Unlike evaporation-induced salt precipitation, displacement does not lead to an increased salt concentration. Thus, as illustrated in figure 3, the brine concentration in the matrix does not increase substantially. The brine is either evaporated or displaced by the gas phase before salt precipitation can occur and stabilise through nucleation and growth. As shown in figure 7(e), under the displacement pattern, the dewetting rate significantly exceeds the salt precipitation rate. Salt precipitation progresses slowly and is more evenly distributed between the near-fracture and inside-matrix regions, as depicted in figure 7( f). The absence of salt clogging facilitates faster dewetting, enabling CO2 to penetrate the matrix more efficiently for sequestration. This dynamic behaviour leads to the formation of the displacement pattern, characterised by open flow channels and effective brine displacement.
In the breakthrough pattern, both the rates of salt precipitation and dewetting are pronounced because
$Pe^{{b}}$
and
$Da_{{s}}$
maintain high values. Although rapid salt precipitation obstructs some flow channels, localised channels remain open, permitting gas entry. As a result, the dewetting rate is not as significantly reduced as in the sealing pattern, as shown in figure 7(e). From figure 7( f), it is evident that salt precipitation initially accumulates extensively in the near-fracture region and then gradually expands into localised areas within the matrix. In the retained flow channels, brine continues to evaporate or be displaced by the gas phase, facilitating the emergence of the breakthrough pattern.
These quantitative findings further emphasise that the interplay between dewetting and salt precipitation rates is critical in determining the resulting salt precipitation patterns. The results highlight fundamental differences in the underlying mechanisms of salt precipitation between experimental conditions and practical engineering injection rates. In real-world reservoirs, the impact of salt precipitation may be more pronounced than that observed experimentally, posing a greater risk to CO2 injection efficiency.
3.2. Regime map according to operation conditions
The analysis above indicates that both the CO2 injection rate and the initial salt concentration are critical factors influencing the resulting salt precipitation pattern. To gain a comprehensive understanding of the mechanisms driving salt precipitation, and to provide practical guidance for CO2 injection engineering applications, simulations were conducted under varying conditions of CO2 injection rate and initial salt concentration. These simulations were performed in a matrix with pore-throat width 85 μm. The synthesised results from these simulations are presented systematically in the regime diagram shown in figure 8, which maps the identified regimes as follows.
Regime diagram of salt precipitation in matrix with porosity 0.48.
-
(i) Low initial concentration, displacement pattern. When the initial salt concentration is low, registering at less than half the saturation concentration, the Damköhler number is also low,
$Da_{{s}}\sim O(10^{-1})$
. Under these conditions, the mass-transfer rate dominates the salt precipitation rate. The increase in salt concentration due to brine evaporation is insufficient to induce substantial salt precipitation, making the formation of salt precipitates less likely, and significantly mitigating the clogging effect within the matrix. In this scenario, the salt precipitation rate remains consistently lower than the dewetting rate. Regardless of the gas injection rate, the gas is able to permeate the matrix, resulting in the formation of a displacement salt precipitation pattern. -
(ii) High initial concentration and low gas injection rate, sealing pattern. With a high initial concentration
$Da_{{s}}\sim O(10^{0})$
, the salt precipitation pattern varied depending on the gas injection rate. At a low gas injection rate
$Pe^{{b}}\sim O(10^{-2})$
, the salt precipitation front advanced gradually from the inlet towards the outlet. Due to the weak advective term, the brine evaporation rate was very low, allowing salt precipitation to dominate dewetting rates. This dominance facilitated extensive and stable salt deposition near the fracture, progressively obstructing the matrix channels, and ultimately resulting in the formation of a sealing pattern. -
(iii) High initial concentration and medium gas injection rate, breakthrough pattern. As the gas injection rate increases to
$O(10^{-2})\lt Pe^{{b}}\lt O(10^{-1})$
, a portion of the gas phase penetrates the brine and enters the matrix, forming localised gas-phase flow channels. Under these conditions, the rates of salt precipitation and dewetting become comparable. Brine evaporation and subsequent salt precipitation occur near these breakthrough channels, leading to the formation of a breakthrough pattern. -
(iv) High initial concentration and high gas injection rate, displacement pattern. With a further increase in the gas injection rate to
$Pe^{{b}}\sim O(10^{-1})$
, the injected gas phase effectively displaces the brine downstream. Unlike evaporation, this displacement process does not significantly elevate the salt concentration, reducing the likelihood of substantial salt precipitation and severe clogging of the matrix. As a result, a displacement precipitation pattern emerges.
From the regime diagram in figure 8, high injection rates typically result in salt precipitation following a displacement mechanism, effectively evaporating or displacing most of the brine. This mechanism enhances injection efficiency and aligns well with experimental observations. Conversely, at lower flow rates, which are more common in practical engineering scenarios, salt precipitation predominantly occurs in a sealing pattern, leading to matrix clogging and a significant reduction in injection efficiency.
3.3. Effect of matrix porosity
Based on the above mechanistic understanding, the gas injection rate and initial salt concentration are critical factors influencing the salt precipitation pattern. To further enhance CO2 injection efficiency, it is essential to identify additional controllable factors. Therefore, we simulated the salt precipitation process under various reservoir conditions. The first factor examined was the effect of reservoir porosity on salt precipitation patterns. Figure 9(a) illustrates the salt precipitation process in a matrix with porosity 0.68. The gas injection rate and initial salt concentration were set at 1.7 × 10−3 m s−1 and 5.75 mol l−1, respectively, corresponding to
$Pe^{{b}}=0.10, \ Da_{{s}}=1.76$
.
(a) Salt concentration evolution during CO2 injection at injection rate U = 1.7 × 10−3 m s−1 and initial concentration 5.75 mol l−1, corresponding to Pe b = 0.10, Da s = 1.76. Porosity of the matrix is 0.68. (b) Localised zoomed-in view of numerical result at t = 800 s to present hydrophilic salt precipitates drawing in neighbouring brines.
Initially, at t = 100 s, brine evaporation led to clogging near the inlet due to salt precipitation. This behaviour closely mirrors the scenario observed with lower matrix porosity 0.48 in supplementary figure S2. By t = 200 s, a portion of the gas began to permeate the matrix. However, owing to the higher porosity and broader pore-throat flow channels, the resulting salt precipitate did not completely obstruct these channels. Consequently, the injected gas continued to displace the brine along the breakthrough flow channels, allowing further brine evaporation. By t = 500 s, the gas phase penetrated deeper into the matrix, leaving many flow channels unobstructed by salt precipitation. Due to the hydrophilic properties of the substrate, a thin brine film persisted on the surfaces of the matrix grains during gas injection. The evaporation of these films facilitated the development of a salt precipitation layer on the matrix surface. These hydrophilic salt deposits attracted neighbouring brine, as shown in figure 9(b), contributing to the gradual thickening of the salt layer, particularly after t = 800 s. The progressive evolution of salt growth and brine film dynamics is captured visually in supplementary movie 4.
This process is consistent with the experimentally observed phenomenon of brine backflow under hydrophilic conditions (He et al. Reference He, Jiang and Xu2019). In high-porosity environments, the wider flow channels provide sufficient space for the formation of salt layers. By contrast, in low-porosity matrices, salt precipitates directly obstruct the narrower flow channels, leading to significant flow impedance. The expansive channel structure in high-porosity matrices facilitates further progression of the gas phase, allowing for continuous brine displacement and evaporation. This results in a transition of the salt precipitation pattern from a breakthrough pattern to a displacement pattern, in stark contrast to the behaviour observed in low-porosity settings, where flow is more readily obstructed by salt precipitation.
To further explore the effect of porosity on salt precipitation patterns, simulations were conducted on matrices with varying porosity values, with the results presented in figure 10. Under the same gas injection rate and initial salt concentration, the salt precipitation pattern transitioned sequentially from sealing to breakthrough, and ultimately to displacement as porosity increased. The thresholds marking these transitions were identified at porosity values 0.45 and 0.55, respectively.
(a) Salt precipitation patterns under different matrix porosity at the same injection rate U = 1.7 × 10−3 m s−1and initial concentration 5.75 mol l−1, corresponding to Pe b = 0.10, Da s = 1.76. (b) Three-dimensional regime diagram of salt precipitation.
As shown in figure 7(e), the dewetting and salt precipitation rates are comparable for the current injection rate and initial salt concentration, which correspond to the breakthrough pattern. The resulting salt precipitation pattern depends on the relative timing of pore-throat blockage by salt, and the advancement of the evaporation front into the matrix. When the porosity is below 0.45, the narrowest pore throats are less than 70 μm wide. In this case, salt precipitation blocks the pore throats before the evaporation front can penetrate the matrix, resulting in complete blockage and the formation of a sealing pattern. Conversely, when the porosity exceeds 0.55, the pore-throat width is greater than 110 μm. Under these conditions, the growth of salt precipitation is insufficient to completely block the flow channels, as illustrated in figure 10(a). This allows the brine evaporation front to advance deeper into the matrix, resulting in a displacement pattern. For porosities ranging between 0.45 and 0.55, with pore-throat widths between 70 and 110 μm, most pore throats near the fracture become blocked by salt precipitation. However, the evaporation front is still able to penetrate at certain locations, leading to the emergence of a breakthrough pattern.
From this analysis, the salt precipitation pattern is determined by the time ratio
$\tau$
, defined as the ratio of the time required for salt precipitation to block the pore throats,
$t_{{s}}$
, to the time required for the brine evaporation front to expand into the matrix,
$t_{{d}}$
, expressed as
\begin{equation} \tau =\frac{t_{{s}}}{t_{{d}}}=\frac{l/J_{{s}}}{L/J_{{d}}}, \end{equation}
where
$l$
denotes the pore-throat width,
$L$
is the fracture width, and
$J_{{s}}$
and
$J_{{d}}$
are the salt precipitation and dewetting rates, respectively, as defined in figure 7(e). According to numerical results, the boundaries between different salt precipitation patterns can be characterised by the time ratio
$\tau$
. When
$\tau$
was less than 0.8, salt precipitation rapidly clogged the flow channels of the matrix, resulting in a sealing pattern. When
$\tau$
exceeded 1.2, the evaporation front advanced into the interior of the matrix, preventing salt precipitation from blocking the flow channels, leading to a displacement pattern. For values between 0.8 and 1.2, salt precipitation followed a breakthrough pattern, where the evaporation front partially expanded into the matrix. This conclusion is further validated by calculations using the data from figure 7(e). For sealing and displacement conditions,
$\tau$
values are 0.2 and 2.7, respectively, both consistent with the expected range. The analysis indicates that higher porosity and wider pore throats favour improved CO2 injection efficiency by enabling displacement patterns. The analysis of
$\tau$
reveals that in addition to
$Pe^{{b}}$
and
$Da_{{s}}$
, changes in pore-throat width
$l$
influence shifts in the salt precipitation pattern. Consequently, the boundaries in the regime diagram shown in figure 8 shift for reservoirs with different porosities. As porosity increases, the displacement pattern region expands, while decreased porosity enlarges the sealing pattern region. More comprehensively,
$Pe^{{b}}$
,
$Da_{{s}}$
and the geometrical parameter
$l$
account for diverse CO2 injection mechanisms affected by both reservoir microstructure and operational conditions. This enables the construction of a three-dimensional regime diagram, as shown in figure 10(b), offering a more accurate determination of salt precipitation patterns compared to the two-dimensional diagram in figure 8.
3.4. Effects of matrix wettability
Prior experimental studies have highlighted the substantial influence of matrix wettability on the salt precipitation process (He et al. Reference He, Jiang and Xu2019). In hydrophilic pore structures, salt precipitation tends to clog pore throats, leading to reduced permeability. Conversely, in hydrophobic matrices, salt primarily forms isolated large crystal structures. To investigate these effects, our preliminary study assessed the role of wettability on salt precipitation dynamics. For the simulations, we employed a medium gas injection rate 1.7 × 10−3 m s−1 and a high initial salt concentration 5.75 mol l−1 in a matrix with porosity 0.48. As shown in figure 8, under hydrophilic conditions, the salt precipitation process followed the breakthrough pattern. To further explore the impact of wettability, we simulated scenarios where the matrix exhibits hydrophobic characteristics, with contact angle
$\theta_m=120^\circ$
.
(a) Salt concentration evolution during CO2 injection at injection rate U = 1.7 × 10−3 m s−1 and initial concentration 5.75 mol l−1 within the hydrophobic matrix, corresponding to Pe b = 0.10, Da s = 1.76. (b) Localised zoomed-in view of salt nucleation, growth and aggregation. (c) Localised zoomed-in view of salt concentration distribution at t = 200 s for both hydrophobic and hydrophilic scenarios at injection rate U = 1.7 × 10−3 m s−1 and initial concentration 5.75 mol l−1. (d) Normalised permeability curves for hydrophilic and hydrophobic conditions at U = 1.7 × 10−3 m s−1, C 0 = 5.75 mol l−1. (e) Injected CO2 saturation curves in the matrix for hydrophilic and hydrophobic conditions.
The numerical results are presented in figure 11(a) and supplementary movie 5. Under hydrophobic conditions, the matrix’s capillary forces more effectively promote the displacement of brine by the gas phase. Consequently, a larger volume of brine is displaced under hydrophobic conditions compared to hydrophilic conditions (see supplementary figure S2). Salt precipitation begins during brine evaporation, as shown at t = 100 s. Hydrophobic conditions facilitate a distinct salt precipitation pattern compared to hydrophilic conditions. Predominantly, the salt forms as large crystals aggregated within macropores, consistent with experimentally observed phenomena (He et al. Reference He, Jiang and Xu2019). This configuration results in a less pronounced clogging effect on the matrix, enabling continued brine displacement and evaporation. Due to the inherently hydrophilic nature of salt, brine is drawn towards pre-existing salt particles. As a result, salt precipitation tends to propagate along these salt clusters, as illustrated at t = 500 s. Figure 11(b) provides a localised magnification of the salt precipitation growth process, while supplementary movie 5(b) captures the dynamic evolution. Initially, salt precipitates nucleate within macropores and begin to grow between 50 s and 100 s. By 200 s, the affinity of salt for brine becomes evident, leading to further growth along brine pathways, and the gradual occupation of large pores. The configuration of the salt precipitate aggregates is shaped primarily by the brine pathways, allowing substantial salt precipitation without causing excessive clogging of matrix flow paths. This facilitates the continuous evaporation and displacement of brine. Compared to hydrophilic conditions, where the precipitation pattern followed the breakthrough mechanism, the hydrophobic matrix transitions to a displacement pattern. This shift significantly enhances injection efficiency, as shown at t = 1000 s. These findings align with experimental observations (He et al. Reference He, Jiang and Xu2019), where salt forms larger crystals under hydrophobic conditions, and brine is more effectively displaced by CO2.
To further elucidate the mechanisms underlying the contrasting salt precipitation behaviours in hydrophilic and hydrophobic matrices, a detailed analysis of the concentration distribution at the gas–water interface was conducted. The results, presented in figure 11(c), highlight key differences in interface configurations and their impact on salt precipitation. On the hydrophobic substrate, the brine interface adopts a convex shape towards the pores. This configuration intensifies evaporation at the liquid surface, resulting in a locally higher salt concentration near the interface. This increased salt concentration promotes nucleation and growth of salt precipitates near the liquid surface, causing salt precipitation to be predominantly concentrated in larger pores. In contrast, in the hydrophilic matrix (figure 11 c), the brine level is located primarily at the pore throats, forming a concave shape towards the interior of the brine. This configuration accelerates evaporation near the interface, particularly at the solid–liquid–gas three-phase junction, where the brine film is thinner, and the salt concentration rises more rapidly. As a result, salt precipitation is more likely to occur near the pore throats, leading to clogging of the matrix flow channels. A comparison of supplementary movies S5(b) and S5(c) effectively illustrates these differences. In hydrophobic matrices, large salt precipitates are more likely to form within macropores, which do not obstruct the primary flow paths of the matrix.
As illustrated in figure 11(d), under hydrophobic conditions, salt precipitation results in a relatively minor reduction in permeability compared to hydrophilic scenarios. At salt saturation level S s = 0.03, the reservoir permeability under hydrophobic conditions is nearly 4–5 times higher than that observed under hydrophilic conditions. Furthermore, figure 11(e) demonstrates that hydrophobic matrices are significantly more conducive to CO2 injection. At 1000 s, the saturation level of injected CO2 in the hydrophobic matrix is approximately four times higher than in the hydrophilic matrix. This indicates that salt precipitation in hydrophobic environments has a substantially reduced clogging effect on the reservoir, thereby facilitating continued brine evaporation and enhancing CO2 injection efficiency.
(a) Curves of injected CO2 saturation into matrix during CO2 injection. Red, blue and green curves represent displacement, breakthrough and sealing modes, respectively. (b) Relative evaporation rate curves for different salt precipitation patterns. (c) Normalised permeability curves for different salt precipitation patterns.
3.5. Implications for CO2 injection
To further quantify injection efficiency across various salt precipitation patterns, the curves of injected gas saturation in a matrix with a porosity 0.48 (
$l=85\,\unicode{x03BC} \text{m}$
) under differing conditions are shown in figure 12(a). The red, blue and green curves represent displacement (D), breakthrough (B) and sealing (S) patterns, respectively. The graph clearly indicates that the rate of water saturation reduction in the matrix is generally the fastest in the displacement mode. Remarkably, even at a slow gas injection rate (U = 5 × 10−4 m s−1), the rate of increase in gas saturation for the displacement pattern is comparable to that observed at a high gas injection rate (U = 8 × 10−3 m s−1) in the breakthrough mode, despite the tenfold difference in injection rates. The displacement precipitation pattern is advantageous as it allows enhanced penetration of the gas phase into the matrix, thereby improving injection efficiency. In contrast, the sealing pattern, characterised by significant blockage of the matrix, impedes gas-phase entry, leading to much slower gas saturation, and making it less favourable for CO2 sequestration. Thus in CO2 injection processes, achieving a displacement precipitation pattern is preferable. Under the current porosity conditions, this can be achieved by identifying or modifying saline aquifers so that the initial salt concentration is less than half of the saturated concentration
$Da_{{s}}\lt O(10^{-1})$
. Alternatively, increasing the gas injection rate to ensure that
$Pe^{{b}}\gt O(10^{-1})$
can also facilitate the displacement pattern.
In assessing the macroscopic impact of salt precipitation patterns on CO2 injection efficiency, we conducted a statistical analysis under three specific conditions: displacement (8 × 10−3 m s−1, 4.0 mol l−1), breakthrough (1.7 × 10−3 m s−1, 4.0 mol l−1) and sealing (5 × 10−4 m s−1, 4.0 mol l−1) patterns. Two key statistical parameters were evaluated: relative drying rate and matrix permeability. The relative drying rate
$\beta$
is defined as the ratio of the brine saturation decrease rate with and without considering the effects of salt precipitation. It can be calculated as
\begin{equation} \beta \left(t\right)=\frac{S_{{b}}\left(t\right)-S_{{b}}\left(t-\Delta t\right)}{\left(S_{{b}}\left(t\right)-S_{{b}}\left(t-\Delta t\right)\right)_{\text{at }{C_{0}}=0\text{ mol}\,\textrm{l}^{-1}}}, \end{equation}
where S b is the brine saturation. The matrix permeability K was calculated by modelling the single-phase flow from the fracture to the matrix. This parameter is expressed as the normalised value K/K 0, where K 0 is the matrix permeability in the absence of salt precipitation.
Figures 12(b) and 12(c) illustrate the calculation results. As shown in figure 12(b), the relative drying rate decreases significantly with increasing salt precipitation saturation, primarily due to the blocking effect of salt precipitation on the matrix. In the sealing pattern, a sharp decline in
$\beta$
is observed, approaching near zero at salt precipitation saturation 0.004. This drastic reduction significantly impacts sequestration efficiency by severely hindering brine removal and gas injection.
The permeability K trends demonstrated a notable initial decrease across all patterns, ultimately stabilising at distinct values: 0.32 for the displacement pattern, 0.25 for the breakthrough pattern, and 0 for the sealing pattern. These results highlight the particularly adverse impact of salt precipitation in the sealing pattern on reservoir permeability. These findings underscore the critical importance of considering salt precipitation patterns in studies of CO2 injection efficiency. A transition from a displacement pattern to a sealing pattern can result in a tenfold decrease in CO2 sequestration capacity, as indicated by the relative drying rate
$\beta$
. Concurrently, the permeability of the reservoir matrix is severely compromised, leading to a substantial reduction in fluid flow capacity. The sealing pattern is therefore highly detrimental to CO2 injection efficiency, and is not conducive to effective sequestration operations.
Beyond optimising gas injection strategies, identifying appropriate reservoir characteristics is equally crucial for enhancing CO2 injection efficiency. As discussed in § 3.4, salt precipitation in hydrophobic reservoirs tends to concentrate in larger pores, causing less blockage and thereby facilitating CO2 injection. Thus identifying hydrophobic reservoirs can significantly improve injection efficiency.
Reservoir porosity is also a critical factor influencing CO2 injection efficiency. In actual reservoirs, the average pore-throat sizes
$l$
typically range from a few microns to several hundred microns (Liu et al. Reference Liu, Yang, Bai, Wu, Zhang, Liu and Xiao2021), with the corresponding
$\tau$
values spanning different salt precipitation regimes. To optimise injection, adjustments must be made according to the corresponding
$Pe^{{b}}$
and
$Da_{{s}}$
boundaries in figure 10(b) under varying reservoir geometrical conditions.
4. Conclusions
In this study, we conducted a pore-scale numerical examination of the mechanisms of salt precipitation during CO2 injection into fractured matrices, using microfluidic experiments as a basis. The models incorporated complex multiphysics processes, including gas–brine multiphase flow with phase change, brine vapour advection and diffusion, salt concentration evolution, and the nucleation and growth of salt precipitates. The numerical findings are qualitatively aligned with experimental results, demonstrating the reliability of the models in predicting salt precipitation characteristics.
We explored salt precipitation under varying CO2 injection rates and initial brine salt concentrations to assess their effects on injection efficiency. Three distinct salt precipitation patterns – displacement, breakthrough, and sealing – were identified from the numerical results. These patterns arose from the interplay between brine dewetting and salt precipitation rates.
-
(i) Displacement pattern. At high CO2 injection rates, which correspond to typical experimental conditions, the numerical results strongly aligned with experimental observations. High injection rates promoted brine evaporation near the inlet, with initial salt deposition near the fracture gradually extending into the matrix. This pattern favoured displacement over precipitation due to the reduced clogging effect of salt in the matrix, enabling a higher dewetting rate than salt precipitation. Consequently, this pattern was the most favourable for gas injection.
-
(ii) Sealing pattern. At slower injection rates, typical of engineering applications, salt precipitation accumulated near the inlet, obstructing further evaporation and brine displacement. This led to matrix clogging near fractures, severely impeding gas injection. In this scenario, salt precipitation dominated dewetting, forming stable deposits that blocked the fracture–matrix interface. The numerical findings indicated that predictions based on experimental CO2 injection rates might be overly optimistic, highlighting the need for mitigation strategies to address salt precipitation.
-
(iii) Breakthrough pattern. At brine salt concentrations near saturation, salt precipitation occurred rapidly near the fracture, even at high gas injection rates. The dewetting and salt precipitation rates were comparable, allowing only limited gas penetration into the matrix. Evaporation and brine displacement occurred primarily in localised regions near the breakthrough gas phases.
A regime diagram delineating these salt precipitation patterns was constructed based on injection rate and initial salt concentration, assuming typical matrix porosity 0.48. At low initial salt concentrations with
$Da_{{s}}\sim O(10^{-1})$
, precipitation consistently followed a displacement pattern, irrespective of the gas injection rate. Conversely, at higher salt concentrations with
$Da_{{s}}\sim O(10^{0})$
, the precipitation pattern transitioned from sealing to breakthrough, and finally to displacement as injection rates increased from
$Pe^{{b}}\sim O(10^{-2})$
to
$Pe^{{b}}\sim O(10^{-1})$
. Analysis of permeability and relative drying rate revealed a significant reduction in both CO2 sequestration capacity and reservoir permeability under sealing conditions, underscoring the importance of avoiding such scenarios in engineering applications.
Additionally, matrix properties significantly influenced precipitation behaviour. In highly porous matrices, wider flow channels reduced the risk of clogging by salt deposits, enabling greater gas penetration, and substantially improving injection efficiency. The time ratio
$\tau$
proved useful for evaluating the impact of porosity on salt precipitation patterns. By incorporating pore-throat width
$l$
, the phase diagram could be extended to three dimensions, accounting for varying reservoir structures, and providing more accurate boundaries between precipitation patterns. Hydrophobic substrates further enhanced injection efficiency by promoting gas displacement over brine, and facilitating the formation of larger salt particles, which reduced clogging risks. Thus hydrophobic and highly porous reservoirs are optimal for CO2 injection.
The numerical simulations in this study were based on an idealised microfluidic pore structure, which introduces certain limitations. In actual saline aquifers, the pore structure is more complex, and further investigation is required using realistic pore structures determined by micro-CT imaging. Additionally, upscaling models need to be developed to extend the pore-scale insights to larger scale. These will be the focus of our future research.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.266.
Acknowledgements
The research is supported by the UK Engineering and Physical Sciences Research Council (EPSRC) under grant nos EP/W026260/1 and EP/V034154/1. ARCHER2 supercomputing resources provided by the EPSRC project ‘UK Consortium on Mesoscale Engineering Sciences (UKCOMES)’ (grant no. EP/X035875/1) are gratefully acknowledged. This work made use of computational support by CoSeC, the Computational Science Centre for Research Communities, through UKCOMES. This work is also supported by the National Natural Science Foundation of China (no. 52206014).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding authors upon reasonable request.
Author contributions
J.Y. and Q.X. conceptualised the methodology and conducted the data analysis. J.Y., T.L., G.W. and J.C. carried out the simulations, while K.H.L. and Q.X. oversaw the programme administration. All authors contributed equally to the writing of the paper.
Open access
