Numerical modeling of the dynamic variation in multiphase CH₄ during CO₂ enhanced gas recovery from depleted shale reservoirs

Abstract

Regarding CO₂ enhanced shale gas recovery, this work focuses on changes in the multiphase (free/adsorbed) CH₄ in the process of CO₂ enhanced shale gas recovery, by utilizing a rigorous numerical model with real geological parameters. This work studies nine injection well (IW) and CH₄ production well (PW) combinations of CO₂ to determine the influence of IW and PW locations on the dynamic interaction of multiphase CH₄ during 10000 d of CO₂ injection. The results indicate that the content of both the adsorbed CH₄ and free CH₄ is strongly variable before (and during) the CO₂-CH₄ displacement. In addition, during the simulation process, the proportion of the adsorbed CH₄ among all extracted CH₄ phases dynamically increases first and then tends to stabilize at 70%–80%. Moreover, the IW-PWs combinations significantly affect the outcomes of CO₂ enhanced shale gas recovery – for both the proportion of adsorbed/free CH₄ and the recovery efficiency. A longer IW-PW distance enables more adsorbed CH₄ to be recovered but results in a lower efficiency of shale gas recovery. Basically, a shorter IW-PWs distance helps recover CH₄ via CO₂ injection if the IW targets the bottom layer of the Wufeng-Longmaxi shale formation. This numerical work expands the knowledge of CO₂ enhanced gas recovery from depleted shale reservoirs.

Graphical abstract

Keywords

CO₂-CH₄ displacement / free gas / Longmaxi shale / CH₄ desorption / numerical simulation

1 Introduction

As a well-known clean energy source, so-called shale gas with CH₄ as the main component is lauded as having a significant supply to meet the increasing energy consumption worldwide (Hou et al., 2018; Cui et al., 2019; Li et al., 2019a; Li et al., 2019b; Gao et al., 2020), making the improvement in gas recovery from tight shale reservoirs a popular topic of study in the petroleum and natural gas exploration domains worldwide (Liu et al., 2016; Zhang et al., 2017; Li et al., 2018; Guo et al., 2020b; Li et al., 2020a; Mazumder et al., 2020). Against this background, the injection of CO₂ into shale reservoirs has received widespread public attention because it enables simultaneous CO₂ sequestration and enhanced gas recovery (CS-EGR) as a result of CO₂/CH₄ interactions in the shale (Chi et al., 2017; Fan et al., 2018; Zhou et al., 2019a; Klewiah et al., 2020; Li et al., 2020b). Therefore, the shale-based CS-EGR technique has been researched in numerous studies, which clarify that it is feasible in theory but difficult in practice because this technique includes strategy optimization, site characterization/monitoring, hazard assessment and management, etc. (Abidoye et al., 2015; Pan et al., 2018; Luo et al., 2019; Zhang et al., 2019; Chen et al., 2020; Iddphonce et al., 2020; Wang et al., 2020). In other words, CS-EGR in shale has been popularized, but full deployment worldwide requires more research to deepen the knowledge of this promising technique.

Basically, in the investigation of shale-based CS-EGR operation, the approaches adopted in the existing literature vary, but some sound achievements are as follows (Iddphonce et al., 2020; Klewiah et al., 2020; Li et al., 2020c). Experiments with low-field NMR scanning have exhibited the desorption of adsorbed CH₄ and CO₂/CH₄ competitive adsorption behavior after CO₂ is injected into shale samples (Liu et al., 2017b; Liu et al., 2019). Another self-designed experiment has confirmed that the preferential adsorption ratio of CO₂ over CH₄ varies between 1.66 and 8.32 (Zhou et al., 2019c). In addition, numerical simulations have been used to accurately characterize the effectiveness of shale-based CS-EGR. For example, a numerical approach was used to show that the proportion of CO₂ in a CO₂-N₂ mixture significantly affects the improved recovery efficiency of shale gas (Li and Elsworth, 2019). For example, the results of grand canonical Monte Carlo simulations have indicated that a higher water content is helpful in trapping more CO₂ during the shale-based CS-EGR process (Zhou et al., 2019b). The majority of interest in CS-EGR in shale has focused on mainly the CH₄ in the adsorbed phase (Wang et al., 2018; Zhang et al., 2020; Liu et al., 2021a; Liu et al., 2021b), causing research on free CH₄ during the CS-EGR process to be lacking to some extent. However, free CH₄ inevitably emerges in the inner spaces (e.g., pores and fractures) of shale during CS-EGR (Fathi and Akkutlu, 2014). Therefore, compared with the considerable attention paid to adsorbed CH₄, the relatively insubstantial attention paid to free CH₄ tends to limit the awareness of the coordinated behavior of CH₄ in multiple phases (adsorbed/free) and thus the comprehensive understanding of shale-based CS-EGR.

Recently, in some investigations on CO₂-enhanced shale gas recovery, how free CH₄ interacts with adsorbed CH₄ in shale has been studied. For example, the proportion of free CH₄ to adsorbed CH₄ during shale-based CS-EGR was explored and the influence of free CH₄ on CO₂ sequestration in shale was investigated by performing numerical work (Tao et al., 2014; Mohagheghian et al., 2019). Nevertheless, existing numerical simulations relating to the interaction between free CH₄ and adsorbed CH₄ are mainly organized using a plain mathematical model with only a single layer or/and isotropic conditions considered, in the process of CO₂ enhanced shale gas recovery. This phenomenon does not conform with the fact that field shale reservoirs is always characterized as multilayer and heterogeneous (Su et al., 2007; Pan et al., 2015; Liu et al., 2017a). Under such circumstances, more attention should be paid to the dynamic interaction of free CH₄ and adsorbed CH₄ in shale after CO₂ is injected, helping promote the development of this CS-EGR technique, which is what this work focuses on.

This work develops an advanced thermal-hydraulic-mechanical (THM) coupled model, where the rock deformation, competitive binary gas sorption, gas and water two-phase flow and thermal expansion in the dual-porosity system are comprehensively considered, using the modeling software COMSOL Multiphysics. This novel model is used to reflect the vertical heterogeneity of the porosity, permeability, Poisson’s ratio, Young’s modulus and Langmuir volume/pressure (CH₄ and CO₂) of a depleted shale formation from the Upper Ordovician Wufeng formation to the lower Silurian Longmaxi formation (WL, for short), ensuring that the model matches the real reservoir conditions as accurately as possible, to focus on the scenarios of most interest. Herein, the WL shale is chosen as the representative example case because of its promising potential for shale gas exploration and exploitation in the Sichuan Basin (Liu et al., 2017a; Liu et al., 2017c; Shan et al., 2017). Accordingly, the dynamic variation in and the outputs of multiphase CH₄ (free/adsorbed) during the CS-EGR process are investigated under nine different arrangements of CO₂ injection well (IW) and CH₄ production well (PW). Furthermore, the implications of this work for CO₂-enhanced shale gas recovery are also suggested. Through these efforts, this work offers a new modeling concept and a novel perspective on CS-EGR in shale and further deepens the understanding of CS-EGR in shale reservoirs.

2 Numerical model description

In this work, the simulated reservoir represents the WL shale as observed from a shale gas exploration well (Well-WQ2) located in north-eastern Sichuan Basin at the intersection of Sichuan Province, Shaanxi Province, Hubei Province and Chongqing municipality (Fig. 1(a)). In this region, the geological background is complicated, as described in previous works in detail (Liu et al., 2016; Liu et al., 2017a; Shan et al., 2017; Zhao et al., 2020). Basically, the thickness of the WL shales is approximately 100 m (in which the Wufeng Formation occupies ~10 m), and the WL shales are located in the burial depth range of 1200–1300 m, according to the drilling and logging information. Accordingly, the depleted pressure is estimated to be approximately 1.7 MPa (Zhang et al., 2015). Moreover, strong vertical heterogeneity of many of the WL shale reservoir parameters (Liu et al., 2017a; Zhang et al., 2018), including permeability, porosity, Poisson’s ratio, Young’s modulus and Langmuir volume/pressure, is observed (Table 1).

Referring to a typical IW-PW arrangement (Li and Elsworth, 2019), the simulation model contains one IW in the center and four PWs in the corners, which is an axisymmetric pattern in geometric space (Fig. 1(b)). A two-dimensional section is selected for this simulation work, in which the anisotropic permeability is considered, and each layer is embedded with the particular parameters shown in Table 1, ensuring that the simulated reservoir reflects the real WL shales as accurately as possible (Fig. 1(b)). To avoid the potential influence from the right boundary (Li et al., 2017), for example, the influence of formation overpressure, the width of the model reservoir is set to 300 m. Herein, the 100 m × 300 m rectangular area is divided into 5085 elements, where the gas pressure at the right boundary is set as the initial gas pressure and the rest of the boundaries do not allow flux. For this simulation work, the temperature is 305 K, the initial pressure of the CO₂ and CH₄ is 0.3 MPa and 1.7 MPa (same as the depleted pressure), respectively. During CS-EGR, a constant pressure of 7 MPa and 0.1 MPa is applied to the IW and the PWs, respectively. Other basic parameters (Table 2) adopted in this numerical work refer to existing work (Sun et al., 2013; Fathi and Akkutlu, 2014; Li and Elsworth, 2015; Fan et al., 2019a; Fan et al., 2019b; Li and Elsworth, 2019; Guo et al., 2020a; Liu et al., 2020; Zhao et al., 2020).

Based on the simulated reservoir, nine IW-PW combinations are studied to determine how the relative locations of IW and PWs affect the dynamic interaction of multiphase CH₄ during CO₂-enhanced gas recovery from depleted WL shales. When the IW targets the lower layer and has a horizontal distance from the PWs of 50 m, this case is labeled L50. Similarly, the remaining cases are L100, L150, M50, M100, M150, H50, H100 and H150 (Fig. 1(b)). Herein, the operation conditions for all cases are identical, except the IW-PWs locations. In this work, the candidate IW is only set in layer 3, the transition zone of layers 5-6, and layer 8 to represent the upper, middle and lower shale reservoirs, respectively. In addition, for each simulated case, the CS-EGR operation is independent, and the continuous run time is set to 10000 days (approximately 30 years).

3 THM coupling process and governing equations

Referring to previous achievements based on numerical works (Vega et al., 2014; Hu et al., 2019; Li and Elsworth, 2019), the simulated shale reservoir is treated as a dual-porosity medium (Fig. 1(c)). In this dual-porosity medium, the injection of CO₂ tends to activate a very complicated process due to a series of feedbacks controlled by the THM coupling effect, as exhibited by Fig. 1(d). Basically, the pivotal equations involved in the numerical modeling of this work refer to previous achievements (Fan et al., 2019a; Fan et al., 2019b; Li and Elsworth, 2019; Zhang et al., 2020; Zhao et al., 2020), mainly describing the following seven aspects

I. Mechanical field

G{u}_{k,ll}+{\displaystyle \frac{G}{\mathrm{1}-\mathrm{2}\nu }}{u}_{l,lk}-({\alpha }_{\mathrm{m}}{p}_{\mathrm{m},k}+{\alpha }_{\mathrm{f}}{p}_{\mathrm{f},k})-K{\alpha }_{\mathrm{T}}{T}_k

(1)-K{\left({\displaystyle \frac{{ϵ}_{L\mathrm{1}}{b}_{\mathrm{1}}{p}_{\mathrm{m}\mathrm{g}\mathrm{1}}+{ϵ}_{L\mathrm{2}}{b}_{\mathrm{2}}{p}_{\mathrm{m}\mathrm{g}\mathrm{2}}}{\mathrm{1}+b_{\mathrm{1}}{p}_{\mathrm{m}\mathrm{g}\mathrm{1}}+b_{\mathrm{2}}{p}_{\mathrm{m}\mathrm{g}\mathrm{2}}}}\right)}_{,k}+f_k=\mathrm{0},

where the subscript i represents the gas component (1 for CH₄ and 2 for CO₂); u_k is the deformation in the k direction; G = E/(2+ 2ν) is the shear modulus, Pa; K = E/3(1-2ν) is the bulk modulus, Pa; E is Young’s modulus, Pa; ν is Poisson’s ratio; α_m and α_f are the Biot effective stress coefficients of the matrix and fracture, respectively; p_m is the total gas pressure in the matrix, Pa; p_mgi is the gas pressure in the pore system, Pa; p_f is the fluid stress in the fracture, Pa; α_T is the thermal expansion coefficient, 1/K; T is the temperature of gas, K; b_i = 1/P_Li; P_Li is the Langmuir pressure constant, Pa; ε_L1 and ε_L2 are the Langmuir strain coefficients of CH₄ and CO₂; and f_k is the volume force component, N.

II. Hydraulic field

{\displaystyle \frac{\partial }{\partial t}}({\phi }_{\mathrm{m}}{\displaystyle \frac{M_{\mathrm{g}i}}{RT}}{p}_{\mathrm{m}\mathrm{g}i}+{\displaystyle \frac{V_{Li}{b}_i{p}_{\mathrm{m}\mathrm{g}i}}{\mathrm{1}+b_{\mathrm{1}}{p}_{\mathrm{m}\mathrm{g}\mathrm{1}}+b_{\mathrm{2}}{p}_{\mathrm{m}\mathrm{g}\mathrm{2}}}}{\rho }_c{\rho }_{\mathrm{g}\mathrm{s}i}

\exp (-{\displaystyle \frac{c_{\mathrm{1}}}{\mathrm{1}+c_{\mathrm{2}}{p}_{\mathrm{m}}}}(T-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\left)\right)+\nabla ·(-D_i{M}_{\mathrm{g}i}\nabla {\phi }_{\mathrm{m}}{\displaystyle \frac{p_{\mathrm{m}\mathrm{g}i}}{RT}})

(2)=-{\displaystyle \frac{\mathrm{1}}{{\tau }_i}}{\displaystyle \frac{M_{\mathrm{g}i}}{RT}}(p_{\mathrm{m}\mathrm{g}i}-p_{\mathrm{f}\mathrm{g}i}),

where M_gi is the gas molar mass, g/mol; ϕ_m is the porosity of the shale matrix; R is the gas molar constant, J/(mol·K); ρ_c is the density of the shale skeleton, kg/m³; ρ_gsi is the gas density under standard conditions, kg/m³; p_fgi is the gas pressure in the fracture system, Pa; V_Li is the Langmuir volume constant, m³/kg; c₁ and c₂ are the thermal coefficients of gas sorption; D_i is the diffusion coefficient, m²/s; T_ref is the reference temperature, K; and τ_i is the desorption time, s.

{\displaystyle \frac{\partial \left(s_{\mathrm{g}}{\phi }_{\mathrm{f}}{\rho }_{\mathrm{f}\mathrm{g}i}\right)}{\partial t}}+\nabla ·(-{\rho }_{\mathrm{f}\mathrm{g}i}{\displaystyle \frac{k{k}_{\mathrm{r}g}}{{\mu }_{\mathrm{g}i}}}(\mathrm{1}+{\displaystyle \frac{b_k}{p_{\mathrm{f}\mathrm{g}i}}})\nabla p_{\mathrm{f}\mathrm{g}i}

+{\displaystyle \frac{\partial \left(s_{\mathrm{w}}{\phi }_{\mathrm{f}}{H}_{\mathrm{g}i}{\rho }_{f\mathrm{g}i}\right)}{\partial t}}+\nabla ·(-H_{\mathrm{g}i}{\rho }_{\mathrm{f}\mathrm{g}i}{\displaystyle \frac{k{k}_{\mathrm{r}\mathrm{w}}}{{\mu }_{\mathrm{w}}}}\nabla p_{\mathrm{f}\mathrm{w}})

(3)={\displaystyle \frac{\mathrm{1}}{{\tau }_i}}{\displaystyle \frac{M_{\mathrm{g}i}}{RT}}(p_{\mathrm{m}\mathrm{g}i}-p_{\mathrm{f}\mathrm{g}i}),

{\displaystyle \frac{\partial \left(s_{\mathrm{w}}{\phi }_{\mathrm{f}}{\rho }_{\mathrm{w}}\right)}{\partial t}}+\nabla ·(-{\rho }_{\mathrm{w}}{\displaystyle \frac{k{k}_{\mathrm{r}\mathrm{w}}}{{\mu }_{\mathrm{w}}}}\nabla p_{\mathrm{f}\mathrm{w}})

+{\displaystyle \frac{\partial }{\partial t}}\left(s_{\mathrm{g}}{\phi }_{\mathrm{f}}{\rho }_{\mathrm{f}\mathrm{v}\mathrm{0}}exp\right({\displaystyle \frac{p_{\mathrm{c}\mathrm{g}\mathrm{w}}}{{\rho }_{\mathrm{w}}{R}_{\mathrm{v}}T}}\left)\right) +

\nabla ·(-{\rho }_{\mathrm{f}\mathrm{v}\mathrm{0}}exp({\displaystyle \frac{p_{\mathrm{c}\mathrm{g}\mathrm{w}}}{{\rho }_{\mathrm{w}}{R}_{\mathrm{v}}T}}){\displaystyle {\sum }_{i=\mathrm{1}}^{\mathrm{2}}}

(4){\displaystyle \frac{k{k}_{\mathrm{r}\mathrm{g}}}{{\mu }_{\mathrm{g}i}}}(\mathrm{1}+{\displaystyle \frac{b_k}{p_{\mathrm{f}\mathrm{g}i}}})\nabla p_{\mathrm{f}\mathrm{g}i}=\mathrm{0},

where ϕ_f is the fracture porosity; s_g is the gas saturation; s_w = 1−s_g is the water saturation; ρ_fgi is the gas density, kg/m³; H_gi is Henry’s coefficient of gas; k is the absolute permeability, m²; k_rg and k_rw are the relative permeabilities of the gas and water; m_gi and m_w are the dynamic viscosities of the gas and water, Pa·s; b_k is the Klinkenberg factor, Pa; p_fw = p_fg−p_cgw is the water pressure in the fracture, Pa; p_fg is the total gas pressure in the fracture, Pa; p_cgw is the capillary pressure, Pa; ρ_w is the density of water, kg/m³; ρ_fv0 is the density of saturated vapor, kg/m³; R_v is the latent heat of vapor, J/(K·kg); s_gr is the residual gas saturation; and s_wr is the irreducible water saturation.

III. Thermal field

{\displaystyle \frac{\partial }{\partial t}}\left({\left(\rho C_p\right)}_{\mathrm{e}\mathrm{f}\mathrm{f}}T\right)+{\eta }_{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla T-\nabla ·({\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla T)+K{\alpha }_{\mathrm{T}}T{\displaystyle \frac{\partial {ϵ}_{\mathrm{V}}}{\partial t}}

(5)+{\displaystyle {\sum }_{i=\mathrm{1}}^{\mathrm{2}}{q}_{\mathrm{s}\mathrm{t}i}{\displaystyle \frac{{\rho }_{\mathrm{s}}{\rho }_{\mathrm{g}\mathrm{s}i}}{M_{\mathrm{g}i}}}{\displaystyle \frac{\partial V_{\mathrm{s}\mathrm{g}i}}{\partial t}}}=\mathrm{0},

where (ρC_p)_eff is the effective heat capacity, J/(m³·K); η_eff is the effective heat convection transfer coefficient, J/(m²·s); l_eff is the effective thermal conductivity, W/(m·K); q_sti is the isosteric heat of gas adsorption, kJ/mol; ε_V is the volume strain; and V_sgi is the absorbed gas content, m³/kg.

IV. Matrix porosity

(6){\phi }_{\mathrm{m}}\mathrm{=}{\phi }_{\mathrm{m}\mathrm{0}}+{\displaystyle \frac{({\alpha }_{\mathrm{m}}-{\phi }_{\mathrm{m}\mathrm{0}})({ϵ}_{\mathrm{e}}-{ϵ}_{\mathrm{e}\mathrm{0}})}{\mathrm{1}+{ϵ}_{\mathrm{e}}}},

where ε_e€=€ε_v€+€p_m/K_s−α_TT−ε_s; and the subscript 0 represents the initial value of the parameters.

V. Fracture porosity

(7){\phi }_{\mathrm{f}}\mathrm{=}{\phi }_{\mathrm{f}\mathrm{0}}(\mathrm{1}+{\displaystyle \frac{{\displaystyle {\sum }_{j=\mathrm{1}}^N\mathrm{\Delta }{b}_j}}{{\displaystyle {\sum }_{j=\mathrm{1}}^N{b}_{j\mathrm{0}}}}})\mathrm{=}{\phi }_{\mathrm{f}\mathrm{0}}(\mathrm{1}+{\displaystyle \frac{{\displaystyle {\sum }_{j=\mathrm{1}}^N(\mathrm{\Delta }{\sigma }_{\mathrm{n}j}-{\alpha }_{\mathrm{f}}\mathrm{\Delta }{p}_{\mathrm{f}})/\left(K_{\mathrm{n}j}{(\mathrm{1}-\frac{{\sigma }_{\mathrm{n}j}^{\prime }}{K_{\mathrm{n}j}\mathrm{\Delta }{v}_{\mathrm{m}\mathrm{a}\mathrm{x}}+{\sigma }_{\mathrm{n}j}^{\prime }})}^{-\mathrm{2}}\right)}}{{\displaystyle {\sum }_{j=\mathrm{1}}^N{b}_{j\mathrm{0}}}}}),

where the subscript j represents the direction; Ds_nj is the normal stress on the fracture, Pa; K_nj is the initial normal stiffness, Pa/m; s’_nj is the normal effective stress, Pa; Dv_max is the maximum fracture aperture as applied normal effective stress goes to infinity, m; b_j0 is the initial fracture width, m; and N represents the dimension – for two dimensions N = 2, and for three dimensions N = 3.

VI. Permeability evolution

(8)k_j\mathrm{=}{k}_{j\mathrm{0}}{(\mathrm{1}+{\displaystyle \frac{\mathrm{\Delta }{b}_j}{b_{j\mathrm{0}}}})}^{\mathrm{3}},

where Db_j is the change of the fracture width, m.

VII. Relative permeability

(9)k_{\mathrm{r}\mathrm{g}}=k_{\mathrm{r}\mathrm{g}\mathrm{0}}{(\mathrm{1}-({\displaystyle \frac{s_{\mathrm{w}}-s_{\mathrm{w}\mathrm{r}}}{\mathrm{1}-s_{\mathrm{g}\mathrm{r}}-s_{\mathrm{w}\mathrm{r}}}}\left)\right)}^{\mathrm{2}}(\mathrm{1}-{\left({\displaystyle \frac{s_{\mathrm{w}}-s_{\mathrm{w}\mathrm{r}}}{\mathrm{1}-s_{\mathrm{w}\mathrm{r}}}}\right)}^{\mathrm{2}}),

(10)k_{\mathrm{r}\mathrm{w}}=k_{\mathrm{r}\mathrm{w}\mathrm{0}}{\left({\displaystyle \frac{s_{\mathrm{w}}-s_{\mathrm{w}\mathrm{r}}}{\mathrm{1}-s_{\mathrm{w}\mathrm{r}}}}\right)}^{\mathrm{4}},

where k_rg0 and k_rw0 are the endpoint relative permeabilities of the gas and water, respectively; s_gr is the residual gas saturation; and s_wr is the irreducible water saturation.

4 Results and discussion

Depending on the rigorous model, this work outputs the multiphase CH₄ content in the depleted WL shales just before external CO₂ is injected into the reservoir. Then, after CO₂ is injected, how the free CH₄ and adsorbed CH₄ dynamically vary during the CS-EGR process is explored. Following this issue, the influence of relative IW-PWs locations on the resulting interaction of multiphase CH₄ is also clarified. On this basis, this work further discusses the implications for CO₂-enhanced gas recovery in depleted WL shales. In addition, it should be noted that the thickness of the simulated reservoir is set to 1 m when the free/adsorbed CH₄ content is calculated.

4.1 Original characteristics of multiphase CH₄ before CO₂ involvement

Controlled by the heterogeneous parameters shown in Table 1, the residual content of CH₄ in each layer of depleted WL shales differs (Fig. 2), where the content of the adsorbed CH₄ is mainly determined by the reservoir pressure and Langmuir volume/pressure, while that of the free CH₄ is calculated based on the reservoir pressure and fracture/matrix porosity. With regard to the total CH₄ content in the depleted WL shales, the greatest content of residual CH₄ occurs in the layer at a depth of ~1275 m (layer 8) and is ~4604 kg, which is roughly two times the residual content (~2334 kg) of CH₄ in the layer at ~1225 m (layer 3). By contrast, the adsorbed CH₄ tends to occupy a greater proportion than the free CH₄ (excluding the situation in layer 1), where the average percentages of adsorbed CH₄ and free CH₄ are ~71% and ~29%, respectively (Fig. 2). Nevertheless, the proportion of multiphase CH₄ is variable among the layers with different burial depths, and the deeper layer generally has a higher proportion of adsorbed CH₄ in the total residual CH₄ (Fig. 2).

In particular, the situation of multiphase CH₄ in the depleted WL shales exhibits the following rules: 1) the Langmuir pressure/volume determines the content of adsorbed CH₄, and 2) the matrix/fracture porosity controls the content of free CH₄, under the condition that all layers in the vertical profile exhibit the same depleted pressure, before CO₂ injection. Therefore, in the depleted WL shales, the original contents of the CH₄ phases in the multiphase CH₄ are independent; in other words, the original content and proportion of adsorbed CH₄ shows no relationship with those of free CH₄ (Fig. 2). It seems almost inevitable that the vertical variability in the content/proportion of multiphase CH₄ among the layers, resulting from the CO₂-CH₄ displacement that occurs when the injected CO₂ enters the layers at different times, tends to increase the complexity of the dynamic interaction between adsorbed CH₄ and free CH₄.

4.2 Dynamic variation in multiphase CH₄ during the CS-EGR process

Current knowledge supports that the CH₄ content in shale reservoirs decreases after CO₂ is injected (Sun et al., 2013; Iddphonce et al., 2020; Klewiah et al., 2020; Zhao et al., 2020). Herein, the dynamic decrease in the free CH₄ and adsorbed CH₄ during 10000 d of CS-EGR is obtained from this simulation work. As shown in Fig. 3, the decrease in multiphase CH₄ is activated by CO₂ injection and increases with time in the process of CS-EGR in WL shales. Meanwhile, the decrease in adsorbed CH₄ is always greater than that in free CH₄ during the CO₂-CH₄ displacement for all IW-PW combinations. Nevertheless, in the whole WL shale, the decreases in adsorbed CH₄ and free CH₄ differ due to the different IW-PWs arrangements. For example, case H50 (Fig. 3(a)) allows more CH₄ to be driven out from the WL shales than case M150 (Fig. 3(b)).

During the decrease in the free CH₄ content in the WL shales after CO₂ injection, the reduction in multiphase CH₄ varies with time. Basically, this variation is j a two-step process, with a dynamic stage (Stage I) and stable stage (Stage II), illustrated by the two representative cases (Fig. 4). In the early phase, the sudden injection of CO₂ forces the CH₄ that is originally in the adsorbed phase to convert to free CH₄, and the CH₄ released from the PWs is limited due to the short extraction time, together enabling the relative proportion of the decrease in adsorbed CH₄ to increase and that of free CH₄ to decrease, i.e., Stage I. Furthermore, the CO₂-CH₄ displacement occurs in a larger area with continued CO₂ migration, exhibiting a dynamic change with the continual decrease in free CH₄ at the PWs, resulting in the decrease in a stable proportion of free/adsorbed CH₄ among the total CH₄ during the continuous CO₂ injection into the WL shales, i.e., Stage II. By contrast, Stage I tends to cover a longer period with a longer IW-PWs distance (Fig. 4). In this work, according to statistics, the resulting proportions of adsorbed CH₄ among all the extracted CH₄ under different IW-PWs combinations are similar (that is, 70%–80%).

With the decrease in CH₄ caused by CO₂ injection, the proportion of the residual multiphase CH₄ in the WL shales varies dynamically. Figure 5 shows that the relative IW-PWs location significantly affects the variation in the proportion of free CH₄ in the modeling reservoir. For the “50” series (H50, M50 and L50), the proportion of free CH₄ first slowly increases and then decreases linearly with CO₂ injection because the insufficient CO₂-CH₄ displacement due to limited CO₂ migration (a short IW-PW distance) is unable to compensate for the free CH₄ released from the PWs. With regard to the “150” series (H150, M150 and L150), CO₂-CH₄ displacement occurs in a large area (a long IW-PW distance) and thus generates massive free CH₄ that is converted from the adsorbed CH₄, overshadowing the free CH₄ extraction from the PWs, ultimately enabling a roughly linear increase in the proportion of residual free CH₄ in the WL shales. For the “100” series (H100, M100 and L100), the results are between those of the “50” series and the “150” series, suggesting a Langmuir-like variation in the proportion of residual free CH_4, which first increases and then tends to be stable (Fig. 5).

In addition, adsorbed CH₄ exists on the surface of the shale skeleton, while free CH₄ appears in both the fractures and the matrix pores in this dual-porosity system of modeling shale reservoirs (Fan et al., 2019a; Li and Elsworth, 2019). Accordingly, it is meaningful to discuss the dynamic variation in free CH₄ in fractures and matrix pores, helping enhance the understanding of the change in free CH₄ during the shale-based CS-EGR process. Generally, CO₂ injection enables a roughly linear decline in the free CH₄ content in both fractures and matrix pores (Fig. 6). Therein, the content of free CH₄ in the matrix pores is always greater than that in the fractures, controlled by the fact that j_f is less than j_m (Table 1). However, Fig. 6 also states that the ratio of the free CH₄ in the fractures to that in the matrix pores differs during CS-EGR operation under variable IW-PWs combinations. According to Fig. 7, the variation range of the relative content of CH₄ in the fractures to that in the matrix pores is not significant and is only between 12.1% and 12.6%. To describe this phenomenon in detail, the variation in this ratio (see Fig. 7) can be divided into three types, corresponding to the “50” series, “100” series and “150” series, where the change law is similar to that shown in Fig. 5.

4.3 Outputs of multiphase CH₄ after 10000 d of CO₂ injection

The vertical nonuniform distribution of multiphase CH₄ in the modeling reservoir, accompanied by the complex/dynamic variation in multiphase CH₄ during CS-EGR operation, varies the outputs of multiphase CH₄ after 10000 d of CO₂ injection under nine IW-PWs combinations. For example, the residual adsorbed CH₄, the residual free CH₄ and their ratio (relative content) resulting from case H50 (Fig. 8(a)) are totally different from those resulting from case M150 (Fig. 8(b)). Nevertheless, regardless of where the IW and PWs are located, the adsorbed CH₄ content is always higher than the free CH₄ content (Fig. 8) throughout the process of CO₂-enhanced shale gas recovery from the WL shales. This phenomenon is determined by the reservoir parameters, even though the local pressure of the modeling reservoir is variable after CO₂ injection (Zhao et al., 2020).

In fact, the residual CH₄ content in the whole reservoir (or in each layer) is the difference between the original CH₄ content and the decreased CH₄ content in the whole reservoir (or in each layer). Thus, in this work, the final decrease in multiphase CH₄ is introduced to analyze the CS-EGR outcomes in the whole modeling reservoir or in each independent layer. The resulting decrease in free/adsorbed CH₄ under different IW-PW combinations differs (Fig. 9), revealing the significance of IW-PWs location in shale-based CS-EGR operation. Briefly, the vertical variability in the decrease in free/adsorbed CH₄ is strong, but a gradual increasing tendency roughly exists with burial depth (Fig. 9). For all the IW-PWs combinations, the bottom layer, with the most favorable sections for in situ shale gas development (Liu et al., 2017a), is the main contributors to the CH₄ recovery during the CS-EGR progress. According to Fig. 9, for every layer, a shorter IW-PW distance enables a greater decrease in CH₄ in the WL shales, in which the proportion of adsorbed CH₄ is lower, and vice versa. In addition, when the PWs are fixed, it seems that the IW location only slightly affects the decrease in CH₄ and the corresponding ratio of the adsorbed/free CH₄.

4.4 Implications for CO₂ enhanced shale gas recovery

From the perspective of the whole modeling reservoir, the injection of CO₂ partly recovers the residual CH₄ from the depleted WL shales (Fig. 10). When the IW location is fixed, the “150” series cases (H150, M150 and L150) have more CH₄ remaining in the reservoirs, lowering the corresponding proportion of adsorbed CH₄ relative to all the residual CH₄ in the WL shales (Fig. 10). That is, the “50” series generates a greater content of recovered CH₄ (only cases H50, M50 and L50 have CH₄ recoveries of more than 5000 kg), while the “150” series is not beneficial for CO₂ shale gas recovery (cases H150, M150 and L150 force only ~4000 kg CH₄ to be recovered) after 10000 d of CS-EGR operation in the WL shales (Fig. 11(a)). Moreover, under different IW-PW combinations, although the total recovered CH₄ contents differ, the proportions of adsorbed CH₄ relative to all the recovered CH₄ are close and vary in the range of 75%–78% (Fig. 11(a)).

Comparing the residual CH₄ contents in the WL shales before and after 10000 d CO₂ injection, the shale gas recovery efficiency is obtained (Fig. 12). Influenced by the variable interaction of multiphase CH₄ after CO2 is injected into the modeling reservoir, the resulting shale gas recovery efficiency is different under different IW-PWs combination. Basically, the efficiency of CO₂-enhanced shale gas recovery varies from 11% to 16%, with an average of 13.3%. Herein, the greatest efficiency of 15.6% emerges for case L50, while the lowest (~11.7%) corresponds to the “150” series (H150, M150 and L150) (Fig. 12). Accordingly, it could be speculated that a shorter IW-PWs distance is needed to obtain a higher efficiency of CO₂-enhanced shale gas recovery from the WL shales if the IW targets the bottom layer. Moreover, this speculation might be suitable for all the WL shales in the Sichuan Basin, assuming that the WL shales in the north-eastern Sichuan Basin are representative (Liu et al., 2017a; Su et al., 2007; Shan et al., 2017).

5 Conclusions

According to a novel numerical model with THM coupling and the measured reservoir parameters of the depleted WL shales located in the north-eastern Sichuan Basin, the dynamic variation in free CH4 and adsorbed CH4 during/after 10000 d of CS-EGR deployment is investigated under nine combinations of IW and PWs. The main conclusions are as follows.

1. As revealed by the numerical simulation, the residual CH₄ content in the depleted WL shales varies in the vertical direction, and the adsorbed CH₄ tends to occupy a greater proportion than the free CH₄. The residual CH₄ content and the corresponding ratio of adsorbed/free CH₄ vary in a complicated manner after CO₂ is injected into the shale reservoir.

2. For the decreased CH₄, this variation in the proportion of the multiphase CH₄ contains two stages, a dynamic stage and stable stage, in which the resulting proportions of adsorbed CH₄ among all extracted CH₄ phases under different IW-PWs combinations are similar (that is, 70%–80%).

3. The relative locations of the IW and PWs significantly affect the predicted CS-EGR outcomes in the WL shales, where the dynamic interaction of adsorbed CH₄ and free CH₄ differs under different IW-PWs combinations, resulting in the variable efficiency of CO₂-enhanced shale gas recovery. Basically, a shorter IW-PWs distance seems to be helpful in increasing the CH₄ recovery efficiency if the IW targets the bottom layer of the WL shales.

References

[1] Abidoye L K, Khudaida K J, Das D B (2015). Geological carbon sequestration in the context of two-phase flow in porous media: a review. Crit Rev Environ Sci Technol, 45(11): 1105–1147

[2] Chen B, Harp D R, Lu Z, Pawar R J (2020). Reducing uncertainty in geologic CO₂ sequestration risk assessment by assimilating monitoring data. Int J Greenh Gas Control, 94: 102926

[3] Chi J, Ju B, Lyu G, Zhang X, Wang J (2017). A computational method of critical well spacing of CO₂ miscible and immiscible concurrent flooding. Pet Explor Dev, 44(5): 771–778

[4] Cui J W, Zhu R K, Mao Z G, Li S X (2019). Accumulation of unconventional petroleum resources and their coexistence characteristics in Chang7 shale formations of Ordos Basin in central China. Front Earth Sci, 13(3): 575–587

[5] Fan C, Elsworth D, Li S, Zhou L, Yang Z, Song Y (2019a). Thermo-hydro-mechanical-chemical couplings controlling CH₄ production and CO₂ sequestration in enhanced coalbed methane recovery. Energy, 173: 1054–1077

[6] Fan C, Luo M, Li S, Zhang H, Yang Z, Liu Z (2019b). A thermo-hydro-mechanical-chemical coupling model and its application in acid fracturing enhanced coalbed methane recovery simulation. Energies, 12(4): 626

[7] Fan Y, Deng C, Zhang X, Li F, Wang X, Qiao L (2018). Numerical study of CO₂-enhanced coalbed methane recovery. Int J Greenh Gas Control, 76: 12–23

[8] Fathi E, Akkutlu I Y (2014). Multi-component gas transport and adsorption effects during CO₂ injection and enhanced shale gas recovery. Int J Coal Geol, 123: 52–61

[9] Gao Z, Fan Y, Xuan Q, Zheng G (2020). A review of shale pore structure evolution characteristics with increasing thermal maturities. Advances in Geo-Energy Research, 4(3): 247–259

[10] Guo C, Li R, Sun J, Wang X, Liu H (2020a). A review of gas transport and adsorption mechanisms in two-component methane-carbon dioxide system. Int J Energy Res, 44(4): 2499–2516

[11] Guo C, Li R, Wang X, Liu H (2020b). Study on two component gas transport in nanopores for enhanced shale gas recovery by using carbon dioxide injection. Energies, 13(5): 1–21

[12] Hou H H, Shao L Y, Li Y H, Li Z, Zhang W L, Wen H J (2018). The pore structure and fractal characteristics of shales with low thermal maturity from the Yuqia Coalfield, northern Qaidam Basin, northwestern China. Front Earth Sci, 12(1): 148–159

[13] Hu S, Hu X, He L, Chen W (2019). A new material balance equation for dual-porosity media shale gas reservoir. Energy Procedia, 158: 5994–6002

[14] Iddphonce R, Wang J, Zhao L (2020). Review of CO₂ injection techniques for enhanced shale gas recovery: prospect and challenges. J Nat Gas Sci Eng, 77: 103240

[15] Klewiah I, Berawala D S, Alexander Walker H C, Andersen P Ø, Nadeau P H (2020). Review of experimental sorption studies of CO₂ and CH₄ in shales. J Nat Gas Sci Eng, 73: 103045

[16] Li D, He Y, Zhang H, Xu W, Jiang X (2017). A numerical study of the impurity effects on CO₂ geological storage in layered formation. Appl Energy, 199: 107–120

[17] Li J, Li B, Wang Z, Ren C, Yang K, Chen S (2020a). An anisotropic permeability model for shale gas recovery considering slippage effect and embedded proppants. Nat Resour Res, 29(5): 3319–3333

[18] Li J, Lu S, Cai J, Zhang P, Xue H, Zhao X (2018). Adsorbed and free oil in lacustrine nanoporous shale: a theoretical model and a case study. Energy Fuels, 32(12): 12247–12258

[19] Li X, Elsworth D (2015). Geomechanics of CO₂ enhanced shale gas recovery. J Nat Gas Sci Eng, 26: 1607–1619

[20] Li Y, Wang Y B, Wang J, Pan Z J (2020b). Variation in permeability during CO₂-CH₄ displacement in coal seams: part 1-experimental insights. Fuel, 263: 116666

[21] Li Y, Wang Z S, Pan Z J, Niu X L, Yu Y, Meng S Z (2019a). Pore structure and its fractal dimensions of transitional shale: A cross-section from east margin of the Ordos Basin, China. Fuel, 241: 417–431

[22] Li Y, Yang J H, Pan Z J, Meng S Z, Wang K, Niu X L (2019b). Unconventional natural gas accumulations in stacked deposits: a discussion of Upper Paleozoic coal-bearing strata in the east margin of the Ordos Basin, China. Acta Geologica Sinica (English Edition), 93(1): 111–129

[23] Li Y, Yang J H, Pan Z J, Tong W S (2020c). Nanoscale pore structure and mechanical property analysis of coal: an insight combining AFM and SEM images. Fuel, 260: 116352

[24] Li Z, Elsworth D (2019). Controls of CO₂–N₂ gas flood ratios on enhanced shale gas recovery and ultimate CO₂ sequestration. J Petrol Sci Eng, 179: 1037–1045

[25] Liu J, Xie L, Elsworth D, Gan Q (2019). CO₂/CH₄ competitive adsorption in shale: implications for enhancement in gas production and reduction in carbon emissions. Environ Sci Technol, 53(15): 9328–9336

[26] Liu J, Xie L, He B, Gan Q, Zhao P (2021a). Influence of anisotropic and heterogeneous permeability coupled with in-situ stress on CO₂ sequestration with simultaneous enhanced gas recovery in shale: Quantitative modeling and case study. Int J Greenh Gas Control, 104: 103208

[27] Liu J, Xie L, He B, Zhao P, Ding H (2021b). Performance of free gases during the recovery enhancement of shale gas by CO₂ injection: a case study on the depleted Wufeng–Longmaxi shale in northeastern Sichuan Basin, China. Petrol Sci, doi:10.1007/s12182-020-00533-y

[28] Liu J, Yao Y, Elsworth D, Liu D, Cai Y, Dong L (2017a). Vertical heterogeneity of the shale reservoir in the Lower Silurian Longmaxi Formation: analogy between the southeastern and northeastern Sichuan Basin, SW China. Minerals (Basel), 7(8): 151

[29] Liu J, Yao Y, Liu D, Elsworth D (2017b). Experimental evaluation of CO₂ enhanced recovery of adsorbed-gas from shale. Int J Coal Geol, 179: 211–218

[30] Liu J, Yao Y, Liu D, Pan Z, Cai Y (2017c). Comparison of three key marine shale reservoirs in the southeastern margin of the Sichuan Basin, SW China. Minerals (Basel), 7(10): 179

[31] Liu J, Yao Y, Zhu Z, Cheng L, Wang G (2016). Experimental investigation of reservoir characteristics of the upper Ordovician Wufeng Formation shale in middle–upper Yangtze region, China. Energy Exploration & Exploitation, 34(4): 527–542

[32] Liu S, Fang H, Sang S, Ashutosh T, Wu J, Zhang S, Zhang B (2020). CO₂ injectability and CH₄ recovery of the engineering test in qinshui Basin, China based on numerical simulation. Int J Greenh Gas Control, 95: 102980

[33] Luo C J, Zhang D F, Lun Z M, Zhao C P, Wang H T, Pan Z J, Li Y H, Zhang J, Jia S Q (2019). Displacement behaviors of adsorbed coalbed methane on coals by injection of SO₂/CO₂ binary mixture. Fuel, 247: 356–367

[34] Mazumder M, Chen L, Xu Q (2020). Integrated ejector-based flare gas recovery and on-site desalination of produced water in shale gas production. Chem Eng Technol, 43(2): 200–210

[35] Mohagheghian E, Hassanzadeh H, Chen Z (2019). CO₂ sequestration coupled with enhanced gas recovery in shale gas reservoirs. Journal of CO₂ Utilization, 34: 646–655

[36] Pan Y, Hui D, Luo P, Zhang Y, Sun L, Wang K (2018). Experimental investigation of the geochemical interactions between supercritical co₂ and shale: implications for co₂ storage in gas-bearing shale formations. Energy Fuels, 32(2): 1963–1978

[37] Pan Z, Ma Y, Connell L D, Down D I, Camilleri M (2015). Measuring anisotropic permeability using a cubic shale sample in a triaxial cell. Journal of Natural Gas ence & Engineering, 26: 336–344

[38] Shan C, Zhang T, Wei Y, Zhang Z (2017). Shale gas reservoir characteristics of Ordovician-Silurian formations in the central Yangtze area, China. Front Earth Sci, 11(1): 184–201

[39] Su W, Li Z, Ettensohn F, Johnson M, Huff W, Wang W, Ma C, Li L, Zhang L, Zhao H (2007). Tectonic and eustatic control on the distribution of black-shale source beds in the Wufeng and Longmaxi formations (Ordovician-Silurian), South China. Front Earth Sci China, 1(4): 470–481

[40] Sun H, Yao J, Gao S H, Fan D Y, Wang C C, Sun Z X (2013). Numerical study of CO₂ enhanced natural gas recovery and sequestration in shale gas reservoirs. Int J Greenh Gas Control, 19: 406–419

[41] Tao Z, Bielicki J M, Clarens A F (2014). Physicochemical factors impacting CO₂ sequestration in depleted shale formations: The case of the Utica shale. Energy Procedia, 63: 5153–5163

[42] Vega B, Dutta A, Kovscek A R (2014). CT imaging of low-permeability, dual-porosity systems using high X-ray contrast gas. Transp Porous Media, 101(1): 81–97

[43] Wang M, Lun Z M, Zhao C P, Wang H T, Luo C J, Fu X X, Li C, Zhang D F (2020). Influences of primary moisture on methane adsorption within Lower Silurian Longmaxi Shales in the Sichuan Basin, China. Energy Fuels, 34(9): 10810–10824

[44] Wang T, Tian S, Li G, Sheng M, Ren W, Liu Q, Zhang S (2018). Molecular simulation of CO₂/CH₄ competitive adsorption on shale kerogen for co₂ sequestration and enhanced gas recovery. J Phys Chem C, 122(30): 17009–17018

[45] Zhang G, Zhou W, Ji S, Liu J, Zhang J, Yang H (2015). Experimental study on CO₂ replacement method used in shale gas exploration. Journal of Chengdu University of Technology, 42(3): 366–371

[46] Zhang H, Diao R, Mostofi M, Evans B (2020). Monte carlo simulation of the adsorption and displacement of CH₄ by CO₂ injection in shale organic carbon slit micropores for CO₂ enhanced shale gas recovery. Energy Fuels, 34(1): 150–163

[47] Zhang K, Jiang H, Qin G (2019). Utilization of zeolite as a potential multi-functional proppant for CO₂ enhanced shale gas recovery and CO₂ sequestration: A molecular simulation study on the competitive adsorption of CH₄ and CO₂ in zeolite and organic matter. Fuel, 249(1): 119–129

[48] Zhang L, Lu S, Jiang S, Xiao D, Chen L, Liu Y, Zhang Y, Li B, Gong C (2018). Effect of shale lithofacies on pore structure of the Wufeng-Longmaxi Shale in southeast Chongqing, China. Energy Fuels, 32(6): 6603–6618

[49] Zhang P, Lu S, Li J, Xue H, Li W, Zhang P (2017). Characterization of shale pore system: a case study of Paleogene Xin’gouzui Formation in the Jianghan basin, China. Mar Pet Geol, 79: 321–334

[50] Zhao P, Xie L, He B, Liu J (2020). Strategy optimization on industrial CO₂ sequestration in the depleted Wufeng-Longmaxi Formation shale in the northeastern Sichuan Basin, SW China: from the perspective of environment and energy. ACS Sustain Chem& Eng, 8(30): 11435–11445

[51] Zhou J, Hu N, Xian X, Zhou L, Tang J, Kang Y, Wang H (2019a). Supercritical CO₂ fracking for enhanced shale gas recovery and CO₂ sequestration: results, status and future challenges. Advances in Geo-Energy Research, 3(2): 207–224

[52] Zhou J, Jin Z, Luo K H (2019b). Effects of Moisture contents on shale gas recovery and CO₂ sequestration. Langmuir, 35(26): 8716–8725

[53] Zhou J, Liu M, Xian X, Jiang Y, Liu Q, Wang X (2019c). Measurements and modelling of CH₄ and CO₂ adsorption behaviors on shales: Implication for CO₂ enhanced shale gas recovery. Fuel, 251: 293–306

RIGHTS & PERMISSIONS

Higher Education Press