THMC Fully Coupled Model of Natural Gas Hydrate under Damage Effect and Parameter Sensitivity Analysis

Abstract

In order to study the influence of damage on the gas production of natural gas hydrate, a multi-physical field theoretical model considering damage effect and coupling thermal-hydraulic-mechanical-chemical (THMC) was established by theoretical analysis and numerical simulation. The THMC model establishes the relationship between the elastic modulus of hydrate sediment and hydrate saturation during the whole process of hydrate decomposition. The THC (thermal-hydraulic-chemical) and THMC fully coupled models not considering or considering the damage effect were compared and analyzed, and the reliability of the THMC fully coupled model was verified. On this basis, the deformation, permeability and damage of hydrate sediments under different initial hydrate saturations and different depressurization amplitudes, as well as the hydrate gas production rate and cumulative gas production, are analyzed. The results showed that higher initial hydrate saturation inhibited the development of damage, maintained stable gas production and increased cumulative gas production. Larger depressurization promoted damage and increased cumulative gas production, but it was easy to cause stability problems.

1. Introduction

Due to population growth, industrial expansion and economic growth, non-renewable fossil fuels cannot meet the growing global energy demand. Therefore, energy diversification is essential and necessary [1]. As a fossil fuel, natural gas is a strong greenhouse gas. Combustion and leakage will cause the greenhouse effect [2]. The release of methane from hydrate-bearing sediments is considered to be slow and chronic, rather than catastrophic [3]. Natural gas hydrate is recognized as a new type of clean energy with broad prospects because of its large reserves, high calorific value and pollution-free products [4]. Unlike the development of conventional oil and gas resources, the development of natural gas hydrates releases natural gas by destroying its thermodynamic equilibrium [5], and the decomposition of natural gas hydrates is a complex process involving thermal-hydraulic-mechanical-chemical (THMC) multi-physical field coupling. The mechanical behavior of hydrate in the process of decomposition is complex [6,7] and will produce mechanical damage, which will affect the stability of the sediments [8,9]. The fracture damage characteristics have an important influence on the strength. Ignoring the decomposition damage effect of hydrates will lead to deviations in the characterization of the natural gas hydrate decomposition process and the prediction of gas production.

Many researchers have carried out a series of studies on the damage effect of gas hydrate sediments. Zhu et al. [10] considered that damage would have a non-negligible effect on the seepage field, the temperature field and the stress field, and proposed the damage evolution equation for rock porous media. Yang et al. [11] considered that the mechanical properties of sediments were the result of the combined effect of hydrate cementation and compaction and established an elastic–plastic damage constitutive model of hydrate sediments that can reflect the volume change characteristics and the influence of confining pressure. Shen et al. [12] established a generalized plastic and elastic damage combination model of nonlinear damage, and studied the effects of hydrate saturation, confining the pressure and void ratio on the stress–strain and volume deformation of sediments. Based on the basic assumption that the strength of hydrate-bearing sediments obeys the Drucker-Prager criterion and Weibull distribution, Li et al. [13] established a damage statistical constitutive model that can simultaneously describe the strain-softening law and strain-hardening law of hydrate-bearing sediments by combining the statistical strength theory with the statistical damage theory. Yan et al. [14] established a statistical damage constitutive model of hydrate-bearing sediments considering the influence of the hydrate occurrence mode and a damage constitutive model considering the influence of temperature and pore pressure and described the evolution of damage factors. Yan et al. [15] introduced a damage factor into the strength of hydrate sediments based on the modified Duncan–Chang model and analyzed the influence of hydrate saturation on strain softening, stiffness and strength. Zhu et al. [16] established an equivalent variable elastic modulus damage constitutive model of natural gas hydrate sediments based on the micromechanical mixing rate theory of composite materials and the rock pore damage theory, considering the influence of hydrate content, confining pressure and internal pore changes on the mechanical properties of natural gas hydrate sediments, combined with the statistical distribution law of micro-element strength of hydrate sediments. Liu et al. [17] established a modified damage constitutive model based on the material parameter equation obtained by discrete element simulation under the Drucker–Prager criterion and triaxial test of hydrate sediments. This model can predict the mechanical behavior of sediments under high confining pressure more accurately. In the above studies on the damage effect of natural gas hydrates, only the mechanical behavior of hydrate sediments was studied, and the coupling analysis of the damage effect and multi-physical field relationship was not carried out.

Many researchers have carried out relevant research on the multi-field coupling problem of hydrate dissociation. Makogon [18] regarded hydrate decomposition as a single-phase gas flow problem, ignoring the decomposition of water production. On this basis, Verigin et al. [19] described the hydrate dissociation problem as gas-water two-phase seepage and considered the mass conservation of water, without considering the effect of temperature change. Then, the gas–water-hydrate three-phase mass conservation, the phase change decomposition of hydrate and the energy conservation in temperature change were gradually considered in the theoretical model [20]. Liang et al. [21] established a thermal-hydraulic-chemical (THC) coupling hydrate depressurization mining model and studied the influence of production pressure drop, initial absolute permeability, water saturation and phase equilibrium parameters on gas production under depressurization mining. Klar et al. [22] established a four-field multi-phase coupling model of thermal-hydraulic-mechanical-chemical (THMC) to simulate the stress-release phenomenon during hydrate decomposition. Sun [23] optimized the Klar A coupling model based on the energy dissipation theory. Wang [24] numerically simulated and analyzed the first trial production of the Shenhu hydrate reservoir in the South China Sea with the established thermal-hydraulic-mechanical-chemical (THMC) multi-field coupling model and studied the mechanical behavior and stability in the trial production engineering. Based on the continuous damage theory, Zhang et al. [25] introduced the three-parameter Weibull distribution and residual strength correction coefficient into the damage statistical constitutive model and established a damage statistical constitutive model of hydrate-bearing sediments considering the influence of damage threshold and residual strength. A multi-field coupling mathematical model of thermal–hydraulic–mechanical–chemical (THMC) considering the damage of hydrate-bearing sediments was established. Based on this model, the influence of structural damage of hydrate-bearing sediments on the deformation, pressure and temperature of sediment reservoirs during hydrate dissociation was discussed. Considering the influence of decomposition damage, load damage and residual strength, Huang et al. [26] established a thermal–hydraulic–mechanical–chemical (THMC) model coupled with damage constitutive and studied the variation of displacement, saturation, average pore pressure and temperature at different mining positions under horizontal well depressurization in permafrost regions. However, the establishment of the above damage constitutive model was mainly based on the Weibull distribution to establish the statistical damage relationship, and then the damage constitutive relationship of the sediment was established. The decomposition damage of hydrates was less considered, and the damage change of the whole process of hydrate decomposition and the influence of hydrate saturation on damage and strength were not fully reflected. At the same time, the damage effect in the process of hydrate decomposition and its influence on gas production needed to be further explored.

This paper holds that hydrate saturation is the main reason for the obvious change in the mechanical properties of sediments. Therefore, the damage relationship directly related to hydrate saturation and mechanical damage is adopted, and this relationship can reflect the complete change of the elastic modulus of sediments during hydrate decomposition. Based on this, a theoretical model of hydrate THMC full coupling considering damage is established, and the evolution of physical parameters of the natural gas hydrate model under different initial hydrate saturations and different depressurization amplitudes and its influence on natural gas production are studied. The research flow chart of this paper is shown in Figure 1.

2. THMC Full Coupling Theoretical Model of Natural Gas Hydrate

2.1. Analysis of THMC Full Coupling

The decomposition of natural gas hydrate is an endothermic reaction and produces methane gas and water, which affects the seepage process of gas and water, changes the temperature and bearing structure of porous media and affects the temperature field and the stress field. The seepage of gas and water will change the pore pressure and effective stress in the medium, make the natural gas hydrate under the phase equilibrium curve, and promote the further decomposition of the hydrate. At the same time, the seepage process is accompanied by the energy transfer process. The seepage of gas and water carries heat and changes the temperature of the porous medium. The temperature of porous media affects the phase equilibrium state and decomposition rate of hydrates, changes the density of percolating gas and produces thermal strain; the change of the stress state in porous media leads to the change of permeability, porosity and pore pressure, changes the phase equilibrium condition of hydrate decomposition, affects the seepage process and causes the change in thermal conductivity, specific heat capacity and other thermodynamic parameters. The coupling analysis of the THMC fully coupled model is shown in Figure 2.

2.2. Basic Assumptions

To establish a fully coupled THMC model, the following basic assumptions are proposed [27]:

• (1). The model hydrate structure is an I-type cubic crystal structure [28], containing 46 water molecules, composed of two small pores and six large pores;

• (2). The hydrate gas is pure methane, with no secondary gas, and the liquid is pure water, regardless of the solubility of methane in water;

• (3). Not considering the effect of temperature change on water density;

• (4). The seepage of methane and water conforms to Darcy’s law;

• (5). Methane and water are a saturated two-phase flow without considering the chemical interaction between sediment and fluid;

• (6). The temperature of sediment particles, hydrate, methane and water in each element of the numerical model is consistent. The hydrate follows the isotropic and is always elastic, and the sediment is elastic-plastic;

• (7). The secondary formation of hydrates in porous media is not considered.

Based on the above assumptions, the saturation formula is proposed;

(1)$S_i=\frac{V_i}{V_{\varphi }},i=g, w, h$

where $g$, $w$ and $h$ represent methane, water and gas-hydrate solids; $V_{\varphi }$ is sediment pore volume.

In this paper, sediments are regarded as skeleton, gas, water and hydrate particles that are filled in the pores of the skeleton, and porosity is defined as:

(2)$\varphi =\frac{V_g+V_w+V_h}{V+V_{\varphi }}$

Therefore, $S_g+S_w+S_h=1$; the porosity does not directly change with the saturation of the natural gas hydrate.

2.3. Decomposition Kinetic Equation of Natural Gas Hydrate

The formation–decomposition chemical equation of natural gas hydrate is:

(3)$C{H}_4\cdot N_h{H}_{2}O\leftrightarrow C{H}_4+N_h{H}_{2}O$

where $N_h$ is the hydration coefficient.

The decomposition can be described by the Kim-Bishoni decomposition kinetic equation, and the decomposition rate $m_h$ is [29,30]:

(4)$m_h=-k_d{M}_h{A}_s(P_e-P_g)$

where $A_s$ is the specific surface area of the hydrate reaction, 1/m; $M_h$ is the molar mass of hydrate, kg/mol; $P_g$ is the gas dynamic pressure, MPa.

$P_e$ is the phase equilibrium pressure, and its formula is [31]:

(5)$P_e=\exp (e_1-\frac{e_2}{T})$

where $T$ is the temperature of the porous medium containing gas hydrate; $e_1$, $e_2$ is the regression coefficient.

According to Equation (5) and parameter values in the literature [31], the phase equilibrium curve is shown in Figure 3.

$k_d$ is the kinetic reaction rate, $\mathrm{Pa}\cdot \mathrm{s}$, which is negatively correlated with temperature. The formula is [30]:

(6)$k_d=k_d^0\exp \left(-\frac{\Delta E_a}{RT}\right)$

where $k_d^0$ is the reaction kinetic constant; mol/(m²$\cdot$Pa$\cdot$s); $\Delta E_a$ is the reaction activation energy; $R$ is the gas constant.

The reaction-specific surface area is related to porosity, hydrate saturation and absolute permeability. The formula is [32]:

(7)$A_s=\varphi S_h\sqrt{\frac{{\varphi }^3{\left(1-S_h\right)}^3}{2k}}$

where $k$ is the absolute permeability. In this paper, a new permeability model is established by combining the Masuda permeability decline model with the cubic theorem in the form of the product [33]:

(8)$k=k_0{\left(1-S_h\right)}^N{\left(\frac{\varphi }{{\varphi }_0}\right)}^3$

where $k_0$ is the initial absolute permeability of hydrate-free porous media, m²; $N$ is the permeability decline coefficient; ${\varphi }_0$ is the porosity of hydrate-bearing sediments when hydrate saturation is 0.

The decomposition process of natural gas hydrate follows the law of conservation of mass, and its equation is:

(9)$\frac{\partial \left(\varphi {\rho }_h{S}_h\right)}{\partial t}={\rho }_h\left(\varphi \frac{\partial S_h}{\partial t}+S_h\frac{\partial \varphi }{\partial t}\right)=m_h$

where ${\rho }_h$ is the density of the natural gas hydrate, kg/m³.

The gas production rate $m_g$ and the water production rate $m_w$ obtained by the chemical reaction formula are:

(10)$m_g=-m_h\frac{M_g}{M_h}$

(11)$m_w=-N_h{m}_h\frac{M_w}{M_h}$

where $M_g$ and $M_w$ are the molar mass of methane and water, kg/mol, respectively.

2.4. Gas-Water Two-Phase Seepage Equation

The mass conservation equations of gas–water two-phase flow in a gas hydrate reservoir are [34]:

(12)$\frac{\partial \left(\varphi {\rho }_g{S}_g\right)}{\partial t}-\nabla \left(\frac{k_{rg}k{\rho }_g}{{\mu }_g}\nabla P_g\right)=m_g$

(13)$\frac{\partial \left(\varphi {\rho }_w{S}_w\right)}{\partial t}-\nabla \left(\frac{k_{rw}k{\rho }_w}{{\mu }_w}\nabla P_w\right)=m_w$

where $k_{rg}$ and $k_{rw}$ are the relative permeability of methane and water, respectively; ${\mu }_g$ and ${\mu }_w$ are the dynamic viscosity of methane and water, respectively,$\mathrm{Pa}\cdot \mathrm{s}$; ${\rho }_g$ and ${\rho }_w$ are the density of methane and water, respectively, kg/m³; $P_w$ is the pressure of reservoir water, MPa.

The methane density varies with temperature and is described by the ideal gas state equation:

(14)${\rho }_g=\frac{P_g{M}_g}{RT}$

The capillary force has a significant effect on gas-water two-phase flow. The magnitude is the difference between the wet phase and the non-wet phase, and the capillary force is related to the effective saturation of water. The calculation formula is [35]:

(15)$P_c=P_g-P_w$

(16)$P_c=p_e{\left(se\right)}^{-1/\lambda }$

where $p_e$ is the reference pressure of air intake value, MPa; $\lambda$ is the coefficient related to pore size distribution; $se$ is the effective saturation, defined as:

(17)$se=\frac{S_w-S_{wr}}{1-S_{wr}-S_{gr}}$

where $S_{wr}$ is the residual saturation of water; $S_{gr}$ is the residual saturation of methane.

The relative permeability is usually a function of the effective saturation of the wet phase. At present, the relative permeability model is based on the displacement experimental data of different fluids in different porous media [36,37,38,39,40,41]. The relative permeability model here uses the formula proposed by Leverett [40]:

(18)$k_{rg}={\left(1-se\right)}^2\left(1-s{e}^2\right)$

(19)$k_{rw}=\sqrt{se}{\left[1-{\left(1-s{e}^{1/m}\right)}^m\right]}^2$

where $m$ is the relative permeability coefficient.

2.5. Energy-Conservation Equation

The heat change during the decomposition of natural gas hydrate can be expressed as [42]:

(20)$C_T\frac{\mathrm{d}T}{\mathrm{d}t}+Q_T+Q_A=Q_h+q$

where $C_T$ is the equivalent specific heat capacity [43]; $C_T\frac{\mathrm{d}T}{\mathrm{d}t}$ is the heat change caused by temperature; $Q_T$ is the heat change caused by heat conduction [42]; $Q_A$ is the heat change caused by thermal convection [43]; $Q_h$ is the latent heat of the hydrate phase change [42]; $q$ is the heat absorbed by hydrate sediments from the external environment [27].

(21)$C_T=\varphi \left({\rho }_w{S}_w{C}_w+{\rho }_g{S}_g{C}_g+{\rho }_h{S}_h{C}_h\right)+\left(1-\varphi \right){\rho }_s{C}_s$

(22)$Q_T=-\nabla \left(\left(\varphi \left(S_w{\lambda }_w+S_g{\lambda }_g+S_h{\lambda }_h\right)+\left(1-\varphi \right){\lambda }_s\right)\cdot \nabla T\right)$

(23)${\mu }_{vg}=-\left(\frac{k_{rg}k}{{\mu }_g}\nabla P_g\right)$

(24)${\mu }_{vw}=-\left(\frac{k_{rw}k}{{\mu }_w}\left(\nabla P_g-\frac{\partial P_c}{\partial \left(se\right)}\cdot \nabla S_w\right)\right)$

(25)$Q_A=\nabla \left(\varphi S_g{C}_g{\rho }_g{\mu }_{vg}T+\varphi S_w{C}_w{\rho }_w{\mu }_{vw}T\right)$

(26)$Q_h=\frac{m_h}{M_h}\left(C_1+C_{2}T\right)$

(27)$q=h(T_c-T)$

where $C_w$, $C_g$, $C_h$ and $C_s$ are the specific heat capacity of the water, gas and hydrate phases and the sediment phase, respectively, J/(kg·K); ${\rho }_s$ is the sediment density, kg/m³; ${\lambda }_w$, ${\lambda }_g$, ${\lambda }_h$ and ${\lambda }_s$ are the thermal conductivity of water, gas and hydrate phase and sediment phase, W/(m·K); ${\mu }_{vg}$ and ${\mu }_{vw}$ are the Darcy velocities of methane and water, respectively, m/s; $C_1$, $C_2$ are regression coefficients; h is the boundary heat conduction coefficient, W/(m²·K); $T_c$ is the environment temperature, K.

2.6. Mechanical-Damage Equation

The influence of pore pressure and temperature change is considered in the equilibrium differential equation expressed in the form of displacement, and its tensor form is expressed as [10]:

(28)$\frac{G}{1-2\nu }\cdot u_{j,ji}+G\cdot u_{i,jj}+\alpha \cdot p_{,i}+K^{\prime }{\alpha }_T{T}_{,i}+F_i=0$

where $G$ is the shear modulus; $\nu$ is Poisson’s ratio; $F_i$ and $u_i$ are the components of body force and displacement in i direction, respectively; $\alpha$ is the Biot coefficient; ${\alpha }_T$ is the bulk thermal expansion coefficient of the whole medium, K⁻¹; $K^{\prime }$ is the drainage bulk modulus;

(29)$K^{\prime }=2G\left(1+\nu \right)/3\left(1-2\nu \right)$

$p$ is equivalent to the pore pressure [44]:

(30)$p=\frac{P_w{S}_w+P_g{S}_g}{S_w+S_g}=\frac{\sigma -{\sigma }^{\prime }}{\alpha }$

where $\sigma$ is the total stress of the element, MPa; ${\sigma }^{\prime }$ is the element effective stress, MPa.

The relationship between porosity and volume strain can be expressed as follows [45]:

(31)$\varphi =\frac{{\varphi }_0\left(1+{\epsilon }_0\right)+\alpha \left(\epsilon -{\epsilon }_0\right)}{1+\epsilon }$

where $\epsilon$ and ${\epsilon }_0$ are the volume strain at any time and the initial time, respectively.

Natural gas hydrate will inevitably lead to sediment damage in the process of decomposition. According to the theory of elastic damage, the elastic modulus of the element can be expressed as [10]:

(32)$E_h=(1-D)E_{t0}$

where $E_h$ is the elastic modulus of the hydrate at any time, GPa; $E_{t0}$ is the elastic modulus of the initial state (undecomposed) of the hydrate sediments, GPa; $D$ is the damage variable. It is assumed that the damage and its evolution are isotropic, so $E_h$, $E_{t0}$ and $D$ are scalars.

The relationship between elastic modulus and hydrate saturation can be obtained from the hydrate experimental data of Hyodo et al. [46,47], which can be expressed by exponential function:

(33)$E_h=E_0\cdot \exp \left(a_E\cdot S_h\right)$

where $E_0$ is the elastic modulus of the sediment when the hydrate saturation is 0, GPa; $a_E$ is the regression coefficient of the equation.

According to this evolution process, the elastic modulus of the initial state of hydrate sediments can be obtained:

(34)$E_{t0}=E_0\cdot \exp \left(a_E\cdot S_{h0}\right)$

where $S_{h0}$ is the initial saturation of the hydrate.

At this time, the elastic modulus of the hydrate at any time can be expressed as:

(35)$E_h=(1-D)E_0\cdot \exp \left(a_E\cdot S_{h0}\right)$

Combining Equations (33) and (35), the expression of the damage variable is obtained:

(36)$D=1-\frac{\exp \left(a_E\cdot S_h\right)}{\exp \left(a_E\cdot S_{h0}\right)}=1-\exp \left(a_E(S_h-S_{h0})\right)$

3. Numerical Simulation Description and Verification

In order to further verify the applicability of the THMC full coupling theoretical model of natural gas hydrate considering damage proposed in this paper, based on the depressurization decomposition test of Berea sandstone core carried out by Masuda [33,48,49], this paper verified the gas production and temperature changes of the THMC full coupling theoretical model of natural gas hydrate considering the damage and the THC full coupling theoretical model without considering the damage, and further compared and analyzed the evolution law of hydrate saturation, the pore pressure and the temperature of the THMC model and the THC model in the process of depressurization decomposition.

3.1. Simulation Conditions

The Berea core model is 300 mm long and 51 mm in diameter. The core is placed in a constant temperature water bath at a temperature of 2.3 °C. The test takes the left end as the outlet, the outlet end maintains a constant pressure of 2.84 MPa, the other boundaries are non-flow boundaries, and a constant confining pressure of 20 MPa is applied to the core. Three monitoring points, A, B and C, were arranged in the model, which were located at 35 mm, 150 mm and 225 mm from the depressurization end on the central axis of the model, respectively. The temperature of the monitoring point and the cumulative gas production at the outlet end were recorded. The model size, initial conditions and boundary conditions are shown in Figure 4. The initial condition parameters of the model are shown in Table 1, and the numerical simulation parameters are shown in Table 2.

The fully coupled model of natural gas hydrate depressurization mining is a non-homogeneous complex partial differential equation set with strong nonlinearity, and its analytical solution cannot be obtained by conventional methods. COMSOL Multiphysics provides a modeling environment based on partial differential equations, in which the numerical solution of partial differential equations can be calculated by the finite element method. The software provides a powerful interactive environment for modeling and solving complex scientific and engineering problems, and users can build models that can be viewed and used in a desktop environment through the model builder. Using this software, it is easy to correlate different physical fields, that is, to solve the scientific problem of multi-physics coupling control. Based on the scientific problems to be solved, after clarifying the basic physical process and constructing a mathematical model that can accurately describe the problem, the user can check whether there are modules related to such problems in the software. If there are, the user can specify the material properties and initial and boundary conditions according to their own needs and establish an applicable physical model. It is worth noting that the multi-module coupling function provided by the software itself enables users to set up and solve multi-physics problems more conveniently and clearly. If the physical field module of the software itself cannot meet the needs of users, users can use the mathematical module to input the mathematical model established by themselves, and then use the solution function of COMSOL to solve the partial differential equations. This paper uses COMSOL Multiphysics 5.6 released in 2020. The new version provides solvers with faster computing speeds and lower memory requirements for multi-core and cluster computing, more efficient CAD assembly processing functions, simulation App layout templates, and a series of image functions including tailoring planes, material rendering and partially transparent views. The THMC fully coupled model in this paper uses COMSOL to solve the coupled partial differential equation, in which the chemical decomposition field adopts the regional ordinary differential interface, the gas seepage field adopts the physical field from Darcy’s law, the water seepage field and the temperature field adopt the general form of ordinary differential interface and the mechanics and damage field adopt the solid mechanics’ physical field. In order to simplify the calculation, a two-dimensional axisymmetric model is used for simulation, as shown in Figure 5.

3.2. Model Verification

Figure 6 showed the cumulative gas production of the THC and THMC models in 0–300 min after depressurization decomposition compared with Masuda’s results. At the initial stage of the hydrate decomposition, the free gas in the model was rapidly discharged from the depressurization end, and the pore pressure was reduced, which destroyed the phase equilibrium, promoted the rapid decomposition of hydrate at the depressurization end and raised the gas production rate. After a period of rapid decomposition, the hydrate near the depressurization end was decomposed, and the pressure at the depressurization end was transmitted to the inside of the model, and the gas continued to be discharged, but the gas production rate gradually slowed down and entered the stable gas production stage. The cumulative gas production of the THMC model was small in the process of gas production because the decomposition of hydrate changed the microstructure of porous media. The hydrate bearing effective stress in the THMC model was decreasing, the mechanical deformation was further generated, the permeability was decreasing and the gas production rate was lower than that of the THC model. When the gas production was completed, the cumulative gas production of the two models was almost consistent.

The cumulative gas production of 0~300 min in Masuda’s experiment was 9091.1 cm³. The cumulative gas production of THC and THMC models was 9260.9 cm³ and 9181.7 cm³, respectively, which were 169.8 cm³ and 90.6 cm³ different from the experimental results. The cumulative gas production of the two numerical models was consistent with the general trend of the cumulative gas production of the experiment, and the THMC model was in good agreement.

Figure 7 showed the comparison of the temperature change curves of the THC model and THMC model with Masuda’s test at monitoring point A. Hydrate decomposition is an endothermic process. The heat sources included the sandstone core itself and the external constant temperature water bath. At the beginning of the experiment, the gas pressure reduction rate and the hydrate decomposition rate were faster. The heat taken away by the fluid flow and the heat absorbed by the decomposition were greater than the heat transferred from the external constant temperature water bath to the core, so the temperature gradually decreased. Measuring point A is close to the depressurization end, and the temperature begins to decrease first. In the experimental results, the temperature of measuring point A decreased to the lowest value of 273.78 K at 50.29 min after depressurization, which was 1.67 K lower than the initial temperature. The temperature of measuring point A in the THC model and THMC model decreased to the lowest value of 273.86 K and 273.89 K, respectively, at 16.68 min and 61.68 min after depressurization, which was 1.59 K and 1.56 K lower than the initial temperature, respectively. As the hydrate further decomposed, the decomposition rate gradually decreased and the gas production slowed down. At this time, the heat taken away by the fluid flow and the heat absorbed by the decomposition began to be less than the heat transferred from the external constant temperature water bath to the core, so the temperature gradually rose with time until it reached the external water bath temperature. In the phase of temperature reduction, the decomposition of the THMC model was relatively slow, and there was a lag phenomenon. In the phase of temperature rise, the heat transfer was faster, the time to reach the water bath temperature was shorter and the temperature rise trend was well-fitted with the test.

3.3. Result Analysis

In this section, the evolution of relevant parameters in the decomposition process of the THC and THMC models was compared and analyzed, and the influence of the damage effect was compared to further explain the importance of the damage effect in hydrate decomposition.

The pore pressure variation of the two groups of models is shown in Figure 8. As the hydrate at the depressurization end decomposed, methane gas and water could flow directly from the port in the initial stage, and the pore pressure of the model decreased rapidly. When further decomposed, the pore pressure gradually transferred to the inside of the model, the gas and water discharge rates decreased, the pore pressure decreased slowly and the gas production pressure was reached at 10,000 s.

Comparing the results of the two groups of models, the pore pressure of the THMC model decreased slowly, and the time to reach the gas production pressure was longer, which was significantly lagging behind the THC model. This was because the porosity of the THC model remained constant, and the pore pressure was transmitted to the model inside at a faster rate. The pore pressure of the THMC model decreased slowly, resulting in a slow increase in the effective stress. The model produced continuous compression deformation, and the porosity and permeability decreased slowly, so the final gas production rate would be slow.

The variation in hydrate saturation after solving the two models is shown in Figure 9, and the left end is the depressurization gas production end. After the start of depressurization, the hydrate state in the model reservoir would be below the phase equilibrium curve. The hydrate would gradually decompose from the depressurization end, producing methane gas and water. Driven by the pressure difference, methane and water gradually flowed to the depressurization end and promoted the further decomposition of the remaining hydrate in the reservoir. The hydrate had a faster decomposition rate in the initial period. When the two groups of models were decomposed to 10,000 s and 5000 s, respectively, the hydrate saturation at the boundary was less than 0.1. With the continuous decomposition process, the boundary was the first to complete the decomposition, and the decomposition front was constantly advancing towards the reservoir center. Finally, at 18,000 s, the hydrate was completely decomposed.

Comparing the two groups of models, the initial decomposition rate of the THC model was faster. As the decomposition process progressed, the THC model took the lead in basically completing the hydrate decomposition at 15,000 s. This was because the model ignored the effect of sediment mechanical behavior on pore volume, and the porosity remained constant throughout the decomposition process. Due to the rapid decomposition, the decomposition of the hydrate in the THMC model produced decomposition damage. The decomposition damage led to an increase in permeability and effective stress. The effective stress would compress the sediment and cause mechanical damage. The mechanical damage of the sediment led to a decrease in porosity and permeability. Therefore, the decomposition rate was slow and the time of complete decomposition of the hydrate was prolonged.

The temperature variation of the two groups of models is shown in Figure 10. Hydrate decomposition would inevitably lead to the redistribution of the temperature field. Since the simulated environment was heated by a water bath at a constant temperature, the boundary remained constant temperature. At the initial stage, the hydrate decomposition absorbed heat, resulting in the temperature of the depressurization end decreasing first. As the decomposition progressed, the temperature near the boundary where the first decomposition was completed gradually increased due to the influence of water bath heating. The decomposition of the undecomposed hydrate would still cause the temperature to decrease. Finally, the model returned to the temperature of the water bath heating at 20,000 s.

Compared with the two groups of models, the temperature variation of THMC in the initial temperature reduction stage was slow, and there was a lag characteristic. The temperature variation trend in the temperature-rising stage was in good agreement with the hydrate decomposition trend. The minimum temperature of point A of the THMC model was 273.89 K, which appeared at 61.68 min and gradually increased after 62 min. The minimum temperature of point A of the THC model was 273.86 K, which appeared at about 16.68 min and then gradually increased. Both models reached the water bath temperature at about 20,000 s.

The Von Mises stress distribution of the THMC model is shown in Figure 11. At the initial stage of depressurization, the overall stress of the model was about 18 MPa. When it was decomposed to 1000 s, the stress of the decomposed part of the depressurization end was significantly reduced, about 16 MPa. The stress concentration occurred in the part that would be decomposed, the stress was about 24 MPa and the stress of the undecomposed part was maintained at 19 MPa. With the decomposition, the stress concentration area gradually expanded along the central axis to the other end. When the decomposition was 8000 s, the high-stress area was merged, the maximum stress was about 23 MPa and then the stress began to decrease gradually. At 18,000 s, the overall stress of the model was maintained at about 20 MPa and remained stable.

The variation of the damage variable of the THMC model is shown in Figure 12. In the initial stage, the damage caused by hydrate decomposition was very small, and the maximum value of the damage variable in the model was still less than 0.1 at 1000 s. As the decomposition continued, obvious damage gradually occurred in the part of the first decomposition at the depressurization end of the model. At 8000 s, the maximum damage variable had reached 0.5 at both ends and the wall of the model, and the damage continued to develop along the central axis to the other end. When the hydrate was decomposed to 10,000 s, the damage variable at the depressurization end reached a maximum value of 0.53. When the hydrate was decomposed to 18,000 s, the overall damage variable of the reservoir reached the maximum value.

4. Influencing Factors and Sensitivity Analysis of Natural Gas hydrate Gas Production under Damage Effect

In the previous section, the reliability of the THMC model under the damage effect was verified, and the difference between the THMC and THC models was compared. In this section, based on the Masuda experimental model, point A was selected as the main monitoring point to discuss the specific influence of different initial hydrate saturation and different depressurization amplitude on the gas production rate, hydrate saturation, permeability and damage variables of natural gas hydrate under the damage effect.

In the simulation of this section, the initial hydrate saturation of the basic model was set to 0.5, the initial absolute permeability was set to 98 mD, and the initial gas pressure was set to 4 MPa. In order to reduce the compression deformation and ensure safe production, the confining pressure was set to 4 MPa, the depressurization range was set to 1 MPa, and the boundary of the depressurization end was set to 60 s to slowly reduce to production pressure, and then the pressure remained constant.

4.1. Effect Analysis of Different Initial Hydrate Saturation on Gas Production

The initial hydrate saturation characterizes the volume ratio of hydrate particles in the pores of hydrate sediments in the model. The larger the initial hydrate saturation is, the lower the effective porosity is, and the stronger the bearing capacity of the model is. In order to study the effect of different initial hydrate saturations on model gas production, three kinds of initial hydrate saturation of 0.3, 0.5 and 0.7 were selected for simulation analysis.

The decomposition of hydrate and the compression of the stress field led to the decrease in the volume strain of the model, and finally reached stability, as shown in Figure 13. When the initial hydrate saturation was larger, the absolute value of the initial volume strain of the model was smaller, because more hydrate particles were cemented in the pores of the sediment and bore the stress together with the sediment, resulting in smaller compression deformation. When the hydrate was completely decomposed, under the action of stress, different initial hydrate saturations produced the same size of compression deformation, but the smaller the initial hydrate saturation, the shorter the time to achieve stable deformation.

After the hydrate began to decompose, the permeability gradually increased from 0, and finally reached the initial absolute permeability of the model. The greater the initial hydrate saturation, the smaller the permeability of the model, because hydrate particles occupied the pore volume in porous media, reducing the seepage velocity; as the hydrate decomposed, the permeability increased continuously. When the hydrate decomposed completely, the permeability reached the initial absolute permeability. Therefore, under the larger initial hydrate saturation, the seepage was slow, and the further decomposition of hydrate was inhibited, the permeability increased slowly and the time to reach the initial absolute permeability was longer.

The variation of damage variable not only reflects the decrease in the elastic modulus of the sediment, but also characterizes the variation of the hydrate saturation. There was symmetry between the variation curve of the hydrate saturation and the variation curve of the damage variable. When the hydrate saturation reached 0, the damage variable reached the maximum and the hydrate was completely decomposed, as shown in Figure 14. The larger the initial hydrate saturation was, the longer the complete decomposition time was. This was because the larger initial hydrate saturation made the effective porosity decrease, and the gas and water produced by the decomposition were not easy to flow out, which inhibited the decomposition. The damage variable was related to the hydrate saturation. The greater the initial hydrate saturation was, the greater the peak value of the damage variable was, but the damage was slower in the initial gas production stage, which again showed that the larger initial hydrate saturation inhibited the initial decomposition and led to the slow development of the damage process. When the hydrate continued to decompose, the damage rate gradually accelerated. The greater the initial hydrate saturation, the faster the damage rate and the faster it finally achieved complete damage. When the initial hydrate saturation was 0.3, 0.5 and 0.7, the peak values of the damage variable under complete decomposition were 0.400, 0.573 and 0.696, respectively.

Since the depressurization end was slowly reduced to the production pressure within 1 min, the gas production rate change was divided into two stages. First, it increased rapidly from 0 to the peak of the gas production rate, and then decreased to a low gas-production rate until the gas production was completed, as shown in Figure 15a. The greater the initial hydrate saturation was, the smaller the peak gas production rate was, the longer the duration of the first stage was and the greater the gas production rate of the second stage was. When the initial hydrate saturation was 0.3, 0.5 and 0.7, the peak gas production rates were 6.72 cm³/s, 4.15 cm³/s and 1.41 cm³/s, respectively. Figure 15b showed the variation of cumulative gas production under different initial hydrate saturations. The greater the initial hydrate saturation was, the greater the cumulative gas production was. This was because the larger the initial hydrate saturation was, the larger the pore volume occupied by hydrate particles was, and the more methane gas was produced by decomposition. When the initial hydrate saturation was 0.3, 0.5 and 0.7, the cumulative gas production was 6452.46 cm³, 10,394.41 cm³ and 13,118.60 cm³, respectively. Consistent with the evolution of permeability, the volume strain and the damage variable, the larger the initial hydrate saturation, the quicker the hydrate decomposition was initially inhibited. Therefore, the cumulative gas production was larger when the hydrate saturation was 0.3 at the initial stage of gas production. With further decomposition, the porosity and permeability gradually increased, the decomposition rate accelerated and the damage continued to occur. Therefore, the gas production rate and gas production at a hydrate saturation of 0.7 gradually exceeded the other initial hydrate saturation conditions.

4.2. Analysis of the Impact of Different Depressurization Amplitude

Gas hydrate particles remain stable in sediment pores. When the pore pressure at the reservoir boundary is changed to provide a depressurization lower than the initial pore pressure, the phase equilibrium state of the hydrate is broken, and the hydrate particles decompose into methane gas and water and flow to the low-pressure end, promoting the continued decomposition of the hydrate. In order to study the influence of different depressurization amplitudes on model gas production, three depressurization amplitudes of 1 MPa, 2 MPa and 3 MPa were selected for simulation analysis.

The decomposition of the hydrate and the compression of the stress field led to the decrease in the volume strain of the model, and finally reached stability, as shown in Figure 16. The greater the depressurization, the greater the absolute value of the final volume strain produced by the model, because the hydrate particles would eventually completely decompose and discharge methane gas and water until the pore pressure at the depressurization end was reached. Under the same confining pressure, the greater the depressurization amplitude, the smaller the pore pressure, the greater the effective stress, the greater the volume strain, and the greater the depressurization amplitude, the shorter the time to reach a stable state.

After the hydrate began to decompose, the permeability gradually increased from 0 and finally reached the initial absolute permeability of the model. The greater the depressurization amplitude was, the greater the permeability of the model was. Because the larger depressurization amplitude moved the hydrate particles further away from the phase equilibrium state and was unable to exist in a solid form, a larger decomposition driving force was generated and continuously transmitted to the inside of the model, which promoted the rapid decomposition of hydrate particles in the model and increased the permeability. When the depressurization amplitude was 1 MPa, the difference with the initial pore pressure was small, and the hydrate solid in the model decomposed slowly. When the depressurization amplitude was 2 MPa and 3 MPa, the hydrate decomposed faster, the permeability increased faster and the time to reach the initial absolute permeability was shorter.

The hydrate saturation variation curve and the damage variable variation curve under different depressurization amplitudes were also symmetrical, as shown in Figure 17. The greater the depressurization was, the shorter the time for the complete decomposition of the hydrate. This was because the larger depressurization produced a larger driving force for decomposition, and the gas and water produced by decomposition flowed out rapidly, which continuously promoted decomposition. The depressurization amplitude had no effect on the peak value of the damage variable, but the larger depressurization amplitude would promote the rapid decomposition of hydrate particles, make the hydrate saturation decrease rapidly, accelerate the damage development process and, finally, achieve complete damage. At this time, the peak value of the damage variable was 0.573.

The gas production rate first increased rapidly from 0 to the peak and then decreased to a low gas production rate until the gas production was completed, as shown in Figure 18a. The greater the depressurization, the greater the peak gas production rate, the shorter the duration of the first stage and the smaller the gas production rate of the second stage. This was because the greater the magnitude of depressurization, the greater the driving force for hydrate decomposition, the faster the decomposition rate and the shorter the time to complete the decomposition, and, therefore, the second stage of continuous gas production time was shorter. The peak values of the gas production rate were 4.15 cm³/s, 7.04 cm³/s and 10.55 cm³/s, respectively, when the depressurization amplitude was 1 MPa, 2 MPa and 3 MPa.

The variation of cumulative gas production under different depressurization amplitudes is shown in Figure 18b. The larger the depressurization amplitude was, the larger the cumulative gas production was. This was because the larger the depressurization amplitude was, the smaller the final pore pressure was, and the larger the volume of decomposed gas and water was. The cumulative gas production at 1 MPa, 2 MPa and 3 MPa was 10,394.41 cm³, 10,529.56 cm³ and 11,150.80 cm³, respectively. Therefore, the larger depressurization amplitude would promote the decomposition of hydrates and contribute to the production of methane gas, but it would also cause large compression deformation. If the decomposition rate could not be controlled, the stability and safety of the model would be affected.

5. Conclusions

In this paper, a fully coupled depressurization production model of THC and THMC for natural gas hydrate was established. Masuda’s experiment was used to verify the reliability of the THMC model, and a sensitivity analysis of relevant parameters for the depressurization production process was conducted with the THMC model. Based on these results, the following conclusions were drawn:

• (1). A multi-physical field theoretical model of natural gas hydrate coupling thermal–hydraulic–mechanical–chemical was established and verified by Masuda’s experimental data. The results showed that the THMC model has better applicability than the THC model;

• (2). In the early stage of gas production, the larger initial hydrate saturation led to smaller initial volume strain and initial permeability, which inhibited the hydrate’s decomposition behavior, and the damage rate was slower. As the hydrate continued to decompose, the greater the damage rate under the larger initial hydrate saturation, the larger the peak value of the damage variable, and the larger the final cumulative gas production;

• (3). The larger depressurization amplitude provided a greater driving force for hydrate decomposition, promoted the rapid decomposition of the hydrate, improved the permeability in the decomposition process, accelerated the development of damage and increased the cumulative gas production. However, the larger pressure amplitude would eventually produce a larger volume strain, which would create stability problems at a faster gas production rate.

Author Contributions

Conceptualization, formal analysis, funding acquisition, writing, and editing, Y.Q., methodology, formal analysis, data curation, writing and editing, X.W., data curation, review, and editing, Z.W., methodology, validation, W.L. and T.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data associated with this research are available and can be obtained by contacting the corresponding author upon reasonable request.

Conflicts of Interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Nomenclature

$S_i$	saturation of phase i
$V_i$	volume of phase i (m3)
$V_{\varphi }$	volume of pores (m3)
$\varphi$	porosity
$V$	volume of sediment skeleton (m3)
$N_h$	hydration coefficient
$m_i$	decomposition or generation rate of phase i (kg·m−3·s−1)
$k_d$	kinetic reaction rate (Pa·s)
$M_i$	molar mass of phase i (kg·mol−1)
$A_s$	reaction specific surface area (m−1)
$P_e$	phase equilibrium pressure (MPa)
$P_g$	gas pore pressure (MPa)
$e_1$	phase equilibrium coefficient
$e_2$	phase equilibrium coefficient
$T$	temperature of porous media (K)
$k_d^0$	reaction kinetic constant (mol·m−2·Pa−1·s−1)
$\Delta E_a$	reaction activation energy (J·mol−1)
$R$	proctor gas constant (J·mol−1·K−1)
$k$	permeability (m−2)
$k_0$	initial absolute permeability (m−2)
$N$	permeability reduction coefficient
${\varphi }_0$	initial porosity
${\rho }_i$	density of phase i (kg·m−3)
$k_{rg}$	relative permeability of gas
${\mu }_g$	dynamic viscosity of gas (Pa·s)
$k_{rw}$	relative permeability of water
${\mu }_w$	dynamic viscosity of water (Pa·s)
$P_w$	pore water pressure (MPa)
$P_c$	capillary force (MPa)
$p_e$	atmos (MPa)
$se$	effective saturation
$\lambda$	aperture distribution coefficient
$S_{wr}$	residual saturation of water
$S_{gr}$	residual saturation of gas
$m$	relative permeability coefficient
$C_T$	equivalent specific heat capacity
$Q_T$	heat change caused by heat conduction
$Q_A$	heat change caused by thermal convection
$Q_h$	latent heat of hydrate phase change
$q$	heat absorbed by hydrate sediments from the external environment
$C_i$	specific heat capacity of phase i (J·kg−1·K−1)
${\lambda }_i$	thermal conductivity of phase i (W·m−1·K−1)
${\mu }_{vg}$	darcy velocity of gas (m·s−1)
${\mu }_{vw}$	darcy velocity of water (m·s−1)
$C_1$	regression coefficients
$C_2$	regression coefficients
$h$	boundary heat conduction coefficient
$T_c$	environment temperature (K)
$G$	shear modulus (GPa)
$\nu$	Poisson’s ratio
$u$	displacement (m)
$\alpha$	Biot coefficient
$p$	equivalent pore pressure (MPa)
$K^{\prime }$	drainage bulk modulus (GPa)
${\alpha }_T$	bulk thermal expansion coefficient (K−1)
$F$	body force (N)
$\sigma$	total stress (MPa)
${\sigma }^{\prime }$	effective stress (MPa)
${\epsilon }_0$	initial volume strain
$\epsilon$	volume strain
$E_h$	elastic modulus of hydrate (GPa)
$D$	damage variable
$E_{t0}$	initial elastic modulus of hydrate (GPa)
$E_0$	elastic modulus of sediment (GPa)
$a_E$	damage reduction factor
$S_{h0}$	initial saturation of hydrate
THMC	thermal-hydraulic-mechanical-chemical
THC	thermal-hydraulic-chemical
g	gas phase
w	water phase
h	solid hydrate phase
s	skeleton sand phase

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Research flow chart.

Figure 2. Coupling analysis of THMC fully coupled model.

Figure 3. Phase equilibrium curve of natural gas hydrate.

Figure 4. Initial and boundary conditions of the model.

Figure 5. Model grid and three-dimensional schematic diagram.

Figure 6. Comparison of cumulative gas production results [29].

Figure 7. Comparison of monitoring point temperature results [29].

Figure 8. Comparison of pore pressure evolution between THC and THMC models. (a) THC model; (b) THMC model.

Figure 9. Comparison of hydrate saturation evolution between THC and THMC models. (a) THC model; (b) THMC model.

Figure 10. Comparison of temperature evolution between THC and THMC models. (a) THC model; (b) THMC model.

Figure 11. Variation of Von Mises stress.

Figure 12. Variation of damage variable.

Figure 13. Variation of volume strain and permeability under different initial hydrate saturation.

Figure 14. Variation of hydrate saturation and damage variable under different initial hydrate saturation.

Figure 15. Gas production law under different initial hydrate saturation. (a) Variation of gas production rate; (b) variation of cumulative gas production.

Figure 16. Variation of volume strain and permeability under different depressurization amplitude.

Figure 17. Variation of hydrate saturation and damage variable under different depressurization amplitude.

Figure 18. Gas production law under different depressurization amplitude. (a) Variation of gas production rate; (b) variation of cumulative gas production.

Table 1

Model initial condition parameters [33].

Definition	Variable	Values and Units
Initial temperature of core	T 0	275.45 K
Water bath temperature of core	Tc	275.45 K
Initial porosity of core	ϕ 0	0.182
Initial gas pressure of core	Pg 0	3.75 MPa
Initial hydrate saturation	Sh 0	0.443
Initial water saturation	Sw 0	0.206
Initial gas saturation	Sg 0	0.351
Initial absolute permeability of core	k 0	98 mD

Table 2

Numerical simulation parameters [27].

Definition	Variable	Values and Units	Definition	Variable	Values and Units
Density of hydrate	ρh	917 kg/m3	Reference pressure value of intake	pe	1 × 105 Pa
Density of sediment	ρs	2150 kg/m3	Hydration coefficient	Nh	6
Gas density under standard conditions	ρg 0	0.652 kg/m3		ΔEa/R	9752.73 K
Density of water	ρw	1000 kg/m3	Aperture distribution coefficient	$\lambda$	2
Molar mass of hydrate	Mh	0.124 kg/mol	Relative permeability coefficient	m	0.6
Molar mass of gas	Mg	0.016 kg/mol	Residual saturation of water	Swr	0.1
Molar mass of water	Mw	0.018 kg/mol	Residual saturation of gas	Sgr	0.05
Specific heat capacity of hydrate	Ch	2220 J/(kg·K)	Latent heat coefficient of phase change	C 1	56,599
Specific heat capacity of sediment	Cs	750 J/(kg·K)	Latent heat coefficient of phase change	C 2	16.744
Specific heat capacity of gas	Cg	2180 J/(kg·K)	Heat conduction coefficient of boundary	h	25 W/(m2·K)
Specific heat capacity of water	Cw	4200 J/(kg·K)	Dynamic viscosity of gas	ug	1.84 × 10−5 Pa·s
Thermal conductivity of hydrate [50]	$\lambda$ h	0.6 W/(m·K)	Dynamic viscosity of water	uw	1.01 × 10−3 Pa·s
Thermal conductivity of sediment [51]	$\lambda$ s	0.8 W/(m·K)	Reaction kinetic constant	$k_d^0$	3.6 × 104 mol/(m2·Pa·s)
Thermal conductivity of gas	$\lambda$ g	0.056 W/(m·K)	phase equilibrium coefficient	e 1	3.98
Thermal conductivity of water	$\lambda$ w	0.5 W/(m·K)	phase equilibrium coefficient	e 2	8.5338
Gas constant	R	8.314 J/(mol·K)	Damage reduction factor	aE	1.7
Elastic modulus of sediment	E 0	1 × 109 Pa	Biot’s coefficient	$\alpha$	0.4
Poisson’s ratio of sediment	$\mathsf{\nu }$	0.2	Thermal expansion coefficient	$\alpha$ T	8 × 10−6 K−1

References

1. Zhao, J.F.; Song, Y.C.; Lim, X.L.; Lam, W.H. Opportunities and challenges of gas hydrate policies with consideration of environmental impacts. Renew. Sustain. Energy Rev.; 2017; 70, pp. 875-885. [DOI: https://dx.doi.org/10.1016/j.rser.2016.11.269]

2. Yu, Y.S.; Zhang, X.W.; Liu, J.W.; Lee, Y.; Li, X.S. Natural gas hydrate resources and hydrate technologies: A review and analysis of the associated energy and global warming challenges. Energy Environ. Sci.; 2021; 14, pp. 5611-5668. [DOI: https://dx.doi.org/10.1039/D1EE02093E]

3. Farahani, M.V.; Hassanpouryouzband, A.; Yang, J.H.; Tohidi, B. Insights into the climate-driven evolution of gas hydrate-bearing permafrost sediments: Implications for prediction of environmental impacts and security of energy in cold regions. RSC Adv.; 2021; 11, pp. 14334-14346. [DOI: https://dx.doi.org/10.1039/D1RA01518D] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35423992]

4. Sloan, E.D. Fundamental principles and applications of natural gas hydrates. Nature; 2003; 426, pp. 353-363. [DOI: https://dx.doi.org/10.1038/nature02135]

5. Wu, P.; Li, Y.H.; Sun, X.; Liu, W.G.; Song, Y.C. Mechanical Characteristics of Hydrate-Bearing Sediment: A Review. Energy Fuels.; 2020; 35, pp. 1041-1057. [DOI: https://dx.doi.org/10.1021/acs.energyfuels.0c03995]

6. Jiang, Y.J.; Li, M.; Luan, H.J.; Shi, Y.C.; Zhang, S.H.; Yan, P.; Li, B.C. Discrete element simulation of the Macro-Meso mechanical behaviors of Gas-Hydrate-Bearing sediments under dynamic loading. J. Mar. Sci. Eng.; 2022; 10, 1042. [DOI: https://dx.doi.org/10.3390/jmse10081042]

7. Jiang, Y.J.; Ma, X.Z.; Luan, H.J.; Liang, W.; Yan, P.; Song, W.Q.; Shan, Q.L. Numerical simulation on the evolution of physical and mechanical characteristics of natural gas hydrate reservoir during depressurization production. J. Nat. Gas Sci. Eng.; 2022; 108, 104803. [DOI: https://dx.doi.org/10.1016/j.jngse.2022.104803]

8. Zhang, N.; Wang, H.N.; Jiang, M.J. A mesoelastic-plastic damage model for hydrate-bearing sediments with various hydrate-growth patterns. Ocean Eng.; 2022; 266, 112919. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.112919]

9. Wang, Z.F.; Qi, C.Z.; Yu, H.; Zhao, F.; Luo, Y. Elastoplastic constitutive model for gas hydrate bearing sediments considering hydrate occurrence habits. Geomech. Geophys. Geo.; 2022; 8, 190. [DOI: https://dx.doi.org/10.1007/s40948-022-00494-z]

10. Zhu, W.C.; Wei, C.H.; Tian, J.; Yang, T.H.; Tang, C.A. Coupled thermal-hydraulic-mechanical model during rock damage and its preliminary application. Rock Soil Mech.; 2009; 30, pp. 3851-3857. [DOI: https://dx.doi.org/10.16285/j.rsm.2009.12.017]

11. Yang, Q.J.; Zhao, C.F. An improved constitutive model coupling elastoplasticity and damage for methane hydrate-bearing sediments. Ocean Eng.; 2014; 32, pp. 57–62, 92. [DOI: https://dx.doi.org/10.16483/j.issn.1005-9865.2014.05.006]

12. Shen, J.; Chiu, A.C.F.; Ng, C.W.W. Modeling the mechanical behavior of methane hydrate-bearing sand using the equivalent granular void ratio. J. Mar. Sci. Eng.; 2022; 10, 1040. [DOI: https://dx.doi.org/10.3390/jmse10081040]

13. Li, Y.L.; Liu, C.L.; Liu, L.L. Damage statistic constitutive model of hydrate-bearing sediments and the determination method of parameters. Acta Pet. Sin.; 2016; 37, pp. 1273-1279. [DOI: https://dx.doi.org/10.7623/syxb201610007]

14. Yan, R.T.; Zhang, B.H.; Yang, D.H.; Li, Y.; Chen, X.X.; Wei, C.F. Damage constitutive model for hydrate-bearing sediment under different temperature and pore pressure conditions. Rock Soil Mech.; 2018; 39, pp. 4421-4428. [DOI: https://dx.doi.org/10.16285/j.rsm.2017.1839]

15. Yan, M.Q.; Yan, R.T.; Yu, H.H. Strain-softening characteristics of hydrate-bearing sediments and modified Duncan-Chang model. Adv. Mater. Sci. Eng.; 2022; 2021, 2809370. [DOI: https://dx.doi.org/10.1155/2021/2809370]

16. Zhu, X.H.; Sun, H.W.; Zhao, J.Z.; Zhang, L.H.; Zhang, S. Damage constitutive model of equivalent variable elastic modulus for gas hydrate sediments. Acta Pet. Sin.; 2019; 40, pp. 1085-1094. [DOI: https://dx.doi.org/10.7623/syxb201909006]

17. Liu, L.; Chen, Y.Q.; Zhou, B.; Xue, S.F. A modified damage constitutive model of gas hydrate-bearing sediments based on discrete element method. Mar. Georesour. Geotechnol.; 2022; pp. 1-11. [DOI: https://dx.doi.org/10.1080/1064119X.2022.2106913]

18. Makogon, Y.F. Hydrates of Hydrocarbons; Penn Well Books: Tulsa, OK, USA, 1997; Volume 610, pp. 50-62.

19. Verigin, N.N.; Khabibullin, I.L.; Khalikov, G.A. Linear problem of the dissociation of the hydrates of a gas in a porous medium. Fluid Dymamics; 1980; 15, pp. 144-147. [DOI: https://dx.doi.org/10.1007/BF01089829]

20. Youslf, M.H.; Abass, H.H.; Selim, M.S.; Sloan, E.D. Experimental and theoretical investigation of methane gas hydrate dissociation in porous media. SPE Reserv. Eng.; 1991; 6, pp. 69-76. [DOI: https://dx.doi.org/10.2118/18320-PA]

21. Liang, W.; Zhao, T.B.; Qiu, Y.; Wang, X.F. Fully coupled numerical model and its application in natural gas hydrate reservoir. Energy Fuels; 2021; 35, pp. 2048-2063. [DOI: https://dx.doi.org/10.1021/acs.energyfuels.0c03465]

22. Klar, A.; Soga, K.; Ng, M.Y.A. Coupled deformation-flow analysis for methane hydrate extraction. Geotechnique; 2010; 60, pp. 765-776. [DOI: https://dx.doi.org/10.1680/geot.9.P.079-3799]

23. Sun, X. Numerical Simulation of Geo-Mechanical Behavior for Hydrate Bearing Sediments Considering the Influence of Hydrate Dissociation. Ph.D. Thesis; Dalian University of Technology: Dalian, China, 2017.

24. Wang, Z.M. Multi-Physics Coupled Numerical Simulation of Natural Gas Hydrate Exploitation. Master’s Thesis; Dalian University of Technology: Dalian, China, 2021; [DOI: https://dx.doi.org/10.26991/d.cnki.gdllu.2021.000846]

25. Zhang, X.L.; Xia, F.; Du, X.L.; Xu, C.S. Study on multi-field coupling model considering damage of hydrate-bearing sediments. Rock Soil Mech.; 2019; 40, pp. 4229–4239, 4305. [DOI: https://dx.doi.org/10.16285/j.rsm.2018.1521]

26. Huang, L.H.; Xu, C.S.; Xu, J.L.; Zhang, X.L.; Xia, F. The depressurization of natural gas hydrate in the multi-physics coupling simulation based on a new developed constitutive model. J. Nat. Gas Sci. Eng.; 2021; 95, 103963. [DOI: https://dx.doi.org/10.1016/j.jngse.2021.103963]

27. Liang, W. Study on Decomposition Characteristics of Natural Gas Hydrate and Multiphysics Coupling. Master’s Thesis; Shandong University of Science and Technology: Qingdao, China, 2020; [DOI: https://dx.doi.org/10.27275/d.cnki.gsdku.2020.000999]

28. Chen, Y.M. Theory and Technology of Natural Gas Hydrate Development; 1st ed. China University of Petroleum Press: Dongying, China, 2011; pp. 36-42.

29. Kim, H.C.; Bishnoi, P.R.; Heidemann, R.A.; Rizvi, S.S.H. Kinetics of methane hydrate decomposition. Chem. Eng. Sci.; 1987; 42, pp. 1645-1653. [DOI: https://dx.doi.org/10.1016/0009-2509(87)80169-0]

30. Clarke, M.; Bishnoi, P.R. Determination of the activation energy and intrinsic rate constant of methane gas hydrate decomposition. Can. J. Chem. Eng.; 2001; 79, pp. 143-147. [DOI: https://dx.doi.org/10.1002/cjce.5450790122]

31. Kamath, V.A. Study of Heat Transfer Characteristics during Dissociation of Gas Hydrates in Porous Media; 1st ed. University of Pittsburgh: Pittsburgh, PA, USA, 1984; pp. 1-24.

32. Sun, X.F.; Mohanty, K.K. Kinetic simulation of methane hydrate formation and dissociation in porous media. Chem. Eng. Sci.; 2006; 61, pp. 3476-3495. [DOI: https://dx.doi.org/10.1016/j.ces.2005.12.017]

33. Masuda, Y.; Fujinaga, Y.; Naganawa, S.; Fujita, K.; Hayashi, Y. Modeling and experiment studies on dissociation of methane gas hydrates in Berea sandstone cores. Proceedings of the Third International Gas Hydrate Conference; Salt Lake City, UT, USA, 18–22 July 1999; pp. 18-22.

34. Sun, X.; Nanchary, N.; Mohanty, K.K. 1-D Modeling of hydrate depressurization in porous media. Transp. Porous Media; 2005; 58, pp. 315-338. [DOI: https://dx.doi.org/10.1007/s11242-004-1410-x]

35. Ma, T.R.; Rutqvist, J.; Oldenburg, C.M.; Liu, W.Q.; Chen, J.G. Fully coupled two-phase flow and poromechanics modeling of coalbed methane recovery: Impact of geomechanics on production rate. J. Nat. Gas Sci. Eng.; 2017; 45, pp. 474-486. [DOI: https://dx.doi.org/10.1016/j.jngse.2017.05.024]

36. Chen, D.; Pan, Z.J.; Liu, J.S.; Connell, L.D. An improved relative permeability model for coal reservoirs. Int. J. Coal Geol.; 2013; 109, pp. 45-57. [DOI: https://dx.doi.org/10.1016/j.coal.2013.02.002]

37. Brooks, R.H.; Corey, A.T. Properties of porous media affecting fluid flow. J. Irrig. Drain. Div.; 1966; 92, pp. 61-88. [DOI: https://dx.doi.org/10.1061/JRCEA4.0000425]

38. Delshad, M.; Lenhard, R.J.; Oostrom, M.; Pope, G.A. A mixed-wet hysteretic relative permeability and capillary pressure model for reservoir simulations. SPE Reserv. Evaluation Eng.; 2003; 6, pp. 328-334. [DOI: https://dx.doi.org/10.2118/86916-PA]

39. Lomeland, F.; Ebeltoft, E. A new versatile capillary pressure correlation. Proceedings of the International Symposium of the Society of Core Analysts; Abu Dhabi, UAE, 9–12 October 2008; Volume 29.

40. Leverett, M.C. Capillary behavior in porous solids. Trans. AIME; 1941; 142, pp. 152-169. [DOI: https://dx.doi.org/10.2118/941152-G]

41. Genuchten, M.T.V. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.; 1980; 44, pp. 892-898. [DOI: https://dx.doi.org/10.2136/sssaj1980.03615995004400050002x]

42. Liang, H.F.; Song, Y.C.; Chen, Y.J. Numerical simulation for laboratory-scale methane hydrate dissociation by depressurization. Energy Convers. Manag.; 2010; 51, pp. 1883-1890. [DOI: https://dx.doi.org/10.1016/j.enconman.2010.02.018]

43. Liu, X.Q.; Sun, Y.; Guo, T.K.; Rabiei, M.; Qu, Z.Q.; Hou, J. Numerical simulations of hydraulic fracturing in methane hydrate reservoirs based on the coupled thermo-hydrologic-mechanical-damage (THMD) model. Energy; 2022; 238, 122054. [DOI: https://dx.doi.org/10.1016/j.energy.2021.122054]

44. Sun, X.; Luo, H.; Soga, K. A coupled thermal–hydraulic–mechanical–chemical (THMC) model for methane hydrate bearing sediments using COMSOL Multiphysics. J. Zhejiang Univ. A; 2018; 19, pp. 600-623. [DOI: https://dx.doi.org/10.1631/jzus.A1700464]

45. Liang, W.; Wang, J.G.; Li, P.B. Gas production analysis for hydrate sediment with compound morphology by a new dynamic permeability model. Appl. Energy; 2022; 322, 119434. [DOI: https://dx.doi.org/10.1016/j.apenergy.2022.119434]

46. Hyodo, M.; Wu, Y.; Nakashima, K.; Kajiyama, S.; Nakata, Y. Influence of fines content on the mechanical behavior of methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth; 2017; 122, pp. 7511-7524. [DOI: https://dx.doi.org/10.1002/2017JB014154]

47. Kajiyama, S.; Hyodo, M.; Nakata, Y.; Yoshimoto, N.; Wu, Y.; Kato, A. Shear behaviour of methane hydrate bearing sand with various particle characteristics and fines. Soils Found.; 2017; 57, pp. 176-193. [DOI: https://dx.doi.org/10.1016/j.sandf.2017.03.002]

48. Wan, Y.Z.; Wu, N.Y.; Hu, G.W.; Xin, X.; Jin, G.R.; Liu, C.L.; Chen, Q. Reservoir stability in the process of natural gas hydrate production by depressurization in the shenhu area of the south China sea. Nat. Gas Ind. B; 2018; 5, pp. 631-643. [DOI: https://dx.doi.org/10.1016/j.ngib.2018.11.012]

49. Yang, S.X.; Zhang, M.; Liang, J.Q.; Lu, J.G.; Zhang, Z.J.; Holland, M.; Schultheiss, P.; Fu, S.Y.; Sha, Z.B. Preliminary results of China’s third gas hydrate drilling expedition: A critical step from discovery to development in the South China Sea. Fire Ice; 2015; 15, pp. 1-5.

50. Waite, W.F.; Santamarina, J.C.; Cortes, D.D.; Dugan, B.; Espinoza, D.N.; Germaine, J.; Jang, J.; Jung, J.W.; Kneafsey, T.J.; Shin, H. et al. Physical properties of hydrate-bearing sediments. Rev. Geophys.; 2009; 47, RG4003. [DOI: https://dx.doi.org/10.1029/2008RG000279]

51. Sun, S.C.; Gu, L.L.; Yang, Z.D.; Lin, H.F.; Zhang, C.X. Gas hydrate dissociation by depressurization along with ice occurrence and sand migration. J. Nat. Gas Sci. Eng.; 2023; 109, 104853. [DOI: https://dx.doi.org/10.1016/j.jngse.2022.104853]