Numerical Simulation of Gas Production and Reservoir Stability during CO 2 Exchange in Natural Gas Hydrate Reservoir

: The prediction of gas productivity and reservoir stability of natural gas hydrate (NGH) reservoirs plays a vital role in the exploitation of NGH. In this study, we developed a THMC (thermal-hydrodynamic-mechanical-chemical) numerical model for the simulation of gas production behavior and the reservoir response. The model can describe the phase change, multiphase ﬂow in porous media, heat transfer, and deformation behavior during the exploitation of NGH reservoirs. Two different production scenarios were employed for the simulation: depressurization and depressur-ization coupled with CO 2 exchange. The simulation results suggested that the injection of CO 2 promotes the dissociation of NGH between the injection well and the production well compared with depressurization only. The cumulative production of gas and water increased by 27.88% and 2.90%, respectively, based on 2000 days of production simulation. In addition, the subsidence of the NGH reservoir was lower in the CO 2 exchange case compared with the single depressurization case for the same amount of cumulative gas production. The simulation results suggested that CO 2 exchange in NGH reservoirs alleviates the issue of reservoir subsidence during production and maintains good reservoir stability. The results of this study can be used to provide guidance on ﬁeld production from marine NGH reservoirs.


Introduction
Natural gas hydrate (NGH), an ice-like substance formed by natural gas (i.e., methane gas) and water under low-temperature and high-pressure conditions, is widely distributed in permafrost locations and sediments below the seafloor around the world, and has a large resource volume [1][2][3].With the ever-increasing global demand for clean energy, the exploitation of NGH reservoirs and novel production technology of NGH have attracted worldwide attention [4][5][6].Methods for gas production from NGH mainly include depressurization, thermal stimulation, injection of inhibitors, and CO 2 exchange [7][8][9].However, production efficiency may not be desirable in every single method.In addition, combining hydrate dissociation and carbon capture & storage would be a potential method, especially regarding marine storage under the goals of carbon peaking and carbon neutrality [10].From an environmental perspective, marine hydrate reservoir development using depressurization and CO 2 exchange may be an attractive way to implement the CO 2 capture and storage (CCS) technique [11].
Energies 2022, 15, 8968 3 of 17 deformation in depressurization and CO 2 exchange cases.The results suggested that the adoption of CO 2 exchange significantly alleviated seafloor subsidence.
However, there are limited studies in the literature with regard to the feasibility of CO 2 exchange for hydrate reservoirs using a coupled THMC model.In addition, one critical question that warrants examination in the depressurization field test is the stability of the NGH reservoir.Our motivations were to investigate the following: (a) the degree of deterioration of NGH reservoir stability; (b) the exact level of NGH reservoir subsidence during production tests; and (c) the effective method to mitigate the reservoir subsidence.Based on these considerations, the method of depressurization and CO 2 exchange with the advantages of enhancing gas production and mitigating reservoir deformation was explored.In this study, the THMC model for gas production from hydrate reservoir was established and implemented in CMG-STARS to simulate fluid production behavior.Moreover, NGH reservoir response were considered using two different production techniques: (a) depressurization only; and (b) depressurization coupled with CO 2 exchange.The influence of CO 2 exchange on the fluid production behavior and the reservoir stability were analyzed in detail.The simulation results of this study can provide guidance on the field production of marine NGH reservoirs using the depressurization coupled with CO 2 exchange technique.

THMC Model and Solution Approach
Fluid production from NGH reservoirs involves complex physical and chemical interactions in geological media, i.e., hydrate phase change, multiphase flow, heat transfer, and reservoir deformation.The THMC model should account for the non-isothermal methane gas release, phase behavior, and fluid flow under the conditions typical for methane hydrate deposits (i.e., in the permafrost and subsea hydrate-bearing sediments) by solving the coupled equations of fluids mass and heat balance associated with such systems.The governing equations for the THMC model are as follows.

Hydrate Reaction Kinetic Model
The dissociation reaction of NGH be described as follows: where M is the gas molecule (CH 4 or CO 2 in this paper), and n represents the hydration number, which is determined by the hydrate structure and the cage occupancy.The Kim-Bishnoi model [43] was used to describe the hydrate phase change behaviour in Equation (2) as follows: where c h represents the hydrate concentration, mol where k 0 = k d 0 when hydrate is dissociating and k 0 = k f 0 when the hydrate is forming; ∆E is the activation energy, J•mol −1 ; R is the universal gas constant, J•mol −1 •K −1 ; T is the temperature, K; A is the hydrate surface area participating in the kinetic reaction, m −1 , where A = A d when hydrate is dissociating and A = A f when the hydrate is forming; p e is the equilibrium pressure of CH 4 or CO 2 hydrate at the system temperature conditions, Pa; p g is the corresponding pressure of gas phase, Pa.Assuming that hydrates in porous media are composed of spherical particles with surface area A HS , the surface area of hydrates during the dissociation and formation process per unit volume of porous media can be approximated using Equations ( 3) and ( 4) following [44]: where S w and S h represent the phase saturation of water and hydrate in the pore space, respectively, and φ is the porosity of the porous medium.Thus, the reaction rate for hydrate dissociation and formation can be described in Equations ( 5) and ( 6): The phase equilibrium relation of NGH is used following the model of Kamath [45]: where P eq represents equilibrium pressure at the system temperature conditions in kPa, T eq represents the equilibrium temperature at the system pressure conditions in K; e 1 and e 2 are the correlation coefficients obtained from the hydrate phase equilibrium measurement.In this study, the phase equilibrium curves of CH 4 and CO 2 hydrate were used as follows in Figure 1 [46].

𝐴 = 𝜙 𝐴 𝑆 𝑆
(3) where Sw and Sh represent the phase saturation of water and hydrate in the pore space, respectively, and ϕ is the porosity of the porous medium.Thus, the reaction rate for hydrate dissociation and formation can be described in Equations ( 5) and ( 6): The phase equilibrium relation of NGH is used following the model of Kamath [45]: where Peq represents equilibrium pressure at the system temperature conditions in kPa, Teq represents the equilibrium temperature at the system pressure conditions in K; e1 and e2 are the correlation coefficients obtained from the hydrate phase equilibrium measurement.
In this study, the phase equilibrium curves of CH4 and CO2 hydrate were used as follows in Figure 1 [46].

Mass Balance Equations
In the numerical model, all three phases (gas, water, and hydrate) were considered, where the hydrate phase as a solid state is not allowed to flow.Considering the water and gas two-phase flow in porous media during the dissociation of hydrate, the mass conservation equations can be represented as: Energies 2022, 15, 8968 where ρ g , ρ w , ρ h represents the density of gas phase, water phase and hydrate phase, respectively, kg•m −3 ; S g is the gas saturation; v g , v w represents the flow rate of gas and water, respectively, m•s −1 ; q g , q w are the source or sink term of gas and water, kg•m −3 •s −1 .dm h /dt, dm g /dt, dm w /dt represents the hydrate dissociation rate and the production rate of dissociated gas and water from hydrate, kg•m −3 •s −1 .

Energy Balance Equation
Hydrate dissociation consumes an abundance of heat, so this important process should be considered in NGH production.The energy conservation equation can be described in Equation (11) considering heat transfer (thermal conduction and heat fluid convection), and the endothermic term caused by hydrate dissociation along with external heat injection: where ρ s is the density of rock skeleton, kg where λ s , λ h , λ g , λ w represents the thermal conductivity of sandy media, hydrate phase, gas phase, and water phase, respectively, W•m −1 •K −1 .

Geomechanics Equation
The three key elements in the description of geomechanical responses are stress, strain, and displacement.Neglecting the change in momentum, the stress equilibrium equation can be expressed as follows [47]: where σ is the total stress tensor, Pa, g is acceleration of gravity, m/s 2 , and ρ r is the density of sediment, kg/m 3 .
Based on the principle of effective stress, the relationship between stress and pore pressure is: where σ is the effective stress tensor, Pa, α is Biot's number, p is the pore pressure, Pa, and I is the unit matrix.The relationship of stress and strain for hydrate-bearing sediment is [47]: where C is the stiffness tensor, Pa, ε is the strain tensor, η is the coefficient of thermal expansion, Pa•K −1 , and ∆T is the temperature difference, K.
Considering the assumption of small strain, the relationship between strain and displacement is as follows [48]: Energies 2022, 15, 8968 6 of 17 where u is the displacement, m.Equations ( 14)-( 16) can be substituted into Equation ( 13) and the displacement can be simplified in Equation ( 17) [47]: where C 1 and C 2 are the correlation parameters.

Dynamic Model of NGH Reservoir Thermophysical Parameters
In addition to the main balance equation, other relationships of the variables should be illustrated.
The summation of all phases equals to unity: The relative permeability model is used as [49]: where k rw,max is the largest relative permeability of the water phase; S wr is the irreducible water saturation; k rg,max is the largest relative permeability of the gas phase; S gr is the critical gas saturation; c w , c g are the fitting parameters.
The latent heat of hydrate phase change is given by the equation [45]: In the coupled THMC numerical model, the change in pressure and temperature results in hydrate dissociation and subsequently yielded changes in stress and deformation of NGH reservoir.The above changes further result in the change of porosity and permeability, which in turn affect the NGH dissociation rate and the fluid seepage behaviors.Thus, in the THMC coupled numerical model, it is necessary to set up a dynamic model for the description of the change of thermophysical parameters, including porosity, permeability, and Young's modulus, etc.
The change of porosity is simultaneously affected by the volumetric strain and hydrate dissociation, and its dynamic change equation is used in Equation ( 22) following [50]: where φ is the porosity at hydrate saturation is S h and φ 0 is the initial porosity.
The permeability reduction model in Equation ( 23) proposed by Moridis et al. [51] is used considering the presence of solid phase, i.e., hydrate or ice, in sandy media: where k 0 is the initial absolute permeability of the formation, µm 2 ; k(φ) is the permeability changed with φ, µm 2 ; φ c is the critical porosity; S s the solid phase saturation, and S s = S h + S I when considering the existence of both hydrate and ice solid phases; N is the index of permeability reduction.The geomechanical properties, such as Young's modulus, Poisson's ratio, internal friction angle etc., of hydrate-bearing sediments, not only depend on the solid skeleton properties, but are also correlated with the amount of the solid hydrate phase existing in the pore space.Thus, the dynamic model of the geomechanical properties can be described in Equations ( 24)-( 26) following [52]: Young's modulus is computed by Poisson's ratio is computed by Friction angle is computed by where m is the empirical coefficient; E b is the Young's modulus of the deposit and the fluid, kPa; E h is the Young's modulus of the solid hydrate, kPa; v b is the Poisson's ratio of the deposit and the fluid; v h is the Poisson's ratio of the solid hydrate; ϕ b is the angle of internal friction of deposit and fluid, degrees; ϕ h is the angle of internal friction of the solid hydrate degrees.

Solution Approach
The iterative solution approach of the coupled THMC numerical model is presented in the flow chart in Figure 2 considering the NGH reservoir seepage and the geomechanical responses.

Results and Discussions
A three-dimensional NGH reservoir geological model with dimensions 500 m × 500 m × 240 m was setup.The NGH reservoir is divided into three layers vertically, including a 150 m layer representing overburden, a 60 m layer representing NGH reservoir, and a 30 m layer representing the underburden.Figure 3 shows a schematic diagram of the geological model.Two vertical wells are designed in the model with one functioning as a depressurization production well and the other as an injection well for CO2.The distance between two wells is 170 m.To avoid the interference of the water intrusion from the underburden, the production section of the two wells is located at the upper 50 m of the NGH layer.

Results and Discussions
A three-dimensional NGH reservoir geological model with dimensions 500 m × 500 m × 240 m was setup.The NGH reservoir is divided into three layers vertically, including a 150 m layer representing overburden, a 60 m layer representing NGH reservoir, and a 30 m layer representing the underburden.Figure 3 shows a schematic diagram of the geological model.Two vertical wells are designed in the model with one functioning as a depressurization production well and the other as an injection well for CO 2 .The distance between two wells is 170 m.To avoid the interference of the water intrusion from the underburden, the production section of the two wells is located at the upper 50 m of the NGH layer.Two simulation cases with different production strategies are conducted for comparison: Case A is designed using depressurization only with a bottom hole pressure (BHP) of 5.0 MPa; Case B is designed using the combination of depressurization and CO2 exchange method.In the initial 500 days of Case B, depressurization using a BHP of 5 MPa was employed, then CO2 was injected at a rate of 5000 m 3 •d −1 .The maximum CO2 injection pressure is 18.0 MPa.The simulation time for both cases is set as 2000 days.Table 1 lists the main parameters of the thermophysical properties.Two simulation cases with different production strategies are conducted for comparison: Case A is designed using depressurization only with a bottom hole pressure (BHP) of 5.0 MPa; Case B is designed using the combination of depressurization and CO 2 exchange method.In the initial 500 days of Case B, depressurization using a BHP of 5 MPa was employed, then CO 2 was injected at a rate of 5000 m 3 •d −1 .The maximum CO 2 injection pressure is 18.0 MPa.The simulation time for both cases is set as 2000 days.Table 1 lists the main parameters of the thermophysical properties.

Analysis of Fluid Production
Figure 4 shows a comparison of the gas and water production (gas production rate Q G , cumulative gas production C G , water production rate Q W , cumulative water production C W ) over time in the afore-mentioned two cases.It was observed that the dynamic of gas and water production shows an identical trend for the initial 500 days due to the same depressurization method applied.After the injection of CO 2 starting at 500 days, the production rate of CH 4 in Case B shows an immediately increasing trend compared with Case A as evidenced by the evolution of Q G and C G .This can be explained by the increase of the pressure near the injection well due to CO 2 injection, which drives the flow of the dissociated CH 4 from NGH to the production well.A significant boost on the gas production rate was observed at about 610 d.This can be explained by the injected CO 2 forming CO 2 hydrate with the free water remaining in the reservoir or the dissociated water from NGH, and the heat released by the formation of CO 2 hydrate promotes the dissociation of NGH further.This resulted in a significant increase in gas production only, but less of an increase in water production.
Figure 4 shows a comparison of the gas and water production (gas production rate QG, cumulative gas production CG, water production rate QW, cumulative water production CW) over time in the afore-mentioned two cases.It was observed that the dynamic of gas and water production shows an identical trend for the initial 500 days due to the same depressurization method applied.After the injection of CO2 starting at 500 days, the production rate of CH4 in Case B shows an immediately increasing trend compared with Case A as evidenced by the evolution of QG and CG.This can be explained by the increase of the pressure near the injection well due to CO2 injection, which drives the flow of the dissociated CH4 from NGH to the production well.A significant boost on the gas production rate was observed at about 610 d.This can be explained by the injected CO2 forming CO2 hydrate with the free water remaining in the reservoir or the dissociated water from NGH, and the heat released by the formation of CO2 hydrate promotes the dissociation of NGH further.This resulted in a significant increase in gas production only, but less of an increase in water production.Based on the above analysis, the injection of CO2 and CO2 exchange promotes the dissociation of NGH.For a production period of 2000 days, CG and CW in Case B (depressurization coupled with CO2 exchange) increased by about 27.88% and 2.90%, respectively, compared with Case A (depressurization only).Based on the above analysis, the injection of CO 2 and CO 2 exchange promotes the dissociation of NGH.For a production period of 2000 days, C G and C W in Case B (depressurization coupled with CO 2 exchange) increased by about 27.88% and 2.90%, respectively, compared with Case A (depressurization only).

Evolution of P, T, and S h
Figure 5 shows the evolution of P at 500 days, 1000 days, and 2000 days of production (left figure panel shows Case A using depressurization method and right figure panel shows Case B using depressurization with CO 2 exchange method).P reduction spreads from the production well to the surroundings gradually.The fastest P drop (DP) located near the well and DP gradually decreases away from the well.With the increase of production time, the degree and the area affected by the pressure reduction gradually increase.Due to the permeability difference between the NGH layer (1100 × 10 −3 µm 2 ) and the overburden and underburden (1 × 10 −3 µm 2 ), pressure drop mainly propagates in the NGH layer.
shows Case B using depressurization with CO2 exchange method).P reduction spreads from the production well to the surroundings gradually.The fastest P drop (DP) located near the well and DP gradually decreases away from the well.With the increase of production time, the degree and the area affected by the pressure reduction gradually increase.Due to the permeability difference between the NGH layer (1100 × 10 −3 μm 2 ) and the overburden and underburden (1 × 10 −3 μm 2 ), pressure drop mainly propagates in the NGH layer.In Case B, it is observed that P increases near the injection well due to CO2 injection.However, the formation of CO2 hydrate accelerates the dissociation of NGH, which possibly increases the porosity and permeability of the NGH reservoir.It further assists the pressure transmission in the NGH reservoir.This can be evidenced in Figure 5f at t = 2000 days, the area with P reduction is much wider and P near the production well is much lower than that in Case A (see Figure 5c).
Figure 6 shows the evolution of T in Case A and Case B at t = 500 days, 1000 days, and 2000 days.Due to the endothermic process of NGH dissociation (DHCH4 = 54.5 kJ/mol), the T near the production well drops rapidly.The T drop is decreasing with the increasing distance away from the production well.The low-T region expands over time with the increasing area of NGH dissociation.In Case B, with the injection of CO2 and the formation of CO2 hydrate releasing heat (DHCO2 = 57.9kJ/mol), T near the injection well increases and is relatively higher in this region.In addition, it was found that the direction of T increase is rightward from the injection well to the production well.This can be explained by the direction of CO2 exchange, which results in the unsymmetric T distribution near the injection well.Around the injection well, T is relatively higher towards the production well.On the other hand, due to the CO2 exchange results in a larger mass of NGH dissociation (an endothermic reaction), the T decrease near the production well is more significant in Case B compared with that in Case A at t = 2000 days.Overall, the T at the lower section of NGH reservoir is higher than that in the upper section based on the geothermal gradient.In Case B, it is observed that P increases near the injection well due to CO 2 injection.However, the formation of CO 2 hydrate accelerates the dissociation of NGH, which possibly increases the porosity and permeability of the NGH reservoir.It further assists the pressure transmission in the NGH reservoir.This can be evidenced in Figure 5f at t = 2000 days, the area with P reduction is much wider and P near the production well is much lower than that in Case A (see Figure 5c).
Figure 6 shows the evolution of T in Case A and Case B at t = 500 days, 1000 days, and 2000 days.Due to the endothermic process of NGH dissociation (DH CH4 = 54.5 kJ/mol), the T near the production well drops rapidly.The T drop is decreasing with the increasing distance away from the production well.The low-T region expands over time with the increasing area of NGH dissociation.In Case B, with the injection of CO 2 and the formation of CO 2 hydrate releasing heat (DH CO2 = 57.9kJ/mol), T near the injection well increases and is relatively higher in this region.In addition, it was found that the direction of T increase is rightward from the injection well to the production well.This can be explained by the direction of CO 2 exchange, which results in the unsymmetric T distribution near the injection well.Around the injection well, T is relatively higher towards the production well.On the other hand, due to the CO 2 exchange results in a larger mass of NGH dissociation (an endothermic reaction), the T decrease near the production well is more significant in Case B compared with that in Case A at t = 2000 days.Overall, the T at the lower section of NGH reservoir is higher than that in the upper section based on the geothermal gradient.
Figure 7 shows the evolution of NGH saturation in Case A and Case B at t = 500 days, 1000 days, and 2000 days.It can be informed that the dissociation of NGH propagates from the production well outward gradually.The region near the production well exhibits a faster dissociation rate compared with that far from the production well.Due to the geothermal gradient effect, NGH located at the lower section of the NGH reservoir practically dissociated more quickly than the upper section.In Case B, the dissociation rate of NGH near the injection well is enhanced with the injected CO 2 forming CO 2 hydrate and releasing heat.With the increasing amount of CO 2 injected and transported, NGH located in between the injection well and the production well dissociated at a much faster rate.It should be noted that in Case A under the depressurization method, NGH dissociation is symmetric and propagates outward based on the center of production well.While in Case B, the majority of NGH dissociation occurs in the region in between two wells.Based on the above analysis, the distance between the injection well and the production well has a profound relationship with the fluid production performance, which need to be further optimized in future studies.Figure 7 shows the evolution of NGH saturation in Case A and Case B at t = 500 days, 1000 days, and 2000 days.It can be informed that the dissociation of NGH propagates from the production well outward gradually.The region near the production well exhibits a faster dissociation rate compared with that far from the production well.Due to the geothermal gradient effect, NGH located at the lower section of the NGH reservoir practically dissociated more quickly than the upper section.In Case B, the dissociation rate of NGH near the injection well is enhanced with the injected CO2 forming CO2 hydrate and releasing heat.With the increasing amount of CO2 injected and transported, NGH located in between the injection well and the production well dissociated at a much faster rate.It should be noted that in Case A under the depressurization method, NGH dissociation is symmetric and propagates outward based on the center of production well.While in Case B, the majority of NGH dissociation occurs in the region in between two wells.Based on the above analysis, the distance between the injection well and the production well has a profound relationship with the fluid production performance, which need to be further optimized in future studies.Figure 8 shows the evolution of CO2 hydrate saturation in Case A and Case B at t = 500 days, 1000 days, and 2000 days.It can be seen in Figure 8 that the formation of CO2 hydrate initially occurred near the injection well and expanded outward from the injection well to the production well.The majority of the CO2 hydrate is located in the region between the injection well and the production well.Due to the geothermal gradient, the saturation of CO2 hydrate in the upper section is practically higher than that in the lower section.Figure 8 shows the evolution of CO 2 hydrate saturation in Case A and Case B at t = 500 days, 1000 days, and 2000 days.It can be seen in Figure 8 that the formation of CO 2 hydrate initially occurred near the injection well and expanded outward from the injection well to the production well.The majority of the CO 2 hydrate is located in the region between the injection well and the production well.Due to the geothermal gradient, the saturation of CO 2 hydrate in the upper section is practically higher than that in the lower section.
Figure 8 shows the evolution of CO2 hydrate saturation in Case A and Case B at t = 500 days, 1000 days, and 2000 days.It can be seen in Figure 8 that the formation of CO2 hydrate initially occurred near the injection well and expanded outward from the injection well to the production well.The majority of the CO2 hydrate is located in the region between the injection well and the production well.Due to the geothermal gradient, the saturation of CO2 hydrate in the upper section is practically higher than that in the lower section.

Effect of CO 2 Exchange on NGH Reservoir Stability
On one hand, the formation of CO 2 hydrate in the pore space increases the NGH reservoir stability.However, the exchange of CO 2 further promotes the dissociation of NGH, which results in the weakening of the NGH reservoir stability.For a fair comparison of the effect of CO 2 exchange on NGH reservoir stability, we compared the displacement of the NGH reservoir under the same C G = 5.81 × 10 7 m 3 in both cases.The production time is 1500 days in Case B compared with 1800 days in Case A, as shown in Figure 4a.
Figure 9 shows the comparison of the Young's Modulus in two cases (left as Case A and right as Case B) with the same C G .In the case of depressurization only, due to the dissociation of NGH, Young's modulus of the formation near the production well reduces, which reduces the reservoir capability to resist stress deformation.In Case B, although the dissociation of NGH reduces the formation strength partially, the formation of CO 2 hydrate further improves the formation strength.The formation of CO 2 hydrate increased Young's modulus of hydrate-bearing sediments significantly between the two wells.This preserves the NGH reservoir's ability to maintain its initial stability.

Effect of CO2 Exchange on NGH Reservoir Stability
On one hand, the formation of CO2 hydrate in the pore space increases the NGH reservoir stability.However, the exchange of CO2 further promotes the dissociation of NGH, which results in the weakening of the NGH reservoir stability.For a fair comparison of the effect of CO2 exchange on NGH reservoir stability, we compared the displacement of the NGH reservoir under the same CG = 5.81 × 10 7 m³ in both cases.The production time is 1500 days in Case B compared with 1800 days in Case A, as shown in Figure 4a.
Figure 9 shows the comparison of the Young's Modulus in two cases (left as Case A and right as Case B) with the same CG.In the case of depressurization only, due to the dissociation of NGH, Young's modulus of the formation near the production well reduces, which reduces the reservoir capability to resist stress deformation.In Case B, although the dissociation of NGH reduces the formation strength partially, the formation of CO2 hydrate further improves the formation strength.The formation of CO2 hydrate increased Young's modulus of hydrate-bearing sediments significantly between the two wells.This preserves the NGH reservoir's ability to maintain its initial stability.It was observed in Figure 10 that the maximum subsidence occurs at the interface of the overburden and NGH layer near the production well.The value of maximum subsidence is 1.06 m in Case B compared with that of 1.08 m in Case A. The maximum uplift occurs at the interface of the underburden and NGH layer near the production well.The uplift distance is 0.058 m in Case B compared with that of 0.057 m in Case A. Overall, the subsidence of the NGH reservoir using the method of depressurization coupled with CO 2 exchange is less than that in the case of depressurization only.The simulation results indicated that NGH dissociation assisted by CO 2 exchange practically offers an advantage for maintaining the reservoir stability.

Conclusions
This paper developed a combined numerical model of THMC (thermal-hydrodynamic-mechanical-chemical) to describe the behavior of gas production and reservoir reaction.Two different production scenarios have been employed for the simulation: depressurization and depressurization coupled with CO2 exchange.The results can be yielded as follows:

Conclusions
This paper developed a combined numerical model of THMC (thermal-hydrodynamicmechanical-chemical) to describe the behavior of gas production and reservoir reaction.Two different production scenarios have been employed for the simulation: depressurization and depressurization coupled with CO 2 exchange.The results can be yielded as follows: (a) A coupled THMC numerical model was developed in this study for the description of the complex phase change behavior, the gas-water two-phase flow, heat transfer, and formation deformation associated with NGH field-scale production; (b) The simulation results suggested that the injection of CO 2 promotes the dissociation of NGH compared with depressurization only.The cumulative production of gas and water increased by 27.88% and 2.90%, respectively based on 2000 days of production simulation; (c) In Case A (depressurization only), NGH dissociation is symmetric and propagates outward based on the center of production well.While in Case B (depressurization coupled with CO 2 exchange), the majority of NGH dissociation occurs in the region between the injection well and the production well; (d) NGH reservoir subsidence propagates outward based on the center of the production well in both cases.The maximum subsidence occurs at the interface of the overburden and NGH layer near the production well.The maximum uplift occurs at the interface of the underburden and NGH layer at the production well; (e) The subsidence of the NGH reservoir is less in Case B (depressurization coupled with CO 2 exchange) compared with that of Case A (depressurization only) at the same cumulative gas production.The results suggest CO 2 exchange in the NGH reservoirs alleviates the issue of reservoir subsidence during production and maintains good reservoir stability.
The results of this study can provide a new insight for high-efficient and safe gas extraction from marine deposits of NGH.

Figure 1 .
Figure 1.Phase equilibrium relation of CO2 and CH4 hydrate in this study.

Figure 2 .
Figure 2. Flow chart of the iterative solution approach of the THMC model.

Figure 2 .
Figure 2. Flow chart of the iterative solution approach of the THMC model.

Figure 3 .
Figure 3. Schematic diagram of the geological model.

Figure 4 .
Figure 4. Comparison of gas (a) and water (b) production in Case A (depressurization only) and Case B (depressurization coupled with CO2 exchange).

Figure 4 .
Figure 4. Comparison of gas (a) and water (b) production in Case A (depressurization only) and Case B (depressurization coupled with CO 2 exchange).

Figure 5 .
Figure 5. Evolution of P after 500 days, 1000 days, 2000 days ((a-c) for Case A and (d-f) for Case B).

Figure 6 .
Figure 6.Evolution of T after 500 days, 1000 days, 2000 days ((a-c) for Case A and (d-f) for Case B).

Figure 7 .
Figure 7. Evolution of NGH saturation after 500 days, 1000 days, 2000 days ((a-c) for Case A and (d-f) for Case B).

Figure 9 .
Figure 9.Comparison of the Young's Modulus in Case A (a) and Case B (b) with the same CG.

Figure 10
Figure 10 shows the comparison of NGH reservoir vertical displacement in both cases at the same CG.The positive value in the figure legend indicates the downward settlement of the formation and the negative value indicates the uplift of the formation.It was observed in Figure 10 that the maximum subsidence occurs at the interface of the overburden and NGH layer near the production well.The value of maximum subsidence is 1.06 m in Case B compared with that of 1.08 m in Case A. The maximum uplift occurs

Figure 9 .
Figure 9.Comparison of the Young's Modulus in Case A (a) and Case B (b) with the same C G .

Figure 10
Figure 10 shows the comparison of NGH reservoir vertical displacement in both cases at the same C G .The positive value in the figure legend indicates the downward

Energies 2022 , 17 Figure 10 .
Figure 10.Comparison of the vertical displacement in Case A (a) and Case B (b) with the same CG.
(a) A coupled THMC numerical model was developed in this study for the description of the complex phase change behavior, the gas-water two-phase flow, heat transfer, and formation deformation associated with NGH field-scale production; (b) The simulation results suggested that the injection of CO2 promotes the dissociation of NGH compared with depressurization only.The cumulative production of gas and water increased by 27.88% and 2.90%, respectively based on 2000 days of production simulation; (c) In Case A (depressurization only), NGH dissociation is symmetric and propagates outward based on the center of production well.While in Case B (depressurization coupled with CO2 exchange), the majority of NGH dissociation occurs in the region between the injection well and the production well; (d) NGH reservoir subsidence propagates outward based on the center of the production well in both cases.The maximum subsidence occurs at the interface of the overburden and NGH layer near the production well.The maximum uplift occurs at the interface of the underburden and NGH layer at the production well; (e) The subsidence of the NGH reservoir is less in Case B (depressurization coupled with CO2 exchange) compared with that of Case A (depressurization only) at the same cumulative gas production.The results suggest CO2 exchange in the NGH reservoirs alleviates the issue of reservoir subsidence during production and maintains good reservoir stability.The results of this study can provide a new insight for high-efficient and safe gas extraction from marine deposits of NGH.Author Contributions: Writing-original draft preparation, Q.L.; writing-review and editing, S.L.; investigation, S.D.; validation, Z.Y.; visualization, L.L. and S.L.All authors have read and agreed to the published version of the manuscript.

Figure 10 .
Figure 10. of the vertical displacement in Case A (a) and Case B (b) with the same C G .
•m −3 ; H s , H h , H g , H w represents the enthalpy of rock, hydrate, methane, and water, respectively, J•kg −1 ; λ e is the effective thermalconductivity, W•m −1 K −1 ; n h is the rate of NGH dissociation per unit volume, mol•m −3 •s −1 ; ∆H h is the dissociation heat of NGH per mole of NGH, J•mol −1 ; Q represents the rate of heat injection from the surrounding per unit volume, J•m −3 •s −1 .The effective thermal conductivity of the hydrate-bearing sediments λ e is used in Equation (12): λ e = (1 − φ)λ s + φS h λ h + φS g λ g + φS w λ w

Table 1 .
The main parameters used in the numerical simulation.Schematic diagram of the geological model.

Table 1 .
The main parameters used in the numerical simulation.