Study on the Seepage and Stability of Carbon Dioxide Geological Storage in Scrap Metal Mine: A Finite Difference Analysis

Abstract

The carbon dioxide geological storage is one of the most effective ways to slow down the earth's climate warming. The ideal location for storage is in the deep scrap mines. Because of the characteristics of the deep underground condition, the long-term mechanical stability must have close relationship with the geostress fields and permeability, whose development can determine and reflect the mechanical behaviors of carbon dioxide geological storage. The permeability of host rock is affected by existing joints and discontinuity surfaces. Furthermore, as the liquid carbon dioxide is injected into the geological storage, the injection pressure and buoyancy force will change the state of geostress obviously. In this paper, based on mechanical properties of the host rock (granite), coupled brittle shear constitutive model was adopted and the analysis module was implemented in software FLAC3D. Then, with the in situ survey of joints in the host rock and the monitoring of in situ geostress, the permeability characteristics and geostress inverse analysis were input in FLAC3D. At last, fluid-solid coupled model was used to analyze the seepage characteristics of each phase after the injection and long-term storage of liquid carbon dioxide. Therefore, the systematic analysis of the storage stability of carbon dioxide geological could be carried out.

1. Introduction

With a growing number of fossil fuels (coal, petroleum, and gas) that are used in industry and daily life, the greenhouse effect and atmospheric pollution are threatening the human survival environment gravely. So the effective methods for capturing and depositing the carbon dioxide have become an extremely urgent problem to be studied. Many scholars in the world have put forward huge constructive methods and technology. Murai and Fujioka discussed the challenges to the carbon dioxide capture and storage (CCS) technology and listed the issues for actual technical application in the engineering [1]. Herzog scaled up the carbon dioxide capture and storage from megatons to gigatons [2]. Wareing et al. raised a composite equation of state for the modeling of sonic carbon dioxide jets in carbon capture and storage scenarios and predicted the thermodynamic physical properties of carbon dioxide by a novel composite two-phase method [3]. Based on different locations and methods, the underground storage can be divided into three categories: ocean storage [4, 5], carbon sink storage [6 –8], and geological storage [9 –12].

Researches show that the most effective method is the carbon dioxide underground geological storage, where the carbon dioxide is injected into a natural pore in strata or engineering cavity in the deep underground. Due to better closure condition, the abandoned oil, gas layer, coal bed, scrap metal mine, salt mine, and other huge mining caverns in deep underground could be regarded as the potential geological storage sites. It is obvious that the key technique of geological storage is long-term safety and stability of sealed structure, to protect the liquid carbon dioxide from permeating into atmosphere or freshwater aquifer. Many scholars adopted huge method to make sure of the stability of structures. Vasco et al. monitored surface deformation to reverse the fluid flow associated with an injected mass of 3 million tons of carbon dioxide storage and obtained a measurable surface displacement of approximately 5 mm/year [13]. Zimmer et al. monitored the carbon dioxide flux at the Ketzin carbon dioxide storage test site and gained adequate long-term baseline data on the carbon dioxide flux variations, carbon dioxide soil gas flux, soil moisture, and temperature measurements by conducting once a month during a 6-year period [14]. Xia et al. used the preparation of sulfur-doped microporous carbons method to protect the stability of hydrogen and carbon dioxide storage [15]. Rucci et al. combined interferometric synthetic aperture radar (InSAR) data to estimate both quasivertical and quasi-east-west displacements of carbon dioxide geologic storage in central Algeria [16]. Galán and Aparicio studied the role of clays as sealing materials in the geological storage of carbon dioxide by experimental method and furthered the knowledge of reactions that occurred between carbon dioxide, diffused or escaped from a geological storage, and the sealing rock [17].

With the fast development of simulation technique, numerical simulation has become an important prediction method in engineering, because of its high accuracy and less cost. Also, the fast development of computer technique makes many scholars use software to solve the stability problems of geological storage of carbon dioxide, such as TOUGH2, TOUGHREACT, and FLAC3D [18 –26]. Although huge efforts have been done on coal bed and salt mine, there are not enough studies on the scrap metal mine. Because of different geological conditions, there must be much difference in host rocks and mechanical characteristics between the metal mine and other disposal structures. For the candidate geological storage site, the studies on stability and mechanical behaviors of metal mine structures are necessary. For the deep metal mines, the long-term mechanical stability must have close relationship with the in situ geostress and permeability, whose developments can determine and reflect the mechanical behaviors of carbon dioxide geological storage. The permeability of host rock is affected by the existing joints and discontinuity surfaces. Also, as the liquid carbon dioxide is injected into the geological storage, the injection pressure and buoyancy force will change the state of geostress obviously. In this paper, based on mechanical property of the host rock (granite), coupled brittle shear constitutive model was adopted and the analysis module was implemented in software FLAC3D. Then, the brittle to shear coupled model with the in situ survey of joints in the host rock, the monitoring of in situ geostress, and the permeability characteristics were input in FLAC3D to predict the geostress inverse fields. At last, fluid-solid coupled model with combining geostress inverse fields was used to analyze the seepage characteristics of each phase after the injection and long-term storage of liquid carbon dioxide. Therefore, the systematic analysis on the seepage fields and the storage stability of carbon dioxide geological could be carried out. Those results could direct the in-situ carbon dioxide injection, upper seepage field analysis, and damage points estimation.

2. Constitutive Model for Simulation

2.1. Elastic Model

According to the generalized Hooke's law

σ_{i j} = 2 G ε_{i j} + 3 λ ε_{k k} δ_{i j},

(1)

where σ_ij is the stresses, ∊_ij is the total strain, ∊_kk is the normal strain, λ = Eν/(1 + ν)(1 − 2ν), G = E/2(1 + ν), and δ_ij is Kronecker's delta, when i = j, the value is 1, and when i ≠ j, its value is 0.

2.2. Brittle to Shear Coupled Constitutive Model for Geostress Analysis

According to numerous studies on hard rocks, there will be typical brittle to shear damage process on the hard rock under different confining pressures [27 –31], such as granite.

Based on the rock fracture mechanics, generalized plastic mechanics, and macroscopic and mesoscopic rock mechanics, Li et al. analyzed the mechanism of brittle to shear under different confining pressures and adopted contrastive analysis method to acquire limitations in Griffith model and Mohr-Coulomb model in brittle to shear damage phenomenon of hard rock and then raised the brittle to shear coupled constitutive model, whose strength formula was [32]

σ_{1} = b (C_{1} σ_{3} + C_{2}) + (1 - b) (D_{1} σ_{3} + D_{2}),

(2)

where C₁, C₂, D₁, and D₂ were the strength parameter of hard rock, C₁, C₂ were used to express the brittle characteristics, and D₁, D₂ showed the shear characteristics, which could be obtained from the rock mechanics tests. b was the representation of the coefficient from the brittleness to shear, and its value was from 0 to 1.

For the material yield criterion, Li et al. gave out the expression form [32]

\begin{matrix} f = σ_{1} - b (C_{1} σ_{3} + C_{2}) - (1 - b) (D_{1} σ_{3} + D_{2}), \\ b = \exp (m σ_{3} + n), \end{matrix}

(3)

where m and n were the coefficients for calculating b. In this paper, according to mine-by tunnel, the coefficients are as follows: m = 28.11 and n = 1 [33, 34].

2.3. Biot Seepage Mechanics Coupling Theory

Biot gave out state variable between the pore pressure p and the change of liquid volume, which contained basic formulas as follows [35]:

\begin{matrix} Equilibrium equation: σ_{i j, j} + ρ X_{j} = 0 (i, j = 1,2, 3), \\ Geometric equation: ε_{i j} = \frac{1}{2} (u_{i, j} + u_{j, i}) ε_{v} = ε_{11} + ε_{22} + ε_{33}, \\ Constitutive equations: σ_{i j}^{'} = σ_{i j} - a p δ_{i j} = λ δ_{i j} ε_{v} + 2 G ε_{i j}, \\ Δ n = \frac{p}{Q} - a ε_{v} = \frac{p}{R} - \frac{σ_{i j}}{3 H}, \\ Seepage equation: K_{i j} \nabla^{2} p = \frac{1}{Q} \frac{\partial p}{\partial t} - a \frac{\partial ε_{v}}{\partial t}, \end{matrix}

(4)

where p was the pore pressure; Δn was the variation of pore; a was the pore pressure coefficient; H, R, Q, and a were Biot constants; δ was Kronecker constant; K_ij was permeability coefficient; σ_ij was the total stress; σ_ij′ was effective stress; ∊_ij was the total strain; G was the shear modulus; λ was coefficient; ρ was the liquid density.

3. Monitoring Data and Parameter Analysis for Simulation

3.1. In Situ Monitoring Geostress

Geostress is a natural stress which lies in earth layer, and it is also noted as rock initiation stress, absoluteness stress, and host rock stress. Geostress is gradually formed in a long geological history, as a result of gravitational field and geotectonic stress. It is an important design parameter for deep underground projects and one of the most important factors for the safety and stability of geotechnical engineering. In this paper, 11 groups’ in situ geostress monitoring data were acquired, and the results were shown in Figure 1.

Figure 1:

The regression curve for the relation between σ_{h, max ⁡}, σ_{h, min ⁡}, σ_v, and depth.

The maximum principal in situ stress was almost horizontal. The directions of the maximum horizontal principal in situ stress of 9 measuring points were all located between NW-SE directions, which were basically the same with those of the maximum principal in situ stress in the regional geotectonic stress field. Among the 9 measuring points, seven maximum principal in situ stress directions were almost horizontal and their angles with horizontal plane were 10° or less. Based on the above analysis, it was suggested that the dominated stress of this golden mine deep crustal stress field was the horizontal geotectonic stress, not the gravitational stress.

3.2. The Characteristics from Brittle to Shear of Host Rock

Based on the in situ monitoring data of geostress, it is easy to learn when the depth is 800 m; the vertical principal stress and minimum horizontal principal stress were about 20 Mpa. So, in order to obtain the mechanical characteristics of granite in the deep underground, the laboratory tests were carried under 0 MPa, 5 MPa, 10 MPa, 15 MPa, and 20 MPa confining pressure. The brittle damage rock specimens were shown in Figure 2, and the shear damage was illustrated in Figure 3. The typical stress-strain and AE accumulative counts curves were shown in Figure 4. Through data analysis, each strength threshold was indicated in Table 1 and the relations were summarized in Figure 5.

Table 1:

The experimental strength of different confining pressure under single loading.

No.	Confining pressure (Mpa)	Crack initiation (Mpa)	Crack damage (Mpa)	Peak strength (Mpa)	Residual strength (Mpa)

1	0	34.72	88.73	109.14	/
2	5	54.12	121.52	170.84	42.67
3	10	56.78	128.32	182.11	60.21
4	15	63.76	132.68	196.89	86.53
5	20	76.78	145.32	232.43	93.21

Figure 2:

The brittle damage under low confining pressure.

Figure 3:

The shear damage under high confining pressure.

Figure 4:

The stress-strain and AE accumulative counts under 10 MPa confining pressure.

Figure 5:

The experimental strength of different confining pressure.

As shown in Table 1 and Figure 5, residual strength was nonexistent and crack damage strength was close to peak strength under low confining pressure. It could conclude that the rock specimen went through crack damage then reached to the peak strength and lost strength rapidly at last. Thus, brittleness properties of granite emerged under low confining pressure were expressed. As confining pressure increased, the rock specimen presented shear damage.

3.3. Joints Statistics for Permeability Coefficient

For the joints statistics, the joints statistics number is 985 groups in total, with average 32.9 cm space. According to the tendency, 4 groups of superior directions were chosen, which had accounted for 84% of the total. The tendency, average tendency, average inclination angle, and the percentage of the total statistics for the four groups were shown in Table 2.

Table 2:

The joint characteristics of 4 groups of superior directions.

No.	Tendency	Average tendency	Average inclination angle	Percentage of the total statistics

1	10°∼50°	23.54°	70.52°	16.6%
2	90°∼160°	121.58°	59.17°	29.6%
3	180°∼240°	209.12°	66.97°	16.6%
4	280°∼350°	313.62°	60.54°	21.2%

According to the trend, 2 groups of superior directions, 25° and 295°, were chosen, whose inclination angles were about 40°∼90°, and the joint distribution characteristics were shown in Figure 6.

Figure 6:

The joints distribution characteristics of host rock.

3.4. Permeability Analysis

Based on the seepage theory of continuous medium, FLAC3D 2005 simulated the fluid flow in the solid medium and calculated the seepage flow equation to express the coupling process [36]. Also, it is based on the following assumptions: (1) the seepage in the host rock meets the Darcy law; (2) the fluid in the model is considered as variable fluid; (3) the solid in the model will follow the effective stresses principle in Terzaghis consolidation theory; (4) the fluid-solid coupled process will abide by Biot equation.

In fluid-solid coupled model of FLAC3D, the fluid conduction simulation met the Darcy law and was expressed as

q_{i} = - k_{i j} \hat{k} (s) {(p - ρ_{f} x_{j} g_{j})}_{{}^{'}l},

(5)

where q_i was the flow vector, $\hat{k} (s)$ was permeability coefficient, p was the pore pressure, ρ_f was fluid density, and g_j was gravitational acceleration vector.

For the small deformation analysis, the equation of hydrostatic was

- q_{i, i} + q_{v} = \frac{\partial ζ}{\partial t},

(6)

where ζ was the change of fluid volume and q_v was the volume flow rate.

From the perspective of momentum balance, the formula was

σ_{i j, j} + ρ g_{i} = ρ \frac{d v_{i}}{d t},

(7)

where ρ = (1 − n)ρ_s + nsρ_w was the bulk density, ρ_w was the fluid density, and ρ_s was the solid density.

So, the mathematics constitutive relation of seepage in FLAC3D was

\frac{1}{M} \frac{\partial p}{\partial t} + \frac{n}{s} \frac{\partial s}{\partial t} = \frac{1}{s} (- q_{i, i} + q_{v}) - α \frac{\partial ε}{\partial t} + β \frac{\partial T}{\partial t} .

(8)

For the host rock, which contained M groups’ joints, the permeability coefficient vector could be expressed by the width of the joints, the distance between them, the dip angle α, and the tendency β.

Wu and Zhang defined the ith permeability coefficient vector by taking the north geographic pole as X axis and the east geographic pole as Y axis [37]. The particular formula was

K_{i} = \frac{ρ g b_{i}^{3}}{12 μ S_{i}} [\begin{bmatrix} K_{x x i} & K_{x y i} & K_{x z i} \\ K_{y x i} & K_{y y i} & K_{y z i} \\ K_{z x i} & K_{z y i} & K_{z z i} \end{bmatrix}],

(9)

where K_xxi = 1 − sin²α_icos ⁡ ²β_i, K_xyi = K_yxi = − sin²α_icos ⁡ β_isinβ_i, K_zxi = K_xzi = − sinα_icos ⁡ α_icos ⁡ β_i, K_yyi = 1 − sin²α_isin²β_i, K_zzi = 1 − cos ⁡ ²α_i = sin²α_i, where α_i was the dip angle; β_i was the tendency; b_i was the width of joints and fissures; S_i was the distance between joints and fissures; ρ was the density of liquid; μ was the dynamic viscosity of liquid.

Therefore, the permeability coefficient vector of M groups’ joints in the host rock was as follows:

K = \sum_{i = 1}^{M} K_{i} = \frac{ρ g}{12 μ} \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} [I - a_{i} a_{i}^{}] = \frac{ρ g}{12 μ} [\begin{bmatrix} \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} (1 - a_{z i}^{2}) & - \sum \frac{b_{i}^{3}}{S_{i}} a_{x i} a_{y i} & - \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} a_{x i} a_{z i} \\ - \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} a_{y i} a_{x i} & \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} (1 - a_{y i}^{2}) & - \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} a_{y i} a_{z i} \\ - \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} a_{z i} a_{x i} & - \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} a_{z i} a_{y i} & \sum_{i = 1}^{M} \frac{b_{i}^{3}}{S_{i}} (1 - a_{z i}^{2}) \end{bmatrix}],

(10)

where a_yi = sina_isinβ_i, a_xi = cos ⁡ β_isina_i, a_zi = cos ⁡ a_i, β_i was the ith tendency of joint fracture, and the value was 0 ≤ β ≤ 360°; a_i was the ith dip angle of joint fracture, and the value was 0 ≤ a_i ≤ 90°.

Based on the forgoing analysis, the permeability coefficient would be changed with the joints. In this paper, the varied permeability can be presented by Rutqvist et al. [38]:

K = K_{0} \exp [A (\frac{ϕ}{ϕ_{0}} - 1)],

(11)

where K₀ was the initial (zero stresses) intrinsic permeability; φ₀ was the initial porosity; A was an empirical (calibration) factor which could be obtained in laboratory tests.

4. Simulation Analysis

4.1. Engineering Background

In this paper, a golden mine would be chosen as case study of carbon dioxide geological storage, and the disposal cavity was located in 850 m depth. For the dimension of the cavity, the length and width were 150 m and the average height was 34 (the maximum height was 45 m in the middle).

Based on the actual engineering geology and structures, a 3D model was built in the finite element analysis software of FLAC3D to accept the constitutive parameters, permeability coefficients, as shown in Figure 7. The geometrical dimensions are length 600 m (X direction), width 600 m (Y direction), and height 529 m (Z direction). The finite element mesh consists of 99835 nodes and 96000 elements. And for the cavity which was between 200 m and 400 m (Y direction), it was meshed as 7224 elements (blue zone) (Figure 7).

Figure 7:

A 3D model of host rock and cavity in FLAC3D.

4.2. Simulation Constitutive Model and Parameters

The model was simplified along with actual layers, and in the 850 depth, the host rock is granite. For the geostress inverse simulation, in the elastic stage, the constitutive model was elastic model, which was expressed in Section 2.1. The property of elastic model (bulk modulus, shear modulus, and density) was shown in Table 3. After the elastic stage, the materials will be presented as plasticity, in this simulation, according to the characteristics of hard rock, the brittle to shear coupled model was chosen as plasticity constitutive model, which was shown in Section 2.2. According to expression forms of this constitutive model, Li et al. programmed the brittle to shear coupled constitutive model by FISH language and C++ language into FLAC3D for simulating the hard rock and obtained the integrating degree between the simulation result and laboratory test data which was 0.9836 [32]. In this model, the basic material property of host rock (bulk modulus, shear modulus, and density), the strength parameter of hard rock (C₁C₂D₁, and D₂), and the coefficients (m and n) for calculating b for the brittle to shear coupled model can be obtained by tests, which was shown in Section 3.2. Also, for the host rock parameters, which were input in simulation, were obtained from the test. The results were shown in Table 3.

Table 3:

Calculation parameters of brittle shear coupled model.

Calculation parameters	Bulk (MPa)	Shear (MPa)	Tension (MPa)	Density (Kg/m³)	m	n	C ₁	C ₂	D ₁	D ₂	Dilation(°)

Value	33000	25000	3	2720	28.11	1	1.8752	38.48	5.86	320	30

For the coupled analysis, the Biot coupled model was adopted in FLAC3D. And based on the in-situ survey data in Sections 3.3 and 3.4, the permeability coefficient and initial moisture content rate of host rock were calculated as 3.35e-7 cm/s and 85.13% respectively.

4.3. Boundary Condition and Initial Conditions for Simulation

For the boundary condition, the four lateral surfaces of model were set as unidirectional normal fixed. For the top surface, it was applied by uniformly distributed load, which was equivalent to overlying rock weight. And for the bottom surface, it was set as three directions’ fixed.

For initial conditions, the gravity was set as 9.8 m/s². The initial principal stresses with depth were input into model, and based on 11 groups’ in situ geostress monitoring data, the site regression equation could be obtained by stress relief method:

\begin{matrix} σ_{h, \max} = 1.433 + 0.043 H, \\ σ_{h, \min} = 1.304 + 0.024 H, \\ σ_{v} = 0.07 + 0.028 H, \end{matrix}

(12)

where σ_{h, max ⁡}, σ_{h, min ⁡}, and σ_v are, respectively, the maximum horizontal principal in situ stress, the minimum horizontal principal in situ stress, and the vertical principal stress and their unit is MPa. H is the depth and its unit is m.

Usually, when the pressure is larger than 7.38 Mpa, the carbon dioxide will be existed as liquid state with 0.4∼0.6 g/cm³ density. In this simulation analysis, the carbon dioxide geological storage was chosen in the 850 m deep underground. With this high pressure situation, in order to keep the liquid state after permeating in the host rock and by referring to the strength of host rock, in this paper, the injecting pressure of carbon dioxide was 10 times of liquefaction critical pressure; it was 73.8 Mpa. Therefore, after excavation, the initial pore pressure 73.8 MPa was applied on the storage wall to simulate the seepage and stability of storage.

4.4. Inverse Simulation of the Geostress

Based on the in situ monitoring data of geostress distribution with depth, which was input into the model, to obtain the inverse geostress fields, and because of different definitions of maximum and minimum principal stress between the software and rock mechanics, the maximum principal stress in rock mass was the minimum principal stress in FLAC3D, and it was also reverse for the minimum principal stress. And the minus value was representation of the compressive stress; the positive value was on behalf of tensile stress, as depicted in Figures 8 and 9.

Figure 8:

The contour of maximum principal stress field.

Figure 9:

The contour of minimum principal stress field.

As shown in Figures 8 and 9, it was easy to obtain that the maximum principal stress would become larger with the depth increase, and the minimum principal stress also followed the same rules. And during the ore excavation period, the geostress also could be affected by excavation. After simulation, the geostress fields were illustrated in Figures 10 and 11.

Figure 10:

Maximum principal stress field after excavation.

Figure 11:

Minimum principal stress field after excavation.

As shown in Figure 10, the maximum principal stress were the compressive stress, which occurred at upper corner of cavity. The scope was 62.5 m width. Contrastive analysis between Figure 8 and Figure 10. The maximum principal stress value in Figure 10 was lower than it in Figure 8, which was changed from 40 MPa∼42.5 MPa to 30 MPa∼35 MPa. And there was a little improved area with 2.5 MPa increased value on the top roof. For the minimum principal stress by contrastive analysis of Figures 9 and 11, there was an obviously raised regions on the lateral bottom wall with 15 m width; the value was increased from 6 MPa∼8 MPa to 18 Mpa∼19.43 Mpa, and it was because that the rock mass was dug out, and the lateral walls have to undertake the supporting capacity of excavation body. Also, because the rock mass was excavated to lead the stress released, there was an obviously compressive stress decreasing region on the bottom floor, and the value was from 8 MPa∼10 MPa to 0 MPa to 4 MPa, in this area, the host rock even presented as tensile damage state. For the top roof, the reduced level of compressive stress was also very huge (from 6 MPa∼8 MPa to 0 MPa, even, tensile state), and it was due to vault construction, which can transmit the compressive stress to lateral walls, and there was also some stress released reason.

The rock mass was excavated, which would cause some damage in the host, as shown in Figure 12.

Figure 12:

The damage zone contour after excavation.

As shown in Figure 12, all the damage zones were the tensile failure, and they mainly distributed on the top of lateral wall and top roof. The damage distance was 12.5 m (width of an element).

4.5. Seepage Analysis

After simulating the seepage stability, the periodical penetrative contours were shown in Figures 13, 14, 15, 16, 17, 18, 19, and 20, and the days and its corresponding center pore pressure were listed in Table 4.

Table 4:

The disposal days and its corresponding center maximum pore pressure.

Disposal days	1	5	10	20	30	50	100	200	300

The max. pressure (MPa)	73.8	73.8	73.742	73.136	71.814	67.584	55.850	39.611	37.243

Figure 13:

The penetrative contours of 1 day.

Figure 14:

The penetrative contours of 5 days.

Figure 15:

The penetrative contours of 10 days.

Figure 16:

The penetrative contours of 30 days.

Figure 17:

The penetrative contours of 50 days.

Figure 18:

The penetrative contours of 100 days.

Figure 19:

The penetrative contours of 200 days.

Figure 20:

The penetrative contours of forever (more than 300 days).

Based on the forgoing coupling simulation results, after 300 days of diffusion and permeation, the pore pressure of liquid carbon dioxide became stable. The largest penetrative distance was about 70 m. The center pressure was 37 MPa, which also could make the carbon dioxide keep in liquid state. As shown in periodical contour of permeation, at the initial stage after liquid carbon dioxide was injected into cavity (from 1 days, 5 days, 10 days, and 30 days to 50 days), the diffusion and permeation were faster, due to high pressure in center and low initial pore pressure of host rock. At 1 day, because of excavation damage, the porosity was enlarged; the liquid carbon dioxide permeated about 25 m just after was injected into the cavity. And till 50 days, the diffused distance was nearly 50 m; the average velocity of permeability was 0.5 m per day. Then the penetrative rate was becoming lower and lower than before. From 50 days to 100 days, the permeated distance was 8 m; the average rate was 0.16 m per day. And for the next 100 days, the diffused distance was 10 m; the rate was 0.1. And for last 100 days, the distance was only 2 m; the rate was 0.02 per day, then, the diffusion became very slow and went to stability after 300 days. For the center pore pressure, because of the high injection pressure in the initial stage, the center pressure was kept about 10 days, and then it declined fast till 200 days. After 200 days, the velocity of pressure decline was becoming lower and lower; it was due to the low permeability. The particular change curve of center pore pressure with days change was shown in Figure 21.

Figure 21:

The change curve of center pore pressure with days.

4.6. Stability Analysis

After the liquid carbon dioxide was injected into the storage, the high pressure also has an effect on the geostress distribution. After coupling simulation, the maximum principle stress and minimum principle stress were shown in Figures 22 and 23, respectively.

Figure 22:

Maximum principal stress field after injection.

Figure 23:

Minimum principal stress field after injection.

As shown in Figure 22, the maximum principal stresses were all the compressive stress. Because of injected pressure which led the stress to be redistributed, the compressive stress was increased obviously on the both top roof and bottom floor. The maximum value (48 MPa) occurred at 30 m depth under bottom storage; the compressive stress value was increased about 8 MPa with comparison of Figure 10. And for the host rock above the storage, the stress was increased in 25 m distance from the top roof; above this region was an decreased area till to the top. For the minimum principal stress by contrastive analysis of Figures 11 and 23, there were obviously raised regions on the bottom floor with 60 m depth, the value was increased from 0 MPa to 25 Mpa∼28 Mpa, and it was because of the pressure and weight of liquid. Also, because the stress was redistributed, there was an obviously compressive stress decreasing region on the lateral walls; the value was from 18 MPa to 0 MPa to 4 MPa, in this area, the host rock even presented as tensile damage state.

After the liquid carbon dioxide was injected into the storage, the high pressure could produce the tensile stress. It was due to the rock tensile strength was much lower than rock compressive strength, there was some tensile damage in the lateral walls. The damage state zone contour after diffused stability was shown in Figure 24.

Figure 24:

The damage zone contour after stability.

As shown in Figure 24, with high pressure liquid, it is easy to learn that the mainly damage zone was enlarged than before. The largest damage range was 75 m on the lateral walls. Also, we could get that most tensile damage occurred in the past, and it means that most damage appeared at the early stage after injection. The damage development and the penetration were a promoting process with each other. This is also a reason for the high permeability during early stage. In Figure 25, it is also due to the high pressure, and the maximal displacement is 4.9242e-6 m, which is much smaller than the real size of storage, so it is negligible.

Figure 25:

The displacement contour after stability.

5. Conclusion

Based on the 985 groups’ joints statistic, the average distance was 32.9 mm, and after calculating, the permeability coefficient was 3.35e-7 cm/s. After the 11 locations in situ geostress monitoring, the relationship between the principal stress and depth could be obtained. Also, through laboratory test, the characteristic and mechanical parameters could be acquired. Those data were necessary and enough accurate for coupling simulation.

By contrastive analysis of inverse geostress between before and after excavation, the maximum principal stress on the corner between the top roof and lateral walls became lower, because rock mass excavation led to released and redistributed stress. For the minimum principal stress, on the lateral bottom wall with 15 m width, the value was increased, and it was because the rock mass was dug out, and the lateral walls have to undertake the supporting capacity of excavation body. Also, there were obviously compressive stress decreasing region on the bottom floor and top roof, in those area, the host rock even presented as tensile damage state.

For the contrastive analysis of geostress between before and after carbon dioxide injection, because of injected pressure to lead the stress redistributed, the compressive stress was increased obviously on the both top roof and bottom floor. For the minimum principal stress, there were obviously raised regions on the bottom floor with 60 m depth, and it was because of the pressure and weight of liquid. Also, because the stress was redistributed, there was an obviously compressive stress decreasing region on the lateral walls.

Through the damage state analysis, after excavating, there were some damage zones, which were mainly distributed on the top of lateral wall and top roof with 12.5 m depth. After the liquid carbon dioxide is injected, the damage zones were enlarged a lot, and they were distributed on the lateral walls with 75 m width. This was also an important reason for the high permeability around the storage structure.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

Murai

and Fujioka

, “Challenges to the Carbon dioxide Capture and Storage (CCS) technology,” IEEJ Transactions on Electrical and Electronic Engineering, vol. 3, no. 1, pp. 37–42, 2008.

Herzog

H. J.

, “Scaling up carbon dioxide capture and storage: from megatons to gigatons,” Energy Economics, vol. 33, no. 4, pp. 597–604, 2011.

Wareing

C. J.

Woolley

R. M.

Fairweather

, and Falle

S. A. E. G.

, “A composite equation of state for the modeling of sonic carbon dioxide jets in carbon capture and storage scenarios,” AIChE Journal, vol. 59, no. 10, pp. 3928–3942, 2013.

Toledo-Cervantes

Morales

Novelo

, and Revah

, “Carbon dioxide fixation and lipid storage by Scenedesmus obtusiusculus,” Bioresource Technology, vol. 130, pp. 652–658, 2013.

Heap

A. D.

Nichol

S. L.

, and Brooke

B. P.

, “Seabed mapping to support geological storage of carbon dioxide in offshore Australia,” Continental Shelf Research, vol. 83, pp. 108–115, 2014.

Eftekhari

A. A.

van der Kooi

, and Bruining

, “Exergy analysis of underground coal gasification with simultaneous storage of carbon dioxide,” Energy, vol. 45, no. 1, pp. 729–745, 2012.

Firouzi

and Wilcox

, “Molecular modeling of carbon dioxide transport and storage in porous carbon-based materials,” Microporous and Mesoporous Materials, vol. 158, pp. 195–203, 2012.

Wee

J. H.

, “A review on carbon dioxide capture and storage technology using coal fly ash,” Applied Energy, vol. 106, pp. 143–151, 2013.

Oladyshkin

Class

Helmig

, and Nowak

, “A concept for data-driven uncertainty quantification and its application to carbon dioxide storage in geological formations,” Advances in Water Resources, vol. 34, no. 11, pp. 1508–1518, 2011.

10.

Rempel

K. U.

Liebscher

Heinrich

, and Schettler

, “An experimental investigation of trace element dissolution in carbon dioxide: applications to the geological storage of CO₂,” Chemical Geology, vol. 289, no. 3–4, pp. 224–234, 2011.

11.

Eshiet

and Sheng

, “Investigation of geomechanical responses of reservoirs induced by carbon dioxide storage,” Environmental Earth Sciences, vol. 71, no. 9, pp. 3999–4020, 2014.

12.

Kim

A.-R.

Cho

G.-C.

, and Kwon

T.-H.

, “Site characterization and geotechnical aspects on geological storage of CO₂ in Korea,” Geosciences Journal, vol. 18, no. 2, pp. 167–179, 2014.

13.

Vasco

D. W.

Rucci

Ferretti

, “Satellite-based measurements of surface deformation reveal fluid flow associated with the geological storage of carbon dioxide,” Geophysical Research Letters, vol. 37, no. 3, Article ID L03303, 2010.

14.

Zimmer

Pilz

, and Erzinger

, “Long-term surface carbon dioxide flux monitoring at the Ketzin carbon dioxide storage test site,” Environmental Geosciences, vol. 18, no. 2, pp. 119–130, 2011.

15.

Xia

Zhu

, and Tang

, “Preparation of sulfur-doped microporous carbons for the storage of hydrogen and carbon dioxide,” Carbon, vol. 50, no. 15, pp. 5543–5553, 2012.

16.

Rucci

Vasco

D. W.

, and Novali

, “Monitoring the geologic storage of carbon dioxide using multicomponent SAR interferometry,” Geophysical Journal International, vol. 193, no. 1, pp. 197–208, 2013.

17.

Galán

and Aparicio

, “Experimental study on the role of clays as sealing materials in the geological storage of carbon dioxide,” Applied Clay Science, vol. 87, pp. 22–27, 2014.

18.

Apps

J. A.

, and Pruess

, “Mineral sequestration of carbon dioxide in a sandstone-shale system,” Chemical Geology, vol. 217, no. 3–4, pp. 295–318, 2005.

19.

Doughty

, “Modeling geologic storage of carbon dioxide: comparison of non-hysteretic and hysteretic characteristic curves,” Energy Conversion and Management, vol. 48, no. 6, pp. 1768–1781, 2007.

20.

Rutqvist

Birkholzer

Cappa

, and Tsang

C.-F.

, “Estimating maximum sustainable injection pressure during geological sequestration of CO₂ using coupled fluid flow and geomechanical fault-slip analysis,” Energy Conversion and Management, vol. 48, no. 6, pp. 1798–1807, 2007.

21.

Celia

M. A.

and Nordbotten

J. M.

, “Practical modeling approaches for geological storage of carbon dioxide,” Ground Water, vol. 47, no. 5, pp. 627–638, 2009.

22.

Laforce

T. C.

, and Blunt

M. J.

, “A three-phase four-component streamline-based simulator to study carbon dioxide storage,” Computational Geosciences, vol. 13, no. 4, pp. 493–509, 2009.

23.

Xiang

Cao

Lan

Wang

, and Broom

D. P.

, “Multiscale simulation and modelling of adsorptive processes for energy gas storage and carbon dioxide capture in porous coordination frameworks,” Energy and Environmental Science, vol. 3, no. 10, pp. 1469–1487, 2010.

24.

Kolditz

Bauer

Böttcher

, “Numerical simulation of two-phase flow in deformable porous media: application to carbon dioxide storage in the subsurface,” Mathematics and Computers in Simulation, vol. 82, no. 10, pp. 1919–1935, 2012.

25.

and Jiang

, “A numerical study of the impurity effects of nitrogen and sulfur dioxide on the solubility trapping of carbon dioxide geological storage,” Applied Energy, vol. 128, pp. 60–74, 2014.

26.

Pan

P. Z.

Rutqvist

Feng

X. T.

, and Yan

, “An approach for modeling rock discontinuous mechanical behavior under multiphase fluid flow conditions,” Rock Mechanics and Rock Engineering, vol. 47, no. 2, pp. 589–603, 2014.

27.

Montési

L. G. J.

, “A constitutive model for layer development in shear zones near the brittle-ductile transition,” Geophysical Research Letters, vol. 34, no. 8, Article ID L08307, 2007.

28.

Zhou

J. W.

W. Y.

, and Yang

X. G.

, “A microcrack damage model for brittle rocks under uniaxial compression,” Mechanics Research Communications, vol. 37, no. 4, pp. 399–405, 2010.

29.

Molladavoodi

and Mortazavi

, “A damage-based numerical analysis of brittle rocks failure mechanism,” Finite Elements in Analysis and Design, vol. 47, no. 9, pp. 991–1003, 2011.

30.

Xie

Zhu

Q. Z.

L. H.

, and Shao

J. F.

, “A micromechanics-based elastoplastic damage model for quasi-brittle rocks,” Computers and Geotechnics, vol. 38, no. 8, pp. 970–977, 2011.

31.

Grady

D. E.

, “Adiabatic shear failure in brittle solids,” International Journal of Impact Engineering, vol. 38, no. 7, pp. 661–667, 2011.

32.

Qiao

Sui

Z.-L.

, and Li

Q.-W.

, “Strength analysis of rock material under the brittle shear failure mode,” Journal of University of Science and Technology Beijing, vol. 34, no. 12, pp. 1364–1370, 2012.

33.

Hoek

Carranza-Tarres

, and Corkum

, “Hoek-Brown failure criterion—2002 edition,” in Proceedings of the 5th North American Rock Mechanics Symposium, Toronto, Canada, 2002.

34.

Cai

and Kaiser

P. K.

, “In-situ rock spalling strength near excavation boundaries,” Rock Mechanics and Rock Engineering, vol. 47, no. 2, pp. 659–675, 2014.

35.

Biot

M. A.

, “General theory of three-dimensional consolidation,” Journal of Applied Physics, vol. 12, no. 2, pp. 155–164, 1941.

36.

Inc

I. C. G.

, FLAC3D (Fast Lagrangian Analysis of Continuain 3 Dimensions) Online Manual Table of Contents, Version 5.0, 2005.

37.

Y. Q.

and Zhang

Z. Y.

, An Introduction to Rock Mass Hydraulics, Southwest Jiaotong University Press, Chengdu, China, 1995.

38.

Rutqvist

Y.-S.

Tsang

C.-F.

, and Bodvarsson

, “A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 39, no. 4, pp. 429–442, 2002.