Next Article in Journal
The Use of a Regulating Transformer for Shaping Power Flow in the Power System
Next Article in Special Issue
An Integrated Approach to Reservoir Characterization for Evaluating Shale Productivity of Duvernary Shale: Insights from Multiple Linear Regression
Previous Article in Journal
Improve Oil Recovery Mechanism of Multi-Layer Cyclic Alternate Injection and Production for Mature Oilfield at Extra-High Water Cut Stage Using Visual Physical Simulation Experiment
Previous Article in Special Issue
Automatic Evaluation of an Interwell-Connected Pattern for Fractured-Vuggy Reservoirs Based on Static and Dynamic Analysis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Micro-Scale Lattice Boltzmann Simulation of Two-Phase CO2–Brine Flow in a Tighter REV Extracted from a Permeable Sandstone Core: Implications for CO2 Storage Efficiency

1
School of Earth and Space Sciences, Peking University, Beijing 100871, China
2
PetroChina Company Ltd., Beijing 100011, China
3
Research Institute of Petroleum Exploration and Development, Beijing 100083, China
*
Author to whom correspondence should be addressed.
Energies 2023, 16(3), 1547; https://doi.org/10.3390/en16031547
Submission received: 13 December 2022 / Revised: 14 January 2023 / Accepted: 23 January 2023 / Published: 3 February 2023

Abstract

:
Deep saline permeable sandstones have the potential to serve as sites for CO2 storage. However, unstable CO2 storage in pores can be costly and harmful to the environment. In this study, we used lattice Boltzmann (LB) simulations to investigate the factors that affect steady-state CO2–brine imbibition flow in sandstone pores, with a focus on improving CO2 storage efficiency in deep saline permeable sandstone aquifers. We extracted three representative element volumes (REVs) from a digital rock image of a sandstone core and selected a tighter REV in the upper subdomain so that its permeability would apparently be lower than that of the other two based on single-phase LB simulation for further analysis. The results of our steady-state LB simulations of CO2–brine imbibition processes in the tighter REV under four differential pressures showed that a threshold pressure gradient of around 0.5 MPa/m exists at a differential pressure of 200 Pa, and that higher differential pressures result in a greater and more linear pressure drop and stronger channelization after the flow are initiated. Furthermore, we conducted simulations over a range of target brine saturations in the tighter REV at the optimal differential pressure of 400 Pa. Our findings showed that the relative permeability of CO2 is greatly reduced as the capillary number falls below a certain threshold, while the viscosity ratio has a smaller but still significant effect on relative permeability and storage efficiency through the lubrication effect. Wettability has a limited effect on the storage efficiency, but it does impact the relative permeability within the initial saturation range when the capillary number is low and the curves have not yet converged. Overall, these results provide micro-scale insights into the factors that affect CO2 storage efficiency in sandstones.

1. Introduction

Climate change and induced greenhouse effects mainly caused by excessive CO2 emissions by human activities pose great challenges to the earth system [1,2,3,4,5]. CO2 emitted in the atmosphere is mainly sourced from the combustion of fossil fuels and electricity generated from the power plants. Amid a great variety of solutions to this problem, CO2 capture and sequestration (CCS) are assumed to be efficient techniques [6,7]. CCS entails the separation of CO2 from other substances produced by power plants or other similar human industrial activities, the compression and transportation of them to the storage sites and the injection of supercritical CO2 into onshore or offshore geological formations that provide ample storage space [8,9]. In this sense, the geological capture and sequestration of CO2 is crucial for mitigating the negative impact of climate change and derived environmental issues [6,10,11]. In addition, most of geological formations that are idealized candidates for CO2 storage are deep saline sandstone reservoirs, which are permeable and widely distributed worldwide [12,13,14]. According to [15], the confidence level of CO2 storage efficiency in sandstone aquifers (90% or higher) is much higher than that for limestone ones (5%). Meanwhile, the preferable mineralogy and geochemistry of sandstone saline aquifers provide exceptional strengths for the selection [16,17,18]. Several experimental investigations [18,19,20,21,22] at the macroscopic scale in the literature have been conducted on sandstones, with results showing the big variance in storage efficiency, depending on sandstone type but also complicated pore structure inside. As a result, the stable storage efficiency of CO2 from sandstones is still a research hotspot and is yet to be fully evaluated at a microscopic scale. Since one can rarely bear the great cost of the unstable CO2 storage after the CCS operation, which could give rise to large economic losses and secondary environmental consequences [23,24,25,26], a better and integrated understanding of the physics of two-phase CO2–brine flow in the porosity of sandstones is extraordinarily essential for preventing CO2 leakage from the formation pores by serving as permanent storage. As such, one must develop reliable models in addition to experimental measurements for both the sandstone formation and the physics of two-phase CO2–brine flow in it.
Multiscale structure of the sandstone formations and the fact that various length scales exist, ranging from pore [9,27,28,29] to core [30,31,32], macroscopic up until reservoir scale [9,29,33]. They necessitate the development of different models at each scale and the integration of all of them into a uniform framework that characterizes the heterogeneities and anisotropy involved [10]. Meanwhile, the injection of supercritical CO2 into a deep saline brine aquifer leads to complicated two-phase flow mechanisms in rock pore systems that are mainly dependent on pore space, pore size distribution and nondimensional parameters such as wettability, viscosity ratio and capillary number, etc. [28,33,34,35]. Even though the numerous traditional experimental methods [30,31,32,36] at core scale were used to study the constitutive relations such as capillary pressure and relative permeability curves, which can then be applied into large-scale continuum-based CFD simulators [37,38,39], the implementation of these experiments is highly complicated, expensive, and time-costly, with great measurement uncertainties that are unquantifiable. Moreover, due to the heterogeneity of a sizeable number of rock samples we encounter in the reality, simply treating the rock sample as a continuum material with uniform macroscopic petrophysical properties is a less effective method for uncovering various micro-scale mechanisms of fluid flow in the pore structures. Given the challenges for experimental studies and continuum-based simulation methods, the pore-scale numerical simulation of fluid flow in porous media is increasingly preferable for understanding phase behaviors and fluid–solid and fluid–fluid interactions in both quantitative and visually intuitive ways. Over the last few years, pore-scale numerical investigations have gained great popularity, especially with respect to their applications in brine aquifers. New techniques, such as X-ray CT scan and imaging, have gradually been applied to pore-scale and REV (representative elementary volume)-scale visualizations of multiphase flow in porous media [40,41,42,43,44].
The numerical methods of simulating flow behaviors in porous media and calculating relevant flow properties can be categorized into direct modeling and pore network modeling methods [44]. For the latter methods, a network representing the rock pore system is extracted and reconstructed based on principles of topology [45,46], and this is then utilized to simulate the flow pattern of single-phase and multiphase flow in it by solving the fluid dynamical governing equations. Although this approach features many remarkable strengths, such as fast computational speeds and finer resolution that satisfy any specific demands for engineering projects, the geometric oversimplification of its pore network with nodes and bonds can cause large, challenging errors in the simulation results [47,48]. The direct modeling approach, by contrast, solves the transport governing equations based on real digital images (micro-CT [40], SEM [49], i.e.,) of rocks, ensuring the physical fidelity of computational results and more importantly that the rock sample itself is not damaged. As a result, one can directly simulate two-phase flow in these images, which renders them better models better than those unnecessarily developed by numerical methods such as the Gaussian field method [50], simulated annealing algorithm [51], Markov chain Monte Carlo [52], etc. For such pore-scale research, several CFD methods have been extensively applied through numerically solving the Navier–Stokes or NS equation that guarantees the momentum conservation, such as the lattice Boltzmann method (LBM) [27], finite element method (FEM) [53], finite volume method (FVM) [54] and finite difference method (FDM) [55]. The NS equation is a partial differential equation of the second order that requires intricate discretization and, therefore, results in demanding programming [56]. In addition, the nonlinear convective term in the NS equation presents numerous numerical challenges [57,58]. Therefore, solving the NS equation in complex porous media, like characterized by highly heterogeneous rocks, either requires simplifying the medium’s geometrical characteristics or using complex grids that necessitate extensive computations. In comparison with other methods, LBM solves the discrete Boltzmann equation with a specific collision model rather than the conventional NS equation. Moreover, it can be effectively applied into pore-scale modelings by associating itself with physics involved in a series of interactions in the pore space. At the same time, it also easily handles complex pore geometries by voxelizing them into discrete nodes that match the lattices in its numerical scheme [58]. With the rapid evolution of computational power and platforms represented by graphical processing unit (GPU)-enhanced parallel computations [59], the large-scale LB-based simulation of multiphase flow in 3D complicated pore space has been realized and also empowered with finer resolution, higher accuracy of results and faster running speed.
To the best of our knowledge, two-phase CO2–brine flow has been extensively reported using the lattice Boltzmann method with its great advantages of modeling irregular flow pathways in complex pore geometries and coupling both fluid–fluid and fluid–solid interactions. Common two-phase LB models are composed of the color–gradient model [60], the pseudopotential model [61], the free energy model [62] and the phase-field model [27]. Though all of these types of LB models have been extensively utilized to simulate multiphase immiscible flow in porous media at pore scale, the color–gradient model and the free energy model show up better performance in handling high viscosity ratios and maintaining a sharp interface between two immiscible fluids [63]. Nevertheless, the color–gradient is more preferred for simulating immiscible multiphase flow in heterogeneous porous media under scenarios of low capillary numbers with good accuracy and efficiency [9]. Most of the literature [28,33,34,48] adopting the color LB model for modeling two-phase CO2–brine flow in rocks is based on pore-scale models, which needs detailed geometrical information and is characterized by a computation domain that is very small and thus inapplicable to field-scale practices [64]. An alternative method to it is the representative elementary volume (REV)-scale method, which is much larger than the pore scale and smaller than the field scale [45]. One of the greatest advantages of the REV-scale method is that some inessential pore structures at the micro scale can be overlooked and statistical flow properties such as permeability, velocity and pressure can be obtained through semi-empirical relations such as the Darcy’s law [65] and Brinkman-extended Darcy models [66], etc.

2. Motivation and Organization for This Study

As far as we are concerned, most of the previous studies regarding micro-scale LB simulation of two-phase CO2–brine flow in sandstone rocks mainly focused on pore-scale insights into the mechanism and flow properties such as relative permeability of the flow processes [19,21,25,26]. Few of them, however, made efforts into quantitatively considering CO2 storage efficiency due to the two-phase immiscible flow at REV scale, especially the imbibition processes in the complex pore space. Meanwhile, some of other studies [50,51,52] in the literature, which apply the generally simplified 2D or 3D pore models generated by computer modeling, can hardly consider the complexity and heterogeneity of sandstone rocks themselves as porous media. Accordingly, we are motivated in this study to understand the CO2 storage efficiency quantitatively due to the two-phase imbibition processes in the heterogeneous sandstones, as illustrated in Figure 1, at REV scale.
In this work, we aim at a micro-scale understanding of CO2 storage efficiency in a tighter REV extracted from a permeable heterogeneous sandstone rock sample based on LB simulation of steady-state CO2–brine imbibition processes, since the tighter subdomain is assumed to have the relative higher storage efficiency with higher capillary threshold for brine breakthrough. Firstly, REV analysis was utilized to extract 3 REV representing the upper, middle, and lower space of the whole core in which pore space was reconstructed using a series of image processing sections, including image preprocessing, filtering, and threshold segmentation. In addition, single-phase LB simulations were conducted on each of these 3 REVs and their calculated results including petrophysical properties and physical fields (pressure; velocity) were further analyzed in order to identify a tighter REV that its central part was then cut as prepared as boundary conditions for subsequent two-phase LB simulations of CO2–brine imbibition flow in it. Two main parts were involved: 1. Modellings were conducted on the tighter subsample under 4 differential pressures (100 Pa; 200 Pa; 300 Pa; 400 Pa) to identify if threshold a pressure gradient [67,68,69,70,71] below which no flow in the computational domain is driven exists and find out the optimized value under which the pore space is fully channelized as characteristics of Darcy flow.2. We conducted LB simulation of the steady-state CO2–brine imbibition processes over a span of target saturations on the tighter REV under 18 scenarios with different combinations of nondimensional parameters including wettability index, capillary number and viscosity ratio under the optimized differential pressure previously obtained to reveal relative permeability curves of CO2 under different combinations of three nondimensional parameters. On this basis, we further evaluated CO2 storage efficiency in these different cases by estimating its residual saturation through setting a lower bound of krCO2 = 0.05.
The organization of the remaining part of this paper is as follows: the next section is mainly concerned with the methodologies involved in this study. In Section 3.2, methodology for digital rock reconstruction and pore space extraction through image scanning and processing techniques are introduced and single-phase LB model, multiphase color–gradient LB model and the derived steady-state LBPM (lattice Boltzmann porous method) are, respectively interpreted in Section 3.2, Section 3.3, and Section 3.4. Results of and discussions on single-phase, two-phase LB simulations and CO2 storage efficiency, based on relative permeability under conditions of various combinations of nondimensional parameters, are analyzed and presented in Section 4, while the main results of the paper are summarized in Section 5 as conclusions. Besides, Section 6 or the last section has been especially created for this paper to demonstrate the expectations and limitations of this study for future work.

3. Methodology

3.1. Digital Rock Reconstruction and REV Analysis

The sample rock in this study is a sandstone extracted from the Daqing Oilfield at a depth of approximately 2047 m. Its porosity and permeability have been measured to be 14.46% and 450 mD, respectively. The main mineralogy of the rock was found to be quartz (75.714%), feldspar (14.65%), and other minerals including clay, calcite, dolomite, and pyrite were also present (less than 10% combined). A X-ray micro-CT imaging device, manufactured by Zeiss GMBH, based in Oberkochen, Germany, at RIPED (Research Institute of Petroleum Exploration and Development), Beijing, was used to scan the rock sample, generating a series of gray-scale scans like Figure 2, in which nearly no fracture was present and pore space is well preserved. These resemble the petrophysical properties of Berea sandstone [28,72]. Some brighter or darker patches are seen in these CT scans, like Figure 2, because various minerals in our rock sample have different densities and then feature different attenuation coefficients for X rays. However, the mineralogical distribution within our sample would not impact our simulation in the following part. A stack of 1400 CT scans like Figure 2 with size of 1998 × 2006 pixels and resolution of 1.98 μm is reconstructed as a 3D digital rock with dimension of 1998 × 2006 × 1400 voxels. To the best of our best knowledge, running efficient LB simulations of fluid flow in porous media like rocks requires graphics processing unit (GPU)-based parallel computing technologies, instead of a computing processing unit (CPU). Therefore, conducting LB simulations of the digital rock image of the whole-core rock sample is hardly realized on personal computers (PC) and, even with HPCC, it would inevitably result in big time consumption [36,64]. As a result, a careful study is needed to identify the size of a representative elementary volume (REV) for different areas in a whole digital rock that is prepared for flow simulation. In order to demonstrate the properties of different regions of the 3D digital rock images, REV analysis in terms of porosity, which is calculated through a series of image processing sections including filtering, threshold-segmentation (Figure 3), combined with central cutting method [73,74], as illustrated in Figure 4, is conducted on the three subdomains (upper, middle, and bottom) of the whole-core computational domain as shown in Figure 5. Figure 6 demonstrates that spatial porosity correlations with REV lengths for REV1 and REV2 converge at 300 voxels, and the convergence is at 350 voxels for REV3. These conclusions are based on the approximation criterion that the porosity difference varies by less than 2% with the REV length since the convergence length. To follow up on the following flow simulations where a tighter REV was needed, we decided to establish 3 REVs with size of 800 × 800 × 400 voxels (Figure 7) in order to sufficiently overlap a tighter REV in the whole-core domain of the rock sample. The calculated porosities for REV1, REV2 and REV3 are 10.654%, 8.642% and 16.721%.

3.2. Single-Phase LB Model

As we all know, the main role that lattice Boltzmann method plays is linking micro-scale mechanisms and macro-scale flow properties possessed as a mesoscopic media such as pressure and velocity [75]. Specifically, the whole computational domain is discretized into many of lattices containing conceptualized groups of fluid particles. In fact, their physical transport is governed by mass and momentum conservative equations [76,77,78]. In addition, the fluid-solid, fluid-fluid interactions can be well characterized using the LB method by utilizing the forces balance between groups of molecules, rather than tracking individual molecules [79] and the Navier–Stokes equation can then be recovered by application of a Taylor series expansion of the lattice Boltzmann equation, followed by the Chapmann–Enskog expansion [80]. The particle distribution function fi (x,t) at position x and time t in the ith velocity direction is assigned at each lattice node. In this chapter, the single-phase LB model is established as follows:
According to the LB equation with the Bhatnagar–Gross–Krook (BGK) approximation, the evolution of fi (x,t) can be defined as follows:
f i ( x + e i Δ t , t + Δ t ) f i ( x , t ) = 1 τ [ f i ( x , t ) f i e q ( x , t ) ]
where fi(x,t) and feq i (x,t) represent the non-equilibrium and equilibrium distribution functions, respectively; Δt is a timestep, τ is the relaxation time as a good measure of how quickly the flow system evolves from non-equilibrium to equilibrium state, related to kinematic viscosity (ν) and written as τ =ν/(c2 s × Δt) + 0.5 with sound speed cs = 1/√3.
In general, the lattice Boltzmann equation is divided into two steps: streaming (Equation (2)) and collision (Equation (3)):
f i ( x + e i Δ t , t + Δ t ) = f t ( x , t + Δ t )
f t i ( x , t + Δ t ) = f i ( x , t ) 1 τ [ f i ( x , t ) f i e q ( x , t ) ]
In the streaming step, fluid particles start to move to neighboring nodes where collision then happens. This can lead the flow system into a dynamic state, that can be revealed on the macroscopic flow behavior until the equilibrium state is achieved in the whole system where equilibrium distribution function feq i (x,t) can be defined as Equation (4), in which wi is the weight of the ith velocity component depending on the selected lattice model, such as D3Q15, D3Q17 and D3Q19 [75,77,78].
f i e q = ω i ρ 1 + 1 c s 2 ( e i u ) + 1 2 c s 4 ( e i u ) 2 1 2 c s 4 u 2
The most widely used D3Q19 (Figure 8) lattice model for 3D flow simulation is adopted in this study, and the lattice discrete velocity ei in D3Q19 lattice model is expressed as follows:
e i = [ 0 1 1 0 0 0 0 1 1 1 1 0 0 1 1 1 1 0 0 0 0 0 1 1 0 0 1 1 0 0 1 1 1 1 0 0 1 1 0 0 0 0 0 1 1 0 0 1 1 1 1 0 0 1 1 1 1 ]
Besides, the weight wi is designated for i = 0, w0 = 1/3, for i = 1–6, wi = 1/18, and for i = 7–18, wi = 1/36.
In order to adjust use of the LB model to fluid flow in porous media such as rocks in this study, one must incorporate porosity 𝜖 into the equilibrium equation (Equation (4)), which can be rewritten as Equation (6). Besides, the external force Fi (x,t) in the ith direction of the lattice model is necessarily considered in the LB equation by modifying Equation (1) into Equation (7). Fi (x,t) can be expressed in Equation 8, where F refers to the total external force of a single predetermined lattice node which is predefined manually and I indicates the identity matrix.
f i e q ( x , t ) = ρ ω i 1 + e i u c s 2 + e i u 2 ϵ c s 4 + u u 2 ϵ c s 2
f i x + e i Δ t , t + δ t f i ( x , t ) = 1 τ f i ( x , t ) f i e q ( x , t ) + Δ t F i ( x , t )
F i ( x , t ) = ρ ω i ( 1 0.5 τ ) e i F c s 2 + ( F u ) e i e i c s 2 I 2 ϕ c s 4
As a result, the macroscopic flow properties such as density, velocity and pressure and can be retrieved through the following equations:
ρ = i f i ( x , t ) , u = i f i ( x , t ) e i ρ + Δ t 2 F , p = ρ c s 2
Furthermore, the Navier–Stokes equation is recovered through the Chapmann–Enskog expansion into the following equation:
ρ t + ( ρ u ) = 0 ρ t + ( ρ u u ) = p + ρ v u + ( u ) T + F
In the end, the Darcy’ law is used to calculate the permeability of rocks:
K = ν < u > Δ p L F
where K is the permeability, <u> is the average velocity of seepage flow system when steady state is realized and ∆p/L denotes the pressure gradient along the flow direction.
Based on the established porosity-adaptive single-phase LB model above, single-phase flow simulation is run on the 3 REVs extracted in Section 3.1 and each of these subdomains is discretized into 800 × 800 × 400 lattices, meaning that one voxel in the 3D rock CT image is equivalent to the size of one lattice in the LB model. The fluid density and viscosity are initialized the same as those of CO2-saturated brine where ρ = 1100 kg/m3 and ν = 0.0011 Pa.s, according to [20]. The pressures for inlet and outlet along the flow direction (upward) are set to 100 Pa and 50 Pa, respectively, or differential pressure is 50 Pa. The main idea of this part of work is to extract a tighter REV from these three based on simulated pressure and velocity fields and calculated permeability for each of them, with relevant results shown in the following Section 4.1.

3.3. Multiphase Color LB Model

The multiphase multi-relaxation time (MRT) color LB model, firstly proposed by Gustensen et al. [81], is one of the most popular multiphase LB models [27,62,78] that can be implemented on complex boundary conditions like pore geometries. In the multiphase color LB model, colors are set to indicate different phases of the flow system where red in general refers to the wetting phase, while the nonwetting phase is indicated by blue.
This work is dependent upon a multiphase color LB model customized for two-phase CO2–brine flow in rocks, in which separate LB equations are utilized to keep track of the mass transport of each phase and total momentum of both. Unlike the single-phase model, two discrete velocity sets are separately considered for the physical evolution of two phases where the D3Q19 lattice model is adopted for momentum transport LB equations, whereas the mass transport LB equations are dependent upon the D3Q7 model, which is aligned with the first seven discrete velocity components of the D3Q19 lattice model. The densities for both fluid phases, ρnw and ρw (CO2 and brine in this case), are solved through mass transport LB equations, respectively formulated, and provide essential information on the distribution of two phases and the location of their interface. Under this circumstance, nonwetting or wetting saturation Snw and Sw can be calculated from the densities of both phases and is introduced and expressed in Equation (12):
S n w = ρ n w ρ n w + ρ w S w = ρ w ρ n w + ρ w
Additionally, Snw can also be related to a phase indicator denoted as φ by the following equation:
φ = 2 S n w 1
Since the maintenance of the two-phase interface is critical to the LB implementation and the diffusion between both fluids is an important mechanism of it, stresses between both phases may occur as a result of interfacial tension and should therefore be incorporated into the momentum transport equation. Logically, the color gradient is defined as follows:
C = φ
In typical LB implementations, the color gradient at the interface can be normalized using the following equation:
n = C C
Additionally, then the D3Q7-based mass transport equations for both phases (wetting and nonwetting) can be established as in Equation (16) based on their momentary densities and velocities until the equilibrium state is matched, where the flux through the interface is minimized and the distributions of both densities and velocities are stabilized.
g n w q = ω q ρ n 1 + e q · u + β e q · n ρ w ρ n / ρ w + ρ n g w q = ω q ρ w 1 + e q · u β e q · n ρ w ρ n / ρ w + ρ n
where β is a parameter controlling the thickness of the interface and gnwq and gwq are defined as the mass transport distribution of nonwetting and wetting phases, respectively, and the weights assigned in the D3Q7 lattice model are w0 = 1/3, for q = 1 − 6, wq = 1/9. Besides, the velocities and time evolution of densities for two phases are determined as follows:
u i ( x , t ) = q   =   0 6 g i q ( x , t ) · e q ρ i ( x , t ) ρ i ( x , t + Δ t ) = q   =   0 6 g i q x e q Δ t , t ( i = w , n w )
A special reminder is given here that a streaming step is directly applied into the equilibrium distribution functions without relaxation time like in the single-phase LB formulation.
At the same time, the momentum transport equations for both phases are characterized by a D3Q19-based MRT-LB scheme:
f q x + e q Δ t , t + Δ t f q ( x , t ) = M q i 1 Λ i i m i e q m i + ω q e q · F c s 2
where q, i = 0, 1, 2, 3, …, 18; and the fq(x,t) is the distribution function of fluid particles, which illustrates the probability that particles at position x travel along the qth direction during the collision step; in the streaming step, fq (x + eqt, t + ∆t) is the particle evolutionary distribution that measures how a given lattice node x influences the neighboring ones over a time step ∆t. In the MRT-LB numerical scheme, Λii is designated as a relaxation diagonal matrix (see [5]) governing the relaxation processes of a total of 19 moment components in the D3Q19 lattice structure and is represented as a function of averaged fluid dynamic viscosity ν, which is calculated from those of both non-wetting and wetting phases (Equation (19)). The 19 moment components are determined through a linear transformation: mi = Mqi fq. Meanwhile, the equilibrium moment components denoted as mieq are reliant on the total flow velocity u and color gradient C. Ultimately, the macroscopic total averaged density and velocity for both fluids can be obtained by Equation (20).
u = ( 1 + φ ) ν n w 2 + ( 1 φ ) ν w 2
ρ = q   =   0 18 f q u = q   =   0 18 f q e q ρ
Moreover, the total fluid pressure p at a given position is calculated from sound speed and total density denoted, given as cs (cs = 1/√3) and ρ, respectively, using the formula p = ρ c s . 2
The non-equilibrium moment components of 19 in total are the following ones:
m 1 e q = j x 2 + j y 2 + j z 2 + α | C | m 9 e q = 2 j x 2 j y 2 + j z 2 + α | C | 2 2 n x 2 n y 2 n z 2 m 11 e q = j y 2 j z 2 + α | C | 2 n y 2 n z 2 m 13 e q = j z j y + α | C | 2 n z n y m 14 e q = j x j y + α | C | 2 n x n y m 15 e q = j x j z + α | C | 2 n x n z
where jx, jy, jz represent the magnitude of components of the total momentum j (j = ρu) along three coordinate axes: x, y, z under Cartesian coordinate system and similarly nx, ny and nz are three components along them as to the normalized color gradient n. α is linearly related to interfacial tension between two phases: σnw, n = 6α. For further information on the concrete formulations of Mqi, m i e q and Λii, please refer to the literature, including but not limited to [5,74,81].

3.4. Two-Phase Flow Parameters and Steady-State LBPM

Relative permeability is the most critical macroscopic flow parameter for two-phase flow in porous media, being a good indication of the flow capacity of either fluid in the multiphase scenarios compared with the single-phase one. Factors controlling it include many rock and fluid properties that can be summarized by nondimensional parameters. Nondimensional flow parameters including capillary number (Ca), viscosity ratio (M) and wettability index (WI), being the first-order or most important ones, play a highly critical role in two-phase flow regimes in complex pore geometries as boundary conditions and therefore have an indispensable impact on macroscopic flow properties like fluid pressure and velocity, ultimately resulting in variation of relative permeability of both nonwetting and wetting fluids.
Regarding multiphase flow in porous media, a multiphase expansion of Darcy’s law (Equation (11)) is widely used to characterize two-phase CO2–brine flow in rocks at reservoir scale. The wetting phase has a stronger affinity to the rock surfaces than nonwetting phase, resulting in the emergence of different flow patterns between them due to thermodynamics. As such, the competitive passage of both fluids through the rocks occurs and the two-phase modified Darcy’s law is constructed to interpret this effect, as in the following equation:
u i = k r i ν i K ( p ρ i g )
where ui (i = nw, w) is the velocity of the nonwetting or wetting fluid phase; kri refers to the relative permeability of both phases; νi indicates their dynamic viscosities and K is the absolute permeability of the rock sample. In addition, ∇p and ρig hereby symbolize pressure gradient along the flow direction, as in the traditional Darcy’s law and the body force of the nonwetting or wetting phase, which is also part of external force applied into the flow system.
Wettability makes a great contribution to the fluid distribution, fluid–fluid interaction and flow pattern in the pore space. It is a good measure of the tendency of either fluid of both phases to compete to stick to and spread across the rock surface. Contact angle θ can be measured at the liquid–solid contact interface in the experiments as shown in Figure 9 and calculated from interfacial tension σnw,n and surface tensions for both fluids (σnw,s and σw,s) based on the Young’s equation [67], as in Equation (23):
cos θ = σ w , s σ n w , s σ n w , n
Therefore, WI can be expressed by the cosine of the contact angle at the rock–fluid interface, which lies in the range between −1 and 1. By referencing [82], wetting conditions for the CO2–brine–sandstone fluid–fluid–solid system in this paper are defined as the following criteria based on WI: WI = 1 or −1 represents the fact that the whole rock system is either completely water-wet or CO2-wet; it is strongly water-wet when WI is between 0.8 and 1, whilst strongly CO2-wet for −1 to 0.8; intermediately water-wet for 0.5 to 0.8 and intermediately CO2-wet for −0.8 to −0.5; weakly water-wet for 0 to 0.5 and weakly CO2-wet for −0.5 to 0. Since extreme cases, such as being completely water or CO2-wet, are seldom seen in real reservoir conditions, most of which, additionally, are water-wet to various extents, therefore only the effects of strongly water-wet, intermediately water-wet, and weakly water-wet scenarios are reported in this paper.
Capillary number or Ca is another essential nondimensional flow parameter that characterizes the domination between viscous and capillary forces in the flow behavior of both phases. It can be defined as the ratio of viscous force and capillary force. The viscous force dominates the flow pattern and enhances the mobility of both wetting and nonwetting fluids at high Ca. On the contrary, when it comes to a low capillary number, which means that there is a more competitive effect of capillary force over the flow system than that of the viscous force, both two phases will reduce their movable contents and be blocked into some local space because of the capillary resistance effect. The formula expressing the capillary number Ca is the following one:
C a =   ν n w u σ n w , w
where νnw is the dynamic viscosity of the nonwetting phase; u is the total velocity and σnw,n is the interfacial tension between two phases.
The viscosity ratio M is the third important controlling factor, except for the two other aforementioned ones for two-phase flow in porous media, and is defined as the ratio of dynamic viscosities between wetting and nonwetting phases as shown in Equation (25). Additionally, its importance level with respect to flow regimes in the pore space is dependent on the capillary number, where at high capillary number (Ca > 10−4, e.g.,), the viscosity ratio has a significant impact on the results, and the smaller contrast is between the viscosities of the two phases, the more favorable flow is for both of them; whereas, there is little or no effect on their flow patterns at low capillary number (Ca < 10−4, e.g.,), since under this condition either wetting or nonwetting phase is blocked due to the capillary effect [5,28,34,82].
M =   ν w   ν n w
The LB method is one of the most suitable numerical methods for digital rock physics represented by multiphase flow processes in rocks, and is used to efficiently deal with complex boundary conditions as pore geometries without losing its characteristics. Meanwhile, it is a mesoscopic method that contains micro-scale physical mechanisms and relates them to nondimensional parameters. For example, local surface wettability can be initialized based on laws accounting for non-equilibrium and equilibrium effects [5,82]. In general cases, the color LB model is the most widely used LB model in simulating immiscible two-phase fluid flow in complex porous media and stays stable for running over a wide range of nondimensional parameters: 10−6 < Ca < 1; 0.01 < M < 100 and, the model is adept at coping with the cases with big density contrast (larger than 100 between two phases).
Experiments for steady-state two-phase flow in rocks are popular for measuring relative permeability in the laboratory by monitoring the pressure drop along the rock sample at a given saturation of one fluid. However, such experiments fail to create rigorous steady-state flow due to spontaneous uncertainties in the fluid saturation [25]. By contrast, LB numerical simulations show uncompetitive strengths by imposing a periodic boundary at both inlet and outlet to secure a constant saturation of one fluid. A color LB model based on the use of digital rock images to simulate the steady-state two-phase flow in porous media is proposed by [5] and is also called steady-state LBPM (lattice Boltzmann models for porous media). Simply put, the successful implementation on the steady-state LBPM relies on two components: 1. a mechanism to drive and supervise the flow pattern of fluids at a constant saturation 2. a mechanism to evolve the flow configuration with changing the target saturation. Notably, during the steady-state LBPM simulation of CO2–brine immiscible flow in digital rock images, both phases are initialized at a target saturation based on morphology adaption [5,28] in the pore space, and the CO2 as the nonwetting phase is typically assigned in the large pores.
In this paper, a tighter 200 × 200 × 200 voxels REV (Figure 10a) with porosity and permeability of, respectively, 8.45% and 16.52 mD, was extracted at the center of REV1 (Figure 10b), and this was significantly less permeable than other the two as demonstrated in the previous part. To follow up, the steady-state LBPM is used to simulate two-phase CO2–brine immiscible flow in this tighter REV to study effects of nondimensional parameters on relative permeability with optimized differential pressure, after which we discussed flow pattern and capacity of CO2 under different cases and gave insights into CO2 storage efficiency based on estimated CO2 residual saturation (SresidualCO2) as discussed in Section 4.3.

4. Results and Discussions

4.1. Single-Phase LB Simulation and Extraction of Tighter REV

Single-phase LB simulations on three REVs (REV1, REV2 and REV3), extracted from the whole-core digital rock image of the rock sample, were conducted based on the model initial inputs as illustrated in Section 3.2 of this paper. Calculated permeabilities of three REVs are 50 mD, 532 mD and 900 mD, seeing big contrast between each other, despite originating from the same macroscopically permeable sandstone core sample. Besides, pressure (Figure 11a–c) and velocity fields (Figure 11d–e) were also visualized after the simulations, which shows that REV1 is remarkably more abnormal than the other two ones where no significant pressure drop is observed in the pressure field and discontinuous streamlines exist across the global domain, with extremely low flow velocity near zero at the outlet. The pressures and velocities are converted into the dimensionless form by dividing the value of a simulated sample at one certain position by the maximum value within the computational domain as the characteristic value under all the scenarios.
Based on the permeability calculations and results from the physical (pressure, velocity) fields of these three REVs, we suggest using a local domain inside REV1 with low permeability to account for the simulation results on the global-field scale. From this point of view, we determined to cut a tighter REV at the center of REV1 with a size of 200 × 200 × 200 voxels from which pore structure is characterized in Figure 10b. Figure 10a points out the cross section of the tighter 200 × 200 × 200 voxels REV extracted at the center of REV1. What can be seen from its pore structure shown in Figure 10b is that plane porosities on bounded facets of the cubic computational domain are limited and that only one connected pore–throat system is across the microstructure. It is then speculated that poor injectivity of fluids into the tighter pore space may cause complicated two-phase CO2–brine flow behavior in it, which directly has an impact on CO2 storage efficiency. Furthermore, we surmise that the effect of this kind of micro-scale CO2 storage mechanism may be upscaled in many of permeable deep saline sandstone aquifers at core scale or even reservoir scale. Based on this inference, a series of simulations of two-phase CO2–brine flow in this tighter REV under different circumstances including various differential pressures and combinations of nondimensional parameters introduced above using the steady-state LBPM. Based on the simulated results obtained under these scenarios, some insights are given into the CO2 storage efficiency in this microstructure.

4.2. Threshold Pressure Gradient and Optimization of Differential Pressure

There exists nearly zero-velocity non-Darcy flow in porous media with low permeability or pore connectivity under a certain pressure gradient across the flow domain which is called as threshold pressure gradient. Besides, non-Darcy low-velocity flow features another interval between the threshold pressure gradient and a point from which the Darcy flow starts [68]. Threshold pressure gradient plays an important role in flow behaviors in ultra-tight reservoirs at the macroscopic scale, which is correlated with many engineering applications such as enhanced oil recovery, ground water discharge and natural gas or hydrate production in these types of reservoirs [68,69]. However, so few reports in the literature were made on the micro-scale effects of threshold pressure gradient on two-phase CO2–brine flow in tight sandstone pore structures. In this section, a couple of LB simulations of two-phase CO2–brine flow in the tighter 200 × 200 × 200 voxels REV were conducted with 4 differential pressures, Pdiff: 100 Pa, 200 Pa, 300 Pa and 400 Pa, and their pressure gradients were converted to 0.25 MPa/m, 0.5 MPa/m, 0.75 MPa/m, 1 MPa/m in order to see if threshold pressure gradient exists in the tighter REV and select the optimized differential pressure where two-phase flow has been fully channelized through the pore structure, characterized by stable and nearly equal flow rate at both inlet and outlet and continuous pressure drop along the flow direction, as typical of the Darcy flow. The reason for the parametric selection was based on the literature [69,70,71], where most of threshold pressure gradient values that are measured in core flooding experiments lie between 0.005 MPa/m and 1.5 MPa/m at the core scale.
Initializations for the steady-state LBPM include the following components: wettability index WI on the rock surface is set to be 0.7, which is typical of intermediately water-wet conditions; viscosity ratio between two phases is 10 where νnw and νw are 1 × 10−7 Pa.s and 1 × 10−6 Pa.s (CO2 is the nonwetting phase and brine is the wetting phase); density for CO2 ρCO2 is 1000 kg/m3 at the supercritical state and brine density ρBrine is designated at 1100 kg/m3 when CO2 is saturated; capillary number Ca in this part is defined as 1 × 10−4; the tighter REV is saturated with brine and CO2 at equal fractions: SCO2 and SBrine are the same as 50%;
Figure 12a–d are simulated results of flow velocity distribution within the domain of tighter REV under the above 4 scenarios of various differential pressures (100 Pa, 200 Pa, 300 Pa, 400 Pa), likewise, the velocity here is converted to the dimensionless form as Figure 11 by dividing the value at any position by characteristic velocity—the greatest value encountered under four runs. These results indicate that for Pdiff =100 Pa (Figure 12a), there is nearly no fluid motion in the whole global domain, symbolizing that the flow system in this scenario is immobilized due to its differential pressure below the threshold one. Flow velocity localization was observed at the inlet when Pdiff reached 200 Pa (Figure 12b), which demonstrates that so few parts of fluid near the entrance under this circumstance have been activated, and the whole computational domain is not yet fully channelized. We therefore assume that the threshold pressure gradients that are frequently encountered in macroscopic experimental investigations also exists at a microscopic scale like in this this study, where threshold pressure gradient is in the vicinity of 0.5 MPa/m, which matches some of the experimental results reported in the literature [55,56,57]. Besides, it is notable in Figure 12c that the flow was fully channelized under Pdiff = 300 Pa and has approximately conformed to Darcy flow behaviors. This discovery can further imply that the linear Darcy flow was realized at a differential pressure P0 equal to or below Pdiff = 300 Pa, and that non-Darcy flow dominated the flow regime in the computational domain when the Pdiff lied in the interval between around 200 Pa and P0. Intuitively, the case with Pdiff = 400 Pa (Figure 12d) can give further proof that flow is not only channelized and flow rates at both sides are closer than Pdiff = 300 Pa. Figure 13 represents the pressure profiles of 4 scenarios with different pressures, from which one can tell that pressure profile keeps nearly constant along the flow direction, accompanied by a drastic plunge in the near-outlet region when differential pressure is 100 Pa and 200 Pa, in sharp contrast to which, significant and continuous pressure drop is observed at Pdiff = 300 Pa and 400 Pa, where the flow has been fully channelized, a typical characteristic of the Darcy flow behaviors; in addition, the linearity of pressure profile gets stronger as Pdiff evolves from 300 Pa to 400 Pa. Based on the demonstrated research results on the pressure profile of the tighter REV under 4 differential pressures, a couple of our insights are given as follows. Variation of pressure distribution along the flow direction in the non-Darcy flow interval (revealed by Figure 13b) and when differential pressure is lower than the threshold pressure (revealed by Figure 13a) is small, based on which we suggest that most of the fluid confined in the pore space can hardly be mobilized under the current differential pressure by overcoming the capillary resistance while the scope of pressure propagation covers the computational domain by equalizing the pressure distribution in it. Besides, the sharp plunge near the end is attributed to pressure continuity between predefined pressures at inlet and outlet from our humble perspective. Moreover, a pressure drop is significantly revealed with stronger linearity, especially as differential pressure gets larger after fluid flow enters the Darcy flow state (Figure 13c,d).
In summary, optimized differential pressure of 400 Pa was selected, at which flow is fully channelized with stably similar flow rates at both inlet and outlet and nearly linear pressure drop is established for the following simulations regarding sensitivity analysis for relative permeability under various combinations of nondimensional parameters in the following part of this paper.

4.3. CO2 Relative Permeability and Storage Efficiency

The two-phase CO2–brine imbibition processes are of great research significance to CO2 storage efficiency in pore spaces of sandstone reservoirs. As discussed before, most of pore surfaces in deep saline aquifers are water-wet to varying degrees. As such, in this section, the steady-state LBPM was utilized to simulate the CO2–brine imbibition processes over a range of target saturations (0.2–0.9) under a total of 18 scenarios (Table 1) constituting different combinations of 3 nondimensional parameters (Ca, M, WI) using the optimized differential pressure (400 Pa), with the highest injectivity of brine in the tighter REV in order to validate the systematic effects of nondimensional parameters on relative permeability for CO2 as the nonwetting phase and therefore discuss CO2 storage efficiency as the residual phase in the pore geometry under these 18 cases (Table 1).
In fact, the use of the steady-state LBPM to precisely and physically predict the CO2–brine imbibition processes, which are naturally unsteady flooding processes, requires an initialized fluid distribution where CO2 almost saturates the whole porosity and the near-inlet region is filled up with brine, as illustrated in Figure 14b, followed by simulations of CO2–brine imbibition processes over a span of target saturations. A better visualized example of the process is presented in Figure 15, from which one can easily see that model initial conditions match scenario 14 at Table 1. Results from Figure 15 show that the wetting-phase breakthrough happened along with brine rapidly dominating the main connected channel in the pore space until the global fluid distribution and flow pattern was stabilized as the target saturation reached a certain level, like SBrine = 0.5 in this case, where relative permeability of CO2 nearly approached zero and it was expelled into isolated big pores and stored as residual phase.
Figure 16 shows simulated relative permeability curves for CO2 (nonwetting) over a span of target saturations ranging from 0.2 to 0.9 under a total of 18 scenarios constituting various combinations of three nondimensional parameters. As a result, a systematic interpretation of effects of these three parameters on relative permeability of CO2 is made, followed by CO2 storage efficiency analysis based on these simulated relative permeability curves through setting up a lower bound (krCO2 = 0.05), below which the ability of CO2 to leave the pore space is infinitely small and can be assumed as zero. Here, in this study, the brine saturation at which the CO2 relative permeability firstly drops below 0.05 is denoted as SBrinemax. The approximated CO2 residual saturation (SresidualCO2) is therefore calculated by the following equation:
S r e s i u a l C O 2 = 1 S B r i n e m a x
As we discussed in the previous sections, capillary number is divided into high or low values based on a boundary value in between, and it is clear from Figure 16 that simulation results at Ca = 1 × 10−5 and 1 × 10−4 are not considerably distinguishable without varying the other two parameters. By contrast, relative permeability for CO2 is significantly higher at Ca = 1 × 10−3, when other conditions are controlled. This obvious comparison can validate the fact that the boundary value of Ca for our study is equal to 1 × 10−4 or lower but much higher than 1 × 10−3, beyond which capillary effect is enough to confine CO2 as nonwetting phase in the pores, especially larger ones, and further effects are obscure with increasing Ca. Additionally, the linearity for relative permeability of CO2 in the case with Ca = 1 × 10−3 is much stronger than those with Ca = 1 × 10−4 or 1 × 10−5 when M = 1, which illustrates that the flow of CO2 at low capillary number is totally dominated by the viscous force; besides, if possibly M is low enough as 1, like in Figure 16a, the CO2 flow pattern is much closer to characteristics of Poiseuille flow than other scenarios, given that brine spreads across the rocks surface as a water film that facilitates the flow of CO2. Nevertheless, the viscosity ratio also plays a crucial role in varying CO2 relative permeabilities at both high and low capillary numbers. A comparison of Figure 16a and Figure 16b offers us an essential piece of information that the slope or the rate of the decline of CO2 relative permeability curves is lower at higher M than at lower the one, while the overall relative permeabilities over the whole saturation range rise to some degree as M decreases, which is mainly caused, in our view, by the lubrication effect where the water film coating the rock surface shares similar viscosity to CO2 so that the viscous coupling effect is less apparent. However, at low a capillary number (Ca = 1 × 10−4, Ca = 1 × 10−5), as illustrated from Figure 16c,d, relative permeabilities are only affected by M for the initial saturation range in which the curves have yet to converge to zero, which means that because of the strong dominance of pore confinement at a low capillary number, the effect of viscosity ratio is hereby limited, with only the convergence rate being alleviated by a finite quantity. Finally, the effect of wettability on the relative permeability performance of CO2 is also discussed in both high- and low-capillary number scenarios, with M being either 1 or 10. Figure 16a and Figure 16b present that relative permeability of CO2 is insensitive to wettability of the rock surface at high capillary number since viscous force is the dominating force controlling the migration of fluids and the capillary force poses little or no resistance to the nonwetting phase flow. Moreover, the capillary force plays an essential role in manipulating the flow system at low capillary, therefore, the effect of wettability on CO2 relative permeability number seems to be magnified. It is suggested from Figure 16c–f that CO2 relative permeability is reduced with water-wetting conditions being weaker and that the effect is much more conspicuous when it transforms from intermediately water-wet to weakly water-wet; as brine saturation rises, the wettability effect is also reduced until none. To the best of our knowledge, residual CO2 is expelled to the larger pores coated by brine and the main porosity is channelized for brine transport through the porous media, with the global target saturation of brine approaching the maximum where no CO2 in the pores can further be pulled out of the computational domain.
The relative permeability of CO2 exerts a great impact on its storage efficiency, as illustrated in the literature, and the higher amount of ultimate safe storage capacity, which is presented as the residual phase in the pore space after the two-phase imbibition processes, means higher storage efficiency at relatively lower risk of leakage. In the rest of the paper, we attempt to analyze CO2 storage efficiency in different scenarios with varying combinations of nondimensional parameters in the tighter REV based on the relative permeability curves of CO2 by means of setting a lower bound value of 0.05 and identifying the target brine saturation with relative permeability that first breaks down the bound, denoted as the maximum brine saturation SBrinemax, which logically derives the CO2 residual saturation SresidualCO2 as expressed in Equation (26). Figure 17 gives a close-up view of relative permeability curves simulated in Figure 17 by narrowing down the brine saturation scope, where relative permeability values are generally below or approach the lower bound value as illustrated in each panel. Apparently, residual CO2 saturations are quite distinguishable between cases at high Ca (Figure 17a,b) and at low Ca (Figure 17c,d). One can easily observe that Figure 17a and Figure 17b both demonstrate much less residual saturation of CO2 than other ones. Estimated CO2 residual saturations are between 0.1 and 0.2 for Ca = 1 × 10−3, M = 10 and 0 and 0.1 for Ca = 1 × 10−3, M = 1, respectively, with relative smaller differences. However, at low Ca, all the cases have obviously lower CO2 residual saturation, with one estimated between 0.5 and 0.6 for Ca = 1 × 10−4 and 1 × 10−5, M = 10 and one between 0.4 and 0.5 for Ca = 1 × 10−4 and 1 × 10−5, M =1. Based on our analysis, since the whole computational domain is controlled by viscous force with rather less pore confinement effect on the flow of fluids at high Ca than low Ca, the flow of CO2 is much more favored with less residual CO2 left in the pores. Consequently, we can primarily conclude that capillary number has an essential role to play in deciding the CO2 storge efficiency of the two-phase CO2–brine flow system as Ca transcends the boundary between low and high Ca values or vice versa, or simply put, higher Ca results in lower CO2 storage efficiency. When the capillary number is controlled at low values (Ca = 1 × 10−4, 1 × 10−5), the viscosity ratio seems to influence the CO2 storage efficiency to some degree. Figure 17b, Figure 17d and Figure 17f present higher CO2 residual saturation with increasing M from 1 to 10, which suggests that decreased similarity in viscosity between CO2 and brine leads to a slight improvement in CO2 storage efficiency due to the mitigated lubrication effect. However, its effect is much less superior to capillary number, for which the pore confinement is responsible. Additionally, it is seen from Figure 17a–f that wettability poses no great changes to the ultimate estimated CO2 residual saturation when the other two parameters, Ca and M, are constant. Accordingly, we assert that the extent of water-wetting conditions of sandstone rock surfaces is less sensitive to CO2 storage efficiency, while it only exerts a somehow measurable influence on the CO2 relative permeabilities before the curves converge to zero at low Ca.

5. Conclusions

In this study, we used the steady-state multiphase color–gradient lattice Boltzmann (LB) model to simulate two-phase CO2–brine flow in a tight rock image volume (REV) extracted from a digital rock image of a whole-core permeable sandstone under four differential pressures (100 Pa, 200 Pa, 300 Pa, and 400 Pa). We identified a threshold pressure gradient less than or equal to 0.5 MPa/m (200 Pa), a non-Darcy flow interval between 200 Pa and a value smaller than 300 Pa, and Darcy flow behavior at Pdiff = 300 Pa and 400 Pa. The channelization of fluids was observed to improve as the differential pressure increased, with closer flow rates at the inlet and outlet. Pressure profiles along the flow direction were presented for the four differential pressures, demonstrating that the immobilized case (Pdiff = 100 Pa) and non-Darcy flow case (Pdiff = 200 Pa) both exhibited nearly constant pressure distributions, with an abrupt drop only seen near the outlet. In contrast, the Darcy flow interval showed a continuous and significant pressure drop across the domain, and the linearity of the pressure drop increased as the differential pressure increased from 300 Pa to 400 Pa. Based on these findings, we posit that differential pressure is a key factor influencing CO2–brine flow behavior in the tight REV, with a threshold pressure gradient. The more significant pressure drops and flow channelization in the REV, the more optimized the differential pressure under which the domain is governed by Darcy flow behavior will be. Therefore, Pdiff = 400 Pa was selected for the subsequent two-phase CO2–brine imbibition process over a range of target saturations in the tight REV.
Furthermore, LB simulations of CO2–brine imbibition processes over a predetermined range of target saturations under the optimized differential pressure (400 Pa) were conducted in the tight REV under a series of scenarios comprising 18 different combinations of nondimensional parameters (Ca, M and WI). These simulations were carried out in order to understand the dependencies of CO2 relative permeability and to analyze storage efficiency, as indicated by the estimated CO2 residual saturation in the pore space. A lower bound of krCO2 = 0.05 was applied to the simulated relative permeability curves. The conclusions in this regard were made based on the simulation results as follows:
  • When examining the effect of the capillary number (Ca) on CO2 relative permeability, it can be divided into high and low values. At high Ca, the CO2 relative permeability is globally higher and the relative permeability curve of CO2 is more linear, resulting in less restricted CO2 flow across the computational domain due to weaker pore confinement. This leads to significantly lower CO2 storage efficiency, as indicated by extremely smaller estimates of CO2 residual saturation compared to those at low Ca.
  • The viscosity ratio (M) is another parameter that influences CO2 relative permeability and facilitates CO2 flow in the pore space through the lubrication effect as it decreases. Lower M leads to more favorable CO2 flow in the imbibition process and thus less residual CO2 in the pores when other conditions are controlled, implying a negative impact on storage efficiency. However, the impact of the viscosity ratio on relative permeability is significantly less than that of the capillary number. At low capillary number, the lubrication effect is overpowered by the pore confinement effect in the relative permeability curves, while at high capillary number, the reduced viscosity contrast between the two phases may result in a finite improvement in CO2 flow capacity and lower storage efficiency. However, this difference is not considerably observable because the viscous force in the pores, already dominant compared to the capillary force, masks it.
  • Wettability only affects CO2 relative permeabilities in the initial saturation range before the curves converge to zero. Under stronger water-wet conditions, CO2 has an easier time passing through the domain due to less affinity to the surface, but its storage efficiency is relatively insensitive to changes in wettability when other conditions are controlled.
The systematic quantification of CO2 storage efficiency in a tighter REV in a permeable and heterogeneous sandstone rock sample gives micro-scale insights into how the mechanisms of two-phase CO2–brine imbibition processes relate to in situ stable CO2 storage.

6. Limitations and Expectations

This study has used the steady-state multiphase color–gradient LB model to investigate the effect of various nondimensional parameters on CO2 storage efficiency, as indicated by CO2 residual saturation, in a tighter REV with a threshold pressure gradient at the micro scale during the two-phase CO2–brine imbibition processes. While some new findings have been obtained, there is still room for improvement in both the LB numerical model and the physics of two-phase CO2–brine seepage flow itself. Some potential avenues for future work include: (1) the development of unsteady-state LB models that consider temporal effects to examine the dynamics of two-phase CO2–brine seepage flow over the continuous saturation range; (2) the incorporation of phase behaviors under pore confinement into the LB numerical scheme; and (3) the use of additional larger REVs or whole cores for modeling through mesh dependence analysis, if computational resources permit.
Although one lattice in the computational domain is represented by one voxel of 3D digital rock images in this study, we have to admit that mesh dependence analysis is necessary in many cases to be conducted on rock images before LB simulations for fine characterization of complex pore structures, fluid–fluid, and fluid–solid interactions, if computational resources are abundant. Furthermore, mesh dependence analysis is typically essential for multiphase flow until an optimum resolution is asserted for defining the effective interface thickness in lattice unit and reducing the spurious velocity magnitude of the LB implementations so as to capture as enough flow details at pore scale as possible. For further and detailed instructions on how to embark on mesh dependence analysis for LB simulations, one is encouraged to refer to [83].

Author Contributions

Conceptualization and methodology, Y.W., W.Z.; LB algorithms and implementation, Y.W. and W.Z.; Data curation and visualization, Y.W., W.Z. and L.J.; formal analysis and investigation, Y.W., W.Z. and Z.C.; writing the manuscript, Y.W. and W.Z.; funding acquisition, and project administration, C.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

This work is supported by CNPC Scientific Research and Technology Development Project “Whole petroleum system theory and unconventional hydrocarbon accumulation mechanism” (2021DJ0101).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Huang, L.; Zhou, W.; Xu, H.; Wang, L.; Zou, J.; Zhou, Q. Dynamic fluid states in organic-inorganic nanocomposite: Implications for shale gas recovery and CO2 sequestration. Chem. Eng. J. 2021, 411, 128423. [Google Scholar] [CrossRef]
  2. Pörtner, H.-O.; Roberts, D.C.; Adams, H.; Adler, C.; Aldunce, P.; Ali, E.; Begum, R.A.; Betts, R.; Kerr, R.B.; Biesbroek, R. Climate Change 2022: Impacts, Adaptation and Vulnerability; IPCC Sixth Assessment Report 2022; IPCC: Geneva, Switzerland, 2022. [Google Scholar]
  3. Matos, S.; Viardot, E.; Sovacool, B.K.; Geels, F.W.; Xiong, Y. Innovation and climate change: A review and introduction to the Special Issue. Technovation 2022, 117, 102612. [Google Scholar] [CrossRef]
  4. Liu, B.; Qin, J.; Shi, J.; Jiang, J.; Wu, X.; He, Z. New perspectives on utilization of CO2 sequestration technologies in cement-based materials. Constr. Build. Mater. 2021, 272, 121660. [Google Scholar] [CrossRef]
  5. McClure, J.E.; Li, Z.; Berrill, M.; Ramstad, T. The LBPM software package for simulating multiphase flow on digital images of porous rocks. Comput. Geosci. 2021, 25, 871–895. [Google Scholar] [CrossRef]
  6. Lichtschlag, A.; Pearce, C.R.; Suominen, M.; Blackford, J.; Borisov, S.M.; Bull, J.M.; de Beer, D.; Dean, M.; Esposito, M.; Flohr, A. Suitability analysis and revised strategies for marine environmental carbon capture and storage (CCS) monitoring. Int. J. Greenh. Gas Control 2021, 112, 103510. [Google Scholar] [CrossRef]
  7. Paltsev, S.; Morris, J.; Kheshgi, H.; Herzog, H. Hard-to-Abate Sectors: The role of industrial carbon capture and storage (CCS) in emission mitigation. Appl. Energy 2021, 300, 117322. [Google Scholar] [CrossRef]
  8. Mouahid, A.; Bombarda, I.; Claeys-Bruno, M.; Amat, S.; Myotte, E.; Nisteron, J.-P.; Crampon, C.; Badens, E. Supercritical CO2 extraction of Moroccan argan (Argania spinosa L.) oil: Extraction kinetics and solubility determination. J. CO2 Util. 2021, 46, 101458. [Google Scholar] [CrossRef]
  9. Bakhshian, S.; Hosseini, S.A. Pore–scale analysis of supercritical CO2–brine immiscible displacement under fractional–wettability conditions. Adv. Water Resour. 2019, 126, 96–107. [Google Scholar] [CrossRef]
  10. Ershadnia, R.; Hajirezaie, S.; Amooie, A.; Wallace, C.D.; Gershenzon, N.I.; Hosseini, S.A.; Sturmer, D.M.; Ritzi, R.W.; Soltanian, M.R. CO2 geological sequestration in multiscale heterogeneous aquifers: Effects of heterogeneity, connectivity, impurity, and hysteresis. Adv. Water Resour. 2021, 151, 103895. [Google Scholar] [CrossRef]
  11. Hou, L.; Yu, Z.; Luo, X.; Wu, S. Self-sealing of caprocks during CO2 geological sequestration. Energy 2022, 252, 124064. [Google Scholar] [CrossRef]
  12. Ali, M.; Al-Yaseri, A.; Awan, F.U.R.; Arif, M.; Keshavarz, A.; Iglauer, S. Effect of water-soluble organic acids on wettability of sandstone formations using streaming zeta potential and NMR techniques: Implications for CO2 geo-sequestration. Fuel 2022, 329, 125449. [Google Scholar] [CrossRef]
  13. Teng, Y.; Wang, P.; Xie, H.; Zhu, J. Capillary trapping characteristics of CO2 sequestration in fractured carbonate rock and sandstone using MRI. J. Nat. Gas Sci. Eng. 2022, 108, 104809. [Google Scholar] [CrossRef]
  14. Abd Rahman, I.; Hasbollah, D.A.; Yunus, N.M.; Kasiman, E.; Mazlan, A. In Carbon Dioxide Storage Potential in Malaysian Sandstone Aquifer: An Overview. IOP Conf. Series: Earth Environ. Sci. 2022, 971, 012022. [Google Scholar] [CrossRef]
  15. Gorecki, C.D.; Sorensen, J.A.; Bremer, J.M.; Knudsen, D.J.; Smith, S.A.; Steadman, E.N.; Harju, J.A. Development of storage coefficients for determining the effective CO2 storage resource in deep saline formations. In Proceedings of the SPE International Conference on CO2 Capture, Storage, and Utilization, San Diego, CA, USA, 2–4 November 2009; OnePetro: Richardson, TX, USA, 2009. [Google Scholar]
  16. Page, S.; Williamson, A.; Mason, I. Carbon capture and storage: Fundamental thermodynamics and current technology. Energy Policy 2009, 37, 3314–3324. [Google Scholar] [CrossRef]
  17. Jin, C.; Liu, L.; Li, Y.; Zeng, R. Capacity assessment of CO2 storage in deep saline aquifers by mineral trapping and the implications for Songliao Basin, Northeast China. Energy Sci. Eng. 2017, 5, 81–89. [Google Scholar] [CrossRef]
  18. Yanzhong, W.; Nianmin, Z.; Xu, C.; Yingchang, C.; Guanghui, Y.; Gluyas, J.G.; Miruo, L. Geologic CO2 storage in arkosic sandstones with CaCl2-rich formation water. Chem. Geol. 2020, 558, 119867. [Google Scholar] [CrossRef]
  19. Koukouzas, N.; Kypritidou, Z.; Purser, G.; Rochelle, C.A.; Vasilatos, C.; Tsoukalas, N. Assessment of the impact of CO2 storage in sandstone formations by experimental studies and geochemical modeling: The case of the Mesohellenic Trough, NW Greece. Int. J. Greenh. Gas Control 2018, 71, 116–132. [Google Scholar] [CrossRef]
  20. Alemu, B.L.; Aagaard, P.; Munz, I.A.; Skurtveit, E. Caprock interaction with CO2: A laboratory study of reactivity of shale with supercritical CO2 and brine. Appl. Geochem. 2011, 26, 1975–1989. [Google Scholar] [CrossRef]
  21. Huq, F.; Blum, P.; Marks, M.A.; Nowak, M.; Haderlein, S.B.; Grathwohl, P. Chemical changes in fluid composition due to CO2 injection in the Altmark gas field: Preliminary results from batch experiments. Environ. Earth Sci. 2012, 67, 385–394. [Google Scholar] [CrossRef]
  22. García-Rios, M.; Luquot, L.; Soler, J.M.; Cama, J. Laboratory-scale interaction between CO2-rich brine and reservoir rocks (limestone and sandstone). Procedia Earth Planet. Sci. 2013, 7, 109–112. [Google Scholar] [CrossRef] [Green Version]
  23. Liu, Q.; Zhu, D.; Jin, Z.; Tian, H.; Zhou, B.; Jiang, P.; Meng, Q.; Wu, X.; Xu, H.; Hu, T. Carbon capture and storage for long-term and safe sealing with constrained natural CO2 analogs. Renew. Sustain. Energy Rev. 2023, 171, 113000. [Google Scholar] [CrossRef]
  24. Ning, Y.; Tura, A. Wellbore Integrity Impact on Carbon Leakage to Ensure Safe Geological Sequestration. In Proceedings of the SPE Annual Technical Conference and Exhibition, Houston, TX, USA, 3–5 October 2022; OnePetro: Richardson, TX, USA, 2022. [Google Scholar]
  25. Ding, K.; Wang, L.; Ren, B.; Li, Z.; Wang, S.; Jiang, C. Experimental study on relative permeability characteristics for CO2 in sandstone under high temperature and overburden pressure. Minerals 2021, 11, 956. [Google Scholar] [CrossRef]
  26. Gholami, R.; Raza, A.; Iglauer, S. Leakage risk assessment of a CO2 storage site: A review. Earth-Sci. Rev. 2021, 223, 103849. [Google Scholar] [CrossRef]
  27. Moradi, B.; Hosseini Moghadam, A.; Rasaei, M.R.; Papi, A. Pore-scale characterization of CO2 front progress through a porous medium using a free energy model based on Phase-Field Lattice Boltzmann Method. J. Pet. Sci. Eng. 2022, 210, 110008. [Google Scholar] [CrossRef]
  28. Kohanpur, A.H.; Rahromostaqim, M.; Valocchi, A.J.; Sahimi, M. Two-phase flow of CO2-brine in a heterogeneous sandstone: Characterization of the rock and comparison of the lattice-Boltzmann, pore-network, and direct numerical simulation methods. Adv. Water Resour. 2020, 135, 103469. [Google Scholar] [CrossRef]
  29. Tahmasebi, P.; Sahimi, M.; Kohanpur, A.H.; Valocchi, A. Pore-scale simulation of flow of CO2 and brine in reconstructed and actual 3D rock cores. J. Pet. Sci. Eng. 2017, 155, 21–33. [Google Scholar] [CrossRef]
  30. Akbarabadi, M.; Piri, M. Relative permeability hysteresis and capillary trapping characteristics of supercritical CO2/brine systems: An experimental study at reservoir conditions. Adv. Water Resour. 2013, 52, 190–206. [Google Scholar] [CrossRef]
  31. Liu, F.; Lu, P.; Griffith, C.; Hedges, S.W.; Soong, Y.; Hellevang, H.; Zhu, C. CO2–brine–caprock interaction: Reactivity experiments on Eau Claire shale and a review of relevant literature. Int. J. Greenh. Gas Control 2012, 7, 153–167. [Google Scholar] [CrossRef]
  32. Tsakiroglou, C.; Avraam, D.; Payatakes, A. Transient and steady-state relative permeabilities from two-phase flow experiments in planar pore networks. Adv. Water Resour. 2007, 30, 1981–1992. [Google Scholar] [CrossRef]
  33. Suwandi, N.; Jiang, F.; Tsuji, T. Relative Permeability Variation Depending on Viscosity Ratio and Capillary Number. Water Resour. Res. 2022, 58, e2021WR031501. [Google Scholar] [CrossRef]
  34. Fan, M.; Dalton, L.E.; McClure, J.; Ripepi, N.; Westman, E.; Crandall, D.; Chen, C. Comprehensive study of the interactions between the critical dimensionless numbers associated with multiphase flow in 3D porous media. Fuel 2019, 252, 522–533. [Google Scholar] [CrossRef]
  35. Liu, Z.; Wu, H. Pore-scale modeling of immiscible two-phase flow in complex porous media. Appl. Therm. Eng. 2016, 93, 1394–1402. [Google Scholar] [CrossRef]
  36. Müller, N. Supercritical CO2-brine relative permeability experiments in reservoir rocks—Literature review and recommendations. Transp. Porous Media 2011, 87, 367–383. [Google Scholar] [CrossRef]
  37. Carrillo, F.J.; Bourg, I.C.; Soulaine, C. Multiphase flow modeling in multiscale porous media: An open-source micro-continuum approach. J. Comput. Phys. X 2020, 8, 100073. [Google Scholar] [CrossRef]
  38. Wang, J. Continuum theory for dense gas-solid flow: A state-of-the-art review. Chem. Eng. Sci. 2020, 215, 115428. [Google Scholar] [CrossRef]
  39. Papp, J.; Wilmoth, R.; Chartrand, C.; Dash, S. Simulation of high-altitude plume flow fields using a hybrid continuum CFD/DSMC approach. In Proceedings of the 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, Sacramento, CA, USA, 9–12 July 2006; p. 4412. [Google Scholar]
  40. Wang, S.; Ruspini, L.C.; Øren, P.E.; Van Offenwert, S.; Bultreys, T. Anchoring Multi-Scale Models to Micron-Scale Imaging of Multiphase Flow in Rocks. Water Resour. Res. 2022, 58, e2021WR030870. [Google Scholar] [CrossRef]
  41. Masalmeh, S.K.; Jing, X.; Roth, S.; Wang, C.; Dong, H.; Blunt, M. Towards predicting multi-phase flow in porous media using digital rock physics: Workflow to test the predictive capability of pore-scale modeling. In Proceedings of the Abu Dhabi International Petroleum Exhibition and Conference, Abu Dhabi, United Arab Emirates, 9–12 November 2015; OnePetro: Richardson, TX, USA, 2015. [Google Scholar]
  42. Dinariev, O.; Evseev, N. Multiphase flow modeling with density functional method. Comput. Geosci. 2016, 20, 835–856. [Google Scholar] [CrossRef]
  43. Koroteev, D.; Dinariev, O.; Evseev, N.; Klemin, D.; Nadeev, A.; Safonov, S.; Gurpinar, O.; Berg, S.; Van Kruijsdijk, C.; Armstrong, R. Direct hydrodynamic simulation of multiphase flow in porous rock. Petrophysics-SPWLA J. Form. Eval. Reserv. Descr. 2014, 55, 294–303. [Google Scholar]
  44. Blunt, M.J.; Bijeljic, B.; Dong, H.; Gharbi, O.; Iglauer, S.; Mostaghimi, P.; Paluszny, A.; Pentland, C. Pore-scale imaging and modelling. Adv. Water Resour. 2013, 51, 197–216. [Google Scholar] [CrossRef]
  45. Øren, P.-E.; Bakke, S. Reconstruction of Berea sandstone and pore-scale modelling of wettability effects. J. Pet. Sci. Eng. 2003, 39, 177–199. [Google Scholar] [CrossRef]
  46. Valvatne, P.H.; Blunt, M.J. Predictive pore-scale modeling of two-phase flow in mixed wet media. Water Resour. Res. 2004, 40, W07406. [Google Scholar] [CrossRef] [Green Version]
  47. Blunt, M.J.; Jackson, M.D.; Piri, M.; Valvatne, P.H. Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv. Water Resour. 2002, 25, 1069–1089. [Google Scholar] [CrossRef]
  48. Ramstad, T.; Øren, P.-E.; Bakke, S. Simulation of two-phase flow in reservoir rocks using a lattice Boltzmann method. SPE J. 2010, 15, 917–927. [Google Scholar] [CrossRef]
  49. Kelly, S.; El-Sobky, H.; Torres-Verdín, C.; Balhoff, M.T. Assessing the utility of FIB-SEM images for shale digital rock physics. Adv. Water Resour. 2016, 95, 302–316. [Google Scholar] [CrossRef]
  50. Liang, Z.; Fernandes, C.; Magnani, F.; Philippi, P. A reconstruction technique for three-dimensional porous media using image analysis and Fourier transforms. J. Pet. Sci. Eng. 1998, 21, 273–283. [Google Scholar] [CrossRef]
  51. Zhou, X.-P.; Xiao, N. 3D numerical reconstruction of porous sandstone using improved simulated annealing algorithms. Rock Mech. Rock Eng. 2018, 51, 2135–2151. [Google Scholar] [CrossRef]
  52. Lin, W.; Li, X.; Yang, Z.; Xiong, S.; Luo, Y.; Zhao, X. Modeling of 3D rock porous media by combining X-ray CT and Markov chain Monte Carlo. J. Energy Resour. Technol. 2020, 142, 013001. [Google Scholar] [CrossRef]
  53. Xie, Y.; Luo, Z.; Cheng, L.; Zhao, L.; Chen, X.; Zhang, N.; Ren, D.; Cao, Y. Numerical modeling and analysis of the matrix acidizing process in fractured sandstone rocks with the Extended–FEM. J. Pet. Sci. Eng. 2023, 220, 111215. [Google Scholar] [CrossRef]
  54. Huang, Y.; Zhou, Z.-f. Simulation of solute transport using a coupling model based on finite volume method in fractured rocks. J. Hydrodyn. 2010, 22, 129–136. [Google Scholar] [CrossRef]
  55. Alpak, F.; Samardžić, A.; Frank, F. A distributed parallel direct simulator for pore-scale two-phase flow on digital rock images using a finite difference implementation of the phase-field method. J. Pet. Sci. Eng. 2018, 166, 806–824. [Google Scholar] [CrossRef]
  56. Gresho, P.M. Some current CFD issues relevant to the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. 1991, 87, 201–252. [Google Scholar] [CrossRef]
  57. Karniadakis, G.E.; Israeli, M.; Orszag, S.A. High-order splitting methods for the incompressible Navier–Stokes equations. J. Comput. Phys. 1991, 97, 414–443. [Google Scholar] [CrossRef]
  58. Soleimani, R.; Norouzi, S.; Rasaei, M.R. Investigation of gas condensate drop-out effect on gas relative permeability by Lattice Boltzmann modelling. Can. J. Chem. Eng. 2019, 97, 1921–1930. [Google Scholar] [CrossRef]
  59. Xu, A.; Shi, L.; Zhao, T. Accelerated lattice Boltzmann simulation using GPU and OpenACC with data management. Int. J. Heat Mass Transf. 2017, 109, 577–588. [Google Scholar] [CrossRef]
  60. Huang, H.; Huang, J.-J.; Lu, X.-Y. Study of immiscible displacements in porous media using a color-gradient-based multiphase lattice Boltzmann method. Comput. Fluids 2014, 93, 164–172. [Google Scholar] [CrossRef]
  61. Chen, L.; Kang, Q.; Mu, Y.; He, Y.-L.; Tao, W.-Q. A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications. Int. J. Heat Mass Transf. 2014, 76, 210–236. [Google Scholar] [CrossRef]
  62. Pooley, C.; Furtado, K. Eliminating spurious velocities in the free-energy lattice Boltzmann method. Phys. Rev. E 2008, 77, 046702. [Google Scholar] [CrossRef]
  63. Huang, H.; Wang, L.; Lu, X.-y. Evaluation of three lattice Boltzmann models for multiphase flows in porous media. Comput. Math. Appl. 2011, 61, 3606–3617. [Google Scholar] [CrossRef]
  64. Tian, Z.; Xing, H.; Tan, Y.; Gao, J. A coupled lattice Boltzmann model for simulating reactive transport in CO2 injection. Phys. A: Stat. Mech. Its Appl. 2014, 403, 155–164. [Google Scholar] [CrossRef]
  65. Cameron, K.C.; Buchan, G.D. Darcys Law. In SW Encyclopedia of Water Science; Trimple: Westminster, CO, USA, 2008; pp. 143–146. [Google Scholar]
  66. Gasmi, H.; Khan, U.; Zaib, A.; Ishak, A.; Eldin, S.M.; Raizah, Z. Analysis of Mixed Convection on Two-Phase Nanofluid Flow Past a Vertical Plate in Brinkman-Extended Darcy Porous Medium with Nield Conditions. Mathematics 2022, 10, 3918. [Google Scholar] [CrossRef]
  67. Zhang, Q.; Choo, J.; Borja, R.I. On the preferential flow patterns induced by transverse isotropy and non-Darcy flow in double porosity media. Comput. Methods Appl. Mech. Eng. 2019, 353, 570–592. [Google Scholar] [CrossRef]
  68. Zhu, W.; Liu, Y.; Shi, Y.; Zou, G.; Zhang, Q.; Kong, D. Effect of dynamic threshold pressure gradient on production performance in water-bearing tight gas reservoir. Adv. Geo-Energy Res. 2022, 6, 286. [Google Scholar] [CrossRef]
  69. Lu, C.; Qin, X.; Ma, C.; Yu, L.; Geng, L.; Bian, H.; Zhang, K.; Zhou, Y. Investigation of the impact of threshold pressure gradient on gas production from hydrate deposits. Fuel 2022, 319, 123569. [Google Scholar] [CrossRef]
  70. Li, D.; Zha, W.; Liu, S.; Wang, L.; Lu, D. Pressure transient analysis of low permeability reservoir with pseudo threshold pressure gradient. J. Pet. Sci. Eng. 2016, 147, 308–316. [Google Scholar] [CrossRef]
  71. Hao, F.; Cheng, L.; Hassan, O.; Hou, J.; Liu, C.; Feng, J. Threshold pressure gradient in ultra-low permeability reservoirs. Pet. Sci. Technol. 2008, 26, 1024–1035. [Google Scholar] [CrossRef]
  72. Christopoulou, M.A.; Koutsovitis, P.; Kostoglou, N.; Paraskevopoulou, C.; Sideridis, A.; Petrounias, P.; Rogkala, A.; Stock, S.; Koukouzas, N. Evaluation of the CO2 Storage Capacity in Sandstone Formations from the Southeast Mesohellenic trough (Greece). Energies 2022, 15, 3491. [Google Scholar] [CrossRef]
  73. Ni, P.; Wang, S.; Wang, C.; Zhang, S. Estimation of REV size for fractured rock mass based on damage coefficient. Rock Mech. Rock Eng. 2017, 50, 555–570. [Google Scholar] [CrossRef]
  74. Min, K.-B.; Jing, L.; Stephansson, O. Determining the equivalent permeability tensor for fractured rock masses using a stochastic REV approach: Method and application to the field data from Sellafield, UK. Hydrogeol. J. 2004, 12, 497–510. [Google Scholar] [CrossRef]
  75. Boek, E.S.; Venturoli, M. Lattice-Boltzmann studies of fluid flow in porous media with realistic rock geometries. Comput. Math. Appl. 2010, 59, 2305–2314. [Google Scholar] [CrossRef]
  76. Xing, H.; Gao, J.; Zhang, J.; Liu, Y. Towards an integrated simulator for enhanced geothermal reservoirs. In Proceedings of the World Geothermal Congress 2010, Bali, Indonesia, 25-29 April 2010. [Google Scholar]
  77. Ziegler, D.P. Boundary conditions for lattice Boltzmann simulations. J. Stat. Phys. 1993, 71, 1171–1177. [Google Scholar] [CrossRef]
  78. Ju, Y.; Wang, J.; Gao, F.; Xie, H. Lattice-Boltzmann simulation of microscale CH4 flow in porous rock subject to force-induced deformation. Chin. Sci. Bull. 2014, 59, 3292–3303. [Google Scholar] [CrossRef]
  79. Zhao, W.; Jia, C.; Zhang, T.; Jiang, L.; Li, X.; Jiang, Z.; Zhang, F. Effects of nanopore geometry on confined water flow: A view of lattice Boltzmann simulation. Chem. Eng. Sci. 2021, 230, 116183. [Google Scholar] [CrossRef]
  80. Bhatnagar, P.L.; Gross, E.P.; Krook, M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 1954, 94, 511. [Google Scholar] [CrossRef]
  81. Gunstensen, A.K.; Rothman, D.H.; Zaleski, S.; Zanetti, G. Lattice Boltzmann model of immiscible fluids. Phys. Rev. A 1991, 43, 4320–4327. [Google Scholar] [CrossRef]
  82. Iglauer, S.; Pentland, C.; Busch, A. CO2 wettability of seal and reservoir rocks and the implications for carbon geo-sequestration. Water Resour. Res. 2015, 51, 729–774. [Google Scholar] [CrossRef]
  83. Vasheghani Farahani, M.; Hassanpouryouzband, A.; Yang, J.; Tohidi, B. Heat transfer in unfrozen and frozen porous media: Experimental measurement and pore-scale modeling. Water Resour. Res. 2020, 56, e2020WR027885. [Google Scholar] [CrossRef]
Figure 1. Schematic showing the CO2 unstable storage due to two-phase CO2–brine imbibition processes where the pore space saturated with CO2 (left) is flooded by brine, resulting in the displacement of in situ CO2 out of the pores (middle), until finally only the residual CO2 and brine are left in the pore space (right).
Figure 1. Schematic showing the CO2 unstable storage due to two-phase CO2–brine imbibition processes where the pore space saturated with CO2 (left) is flooded by brine, resulting in the displacement of in situ CO2 out of the pores (middle), until finally only the residual CO2 and brine are left in the pore space (right).
Energies 16 01547 g001
Figure 2. Example of gray-scale micro-CT scan of sandstone rock sample in this study.
Figure 2. Example of gray-scale micro-CT scan of sandstone rock sample in this study.
Energies 16 01547 g002
Figure 3. Workflow for imaging processing of raw CT image (left, e.g.,) through filtering, obtaining the noise-free image which is further threshold-segmented into separated pore space (blue) and rock matrix (gray) on the right.
Figure 3. Workflow for imaging processing of raw CT image (left, e.g.,) through filtering, obtaining the noise-free image which is further threshold-segmented into separated pore space (blue) and rock matrix (gray) on the right.
Energies 16 01547 g003
Figure 4. Schematic diagram showing REV central cutting.
Figure 4. Schematic diagram showing REV central cutting.
Energies 16 01547 g004
Figure 5. Schematic diagram of REVs oriented at three different subdomains (upper (orange), middle (gray) and bottom (green)) divided from a whole-core computational domain of CT-scanned rock sample.
Figure 5. Schematic diagram of REVs oriented at three different subdomains (upper (orange), middle (gray) and bottom (green)) divided from a whole-core computational domain of CT-scanned rock sample.
Energies 16 01547 g005
Figure 6. Results of REV analysis in terms of porosity spatial correlation (right), performed on three parts of the whole core domain defined in Figure 5, based on the REV central cutting method explained in Figure 5 with top view presented on the left.
Figure 6. Results of REV analysis in terms of porosity spatial correlation (right), performed on three parts of the whole core domain defined in Figure 5, based on the REV central cutting method explained in Figure 5 with top view presented on the left.
Energies 16 01547 g006
Figure 7. Schematic plot showing spatial position (left), gray-scale CT image (middle) and extracted pore space (right) of REVs characterizing upper, middle, and bottom subdomains of the whole-core domain.
Figure 7. Schematic plot showing spatial position (left), gray-scale CT image (middle) and extracted pore space (right) of REVs characterizing upper, middle, and bottom subdomains of the whole-core domain.
Energies 16 01547 g007
Figure 8. Schematic diagram of D3Q19 model.
Figure 8. Schematic diagram of D3Q19 model.
Energies 16 01547 g008
Figure 9. Schematic diagram showing the contact angle where a fluid droplet drops onto the rock surface [5].
Figure 9. Schematic diagram showing the contact angle where a fluid droplet drops onto the rock surface [5].
Energies 16 01547 g009
Figure 10. Cross section of the pressure field of REV1 with a square bounded by black dash lines, indicating the scope of a tighter 200 × 200 × 200 voxels REV to be extracted (a) and the pore structure(b) of it denoted in green bounded by a cube as the computational domain.
Figure 10. Cross section of the pressure field of REV1 with a square bounded by black dash lines, indicating the scope of a tighter 200 × 200 × 200 voxels REV to be extracted (a) and the pore structure(b) of it denoted in green bounded by a cube as the computational domain.
Energies 16 01547 g010
Figure 11. Pressure (ac) and velocity (df) distributions of REV1, REV2 and REV3, two parameters are presented in the dimensionless form ranging from 0 to 1, aiming to illustrate the big contrast.
Figure 11. Pressure (ac) and velocity (df) distributions of REV1, REV2 and REV3, two parameters are presented in the dimensionless form ranging from 0 to 1, aiming to illustrate the big contrast.
Energies 16 01547 g011
Figure 12. Velocity distributions of two-phase flow (nonwetting phase–CO2, wetting phase–brine) in the tighter REV when Pdiff is equal to 100 Pa (a),200 Pa (b),300 Pa (c),400 Pa (d) at SCO2 = 0.5, Ca = 1 × 10−4, M = 10 and WI = 0.7.
Figure 12. Velocity distributions of two-phase flow (nonwetting phase–CO2, wetting phase–brine) in the tighter REV when Pdiff is equal to 100 Pa (a),200 Pa (b),300 Pa (c),400 Pa (d) at SCO2 = 0.5, Ca = 1 × 10−4, M = 10 and WI = 0.7.
Energies 16 01547 g012
Figure 13. Pressure profiles of the tighter REV versus CT slice when steady state is reached as Pdiff is equal to 100 Pa (a), 200 Pa (b), 300 Pa (c), 400 Pa (d) with the same initial conditions as Figure 12.
Figure 13. Pressure profiles of the tighter REV versus CT slice when steady state is reached as Pdiff is equal to 100 Pa (a), 200 Pa (b), 300 Pa (c), 400 Pa (d) with the same initial conditions as Figure 12.
Energies 16 01547 g013
Figure 14. The tighter REV saturated with CO2 (a) was initially saturated with brine in the near-inlet region with global target brine saturation of 0.05, (b) as prepared for the CO2–brine imbibition process over a range of target saturation from 0.2 to 0.9 with initial model conditions the same as in Section 4.2.
Figure 14. The tighter REV saturated with CO2 (a) was initially saturated with brine in the near-inlet region with global target brine saturation of 0.05, (b) as prepared for the CO2–brine imbibition process over a range of target saturation from 0.2 to 0.9 with initial model conditions the same as in Section 4.2.
Energies 16 01547 g014
Figure 15. Steady-state brine distribution of the tighter REV, during the CO2–brine imbibition process over a range of target saturation from 0.2 to 0.5., the flow pattern and fluid distribution are stabilized when the global target saturation reaches 0.5 where the flow is channelized and dominated by brine.
Figure 15. Steady-state brine distribution of the tighter REV, during the CO2–brine imbibition process over a range of target saturation from 0.2 to 0.5., the flow pattern and fluid distribution are stabilized when the global target saturation reaches 0.5 where the flow is channelized and dominated by brine.
Energies 16 01547 g015
Figure 16. Simulated relative permeability curves for CO2 (nonwetting) under the total 18 scenarios, with different combinations of three nondimensional parameters including capillary number (Ca), wettability index (WI) and viscosity ratio (M) described in Table 1. WI = 0.3, WI = 0.7 and WI = 0.9 respectively weakly water-wet, intermediately water-wet, and strongly water-wet. (Note: (a) and (b) represent cases under high capillary number (Ca = 1 × 10−3), while those under low capillary numbers (Ca = 1 × 10−4, Ca = 1 × 10−5) are done by (cf).)
Figure 16. Simulated relative permeability curves for CO2 (nonwetting) under the total 18 scenarios, with different combinations of three nondimensional parameters including capillary number (Ca), wettability index (WI) and viscosity ratio (M) described in Table 1. WI = 0.3, WI = 0.7 and WI = 0.9 respectively weakly water-wet, intermediately water-wet, and strongly water-wet. (Note: (a) and (b) represent cases under high capillary number (Ca = 1 × 10−3), while those under low capillary numbers (Ca = 1 × 10−4, Ca = 1 × 10−5) are done by (cf).)
Energies 16 01547 g016
Figure 17. A close-up view of CO2 relative permeability curves under 18 scenarios at Table 1 by narrowing the scope of brine saturation within which most of the values have fallen into the scale of 0.1. A lower bound value of krCO2 = 0.05, as the minimum value for assuming the flow of CO2 through the computational domain is present, is defined; therefore, brine saturations at which the CO2 relative permeability firstly falls below the 0.05 are considered approximate maximum brine saturation (SBrinemax) in the computational domain where CO2 residual saturation (SresidualCO2) is consequently estimated using Equation (26). (Note: (a) and (b) represent cases under high capillary number (Ca = 1 × 10−3), while those under low capillary numbers (Ca = 1 × 10−4, Ca = 1 × 10−5) are done by (cf).)
Figure 17. A close-up view of CO2 relative permeability curves under 18 scenarios at Table 1 by narrowing the scope of brine saturation within which most of the values have fallen into the scale of 0.1. A lower bound value of krCO2 = 0.05, as the minimum value for assuming the flow of CO2 through the computational domain is present, is defined; therefore, brine saturations at which the CO2 relative permeability firstly falls below the 0.05 are considered approximate maximum brine saturation (SBrinemax) in the computational domain where CO2 residual saturation (SresidualCO2) is consequently estimated using Equation (26). (Note: (a) and (b) represent cases under high capillary number (Ca = 1 × 10−3), while those under low capillary numbers (Ca = 1 × 10−4, Ca = 1 × 10−5) are done by (cf).)
Energies 16 01547 g017
Table 1. All 18 scenarios constituting different combinations of three nondimensional parameters (capillary number (Ca), viscosity ratio (M) and wettability (WI)) for two-phase CO2–brine imbibition processes.
Table 1. All 18 scenarios constituting different combinations of three nondimensional parameters (capillary number (Ca), viscosity ratio (M) and wettability (WI)) for two-phase CO2–brine imbibition processes.
Scenario No.Capillary Number (Ca)Viscosity Ratio (M)Wettability (WI)
11 × 10−310.3
21 × 10−410.3
31 × 10−510.3
41 × 10−310.7
51 × 10−410.7
61 × 10−510.7
71 × 10−310.9
81 × 10−410.9
91 × 10−510.9
101 × 10−3100.3
111 × 10−4100.3
121 × 10−5100.3
131 × 10−3100.7
141 × 10−4100.7
151 × 10−5100.7
161 × 10−3100.9
171 × 10−4100.9
181 × 10−5100.9
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Wan, Y.; Jia, C.; Zhao, W.; Jiang, L.; Chen, Z. Micro-Scale Lattice Boltzmann Simulation of Two-Phase CO2–Brine Flow in a Tighter REV Extracted from a Permeable Sandstone Core: Implications for CO2 Storage Efficiency. Energies 2023, 16, 1547. https://doi.org/10.3390/en16031547

AMA Style

Wan Y, Jia C, Zhao W, Jiang L, Chen Z. Micro-Scale Lattice Boltzmann Simulation of Two-Phase CO2–Brine Flow in a Tighter REV Extracted from a Permeable Sandstone Core: Implications for CO2 Storage Efficiency. Energies. 2023; 16(3):1547. https://doi.org/10.3390/en16031547

Chicago/Turabian Style

Wan, Yidi, Chengzao Jia, Wen Zhao, Lin Jiang, and Zhuxin Chen. 2023. "Micro-Scale Lattice Boltzmann Simulation of Two-Phase CO2–Brine Flow in a Tighter REV Extracted from a Permeable Sandstone Core: Implications for CO2 Storage Efficiency" Energies 16, no. 3: 1547. https://doi.org/10.3390/en16031547

APA Style

Wan, Y., Jia, C., Zhao, W., Jiang, L., & Chen, Z. (2023). Micro-Scale Lattice Boltzmann Simulation of Two-Phase CO2–Brine Flow in a Tighter REV Extracted from a Permeable Sandstone Core: Implications for CO2 Storage Efficiency. Energies, 16(3), 1547. https://doi.org/10.3390/en16031547

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop