Numerical Simulation of Improved Gas Production from Oceanic Gas Hydrate Accumulation by Permeability Enhancement Associated with Geomechanical Response

Abstract

In the Shenhu Area of the South China Sea, although some numerical studies are conducted on the gas production at well SHSC-4, the geomechanical responses have not been taken into account, and the associated impact of permeability enhancement on gas production has not been thoroughly investigated. In this study, pTOUGH+HYDRATE V1.5 coupled with the RGMS is applied to account for geomechanical responses. Based on actual geological conditions, the reservoir model has five layers: the hydrate-bearing layer (HBL), the three-phase layer (TPL), the free gas layer (FGL), the overburden, and the underburden. The numerical results match the trial production data, validating the numerical model. The analysis shows that gas production from the FGL contributed the most (72.17%) to the cumulative gas production (V_g), followed by the TPL (23.54%) and the HBL (4.29%). The cumulative water-to-gas ratio (R_wgT) gradually decreased during gas production, with the HBL exhibiting the highest value. Permeability enhancement can improve gas production, with the FGL being the most responsive to such enhancement. It increased V_g by 87% and reduced R_wgT to 85%. To achieve more realistic production schemes and better enhance energy recovery, it is advisable to conduct numerical investigations that incorporate geomechanical considerations due to the intricate nature of hydrate-bearing sediments.

1. Introduction

1.1. Background

Natural gas hydrates (NGHs) are commonly referred to as combustible ice, which is a non-stoichiometric cage crystalline compound [1], because of their ice-like appearance and ability to be burned. NGH, being an unconventional source of clean energy, produces the least amount of CO₂ per unit of energy [2]. The potential reserves of hydrated gas are over 1.5 × 10¹⁶ m³ and are widely distributed throughout the earth [3], with over 230 hydrate deposits discovered globally in ocean floors and permafrost zones. To extract methane gas from hydrate reservoirs, the in-situ equilibrium condition (high pressure and low temperature) of NGH must be broken, allowing it to decompose and be produced as fluid. Four methods can be used for gas recovery from hydrate-bearing sediments, including depressurization [4,5], thermal stimulation [6], inhibitor injection [7], and CO₂–CH₄ replacement [8,9,10,11]. National programs exist in many countries to research and produce natural gas from gas hydrate deposits in order to discover the commercialization possibility of methane hydrate resources, leading to various studies on the Alaska North Slope [12,13,14,15], the Mallik site in Canada [4,16,17], the Black Sea [18,19], the Krishna–Godavari basin in India [20,21], the Ulleung basin in Korea [22,23], the Nankai Trough in Japan [24,25,26], and the South China Sea [27,28,29] and Qilian Mountain [30] in China.

In general, both marine and terrestrial gas hydrate accumulations are targeted for exploitation. Extracting methane from hydrates can provide a significant global energy supply as long as the operation is controlled to prevent any leakage. Concerns about potential environmental risks are increasing, particularly regarding methane escape to the seafloor during hydrate exploitation in marine region [31]. In addition to production leakage, the warming atmosphere can disrupt the stability of hydrate fields, leading to the release of sequestered methane into the sediments and soils above. The methane, along with methane-derived carbon, can further contribute to greenhouse warming when it reaches the atmosphere [32]. If this happens in the terrestrial region, it has the potential to cause a catastrophic disaster, even with only around 1% of the global gas-in-place estimated to be in the permafrost-associated NGH [33]. This topic has sparked interest in studies that aim to quantify natural gas hydrates from ice [34] and understand how gas hydrate-bearing permafrost sediments respond to changes in environmental temperature [35]. While both NGH exploitation and global warming can have environmental effects, the focus of this study is to evaluate the production potential of offshore hydrate accumulations in the South China Sea.

The main focus of this study is on the Shenhu area located in the northern part of the South China Sea. Between the years 2007 and 2016, the Guangzhou Marine Geological Survey (GMGS) conducted three gas hydrate drilling expeditions, namely, GMGS1, GMGS3, and GMGS4 [36]. During GMGS3, the team drilled a total of 23 sites for logging-while-drilling (LWD) and 4 sites for core drilling. The results indicated that the average hydrate saturation ranged from 13.7% to 45.2% [36]. Among these drilling sites, the first gas hydrate production site, W17, was selected [29]. This was a significant accomplishment since hydrates in the Shenhu area exist in clayey silt sediments, which pose production challenges due to their low permeability and high levels of clays. The W17 well test demonstrated the possibility of production in such a location, which accounts for 90% of total hydrate reservoirs [37]. In 2017, the China Geological Survey set a new world record for gas hydrate production by performing a successful offshore methane hydrate production test in the Shenhu Area. The team used a single vertical well and depressurization for a production period of 60 days, resulting in a total gas production of 3.09 × 10⁵ ST m³, averaging 5.15 × 10³ ST m³/d [29]. This achievement marked the most significant amount of gas production and the longest production period in history. Recently, the China Geological Survey performed another production test using a horizontal well for the first time in the Shenhu Area, achieving a 30-day continuous gas production process. The total gas production reached 8.614 × 10⁵ ST m³, averaging 2.87 × 10⁴ ST m³/d [27]. The gas production rates in both production tests remain significantly below 5.00 × 10⁵ ST m³/d, which is necessary for the commercial exploitation of NGHs [1].

The commercialization of methane hydrate resources is a crucial challenge, necessitating the need for enhanced gas recovery from methane hydrate deposits and leading to the exploration of various methods to increase gas production. The second production test in the Shenhu area indicates that horizontal wells have been found to enhance gas production by increasing the contact area between the well and the reservoir. However, the cost of horizontal wells is 1.5 to 2.5 times more than vertical wells [38]. In addition to horizontal wells, permeability enhancement (e.g., hydraulic fracturing), allowing gas to flow more freely towards the wellbore, might be a cost-effective option for methane hydrate reservoirs to stimulate gas production, resulting in a significant increase in gas recovery. This enhanced gas recovery can maximize the productivity of wells and boost overall production rates, ultimately paving the way for the commercialization of methane hydrate resources. Previous studies in the eastern Nankai Trough have demonstrated that permeability anisotropy and permeability improvement can enhance gas recovery [24,39]. Therefore, it appears feasible to improve gas recovery in the South China Sea via permeability enhancement as well.

Conducting field tests to determine gas recovery via different production schemes is impossible, thus necessitating the use of a numerical simulator to simulate production performance. The most commonly used simulator for studying production behavior in gas-hydrate-bearing sediments is TOUGH+HYDRATE (T+H). T+H is capable of describing mass and heat flows, hydrate formation and dissociation, and the inhibitor effect [40]. When compared to other simulators available for hydrate simulations, T+H demonstrates superior predictive results when assessing production feasibility for permafrost and marine hydrates [41]. T+H has been utilized to identify gas recovery from low-permeability hydrate reservoirs via depressurization [42], as well as natural gas hydrate reservoirs in the eastern Nankai Trough [26]. Some studies have also incorporated a geomechanical simulator to analyze geomechanical responses [43,44,45]. By utilizing these simulators, the effectiveness of permeability enhancement can be determined.

In order to evaluate the potential for gas production resulting from permeability enhancement in the South China Sea, various numerical simulations have been carried out. In the Liwan 3 Area, Zhang et al. [46] focused on a methane hydrate reservoir with the conditions of low permeability located at site LW3-H4-2. Their approach involved the use of an enlarged, highly permeable well wall, which was proposed as a means of promoting gas production. At site SH2 in the Shenhu Area, Li et al. [47] and Sun et al. [48] suggested the stratification split grouting foam mortar method and hydraulic fracturing, respectively, as methods to improve permeability for production enhancement. The latter study demonstrated the influences of horizontal and vertical fractures. Despite the different methods used for permeability enhancement, Yu et al. [49] conducted a study showing that gas production can be significantly enhanced by the complex mechanisms associated with permeability enhancement at well SHSC-4. These models were based on various geological conditions in the South China Sea and explored the impact of different parameters, including permeability, on gas recovery.

Despite the progress made in enhancing gas production from methane hydrate deposits, the effect of permeability associated with geomechanical responses on gas production enhancement has not been thoroughly investigated. Previous studies have predominantly relied on simplistic adjustments of formation porosity and permeability based on pressure and temperature, with a primary focus on understanding the interplay between system flow and thermal aspects. Unfortunately, these studies have overlooked the crucial factors of media deformations and changes in stress fields. Geomechanics is a critical factor in understanding the behavior of hydrate behavior because of the state and phase changes, leading to changes in its porosity, permeability, and flow characteristics [4,50,51]. Specifically, changes in stress and pressure can cause the rock to compact or expand, affecting the pore structure and connectivity. These changes in permeability can significantly impact the flow of gas within the reservoir, influencing production rates and overall productivity.

Before commercial gas production from hydrate deposits can be developed, it is essential to address and understand the geomechanical response of hydrate-bearing sediments. Particularly, potential wellbore instability and casing deformation are significant concerns that must be tackled [52]. The application of mechanical loads to sediments containing hydrates can cause hydrate dissociation and result in a substantial increase in pressure. This pressure increase can negatively impact the wellbore assembly, the hydrate-bearing sediments, and the surrounding formations, posing risks to their integrity and stability [50]. To prevent failures like tensile or buckling failure in the well assembly, it becomes crucial to implement appropriate well design based on numerical simulations accounting for the geomechanical response [53].

To evaluate the production potential in the South China Sea, it is crucial to have a comprehensive understanding of the intricate interactions between permeability, geomechanical responses, and gas production enhancement. To gain a more accurate and holistic understanding of these phenomena, it is imperative to consider the influence of media deformations and stress field changes. Only a geomechanical model can accurately capture the mechanical behavior of subsurface rocks, including (a) the deformation and potential failure of the reservoir media and well assembly; (b) the effects of changes in stresses, fluid pressures, and temperatures on flow properties like porosity and permeability that control production; and (c) the interdependence of system flow, thermal, and geomechanical properties. Hence, the accurate determination of (a) the impact of permeability enhancement associated with geomechanical responses on hydrate development and (b) the effectiveness of permeability enhancement can be achieved.

1.2. Targeted Accumulation

This study focuses on the oceanic hydrate accumulation at well SHSC-4 in the Shenhu area (Figure 1), which is located in the northern South China Sea between the Xisha Trough and Dongsha Islands and is structurally part of the Baiyun Sag of the Zhuer Depression [29]. The seafloor exhibits a complex topography, including various features such as eroded channels, sea valleys, sea mounts, steep slopes, reverse slopes, plateaus, alluvial fans, and slide fans [54]. The continental slope in the northern South China Sea displays both passive and active continental features due to interactions between the Eurasian, Pacific, and Indochina plates. The reservoir lithology is composed of clayey silt, with montmorillonite and illite constituting over a quarter of the total minerals [29]. Based on the logging interpretation and core analysis of SHSC-4, the overlying formation extends from the sea bottom to a depth of 1495 m (201 m below seafloor (mbsf), with a water depth of 1266 m). The NGH system at the test site comprises three intervals, with the first interval ranging from 1495 to 1530 m (201–236 mbsf), and its pore space is filled with solid NGH and liquid water. The second interval spans from 1530 to 1545 m (236–251 mbsf), and its pore space is filled with three phases: solid NGH, free hydrocarbon gas, and liquid water. The third interval spans from 1545 to 1572 m (251–278 mbsf), and its pore space is filled with two phases: free hydrocarbon gas and liquid water.

1.3. Objectives

The objective of this study is to assess the technical feasibility of gas production and production enhancement by the numerical simulation of oceanic hydrate accumulation in the Shenhu Area of the South China Sea. Specifically, the study evaluates the impact of increased permeabilities in the near-wellbore region on gas production associated with the geomechanical responses. The analysis considers the properties and conditions obtained at well SHSC-4; studies the gas production behavior at well SHSC-4 using a single vertical well by depressurization; assesses the contribution of different layers to gas production; and investigates the feasibility of permeability enhancement, the potential of increased gas production, and the geomechanical response of the geological system during production. To validate the accuracy of the numerical model used in the study, a 60-day simulation result is compared against the actual trial production test conducted at well SHSC-4. The impact of permeability enhancement is evaluated by quantifying the ratios of the cumulative gas production and cumulative water-to-gas ratio, comparing them to the original conditions.

2. Methodology

2.1. Coupled Numerical Simulators

Two parallel simulators, pTOUGH+HYDRATE (pT+H) V1.5 and the Reservoir GeoMechanics Simulator (RGMS) [55,56], based on an MPI (message passing interface) were utilized in this study to simulate the coupled processes of flow, thermal, chemical, and geomechanical processes associated with gas production induced by depressurization. pT+H V1.5 is a parallel version of TOUGH+HYDRATE V1.5 used for conducting numerical investigations of flow and thermal behavior in hydrate-bearing geologic media. The code uses the integral finite difference method (IFDM) [57,58] for space discretization and the Newton–Raphson (NR) iteration for a fully implicit solution. The computational domain is decomposed into subdomains for efficient workload distribution to multiple processes. The RGMS is a parallel geomechanical simulator that accurately characterizes deformations and stresses in subsurface systems by employing parallel strategies in conjunction with the finite element method (FEM). It has the capability to handle grids with Cartesian (2D and 3D) and cylindrical (2D) coordinates. Additionally, the simulator can be used either as a standalone or coupled with T+H, pT+H, and other simulators describing flow and thermal behavior. To explore the geomechanical responses to gas production from hydrate reservoirs, pT+H V1.5 and the RGMS are coupled via the fixed-stress split iterative scheme. The coupling scheme was validated in previous studies [55,56].

2.2. Governing Equations

2.2.1. Flows of Fluid and Heat

The mass and energy balance equation includes the accumulation, flux, and source/sink terms through a control volume in a porous medium [40], which is expressed as

$$\frac{d}{dt}{\displaystyle {\int }_{V_n}{M}^{\kappa }\mathrm{d}V=}{\displaystyle {\int }_{{\Gamma }_n}{F}^{\kappa }\cdot n\mathrm{d}A+}{\displaystyle {\int }_{V_n}{q}^{\kappa }\mathrm{d}V},$$

(1)

where $t$ is the time, $V_n$ is the volume of subdomain n, $\mathrm{d}V$ is the differential volume, $M^{\kappa }$ is the accumulation of component $\kappa$ in terms of mass or heat, ${\Gamma }_n$ is the surface area of subdomain n, $\mathrm{d}A$ is the differential surface area, $F^{\kappa }$ is the flux vector of component $\kappa$ in terms of mass or heat, $n$ denotes the inward unit normal vector, and $q^{\kappa }$ is the source/sink term of component $\kappa$.

For component $\kappa$, the mass accumulation $M^{\kappa }$ is calculated by

$$M^{\kappa }={\displaystyle \sum _{\beta =A,G,I,H}\varphi S_{\beta }}{\rho }_{\beta }{X}_{\beta }^{\kappa }, \kappa =w,m,i,$$

(2)

where $\beta$ indicates phase (A refers to aqueous, G refers to gaseous, I refers to solid ice, and H refers to solid hydrate), $\kappa$ indicates component (w refers to H₂O, m refers to CH₄, and i refers to water-soluble inhibitor), $\varphi$ represents the porosity, $S_{\beta }$ represents the saturation of phase $\beta$, ${\rho }_{\beta }$ represents the density of phase $\beta$, and $X_{\beta }^{\kappa }$ represents the mass fraction of component $\kappa$ in phase $\beta$.

The heat accumulation $M^{\theta }$ comprises contributions from both the rock matrix and all the phases and is expressed by

$$M^{\theta }={\displaystyle {\int }_{T_0}^T(1-\varphi ){\rho }_R{C}_{R}dT+{\displaystyle {\sum }_{\beta =1,\cdots ,N_{\beta }}\varphi S_{\beta }{\rho }_{\beta }{U}_{\beta }}+Q_{diss}},$$

(3)

where ${\rho }_R$ is the density of the rock, C_R is the heat capacity of the dry rock, T is the temperature, and $U_{\beta }$ is the specific internal energy of phase $\beta$. The energy change of hydrate dissociation Q_diss is

$$Q_{diss}=\begin{array}{c}\\ \end{array}\{\begin{array}{c}\Delta \left(\varphi {\rho }_H{S}_H\Delta H^0\right) \mathrm{for} \mathrm{equilibrium} \mathrm{dissociation}\\ Q_H\Delta H^0 \hspace{1em}\hspace{1em}\mathrm{for} \mathrm{kinetic} \mathrm{dissociation}\end{array},$$

(4)

where $\Delta ()$ denotes the quantity change during a given timestep, and $\Delta H^0$ is the specific enthalpy of hydrate dissociation, which can be calculated by Equation (5) from [59]:

$$\Delta H^0=C_f\left(C_1+C_2/T\right),$$

(5)

where the conversion factor C_f is 33.72995 J∙gmol∙kg⁻¹∙cal⁻¹,

$$C_1=\{\begin{array}{l}13521\\ 6534\end{array}, \mathrm{and} C_2=\{\begin{array}{l}-4.02\hspace{1em}\mathrm{for}\hspace{1em} 0 °\mathrm{C}<T_c\le 25 °\mathrm{C}\\ -11.97 \mathrm{for} -25 °\mathrm{C}<T_c\le 0 °\mathrm{C}\end{array}.$$

(6)

It is important to note that hydrate dissociation can be treated in two ways: (a) a chemical equilibrium reaction (the hydrate is considered a thermodynamic state of the CH₄ and H₂O system) or (b) a kinetic reaction (the hydrate is regarded as a distinct component) [60]. The former was employed in this study. Interested readers can refer to Moridis [40,61] for in-depth information on the specifics of the two models and the associated thermodynamics.

The mass flux of component $\kappa$ (e.g., H₂O, CH₄, and inhibitor) contributed by the aqueous and gaseous phases is defined as

$$F^{\kappa }={\displaystyle \sum _{\beta =A,G}{F}_{\beta }^{\kappa }}={\displaystyle \sum _{\beta =A,G}{X}_{\beta }^{\kappa }{F}_{\beta }}, \kappa =w,m,i.$$

(7)

Following Darcy’s law, the mass flux of phase $\beta$ is defined as

$$F_{\beta }=-k\frac{k_{r\beta }{\rho }_{\beta }}{{\mu }_{\beta }}\left(\nabla P_{\beta }-{\rho }_{\beta }g\right),$$

(8)

in which, for phase β, $k_{r\beta }$ is the relative permeability, ${\mu }_{\beta }$ is the viscosity, and $P_{\beta }$ is the fluid pressure. Additionally, k is the absolute permeability tensor, and g is the gravity vector.

The heat flux is defined as

$$F^{\theta }=-\overline{k_{\theta }}\nabla T+{\displaystyle \sum _{\beta =A,G}{h}_{\beta }{F}_{\beta }},$$

(9)

where $h_{\beta }$ is the specific enthalpy of phase $\beta$, and $\overline{k_{\theta }}$ is the composite thermal conductivity considering the properties of medium and fluid.

2.2.2. Geomechanics

The quasi-static momentum conservation equation is based on the underlying assumption that there is always equilibrium between the fluid and rock, which is calculated by

$$\nabla \cdot \sigma +{\rho }_{b}g=0,$$

(10)

where $\sigma$ is the total stress tensor. The bulk density ${\rho }_b$ is found using

$${\rho }_b=\left(1-\varphi \right){\rho }_R+\varphi {\rho }_f,$$

(11)

in which ρ_f is the saturation-weighted fluid density, which is calculated using

$${\rho }_f={\displaystyle \sum _{\beta }{S}_{\beta }}{\rho }_{\beta }.$$

(12)

In accordance with convention, the equations below always consider tensile stress as positive. The relationship that defines the stresses associated with the rock skeleton is expressed as

$$\sigma ={\sigma }^{\prime }-\alpha I{P}_t=C:\epsilon -\alpha I{P}_t,$$

(13)

where C is the elasticity tensor, ${\sigma }^{\prime }$ is the effective stress tensor, I is the identity matrix, and $\epsilon$ is the strain tensor. Biot’s coefficient α [62] is defined as

$$\alpha =1-\frac{K_{dr}}{K_s},$$

(14)

where K_s is the skeletal grain modulus, and K_dr is the drained bulk modulus. The average mobile fluid pressure P_t is calculated using

$$P_t=\frac{{\displaystyle \sum _{\beta }{S}_{\beta }{P}_{\beta }}}{{\displaystyle \sum _{\beta }{S}_{\beta }}}.$$

(15)

By utilizing the assumption of infinitesimal deformation, the strain tensor is determined using

$$\epsilon =\frac{1}{2}\left[\nabla u+{\left(\nabla u\right)}^T\right],$$

(16)

in which u is the displacement vector.

2.2.3. The Coupling Method between Geomechanics and Flows of Fluid and Heat

To investigate the geomechanical response during production, pT+H and the RGMS were coupled with the fixed-stress split iterative scheme so that the changes in fluid pressure, temperature, phase saturations, and associated stresses can be found. The accuracy and reliability of the scheme were verified by comparing its numerical results and the analytical solutions of Terzaghi’s problem [63] and McNamee–Gibson’s problem [64,65] in prior research [55].

The scheme operates in a sequential manner as shown in Figure 2, where pT+H solves the problem when the stress field is frozen. The strongly nonlinear equations are solved via the NR iteration, and the porosity is calculated using

$${\varphi }^k={\varphi }^{k-1}+\frac{\alpha \left(1+{\epsilon }_v\right)-{\varphi }^{k-1}}{K_{dr}}\delta P_t{}^{k-1},$$

(17)

in which k is the number of the NR iterations, K_dr is the drained modulus, ε_v is the volumetric strain, and $\delta P_t{}^{k-1}$ is the difference at the NR iteration k and k − 1.

Subsequently, the porosity in the RGMS is calculated using

$${\varphi }^n={\varphi }_0^{}+\alpha {\epsilon }_v+\frac{\left(\alpha -{\varphi }_0^{}\right)\left(1-\alpha \right)}{K_{dr}}\left(P_t-P_{t0}\right),$$

(18)

where subscript n denotes the initial states.

pT+H and the RGMS solve the problem sequentially; do not proceed to the next time step until $\left({\varphi }^n-{\varphi }^k\right)/{\varphi }^n$ is lower than a predetermined tolerance level.

The mechanical properties of the hydrate are influenced by the existence of hydrate, which is more pronounced when S_H has a higher value. However, there is a scarcity of research on how the properties are related to the hydrate phase. Rutqvist and Moridis [50] proposed a standard approach that utilizes linear interpolation equations as follows:

$$K_{dr}=K_{dr0}\left(1-S_H\right)+K_{dr1}{S}_H,$$

(19)

and

$$G=G_0\left(1-S_H\right)+G_1{S}_H,$$

(20)

where G is the shear modulus, subscript 0 denotes S_H = 0, and subscript 1 denotes S_H = 1.

3. Numerical Model

3.1. The Geologic Model

The geological model was developed from a combination of geophysical surveys and an analysis of core samples retrieved from the well. The model consists of five units: the overburden (OB), the hydrate-bearing layer (HBL), the three-phase layer (TPL), the free gas layer (FGL), and the underburden (UB). The computational domain used in both the production and geomechanical studies accurately represents the geological model, encompassing the entire OB to the ocean floor and a portion of the UB that was thick enough to fully account for the necessary heat and water exchange with the reservoir, and provides a true zero-displacement boundary for geomechanical computations.

Combining geophysical surveys and core sample analysis from the well, the geological model was developed (Figure 3), comprising five units: the overburden (OB), the hydrate-bearing layer (HBL), the three-phase layer (TPL), the free gas layer (FGL), and the underburden (UB). In pT+H V1.5 and the RGMS, the same computational domain was utilized to represent the geological model, covering the entire OB to the ocean floor and a thick enough portion of the UB to account for water and heat exchange within the reservoir, and to make sure there was no displacement occurring at the bottom boundary for geomechanical computations.

3.2. Domain Discretization

A 2D cylindrical domain with an outer radius of 300 m, as presented in Figure 4, was utilized to investigate gas production with a single vertical well. The domain was divided into 412 segments radially and 232 segments vertically, leading to 95,584 gridblocks. Gridblocks with different lengths were produced due to the high-resolution radial discretization in the vicinity of the well; specifically, the segment length in the radial direction (Δr) of 0.10 m was used for r_w < r ≤ 1 m, while that of 0.20 m was used for 1 m ≤ r ≤ 21 m. For distances greater than 21 m but less than 300 m (r_max), Δr increased logarithmically for r > 0.20 m. The segment length in the vertical direction (Δz) was 0.5 m in the hydrate accumulation area and was larger in OB and UB. A mesh representation of the domain used in this study is presented in Figure 5a, while a more detailed representation of the grid near the wellbore is given in Figure 5b.

For each element, four equations were designed to account for the mass balance of three components (H₂O, CH₄, and NaCl) as well as the heat balance of the system. Previous research has demonstrated that a chemical equilibrium reaction provides an accurate result during hydrate formation/dissociation in gas production [60]. To account for the scale of the problem, a total of 560,000 equations were formulated, consisting of approximately 380,000 equations in pT+H V1.5 and 180,000 equations in the RGMS. Consequently, the problem sizes necessitated the use of pT+H V1.5 and the RGMS to provide practical solutions.

3.3. Well Description

At the center of the cylindrical domain, the vertical production well was perforated from 201 mbsf to 268 mbsf, covering the HBL, the TPL, and a portion of the FGL (Figure 4). The well was treated as a pseudo-porous medium to use Darcy’s flow concepts within the wellbore, of which the properties are as follows: the vertical permeability was 5 × 10⁻⁹ m², the porosity was 1, the capillary pressure was 0, the irreducible gas saturation was 0.005, and the relative permeabilities had linear relationships with phase saturations. This approach was able to simulate the pressure drop in a steel wellbore, which was validated in a previous study [66]. The bottomhole pressure (P_bh) was 3 MPa [67] at a gridblock above the topmost well gridblock.

3.4. System Properties

Table 1. Properties, conditions, and models used in pT+H V1.5.

Properties, Conditions, Models	Values
Initial pressure at the bottom of TPL	14.93 MPa
Initial temperature at the bottom of TPL	14.82 °C
Gas composition	100% CH4
Initial saturation of HBL	SH = 0.34
Intrinsic permeabilities of HBL	kh = 2.86 × 10−15 m2 = 2.9 mD; kz = kh
$\mathrm{Porosity} \varphi$ of HBL	0.35
Initial saturation of TPL	SH = 0.31, SG = 0.078
Intrinsic permeabilities of TPL	kh = 1.48 × 10−15 m2 = 1.5 mD; kz = kh
$\mathrm{Porosity} \varphi$ of TPL	0.33
Initial saturation of FGL	SG = 0.078
Intrinsic permeabilities of FGL	kh = 7.30 × 10−15 m2 = 7.4 mD; kz = kh
$\mathrm{Porosity} \varphi$ of FGL	0.32
Intrinsic permeabilities of OB	kh = 9.87 × 10−18 m2 = 0.01 mD; kz = kh
$\mathrm{Porosity} \varphi$ of OB	0.10
Intrinsic permeabilities of UB	kh = 9.87 × 10−18 m2 = 0.01 mD; kz = kh
$\mathrm{Porosity} \varphi$ of UB	0.10
Dry thermal conductivity	kθd = 1 W∙m−1∙K−1
Specific heat CR	1000 J kg−1∙K−1
Grain density ρR	2650 kg∙m−3
Composite thermal conductivity model [40]	$\begin{array}{c}\overline{k_{\theta }}=k_{\theta d}+\left(\sqrt{S_A}+\sqrt{S_H}\right)\left(k_{\theta w}-k_{\theta d}\right)\\ \hspace{1em}+\varphi S_I{k}_{\theta I}\hfill \end{array}$
Relative permeability model EPM#2 [40]	$k_{rA}=\max \left\{0,\min \left\{{\left[\frac{S_A-S_{irA}}{1-S_{irA}}\right]}^n,1\right\}\right\}$; $k_{rG}=\max \left\{0,\min \left\{{\left[\frac{S_G-S_{irG}}{1-S_{irA}}\right]}^{n_G},1\right\}\right\}$; $k_{rH}=0$
SirA, SirG, n, nG [69]	0.65; 0.03; 3.50; 2.50
Capillary pressure model [70]	$P_{cap}=-P_0{\left[{\left(S^{*}\right)}^{-\frac{1}{\lambda }}-1\right]}^{1-\lambda }$ $S^{*}=\frac{S_A-S_{irA}}{S_{mxA}-S_{irA}}$
λ, P0, SirA, SmxA of HBLs	0.45; 104 Pa; 0.65; 1.0
Porosity–permeability relationship [71]	$\frac{k}{k_0}=\exp \left[\gamma \left(\frac{\varphi }{{\varphi }_0}-1\right)\right]$
Empirical permeability reduction factor γ [71]	29.0

provides the properties used in pT+H V1.5 that are based on the previous studies [29,49,67,68]. The properties used in the RGMS are based on the previous studies [53,69] as shown in

Table 2. Properties used in RGMS.

Properties	Values
Young’s modulus of HBL	E = 200 MPa at SH = 0; E = 1.4 GPa at SH = 1
Young’s modulus of TPL	E = 200 MPa at SH = 0; E = 1.4 GPa at SH = 1
Young’s modulus of FGL	E = 200 MPa
Young’s modulus of OB	E = 70 MPa
Young’s modulus of UB	E = 200 MPa
Poisson’s ratio of HBL	ν = 0.15
Poisson’s ratio of TPL	ν = 0.15
Poisson’s ratio of FGL	ν = 0.45
Poisson’s ratio of OB	ν = 0.45
Poisson’s ratio of UB	ν = 0.45
Biot’s coefficient	α = 0.99

.

3.5. Initial Conditions

Table 1 presents the temperature and pressure values at the bottom boundary of the TPL, which conform to the geothermal and hydrostatic gradients, respectively, as is typical of most hydrate deposits. Following the initialization process outlined by Moridis and Reagan [66], the simulation with the initial conditions continued until the entire domain reached a steady state without any change in various distributions. For geomechanical simulation, the initial total stress field was obtained by assuming that the initial effective stress field was zero, without requiring the application of the overburden pressure.

3.6. Model Validation

In order to guarantee the prediction accuracy of gas production at well SHSC-4, it is imperative to first validate the numerical model built earlier. As stated in the Introduction section, a trial production test was conducted at well SHSC-4, lasting 60 days, with a total gas production of 3.09 × 10⁵ ST m³ [29]. Using the constructed model, a simulation was performed to replicate the 60-day production, taking into account the geomechanical responses. The simulation resulted in a total gas production of 3.08 × 10⁵ ST m³, as depicted in Figure 6. This successful replication serves as validation for the constructed model.

3.7. Simulations Cases

The base case was based on the original formation information obtained from the first test production at well SHSC-4 (Figure 4). In addition to the base case, this study investigates whether the cumulative gas production can be improved by increasing the permeabilities of the HBL, TPL, and FGL, individually. As illustrated in Figure 7, the permeability was increased to k_s within the region from 0 to the stimulated radius (r_s) within a specific layer. To assess the effectiveness of permeability enhancement, a permeability enhancement ratio (f_k = k_s/k₀) was proposed, in which k₀ is the original permeability of the layer. There are four different values of r_s (0.3 m, 0.5 m, 1 m, and 2 m) and three different values of f_k (2, 4, and 8), combined with three layers, resulting in a total of thirty-six cases.

4. Results and Discussion

In the analysis of the production potential associated with the geomechanical response at well SHSC-4, the monitoring included various parameters related to flow and geomechanics both with and without permeability enhancement, which were obtained from pT+H and the RGMS. The monitored flow-related parameters were pressure (P), temperature (T), gas saturation (S_G), hydrate saturation (S_H), the production rates of CH₄ and H₂O (Q_g and Q_w, respectively), and the cumulative production of CH₄ and H₂O (V_g and M_w, respectively), The water-to-gas ratio was also monitored, both instantaneously (R_wg = Q_w/Q_g) and cumulatively (R_wgT = M_w/V_g). Geomechanics-related parameters monitored were radial and vertical displacements (u_r and u_z, respectively) at key locations. To evaluate the influence of permeability enhancement, the key parameters are V_g and R_wgT. Specifically, more gas and less water are desired after permeability enhancement, so larger V_g and smaller R_wgT values are better.

4.1. Base Case

4.1.1. Fluid Production

Figure 8a shows Q_g produced at the well from the HBL, TPL, and FGL and all layers in the base case. The value of Q_g has an initial peak after production begins, followed by a decline and minor fluctuations within a certain range in the subsequent production period. The initial peak is caused by the rapid dissociation of hydrates near the wellbore region and a subsequent surge in gas production rate after the bottomhole pressure drops. In addition, the free gas in the TPL and FGL contributes to the initial peak of Q_g. The average gas production over the entire production period is 0.074 ST m³/s, which is far below the gas production rate of 0.579 ST m³/s (=5.00 × 10⁵ ST m³/d) required for the commercial exploitation of NGHs [1]. Compared with hydrate deposits in Mount Elbert, Alaska North Slope, where there exists a lag time before substantial gas production [14], the hydrate deposit at well SHSC-4 does not exhibit such a lag phenomenon but instead has the highest gas production rate in the early stage of production, indicating that this class of hydrate deposit is conducive to exploitation. The contribution of each layer to the total gas production rate was ranked from highest to lowest as FGL, TPL, and HBL, indicating that the FGL is the primary source of gas production. This suggests that the FGL is the most important layer for gas production in the studied area.

Figure 8b shows V_g produced at the well from the HBL, TPL, and FGL and all layers in the base case. As shown in the figure, in the later production period, since Q_g fluctuates within a certain range, V_g, which is the integral of gas production over time, shows a nearly linear relationship with time. After 120 days of production, the FGL, TPL, and HBL accounted for 72.17%, 23.54%, and 4.29% of the total cumulative gas production, respectively. This also indicates that the FGL is the most important layer for gas production because it has the highest contribution to the total cumulative gas production. It is anticipated that the FGL will exhibit the most pronounced response to permeability enhancement.

Figure 9 shows R_wg and R_wgT produced at the well from the HBL, TPL, and FGL and all layers in the base case. Apart from directly evaluating production via Q_g and V_g, R_wg and R_wgT can also be used to indirectly characterize production performance. In real practice, more gas and less water are desired, so smaller R_wg and R_wgT values are better. R_wg produced from the HBL, TPL, and FGL and all layers decreased gradually during production. Among the three layers, the HBL has the highest R_wg with the smallest contribution to gas production observed in Figure 8. Moreover, R_wg produced from the HBL is tens of times higher than those from TPL and FGL. If permeability enhancement is carried out within HBL, gas production may increase, while water production may also increase. R_wgT reaches a short-term peak in the first two days of production and then shows a decreasing trend throughout the entire production period. R_wgT produced from the FGL and all layers tends to stabilize in the later period of production. Due to the large R_wgT produced from HBL, the total R_wgT was far higher than those from TPL and FGL.

Figure 10 shows u_r and u_z at key locations in the base case. As P_bh is lower than the pressure of the formation, the reservoir is “squeezed” and moves toward the vertical well in the radial direction, but the compaction is not significant. Although the location with the largest radial displacement occurs at (r, z) = (1 m, −201 m), the absolute value does not exceed 0.01 m when the simulation ends. In the vertical direction, the subsidence at the top of the HBL and the uplift at the bottom are observed. As the gridblock that was set to the bottomhole pressure is closer to the top of the HBL, the subsidence at the top of HBL is more obvious, with a maximum level of no more than 0.08 m. Overall, the displacement within the formation is not significant.

4.1.2. Spatial Distributions

Figure 11 illustrates the spatial evolution of P in the base case. The cylindrical region of pressure drop gradually expands over production time, with those in the HBL and FGL being larger than that in the TPL. This is mainly due to the higher permeability of the HBL and FGL compared to that of the TPL. Additionally, as the hydrate saturation is highest in the hydrate layer and the effective permeability is lower, fluid replenishment becomes difficult, resulting in the most significant pressure drop. Although the cylindrical region of the pressure drops in the TPL and FGL is relatively small during production, the rates of pressure drop are faster than that in the HBL. Particularly, the pressure drops in the FGL and the lower part of the TPL are very significant. In previous studies [49,67], the pressure drop area obtained using T+H was only not significant in the FGL, which differs from the results in this study. It is evident that considering mechanical response leads to different hydrate production mechanisms.

Figure 12 illustrates the spatial evolution of T in the base case. As hydrate dissociation is endothermic, the low-temperature region generally indicates the location of hydrate dissociation. In this figure, the low-temperature region is found in the FGL and the lower part of the TPL, coinciding with the significant pressure drop area and indicating a large amount of hydrate dissociates. Unlike production lasting over 1000 days in previous studies [49,67], the numerical simulation in this study was halted due to the temperature in the low-temperature region reaching 0.01 °C after 120 days of production.

Figure 13 illustrates the spatial evolution of S_H in the base case. Hydrate dissociation occurs in the area where the significant pressure drop is shown in Figure 11 and the low temperature is presented in Figure 12. Hydrates gradually dissociate during production, but the unevenness of hydrate dissociation progress in each layer becomes apparent. The dissociation rate of the HBL is relatively uniform, while the lower part of the TPL undergoes hydrate dissociation, followed by hydrate formation. Moreover, in the FGL, the hydrate forms and dissociates alternatively. This phenomenon may be caused by the Joule–Thomson cooling effect, the capillary effect, the “upstream weighting” approach applied in the simulator [72], and the equilibrium model used in this study. The cyclic process of hydrate formation and dissociation in the FGL ultimately led to the temperature reaching 0.01 °C, resulting in the simulation stopping.

Figure 14 illustrates the spatial evolution of S_G in the base case. The evolution of gas saturation in the base case is also presented. The gas saturation in the HBL and TPL gradually expands, and some gas migrates from the TPL to the HBL. Gas dissociated from hydrate in the lower part of the TPL migrates toward the wellbore radially and toward the HBL vertically at a very slow rate due to the low permeability of the TPL (1.5 mD), resulting in gas accumulation in the lower part of the TPL. In the FGL, a large amount of gas flows into the wellbore because of the significantly lower wellbore pressure compared to the formation pressure and higher formation permeability. The Joule–Thomson effect caused by the rapidly migrating gas may have caused the low temperature in the FGL, meeting the conditions for hydrate generation. Thus, hydrates are formed, hindering the radial gas migration in the FGL and causing the gas to accumulate on the side away from the wellbore.

4.2. Effect of Permeability Enhancement

4.2.1. Fluid Production

In the following tables, the values marked in green are the minimum, while the values marked in red are the maximum.

Table 3. The cumulative gas production with permeability enhancement (V_g) after 120-day production.

Cumulative Gas Production (ST m3)
			rs (m)
			0.3	0.5	1	2
kf	HBL	2	631,851	637,193	642,373	645,883
		4	639,688	648,514	661,047	679,954
		8	646,531	657,912	677,975	711,590
	TPL	2	641,762	648,207	657,112	667,160
		4	654,888	667,240	685,330	702,448
		8	660,201	676,255	698,987	706,541
	FGL	2	687,226	696,737	714,312	731,513
		4	712,884	740,743	811,315	855,335
		8	758,690	788,555	924,427	1,160,649

and

Table 4. The ratios of cumulative gas production with permeability enhancement (V_g) to that in the base case (V_g,0) after 120-day production.

Vg/Vg,0
			rs (m)
			0.3	0.5	1	2
kf	HBL	2	1.019	1.028	1.036	1.042
		4	1.032	1.046	1.066	1.097
		8	1.043	1.061	1.094	1.148
	TPL	2	1.035	1.046	1.060	1.076
		4	1.056	1.076	1.106	1.133
		8	1.065	1.091	1.128	1.140
	FGL	2	1.109	1.124	1.152	1.180
		4	1.150	1.195	1.309	1.380
		8	1.224	1.272	1.491	1.872

present V_g with permeability enhancement and its ratio relative to that in the base case after 120-day production, respectively. For reference, the V_g in the base case (V_g,0) is 620,668 ST m³ following 120 days of production. The permeability enhancement of three different layers (i.e., the HBL, TPL, and FGL) shows that the larger the values of r_s and k_f, the more significant the increase in gas production. When (k_f, r_s) = (8, 2 m), the V_g values, predicted by improving the permeabilities of the HBL, TPL, and FGL, are 711,590, 706,541, and 1,160,649 ST m³, respectively. Compared to the base case, the production was increased by 15%, 15%, and 87% corresponding to the modification in the HBL, TPL, and FGL, respectively. In order to reach higher production, the permeability of the FGL should be enhanced.

Table 5. The cumulative water-to-gas ratio with permeability enhancement (R_wg) after 120-day production.

Cumulative Water-to-Gas Ratio (kg H2O/m3 CH4)
			rs (m)
			0.3	0.5	1	2
kf	HBL	2	3.599	3.624	3.676	3.760
		4	3.683	3.739	3.834	3.939
		8	3.725	3.800	3.920	4.029
	TPL	2	3.373	3.349	3.318	3.282
		4	3.325	3.283	3.226	3.172
		8	3.308	3.256	3.192	3.173
	FGL	2	3.308	3.307	3.289	3.221
		4	3.294	3.268	3.159	2.969
		8	3.234	3.236	3.066	2.406

and

Table 6. The ratios of cumulative water-to-gas ratio with permeability enhancement (R_wgT) to that in the base case (R_wgT,0) after 120-day production.

RwgT/RwgT,0
			rs (m)
			0.3	0.5	1	2
kf	HBL	2	1.268	1.277	1.296	1.325
		4	1.298	1.318	1.351	1.388
		8	1.313	1.339	1.382	1.420
	TPL	2	1.189	1.180	1.169	1.157
		4	1.172	1.157	1.137	1.118
		8	1.166	1.148	1.125	1.118
	FGL	2	1.166	1.165	1.159	1.135
		4	1.161	1.152	1.113	1.046
		8	1.140	1.140	1.080	0.848

present R_wgT with permeability enhancement and its ratio relative to that in the base case after 120-day production, respectively. For reference, the cumulative water–gas ratio (R_wgT,0) is 2.84 kg H₂O/m³ CH₄ following 120 days of production. The results of permeability enhancement in the HBL show that the larger the values of r_s and k_f, the larger the R_wgT. The increase in permeability near the wellbore area in the HBL results in a greater increase in water production than gas production, leading to a larger R_wgT. The results of permeability enhancement in the TPL and FGL show that the larger the values of r_s and k_f, the smaller the R_wgT. When (k_f, r_s) = (8, 2 m), the R_wgT values, determined by permeability enhancement in the HBL, TPL, and FGL are 4.03, 3.17, and 2.41 kg H₂O/m³ CH₄, respectively. The resulting ratios of R_wgT to R_wgT,0 are 1.42, 1.12, and 0.85 when the permeabilities of the HBL, TPL, and FGL are increased, respectively. To reduce the amount of separated water required for unit gas production, the FGL should be treated to enhance its permeability.

The radial and vertical displacements with permeability enhancement are not significant, and therefore, further discussion is not included in this section.

4.2.2. Spatial Distributions

In this section, the figures depicting the spatial distributions of P, T, S_H, and S_G with permeability enhancement after 120 days of production are arranged in a manner where the r_s values increase from left to right and the k_f values increase from top to bottom. The highest V_g and the lowest R_wgT are achieved after increasing the permeabilities of the FGL, and only the spatial distributions with permeability enhancement in the FGL are discussed.

Figure 15 shows the spatial distributions of P with permeability enhancement in the FGL after 120 days of production.

In cases where the value of k_f is small, augmenting r_s has an insignificant effect on the area of pressure drop, and the spatial distributions of P remain relatively unchanged compared to the base case. Conversely, when the value of k_f is large, the increasing r_s results in a narrower area of pressure drop in the upper section of the FGL and the lower section of the TPL, which is closer to the wellbore.

Figure 16 shows the spatial distributions of T with permeability enhancement in the FGL after 120 days of production. When the value of k_f is small, the varying r_s has a negligible effect on the spatial distributions of T compared to the base case, which is comparable to the area of pressure drop presented in Figure 15. However, when k_f equals 4, an increase in r_s results in a shrinkage of the low-temperature area within the FGL, with a gradual shift of the lowest temperature from the upper to the middle section of the FGL. Furthermore, when (k_f, r_s) = (8, 2 m), the low-temperature region within the FGL becomes exceedingly small.

Figure 17 and Figure 18 show the spatial distributions of S_H and S_G with permeability enhancement in the FGL after 120 days of production, respectively. In Figure 17, as the k_f and r_s values increase, the formation and dissociation of gas hydrate transpire in closer proximity to the wellbore, thereby facilitating the production of gas dissociated from gas hydrate. Furthermore, Figure 18 demonstrates that the formation of gas hydrate is less likely to obstruct the flow of gas, resulting in less gas accumulating on the side of the gas hydrate that is farther from the wellbore. These two figures collectively suggest that augmenting k_f and r_s values is more conducive to gas production. This assertion is supported by Table 3 and Table 4, which indicate that larger k_f and r_s values yield higher cumulative gas production.

5. Conclusions

The following conclusions can be drawn from this study:

• To evaluate the effectiveness of permeability enhancement considering the geomechanical responses in the Shenhu area, a coupled simulation using pTOUGH+HYDRATE V1.5 and the RGMS (Reservoir Geomechanics Simulator) is implemented.
• Based on the geophysical surveys and analysis of core samples at well SHSC-4 located in the Shenhu area of the northern South China Sea, the established numerical simulation model is accurate, and the simulation results are highly consistent with the trial production data, ensuring the reliability of the outcomes obtained in this study.
• In the base case, the formation and dissociation of gas hydrates in the free gas layer (FGL) alternate, ultimately resulting in a low-temperature region near 0 °C and leading to the cessation of the simulation after 120 days of production. The cumulative gas production reached 6.2 × 10⁵ ST m³.
• In the base case, the FGL contributes the most to gas production, accounting for 72.17% of the cumulative gas production (V_g), followed by the three-phase layer (TPL), accounting for 23.54% of the cumulative gas production, and the hydrate-bearing layer (HBL) contributes the least, accounting for only 4.29% of the cumulative gas production.
• In the base case, the cumulative water-to-gas ratio (R_wg) from the HBL, TPL, and FGL gradually decreases during the production of gas hydrates. R_wgT from the HBL, which contributes the least to gas production, is the highest, with a value several times those from TPL and FGL.
• In the base case, the gas production obtained without permeability enhancement is insufficient for commercial production. Permeability enhancement can be an option used to increase gas production.
• After increasing the permeabilities of the HBL, TPL, and FGL with the same permeability enhancement ratio (f_k) and the same simulated radius (r_s), the improvement effect of modifying the FGL is the best, with a maximum increase of 87%. The required mass of water separated from a unit of gas is the lowest when applying permeability enhancement in the FGL, with a minimum value of 85% of the original separation mass.
• The results of modifying the FGL show that the higher the degree of permeability enhancement, the deeper the impact of permeability enhancement and the closer the formation and dissociation of gas hydrates are to the wellbore, making it more difficult for gas to be obstructed by the formation of gas hydrates, which is more conducive to production.
• Although permeability enhancement is attempted in this study, it did not extend the production period as the simulation still ends due to low temperature in the FGL. Future research should focus on exploring methods to prevent such low temperatures from occurring in the FGL.
• The results obtained by considering geomechanical responses differ from previous numerical studies that only considered flow and thermal behaviors. This indicates that neglecting geomechanical responses may result in an incorrect natural gas hydrate production scheme. Therefore, future numerical studies should take geomechanical responses into consideration to obtain more realistic results.
• In future work, it is imperative to discover production schemes that effectively mitigate the occurrence of a low-temperature region after 120 days of production, which currently causes disruptions in numerical simulations, thus enabling the extension of the observation period. Moreover, new production schemes combined with permeability enhancement should be explored to facilitate the achievement of production rates that meet the necessary threshold for the commercial exploitation of natural gas hydrates.