Gas Production and Reservoir Settlement in NGH Deposits under Horizontal-Well Depressurization

1 Introduction

Nature gas hydrate (NGH) is an ice-like crystal that mainly exists in Oceanic continental margin and permafrost region [1]. Global reserves of NGH are estimated to be twice those of conventional fossil fuels, making them a promising alternative energy source for the future [2]. Various techniques have been proposed for gas production from NGH reservoirs, including depressurization [3], thermal stimulation [4], chemical inhibitor injection [5], and gas replacement [6]. Among these methods, chemical inhibitor injection has been extensively studied for natural gas pipeline transportation and well drilling, but has received relatively limited attention for NGH production due to its potential environmental risks [7]. The gas replacement method offers dual benefits of CO2 storage and NGH production, along with the added advantage of reservoir repair [8]. CO2 can also enhance production by forming a hydrate cap or through displacement as a mixed gas with N2 [9]. However, the low replacement rate makes it difficult to meet the commercial production demands [10]. Thermal stimulation methods [11,12,13] are designed to supply the heat required to maintain equilibrium temperature during the rapid decomposition of NGH. Nevertheless, field tests have revealed significant energy losses associated with thermal stimulation methods due to the low permeability [14] and thermal conductivity [15] of the NGH reservoirs. In general, depressurization is widely regarded as the most economical method for NGH production, as it does not necessitate additional energy input. Both numerical simulation at the field scale and field trials have shown their potential commercial application [16].

Previously, field trials primarily employed vertical wells to test gas production from the NGH reservoir. To obtain a high gas production rate from vertical wells, a significant pressure differential is applied to accelerate the NGH decomposition rate [17]. However, large pressure drawdown can trigger Joule-Thomson effect combined with the endothermic of NGH decomposition, causing a rapid decline in reservoir temperature and resulting in the reformation of NGH and ice [15]. In comparison, horizontal wells can expand the pressure relief zone in the reservoir, which facilitates gas production from individual wells even at a lower pressure differential. This approach has shown promise, as demonstrated by the second NGH production field trial in the South China Sea [18].

Beyond economic considerations, safety risks represent a critical constraint on NGH exploitation. Marine NGH reservoirs are typically characterized by shallow burial depth and weak diagenesis properties [19]. In porous sediments, NGH often serves to cement skeletal particles as the primary binding agent [20]. The dissociation of NGH weakens the mechanic properties of the reservoir, leading to reservoir settlement, which may cause wellbore damage, submarine landslide [21], and methane leakage incidents [22]. Therefore, it is necessary to understand the geomechanical response characteristics during NGH production to evaluate the stability of disturbed formations. The mechanical properties of NGH reservoirs have been widely studied, revealing that their strength and stiffness depend on the NGH saturation [23], distribution morphology [19] and the characteristic of the sediment skeleton particles [24]. Numerous studies show that hydrates have less influence on the strength of NGH bearing sediments when the NGH saturation is less than 10% [25], then its strength increases significantly with rising hydrate saturation [26]. Microscopic studies indicate that this phenomenon is attributed to the distribution morphology of the hydrate in porous media. Hydrate is more likely to fill the pores space when its saturation is below 10% [21], then gradually cements the sediment particles with continues increase in saturation [27]. The earth pressure coefficient of NGH is mainly controlled by the sediment skeleton particles at high stresses [28]. Hydrate exploitation leads to the deterioration of the mechanical properties of the reservoir, and the reservoir is bound to settle under the stress of the overlying formation. Several Studies show that the maximum strain occurs at the reservoir boundary directly below the well, with values reaching up to 0.6% [29]. These deformations and settlements cause uneven horizontal in-situ stress on hydrate-bearing reservoirs, which will cause wellbore wall instability during hydrate production [30].

Several software packages are available for the numerical simulation of hydrate reservoirs, such as TOUGH+Hydrate [31] and COMSOL [32]. TOUGH+Hydrate can be coupled with FLAC3D to simulate the geomechanical behaviors of hydrate reservoirs. However, mesh generation remains complex, the numerical simulation of hydrate reservoirs, most studies focused on the geomechanical response characteristics in the NGH production process by vertical wells with field scale [33,34,35]. Seldom considers the impact of heterogeneous hydrate distribution on reservoir settlement during NGH production when using horizontal wells. COMSOL Multiphysics is a simulation software that enables the formulation of partial differential equations to model hydrate formation and decomposition [36]. But the numerical instability issues may arise in long-term hydrate production simulations.

In this study, we investigate the settlement characteristics of the NGH reservoir with heterogeneous hydrate distribution during gas production via horizontal well depressurization. A fully coupled thermal-hydraulic-mechanical-chemical (THMC) numerical simulation approach is employed to analyze the process. Additionally, the effects of gas production pressure and well location on reservoir settlement behavior are systematically examined.

2 Numerical Simulation Method

In this study, the hydrate reservoir is constructed based on the two-phase flow in porous media. We enhance the chemical processes and engineering thermophysical processes to describe dissociation of the hydrate reservoir. The reliability of the model has been verified in our previous paper [25].

2.1 Continuity Equations for Two Phases Flow

The Continuity Equations of gas and water flow in porous media can be written as: ρwCp,w∂Pw∂t+∇ρwvw=m⋅w+ρwCp,w∂Pg∂t(1) ρgS∂Pg∂t+∇ρgvg=m⋅g+ρgCp,w∂Pw∂t(2) vi=−Khkr,iμi∇Pi+ρig∇D,i=g,w(3) where, i = g , w means non-wetting and wetting, respectively, ρ i is density, v i is Darcy’s velocity, μ i is dynamic viscosity, g is the gravitational acceleration, P i is fluid pressure, C p , w is a function related to pore pressure and saturation and can be defined as: Cp,w=−nwg∂Sw∂Pc(4) nwg=n1−Sh−nr(5) Pc=Pg−Pw=P0 Sw−1/m−11−m(6) where, n and n w g are initial porosity and effective porosity, respectively, nr is residual volume fraction, nr = 0.248 [37]; P c is pore pressure, P g and P w are the pressure of the non-wetting phase and wetting phase, respectively, S w and S h are wetting phase saturation and hydrate saturation.

In Eq. (7), S is the storage coefficient which can be defined as: S=nwgSgχf+Cp,w(7) χf=1ρg∂ρg∂Pg(8) where, χ f is the gas compression coefficient. According to the gas state equation, the relationship between the density of methane gas, temperature, and pressure can be expressed as: ρg=PgMgRT1+9128PgPg,cTcT1−6Tc2T2−1(9) ∂ρg∂Pg=MgRT1+9128PgPg,cTcT1−6Tc2T2−2(10) where, P g , c and T c are the critical pressure and critical temperature, and the corresponding values are 4640 kPa and 190.8 K, respectively.

The coefficient K h and k r , i in Eq. (11) are absolute permeability and relative permeability, respectively. Based on the Amyx empirical formulation [38,39], the expression of the absolute permeability is as follows:

The relative permeability of fluid and gas was calculated by Formulas (12) and (13). krw=SwL1−1−Sw1/mm2(12) krg=1−SwL1−Sw1/m2m(13) where m and L are coefficient taken as 0.5 and 0.7, respectively [38,39].

2.2 Kinetic Model of Hydrate Decomposition

NGH is composed of water molecules and methane molecules. The NGH decomposition and synthesis reaction is as follows:

Based on the Kim-Bishnoi empirical formulation [40], the production rate of methane gas is calculated by Formula (16): m˙=kdMgAsPe−Pg(15) where, M g is the molar mass of methane gas, P e is the three-phase equilibrium pressure of NGH, P e = 1.15 exp ( 49.3185 − 9459 T ) , k d is the decomposition constant of NGH, which is described by the following empirical formula: kd=k0e−ΔEaRT(16) where, k 0 is taken as 3.6 × 104 mol/(m2·Pa·s) [41], R is the ideal gas constant, T is temperature, Δ E a is the kinetic energy, J, Δ E a / R = 9753.73 K [42], A s in Eq. (18) is the reaction area of NGH particles [43], which is described by the following empirical formula: As=nShAgeo(17) where, A g e o is the surface area ratio of hydrate particles per unit volume, assuming that the natural gas hydrate particles are spherical particles, the calculation is 7.5 × 105 m−1 [44].

2.3 Energy Conservation Equation

The decomposition of natural gas hydrate sediments is a process of absorbing heat from the surrounding environment, and the heat exchange process can be controlled by heat conduction and heat convection equations. ρCeq∂T∂t+ρgCgvg+ρwCwvw·∇T−∇·λeq∇T=Q(18) ρCeq=nρgSgCg+nρwSwCw+nρhShCh+1−nρsCs(19) λeq=nSgλg+nSwλw+nShλh+1−nρsλs(20) Q=−m˙hΔHD−nρSgσg∂Pg∂t−ρgvgσg∇Pg(21) ΔHD=c1+c2T(22) where, T represents the temperature, which is the dependent variable of the energy governing equation. Cg, Cw, Ch, Cs are the specific heat capacities of methane non-wetting phase, wetting phase, NGH, and skeleton particles, respectively. λ g , λ w , λ h , λ s are thermal conductivity of methane non-wetting phase, wetting phase, NGH, and skeleton particles, respectively. Δ H D is the enthalpy of the methane decomposition process. v g , v w are Darcy velocities of gas and water. σ g is the Joule–Thomson throttling coefficient [44,45] (Appendix A).

2.4 Mechanical Models

We have assumed that the solid particles cannot be compressed, and the settlement of formation only causes the compaction of pore volume. Then the real porosity of the formation can be expressed as: np=1−(1−n0)/J(23) J=det(F)(24) F=I+∂u∂x(25) where n p and n 0 are the real and initial porosities, J is the determinant of the deformation gradient F, I is identity tensor. u represents the dependent variable of the mechanical equilibrium equation, and the main governing equations are described as follows [46,47,48]: ρc∂2u∂t2=∇⋅S+F⋅v(26) where, ρ c is the equivalent density of hydrate-bearing sediment, and its calculation formula is ρ c = ( 1 − n ) ρ s + n S h ρ h , n is the porosity. ρ s , ρ h are density of skeleton particles and hydrate, respectively. σ is the stress tensor which can be expressed as: σ=σ0+C:ε(27) where, S 0 is the initial stress tensor. C is the elastic coefficient matrix. ε is the total strain tensor, which can be expressed by displacement gradient. The calculation formula of displacement gradient is as follows:

According to the incremental theory:

The ε i j e and ε i j p are elastic and the plastic strain, respectively.

The plastic strain increment is calculated using the plastic potential energy function:

The potential Q is written in terms of at most three invariants of Cauchy’s stress tensor. Q = Q ( I 1 ( σ ′ ) , J 2 ( σ ′ ) , J 3 ( σ ′ ) ) , I 1 σ ′ , J 2 σ ′ , J 3 σ ′ are the invariants of stress tensor. λ is a scalar that is determined from the consistency condition to keep the stress point on the yield surface, f σ ′ = 0 .

In this paper, Mohr-Coulomb yield criterion is used to judge the yield conditions: fσ1′,σ3′,Sh=σ1′−σ3′Nφ−2c′ShNφ(31) where, σ 1 ′ and σ 3 ′ are the maximum and the minimum principal stress, respectively. N φ = 1 + sin φ ′ / ( 1 − s i n φ ′ ) ; φ ′ is the internal friction angle; c ′ is the cohesion.

Mechanical Properties of NGH Bearing Sediment

The relationship between the cohesion (c), internal friction angle (φ), secant modulus (Ei) with hydrate saturation (Sh) have been tested with triaxial compression experiments and reported in our previous paper [25]. The fitting mathematical relationships are shown as follows:

3 Model Construction

Fig. 1 presents a schematic diagram of the hydrate reservoir in the Shenhu area of the South China Sea. The geological system simulated in this study is situated at site SH02, located on the northern slope of Shenhu area. NGH deposits were confirmed in this region during the gas hydrate drilling expedition (GMGS6) conducted by the Guangzhou Marine Geological Survey [49]. Logging results indicate that the water depth ranges from 800 m to1500 m, with NGH deposits found at depths between 207 and 253.4 m below the seafloor (mbsf). The saturation levels of gas hydrate and free gas vary with depth. Within the interval of 207.8–253.4 mbsf, no free gas is observed, while the saturation of gas hydrate ranges from 0%–54.5%, with an average saturation of 31%. Free gas was detected at depths range between 253.4 and 278 mbsf, with an average saturation of 13%.

Figure 1: Schematic diagram of hydrate reservoir in Shenhu area, South China Sea [50].

3.1 Domain Discretization and Boundary Condition

Fig. 2 shows mapped mesh for hydrate bearing sediment, boundary condition and production well location schematic built in COMSOL. The geometric model is simplified to a two-dimensional model without considering the rate of pressure decay along the horizontal wellbore during production. The model is 73 m in the Y direction, ±25 m in the X direction, and 1 m in thickness, respectively. The diameter of horizontal wellbore is 0.16 m. The total number of grids are discretized into 32,488 elements. Considering phase change, mass transportation and other complex condition around the wellbore, we chose the triangular grid structure while increasing the mesh grid density near the borehole.

Figure 2: Mapped mesh of the NGH reservoir, boundary condition and production well location setting in the simulation model.

Since the wellhead is far from the model boundary, a constant pressure and temperature boundary with initial condition is used. And a constant production pressure is set to simulate the depressurization production process. To prevent reformation of NGH and ice formation, the temperature of the production well is set as constant value equal to initial temperature. Loading pressure of 12.07 MPa is applied to the upper boundary. TOP1 point located at the top boundary of model is set to monitor the reservoir settlement. The lower boundary is subject to a fixed constraint, while left and right boundaries are subject to roller constraints. Other parameters employed in our model are listed in Table 1.

3.2 Initial Conditions and Simulation Parameters

Prior to production of NGH, it is essential to establish the equilibrium pressure-temperature initial conditions to ensure that the simulation model’s temperature, pressure and saturation align with the actual geological conditions. Well logging results from the SH02 site indicate that the NGH deposits are located at a depth of 208 m below sea floor (mbsf) in a water column of 800 m. Due to the absence of in-situ temperature measurements, the following empirical formula is used to calculate the initial formation temperature. T=T0+grad×Z×10−3(38) T=T0+grad×Z×10−3(39) where, T 0 is the sea floor temperature of 275.15 K, grad is geothermal gradient of 44 K/km, Z is depth of formation in meters ranging from 208 m to 281 m. Then the pore water pressure can be calculated using the hydrate pressure-temperature equilibrium curve. Initial temperature and pressure of the NGH reservoir are shown in Fig. 3. Other initial parameters used in this study are listed in Table 2. And the distribution of NGH saturation is set based on the log data obtained from GMGS6-SH02, as shown in Fig. 3c.

Figure 3: Initial pressure (a), temperature (b), NGH saturation (c), and logging data of SH02 [49].

4 Simulation Results and Discussion

4.1 Effect of Outlet Pressure on Gas Production Process

Fig. 4 shows the pore pressure and temperature change with production time. A significant pressure drop is observed nearby the production well in the early stage of production, and then it remains in a relatively stable state. The pressure drop area shows as a funnel shape, and its decomposition range is closely related to the production pressure; specially, a lower production pressure results in a larger decomposition area. In general, due to the low permeability of the reservoir, the decomposition area from the pressure drop remains limited even after six months of production. It is well known that the production pressure is the main factor to control the hydrate decomposition rate in depressurization method. After half a year of production, the radius of decomposition region is about 2 m with production pressure of 6.5 MPa and it is 4 m with production pressure of 4.5 MPa. Since hydrate decomposition and free gas production is an endothermic process, a decrease in production pressure accelerates both hydrate decomposition and gas production, leading to a more pronounced temperature drop near the wellbore.

Figure 4: Formation pore pressure and temperature change with production time.

Fig. 5 shows the gas and water production rate for per unit well length with different production pressure. Since there is free gas in the NGH reservoir, and no drainage depressurization process in the numerical simulation, water and gas can be produced instantly at the beginning of production. It shows a high initial gas production rate, then rapid decayed for the consumption of free gas. During the initial production phase, the migration of free gas and water near the wellbore towards the wellbore leads to an enhanced production rate; however, subsequent migration of free gas and water from low permeability formations is hindered, resulting in a rapid decline in production rate [51]. As the production progresses, the system pressure gradually decreases due to the discharge of free gas and water. Once it reaches the equilibrium pressure of natural gas hydrates (NGH), NGH decomposition occurs to supplement the gas source, resulting in enhancement in production rate. It has been reported that the average gas production rate of the second hydrate production test in the Shenhu area of the South China Sea is 2.9 × 104 m3/day [49]. Under similar conditions with a horizontal well length, simulation results indicate a production rate of 3.8 × 104 m3/day at an outlet pressure of 6.5 MPa, which aligns closely with field data. As the production pressure is reduced to 4.5 MPa, the gas and water production rates increase to 2 times and the produced cumulative gas increase to 3 times. Therefore, reducing production pressure during the production can improve gas production effectively.

Figure 5: Gas and water production rate (a), accumulation of produced gas and water (b) 1 m well length with different production pressure.

Fig. 6 shows NGH saturation change with production time. Since the NGH saturation reservoir is stratified, the degree of hydrate saturation within the pore increases inversely with respect to the permeability, the permeability is diminished as a result of the presence of NGH in the pore. The production pressure drop can be transmitted quickly into depth of the low Sh area, which shows a significant ahead of NGH decomposition front in low Sh layers than that in high Sh region. Since production pressure is the main factor to control NGH decomposition rate in the production process, the radius of decomposition region is about 2 m with production pressure of 6.5 MPa after half a year of production, and it is 6 m with production pressure of 4.5 MPa. Additionally, there are a few secondary hydrates occurring near the hydrate dissociation front. The heat absorption associated with NGH dissociation alters the pressure-temperature conditions of phase equilibrium, leading to the reformation of NGH.

Figure 6: NGH saturation variation with production time and different production pressure.

Fig. 7 shows the reservoir settlement characters with different production pressure. The primary settlement is concentrated in the upper deposit around the wellbore, with wellbore settlement values of 0.032 m and 0.018 m observed at pressures of 4.5 MPa and 6.5 MPa, respectively. Meanwhile, the settlement of TOP1 point in Run1 reaches to nearly 0.1 m with production pressure of 4.5 MPa. And the settlement curves show that settlement displacement will continue to increase with the NGH production. The top boundary exhibits a U-shaped sag over a production period of 180 days. The settlement rate of formation refers to subsidence in per unit time. Its maximum value occurs at the beginning of production and decreasing with production time, which is very similar to the gas production rate curve. Lower production pressures induce greater settlement and sag in the top boundary, increasing likelihood of effective stress concentration and geological landslide. This is corroborate by the settlement rate observed at the upper boundary. With the free gas and water discharged, the mechanical weakness of reservoir caused by NGH decomposition become the main reason for reservoir settlement. And the settlement rate shows a gradual decline with the production time.

Figure 7: Reservoir settlement characteristics with different production pressure. (a) Settlement, (b) Sedimentation rate (c) Top boundary settlement.

4.2 Influence of Horizontal Well Location

Fig. 8 shows NGH saturation distribution contours with different well location.

Figure 8: NGH saturation distribution contours with different well location.

The NGH saturation exerts a significant influence on the permeability, induces regional distribution of hydrate decomposition. The decomposition shape of the Run6 model starts with the wellbore and diffuses along the low saturation hydrate layer. The decomposition area forms a diamond-shaped distributed after approximately 180 days, and the movement of the decomposition interface in the lower layer of the well site is slower due to higher NGH saturation of formation. The decomposition interface of the Run3 exhibits an approximate circular shape with a radius ranging from 4 m after a period of 180 days. The occurrence of uniformly distributed, highly saturated hydrate layers surrounding the wellbore is responsible for this phenomenon. The well location of the Run7 model is located within the free gas layer, thereby facilitating the most rapid decomposition of hydrates and resulting in a largest decomposition area.

Fig. 9 illustrates gas and water production rates and accumulated production for different well location with a well length of 1 m. Initially, production rates are high at the onset of production, rapidly decreasing to a lower value, then gradually rising and stabilizing. During the decline stage of production rate, Run3 exhibits the lowest rate while Run7 is the highest. At the production enhancement stage, Run3’s rate gradually aligns with that of Run6, whereas Run7’s rate remains relatively constant. The wellbore of Run3, located in a high-saturation hydrate formation, impedes free gas and water migration due to elevated saturation. However, the wellbores of Run6 and Run7, located in lower-saturation hydrate and free gas zone, respectively, indicates the greater permeability, facilitating fluid migration to the wellbores. In the mid-to-late production stages, gas from hydrate dissociation sustains decomposition rate, resulting in a higher recovery rate for Run3 compared to Run6 in later stages. The cumulative gas production of Run7 is highest, as it includes free gas, hydrate-derived gas and dissolved gas from pore liquid.

Figure 9: Gas and water production rate (a) and accumulated production (b) with different well location for 1 m well length.

Fig. 10 shows reservoir settlement characteristics under different well locations. It is shown that after 30 days production, the settlement and settlement rate of TOP1 point in the Run3 surpasses that of the Run6. The production well in Run6 is located in the low Sh layer where with relatively high permeability, and the pore pressure drops significantly quicker than Run3, then leading to a faster settlement in Run6. With the continuous gas production, the hydrate decomposition within the pressure drop influence range becomes the main gas source. The pressure drop influence range in Run6 is significantly larger than that in Run3, and the pressure drop gradient at local points is smaller than that in Run3, resulting in a smaller amount of hydrate decomposition at local points in Run6 and a smaller local settlement for formation weakening caused by hydrate decomposition. The settlement of Run7 is the largest among all cases after 180 days gas production. Since the pore pressure drop gradient is more evenly distributed due to the wellbore located in the free gas region in Run7, it shows a minimum sag in its top boundary, which is benefit to prevent the formation of tension fractures caused by uneven settlement. The settlement of Run6 is the smallest among all, while exhibiting the largest sag in its top boundary. And the settlement of Run3 is comparatively lower than that of Run6 due to its location in a highly saturated formation with enhanced strength. Above results show that well location would significantly affect the settlement and gas characteristics of the NGH reservoir. Although locating the production well in the overlying free gas layer would lead to a largest settlement, the formation uneven settlement is minimal in this way, which can effectively reduce the tensile fractures development caused by the uneven settlement. Meanwhile, the gas production rate of this wellbore arrangement can significantly increase the gas production rate.

Figure 10: Reservoir settlement characteristics under different well location (a) Displacement of point TOP1, (b) Sedimentation rate of point TOP1, (c) Displacement of upper boundary.

5 Conclusion

This study investigates the settlement behavior of horizontal wellbores in natural gas hydrate (NGH) reservoirs with heterogeneous hydrate distributions. The influences of well placement under various production pressures and single-well production conditions on both gas production and formation settlement are analyzed. Results shows that although lower production pressures enhance productivity, but also lead to increased formation settlement and a higher risk of tensile fractures. Previous studies have also indicated that excessively low production pressures may promote secondary hydrate formation, thereby impeding gas flow [52]. Thus, it is crucial to optimize production pressure based on specific NGH reservoir characteristics. Adjusting the horizontal well position can help compensate for the production deficiency caused by high outlet pressure. The optimal configuration appears to place the well near free gas zones within low-saturation hydrate formations. This strategy not only improves gas production but also mitigates uneven settlement, thereby reducing the potential for tensile fracture development.