Numerical Simulation and Modeling on CO2 Sequestration Coupled with Enhanced Gas Recovery in Shale Gas Reservoirs

CO2 geological sequestration in shale is a promising method to mitigate global warming caused by greenhouse gas emissions as well as to enhance the gas recovery to some degree, which effectively addresses the problems related to energy demand and climate change. With the data from the New Albany Shale in the Illinois Basin in the United States, the CMG-GEM simulator is applied to establish a numerical model to evaluate the feasibility of CO2 sequestration in shale gas reservoirs with potential enhanced gas recovery (EGR). To represent the matrix, natural fractures, and hydraulic fractures in shale gas reservoirs, a multicontinua porous medium model will be developed. Darcy’s and Forchheimer’s models and desorption-adsorption models with a mixing rule will be incorporated into the multicontinua numerical model to depict the three-stage flow mechanism, including convective gas flow mainly in fractures, dispersive gas transport in macropores, and CH4-CO2 competitive sorption phenomenon in micropores. With the established shale reservoir model, different CO2 injection schemes (continuous injection vs. pulse injection) for CO2 sequestration in shale gas reservoirs are investigated. Meanwhile, a sensitivity analysis of the reservoir permeability between the hydraulic fractures of production and injection wells is conducted to quantify its influence on reservoir performance. The permeability multipliers are 10, 100, and 1,000 for the sensitivity study. The results indicate that CO2 can be effectively sequestered in shale reservoirs. But the EGR of both injection schemes does not perform well as expected. In the field application, it is necessary to take the efficiency of supplemental energy utilization, the CO2 sequestration ratio, and the effect of injected CO2 on the purity of produced methane into consideration to design an optimal execution plan. The case with a permeability multiplier of 1,000 meets the demand for both CO2 sequestration and EGR, which indicates that a moderate secondary stimulation zone needs to be formed between the primary hydraulic fractures of injection and production wells to facilitate the efficient energy transfer between interwell as well as to prevent CO2 from channeling. To meet the demand for CO2 sequestration in shale gas reservoirs with EGR, advanced and effective fracking is essential.


Introduction
At present, fossil fuels are the primary source of energy consumption in the world [1][2][3], which leads to abundant CO 2 being released into the atmosphere via the combustion of fossil fuels (oil, natural gas, and coal). Before the industrial revolution, the average concentration of CO 2 in the atmosphere was 0.03%, which has increased to 0.04% in 2005 and is expected to reach 0.01% by 2100 without any intervention [4][5][6]. The release of huge amounts of CO 2 into the atmosphere leads to global warming and ocean acidification, which will be harmful to the whole world. Therefore, mitigating the contribution of CO 2 emissions to global warming has become a common problem faced by all the countries. It is expected that China's total CO 2 emissions will reach 67 × 10 8 t in 2030 and surpass the United States, to become the world's largest CO 2 emitter [7,8]. CCS, as an emerging technology, is expected to play a key role in CO 2 emission reduction [9,10]. From 2010 to 2050, 14% of the CO 2 emission reduction arises from the application of CCS, making it the largest contributor to GHG emission reduction [11].
With the booming of commercial shale gas development, many researchers focus on the feasibility study of CO 2 sequestration in shale gas reservoirs with enhanced gas recovery [12][13][14]. Based on different shale samples, several researchers have extensively studied the interaction between CO 2 and CH 4 in shales. Nuttal et al. [15] systematically investigated the sorption capacity of CH 4 and CO 2 in the Ohio Shale of the Upper Devonian in eastern Kentucky. And the results indicated that the organic matter of the Ohio black shale has a complex microporous structure similar to the coal, which can facilitate the adsorption of large amounts of gas. With organic-rich shale samples from the Fort Worth Basin, Kang et al. [16] tested the sorption capacity for CO 2 and CH 4 . The lab data demonstrated that 40% more CO 2 was preferentially adsorbed. Compared to CH 4 , the preferential sorption of CO 2 by shale is beneficial for the enhanced recovery of shale gas with the CO 2 injection while effectively sequestering a certain amount of CO 2 . But the actual reservoir performance is highly related to the reservoir geological conditions and engineering parameters. Godec et al. [17] showed that the primary recovery of shale gas in the Marcellus Shale reservoir or shale reservoirs similar to the Marcellus Shale is 20%-35%, with an average recovery of 25%. With an appropriate well spacing between injector and producer, 7% of recovery increment can be obtained with the CO 2 injection. Zhang [18] conducted in-house experiments on CO 2 injection enhancing shale gas recovery for terrestrial shale samples from Fuxian of Ordos Basin. The lab results showed that the recovery of CO 2 injection is 80.29%, 7.66% higher compared to the depleted development.
However, as to the Devonian and Mississippian Shales in the Illinois Basin and the Silurian Ohio Shale in eastern Kentucky, the EGR with CO 2 injection does not perform well as expected, with a recovery increment less than 1%. In the case of Devonian and Mississippian Shales, the energy transfer efficiency is impeded by the tight and unstimulated matrix between the hydraulic fractures of the injector and producer. In the case of Ohio Shale, the reservoir pressure is still high while injecting CO 2 , limiting the amount of CO 2 injection and the corresponding shale gas recovery. In China, there are some CCS projects on CO 2 sequestration in depleted oil reservoirs or saline aquifers, such as the Shenhua Ordos pilot-scale project for CO 2 deep saline aquifer storage [19]. The CS-EGR in shale reservoirs is still in the preliminary stage even with abundant shale resources [20]. But the CS-EGR in shale reservoirs has attracted extensive attention recently in China [21,22]. There are still many challenges to prove the viability of sequestration and enhanced recovery in shale reservoirs because of the complicated mechanisms of the process and the engineering complexity of CO 2 injection in shale reservoirs.
In the paper, to objectively evaluate the feasibility of CO 2 sequestration in shale gas reservoirs with potential enhanced gas recovery (EGR), numerical simulation studies are carried out to investigate the mechanism of the process and the effects of several dominating engineering parameters on the reservoir performance to explore the engineering complexity of CO 2 injection in shale reservoirs, taking the New Albany Shale as an example. A multicontinua porous medium model will be developed to represent the domains (matrix, natural fractures, and hydraulic fractures) in shale gas reservoirs. A different domain has its own scale and corresponding transport mechanism. Darcy's and Forchheimer's models and desorption-adsorption models with a mixing rule will be incorporated into the multicontinua numerical model to mimic gas migration in different domains. The gridding scheme, local grid refinement with logarithmic spacing, is employed to accurately simulate the detailed transient gas flow phenomenon around hydraulic fractures. With the established shale reservoir model, different CO 2 injection schemes (continuous injection vs. pulse injection) for CO 2 sequestration in shale gas reservoirs are investigated. Meanwhile, a sensitivity analysis of the reservoir permeability between the hydraulic fractures of the producer and injector is conducted to quantify its influence on reservoir performance. The insights obtained from the study will not only improve the understanding of CO 2 sequestration in shale reservoirs but also provide new methods for enhancing shale gas recovery, which effectively promotes the technology development and wide application of CCS in China and across the world.

Numerical Simulation Methodology.
With the data from the New Albany Shale, the CMG-GEM simulator is implemented to establish a reservoir model. The governing equations employed in the CMG's general EOS compositional simulator (GEM), which depicts the total mass balance for each component including accumulation term as well as convection term and sink/source term, are expressed by the continuity equations below [24]: where ρ k=o, g,w denotes the density of phase k, where k represents the phases (o is oil phase, g is gas phase, and w is water phase); v k=o, g,w is Darcy's flow velocity of each phase; s k=o, g,w is the saturation of each phase; y i is the mole fraction of 3 Geofluids component i in the gas phase; x i is the mole fraction of component i in the oil phase; ∅ is the porosity; and q i denotes the injection/production of component i.
To meet the needs of thermodynamic equilibrium, the Peng-Robinson equation of state is generally applied in the GEM to determine the component composition and compressibility factor for each phase.
The Langmuir isotherm has been widely employed to simulate single component adsorption: where VðPÞ is the gas volume of adsorption at pressure P; V L is the Langmuir volume, referred to as the maximum adsorbed gas volume at the infinite pressure; and P L is the Langmuir pressure, representing the pressure corresponding to a one-half Langmuir volume. For modeling the competitive multicomponent adsorption-desorption process, an extended Langmuir isotherm is implemented [25]: where w i is the moles of adsorbed component i per unit mass or rock; w i,max is the maximum moles of adsorbed component i per unit mass or rock; B i is the parameter for Langmuir isotherm relation; P is the pressure; and y ig is the molar fraction of adsorbed component i in the gas phase. Both B i and ω i,max are parameters of the Langmuir isotherm for single component i (CH 4 and CO 2 ), which are determined in the lab with the New Albany Shale samples.

Geofluids
For the simulation of shale reservoirs developed with multistage hydraulic fractured horizontal wells, some more equations or models are applied to address the specialty. By implementing the correlation proposed by Evan and Civan, the Forchheimer model with the non-Darcy beta factor can be utilized to simulate a turbulent gas flow within hydraulic fractures, accounting for the inertial effects on the flow characteristics [26]: where v is the velocity, K is the permeability, μ is the viscosity, ρ is the density, P is the pressure, and β is the non-Darcy beta factor, determined by the correlation proposed by Evan and Civan. Local grid refinement with logarithmic spacing, which discretizes the reservoir to a finer degree region around hydraulic fractures and more coarsely further away from the hydraulic fractures, is implemented to accurately depict the detailed transient gas flow phenomenon around the hydraulic fractures. A dual-permeability model is employed to take natural fractures acting as boundaries to matrix ele-ments in three directions into consideration, where the governing equations of the dual-permeability model are an extension of the equations for single porosity systems. There are two sets of mass balance equations, with one for the matrix system and the other one for the natural fracture system. Meanwhile, new terms, accounting for the matrix-  5 Geofluids fracture transfer in each phase for each component, are included in the mass balance equations of the dualpermeability model. The logarithmically spaced, locally refined, and dual permeability (LS-LR-DK) methodology has been widely applied to simulate gas flows in hydraulically fractured shale gas reservoirs, which has been validated by previous work to both accurately and efficiently simulate stimulated fractured shale reservoirs [23,27,28].

Reservoir Model.
Based on the CMG-GEM simulator, a homogeneous 3D multicontinua porous medium model is developed to evaluate the feasibility of CO 2 sequestration in shale gas reservoirs with potential enhanced gas recovery (EGR). The dimensions of the numerical model are 1445 m × 914 m × 30 m, corresponding to the length, width, and thickness of the shale gas reservoir, respectively, as shown in Figure 2. Two horizontal wells are simulated with four fracturing stages each, upon which each fracturing stage has a single perforated interval. The local grid refinement with logarithmic spacing is employed to model hydraulic fractures explicitly in the matrix portion by defining high permeability values for the hydraulic fractures and low permeability values for the shale matrix. The competitive CH 4 -CO 2 adsorptiondesorption is simulated based on the extended Langmuir model.
Due to expensive computational time and cost for the entire field case, in the study, a submodel with one fracking stage for each horizontal well is extracted from the whole reservoir model. The dimensions of the submodel are 164:6 m × 914:4 m × 30:5 m, as shown in Figure 3. The entire simulation period is 30 years. First of all, the shale reservoir is depleted for five years. Methane is recovered from the CH 4 producer drilled in the shale reservoir with a maximum gas rate at a surface condition (STG) of 2:8 × 10 4 m 3 /day and minimum bottom-hole pressure (BHP) of 1379 kPa. With a constant daily injection rate of 1:1 × 10 3 m 3 /day, CO 2 is injected by continuous and pulsed injection, respectively. The injection time is 5 years. The specific parameters used in the numerical model are listed in Table 1.

Simulation Schemes.
With the established shale reservoir model, different CO 2 injection schemes (continuous injection vs. pulse injection) for CO 2 sequestration in shale gas reservoirs are investigated. The following three schemes are simulated. Scheme 1 is a depletion development scheme without CO 2 injection. In scheme 2 and scheme 3, CO 2 is injected from the 5th year to the 10th year. Scheme 2 is continuous injection, and scheme 3 is pulse injection.
Scheme 1: depletion development. Based on this scheme, 9:52 × 10 6 m 3 of CH 4 is produced over 30 years, as shown in Figure 4.
Scheme 2: continuous injection. In scheme 2, CO 2 is injected continuously from the 5th year to the 10th year. During the five-year injection, the cumulative injection of CO 2 is 2:07 × 10 6 m 3 through the one-stage hydraulic fracture. And the total gas production is 9:56 × 10 6 m 3 at the end of the simulation. The gas production of each component is CH 4 (9:31 × 10 6 m 3 ) and CO 2 (0:25 × 10 6 m 3 ), respectively, as shown in Figure 5. The reservoir pressure is improved by CO 2 injection, resulting in higher total gas production. Compared to scheme 1, the decrease of CH 4 production of scheme 2 is due to the change of composition by CO 2 injection, reducing the purity of produced methane. The EGR of the scheme does not perform well as expected. Most of the supplemental energy is trapped around an injector due to a tight formation impeding the effective pressure communications between an injector and a producer, dominating the success of the EGR, as shown in Figure 6 6 Geofluids why the reservoir productivity is enhanced by the injection process. But the increment is not substantial. Overall, the supplemental energy is not effectively utilized to offset the impact of CO 2 injection on the purity of produced methane. But the injected CO 2 is effectively sequestered. At the end of the 30-year simulation, 87.9% of CO 2 is still effectively sequestered in shale reservoirs. Scheme 3: pulse injection. In scheme 3, CO 2 is injected by pulse injection from the 5th year to the 10th year, pulse injection starting from the first month of the 5th year, with CO 2  7 Geofluids injection for one month and then shutting in for one month repeatedly throughout the 5 years. During the five years, 1:04 × 10 6 m 3 of CO 2 is injected through the hydraulic fracture. And the total gas production is 9:53 × 10 6 m 3 . The gas production of each component is CH 4 (9:33 × 10 6 m 3 ) and CO 2 (0:20 × 10 6 m 3 ), as shown in Figure 7. Like scheme 2, injected CO 2 is effectively sequestered. At the end of the 30-year simulation, 80.8% of CO 2 is still effectively sequestered in the shale reservoir. With the CO 2 injection at half the amount of scheme 2, the total gas production of scheme 3 is basically the same as scheme 2 with the difference of 0:03 × 10 6 m 3 . The difference in total gas production between scheme 2 and scheme 3 is small, mainly due to the low efficiency of supplemental energy utilization. Most of the  8 Geofluids supplemental energy is trapped around an injector, which cannot be fully harnessed for efficient development, as shown in Figure 6(b). Meanwhile, the cumulative injection of CO 2 is halved for the pulse injection, which reduces its impact on produced natural gas purity, resulting in higher CH 4 produc-tion in scheme 3 compared to scheme 2. The simulation results of the three schemes are summarized in Table 2 Figure 7: (a) Simulation results of scheme 3: cumulative gas production, CO 2 injection; (b) cumulative production of each component.  10 Geofluids on the primary factors affecting shale gas development with CO 2 injection. Among them, reservoir permeability has the greatest impact on reservoir performance. Based on the above numerical simulation results, it can be observed that the tight and unstimulated shale matrix between the hydraulic fractures of the producer and injector seriously affects the energy transfer efficiency, resulting in the inefficient development of shale gas. Therefore, the sensitivity analysis of inter-fracture reservoir permeability (reservoir permeability between the hydraulic fractures of the producer and injector) is carried out to further explore and quantify the potential of the shale gas reservoir performance with the CO 2 injection. The spacing between the hydraulic fractures of the producer and injector is 55 m. On the basis of scheme 2, three subschemes are established: scheme 2.1 increases the interfracture permeability to 10 times of the original permeability, scheme 2.2 increases to 100 times, and scheme 2.3 increases to 1000 times. Here, scheme 2 is set as a benchmark for the three subschemes.

Geofluids
With the outcomes of simulation (Figures 8-10), it is found that under the same gas injection volume as scheme 2, the total gas production of scheme 2.1 is 9:58 × 10 6 m 3 , which is increased by 0.21% compared with scheme 2. CH 4 production is 9:33 × 10 6 m 3 , increased by 0.21%; CO 2 production is 0:25 × 10 6 m 3 , increased by 1.08%. Scheme 2.2 total gas production is 9:71 × 10 6 m 3 , increased by 1.57%; the CH 4 production is 9:42 × 10 6 m 3 , increased by 1.18%. CO 2 production is 0:29 × 10 6 m 3 , increased by 16%. Scheme 2.3 total gas production is 10:25 × 10 6 m 3 , increased by 7.21%; CH 4 production is 9:75 × 10 6 m 3 , increased by 4.73%; CO 2 production is 0:50 × 10 6 m 3 , increased by 100.56%. Based on the above numerical experiments, in a word, the bigger 12 Geofluids the interfracture permeability is, the better the reservoir performance is. Due to the increase of interfracture permeability, the pressure transmits easier from the injector to the producer ( Figure 11). The efficiency of supplemental energy utilization is improved, and the average pressure of the whole reservoir declines faster, which is the main reason for the increase in CH 4 production. The finding indicates that energy utilization efficiency plays a key role in improving shale gas recovery by CO 2 injection. Otherwise, the supplemental energy will not be effectively harnessed to benefit the reservoir development. The pressure status of the grid blocks, which is highlighted in Figures 8-10, also indicates that the producer of the case with bigger interfracture permeability exhibits higher energy utilization efficiency in the drainage area. In other words, an effective pressure-driven system between the injector and producer is established and more supplemental energy is utilized to benefit the recovery process for the case with bigger interfracture permeability. Meanwhile, not only CH 4 production has increased, but also CO 2 production has increased significantly with the increase of interfracture permeability (Table 3), which is harmful to the sequestration of CO 2 in shale reservoirs. In order to meet the demand for both CO 2 sequestration and EGR, it is essential to establish effective communication between the injector and producer via fracking, which is beneficial to efficient energy transmission. Meanwhile, the communication between the injector and producer should be appropriate to prevent CO 2 from channeling, which is beneficial to the effective sequestration of CO 2 . The successful CO 2 interwell flooding strategy in shale reservoirs for the sequestration process sets higher requirements for on-site fracking operations. A moderate secondary stimulation zone needs to be formed between the primary hydraulic fractures of the injector and producer to facilitate the efficient energy transfer between interwell as well as to prevent CO 2 from channeling, such as scheme 2.3. Before the field pilot, it is necessary to take the efficiency of supplemental energy utilization, the CO 2 sequestration ratio, and the effect of injected CO 2 on the purity of produced methane into consideration to design an optimal execution plan.

Conclusions
In this paper, with the data from the New Albany Shale reservoir in the Illinois Basin, the CMG-GEM simulator is implemented to establish a numerical model to evaluate the feasibility of CO 2 sequestration in shale gas reservoirs with potential enhanced gas recovery (EGR   13 Geofluids shale reservoir model, different CO 2 injection schemes (continuous injection vs. pulse injection) for CO 2 sequestration in shale gas reservoirs are investigated. Meanwhile, a sensitivity analysis of the reservoir permeability between the hydraulic fractures of production and injection wells is conducted to quantify its influence on reservoir performance. Based on the above research, the following conclusions are obtained: (1) CO 2 sequestration in shale gas reservoirs is technically feasible. With appropriate well spacing and effective stimulation, the gas production is improved by the injection of carbon dioxide which also lowers the purity of produced natural gas (2) Due to the tight and unstimulated matrix between the hydraulic fractures of the injector and producer, the pressure transfer efficiency of both continuous injection and pulse injection is low, which leads to the inefficient development of shale gas. But the injected CO 2 is effectively sequestered. The tight interfracture formation acts as a double-edged sword here. The tight interfracture formation reduces the energy transfer efficiency, which dominates the success of the EGR. Meanwhile, the tight interfracture formation effectively prevents CO 2 from channeling, which benefits CO 2 sequestration in shale. In the field application, it is necessary to take the efficiency of supplemental energy utilization, the CO 2 sequestration ratio and the effect of injected CO 2 on the purity of produced methane into consideration to design an optimal execution plan (3) Based on the sensitivity analysis of interfracture reservoir permeability conducted to quantify its influence on reservoir performance, the success of the EGR is determined by the energy transfer efficiency between the injector and producer. Meanwhile, with the increase of the interfracture reservoir permeability, the CO 2 sequestration ratio decreases. In order to meet the demand for CO 2 sequestration in shale gas reservoirs with EGR, advanced and effective fracking is essential, which means that a moderate secondary stimulation zone needs to be formed between the primary hydraulic fractures of injection and production wells to facilitate efficient energy transfer between interwell as well as to prevent CO 2 from channeling, such as scheme 2.3

Data Availability
Data is available upon request.

Conflicts of Interest
The authors declare no conflict of interest.