A Case Study on the CO₂ Sequestration in Shenhua Block Reservoir: The Impacts of Injection Rates and Modes

Abstract

Carbon capture and storage (CCS) is the most promising method of curbing atmospheric carbon dioxide levels from 2020 to 2050. Accurate predictions of geology and sealing capabilities play a key role in the safe execution of CCS projects. However, popular forecasting methods often oversimplify the process and fail to guide actual CCS projects in the right direction. This study takes a specific block in Shenhua, China as an example. The relative permeability of CO₂ and brine is measured experimentally, and a multi-field coupling CO₂ storage prediction model is constructed, focusing on analyzing the sealing ability of the block from the perspective of injection modes. The results show that when injected at a constant speed, the average formation pressure and wellbore pressure are positively correlated with the CO₂ injection rate and time; when the injection rate is 0.5 kg/s for 50 years, the average formation pressure increases by 38% and the wellbore pressure increases by 68%. For different injection modes, the average formation pressures of various injection methods are similar during injection. Among them, the pressure increases around the well in the decreasing injection mode is the smallest. The CO₂ concentration around the wellbore is the largest, and the CO₂ diffusion range continues to expand with injection time. In summary, formation pressure increases with the increase in injection rate and injection time, and the decreasing injection mode has the least impact on the increase in formation pressure. The CO₂ concentration is the largest around the well, and the CO₂ concentration gradually decreases. The conclusion helps determine the geological carrying capacity of injection volumes and provides insights into the selection of more appropriate injection modes. Accurate predictions of CO₂ storage capacity are critical to ensuring project safety and monitoring potentially hazardous sites based on reservoir characteristics.

1. Introduction

CCS aims to mitigate human-induced carbon dioxide emissions by injecting and storing carbon dioxide in specific geological structures [1,2]. In the pursuit of achieving carbon neutrality by the mid-21st century, CCS stands out as a pivotal carbon-negative technology, garnering significant attention and interest from countries globally [3,4]. Before CO₂ injection can proceed, a proper assessment of the risk of CO₂ leakage from injection wells and geological storage sites must be conducted [1,5]. The Shenhua Carbon Capture and Storage Demonstration Project in China’s Ordos Basin stands as Asia’s first and largest full-chain saline aquifer carbon dioxide storage project. There is a lot of engineering and research going on there. These studies include stress and deformation changes induced by injection, potential damage modes and safety factors, interactions between coal mining and carbon dioxide geology storage, and determination of injection pressure limits, and the upper limit of wellhead pressure is 18 MPa, which is reliable [6,7]. Prior to project implementation, a rigorous consideration of the impact of fluid flow in the formation and an accurate prediction of the formation’s storage capacity are essential to ensure the safe development of the project [8,9]. CCS is a complex process that requires focus on its impact on formation pressure and CO₂ distribution. Different injection modes have important effects on the formation.

Geological formations such as basalts, coal seams, depleted oil reserves, soils, deep saline aquifers, and sedimentary basins exhibit vast potential for carbon dioxide storage [10]. There is potential for carbon dioxide (CO₂) recovery in ultra-deep water subsalt carbonate reservoirs for carbon capture and storage [11]. To facilitate the geological storage of CO₂, the pressure must exceed 7.38 MPa and the temperature must surpass 31.1 °C, indicating a theoretical storage depth exceeding 800 m. The Shenhua saline aquifer block in China satisfies these requirements. The profound geological structure’s complexity, exploration extent, and limitations in indoor physical simulation experimental conditions pose challenges to reservoir characterization and geological modeling. Simultaneously, accurately predicting CO₂ migration patterns and ensuring reservoir safety are critical issues in current research. However, there is no very definite conclusion about the storage capacity of the reservoir. In particular, the research on the storage capacity under different injection modes is still blank.

Countries such as the United States, China, Russia, the United Kingdom, Croatia, and India are actively accelerating the global deployment of CCS, making significant contributions to the reduction in global greenhouse gas emissions [12,13,14,15,16,17,18]. Challenges such as carbon dioxide leakage, energy inefficiency, and high implementation costs pose significant obstacles to the development of CCS. The safety assessment of CCS represents one of the greatest challenges; accurate predictions of geological carrying capacity and storage capacity are essential prerequisites for formulating and implementing a viable plan [19,20,21,22,23,24].

Berrezueta conducted laboratory studies on carbon dioxide–brine–rock interactions and performed some sensitivity analyses [25]. Xie investigated the influence of geological and engineering parameters on CO₂ migration and flow characteristics through indoor injection experiments, supplemented by X-ray computed tomography (CT) and scanning electron microscopy (SEM) experiments and computational fluid dynamics (CFD) numerical simulations [26]. These experimental studies did not conduct real simulations under formation conditions and could not directly guide the project. Hu integrated CFD simulation technology into the experimental study of atmospheric CO₂ diffusion in full-scale blasting emission tests of high-pressure supercritical phase CO₂ pipelines, and quantitatively analyzed the relationship between supercritical CO₂ leakage diameter and dangerous distance [27]. However, this study did not directly analyze the formation CO₂ distribution. Tutolo used high-performance computing techniques to study the coupled effects of cold CO₂ injection and background hydraulic head gradients on reservoir-scale mineral volume changes. Research has found that the migration and flow characteristics of CO₂ in sandstone during the geological storage process have a significant impact on the physical and mechanical properties of the rock [28]. Therefore, research on CO₂ storage must consider the physical and mechanical properties of rocks and fluid flow characteristics. Yang utilized the VOF (Volume of Fluid) method, capable of tracking dynamic changes in the two-phase interface, to establish two-dimensional and three-dimensional models and numerically simulate a supercritical CO₂–brine two-phase flow [29]. Without considering the influence of reservoir mechanical properties on seepage characteristics, it is inaccurate to simply study two-phase flow. In water–mechanical–chemical-coupled simulations, simplified flow mechanisms can lead to significant deviations in predicted throughput and storage performance [30]. Ratnakar and Omosebi et al. developed a machine learning-based workflow to inject single-phase supercritical carbon dioxide into deep saline aquifers to assess leakage risks [31,32,33,34,35,36]. The shortcoming is that these studies did not conduct sufficient and effective analysis and research on formation pressure changes.

Given the limitations of the current body of research, this study addresses the relatively singular factors considered and explores other issues. It entails experimental measurements of CO₂ migration and brine flow characteristics under different driving pressures. Additionally, a comprehensive large-scale multi-field coupling model of the reservoir (encompassing seepage, chemical diffusion, and solid mechanics fields) was established. Experimental data were incorporated into the model, and the influence of formation pressure on rock permeability characteristics was thoroughly examined. The primary focus of this study is the innovative exploration of various CO₂ injection modes, with an evaluation of reservoir storage capacity and risk conducted through the analysis of changes in pressure around the well, average formation pressure alterations, and CO₂ distribution. Taking a specific CO₂ geological storage project in the Shenhua saline aquifer as the research subject, the study integrates experiments and simulations, aligning with the actual engineering background and conditions. This approach aims to elucidate the CO₂ geological storage mechanism in saline aquifers. The research methods and conclusions derived from this study provide valuable insights into the geological storage mechanism and seepage laws of CO₂ in saline aquifers, playing a pivotal role in informing the scientific and safe implementation of storage projects.

The framework of this article is structured as follows: Section 2 delves into core methods, including mathematical models, physical models, physical properties, seepage characteristic parameters, introduction to injection methods, and an overview of the block. Section 3 engages in a discussion of the results, covering model verification, reservoir pressure comparison, and CO₂ distribution. Finally, the article concludes with a summary.

2. Methodology

This study fully considered the physical properties of the reservoir and fluid. Multi-fields mainly include multiphase transfer field, Darcy seepage field, and solid mechanics field. The coupling method is introduced in detail in the mathematical model. Figure 1 is the flow of the process.

2.1. Mathematical Model

This study takes into account multi-field coupling (encompassing seepage, chemical diffusion, and solid mechanics fields), diffusion effects, and effective stress. The mathematical model comprises the multiphase fluid flow mass conservation theory, seepage mechanics momentum equation, solid mechanics stress balance differential equation, constitutive equation, geometric equation, Terzaghi effective stress principle, and diffusion equation. This section primarily introduces each equation and its physical meaning, elucidating how they are coupled to establish connections.

The continuity equation of multiphase seepage delineates the mass conservation of multiphase mixed fluids. This equation is articulated in terms of the volume fraction of each phase [37]:

$$\frac{\partial {ϵ}_p{\rho }_{s_i}{s}_i}{\partial t}+\nabla ·N_i=0$$

(1)

Among them, ${ϵ}_p$ is the porosity, ${\rho }_{s_i}$ is the fluid density, and $s_i$ is the volume fraction.

$$N_i={\rho }_{Si}{u}_i$$

(2)

$u_i$ is the fluid velocity.

Darcy’s seepage flux equation is a constitutive equation that characterizes the flow of liquid through porous media. This equation finds extensive application in petroleum engineering and groundwater engineering:

$$u_i=-\frac{{\kappa }_{{rs}_i}}{{\mu }_{s_i}}\kappa (\nabla p-{\rho }_{i}g)$$

(3)

i = 1, 2, 3. $s_1$ is brine saturation. $s_2$ is carbon dioxide saturation. $s_3$ is bound brine saturation.

$$s_1+s_2+s_3=1$$

(4)

If Darcy’s multiphase seepage equation considers diffusion effects, then:

$$u_i=-\frac{{\kappa }_{{rs}_i}}{{\mu }_{s_i}}\kappa (\nabla p-{\rho }_{i}g)-D_{{cs}_i}\nabla s_i$$

(5)

$D_{{cs}_i}$ is the diffusion coefficient, and the value here is 6 × 10⁻⁹ m²/s from the literature [38].

The stress balance equation of a solid elucidates the equilibrium of forces at each point within a stationary solid. In a three-dimensional space, for a solid within a volume element and considering three directions (x, y, and z), the stress balance equation can be expressed as [39]:

$$\left.\begin{array}{c}\frac{\partial {\sigma }_x}{\partial x}+\frac{\partial {\tau }_{yx}}{\partial y}+\frac{\partial {\tau }_{zx}}{\partial z}+f_x=0\\ \frac{\partial {\sigma }_y}{\partial y}+\frac{\partial {\tau }_{zy}}{\partial z}+\frac{\partial {\tau }_{xy}}{\partial x}+f_y=0\\ \frac{\partial {\sigma }_z}{\partial z}+\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\tau }_{yz}}{\partial y}+f_z=0\end{array}\right\}$$

(6)

${\sigma }_x$, ${\sigma }_y$, ${\sigma }_z$, ${\tau }_{xy}$, ${\tau }_{yz}$, ${\tau }_{zx}$ are the stress components, and $f_x$, $f_y$, $f_z$ are the body force components. Here, we only consider gravity:

$$\left.\begin{array}{c}{f}_x=0 \\ f_y=0 \\ f_z=\rho g\end{array}\right\}$$

(7)

The relationship between the shear stress components and displacement components:

$$\left.\begin{array}{c}{\tau }_{yx}={\tau }_{xy}=\frac{E}{2\left(1+\nu \right)}\left(\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}\right)\\ {\tau }_{xz}={\tau }_{zx}=\frac{E}{2\left(1+\nu \right)}\left(\frac{\partial u_z}{\partial x}+\frac{\partial u_x}{\partial z}\right)\\ {\tau }_{zy}={\tau }_{yz}=\frac{E}{2\left(1+\nu \right)}\left(\frac{\partial u_y}{\partial z}+\frac{\partial u_z}{\partial y}\right)\end{array}\right\}$$

(8)

$E$ is the elastic modulus; $\nu$ is the Poisson’s ratio, and $u_x$, $u_y$, $u_z$ are the displacement components.

Terzaghi’s effective stress principle asserts that while the stress in the soil is borne by both the soil skeleton and the water vapor in the soil, only the effective stress transmitted through the soil particles induces soil deformation. The pressure transmitted through water vapor in the pores does not contribute to the strength of the soil [39]:

$$\left.\begin{array}{c}{\sigma }_x={\sigma }_x^{\prime }+\alpha p\\ {\sigma }_y={\sigma }_y^{\prime }+\alpha p\\ {\sigma }_z={\sigma }_z^{\prime }+\alpha p\end{array}\right\}$$

(9)

${\sigma }_x^{\prime }$, ${\sigma }_y^{\prime }$, ${\sigma }_z^{\prime }$ are the effective stress components; $\alpha$ is the Biot coefficient, and $p$ is the reservoir pressure. Equilibrium Equation (10) by bagging (7), (8), and (9) into (6):

$$\left.\begin{array}{c}\frac{\partial {\sigma }_x^{\prime }}{\partial x}+\alpha \frac{\partial p}{\partial x}+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_x}{\partial y^2}+\frac{{\partial }^2{u}_y}{\partial x\partial y}\right)+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_x}{\partial z^2}+\frac{{\partial }^2{u}_z}{\partial x\partial z}\right)=0\\ \frac{\partial {\sigma }_y^{\prime }}{\partial y}+\alpha \frac{\partial p}{\partial y}+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_y}{\partial z^2}+\frac{{\partial }^2{u}_z}{\partial y\partial z}\right)+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_y}{\partial x^2}+\frac{{\partial }^2{u}_x}{\partial y\partial x}\right)=0\\ \frac{\partial {\sigma }_z^{\prime }}{\partial z}+\alpha \frac{\partial p}{\partial z}+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_z}{\partial x^2}+\frac{{\partial }^2{u}_x}{\partial z\partial x}\right)+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_z}{\partial y^2}+\frac{{\partial }^2{u}_y}{\partial z\partial y}\right)+\rho g=0\end{array}\right\}$$

(10)

And three-dimensional partial differential equilibrium equation, denoted as Equation (11):

$$\left.\begin{array}{c}\frac{E\upsilon }{\left(1+\upsilon \right)\left(1-2\upsilon \right)}\left(\frac{{\partial }^2{u}_x}{\partial x^2}+\frac{{\partial }^2{u}_y}{\partial x\partial y}+\frac{{\partial }^2{u}_z}{\partial x\partial z}\right)+\frac{E}{\left(1+\upsilon \right)}\frac{{\partial }^2{u}_x}{\partial x^2}+\alpha \frac{\partial p}{\partial x}+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_x}{\partial y^2}+\frac{{\partial }^2{u}_y}{\partial x\partial y}\right)+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_x}{\partial z^2}+\frac{{\partial }^2{u}_z}{\partial x\partial z}\right)=0\\ \frac{E\upsilon }{\left(1+\upsilon \right)\left(1-2\upsilon \right)}\left(\frac{{\partial }^2{u}_x}{\partial y\partial x}+\frac{{\partial }^2{u}_y}{\partial y^2}+\frac{{\partial }^2{u}_z}{\partial y\partial z}\right)+\frac{E}{\left(1+\upsilon \right)}\frac{{\partial }^2{u}_y}{\partial y^2}+\alpha \frac{\partial p}{\partial y}+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_y}{\partial z^2}+\frac{{\partial }^2{u}_z}{\partial y\partial z}\right)+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_y}{\partial x^2}+\frac{{\partial }^2{u}_x}{\partial y\partial x}\right)=0\\ \frac{E\upsilon }{\left(1+\upsilon \right)\left(1-2\upsilon \right)}\left(\frac{{\partial }^2{u}_x}{\partial z\partial x}+\frac{{\partial }^2{u}_y}{\partial z\partial y}+\frac{{\partial }^2{u}_z}{\partial z^2}\right)+\frac{E}{\left(1+\upsilon \right)}\frac{{\partial }^2{u}_z}{\partial z^2}+\alpha \frac{\partial p}{\partial z}+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_z}{\partial x^2}+\frac{{\partial }^2{u}_x}{\partial z\partial x}\right)+\frac{E}{2\left(1+\upsilon \right)}\left(\frac{{\partial }^2{u}_z}{\partial y^2}+\frac{{\partial }^2{u}_y}{\partial z\partial y}\right)+\rho g=0\end{array}\right\}$$

(11)

We bring ${\nabla }^2=\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}$ into (11) and simplify the solid (stress field) Equation (12) for carbon dioxide reservoir calculation considering the effect of effective stress:

$$\left.\begin{array}{c}\frac{E}{2\left(1+\upsilon \right)\left(1-2\upsilon \right)}\frac{\partial {\epsilon }_V}{\partial x}+\frac{E}{2\left(1+\upsilon \right)}{\nabla }^2{u}_x+\alpha \frac{\partial p}{\partial x}=0\\ \frac{E}{2\left(1+\upsilon \right)\left(1-2\upsilon \right)}\frac{\partial {\epsilon }_V}{\partial y}+\frac{E}{2\left(1+\upsilon \right)}{\nabla }^2{u}_y+\alpha \frac{\partial p}{\partial y}=0\\ \frac{E}{2\left(1+\upsilon \right)\left(1-2\upsilon \right)}\frac{\partial {\epsilon }_V}{\partial z}+\frac{E}{2\left(1+\upsilon \right)}{\nabla }^2{u}_z+\alpha \frac{\partial p}{\partial z}+\rho g=0\end{array}\right\}$$

(12)

${\epsilon }_x$, ${\epsilon }_y$ and ${\epsilon }_z$ are the strain components, and ${\epsilon }_V$ is the volume component:

$$\left.\begin{array}{c}{\epsilon }_x=\frac{\partial u_x}{\partial x}\\ {\epsilon }_y=\frac{\partial u_y}{\partial y}\\ {\epsilon }_z=\frac{\partial u_z}{\partial z}\\ {\epsilon }_V={\epsilon }_x+{\epsilon }_y+{\epsilon }_z\end{array}\right\}$$

(13)

α takes 1, which is experience from the engineering site, then:

$$\left.\begin{array}{c}\frac{E}{2\left(1+\upsilon \right)\left(1-2\upsilon \right)}\frac{\partial {\epsilon }_V}{\partial x}+\frac{E}{2\left(1+\upsilon \right)}{\nabla }^2{u}_x+\frac{\partial p}{\partial x}=0\\ \frac{E}{2\left(1+\upsilon \right)\left(1-2\upsilon \right)}\frac{\partial {\epsilon }_V}{\partial y}+\frac{E}{2\left(1+\upsilon \right)}{\nabla }^2{u}_y+\frac{\partial p}{\partial y}=0\\ \frac{E}{2\left(1+\upsilon \right)\left(1-2\upsilon \right)}\frac{\partial {\epsilon }_V}{\partial z}+\frac{E}{2\left(1+\upsilon \right)}{\nabla }^2{u}_z+\frac{\partial p}{\partial z}+\rho g=0\end{array}\right\}$$

(14)

Subsequently, the connection between the multiphase seepage field and the solid mechanical field can be established through Equations (5) and (14), considering both the effective stress principle and diffusion principle.

2.2. Physical Models and Numerical Methods

This multi-field coupling model employs the finite element method (FEM, COMSOL) to numerically solve the aforementioned equations. As illustrated in Figure 2, the geometric model dimensions are 2400 m × 2400 m × 300 m. A 1/4 symmetrical structure was utilized for the study. Our research hypothesis involves incorporating roller supports around and at the bottom of the reservoir to restrict normal movement. The upper surface of the reservoir is free, and the outer boundary serves for outflow. No flux is present at the boundary, the internal boundary of the reservoir, and the upper and lower boundaries. FEM calculations were executed based on the mathematical and physical models, with grid divisions as depicted in Figure 3. The total number of units is 242,895, and the grid around the well is dense. The maximum element size is 31.8 m; the minimum element size is 6 m; the maximum element growth rate is 1.13; the curvature factor is 0.5; the minimum element quality is 0.1842; the average element quality is 0.6632; and the overall quality of the grid is good.

According to engineering site data, the Young’s modulus of the entire reservoir is 1.3 GPa, the Poisson’s ratio is 0.23, and the density is 2560 kg/m³. The salinity of brine is 23.5 g/L with a density of 0.984 g·cm⁻³. The main ionic components are Na⁺ and Cl⁻. The relationship between porosity and permeability and reservoir pressure is as follows:

$$\left.\begin{array}{c}epsilon=0.2984\times (1+(2\times 10^\left(-2.5\right))\times p/1[MPa\left]\right)\\ kappa=1\times 10^\left(-12\right)\times {\left(\frac{epsilon}{0.2984}\right)}^3(m^2)\end{array}\right\}$$

(15)

Epsilon is porosity, and kappa is permeability. And the empirical formula is given by the project.

2.3. Physical Properties and Seepage Characteristic Parameters

Figure 4 illustrates the experimental system designed for testing the relative permeability of supercritical CO₂–brine. The system comprises a CO₂ storage tank, booster pump, pressure gauge, incubator, core holder, confining pressure pump, back pressure valve, back pressure pump, gas–liquid separator, and gas–water metering device (Appendix A shows the experimental equipment and procedure). The experiment simulates a temperature of 70 °C (temperature of the reservoir), with inlet and outlet pressures set at 10 MPa and 8 MPa, respectively. Under these conditions, CO₂ attains a supercritical state. Figure 5 illustrates the result of relative permeability.

The relative permeability of brine and CO₂ is obtained:

$$\left.\begin{array}{c}{K}_{rw}=\left(\frac{\left(s_w-0.1706\right)}{0.8294}\right)^7.661\\ K_{rg}=0.1909\times \left(1-\frac{\left(s_w-0.1706\right)}{0.8294}\right)^3.502\end{array}\right\}$$

(16)

The density and viscosity change curves of supercritical carbon dioxide with pressure at 70 °C are plotted based on the thermophysical parameters from the National Institute of Standards and Technology (NIST) in Figure 6 and Figure 7.

A magnetic levitation balance with precision 1 $\mathsf{\mu }\mathrm{g}$ was employed for measuring the density of saltwater in a measuring cylinder. This electromagnetic levitation balance utilizes a combination of an electromagnet positioned outside the measuring container and a permanent magnet inside the measuring container to directly measure the absolute density of the fluid within the isolated and closed measuring container.

In the investigation of the impact of injection pressure on the density of saltwater under CO₂ storage conditions within the saltwater layer, the baseline temperature was set to 70 °C. Each set of experiments was conducted over a time period of 120 h, and the density under various pressures is illustrated in Figure 8.

Under identical conditions, the viscosity of the carbon dioxide aqueous solution, as measured with a viscometer, is presented in Figure 9 [40].

2.4. Introduction to Injection Methods

This study delves into the carbon sequestration capacity of the reservoir from the perspectives of injection amount and injection mode.

Table 1. Injection modes.

	0–10 Years	10–20 Years	20–25 Years	25–30 Years	30–40 Years	40–50 Years
Mode 1	0.1	0.1	0.1	0.1	0.1	0.1
Mode 2	0.2	0.2	0.2	0.2	0.2	0.2
Mode 3	0.3	0.3	0.3	0.3	0.3	0.3
Mode 4	0.4	0.4	0.4	0.4	0.4	0.4
Mode 5	0.5	0.5	0.5	0.5	0.5	0.5
Mode 6	0.5	0.4	0.3	0.3	0.2	0.1
Mode 7	0.1	0.2	0.3	0.3	0.4	0.5
Mode 8	0.5	0.5	0.5	0.1	0.1	0.1
Mode 9	0.4	0.4	0.4	0.2	0.2	0.2
Mode 10	0.1	0.1	0.1	0.5	0.5	0.5
Mode 11	0.2	0.2	0.2	0.4	0.4	0.4

shows all injection modes. Five injection rates modes: 0.1, 0.2, 0.3, 0.4, and 0.5 kg/s. Seven injection modes: 10 years at intervals of 0.5–0.1 kg/s; 10 years at intervals of 0.1–0.5 kg/s; 25 years at 0.5 kg/s, 25 years at 0.1 kg/s; 25 years at 0.4 kg/s and 25 years at 0.2 kg/s; 0.1 kg/s in 25 years, 0.5 kg/s in 25 years; 0.2 kg/s in 25 years and 0.4 kg/s in 25 years, and Mode 3. We compared the changes in pressure over time: pressure around the well, average formation pressure, and changes in CO₂ distribution over time.

2.5. Block Introduction

The Shenhua Group is currently executing China’s inaugural full-chain carbon dioxide capture and geological storage demonstration project, situated in the Ordos Basin in the eastern part of northwest China. As depicted in Figure 10, it spans five provinces (autonomous regions), including Shaanxi, Shanxi, and Inner Mongolia, covering a total area of more than 27.6 × 10⁴ km². The Ordos Basin can be divided into six primary tectonic unit structures based on the history of geological structural changes. These include the Yishaan slope, the western margin thrust belt, the Shanxi burned skirt belt, the Tianhuan depression, the Yimeng uplift in the north, and the Yimeng uplift in the south. The Ordos Basin stands as one of the largest terrestrial sedimentary basins in China, characterized as a craton sedimentary basin. It lacks major fault zones traversing the entire basin, exhibits geological stability, even stress distribution, and boasts a thick sedimentary layer (with an average thickness of about 6000 m). Given these geological characteristics, it can conservatively be inferred that the Ordos formation possesses a significant geological storage capacity for CO₂ [39].

In this project, carbon dioxide from coal tail gases is captured through liquid processing and stored in deep brine aquifers. The primary target layer for CO₂ injection is the saline aquifer beneath the mined coal seam. The formation receiving the carbon dioxide injection is characterized by low porosity, low permeability, and high heterogeneity. The project is currently operating successfully, with no reported CO₂ leaks or associated environmental hazards, and only minor pressure build-up has been observed.

3. Results and Discussion

3.1. Model Validation

To validate the accuracy of the multi-field coupling simulation, we conducted permeability experiments and simulations on cores, comparing them based on the relationship between permeability and reservoir pressure at the engineering site. Water injection experiments were carried out using on-site provided cores, with a core length of 6 cm and a diameter of 2.5 cm. The outlet pressure was set to 8 MPa, and the inlet pressure varied at 10 MPa, 12 MPa, 14 MPa, 16 MPa, 18 MPa, and 20 MPa.

Figure 11 illustrates the relationship between the average core flow velocity and pressure for three experiment cases: one without considering solid mechanics, one considering solid mechanics in simulation, and one with experiments. We observed that simulations considering solid mechanics align closely with experimental results. However, simulations neglecting solid mechanics introduce increasing errors as the pressure rises. Therefore, to ensure simulation accuracy, accounting for the influence of rock mass solid mechanics is essential. According to Equation (15), as the pressure increases, the permeability increases.

3.2. Reservoir Pressure Comparison

In the simulations, the initial formation pressure is 8 MPa, and the temperature is 70 °C, placing CO₂ in a supercritical state. Injection was conducted over 50 years at five rates (0.1–0.5 kg/s) with intervals of 0.1 kg/s. Figure 12a,b illustrate the changes in pressure around the well and in the formation at different injection rates. It is evident that both the pressure around the well and the average pressure in the formation increase with CO₂ injection. The higher the injection amount per unit time, the greater the pressure change.

Figure 12c,d depict pressure cloud diagrams with injection rates of 0.1 kg/s and 0.5 kg/s, respectively. The pressure in the reservoir rises annually with injection, and the pressure around the well is notably higher than in other locations. Due to the influence of gravity, the pressure value in the lower layer of the reservoir is higher. Key findings for different injection rates after 50 years include:

• 0.1 k g/s: Average max formation pressure is about 8.6 MPa (7% higher), and max wellbore pressure is about 9.1 MPa (14% higher, the average pressure value of the wellbore).
• 0.2 kg/s: Average max formation pressure is about 9.2 MPa (15% higher), and max wellbore pressure is about 10.4 MPa (30% higher).
• 0.3 kg/s: Average max formation pressure is about 9.8 MPa (23% higher), and max wellbore pressure is about 11.4 MPa (42% higher).
• 0.4 kg/s: Average max formation pressure is about 10.6 MPa (32% higher), and max wellbore pressure is about 12.6 MPa (57% higher).
• 0.5 kg/s: Average max formation pressure is about 11.0 MPa (38% higher), and max wellbore pressure is about 13.5 MPa (68% higher).

Therefore, in CCS projects, considering the specific working conditions of the reservoir is crucial to estimating the maximum injection rate. This consideration becomes particularly important for controlling the injection rate and determining the appropriate injection time.

Similarly, for the other six injection modes mentioned in Section 2.4, the average injection rate is 0.3 kg/s: 10 years is an interval of 0.5–0.1 kg/s; 10 years is an interval of 0.1–0.5 kg/s; 25 years is 0.5 kg/s, 0.1 kg/s in 25 years; 0.4 kg/s in 25 years and 0.2 kg/s in 25 years; 0.1 kg/s in 25 years, 0.5 kg/s in 25 years; 0.2 kg/s in 25 years, 0.4 kg/s in 25 years. Analysis of the results. Examining Figure 13, which shows the average pressure changes around the well and in the formation with different injection modes, several observations can be made:

• Average Formation Pressure: It increases with injection time, peaking at 50 years for various injection modes. The maximum values at 50 years are 9.74, 9.89, 9.73, 9.77, 9.89, and 9.86 MPa. These values are relatively close to the case of a constant injection rate of 0.3 kg/s, which reaches 9.81 MPa after 50 years.
• Wellbore Pressure: The pressure around the well does not exhibit a simple monotonic change over time. The maximum value occurs at different times for various injection modes, and there is a considerable gap between these maximum values. It is worth noting that Mode 6, Mode 8, and Mode 9 each experienced a decrease in wellbore pressure in different years, which was due to their reduced injection rates. The pressure around the well is affected by both the injection time and injection rate. The pressure around the well becomes higher as the injection time and injection rate increase.

Considering these findings, the injection mode that decreases year by year seems to be the most suitable for this model. This mode results in a maximum wellbore pressure of only 10.93 MPa, making it more conducive to the safe development of the project. This information is valuable for optimizing injection strategies and ensuring the project’s safety and efficacy.

3.3. CO₂ Distribution Analysis

In this section, we discuss the CO₂ distribution. Figure 14 shows the distribution of CO₂ injection volume and rate. Figure 14a shows the distribution of injection rate 0.1 kg/s CO₂ over time. Figure 14b shows the distribution of injection rate 0.5 kg/s CO₂ over time. Figure 14c shows the distribution cloud diagram of CO₂ with an injection volume of 0.1 kg/s; Figure 14d is the distribution cloud diagram of CO₂ with an injection volume of 0.5 kg/s. We can find that as the injection time increases, the diffusion range of CO₂ becomes larger and larger, and the concentration of CO₂ on the diffusion path becomes larger, and the maximum concentration is around the well. And, the greater the injection rate, the wider the diffusion range and the higher the concentration. As the diffusion range increases, the concentration of carbon dioxide becomes smaller and smaller. When the volume fraction of carbon dioxide is less than 1%, we consider it to be no longer diffusing. Due to the effect of gravity, the CO₂ concentration in the upper layer of the reservoir is greater than that in the lower layer and is distributed in a circular cone. For the five injection rates of 0.1–0.5 kg/s, the maximum diffusion ranges are 596 m, 608 m, 622 m, 621 m, and 640 m, respectively, in 50 years (The distance between the uppermost layer of the reservoir and the well is marked by an orange double arrow). For the other six injection modes, the CO₂ distribution is shown in Figure 15. As the injection time increases, the CO₂ diffusion range increases year by year. The maximum concentration is also around the well and occurs in the year of maximum injection volume. The diffusion ranges under these six working conditions are 683 m, 576 m, 696 m, 690 m, 564 m, and 557 m.

In summary, we should pay more attention to whether CO₂ leakage occurs around the well, and the greater the injection rate, the higher the frequency of attention. For the reservoir after 690 m, monitoring can be relatively reduced.

Although this study helps determine the geological carrying capacity of the injection volume, provides insights into selecting a more appropriate injection mode, and has a good guiding role for engineering, there are still some limitations: Failure to consider the impact of changes in reservoir temperature. The boundary conditions and parameters we used were all from the site and meet the engineering requirements to the greatest extent. The calculation results have good convergence and are in line with the site’s basic understanding of the pressure and CO₂ distribution around the well. It has good engineering guidance analysis. Subsequent engineering development will also be closely integrated with simulation.

4. Conclusions

To precisely anticipate the CO₂ migration pattern, assess the CO₂ storage capacity and formation safety, and ensure the project’s seamless advancement, we amalgamated experiments and simulations to formulate a multi-field coupling CO₂ storage prediction model for a specific block in Shenhua. Our study focused on the carbon storage capacity of the reservoir concerning injection volume and mode. The key findings are as follows:

The average pressure in the formation and around the well rises proportionally to the total volume of injected CO₂, with the pressure around the well within the reservoir exhibiting the most significant increase. Additionally, higher injection rates correlate with elevated reservoir pressures. For instance, injecting at a rate of 0.1 kg/s for 50 years resulted in an approximately 7% increase in the average formation pressure compared to pre-injection levels, accompanied by a 14% increase in the maximum pressure around the well. In contrast, injecting at a rate of 0.5 kg/s for the same duration led to a roughly 38% surge in the average maximum formation pressure and a 68% increase in the maximum pressure around the well compared to pre-injection levels.

In the case of various injection modes, the average formation pressure rises with the total injection volume. After 50 years of injection, the maximum average pressure values in the formation become quite similar. Among the modes, the decreasing injection mode with a 10-year interval results in the smallest maximum pressure value around the well, measuring only 10.93 MPa.

The maximum concentration of CO₂ within the reservoir is concentrated around the well, and the extent of CO₂ diffusion expands with the cumulative injection volume. Larger injection rates per unit time led to higher maximum concentrations of CO₂ around the well, increased concentrations along the diffusion path, and broader diffusion ranges. The maximum diffusion range remains under 690 m. Enhanced CO₂ leakage monitoring is recommended around the well and within a 690 m radius from the well.