Performance of free gases during the recovery enhancement of shale gas by CO2 injection: a case study on the depleted Wufeng–Longmaxi shale in northeastern Sichuan Basin, China

In this work, a novel thermal–hydraulic–mechanical (THM) coupling model is developed, where the real geological parameters of the reservoir properties are embedded. Accordingly, nine schemes of CO2 injection well (IW) and CH4 production well (PW) are established, aiming to explore the behavior of free gases after CO2 is injected into the depleted Wufeng–Longmaxi shale. The results indicate the free CH4 or CO2 content in the shale fractures/matrix is invariably heterogeneous. The CO2 involvement facilitates the ratio of free CH4/CO2 in the matrix to that in the fractures declines and tends to be stable with time. Different combinations of IW–PWs induce a difference in the ratio of the free CH4 to the free CO2, in the ratio of the free CH4/CO2 in the matrix to that in the fractures, in the content of the recovered free CH4, and in the content of the trapped free CO2. Basically, when the IW locates at the bottom Wufeng–Longmaxi shale, a farther IW–PWs distance allows more CO2 in the free phase to be trapped; furthermore, no matter where the IW is, a shorter IW–PWs distance benefits by getting more CH4 in the free phase recovered from the depleted Wufeng–Longmaxi shale. Hopefully, this work is helpful in gaining knowledge about the shale-based CO2 injection technique.


Introduction
Current knowledge strongly supports CO 2 sequestration in geological formations (e.g., shale and coalbed) as a promising approach to respond to the issue of excessive anthropogenic CO 2 emissions into the atmosphere, which is treated as carbon-negative technology and thus has been receiving growing attention (Abidoye et al. 2015;Zhang et al. 2015;Wan and Liu 2018;Liu et al. 2020;Rani et al. 2020;Řimnáčová et al. 2020). Therein, injecting CO 2 into shale reservoirs is becoming increasingly attractive and has been lauded as a win-win solution for meeting geological CO 2 storage and simultaneous enhanced gas recovery (CS-EGR) from shales (Liu et al. 2017b(Liu et al. , 2019aKalra et al. 2018;Pan et al. 2018;Yang et al. 2019;Zhang et al. 2020). However, understanding of shale-based CS-EGR is limited, with only few pilot field studies undertaken, in which strategy optimization, site characterization/monitoring, and hazard assessment & management need to be taken into account in a comprehensive manner (Chi et al. 2017;Ajayi et al. 2019;Zhang et al. 2019;Chen et al. 2020;Iddphonce et al. 2020). In other words, despite being proven to be a feasible technique, CS-EGR in the shale is not yet mature enough for full national or worldwide deployment.
To promote the development of the shale-based CS-EGR technique, considerable work has been organized, in which significant achievements have been drawn (Iddphonce et al. 2020). For example, it has been proven that the CO 2 -CH 4 displacement behavior differs in the shales from different sedimentary environments (Du et al. 2019). Another example is that the CS-EGR efficiency is considerably affected by Edited by Yan-Hua Sun the CO 2 proportion in the injected CO 2 -N 2 mixture (Li and Elsworth 2019). In one other example, NMR-based investigation indicated that the properties (e.g., mineral composition) of shale reservoir can interfere with the CS-EGR outputs (Liu et al. 2019b). However, the examples mentioned above mainly pay attention to the gases (CO 2 /CH 4 ) in the adsorbed phase during the CO 2 -CH 4 interaction, in which the performance of the free gases (CO 2 /CH 4 ) tends to be ignored. Unfortunately, neglecting the free gases in the process of shale-based CS-EGR does not exist in isolation and even is a common phenomenon reported in many literatures (Liu et al. 2017b;Zhou et al. 2019;Řimnáčová et al. 2020).
Although the competitive adsorption of CO 2 /CH 4 is the primary element enabling CS-EGR in shales (Wang et al. 2018b;Liu et al. 2019a), the gases in free phase inevitably appear to the spaces of the pores and fractures in shale reservoirs during the CO 2 -CH 4 interaction (Fathi and Akkutlu, 2014). Basically, free gas is crucial for shale gas exploration/exploitation and even occupies over 50% of in situ gas content in shales (Liu et al. 2016b). To separately study the multiphase (free/adsorbed) gas in shales, experimental approaches have been conducted, and sound achievements have been made. For instance, isotopic geochemistry (Liu et al. 2016b) and low-field NMR (Yao et al. 2019) were successfully introduced to quantitatively identify the CH 4 in the free and adsorbed phase. Nevertheless, it was found that a single experimental methodology was usually performed at a restricted scale, and could hardly recognize the free CO 2 and free CH 4 from the CO 2 -CH 4 mixture at the same time during the dynamic CO 2 -CH 4 interaction in shales. Given this situation, physics-based numerical simulations can be employed to address concerns about the free gases during the shale-based CS-EGR process, in which multi-factors can be investigated in a clear and transparent manner.
Up to a point, numerical attempts have given rise to knowledge improvement on the performance of free gases during the CS-EGR operation in shales. There are two examples; one is that numerical work advocated that the portion of free gas is always less than that of adsorbed gas in the process of CO 2 -CH 4 interplay in shales, under variable reservoir permeability and temperature/pressure (Mohagheghian et al. 2019). Another modeling study has suggested that an increase in the ratio of the free CH 4 to adsorbed CH 4 would decrease the CO 2 sequestration potential of the Utica shale formation (Tao et al. 2014). Unfortunately, current simulations about free CH 4 /CO 2 in shales are usually run with a plain mathematical model with only a single layer or/and isotropic conditions being simulated, which is contrary to field fact that multilayer and heterogeneous characteristics exist extensively in shale reservoirs (Tang et al. 2016;Liu et al. 2017a). For this reason, limited existing numerical achievements can be directly generalized and carried out in a real shale formation that usually has a thickness of tens to several hundred meters, where vertical heterogeneity emerges in many reservoir parameters, such as the porosity/ permeability, the Langmuir volume, and Young's modulus (Tang et al. 2016;Liu et al. 2017a). Clearly, the development of a comprehensive model that contains as many realistic conditions as possible is, therefore, the need of the hour, aiming to enhance understanding of the free gases behavior during the CS-EGR process in shales.
In this work, because the Wufeng-Longmaxi Formation in the northeastern Sichuan Basin (SW China) has huge potential resources of shale gas (Liu et al. 2017a), it is selected for simulation where the CS-EGR technique has the opportunity to be implemented. Accordingly, a thermal-hydraulic-mechanical (THM) coupled model is developed using the COMSOL Multiphysics software, attempting to simulate the depleted condition of the Wufeng-Longmaxi shale. Herein, this integrated model not only contains the responses of the rock deformation, competitive CH 4 /CO 2 sorption, gas/water two-phase flow, and thermal expansion in the dual-porosity system but also involves the real heterogeneity of the reservoir properties of the Wufeng-Longmaxi shale. Using a mathematical approach, this work focuses on the dynamic interaction of free CO 2 and free CH 4 during the continuous CO 2 injection into the depleted Wufeng-Longmaxi shale, with a run time of roughly 30 years (10,000 d). Accordingly, in a quantitative manner, a series of simulations are set to clarify the anisotropic content variation of free gases (CO 2 and CH 4 ) in the modeling reservoir, under different relative locations of the CO 2 injection well (IW) and the CH 4 production well (PW). Furthermore, the strengths and weaknesses of each IW-PWs combo are also discussed, in terms of the different desired objectives. Hopefully, the peculiar perspectives from the gases in the free phase during the CO 2 -CH 4 interaction in shales will help to deepen awareness on the potential field deployment of the CS-EGR project in the shale reservoir.

Geological background
In this work, the case Wufeng-Longmaxi shale is from the Well-WQ2 that is a shale gas exploration well and is located in the northeastern Sichuan Basin, Southwest (SW) China (Fig. 1a). In this region, an arc tectonic belt system with a thrust belt trending to the southwest was formed by a series of tectonic activities, in which the Chengkou-Fangxian and Shashi Buried fault zones are quintessential examples (Li et al. 2013;Ji et al. 2015;Liu et al. 2017a). Figure 1b shows the continuous sedimentary succession in the study area from the late Sinian to the Cenozoic Era, except for the missing strata from the Devonian to the Carboniferous period (Meng and Zhang 2000;Shi et al. 2012;Li et al. 2013). Under this sedimentary background, previous studies have identified two conventional hydrocarbon source systems ( Fig. 1b): the bottom one (from the Cambrian Guojiaba Formation to the Silurian Cuijiagou Formation) and the upper one (from the Permian Qixia Formation to the Triassic Leikoupo Formation) (Li et al. 2013;Liu et al. 2017a). Therein, during the shale gas exploration and exploitation, the Wufeng-Longmaxi shale is usually treated as a close combination due to the conformable contact relationship between the upper Ordovician Wufeng Formation and the lower Silurian Longmaxi Formation (Zhu et al. 2016;Ye et al. 2017;Wang et al. 2018a).
The Wufeng-Longmaxi shale in the Well-WQ2 has a thickness of roughly 100 m with a buried depth of 1200-1300 m, where the bottom section (about 10 m) belongs to the Wufeng shale (Liu et al. 2016a(Liu et al. , 2017a. At such a buried depth, the depleted pressure of the Wufeng-Longmaxi shale is estimated at about 1.7 MPa (Zhang et al. 2015). Statistics suggest that the vertical heterogeneity is strong in the involved Wufeng-Longmaxi shale which also has anisotropic characteristics in the reservoir properties (e.g., porosity/permeability) (Liu et al. 2016a(Liu et al. , 2017a (Fig. 1c). In addition, the well logs of Well-WQ2 indicate the temperature of the Wufeng-Longmaxi shale is 36 °C (310 K), the water saturation is 0.3, and the in-situ stresses in the horizontal and vertical directions are 30 MPa and 36 MPa, respectively . Such field data provide a basis for the model development in this numerical work, ensuring this model is as close to the real formation conditions as possible.   . This pattern is axisymmetric, so a two-dimensional representative is selected for the numerical simulation (Fig. 2b). According to the geological parameters, the model height is determined as 100 m and is equally divided into ten sections to reveal the vertical heterogeneity of the reservoir properties. The model width is set 300 m to avoid possible influence of the right boundary condition (e.g., formation overpressure) (Li et al. 2017). Accordingly, the simulation area of 100 m × 300 m is divided into 5085 elements, and all boundaries are airtight for flux, except for the right boundary with the initial gas pressure being set. As the CS-EGR operation is conducted in the depleted conditions of the Wufeng-Longmaxi shale, the initial pressure of the simulative reservoir is set as 2 MPa (that is, 0.3 MPa for CO 2 and 1.7 MPa for CH 4 ). In this modeling process, under a temperature of 305 K, constant pressures of 0.1 MPa and 7 MPa are applied to the PWs and IW during the CS-EGR process, respectively. To develop the mathematical model, the primary parameters (see "Appendix") are obtained from the existing literatures (Dahaghi 2010;Sun et al. 2013;Fathi and Akkutlu 2014;Elsworth 2015, 2019;Fan et al. 2019a, b, c;Ma et al. 2020;Zhao et al. 2020), while the key heterogeneous information on the Wufeng-Longmaxi shale is derived from real fieldwork ( Fig. 1).

Model description for the Wufeng-Longmaxi shale
As shown in Fig. 2b, three IW candidates and three PWs candidates are organized to investigate how the relative locations of IW and PWs affect the free gas behavior of the CS-EGR operation in shales. Among these nine IW-PWs combos, when the IW is located at the lower section and is 50 m far away from the PWs in the horizontal direction, the modeling case is labeled as L50. In a similar fashion, the rest of the cases are marked as L100, L150, M50, M100, M150, H50, H100, and H150, respectively (Fig. 2). Herein, the layer 3, interface of layers 5-6, and layer 8, where the IW candidate is set, represent the upper, middle, and lower sections of the simulative shale reservoir, respectively.

THM coupling process and governing equations
Since many parameters are involved in this simulation work (see "Appendix"), CS-EGR in shales is identified as a complicated process and involves integral feedback in the THM coupling phenomenon. Basically, the hydraulic field contains complex mass transport of binary gases (CO 2 /CH 4 ) with the gas-water represented as two-phase flows superposed on competitive non-isothermal adsorption (Fig. 2c). Once the shale reservoir captures the nonthermally equilibrated CO 2 , heat transfer, such as thermal conduction/convection, occurs among the CH 4 , CO 2 , water, and shale skeleton, which is accompanied by energy release/adsorption associated with gas adsorption/desorption, in turn impacting on the thermal field (Fan et al. 2019a). Furthermore, the dynamic variation in the hydraulic and thermal fields can interfere with the mechanical field (e.g., shrink/swelling of the shale matrix), influencing the anisotropic porosity/permeability and thus affecting the convective fluxes of the water, gas, and energy (Fan et al. 2019a). The shale reservoir in this work is described as a dualporosity media (Fig. 2c), which is used extensively in previous numerical works (Bandis et al. 1983;Nassir et al. 2014;Elsworth 2015, 2019;Zhou et al. 2018;Zhao et al. 2020). According to these works above, a series of governing equations representing rock deformation, competitive CH 4 /CO 2 sorption, gas/water two-phase flow, and thermal expansion are used for the CS-EGR model development in this work.

Results and discussion
In the dual-porosity media, fractures and matrix pores are the internal spaces of shale for trapping free gases (CO 2 / CH 4 ). In this work, four types of free gases are classified to investigate the behavior of free CO 2 /CH 4 during the CS-EGR process in the depleted Wufeng-Longmaxi shale-CH 4 in the fracture, CO 2 in the fracture, CH 4 in the matrix and CO 2 in the matrix. Accordingly, this work separately studies the performance of free CH 4 or free CO 2 , and then the dynamic interaction between free CH 4 and free CO 2 during the recovery enhancement of shale gas by CO 2 injection into the depleted Wufeng-Longmaxi shale, under variable IW-PWs locations. Note that this numerical model's thickness value is 1 m when the gas content is calculated, where the involved gas content indicates the content of gas in free phase.

Content variation of free CH 4 during shale-based CS-EGR process
Prior to the involvement of the injected CO 2 , the original content of free CH 4 both in the fracture and in the matrix is strongly heterogeneous in the vertical direction, where the free CH 4 in the matrix is visibly more than that in fracture (Fig. 3). This heterogeneous tendency is similar to the vertical variation of the fracture/matrix porosity of the Wufeng-Longmaxi shale shown in Fig. 1, owing to the free gas content being directly determined by the free space under the depleted reservoir pressure. Once the CO 2 injection starts, the content variation of the free CH 4 in the fracture/matrix occurs. Taking M100 as an example, it is found that the free CH 4 in the fracture continuously varies during the CS-EGR process, intuitively reflecting the accumulation of free CH 4 in the left side of the PWs and the decrease in that in the right side of the PWs (Fig. 4a). For the same case, the cloud pictures qualitatively indicate an overall reduction of free CH 4 in the matrix over time after CO 2 is injected into the shale reservoir (Fig. 4b). These variations address the concept that the CO 2 injection can facilitate the migration of free CH 4 in the depleted Wufeng-Longmaxi shale reservoir.
After a production period of 10,000 d, the CS-EGR operation enables the content difference of free CH 4 in the fracture/matrix under variable IW-PWs combos. According to the cloud pictures, at a fixed IW location, a farther IW-PWs distance in the horizontal direction makes more CH 4 in the free phase trapped in the shale fracture and matrix after 10,000 d of CS-EGR operation (Fig. 5a, b). Comparatively, when PWs are located as fixed, there are less obvious changes for the free CH 4 in the shale fracture/matrix along with the variation of the IW location. Anyhow, after 10,000 d of CS-EGR production, the vertical heterogeneity emerges for all outputs of free CH 4 under different IW-PWs relative locations (Fig. 5a, b).
To further display the heterogeneous characteristics of free CH 4 in the Wufeng-Longmaxi shale, the content variation of free CH 4 in each layer is concluded at different points in time during the CS-EGR process. As shown in Fig. 6, the free CH 4 content in each layer is dramatically affected by the relative locations of IW-PWs. A shorter IW-PWs distance (fixed IW location) makes the CO 2 injection-induced decrement of the free CH 4 content greater, which is more obvious for the layers near the IW. Besides, when the PWs location is fixed, a smaller decrement goes to the free CH 4 content in the layers close to the IW. Moreover, in some operating cases, the content decrement is negative for the free CH 4 in the layers near IW (e.g., Fig. 6a-c, e, f, h, i), revealing that the CO 2 injection facilitates the increment of free CH 4 in shale fractures at the initial CS-EGR phase. Taking H150 as an example, the free CH 4 content in the fracture of layer 1 first increases and then decreases (Fig. 7a), while that of layer 8 decreases monotonously (Fig. 7b), after CO 2 is injected. This phenomenon is caused by the content variation of free CH 4 in the matrix. Using the same example of H150, the massive accumulation of free CH 4 in the matrix, caused by the displacement of the adsorbed CH 4 by CO 2 injection, allows the release of the free CH 4 in the matrix into the fracture, regarding layer 1 which is close to the IW location (Fig. 7a). Comparatively, for layer 8 of H150, the decrease of free CH 4 in the fracture/matrix is monotonic and mainly controlled by the release of CH 4 pressure in the shale reservoir (Fig. 7b).
With regard to the mass of CH 4 in the matrix, its decrement varies in a heterogeneous way for all layers, under  all IW-PWs combos (Fig. 8). Basically, the CH 4 content in the matrix at a depth of 1275 m (layer 8) decreases more than that in the rest of the layers, probably because the highest content of the original free CH 4 in the matrix (Fig. 3) and the highest permeability (Fig. 1) promotes the CH 4 release from the bottom PW in layer 8. When IW  has a fixed location, the content decrement of free CH 4 in the matrix at a shorter horizontal distance of IW-PWs is greater than that for a longer one. By comparison, the IW location only marginally affects the performance of free CH 4 in the matrix in a fixed situation of PWs, where an IW at the bottom layer makes the decrement of CH 4 content become slightly less than that at the upper layer (Fig. 8).
Besides, no matter where the IW and PWs are, the content of free CH 4 in the matrix is invariably greater than that in the fracture during the CS-EGR process, in that the matrix porosity φ m is considerably higher than the fracture porosity φ f . This phenomenon is illustrated by the examples L150 and M100, where the ratio of CH 4 in the matrix relative to that in the fracture decrease rapidly, then slow by after CO 2 involvement, which however is always greater than 1 (Fig. 9). Figure 9 also indicates this dynamic ratio of CH 4 in the matrix to CH 4 in the fracture differs under different IW-PWs patterns, in which the detailed mechanism will be described in future work.   Fig. 9 CO 2 injection-induced variation of the ratio of CH 4 in the matrix relative to that in the fracture during the CS-EGR process

Content variation of free CO 2 during shale-based CS-EGR process
As exhibited by the representative case M100, the mass content and distribution area of the free CO 2 in the reservoir gradually become greater with time after CO 2 is injected into the Wufeng-Longmaxi shale (Fig. 4c, d). Here, the variation is characterized as heterogeneous in the vertical direction, in which the performance of free CO 2 in the fracture differs from that in the matrix. This heterogeneous variation leads to the variable outcomes of free CO 2 content, either in the fracture or in the matrix during 10,000 d of CS-EGR operation, under different relative locations of IW-PWs (Fig. 5c,  d). From the qualitative perspective, after 10,000 d of CO 2 injection, the free CO 2 tends to be trapped more at a longer IW-PWs distance with a fixed IW location, while a bottom IW location enables more CO 2 in the free phase to be trapped than an upper one. Herein, this phenomenon is general and is suitable for the free CO 2 , both in the fracture and in the matrix (Fig. 5c, d).
During the CS-EGR process, the fracture is the first space in the shale reservoir to encounter the injected CO 2 ; therefore, the relative locations of IW and PWs significantly affect the performance of free CO 2 in the fracture, as shown by the numerical outputs of every single layer (Fig. 10). For all IW-PWs combos, the layers close to the IW meet the CO 2 first and tend to trap more CO 2 in the fracture. For example, a 2000 d of CS-EGR operation enables a considerable amount of free CO 2 to be trapped in layers 8-10 but less free CO 2 to be trapped in layers 1-3 in the case L50; furthermore, after CO 2 is injected for 10,000 d, the content of free CO 2 in the fracture of layers 8-10 is obviously more than that of rest layers, for the example of case L50 (Fig. 10g). Besides, it is also noted that a longer IW-PWs distance allows more free CO 2 to be trapped in the fracture of each layer than a shorter one, when the IW is fixed (Fig. 10). This is due to a longer path for the CO 2 migration that usually corresponds to a greater area for CO 2 accumulation in the fracture.
As for the free CO 2 in the matrix, its content in each layer is variable during the process of CO 2 injection, under all H50 H100 H150 Fig. 10 Accumulation of free CO 2 in the fracture during the CS-EGR process for different IW-PWs combos IW-PWs combos (Fig. 11). At first glance, for all IW-PWs situations, the matrix of layer 8 (depth of about 1275 m) performs effectively in the CO 2 trapping among all layers (Fig. 11), a result of the significant matrix porosity φ m of layer 8 when compared with other layers (Fig. 1). In addition, the IW location affects the accumulation of free CO 2 in the matrix at a fixed situation of PWs, where the layers near the IW are likely to trap more free CO 2 in the matrix-similar to the effect for free CO 2 in the fracture. While for the situation of the fixed IW location, a longer IW-PWs distance makes more free CO 2 accumulated in the matrix, which is similar to the performance of free CO 2 in the fracture, with similar reasoning (Fig. 11).
Because the matrix porosity φ m is significantly greater than the fracture porosity φ f , for each layer, the content of free CO 2 in the matrix is higher than that in the fracture, revealed by Figs. 10 and 11. Quantitatively, the examples of H50 and M100 indicate that the ratio of free CO 2 in the matrix relative to that in the fracture is invariably superior during the whole CS-EGR process (Fig. 12). Nonetheless, Ratio of CO 2 in matrix to CO 2 in fracture H50 M100 Fig. 12 Ratio of CO 2 in the matrix to that in the fracture after CO 2 injection during the CS-EGR process in Fig. 12, a sharp decline occurs in the very initial period of CO 2 injection, resulting from the injected CO 2 mainly staying in the fracture before arriving in the matrix pores; afterwards, this ratio has a small increase after CO 2 enters into the matrix pores of shale reservoir, and then tends to be stable.

Interactive behavior of free CO 2 -CH 4 during the shale-based CS-EGR process
In the process of CS-EGR operation, the behavior of free gases is complex. For example, the free CH 4 in the matrix after CO 2 injection contains the original free CH 4 in the matrix and the desorbed CH 4 (originally in adsorbed) by CO 2 displacement. So, the interactive behavior between the free CO 2 and the free CH 4 is dynamic and complicated during the CS-EGR process in the shale reservoir. In Fig. 13, the content variation of free CH 4 and CO 2 echoes the outputs shown by Fig. 5, implicating the intensified impact on the performance of free gases from the relative locations of IW and PWs. Overall, the injection of CO 2 into the Wufeng-Longmaxi shale allows the content of free CH 4 and free CO 2 in the reservoir to decrease and increase, respectively (Fig. 13). Although the tendency of content variation for free CH 4 and free CO 2 is opposite, there are some common elements; for instance, the residual free CH 4 and trapped free CO 2 for case L150 are the greatest, while those for H50/M50 are the lowest (Fig. 13).
According to the statistics, the content fraction of free CH 4 among all free gases both in the fracture and in the matrix gets a continuous decrease with time, during the CS-EGR operation in the depleted Wufeng-Longmaxi shale (Fig. 14). Therein, the fraction of free CH 4 content in the fracture is consistently lower than that in the matrix, on the basis of two representative examples (Fig. 14). This phenomenon possibly is due to the extracted CH 4 from the PWs containing the free CH 4 in the fracture in the right side of the PWs (Fig. 4a), while the CH 4 -CO 2 displacement that partly supplies the free CH 4 in the matrix mainly exists in the left side of PWs (Fig. 4). However, this hypothesis needs more attention. In addition, the variable fraction of free CH 4 in the fracture and that in the matrix together form the content fraction of free CH 4 relative to all free gases in the whole shale reservoir. Herein, reflected by the slope of change curves, the decreasing tendencies in Fig. 14 are of the "fast followed by slow" type, which suggests the proportion of free CH 4 among all free gases has a tendency to be constant after a sufficient time of CO 2 injection. Therefore, it can be speculated that the CS-EGR operation probably ends when the fraction of free CH 4 among all free gases changes insignificantly with time. Furthermore, after a 10,000 d of CS-EGR production, the resulting proportion of free CH 4 and free CO 2 differs, under different IW-PWs combos (Fig. 15). For all cases, no matter whether in the fracture or the matrix, the content of free CO 2 is primarily greater than that of free CH 4 . After the CS-EGR operation runs for 10,000 d, for all IW-PWs combos, the content proportions of free CO 2 among all free gases in the matrix, the fracture and the whole reservoir are 65.8%, 69.5%, and 66.6% on average, respectively. Comparatively, it seems that a shorter IW-PWs distance (cases H50, M50 and L50) tends to enable a relatively higher proportion of free CH 4 both in the fracture and in the matrix, and thus in the whole reservoir (Fig. 15).

Location selection of IW-PWs for desired performance of free CO 2 -CH 4
Since different IW-PWs locations facilitate variable interactive behavior of free CO 2 and free CH 4 during the shalebased CS-EGR process, it is necessary to make an appropriate selection on the relative locations of IW and PWs to achieve the desired purpose. For this selection, a parameter called the recovery efficiency of free CH 4 ( f free-CH 4 , for short) is defined, where C o and C r are the content of the original free CH 4 (before CO 2 involvement) and the residual free CH 4 (after CO 2 involvement) in the fracture/matrix, respectively, kg.
Herein, a higher f free-CH 4 value indicates that more CH 4 in the free phase is recovered from PWs, and vice versa. As exhibited in Fig. 16, the vertical heterogeneity is shown in the f free-CH 4 value under each IW-PWs combo, during the CS-EGR operation in the Wufeng-Longmaxi shale, which is codetermined by the reservoir properties and the IW-PWs strategy. For all IW-PWs cases, the vertical heterogeneity of f free-CH 4 has a similar relation to that of the free CH 4 in the matrix; that is, the f free-CH 4 for bottom layers is higher than that for the upper layers (Fig. 16). With regard to the free CH 4 in the fracture, the f free-CH 4 is affected by the IW location, and when the IW locates at the upper layers (or bottom layers), the f free-CH 4 of bottom layers (or upper layers) becomes higher (Fig. 16).
The performance of free CH 4 /CO 2 in the whole reservoir consists of that in the fracture and in the matrix of each layer. Basically, a longer IW-PWs distance generates Free gases in matrix Fig. 15 Proportion of free CH 4 and free CO 2 in the Wufeng-Longmaxi shale after 10,000 d of CS-EGR operation a lower recovered content of free CH 4 (Fig. 17a), accompanied with a higher trapped content of free CO 2 (Fig. 17b), both in the fracture and the matrix. Accordingly, an appropriate IW-PWs strategy can be selected for different expected targets. For example, if the CS-EGR operation aims to trap more CO 2 in the free phase in the depleted Wufeng-Longmaxi shale, the IW should be located in the bottom layers and have a longer horizontal distance with the PWs (like L150 in this work) (Fig. 17). One more example, if the CS-EGR operation is designed to get more CH 4 in the free phase recovered from the depleted Wufeng-Longmaxi shale, the IW location is flexible and only a shorter IW-PWs distance is needed, such as H50, M50, and L50 in this work (Fig. 17).

Conclusions
In developing a novel THM coupling model, the performance of free CH 4 and free CO 2 during the CS-EGR process in the Wufeng-Longmaxi shale is clearly obtained. The main conclusions are.
Vertical heterogeneity exists in the content of free CH 4 or free CO 2 in the fracture/matrix throughout the whole process of CS-EGR operation, codetermined by the reservoir properties and the IW-PWs strategy. Because the matrix porosity φ m is significantly higher than the fracture porosity φ f , the free CH 4 /CO 2 in the matrix is much higher than that in the fracture for either layer or the whole   Trapped free CO 2 , 10 3 kg H50 H100 H150 M50 M100 M150 L50 L100 L150 H50 H100 H150 M50 M100 M150 L50 L100 L150

Fig. 17
Recovered free CH 4 and trapped free CO 2 in the Wufeng-Longmaxi shale after 10,000 d of CS-EGR operation reservoir. After CO 2 involvement, the ratio of free CH 4 / CO 2 in the matrix relative to that in the fracture declines and tends to be stable with time, where the change behavior is different for the free CH 4 and free CO 2 . For the free CH 4 in the fracture/matrix, its recovery is lower at a longer IW-PWs distance (fixed IW location) and is insignificantly affected by the variation of IW location at a PW location during the CS-EGR operation. For the free CO 2 in the fracture/matrix, it is trapped more at a longer IW-PWs distance (fixed IW location) and tends to be more trapped when the IW locates at bottom layers (fixed location of PWs). After CO 2 is injected into the Wufeng-Longmaxi shale, the free CH 4 content in the fracture/matrix of the layers near the IW location increases first and decreases later, while that of the layers far away from the IW location suffers a monotonic decrease.
During the CS-EGR operation in the Wufeng-Longmaxi shale, the content of free CH 4 among all free gases in the fracture/matrix has a continuous decline with time-in a "fast followed by slow" way. A 10,000 d of CO 2 injection enables the content of free CO 2 to be greater than that of free CH 4 in the fracture/matrix, in which a shorter IW-PWs distance results in a relatively higher proportion of free CH 4 . In addition, when the IW locates at the bottom layers and has a farther distance to PWs, more CO 2 in the free phase tends to be trapped in the depleted Wufeng-Longmaxi shale; furthermore, no matter where the IW is, a shorter IW-PWs distance is helpful for getting more CH 4 in the free phase recovered from the depleted Wufeng-Longmaxi shale.