Investigating the Interaction Effects between Reservoir Deformation and Hydrate Dissociation in Hydrate-Bearing Sediment by Depressurization Method

the Interaction Effects between Reservoir Deformation Hydrate Hydrate-Bearing Abstract: Natural gas hydrate (NGH) has been widely focused on having great potential for alter-native energy. Numerous studies on gas production from hydrate-bearing sediments have been conducted in both laboratory and ﬁeld. Since the strength of hydrate-bearing sediments depends on the saturation of NGH, the decomposition of NGH may cause the failure of sediments, then leading to reservoir deformation and other geological hazards. Plenty of research has shown that the reservoir deformation caused by hydrate decomposition is considerable. In order to investigate this, the inﬂuence of sediment deformation on the production of NGH, a fully coupled thermo-hydro-chemo-mechanical (THMC) model is established in this study. The interaction effects between reservoir deformation and hydrate dissociation are discussed by comparing the simulation results of the mechanical coupling and uncoupled models on the laboratory scale. Results show that obvious differences in behaviors between gas and water production are observed among these two models. Compared to the mechanical uncoupled model, the mechanical coupling model shows a signiﬁcant compaction process when given a load equal to the initial pore pressure, which leads to a remarkable decrease of effective porosity and reservoir permeability, then delays the pore pressure drop rate and reduces the maximum gas production rate. It takes a longer time for gas production in the mechanical coupling model. Since the reservoir temperature is impacted by the comprehensive effects of the heat transfer from the boundary and the heat consumption of hydrate decomposition, the reduced maximum gas production rate and extended gas production process for the mechanical coupling model lead to the minimum reservoir temperature in the mechanical coupling model larger than that of the mechanical uncoupled model. The reduction of the effective porosity for the mechanical coupling model causes a larger cumulative water production. The results of this paper indicate that the reservoir deformation in the gas production process should be taken into account by laboratory and numerical methods to accurately predict the behaviors of gas production on the ﬁeld scale.


Introduction
Natural gas hydrate (NGH) is an ice-like crystal in which gas molecules, most often methane, are trapped by water molecules. Naturally existing in marine and permafrost, NGH is always found at low temperatures and high pressures along with gas and water, simultaneously [1,2]. Its properties and a large number of reserves in nature make it an attractive future energy source as well as a potential hazard [3][4][5]. Laboratory and short-term field studies have shown that gas production from hydrate-bearing sediment is technically feasible using depressurization [6,7], but there are still challenges to produce gas from hydrate-bearing sediment safely, efficiently, and economically due to the complex coupled thermo-hydro-chemo-mechanical (THMC) process [8]. Among them, the ground deformation induced by gas hydrate dissociation has been widely focused on for its risks of wellbore stability, gas leakage and submarine landslide [9,10].
A vast number of triaxial test results based on synthetic and nature hydrate samples have indicated that the mechanic properties of hydrate-bearing sediment are closely related to NGH saturation and distribution morphology [11][12][13][14][15][16][17][18][19]. With increasing NGH saturation, the bonding strength would be significantly strengthened, leading to the improvement of hydrate-bearing sediment strength as well as the dilation of hydrate-bearing sediment samples [20,21]. Therefore, the decomposition of NGH would undoubtedly weaken the strength of the hydrate-bearing sediment [22]. Changes in in situ geotechnical parameters of hydrate-bearing sediments may cause the heterogeneous deformation of the seabed with constant upper loads [23,24]. Meanwhile, pore pressure variations induced by the changes of effective stress state during hydrate dissociation could also adversely affect seabed stability. [25]. These pore pressure changes caused by NGH dissociation vary according to the initial effective permeability of hydrate-bearing sediments, and the overpressure is often generated in low permeability layers with rapid dissociation of NGH, thereby decreasing the effective pressure and shear strength of unconsolidated sediments [26,27].
In the limited number of field NGH production tests, the seafloor stability during methane production was closely monitored. Series of marine geophysical equipment based on seafloor environment monitor systems were developed, including the 3-component servo-accelerometer system and the tiltmeters and pressure sensors based on seafloor deformation monitoring system developed by Japan [28,29], the Gas Hydrate Sea-Floor Observatory developed by the USA [30], the three-dimensional environmental monitoring system and "four-in-one" environmental monitoring system developed by China [31,32]. The advantages of quick response, high-efficiency, and reflecting the real situation of sediment deformation accompanied by NGH production are obvious by these systems. However, as a time and money cost program, the field gas production test results obtained by these systems may be influenced by the normal background changes in the seafloor [33]. Therefore, it is still difficult to understand the mechanism of sediment deformation accompanied by hydrate dissociation for safety production in field tests.
Numerical simulation is considered as a feasible method to predict production characteristics of gas and water, reservoir temperature and pressure distribution, and the mechanical behavior of the reservoir. TOUGH+HYDRATE is a typical successful simulator for the simulation of system behavior in hydrate-bearing stratum by coupling of hydraulic, thermal, and thermodynamic governing equations [34]. It has been widely used to simulate the gas production behavior in hydrate-bearing sediment by different production methods from laboratory to reservoir scale [35][36][37][38]. It can also predict geomechanical responses during hydrate production by linking FLAC3D through a coupling code [25]. Implementing the governing equations related to hydraulic, thermal, and thermodynamic behavior in hydrate-bearing sediment directly into FLAC3D is also a practical solution [39]. Other simulators, including Hydrsim simulator, MH21 code, STOMP-HYD, CMG-STARS, are also developed to simulate the behavior of gas production from hydrate-bearing sediment [40][41][42][43]. The most recent report conducted by White et al. [44] has compared mainstream NGH production simulation codes.
The purpose of this work is to study the influence of sediment deformation on the gas production from hydrate-bearing sediment by depressurization method with Comsol software that incorporates hydraulic, thermal, thermodynamic, and mechanics. The interaction effects between reservoir deformation and hydrate dissociation are considered in Energies 2021, 14, 548 3 of 16 the gas production process by coupling the formation strain into the two-phase seepage governing equations. Although the present work is based on a laboratory-scale model, the methodology behind this work could also apply for field-scale modeling with detailed stratum parameters.

Models for Two-Phase Flow in Porous Media
The two-phase flow model based on Darcy's theory is used to describe the methane and water flow in the porous media [45].
where ρ i , P i , v i , µ i are density, pressure, velocity, and viscosity, respectively. The subscripts are non-wetting and wetting phases. g is the gravitational acceleration. ε v is the volumetric strain of sediment. In order to establish a two-phase Darcy equation, Equation (1) is converted into pressure form.
The dependent variables of the above equations are P w and P g , respectively. It is a function related to pore pressure and saturation and can be defined as: where n and n wg are porosity and effective porosity, respectively. n 0 is initial porosity. ε v is the volumetric strain of sediment.P c is pore pressure.S w and S h are wetting phase saturation and hydrate saturation. P g and P w are the pressure of the non-wetting phase and wetting phase, respectively. In Equation (5), S is the storage coefficient that is defined as: where χ f is the gas compression coefficient. According to the gas state equation, the relationship between the density of methane gas, temperature, and pressure can be expressed as [46]: where P g,c and T c are the critical pressure and critical temperature, and the corresponding values are 4640 kPa and 190.8 K, respectively. The coefficient K h and k r,i in Equation (6) are absolute permeability and relative permeability, respectively. The mathematical formula describing the absolute permeability of hydrate deposits is generally related to the hydrate saturation. In this paper, the mathematical model proposed by Amyx is used to describe the absolute permeability of hydrate sediments [47]: K f × 10 8 × n wg 9.147 , n wg ≥ 0.11 (15) where K i is the initial permeability, K f is the permeability without hydrate. This mathematical model has been successfully used by Yousif et al. to describe the permeability of gas hydrate sediments during decomposition [48]. Relative permeability can be described and calculated by the model set up by Stone [49]: Hydrate saturation, non-wetting phase saturation and wetting phase saturation meet the following conditions:

Kinetic Model of Hydrate Decomposition
NGH is composed of water molecules and methane molecules. The expression of decomposition and synthesis reaction is as follows: Based on the Kim-Bishnoi empirical formulation [50], the production rate of methane gas was calculated by Formula (20).
Then the water production rate and the decomposition and synthesis rate of methane can be calculated as: .
where M g , M w , M h are the molar mass of methane, the molar mass of water and the molar mass of hydrate, respectively. N h is the hydration number, with a value of 5.75. K d is the decomposition constant of NGH; A s is the reaction area of NGH particles [51], which is described by the following empirical formula: The value in Equation (23) is the natural gas hydrate decomposition rate constant, which is generally taken as 3.6 × 10 4 mol/ m 2 · Pa · s [52]; R is the ideal gas constant; T is the temperature, K; ∆E a is the kinetic energy, J; ∆E a /R = 9753.73K. A geo is the surface area ratio of hydrate particles per unit volume, assuming that the natural gas hydrate particles are spherical particles, the calculation can be obtained as A geo is 7.5e5 m −1 . P e is the three-phase equilibrium pressure of natural gas hydrate, and the empirical formula P e = 1.15 exp(49.3185 − 9459 T ) is used [53].

Energy Conservation Equation
The decomposition of natural gas hydrate sediments is a process of absorbing heat from the surrounding environment, and the heat exchange process can be controlled by heat conduction and heat convection equations.
(ρC) eq ∂T ∂t The dependent variable of the energy governing equation is temperature T. C g , C w , C h , C s are the specific heat capacities of methane non-wetting phase, wetting phase, NGH, and skeleton particles, respectively. λ g , λ w , λ h , λ s are thermal conductivity of methane non-wetting phase, wetting phase, NGH, and skeleton particles, respectively. ∆H D is the enthalpy of the methane decomposition process. v g , v w are Darcy velocities of gas and water. σ g is the Joule-Thomson throttling coefficient [54].

Mechanical Models
The dependent variable of the mechanical equilibrium equation is displacement u, and the main governing equations are described as follows: where ρ c is the equivalent density of hydrate-bearing sediment, and its calculation formula is ρ c = (1 − n)ρ s +nS h ρ h , n is the porosity. ρ s , ρ h are the density of skeleton particles and hydrate, respectively. S is the stress tensor, which can be expressed as: where S 0 is the initial stress tensor. C is the elastic coefficient matrix. Assuming that the soil mass model is isotropic, the elastic modulus is defined as a function related to the gas hydrate saturation E = (300 + 1350S h ), MPa. ε is the total strain tensor, which can be expressed by displacement gradient. The calculation formula of the displacement gradient is as follows: According to the incremental theory: The ε e ij and ε p ij are elastic and the plastic strain, respectively. The plastic strain increment is calculated using the plastic potential energy function [55]: The potential Q is written in terms of at most three invariants of Cauchy's stress tensor.
are the invariants of the stress tensor. λ is a scalar that is determined from the consistency condition to keep the stress point on the yield surface, f (S ) = 0.
In this paper, the Mohr-Coulomb yield criterion is used to judge the yield conditions [56]: where S 1 and S 3 are the maximum and the minimum principal stress, respectively. N ϕ = (1 + sin ϕ )/(1 − sin ϕ ); ϕ is the internal friction angle; c is the cohesion, which is a function related to the NGH saturation and expressed as C= 3.3S h × 10 6 Pa [57].

Model Description
Since Masuda's water-bearing sandstone core decompression experiment has been widely used to verify numerical models [58], our simulation model is also verified by this experimental data without considering the mechanical model. After the verification, the mechanical model is added to the former model with the same geometric model. An axisymmetric hydrate-bearing sediment model with a diameter of 51 cm and a length of 30 cm is established. The initial and boundary conditions are the same as the experiment conducted by Masuda. As shown in Figure 1, six points numbered C1 to C6 are fixed at the symmetry axis to monitor the system pressure, and the other six points numbered B1 to B6 are located at the top of the model to monitor the system temperature. Both temperature and pressure are monitored at points named T1 to T6, which lie at the bottom of the model.

Boundary Condition and Initial Values
The initial parameters used in this study are similar to the experiment conducted by Masuda and listed in Table 1. Table 1. Parameters adopted in this study [54,[58][59][60].

Parameters
Value Parameters Value

Boundary Condition and Initial Values
The initial parameters used in this study are similar to the experiment conducted by Masuda and listed in Table 1. Table 1. Parameters adopted in this study [54,[58][59][60]. P c can be calculated by Equation (10) and the initial saturation of the water. The mass flux on the cylindrical surface and the lower end of the core is set as zero, and the upper end is defined as the pressure outlet (2.99 MPa). All the walls are set with a constant temperature (275.25 K) since the reactor is put into a constant temperature bath in the experiment. Then the pressure and heat transfer boundary are given by:

Parameters
where λ heat is the heat transfer coefficient of the boundary, and T 0 is ambient temperature.
In the sediment deformation considered model, the boundary conditions of the added mechanics module are the same as those of the conventional lateral confined triaxial compression test. A fixed constraint is applied to the lower boundary of the model. A constant axial pressure (3.75 MPa), which is equal to the initial pore pressure, is applied to the pressure outlet boundary, and lateral constraint conditions are applied to the cylinder:

Comparison of Simulation and Experimental Results
Figures 2-4 are comparative results of pore pressure, system temperature, accumulated production gas for simulation with Masuda's experiment. Figure 2 shows the pore pressure changes with the production time at the far end of the pressure outlet without the consideration of sediment deformation. It is obvious that the pressure drop processes are very similar to simulation and experiment. In general, as the mining progresses, the pressure gradually decreases until it reaches the outlet pressure. Another worth noting phenomenon is an obvious slowdown in the pressure drop at about 40 min due to the fact that the methane gas released from the NGH decomposition buffers the pressure drop speed. are very similar to simulation and experiment. In general, as the mining progresses, the pressure gradually decreases until it reaches the outlet pressure. Another worth noting phenomenon is an obvious slowdown in the pressure drop at about 40 min due to the fact that the methane gas released from the NGH decomposition buffers the pressure drop speed.   Figure 3 shows the temperature change of the core at point T3. It can be seen from the figure that the results of the numerical analysis are basically consistent with the trends of the experimental results. The decomposition of natural gas hydrate is an endothermic process. As the mining progresses, the hydrate gradually decomposes and absorbs heat in the model so that the temperature at point T3 gradually decreases. At the same time, the ambient temperature at the wall boundary transfers heat to the inside of the reactor. When the ambient temperature can provide enough heat for the dissociation reaction of NGH, the temperature starts to increase until it reaches the ambient temperature of 275.25 K. Therefore, the temperature curve shows a downward trend followed by an upward trend. Figure 4 shows the accumulated production gas of simulation and experiment. It shows that in the initial stage of decomposition, the gas production rate of the simulation result is slightly lower than the experimental measurement value. After 50 min, the gas production rate of the simulation result gradually increases, and the total gas production of the final simulation result is basically equal to the total gas production of the experiment result. The comparison has shown a consistent result of simulation and Masuda's experiment, which indicates that the simulation model used in this study is validated.
Energies 2021, 14, x FOR PEER REVIEW 10 of 19 Figure 3 shows the temperature change of the core at point T3. It can be seen from the figure that the results of the numerical analysis are basically consistent with the trends of the experimental results. The decomposition of natural gas hydrate is an endothermic process. As the mining progresses, the hydrate gradually decomposes and absorbs heat in the model so that the temperature at point T3 gradually decreases. At the same time, the ambient temperature at the wall boundary transfers heat to the inside of the reactor. When the ambient temperature can provide enough heat for the dissociation reaction of NGH, the temperature starts to increase until it reaches the ambient temperature of 275.25 K. Therefore, the temperature curve shows a downward trend followed by an upward trend.  Figure 4 shows the accumulated production gas of simulation and experiment. It shows that in the initial stage of decomposition, the gas production rate of the simulation result is slightly lower than the experimental measurement value. After 50 min, the gas production rate of the simulation result gradually increases, and the total gas production

Interaction Effects between Reservoir Deformation and Hydrate Dissociation
When performing the depressurization process, with fluids discharging out of the reactor, the pore pressure gradually decreases, resulting in the increment of effective stress for the skeleton particles. Because the cylindrical boundary of the model is set with no lateral displacement, the sediment can only be compressed in the axial direction. Generally, with the gradual compression of skeleton particles caused by hydrate decomposition, the internal pore structure will be changed, and the absolute porosity of the sample will be reduced, thereby controlling the NGH decomposition indirectly. Figure 5a shows the cumulative gas and gas production rate variation with production time. Since there are free gas and water in the initial status of both the mechanical uncoupled model and mechanical coupling model, the gas production rate shows an increment at the beginning of production, then reduces for the rapid consumption of free gas. The system pressure decreases gradually with the discharge of free gas and water. When it reaches the NGH equilibrium pressure, the NGH begins to decompose to supplement the gas source. Then the gas production rate keeps a constant growth over a long production time with nearly 80 min by the mechanical coupling model and 110 min by the mechanical uncoupled model. A significant difference in gas production rate is shown for these two models; the growth and decline magnitude of gas production rate for the loading model are relatively small compared to that of the mechanical uncoupled model, which indicates that the sediment deformation would restrain the decomposition of NGH and leads to the extension of production time. Figure 5b shows the cumulative water and water production rate variation with production time. The water production rate of the mechanical uncoupled model varies with production time rapidly, while the water production rate for the mechanical coupling model keeps relatively smooth in the long term of production time. Excepting that the water production rate of the mechanical uncoupled model has a little advantage between 50 and 150 min, the water production rate and the cumulative water for the mechanical coupling model are significantly high.

Interaction Effects between Reservoir Deformation and Hydrate Dissociation
When performing the depressurization process, with fluids discharging out of the reactor, the pore pressure gradually decreases, resulting in the increment of effective stress for the skeleton particles. Because the cylindrical boundary of the model is set with no lateral displacement, the sediment can only be compressed in the axial direction. Generally, with the gradual compression of skeleton particles caused by hydrate decomposition, the internal pore structure will be changed, and the absolute porosity of the sample will be reduced, thereby controlling the NGH decomposition indirectly. Figure 5a shows the cumulative gas and gas production rate variation with production time. Since there are free gas and water in the initial status of both the mechanical uncoupled model and mechanical coupling model, the gas production rate shows an increment at the beginning of production, then reduces for the rapid consumption of free gas. The system pressure decreases gradually with the discharge of free gas and water. When it reaches the NGH equilibrium pressure, the NGH begins to decompose to supplement the gas source. Then the gas production rate keeps a constant growth over a long production time with nearly 80 min by the mechanical coupling model and 110 min by the mechanical uncoupled model. A significant difference in gas production rate is shown for these two models; the growth and decline magnitude of gas production rate for the loading model are relatively small compared to that of the mechanical uncoupled model, which indicates that the sediment deformation would restrain the decomposition of NGH and leads to the extension of production time. Figure 5b shows the cumulative water and water production rate variation with production time. The water production rate of the mechanical uncoupled model varies with production time rapidly, while the water production rate for the mechanical coupling model keeps relatively smooth in the long term of production time. Excepting that the water production rate of the mechanical uncoupled model has a little advantage between 50 and 150 min, the water production rate and the cumulative water for the mechanical coupling model are significantly high.  Figure 6a,b shows the pore pressure variation with production time for the mechanical coupling and uncoupled models at C1 to C6, respectively. The pore pressure of C1 is equal to outlet pressure since it locates at the pressure outlet boundary. For the initial stage of the gas production, the free gas and water are discharged out of the reactor rapidly, leading to the pore pressure dropping significantly faster than other stages for both mechanical coupling and uncoupled models. There is an obvious pore pressure rebound point after the rapid dropping stage, which is caused by the balance of NGH decomposition and fluids discharge. The pore pressure rebound point for the mechanical uncoupled model appears about 15 min after depressurization, about 10 min earlier than that of the loading model. Meanwhile, the pore pressure at the pressure rebound point is related to the location of the monitory point for the same model. The further the monitory point to the pressure outlet and the higher the pore pressure at the pressure rebound point are, the longer the rebound duration is. Compared to the mechanical uncoupled model, the pore pressure of all monitory points at the pressure rebound points is higher, and the pore pressure drop rate is significantly lower in the mechanical coupling model. It costs about 300 and 400 min for pore pressure to drop to the outlet pressure at the bottom of the cylinder for the mechanical uncoupled and coupling models, respectively. It is shown that the sediment deformation may promote the occurrence of overpressure in sediment pores.  Figure 6a,b shows the pore pressure variation with production time for the mechanical coupling and uncoupled models at C1 to C6, respectively. The pore pressure of C1 is equal to outlet pressure since it locates at the pressure outlet boundary. For the initial stage of the gas production, the free gas and water are discharged out of the reactor rapidly, leading to the pore pressure dropping significantly faster than other stages for both mechanical coupling and uncoupled models. There is an obvious pore pressure rebound point after the rapid dropping stage, which is caused by the balance of NGH decomposition and fluids discharge. The pore pressure rebound point for the mechanical uncoupled model appears about 15 min after depressurization, about 10 min earlier than that of the loading model. Meanwhile, the pore pressure at the pressure rebound point is related to the location of the monitory point for the same model. The further the monitory point to the pressure outlet and the higher the pore pressure at the pressure rebound point are, the longer the rebound duration is. Compared to the mechanical uncoupled model, the pore pressure of all monitory points at the pressure rebound points is higher, and the pore pressure drop rate is significantly lower in the mechanical coupling model. It costs about 300 and 400 min for pore pressure to drop to the outlet pressure at the bottom of the cylinder for the mechanical uncoupled and coupling models, respectively. It is shown that the sediment deformation may promote the occurrence of overpressure in sediment pores.  Figure 7 shows the temperature variation with a production time at points T1 to T6 for the mechanical coupling and uncoupled models at monitoring points T1 to T6. With  Figure 7 shows the temperature variation with a production time at points T1 to T6 for the mechanical coupling and uncoupled models at monitoring points T1 to T6. With the consumption of hydrate-bearing sediment sensible heat by NGH decomposition, the reservoir temperature reduces rapidly. Since the gas production rate for the mechanical uncoupled model is significantly larger than that of the mechanical coupling model at the initial production stage, the system temperature reaches the minimum temperature of 274 K at a production time of 150 min for the mechanical uncoupled model, which is ahead of 150 min than that of the mechanical coupling model with the lowest temperature of 274.3 K. With the ambient heat transfers into the model, there is a slow rise in system temperature until it reaches to the ambient temperature. It indicates that the gas production process is significantly influenced by ambient temperature. Although this influence exists widely in NGH production tests at both the laboratory and field scale, it is more dramatic for small scale models. However, most gas production tests on a laboratory scale neglect this influence, leading to the overestimation of gas production effects.

System Temperature Changing Characteristics
(a) (b) Figure 6. Pore pressure variation with production time for the mechanical (a) uncoupled and (b) coupled model. Figure 7 shows the temperature variation with a production time at points T1 to T6 for the mechanical coupling and uncoupled models at monitoring points T1 to T6. With the consumption of hydrate-bearing sediment sensible heat by NGH decomposition, the reservoir temperature reduces rapidly. Since the gas production rate for the mechanical uncoupled model is significantly larger than that of the mechanical coupling model at the initial production stage, the system temperature reaches the minimum temperature of 274 K at a production time of 150 min for the mechanical uncoupled model, which is ahead of 150 min than that of the mechanical coupling model with the lowest temperature of 274.3 K. With the ambient heat transfers into the model, there is a slow rise in system temperature until it reaches to the ambient temperature. It indicates that the gas production process is significantly influenced by ambient temperature. Although this influence exists widely in NGH production tests at both the laboratory and field scale, it is more dramatic for small scale models. However, most gas production tests on a laboratory scale neglect this influence, leading to the overestimation of gas production effects.   Figure 8 shows NGH decomposition rate variation with production time for the mechanical Figure 8a uncoupled and Figure 8b coupled models at monitoring points B1 to B6 and T1 to T6, B1 to B6 corresponding to the top of the model and T1 to T6 corresponding to the bottom of the model. B1 to B6, which are closer to the pressure outlet, have a high NGH decomposition rate for the rapid drop of pore pressure near the outlet. Since the equilibrium pressure of NGH is temperature-dependent, as mentioned former, with the decrement of model temperature caused by sensible heat consumption, the equilibrium pressure of NGH also drops, which leads to the slowdown of NGH decomposition rate for both mechanical coupling and uncoupled models. In the later stage, the model temperature rises for heat absorption from around the environment, giving rise to the improvement of NGH equilibrium pressure and corresponding to the NGH decomposition rate increment.

NGH Decomposition Rate Characteristics
the equilibrium pressure of NGH is temperature-dependent, as mentioned former, with the decrement of model temperature caused by sensible heat consumption, the equilibrium pressure of NGH also drops, which leads to the slowdown of NGH decomposition rate for both mechanical coupling and uncoupled models. In the later stage, the model temperature rises for heat absorption from around the environment, giving rise to the improvement of NGH equilibrium pressure and corresponding to the NGH decomposition rate increment.  Figure 9 shows the effective porosity variation with production time for the mechanical (a) uncoupled and (b) coupled models. Generally, the effective porosity of hydratebearing sediment is considered as a hydrate saturation related parameter, which is increased with the reduction of hydrate saturation, as shown in Equation (8). With the help of ambient heat, the NGH in the mechanical coupled model decomposes completely,  Figure 9 shows the effective porosity variation with production time for the mechanical (a) uncoupled and (b) coupled models. Generally, the effective porosity of hydrate-bearing sediment is considered as a hydrate saturation related parameter, which is increased with the reduction of hydrate saturation, as shown in Equation (8). With the help of ambient heat, the NGH in the mechanical coupled model decomposes completely, which leads to the effective porosity reaching the initial defined absolute porosity of 0.18. The sediment deformation causes a significant decrease of the effective porosity, as shown in Figure 9b, and the final effective porosity of the mechanical coupling model is 0.174.

Effective Porosity and Absolute Permeability Changing Characteristics
The deep-rooted reasons of sediment deformation restrained on the gas production of hydrate-bearing sediment can be obtained from the information contained in Figure 9, which shows the absolute permeability variation with production time for the mechanical (a) uncoupled and (b) coupled models. As shown in Figure 10, even the initial porosity is the same with the mechanical coupling and uncoupled models, the final absolute permeability of the mechanical coupling model is reduced nearly 30% from 8.03 × 10 −14 m 2 to 5.57 × 10 −14 m 2 . Meanwhile, the growth of absolute permeability is lagged significantly in the mechanical coupling model. The pressure drop rate is controlled by the discharging rate of free gas and water, and the lowered absolute permeability restrains the discharge of decomposed gas and water, thus limiting the pressure drop rate and NGH decomposition. which leads to the effective porosity reaching the initial defined absolute porosity of 0.18. The sediment deformation causes a significant decrease of the effective porosity, as shown in Figure 9b, and the final effective porosity of the mechanical coupling model is 0.174.
(a) (b) Figure 9. Effective porosity variation with production time for the mechanical (a) uncoupled and (b) coupled models.
The deep-rooted reasons of sediment deformation restrained on the gas production of hydrate-bearing sediment can be obtained from the information contained in Figure 9, which shows the absolute permeability variation with production time for the mechanical (a) uncoupled and (b) coupled models. As shown in Figure 10, even the initial porosity is the same with the mechanical coupling and uncoupled models, the final absolute permeability of the mechanical coupling model is reduced nearly 30% from 8.03 × 10 −14 m 2 to 5.57 × 10 −14 m 2 . Meanwhile, the growth of absolute permeability is lagged significantly in the mechanical coupling model. The pressure drop rate is controlled by the discharging rate of free gas and water, and the lowered absolute permeability restrains the discharge of decomposed gas and water, thus limiting the pressure drop rate and NGH decomposition. As the field, NGH production test conducted at Nankai Trough by Japan in 2013 showed that there is a significant difference between the numerical predicting and field test results, which is attributed to the heterogeneity distribution of NGH. However, the comparative results of the NGH production characteristics for the mechanical coupling and uncoupled models show that there is a significant influence of sediment deformation  The deep-rooted reasons of sediment deformation restrained on the gas production of hydrate-bearing sediment can be obtained from the information contained in Figure 9, which shows the absolute permeability variation with production time for the mechanical (a) uncoupled and (b) coupled models. As shown in Figure 10, even the initial porosity is the same with the mechanical coupling and uncoupled models, the final absolute permeability of the mechanical coupling model is reduced nearly 30% from 8.03 × 10 −14 m 2 to 5.57 × 10 −14 m 2 . Meanwhile, the growth of absolute permeability is lagged significantly in the mechanical coupling model. The pressure drop rate is controlled by the discharging rate of free gas and water, and the lowered absolute permeability restrains the discharge of decomposed gas and water, thus limiting the pressure drop rate and NGH decomposition. As the field, NGH production test conducted at Nankai Trough by Japan in 2013 showed that there is a significant difference between the numerical predicting and field test results, which is attributed to the heterogeneity distribution of NGH. However, the comparative results of the NGH production characteristics for the mechanical coupling and uncoupled models show that there is a significant influence of sediment deformation As the field, NGH production test conducted at Nankai Trough by Japan in 2013 showed that there is a significant difference between the numerical predicting and field test results, which is attributed to the heterogeneity distribution of NGH. However, the comparative results of the NGH production characteristics for the mechanical coupling and uncoupled models show that there is a significant influence of sediment deformation on the production behavior of NGH in hydrate-bearing sediment. Therefore, the lack of accurate prediction of the sediment deformation may be one of the reasons causing the difference between numerical predicting and field test results, and it is necessary to consider the sediment deformation when employing numerical simulation to predict NGH production test at the field scale to ensure the prediction accuracy.

Conclusions
This paper has discussed the interaction effects between reservoir deformation and hydrate dissociation in hydrate-bearing sediment by the depressurization method. A fully coupled thermo-hydro-chemo-mechanical (THCM) model is established to simulate the process of gas production from hydrate-bearing sediment. The gas and water production behavior, pore pressure change characteristics, effective porosity and absolute permeability change characteristics for the mechanical coupling and uncoupled models with the laboratory scale are compared. Results show that obvious differences between production Energies 2021, 14, 548 14 of 16 behaviors of gas and water in the two models have been observed. Comparing to the mechanical uncoupled model, the mechanical coupling model shows a significant compaction process when given a load equal to initial pore pressure, which leads to a remarkable decrease of effective porosity and reservoir permeability, then delays the pore pressure drop rate and reduces the maximum gas production rate. Gas production takes a longer time in the mechanical coupling model. Since the reservoir temperature is influenced by the comprehensive effect of heat transfer from the boundary and heat consumption for hydrate decomposition, the reduced maximum gas production rate and extended gas production duration lead to the minimum reservoir temperature in the mechanical coupling model, which is larger than that of the mechanical uncoupled model. The reduction of effective porosity for the mechanical coupling model causes a larger cumulative water production. Since reservoir deformation has not been considered by most laboratory apparatus in evaluation studies of production methods, it is difficult to make an accurate conclusion for these methods. Therefore, the reservoir deformation should be taken into account in the gas production process by laboratory and numerical methods to predict gas production behaviors at the field scale. Data Availability Statement: All Data has been provided in the paper.