Study on the Evolution Law of Temperature, Pressure, and Productivity near the Well for Gas Hydrate Exploitation by Depressurization

Abstract

In this paper, a one-dimensional model of gas–water two-phase productivity for hydrate depressurization is established, which takes into account permeability variation and gas–water two-phase flow. By solving the coupled algebraic equations of dissociation front position, equilibrium temperature, and pressure in an iterative scheme, the movement law of the hydrate dissociation front and the evolution process of temperature and pressure near the well were obtained, and the effects of bottom hole pressure, reservoir temperature, and hydrate saturation on productivity were analyzed. The results show that the hydrate reservoir is divided into a decomposed zone and an undecomposed zone by the dissociation front, and the temperature and pressure gradients of the former are greater than those of the latter. Reducing bottom hole pressure, increasing reservoir temperature, and increasing hydrate saturation all lead to an increase in temperature and pressure gradient in the decomposed zone. Methane gas production is a sensitive function of bottom hole pressure, reservoir temperature, and hydrate saturation. The lower the bottom hole pressure, the higher the reservoir temperature, the lower the hydrate saturation (within a certain range), and the higher the gas production rate. The trend of the water production curve is the same as that of gas, but the value is 3–4 orders of magnitude smaller, which may be due to the large difference in the viscosity of gas and water, and the gas seepage speed is much larger than that of water.

1. Introduction

Natural gas hydrate (NGH) is a crystalline compound with a cage structure formed by hydrocarbon gas and water. NGH is a highly compressed natural gas resource, and 1 m³ of NGH can be decomposed to obtain 160–180 m³ (standard) of natural gas [1]. NGH exists widely in marine sediments and permafrost. The global methane carbon content in gas hydrates is estimated to be 10¹⁶ kg [2]. Therefore, NGH is recognized as a promising future energy source in the 21st century.

The principle of hydrate dissociation is based on changing the thermodynamic equilibrium of the three-phase system (water–hydrate–gas) [3,4]. In this study, we focus on the depressurization method for hydrate reservoirs. By reducing the bottom hole pressure, the equilibrium of the hydrate phase is broken, resulting in hydrate dissociation and the release of natural gas. Makogon (1997) used the classical Stefan problem to describe the process of hydrate dissociation and obtained a self-similar solution of the pressure distribution [5]. However, the model does not take into account the effect of water release. Based on the research of Makogon, Verigin et al. (1980) proposed a research model including the effect of water release, which considered the mass balance of gas and water at the dissociation front and assumed that the water produced by hydrate dissociation was stationary and did not affect the flow of gas [6]. Yousif and Sloan (1991) regarded the dissociation process in porous media as the Kim–Bishnoi dynamic isothermal process [7]. All the above models simplify the hydrate dissociation into an isothermal process, ignoring the thermal effect.

Holder and Angert (1982) evaluated the temperature distribution in the hydrate layer by using the conduction and heat transfer equation, considering the temperature change during the hydrate dissociation process [8]. Burshears et al. (1986) extended Holder and Angert’s model, taking into account the influence of water transport in the formation [9]. In their model, however, only heat conduction is included in the energy equation. Selim and Sloan (1989) considered convective heat transfer in a one-dimensional model and obtained analytical expressions of reservoir temperature and pressure distribution under the assumption that reservoir water does not flow and the well temperature remains unchanged [10]. Makogon (1997) considered the influence of heat conduction, heat convection, and throttling processes in revealing the interface heat transfer process of hydrate dissociation, and finally obtained the analytical expression of temperature and pressure in a one-dimensional model [5]. Tsypkin (2000) described the movement of gas and water in hydrate reservoirs, respectively, and obtained self-similar solutions of temperature and pressure distributions by using the micro-perturbation method [11]. Zhao and Shang (2010) further considered the influence of hydrate dissociation heat and added the energy conservation equation on the dissociation front, proving that reducing bottom hole pressure can significantly increase gas production [12]. Hao et al. (2020) established a semi-analytical productivity model of NGH by depressurization and analyzed the factors affecting productivity [13]. Based on the classical analytical model of hydrate dissociation proposed by Selim and Sloan, Li et al. (2023) further considered geomechanics and established a new analytical model [14]. With the gradual deepening of research, scholars have gradually adopted numerical models to characterize hydrate dissociation and production. Yu et al. (2023) analyzed the production and dissociation characteristics of NGH by numerical method and considered the influence of reservoir stress in the model [15]. Liu et al. (2023) further used a multi-scale numerical model to simulate the dissociation process of methane hydrate in porous media, including two-phase flow, heat and mass transfer, dissociation kinetics, and hydrate structure evolution [16].

Both field data and experimental studies show that hydrate dissociation causes changes in reservoir permeability, and gas–water two-phase seepage after decomposition is an important factor affecting productivity [17,18]. Based on this, a one-dimensional gas–water, two-phase productivity model for hydrate depressurization is proposed in this paper. Considering hydrate saturation variation and water phase flow, a set of approximate self-similar formats is used to obtain the variational solutions of pressure, temperature, and productivity in the reservoir. This method is solved in an iterative format. The distribution of pressure and temperature near the well, the location of the front, and the evolution of gas/water production over time were obtained.

This paper mainly consists of three parts. Among them, Section 2 is the model building; Section 3 is the model solving; Section 4.1 is the dynamic analysis of hydrate dissociation; and Section 4.2 is the analysis of influencing factors in hydrate productivity.

2. Model Building

2.1. Physical Model and Basic Assumptions

Assume that there is a large pressurized methane hydrate reservoir underground, the reservoir pressure is pe, the reservoir temperature is T_e, and there are solid hydrates and natural gas in the reservoir pores. The hydrate is in a stable state at the initial temperature and pressure. When the bottom hole pressure p_G drops to a certain value, that is, p_G < p_D < p_e, where p_D is the dissociation pressure of hydrate at the dissociation temperature T_D. At this point, the balance of the hydrate phase near the well is disrupted, and it begins to decompose into natural gas and water. Over time, the process of hydrate dissociation spreads outward. It is assumed that hydrate dissociation in porous media occurs in a narrow region, which can be regarded as a surface, the so-called dissociation front. This moving front divides the volume of the reservoir into two distinct zones. The area near the well is the hydrate decomposed zone, where hydrate, water, and gas coexist in three phases, while the area far away from the well is the undecomposed zone, where only solid hydrate and natural gas exist. Because of the pressure gradient, the gas and water move inside the well, while the dissociation front moves in the opposite direction. The model adopts three important assumptions: (1) The pressure and temperature at any point of the dissociation front are the equilibrium pressure p_D and the equilibrium temperature T_D; (2) The hydrate reservoir is assumed to be porous and contain natural gas. When the dissociation front moves towards the hydrate zone, heat must be supplied to the front edge due to the endothermic property of the hydrate dissociation process. In this case, heat conduction is ignored compared with thermal convection. (3) The decomposed zone contains three phases of hydrate, methane gas and water, in which only two phases of gas and water flow.

2.2. Basic Mathematical Model

In this section, the basic mathematical model adopted by Makogon (1997) [5] and Ji et al. (2001) [19] is followed. The difference is that gas–water two-phase seepage is considered in the decomposed zone. Assuming a hydrate reservoir, the well location is shown in Figure 1. Using a one-dimensional model, for the decomposed zone, the gas–water two-phase pressure equation can be seen according to the mass conservation equation as follows:

$${\displaystyle \frac{2{\varphi }_1{\mu }_g}{K_1{k}_{rg1}}}{\displaystyle \frac{\partial p_1}{\partial t}}={\displaystyle \frac{{\partial }^2{p}_1^2}{\partial x^2}}$$

(1)

$${\displaystyle \frac{2{\varphi }_1{\mu }_w}{K_1{k}_{rw1}}}{\displaystyle \frac{\partial p_1}{\partial t}}={\displaystyle \frac{{\partial }^2{p}_1^2}{\partial x^2}}$$

(2)

where

$${\varphi }_1=(1-S_H-S_w)\varphi$$

(3)

Only free methane gas flows in the undecomposed zone. Similarly, the pressure expression in the undecomposed region can be expressed as:

$${\displaystyle \frac{2{\varphi }_2{\mu }_g}{K_2{k}_{rg2}}}{\displaystyle \frac{\partial p_2}{\partial t}}={\displaystyle \frac{{\partial }^2{p}_2^2}{\partial x^2}}$$

(4)

where

$${\varphi }_2=(1-S_{w0}-S_{H0})\varphi$$

(5)

In the above equation, ϕ₁ and ϕ₂ are the porosities of the decomposed and undecomposed zones, respectively; μ_g and μ_w are the viscosities of gas and water, respectively; K₁ is the absolute permeability of the decomposed zone, and K₂ is the initial absolute permeability of the reservoir; k_rg1 and k_rg2 are gas relative permeability in the decomposed and undecomposed zones, respectively; k_rw1 is the relative permeability of water in the decomposed zone; k_rg0 is the initial gas relative permeability of the reservoir; p₁ and p₂ represent the pressure in the decomposed and undecomposed zones, respectively; S_H is the hydrate saturation in the decomposed zone; S_w is the water saturation of the decomposed zone; S_w0 is the original water saturation of the reservoir; and S_H0 is the initial hydrate saturation of the reservoir.

According to the hypothesis of the model, conduction heat transfer is much smaller than convective heat transfer in the case of natural gas in the hydrate reservoir. Therefore, the heat convection effect of fluid is considered in this paper, and the heat conduction effect is ignored. The temperature distribution equation of the hydrate reservoir is as follows:

$${\displaystyle \frac{\partial T_1}{\partial t}}={\displaystyle \frac{c_v{K}_1}{c_1{\mu }_1}}{\displaystyle \frac{\partial p_1}{\partial x}}({\displaystyle \frac{\partial T_1}{\partial x}}-\delta {\displaystyle \frac{\partial p_1}{\partial x}})+\eta {\displaystyle \frac{{\varphi }_1{c}_v}{c_1}}{\displaystyle \frac{\partial p_1}{\partial t}}$$

(6)

$${\displaystyle \frac{\partial T_2}{\partial t}}={\displaystyle \frac{c_v{K}_2}{c_2{\mu }_2}}{\displaystyle \frac{\partial p_2}{\partial x}}({\displaystyle \frac{\partial T_2}{\partial x}}-\delta {\displaystyle \frac{\partial p_2}{\partial x}})+\eta {\displaystyle \frac{{\varphi }_2{c}_v}{c_2}}{\displaystyle \frac{\partial p_2}{\partial t}}$$

(7)

where T₁ and T₂ are the temperatures of the decomposed and undecomposed zones, respectively; c_v is the volume heat capacity of gas; c₁ and c₂ are the comprehensive specific heat capacities of the decomposed and undecomposed zones, respectively; δ is the throttle coefficient of gas; η is the adiabatic coefficient of gas.

$$\begin{array}{l}{c}_1=\varphi (S_w{c}_w+S_g{c}_g+S_H{c}_H)+(1-\varphi )c_r\\ c_2=\varphi (S_{w0}{c}_w+S_{g0}{c}_g+S_{H0}{c}_H)+(1-\varphi )c_r\\ {\mu }_1=S_w{\mu }_w+(1-S_w){\mu }_g\\ {\mu }_2=S_{w0}{\mu }_w+(1-S_{w0}){\mu }_g\end{array}$$

(8)

where c_w, c_g, c_H, and c_r represent the specific heat capacities of water, gas, hydrate, and rocks, respectively; S_g is the gas saturation in the decomposed zone; S_g0 is the original gas saturation of the reservoir.

The auxiliary equation is:

$$S_g+S_w+S_H=1$$

(9)

The initial and boundary conditions are:

$$\begin{array}{l}{p}_1(0,t)=p_G\\ p_2(x,0)=p_2(\infty ,t)=p_e\end{array}$$

(10)

$$p_1(l(t),t)=p_2(l(t),t)=p_D(T_D)$$

(11)

$$T_2(x,0)=T_2(\infty ,t)=T_e$$

(12)

$$T_1(l(t),t)=T_2(l(t),t)=T_D$$

(13)

where p_G is the bottom hole pressure, p_e is the initial reservoir pressure, T_e is the initial reservoir temperature, p_D is the dissociation front pressure, and T_D is the dissociation front temperature. According to the hypothesis of the model, p_D and T_D are the equilibrium pressure and temperature. l(t) is the distance between the dissociation front and the bottom hole.

Constant bottom hole pressure and constant reservoir pressure are, respectively, used to approximate the pressure treatment, as follows:

$$\begin{array}{l}{\displaystyle \frac{\partial p_1^2}{\partial t}}\approx 2p_G{\displaystyle \frac{\partial p_1}{\partial t}}\\ {\displaystyle \frac{\partial p_2^2}{\partial t}}\approx 2p_e{\displaystyle \frac{\partial p_2}{\partial t}}\end{array}$$

(14)

Equations (1), (2), and (4) are linearized to:

$${\displaystyle \frac{\partial p_1^2}{\partial t}}={\chi }_{1i}{\displaystyle \frac{{\partial }^2{p}_1^2}{\partial x^2}}$$

(15)

$${\displaystyle \frac{\partial p_2^2}{\partial t}}={\chi }_2{\displaystyle \frac{{\partial }^2{p}_2^2}{\partial x^2}}$$

(16)

where the subscripts i = g, w, respectively, represent the gas phase and the water phase.

$$\begin{array}{l}{\chi }_{1i}={\displaystyle \frac{K_1{k}_{ri1}{p}_G}{\varphi (1-S_w-S_H){\mu }_i}}\\ {\chi }_2={\displaystyle \frac{K_2{k}_{rg2}{p}_e}{\varphi (1-S_{w0}-S_{H0}){\mu }_g}}\end{array}$$

(17)

From the boundary condition relation (10)–(13), self-similar solutions of Equations (15) and (16) can be obtained according to Makogon’s (1997) method [5]:

$$p_1^2=p_G^2-(p_G^2-p_D^2){\displaystyle \frac{erf{\lambda }_{1i}}{erf{\alpha }_{1i}}}$$

(18)

$$p_2^2=p_e^2-(p_e^2-p_D^2){\displaystyle \frac{erf{\lambda }_2}{erf{\alpha }_2}}$$

(19)

where

$$\begin{array}{l}{\lambda }_{1i}={\displaystyle \frac{x}{2\sqrt{{\chi }_{1i}t}}}\\ {\lambda }_2={\displaystyle \frac{x}{2\sqrt{{\chi }_{2}t}}}\end{array}$$

(20)

$$\begin{array}{l}{\alpha }_{1i}=\sqrt{{\displaystyle \frac{\gamma }{4{\chi }_{1i}}}}\\ {\alpha }_2=\sqrt{{\displaystyle \frac{\gamma }{4{\chi }_2}}}\end{array}$$

(21)

$$l(t)=\sqrt{\gamma t}$$

(22)

In the above equation, γ represents the constant of the dissociation front movement, which needs to be solved iteratively in this paper. The error function and complementary error function are defined as:

$$\begin{array}{l}erf(\xi )={\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle {\int }_0^{\xi }{e}^{-{\eta }^2}}d\eta \\ erfc(\xi )=1-erf(\xi )\end{array}$$

(23)

Similarly, self-similar solutions to Equations (6) and (7) can be obtained:

$$T_1=T_D+A_{1i}\delta \left[erf{\lambda }_{1i}-erf{\alpha }_{1i}+({\displaystyle \frac{\eta }{\delta }}{B}_1-1)({\Psi }_1({\lambda }_{1i})-{\Psi }_1({\alpha }_{1i}))\right]$$

(24)

$$T_2=T_e-A_2\delta \left[erfc{\lambda }_2+({\displaystyle \frac{\eta }{\delta }}{B}_2-1){\Psi }_2({\lambda }_2)\right]$$

(25)

where

$$\begin{array}{l}{\Psi }_1({\xi }_1)={\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle {\int }_0^{{\xi }_1}\frac{\eta e^{-{\eta }^2}}{\eta +C_{1i}{e}^{-{\eta }^2}}}d\eta \\ {\Psi }_2({\xi }_2)={\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle {\int }_{{\xi }_2}^{\infty }\frac{\eta e^{-{\eta }^2}}{\eta +C_2{e}^{-{\eta }^2}}}d\eta \end{array}$$

(26)

$$\begin{array}{l}{A}_{1i}={\displaystyle \frac{1}{2erf{\alpha }_{1i}}}{\displaystyle \frac{p_D^2-p_G^2}{p_G}}\\ A_2={\displaystyle \frac{1}{2erfc{\alpha }_2}}{\displaystyle \frac{p_e^2-p_D^2}{p_e}}\end{array}$$

(27)

$$\begin{array}{l}{B}_1={\displaystyle \frac{{\varphi }_1{c}_v}{c_1}}\\ B_2={\displaystyle \frac{{\varphi }_2{c}_v}{c_2}}\end{array}$$

(28)

$$\begin{array}{l}{C}_{1i}={\displaystyle \frac{p_D^2-p_G^2}{p_G}}{\displaystyle \frac{c_v}{c_1}}{\displaystyle \frac{1}{2\sqrt{\pi }erf{\alpha }_{1i}}}{\displaystyle \frac{K_1{k}_{ri1}}{{\mu }_1{\chi }_{1i}}}\\ C_2={\displaystyle \frac{p_e^2-p_D^2}{p_e}}{\displaystyle \frac{c_v}{c_2}}{\displaystyle \frac{1}{2\sqrt{\pi }erfc{\alpha }_2}}{\displaystyle \frac{K_2{k}_{rg2}}{{\mu }_2{\chi }_2}}\end{array}$$

(29)

2.3. Dissociation Front Motion

At the dissociation front, according to the phase equilibrium of natural gas hydrates, the relationship between T_D and p_D can be expressed as:

$${\log }_{10}{p}_D=a(T_D-T_0)+b{(T_D-T_0)}^2+c$$

(30)

where T₀ is the reference temperature, which is 273.15 K in this paper; a, b, and c are empirical constants related to hydrate composition. Makogon (1997) [5] used the least squares method to obtain three values of 0.0342/K, 0.0005/K², and 6.4804 based on the equilibrium temperature and pressure data of methane hydrate.

The mass balance equation of gas in the dissociation front is obtained by Verigin et al. (1980) [6]:

$${\rho }_1{v}_1-{\rho }_2{v}_2=-\left[S_H\epsilon {\rho }_H-(1-S_w){\rho }_1+(1-S_H){\rho }_2\right]\varphi {\displaystyle \frac{dl}{dt}}$$

(31)

In the above equation, ρ₁ and ρ₂ are methane gas densities in the decomposed and undecomposed zones, respectively; ρ_H is the hydrate density; ε is the mass fraction of gas in methane hydrate; and v₁ and v₂ are methane gas velocities in the decomposed and undecomposed zones, respectively.

At the dissociation front,

$${\rho }_1(l,t)={\rho }_2(l,t)={\rho }_0{\displaystyle \frac{p_D{T}_0}{z{p}_0{T}_D}}$$

(32)

z is the compression factor of methane gas; ρ₀ is the density of methane gas in the standard state (pressure p₀, temperature T₀).

Substituting Equation (32) into Equation (31) yields:

$$v_1(l,t)-v_2(l,t)=-\left[\epsilon S_H{\displaystyle \frac{{\rho }_H}{{\rho }_0}}{\displaystyle \frac{p_0}{p_D}}{\displaystyle \frac{T_D}{T_0}}z-(S_H-S_w)\right]\varphi {\displaystyle \frac{dl}{dt}}$$

(33)

The mass balance equation of water is:

$${\rho }_w{\varphi }_1{S}_w=(1-\epsilon ){\rho }_H{\varphi }_2{S}_H$$

(34)

where ρ_w is the density of water.

Substituting Equations (18) and (19) into Equation (33) yields an equation with the determinative constant γ:

$$K_1{k}_{ri1}{\displaystyle \frac{p_D^2-p_G^2}{\sqrt{\pi {\chi }_{1i}}}}{\displaystyle \frac{e^{-{\alpha }_{1i}^2}}{erf{\alpha }_{1i}}}-K_2{k}_{rg2}{\displaystyle \frac{p_e^2-p_D^2}{\sqrt{\pi {\chi }_2}}}{\displaystyle \frac{e^{-{\alpha }_2^2}}{erf{\alpha }_2}}=A\sqrt{\gamma }$$

(35)

$$A=\left[\epsilon S_H{\displaystyle \frac{{\rho }_H{p}_0{T}_D}{{\rho }_0{T}_0}}z-(S_H-S_w)p_D\right]\varphi {\mu }_1$$

(36)

In addition, at the dissociation front, λ₂ = α₂, Equation (25) becomes:

$$T_D=T_e-A_2\delta \left[erfc{\alpha }_2+({\displaystyle \frac{\eta }{\delta }}{B}_2-1){\Psi }_2({\alpha }_2)\right]$$

(37)

2.4. Permeability Model of Decomposition Zone

The hydrate dissociation zone is a three-phase flow zone of gas–water–hydrate. With the dissociation of hydrate, the absolute permeability of the decomposed zone increases gradually, and the relative permeability of gas and water changes constantly. Masuda et al. (1999) [20] gave the empirical equation for the change in absolute permeability with hydrate saturation:

$$K=K_0{(1-S_H)}^N$$

(38)

In the above equation, K₀ is the absolute reservoir permeability when hydrate saturation equals 0; N is the permeability decline index. The Masuda model has been widely cited and improved by a large number of scholars, but different hydrate reservoirs have different N values, which need to be obtained by fitting the seepage experiment data. Minagawa et al. (2005) used an exponential relationship to fit the permeability data under different hydrate saturations and obtained a decline index N ranging from 2.5 to 9.8 [21]. Liang et al. (2011) believe that the value of N is between 2 and 15 [22]. In this paper, the N value is 5.

The gas–water relative permeability is based on the algorithm proposed by Hong and Pooladi-Darvish [23,24]:

$$k_{rg}=k_{rg0}{\overline{S_g}}^{1/2}{(1-{\overline{S_{wH}}}^{1/m})}^{2m}$$

(39)

$$k_{rw}=k_{rw0}{\overline{S_w}}^{1/2}{\left[1-{(1-{\overline{S_w}}^{1/m})}^m\right]}^2$$

(40)

where

$$\begin{array}{l}\overline{S_g}={\displaystyle \frac{1-S_w-S_H-S_{gr}}{1-S_{wr}-S_{gr}}}\\ \overline{S_{wH}}={\displaystyle \frac{S_w+S_H-S_{wr}}{1-S_{wr}-S_{gr}}}\\ \overline{S_w}={\displaystyle \frac{S_w-S_{wr}}{1-S_{wr}-S_{gr}}}\end{array}$$

In the equation, m = 0.45, S_wr = 0.3, S_gr = 0.05, k_rw0 = 0.5, and k_rg0 = 1.0.

3. Model Solving

In the decomposed zone, given a hydrate saturation S_H, the corresponding gas saturation S_g and water saturation S_w can be obtained by combining Equations (9) and (34), and then the absolute permeability and relative permeability of gas and water under different hydrate saturations can be calculated. The nonlinear Equations (30), (35) and (37) are solved iteratively to obtain T_D, p_D, and γ values under different conditions. Methane gas production (per unit length) can be obtained by the following equation:

$$Q_g={\displaystyle \frac{K_1{k}_{rg1}}{{\mu }_1}}{\displaystyle \frac{\partial p_1(0,t)}{\partial x}}={\displaystyle \frac{K_1{k}_{rg1}}{{\mu }_1}}{\displaystyle \frac{p_D^2-p_G^2}{p_G}}{\displaystyle \frac{1}{erf{\alpha }_{1g}}}{\displaystyle \frac{1}{2\sqrt{\pi {\chi }_{1g}t}}}$$

(41)

The water production per unit length is:

$$Q_w={\displaystyle \frac{K_1{k}_{rw1}}{{\mu }_1}}{\displaystyle \frac{\partial p_1(0,t)}{\partial x}}={\displaystyle \frac{K_1{k}_{rw1}}{{\mu }_1}}{\displaystyle \frac{p_D^2-p_G^2}{p_G}}{\displaystyle \frac{1}{erf{\alpha }_{1w}}}{\displaystyle \frac{1}{2\sqrt{\pi {\chi }_{1w}t}}}$$

(42)

4. Result Analysis

A one-dimensional dissociation model is established, and the basic parameters of the reservoir are shown in

Table 1. Parameters of hydrate reservoir.

Reservoir Parameter	Value	Reservoir Parameter	Value
Initial reservoir pressure pe (MPa)	15	Density of hydrate ρH (kg·m−3)	900
Porosity ϕ	0.4	Methane gas density under standard conditions ρ0 (kg·m−3)	0.7
Gas compression factor z	0.8	Density of water ρw (kg·m−3)	1000
Volume heat capacity of gas cv (J·kg−1·K−1)	3000	Initial hydrate saturation SH0	0.6
Specific heat capacity of water cw (J·kg−1·K−1)	4200	Initial water saturation Sw0	0.2
Specific heat capacity of gas cg (J·kg−1·K−1)	3000	Initial gas saturation Sg0	0.2
Specific heat capacity of hydrate cH (J·kg−1·K−1)	2200	Mass fraction of gas in hydrate ε	0.129
Specific heat capacity of rock cr (J·kg−1·K−1)	1800	Viscosity of water μw (mPa·s)	1.2
		Viscosity of gas μg (mPa·s)	0.02
		Throttling coefficient of gas δ (K·Pa−1)	8 × 10−7
		Adiabatic coefficient of gas η (K·Pa−1)	3.2 × 10−7

. In this section, numerical results of the time evolution of hydrate reservoir pressure and temperature profiles under different conditions are given. In addition, changes in methane gas/water production over time and the location of the dissociation front were evaluated. The sensitivity of natural gas production to different reservoir parameters was discussed.

4.1. Dynamic Analysis of Hydrate Depressurization

Based on the given reservoir temperature T_e (15 °C, 12 °C, 10 °C), bottom hole pressure p_G (10 MPa, 8 MPa, 5 MPa), and hydrate saturation S_H (0.5, 0.3, 0.2), the Equations (30), (35) and (37) are solved iteratively. Hydrate dissociation temperature T_D, dissociation pressure p_D, and values representing the movement of the hydrate dissociation front were obtained (

Table 2. Values of dissociation temperature, dissociation pressure, parameter γ, and the distance l from the dissociation front to the wellbore for a given reservoir and well condition.

Te (°C)	pG (MPa)	SH (MPa)	TD (°C)	pD (MPa)	γ (m2/s)	l (m)
15	10	0.5	14.92	12.62	1.26 × 10−4	18.09
15	8	0.5	14.92	12.62	1.68 × 10−4	20.88
15	5	0.5	14.92	12.62	2.14 × 10−4	23.56
12	10	0.5	11.87	10.94	1.02 × 10−4	16.25
10	10	0.5	9.85	10.44	7.69 × 10−5	14.12
15	10	0.3	14.92	12.65	3.90 × 10−3	100.09
15	10	0.2	14.92	12.65	1.64 × 10−2	206.02

).

As can be seen from the table, when the bottom hole pressure gradually decreases from 10 MPa to 8 MPa and 5 MPa, the dissociation temperature and pressure remain unchanged, which are 14.92 °C and 12.62 MPa. Similarly, as the hydrate saturation decreases, the dissociation temperature and pressure remain almost constant. Reservoir pressure, however, has a great influence on the dissociation temperature and pressure. With the decrease in reservoir temperature, the dissociation temperature and pressure of hydrate decrease sharply. Temperature is particularly affected; for every 1 °C decrease in reservoir temperature, the dissociation temperature decreases by about 1 °C. The above conclusion, that hydrate dissociation temperature and pressure are sensitive functions of reservoir temperature, is consistent with the research of Ji et al. (2001) [19].

Table 2 also shows that under the premise of keeping other conditions unchanged, the lower the bottom hole pressure, the faster the hydrate dissociation front moves, and the farther the dissociation front spreads to the depth of the reservoir. Additionally, under the same conditions, the lower reservoir temperature provides less energy for hydrate dissociation, resulting in slower movement of the dissociation front. In addition, the decrease in hydrate saturation leads to a sharp increase in the movement velocity of the dissociation front. When the hydrate saturation gradually decreases from 0.5 to 0.3 and 0.2, the movement distance of the dissociation front rapidly increases from 14.12 m to 100.09 m and 206.02 m.

Figure 2 shows the variation curve of the dissociation front with time under different conditions. It can be seen that the front advances with time, and the advancing speed gradually slows down. It may be that in the early stage, hydrate dissociation occurs near the bottom of the well, the dissociation rate is fast, and the corresponding dissociation front moves faster. As the dissociation continues, the front moves deeper into the reservoir, the dissociation speed slows down, and the corresponding movement speed of the front decreases. The variation trend of the dissociation front with bottom hole pressure, reservoir temperature, and hydrate saturation is consistent with the conclusion in Table 2. With the decrease in bottom hole pressure and the increase in reservoir temperature, the motion speed of the dissociation front gradually increases (Figure 2a,b). The dissociation front, especially the sensitivity function of hydrate saturation, moves rapidly with the decrease in hydrate saturation (Figure 2c). It should be pointed out that when we discussed the relationship between the dissociation front and time, we assumed that the hydrate saturation remained unchanged, which is contrary to the fact, and this is also a limitation of the model in this paper. Consideration of dissociation dynamics will be added to the model in the subsequent research.

Figure 3 shows the pressure distribution curves near the well under different reservoir and well conditions. As shown in the figure, the hydrate reservoir is divided into two zones by the dissociation front, the decomposed zone and the undecomposed zone, and the pressure distribution in the two zones is quite different. In the decomposed zone, the pressure decreases rapidly from the dissociation pressure to the bottom hole pressure, and the pressure drop funnel is deep. In the undecomposed zone, the pressure slowly recovers to the original reservoir pressure, and the pressure drop funnel is relatively gentle. The pressure distribution curves of different reservoirs and well conditions have certain differences. The pressure drop funnel in the decomposed zone deepens as the bottom hole pressure decreases, while the pressure distribution curve in the undecomposed zone almost coincides (Figure 3a). This is due to the gradual increase in the production pressure differential as the bottom hole pressure decreases, which in turn increases the driving force of hydrate dissociation. As a result, the rate of hydrate dissociation is accelerated, more fluid is extracted from the reservoir, and the rate of pressure reduction in the decomposed zone is accelerated. Reservoir temperature and hydrate saturation have a significant influence on the pressure distribution of the hydrate reservoir. As the reservoir temperature increases, the pressure drop funnel in the decomposed zone deepens, while the pressure drop funnel in the undecomposed zone decreases (Figure 3b). This is because the reservoir temperature provides an energy supply. The higher the reservoir temperature, the faster the hydrate dissociation in the decomposed zone, and the deeper the pressure drop funnel formed. In addition, it can be seen from Table 2 that the dissociation pressure increases with an increase in reservoir temperature. Therefore, the higher the reservoir temperature, the closer the dissociation pressure is to the original reservoir pressure, and the shallower the pressure drop funnel in the undecomposed zone. On the contrary, as hydrate saturation decreases, less hydrate is available for dissociation, the pressure drop funnel in the decomposed zone decreases, and the corresponding pressure drop funnel in the undecomposed zone increases (Figure 3c). This is because the lower the hydrate saturation in the decomposed zone, the farther the dissociation front moves (Table 2). At this time, under the same pressure difference (the difference between the dissociation pressure and the bottom hole pressure), a wider decomposed zone corresponds to a smaller pressure drop funnel. Correspondingly, when the hydrate saturation is lower, the disturbance in the undecomposed zone becomes larger, and the pressure drop funnel increases.

Figure 4 shows the temperature distribution curve near the well under different conditions. As mentioned above, the hydrate reservoir is divided into two zones by the dissociation front, and the temperature changes in the two zones are quite different. In the decomposed zone, the lowest temperature appears at the wellbore, which may be due to the fact that only heat convection is considered in this paper and the heat from deep formation cannot be transferred to the well circumference in time. From the wellbore to the dissociation front, the temperature rises rapidly, and the temperature gradient is greatest at the dissociation front. This is because hydrate dissociation is an endothermic reaction, and the dissociation front is where the decomposition is most intense and absorbs a large amount of heat. In the undecomposed zone, the temperature slowly rises from the dissociation temperature to the original reservoir temperature.

As can be seen from Figure 4a, bottom hole pressure has a great influence on the temperature change in the decomposed zone. When the bottom hole pressure is 10 MPa, 8 MPa, and 5 MPa, the temperature gradient in the decomposition zone is about 0.087 °C/m, 0.169 °C/m, and 0.378 °C/m, respectively, showing a gradually increasing trend. As noted before, the lower bottom hole pressure causes more hydrates to reach equilibrium conditions and begin to decompose, thus absorbing more heat, resulting in a faster cooling of the decomposed zone. Similarly, the variation in reservoir temperature will also affect the temperature distribution in the decomposed zone. When the reservoir temperature is 15 °C, 12 °C, and 10 °C, respectively, the temperature gradient in the decomposition region shows a decreasing trend (Figure 4b). This may be due to the fact that when the initial reservoir temperature is low, hydrate dissociation absorbs a large amount of heat, resulting in a substantial decrease in local reservoir temperature, causing hydrate self-protection effects, secondary formation, and other problems, thus slowing down hydrate dissociation [13]. The change in hydrate saturation also affects the temperature distribution around the well. As hydrate saturation gradually decreases, less hydrate is available for dissociation, and the heat consumed gradually decreases, thus reducing the temperature gradient in the decomposed zone (Figure 4c).

For the undecomposed zone, the temperature distribution varies under different bottom hole pressures (Figure 4a). It can be seen from Table 2 that the change in bottom hole pressure has no effect on the dissociation temperature. Whereas, the lower the bottom hole pressure, the farther the dissociation front moves. This indicates that the lower the bottom hole pressure, the more dissociation heat is required from the reservoir and thus, the lower the temperature in the undecomposed zone. Under different reservoir temperatures, the temperature distribution curves in the undecomposed zone are almost parallel, with differences only in numerical values. All curves converge from dissociation temperature to reservoir temperature (Figure 4b). The effect of hydrate saturation on the temperature distribution in the undecomposed zone is similar to that of bottom hole pressure. The dissociation temperature does not change with hydrate saturation. The lower the hydrate saturation, the farther the dissociation front moves (Table 2). This indicates that when hydrate saturation is low, the undecomposed zone is disturbed more (more heat needs to be provided for hydrate dissociation), thus the temperature drop gradient is larger.

4.2. Analysis of Factors Affecting Production Capacity of NGH Depressurization

In this paper, the simulated natural gas production was verified by the data reported in the literature. Ye et al. reported that in the Shenhu area of the South China Sea, the natural gas production was about 2.87 × 10⁴ m³/day. The thickness of the gas-bearing zone in this area was about 89 m, and the length of the horizontal well was about 300 m [25]. As our research shows, the predicted gas production is of the order of 10⁻⁵ m³/s per 1 m² of well. This is about 2.31 × 10⁴ m³/day, based on reservoir and well data in the Shenhu area. Considering the difference in reservoir permeability and mining conditions, the simulated gas production is within an acceptable range.

4.2.1. Influence of Bottom Hole Pressure on Productivity

The variation in methane gas production (unit length) with time under different bottom hole pressures is shown in Figure 5a. As expected, gas production decreases with the square root of time. Methane gas production is a sensitive function of bottom hole pressure. As bottom hole pressure decreases, the pressure gradient near the well increases (Figure 3a). Accordingly, the driving force of hydrate dissociation increases, and the movement speed of the hydrate dissociation front accelerates (Figure 2a). Therefore, the higher the peak gas production of production wells. Although lowering the bottom hole pressure can improve the gas production rate, too low a bottom hole pressure drop will lead to a large and rapid decomposition of hydrate, and the temperature near the wellbore will drop significantly (Figure 4a), resulting in problems such as the self-protection effect and secondary formation of hydrate, which hinder further decomposition of hydrate. Therefore, it is necessary to calculate and select the optimal bottom hole pressure in actual depression-mining so that gas production can meet industrial demand without causing the self-protection effect and secondary formation of hydrate due to excessive gas production [13]. Figure 5b shows the variation in the water production (unit length) with time under different bottom hole pressures. The trend of the water production curve is the same as that of gas, but it is 4 orders of magnitude less numerically. There may be two reasons for the low water yield. One is that there is a big difference in the viscosity between the two phases of gas and water, and the gas seepage is much higher than that of water. The other is that most of the hydrate decomposing water is trapped in the formation.

4.2.2. Influence of Reservoir Temperature on Productivity

Figure 6 shows the relationship between gas and water production (unit length) over time at different reservoir temperatures. As shown in the figure, both gas and water production decrease with the square root of time, and both are sensitive functions of reservoir temperature. The higher the reservoir temperature, the more gas and water production gradually increases. In the one-dimensional model presented here, the heat required for hydrate dissociation must be supplied by the hydrate reservoir. A higher initial temperature means that more heat is available for hydrate dissociation. Consequently, more hydrates begin to decompose, and thus more decomposition products are produced. Therefore, reservoir temperature becomes an important control parameter. It should be emphasized that for real gas hydrate reservoirs, heat can also be transferred from the side, which will significantly affect the production process of natural gas.

4.2.3. Influence of Hydrate Saturation on Productivity

Figure 7a shows the variation curve of gas production (unit length) when hydrate saturation values are 0.5, 0.3, and 0.2. As can be seen from the figure, gas production increases with a decrease in hydrate saturation. This is because when the hydrate saturation is small, the absolute permeability of the reservoir is relatively large, and thus the corresponding gas production is large. Unlike gas production, water production does not always increase with a decrease in hydrate saturation. As can be seen from Figure 7b, as the hydrate saturation decreases, the water production increases first and then decreases. It should be pointed out that when the hydrate saturation is relatively high (S_H = 0.5), the water production is almost zero. As mentioned above, when the hydrate saturation reaches a certain value, such as 0.5, the absolute permeability of the reservoir is limited. Additionally, the seepage velocity of gas is much higher than that of water. Thus, in the limited permeability, gas preferentially produces, and most of the water remains in the reservoir, resulting in almost zero water production. With the decrease in hydrate saturation (S_H = 0.3), the absolute permeability of the reservoir increases, resulting in an increase in water production. When the hydrate saturation decreases further (S_H = 0.2), although the absolute permeability of the reservoir is relatively large, the water production, however, is reduced due to the relatively small amount of final decomposed water.

In the exploitation of natural gas hydrates, methane gas production has always been the focus of attention. According to the above results, methane gas production is a sensitive function of bottom hole pressure, reservoir temperature, and hydrate saturation. By lowering the bottom hole pressure and increasing the reservoir temperature, the pressure and temperature gradient in the decomposed zone will increase. Both of these increases will improve the driving force of hydrate dissociation and push the dissociation front to move farther, thus promoting the development of methane gas. Nevertheless, with the increase in hydrate saturation in a certain range, the gas production of methane decreases, although the pressure and temperature gradient in the decomposed zone increase. This may be due to the fact that methane gas production is also affected by the absolute permeability of the reservoir. The increase in hydrate saturation will increase the driving force of hydrate dissociation, but it will also lead to a decrease in the absolute permeability of the reservoir, and the latter has a more obvious inhibitory effect on methane gas extraction. As a result, methane gas production is reduced.

5. Conclusions

In this paper, a one-dimensional numerical model of gas hydrate depressurization was established, in which convective heat transfer and gas–water two-phase seepage were taken into account. The movement law of the hydrate dissociation front and the evolution process of temperature and pressure near the well were analyzed, and the influence mechanisms of bottom hole pressure, reservoir temperature, and hydrate saturation on hydrate productivity were also discussed. According to the results obtained by the model, the following understandings were obtained:

(1) The dissociation front moves to the depth with time, and the moving speed slows down gradually. The dissociation front movement speed increases slowly with the decrease in bottom hole pressure and reservoir temperature, yet increases sharply with the decrease in hydrate saturation.

(2) The hydrate reservoir is divided into two zones by the dissociation front: the decomposed zone and the undecomposed zone. In the decomposed zone, the pressure decreases rapidly from the dissociation pressure to the bottom hole pressure, and the pressure drop funnel is deeper. In the undecomposed zone, the pressure slowly recovers from the dissociation pressure to the original reservoir pressure, and the pressure drop funnel is relatively gentle. Lower bottom hole pressure, higher reservoir temperature, and higher hydrate saturation all lead to an increased pressure gradient in the decomposed zone.

(3) Bottom hole pressure, reservoir temperature, and hydrate saturation all affect the temperature distribution in the decomposed zone. The temperature gradient in the decomposed zone is proportional to the reservoir temperature and hydrate saturation and inversely proportional to the bottom hole pressure.

(4) Gas and water production both increase with the decrease in bottom hole pressure, yet the latter is 4 orders of magnitude less than the former in numerical value, which may be due to the large difference in the viscosity of gas–water two phases, and the gas seepage speed is much larger than that of water.

(5) Both gas and water production are sensitive functions of reservoir temperature. A higher initial reservoir temperature can provide more heat for hydrate dissociation, and the corresponding decomposition products are higher.

(6) The gas production increases with the decrease in hydrate saturation, while the water production increases first and then decreases with the decrease in hydrate saturation. This indicates that when the hydrate saturation is low, although the absolute permeability of the reservoir is relatively large, the water production is still low even if the hydrate is completely decomposed due to the small hydrate content.

Based on the previous model, this paper has made some progress in considering the gas–water two-phase seepage in the decomposed zone, whereas the model itself still has some shortcomings. For instance, the model in this paper is a simple one-dimensional model, and only heat convection is considered in the thermal model. Additionally, the hydrate saturation in the decomposed zone is regarded as a fixed value when calculating the permeability, and its variation with time is not considered. These model defects need to be further improved in future studies.