Coupling Submarine Slope Stability and Wellbore Stability Analysis with Natural Gas Hydrate Drilling and Production in Submarine Slope Strata in the South China Sea

Abstract

This research explores the geomechanical challenges associated with gas hydrate extraction in submarine slope zones, a setting posing a high risk of significant geological calamities. We investigate slope and wellbore deformations driven by hydrate decomposition within a subsea environment. Utilizing Abaqus, a fluid-solid-thermal multi-field coupling model for gas hydrate reservoirs was created. Hydrate decomposition during drilling is minimal, resulting in minor formation deformation near the wellbore. However, a year of hydrate production caused a maximum displacement of up to 7 m in the wellbore and formation, highlighting the risk of submarine landslides. This indicates the need for meticulous surveillance of formation subsidence and wellhead equipment displacement. In the aftermath of a hydrate-induced submarine landslide, both the hydrate layer and the overlying strata descend together, inflicting considerable damage on the formation and wellbore. Our study presents a holistic examination of the interplay between environmental geomechanics risks and engineering structure risks for submarine slope instability and wellbore stability during hydrate development, providing crucial insights for enhancing safety measures in hydrate drilling and production, and ensuring wellbore stability.

1. Introduction

Natural gas hydrates (NGHs) resemble ice-like crystalline solids formed from natural gas and water under conditions of high pressure and low temperature, prevalently distributed in permafrost zones and shallow sediments in the deep sea [1,2]. NGH’s storage capacity is substantial, and its carbon content is nearly double the total carbon content in conventional fossil fuels, making it a promising alternative energy source [3]. However, the harsh environments of NGHs accumulation, their unique occurrence conditions, and their complex nature pose considerable challenges for current exploitation technologies, which have yet to meet the prerequisites for large-scale commercial hydrate development [4].

Apart from the issues of slow decomposition and low efficiency confronted by hydrate development, the potential geological and engineering problems arising from large-scale hydrate development also significantly impede its progression [5]. Existing studies suggest that natural hydrate decomposition could trigger extensive seabed landslides, such as the second phase of the Storegga seabed landslide, which was closely tied to hydrate decomposition [6]. Man-made hydrate extraction activities can similarly induce geological and engineering problems, as evidenced by the notable subsidence of the seabed surface in the hydrate decomposition area during a short-term hydrate trial in Japan’s Nankai sea area [7]. Secondary issues, such as engineering safety concerns caused by seabed landslides and subsidence, and greenhouse gas leaks, further complicate hydrate development [8,9].

Simultaneously, due to the control exerted by hydrate phase equilibrium pressure, the proven water depth range for marine hydrates generally falls between 800 and 2000 m [10,11]. This depth corresponds predominantly to continental shelves. Although the overall inclination of continental shelves is minimal, numerous secondary slopes occupy extensive areas. Unlike terrestrial landslides, submarine landslides can occur at much lower angles, in some cases, even less than 2° [12,13]. Current research suggests that the second phase of these landslides is triggered by the natural decomposition of hydrates. Therefore, studying the impact of large-scale hydrate decomposition on the stability of submarine slopes is of critical importance.

The peculiar mechanical properties of offshore hydrate reservoirs explain why offshore hydrate development can result in seabed landslides and subsidence. Unlike conventional oil and gas reservoirs, hydrate reservoirs are semi−consolidated or unconsolidated, weakly cemented, and rich in mud content. Hydrates, prior to decomposition, act as bonding and supporting agents within the reservoir [14]. At saturations generally exceeding 40%, pore-filling hydrates transition to a load-bearing state, where the hydrates bridge neighboring grains and become integral components of the sediment skeleton, thereby contributing to the mechanical stability of the sediment [15]. Post decomposition, the reservoir framework collapses, and the water content in the reservoir dramatically increases, leading to seabed instability and even triggering geological catastrophes like seabed landslides and subsidence [16].

Current strategies for hydrate trial production include pressure reduction, heating injection, chemical injection, and CO₂ exchange methods [17,18,19,20]. Regardless of the technique employed, the primary process for hydrate production still involves the drilling and development methods utilized in the oil and gas industry. Given that the wellbore is the sole conduit linking the hydrate reservoir to the exterior, maintaining its integrity during hydrate development is pivotal for efficient hydrate production [21]. Ensuring wellbore stability during hydrate development has significant implications for the feasibility evaluation of regional hydrate development, disaster prevention, and well structure design [22,23,24].

Nevertheless, there is a current lack of research into the effects of hydrate decomposition on wellbore stability, along with relevant landslides and seabed subsidence. The issue of hydrate−induced underwater landslides involves complex multiphase fluid mechanics, soil mechanics, hydrate decomposition kinetics, and heat and mass transfer processes. Owing to the vast scale of underwater landslides and the complexities of the hydrate simulation development system, along with the poor adaptability of the hydrate drilling simulation system, studying this issue via indoor experiments proves challenging. Comparatively, numerical simulation can more effectively represent the coupling of multiple fields in the hydrate development process, thus replicating the full process of drilling, recovering, and landslide in hydrate−rich underwater slope regions [11,25]. Numerical simulation research on hydrates has also evolved from early fluid–solid coupling analysis to heat–flow–chemical–solid multifield coupling analysis. Currently, scholars use finite element software such as Comsol, Hydrate+Tough, Abaqus, and Flac3D to conduct mature multiphase coupling numerical simulation of hydrate formations [26,27,28,29,30]. Therefore, this paper uses numerical simulation to evaluate the stability of the wellbore during the drilling of hydrates.

2. Methodology

This part introduces the core tool and method of this research, which is a numerical model for describing the thermal-fluid-solid-phase change coupling of the hydrate layer, developed based on Abaqus secondary development and subroutine establishment. This model is combined with the finite element strength reduction method. Based on this method, the feasibility of the method is proved in this chapter by verifying existing hydrate decomposition experimental results and numerical simulation results and analyzing benchmark cases.

2.1. Description of Mechanical Properties of Natural Gas Hydrate Reservoirs and Wellbore

The properties of hydrate reservoirs in low-temperature permafrost regions differ from those of hydrate reservoirs under high-pressure conditions on the seafloor. The hydrate reservoir properties discussed in this paper pertain exclusively to those found in seafloor conditions [15]. Hydrate decomposition kinetics describe the rate of decomposition and production of decomposition products of hydrates under changing external conditions. It is a key factor in determining the saturation field of hydrate formations.

This paper implements the dynamic coupling of Abaqus porous media and the hydrate saturation field by compiling related hydrate decomposition kinetics equations into subroutines using the USDFLD subroutine. Abaqus has advantages over other general-purpose finite element software in the fields of solid mechanics and geotechnics because it offers a wide array of constitutive models, more diverse contact mechanisms, and allows users to customize materials through subroutines. This ensures a comprehensive understanding of the constitutive models and mechanical behavior of hydrate reservoirs. In terms of multi−physics coupling, Abaqus can focus computational problems on the mechanical properties of the soil and wellbore while ensuring multi-field coupling. This enables us to better grasp the overall stability of the slope after hydrate dissociation, while controlling the hydrate decomposition process through multiple physical fields. Based on hydrate saturation, the dynamic changes of hydrate reservoir mechanical strength and seepage properties are further realized, and the dynamic update and coupling solution of natural gas hydrate reservoir properties in Abaqus is finally completed.

2.2. Hydrate Decomposition Kinetics Model

Natural gas hydrates are formed by alkanes, primarily methane, and water under low-temperature and high-pressure conditions. Currently, methane hydrates are the most widely distributed and studied, with methane accounting for more than 95% of the gas products, and its decomposition equation is (Moridis, 2003):

$$\mathrm{C}{\mathrm{H}}_4\cdot n{H}_2\mathrm{O}+\mathrm{heat}\leftrightarrow \mathrm{C}{\mathrm{H}}_4+n{H}_{2}O$$

(1)

where, n has a value range of 5.5–6, and in this work, it is taken as 6 [31]. The hydration number n primarily determines the ratio of water molecules to gas in the hydrate. Upon hydrate dissociation, both water and gas are expelled through the formation pores and wellbore, without inhibiting the rate of hydrate decomposition. Consequently, this does not have an impact on the mechanical strength of the formation.

The decomposition of natural gas hydrates is an endothermic reaction [32], with a reaction heat of 53 kJ·mol⁻¹ to 58 kJ·mol⁻¹.

The phase equilibrium conditions for methane hydrate are shown in the following equation [33,34]:

$$P_{\mathrm{eq}}=\mathit{exp}\left(e_1-\frac{e_2}{T}\right)$$

(2)

where $P_{eq}$ is the phase equilibrium pressure at temperature T, Pa, $T$ is temperature, °C, and $e_1$ and $e_2$ are the decomposition coefficients, with $e_1$ = 39.08 and $e_2$ = 3533 when $T$ > 0 °C.

The decomposition rate of hydrates is influenced by various factors. Based on the research results of Sun et al., the hydrate decomposition kinetics can be simplified as a function of pressure, temperature, permeability, effective porosity, and time [34,35]. The decomposition rate, gas production rate, water production rate, and reaction heat of hydrates can be described as:

$$\begin{gathered} -g^h=M_h{k}_{\mathrm{reac}}{\Gamma }_r{A}_s\left(P_{\mathrm{eq}}-f_g\right) \\ {g}^g=g^h\frac{M_g}{M_h} \\ {g}^w=n\ast g^h\frac{M_w}{M_h} \\ {A}_s=\sqrt{\frac{{\varphi }_e^3}{2k}}{\Gamma }_r={\left(s_g{s}_w{s}_h\right)}^{2/3}\hspace{1em} \\ {h}_a=h_r\frac{g^h}{M_h} \\ k\left(s_h\right)=k_0{(1-s_h)}^N \end{gathered}$$

(3)

• where $g^h$ is the rate of decomposition of the natural gas hydrate, kg·m⁻³·s⁻¹, $g^g$ is the rate of generation of natural gas, kg·m⁻³·s⁻¹, $g^w$ is the rate of water generation, kg·m⁻³·s⁻¹, $M_g$ is the relative molecular mass of methane gas, $M_h$ is the relative molecular mass of natural gas hydrate, $M_w$ is the relative molecular mass of water, $k_{\mathrm{reac}}$ is the kinetic constant of hydrate dissociation, mol·m⁻²·s⁻¹·Pa⁻¹, $f_g$ is the fugacity of methane gas obtained from the Peng–Robinson equation, Pa, $A_s$ is the total surface area between the hydrate particles and the surrounding phase, m²·m⁻³, ${\Gamma }_r$ is te ratio of the active reaction surface to the total surface area, $k$ is the permeability, $k_0$ is the intrinsic permeability, and $N$ is the permeability component;
• where $A_s$ and ${\Gamma }_r$ are determined based on the relationship proposed by Sun and Mohanty (Sun and Mohanty, 2006). Additionally, the theory suggests that the generation of the hydrate phase takes place on the reservoir framework rather than directly within the reservoir pore. The decomposition process is the same;
• where $k$ is the permeability of the hydrate layer, which is a function of the initial and current saturation of the hydrate $\mathrm{m}\mathrm{D}$, $s_g$ is the saturation of methane gas in the reservoir pore, $s_w$ is the water saturation in the reservoir, $s_h$ is the saturation of the natural gas hydrate, $h_a$ is the total heat absorbed by the reaction, 10³ kJ·m⁻³·s⁻¹, and $h_r$ is the reaction heat of the hydrate decomposition per unit amount of substance, 53 kJ·mol⁻¹.

2.3. Hydrate Reservoir Mechanics Model

The pore-scale habit of gas hydrates is of critical importance and warrants thorough examination. Various morphological forms of gas hydrates can occur within sedimentary structures, including pore-filling, grain-coating, and vein-like structures [36,37]. Each of these forms carries unique implications for wellbore stability, relevant to both drilling and production phases. This paper focuses on pore-filling hydrates, which occupy the interstitial spaces in the sediment and serve as a cementing agent, thereby bolstering the overall mechanical strength of the formation.

The experiment shows that the strength of hydrate deposits increases linearly with the increase in hydrate saturation and confining pressure. Under fixed confining pressure, the bulk modulus, shear modulus, and cohesion are represented as [11]:

$$\begin{gathered} K=\left(K_{\mathrm{sh}1}-K_{\mathrm{sh}0}\right)\times s_h+K_{\mathrm{sh}0} \\ G=\left(G_{\mathrm{sh}1}-G_{\mathrm{sh}0}\right)\times s_h+G_{\mathrm{sh}0} \\ c=\left(c_{\mathrm{sh}1}-c_{\mathrm{sh}0}\right)\times s_h+c_{\mathrm{sh}0} \end{gathered}$$

(4)

where: $K$ is the layer volume modulus of the natural gas hydrates, Pa, $K_{\mathrm{sh}1}$ is the layer volume modulus of the natural gas hydrates when saturation is 1, Pa, $K_{\mathrm{sh}0}$ is the layer volume modulus of the natural gas hydrates when saturation is 0, Pa, $G$ is the layer shear modulus of the natural gas hydrates, Pa, $G_{\mathrm{sh}1}$ is the layer shear modulus of the natural gas hydrates when saturation is 1, Pa, $G_{\mathrm{sh}0}$ is the layer shear modulus of the natural gas hydrates when saturation is 0, Pa, $c$ is the internal cohesion of the natural gas hydrates layer, Pa, $c_{\mathrm{sh}1}$ is the internal cohesion of the natural gas hydrates layer when saturation is 1, Pa, and $c_{\mathrm{sh}0}$ is the internal cohesion of the natural gas hydrates layer when saturation is 0, Pa.

2.4. Analysis of Porous Media in Abaqus

The effective stress principal is adopted in Abaqus to express the basic stress determining rules for porous media, shown as:

$$\begin{array}{c}{\overline{\mathit{\sigma }}}^{\mathrm{*}}=(1-n_t)\overline{\mathit{\sigma }}-n_t{\overline{p}}_t\mathbf{I},\end{array}$$

(5)

where ${\overline{\mathit{\sigma }}}^{\mathrm{*}}$ is the effective stress, $n_t$ is the proportion of wet, $\overline{\mathit{\sigma }}$ denotes the skeleton effective stress, and ${\overline{p}}_t$ represents the average compressive stress of wet liquid.

The solution strategy based on the effective stress principle involves the following steps:

(1) Discretizing the porous media to establish the stress equilibrium statement by using the virtual work principal:

$${\int }_V\mathit{\sigma }:\delta \mathit{\epsilon }dV={\int }_S\mathit{t}\cdot \delta \mathit{v}dS+{\int }_V{\mathit{f}}_1\cdot \delta \mathit{v}dV+{\int }_V(sn+n_t){\rho }_w\mathit{g}\cdot \delta \mathit{v}dV$$

(6)

where $\mathit{\sigma }$ and ${\mathit{f}}_1$ represent the total stress matrix and the body-force vector, respectively, $\mathit{t}$ is the surface traction vector in unit area, $\delta \mathit{v}$ and $\delta \mathit{\epsilon }$ represent the virtual velocity field and the matrix of virtual strain rate, respectively, $s$ is the liquid saturation, $n$ is the liquid volume ratio to pore volume, and ${\rho }_w$ is the density of liquid;

(2) Determining the constitutive behavior of both the fluid and skeleton structure within the porous medium [8];

(3) Ensuring the continuity of the wetting liquid phase throughout the porous medium, as described [19]:

$${\int }_V{\rho }_w(d{V}_w+d{V}_t)={\int }_V{\rho }_w(n_w+n_t)dV$$

(7)

where $V_w$ and $V_t$ separately denote the volumes of unbound liquid and bound liquid and $n_w$ is the proportion of unbound liquid in a unit volume.

2.5. Flowchart of the Methodology

Figure 1 presents a flowchart outlining the methodology. The yellow section illustrates how hydrate saturation is updated within the computational domain, a crucial aspect of the simulation. The green box denotes the primary theoretical framework employed in this method. The blue box highlights the macroscopic objects, loads, and ultimate research objectives.

2.6. Model Verification

Gupta et al. conducted experiments on hydrate mechanical properties and conducted numerical simulations combining heat, chemical, water, and geomechanical fields [22]. The mechanical experiments were carried out using the NESSI high-pressure equipment (Seabed Interaction Natural Environment Simulator) [26]. The stainless-steel container has a volume of 40 L and is equipped with a three−axis system. Hydrate sediment was first generated in situ in the container and then decomposed by pressure reduction method. Based on the hydrate core experiment conditions, Gupta et al. established a multi-field coupling numerical model based on C++. The experimental conditions and boundary conditions for numerical simulations are shown in Figure 2, and some of the parameters used in the numerical model are listed in

Table 1. Parameters used in the experiment.

Parameter	Value	Unit	Parameter	Value	Unit
$k_s^c$	$0.383\cdot \mathit{ln}(T)-1.581$	W·m−1·K−1	${\rho }_s$	2100	kg·m−3
$k_h^c$	2.1	W·m−1·K−1	$k_{\mathrm{reoc}}$	3.2 × 10−10	mol·m−2
$k_s^c$	1.9	W·m−1·K−1	$n$	5.75	Pa−1·s−1
${\rho }_h$	$\frac{{\rho }_g}{z{R}_{g}T}$	kg·m−3	$E_s$	160	$\mathrm{MPa}$
${\rho }_w$	$\mathrm{Vapor}:0.0022\frac{{\rho }_g}{T}$ Liquid: 1000	kg·m−3	$E_h$	360	$\mathrm{MPa}$
$c$	3	kPa

.

In assessing the reliability of our model, we reconstructed a multiphase coupling model of a hydrate reservoir based on Abaqus, in line with the theoretical premises and the experimental and numerical conditions proposed by Gupta et al. Figure 3 presents the experimental and numerical results from Gupta et al., detailing data on gas production and volume strain resulting from hydrate decomposition. In their experiment, hydrate decompression via pressure reduction was initiated on the 12th day. Figure 3a depicts the temporal evolution of gas production from hydrates, exhibiting a generally linear trend over time. As illustrated in Figure 3a, aside from the final experimental outcomes, the cumulative gas production from both numerical simulations closely aligns with the experimental data and provides a good fit. Figure 3b illustrates the total volume strain of the samples, serving as an indicator of the reliability of both the hydrate sediment geomechanics model and hydrate kinetics. In Gupta’s work, the curves for both the experimental and numerical methods display a step-like upward trend for the initial 0.2 days before the start of hydrate decomposition. Conversely, during this period, the numerical simulation results generated in our study indicate a steady increase in volume strain. Moreover, during the subsequent hydrate decomposition process, the rate of increase in volume strain significantly diminishes and essentially stabilizes.

3. Results

3.1. Overview of the Model

The overall gradient of the northern continental slope in the South China Sea is approximately 2 degrees, though there exist secondary slopes on the continental slope that feature steeper angles, of around 10 degrees. The expanse of these seafloor slopes is generally substantial, with some spanning several tens of kilometers, as shown in Figure 4 [38]. However, a single secondary slope typically spans around 1 km. Given the narrower span of these secondary slopes, their incline is steeper, and they pose a higher risk of landslides when developing gas hydrates. Therefore, the subject of this study is a single secondary slope.

During the experimental development process of gas hydrates, vertical wells are commonly employed. Owing to the nearly perpendicular angle between the vertical wellbore and the seafloor landslide, the vertical wellbore experiences a stronger shear force. Conversely, while the use of horizontal wells can significantly boost the speed of gas hydrate development, the spatial relationship between the structure and location of the wellbore in horizontal wells and the slopes is challenging to predict. Therefore, a vertical well is adopted in this study.

The completion engineering design of NGH involves planning and researching wellbore types, completion methods, sand control measures, and cementing techniques. The actual wellbore structure upon completion of gas hydrate drilling is intricate. However, from a finite element perspective, an excessively detailed wellbore structure leads to an overly complicated model that hampers calculation convergence.

In response to this challenge, the present study greatly simplifies the relevant structures. As depicted in Figure 5, the model assumes a vertical wellbore positioned at the center of the slope, piercing through both the overlying rock layer and the gas hydrate layer. Given that the formation near the gas hydrate is the most fragile, the bit size and casing size for the hydrate layer, along with the corresponding cement ring thickness, are selected as key parameters for the wellbore in the numerical model. The wellbore structure of hydrate wells is similar to that of conventional oil and gas wells. However, due to the generally shallower burial depth of hydrates, there are fewer stratigraphic layers. From top to bottom, the wellbore comprises casing and tubing, followed by open−hole completion. From the inner to the outer layers, it includes the casing and cement sheath. The surface casing section may consist of two layers of casing and a cement layer. Therefore, the weakest part is the lower cement layer. Referencing previous structural designs of vertical and horizontal gas hydrate wells conducted by Japan and China, the wellbore diameter within the gas hydrate layer in our model is set to 311.15 mm, with a casing diameter of 244.5 mm [39,40,41]. Using this as a basis, in this model, the wellbore is a single cement ring with a uniform thickness, an outer diameter of 311 mm, an inner diameter of 245 mm, and a thickness of 33 mm. In the gas hydrate layer, the wellbore is set as a permeable interface to simulate the wellbore seepage of perforation completion. In the sedimentary layer, the wellbore is set as an impermeable interface, and the pressure inside the wellbore will not affect the formation pressure.

Based on the preceding analysis, this study uses Abaqus software to construct a series of 3D models of hydrate-bearing slopes. The total length of the numerical model is 1200 m, incorporating an 800 m span for the slope, and a horizontal seabed plane of 200 m on each side. The model’s maximum height is set at 600 m, and the minimum height at 400 m, resulting in a slope angle of 14.04°. The hydrate zone is situated directly beneath the slope’s highest point, buried at a depth of 200 m and running parallel to the slope’s surface. At the model’s highest point, the height of the overlying seawater is 1300 m. This model allows us to accurately simulate and analyze the thermal-fluid-solid-phase change interactions and the structural stability of the hydrate-bearing slopes.

3.2. Model Mesh and Boundary Conditions

Figure 5 illustrates the numerical model’s mesh scheme. The blue region symbolizes the hydrate layer, while the gray zone denotes the sedimentary layer devoid of hydrates, and the red area signifies the wellbore section. To accurately capture the coupled interactions between stress and pore pressure, the C3D8P mesh type is employed. Specifically designed for porous media in Abaqus, this 8-node trilinear mesh accommodates both displacement and pore pressure in soil, making it apt for simulating the behaviour of sand on the slope and contact behaviour with the wellbore. The mesh is refined around the wellbore and the neighboring strata, with 16 seeds per half circumference of the wellbore and three layers radially. In the x-direction, the mesh size gradually increases on both sides from the wellbore to efficiently characterize the layers near the hydrate wellbore while also reducing the overall computational load. In the y-direction, the mesh seed number is 16. In the z-direction, the number of meshes for the overlying and underlying strata is 18, while for the hydrate layer, it is 8. The model comprises a total of 83,574 mesh elements. In this instance, the stress field, seepage field, temperature field, and a custom field variable (the hydrate saturation field) can be simultaneously coupled and computed. This grid design allows for a detailed understanding of the interactions between these different factors within the model.

3.3. Parameters and Initial States

In the study, it is assumed that the temperature at the seafloor stays constant at 2 °C and the temperature gradient of the seafloor is 45 °C/km. The densities of both the seawater and the pore water are presumed to be unchanging. The model simulates the decomposition of gas hydrates brought about by variations in the wellbore’s pore pressure. Additionally, the model emulates the propagation of temperature and the transport of fluid during the depressurization and production of gas hydrates in the formation. The thermodynamic and hydraulic parameters used in the model are detailed in

Table 2. Parameter for the numerical model.

Parameters	Value	Parameters	Value
Porosity	0.4	Sea water density	1040 kg·m−3
Hydrate density	920 kg·m−3	Sea water thermal conductivity	3.0 W·m−1 K−1
Hydrate thermal conductivity	0.5 W·m−1 K−1	Sea water specific heat capacity	4.2 kJ·kg−1 K−1
Hydrate specific heat capacity	2.1 kJ·kg−1 K−1	Intrinsic water saturation	1
Formation skeleton density	2600 kg·m−3	Intrinsic water saturation with hydrate	$1-s_h$
Formation skeleton thermal conductivity	1.0 W·m−1 K−1	$\mathrm{Intrinsic}\mathrm{permeability}{s}_h=0$	0.6 mD
Formation skeleton specific heat capacity	1.0 kJ·kg−1 K−1	$\mathrm{Intrinsic}\mathrm{permeability}{s}_h=1$	0.2 mD
Seabed seawater temperature	2 °C	Geothermal gradient	0.035 °C·m−1
$\mathrm{Formation}\mathrm{cohesive}\mathrm{strength}{s}_h=0$	0.1 MPa	$\mathrm{Volumetric}\mathrm{modulus}{s}_h=0$	2.0 MPa
$\mathrm{Formation}\mathrm{cohesive}\mathrm{strength}{s}_h=1$	2.0 MPa	$\mathrm{Volumetric}\mathrm{modulus}{s}_h=1$	32 MPa
$\mathrm{Formation}\mathrm{shear}\mathrm{modulus}{s}_h=0$	1.5 MPa	Internal friction angle	5°
$\mathrm{Formation}\mathrm{shear}\mathrm{modulus}{s}_h=1$	29 MPa	Effective stress coefficient	1
Wellbore cement cohesive strength	60 MPa	Friction angle of wellbore cement	15°
Wellbore cement Young’s modulus	500 MPa	Poisson’s ratio of wellbore cement	0.3

. This simulation offers an accurate representation of the thermal and fluid dynamics involved in the process of gas hydrate extraction. It also offers a framework for understanding how different variables might affect the extraction process, providing critical insights into the optimization of gas hydrate extraction techniques [18,42,43].

Having incorporated the given parameters, boundary and initial conditions, and subroutines into the software, the initial state of the model can be established. As illustrated in Figure 6, from top to bottom, the graphic represents the initial distribution of pore pressure, temperature, and maximum principal stress, respectively.

The distribution of various fields in the formation, after initialization, aligns with the theoretical understanding and demonstrates an appropriate degree of non-uniformity, in accordance with the material properties. This accurately mirrors the real conditions of the formation. Such visual representations of the initialization process are essential in comprehending the baseline or starting conditions of the simulation, thus, offering a benchmark against which to measure any subsequent changes or variations.

4. Results

4.1. Slope and Wellbore Deformation during Hydrate Drilling Process

Deep-water hydrate reservoirs are typically located at relatively shallow depths and possess the characteristics of unconsolidated or semi-consolidated formations. This means that they do not usually encounter complex structures during the drilling operations, which significantly shortens the time required from drilling to the completion of the hydrate layer. During hydrate drilling operations, precautions are taken to prevent wellbore instability, like balanced drilling and the use of low-temperature drilling fluids. Additionally, hydrate inhibitors are strategically added to the drilling fluid. Due to these measures, the degree of hydrate decomposition during the drilling process remains very low.

The deformation of the wellbore and the formation is primarily dependent on the impact of drilling on local in situ stress, rather than large-scale changes. This highlights the importance of maintaining stable drilling conditions and using appropriate drilling fluids to mitigate the risks associated with the decomposition of hydrates, such as wellbore instability or potential changes to the formation’s mechanical properties.

As shown in Figure 7a, when drilling is completed, only a small range of formation around the wellbore undergoes a slight deformation, of less than 10⁻³ m, and the displacement value is negligible compared to the overall scale of the formation. The range of influence on the formation gradually increases from top to bottom, which is related to the increasing overlying pressure. Figure 7 shows the distribution of maximum principal strain, maximum principal stress, and deformation of the wellbore. In Abaqus, compressive stress and strain are positive values. As can be seen from Figure 7b,c, the deformation of the wellbore is mainly radial towards the wellbore, and the direction of the maximum principal stress is close to the horizontal direction. The maximum displacement of the wellbore is about 4 × 10⁻⁴ m. Meanwhile, results show that the wellbore did not undergo plastic deformation during the drilling stage. Therefore, it can be concluded that during the hydrate drilling stage, significant hydrate decomposition and large-scale deformation and damage of the formation did not occur, and the wellbore remained stable.

4.2. Formation Subsidence and Wellbore Integrity

In the process of hydrate decompression and production at a pressure of 6 MPa over a year, the saturation distribution of the hydrate can be seen in Figure 8a. Due to the thermal and pressure conditions of the formation, the hydrate dissociation zone emerges as a triangle in the plane view and as a cone in the three-dimensional perspective. During the process of hydrate dissociation involving the slope area, the degree of hydrate dissociation is not uniform on both sides of the wellbore. This discrepancy can be attributed to the initially low pore pressure in the relatively elevated position.

Figure 8b displays the distribution of pore pressure within the formation. The lower the initial pore pressure of the formation, the more it facilitates the propagation of wellbore pressure to further formations, thereby resulting in a larger hydrate dissociation zone. Due to the low permeability of the hydrate area, the propagation distance of the pressure drop within the hydrate dissociation zone is significantly lower than that within the sedimentary layer. This highlights the role of the permeability in the hydrate dissociation process, which influences the pressure propagation and, subsequently, the size of the hydrate dissociation zone.

Figure 8c illustrates the distribution of maximum principal stress in the formation after depressurization. By comparing Figure 6c and Figure 8c, it becomes evident that the decreased pore pressure increases the effective stress in the formation. Concurrently, the decomposition of hydrates undermines the strength of the formation. The combined impact of these two factors results in significant damage to the formation skeleton within the hydrate decomposition zone. There is evident deformation both in the decomposition zone and the formation affected by the pressure drop.

Figure 8d presents the distribution of principal strain in the formation post-depressurization. In the core area near the wellbore, the direction of the maximum principal strain approximates the direction of the vertical slope. Nonetheless, there is no noticeable deformation in other areas. This distribution implies that the most substantial changes and potential for damage occur near the wellbore, where the depressurization and hydrate decomposition effects are strongest.

Figure 9a illustrates the significant displacement that has occurred within the formation due to gas hydrate dissociation. The displacement is up to approximately 7 m within the gas hydrate dissociation zone, with the affected area extending across the entire formation. These numerical simulation results align closely with those from previous studies on seafloor subsidence that did not involve slopes. However, in sloping formations, the direction of deformation is not vertically downward, but, rather, roughly perpendicular to the direction of the slope. The simulation results suggest that the angle between the direction of formation displacement and the wellbore is a key determinant in causing wellbore deformation.

This means that the structural integrity of the wellbore can be significantly impacted by the inclination of the slope in the formation. In such scenarios, the wellbore might suffer from deformation or even damage, which can cause issues with drilling or production operations, and, potentially, pose risks to safety and environment.

The wellbore exhibits significant deformation, as illustrated in Figure 9. The wellbore displacement aligns closely with the direction of ground subsidence, predominantly characterized by a downward shift, complemented by substantial lateral displacement in the positive x−direction. As Figure 9b demonstrates, displacement gradually decreases from the wellbore’s top to its bottom. The top of the wellbore sees minimal lateral displacement, which then escalates with increasing depth. Around the 60 m mark, the wellbore’s maximum lateral displacement hits approximately 1 m. As the depth continues to grow, the displacement in both directions tapers off substantially, ultimately reaching 0.49 m horizontally and 2.38 m vertically at the wellbore’s base.

Figure 9c illustrates the wellbore’s equivalent plastic strain, effectively portraying the cumulative plastic strain throughout the entire deformation process. In the context of the wellbore, it indicates the extent and magnitude of the plastic deformation that takes place.

The entire wellbore undergoes plastic strain, however, for most of the wellbore, the plastic strain value remains below 0.1. At this juncture, the ground itself has not yet exhibited any signs of plastic deformation.

4.3. Submarine Landslides and Wellbore Integrity

As the extraction of gas hydrates proceeds, there is a risk of submarine landslides in the geological layers, with the displacement distribution and vectors shown in the corresponding Figure 10. Figure 10a shows the distribution of equivalent plastic strains within the geological strata. The highest plastic strain is seen in the gas hydrate dissociation area. The lower boundary of plastic deformation runs along the gas hydrate layer. At this stage, a large-scale landslide transpiring across the entire slope of the geological strata takes place, with the primary displacement direction leaning to the left, as illustrated in Figure 10b. The wellbore undergoes significant shear deformation at the bottom, while the upper part of the wellbore displaces concurrently with the landslide.

Figure 11 shows the direction of displacement, stress-strain distribution, and the deformation experienced by the entire wellbore. As presented in Figure 11a, the full length of the wellbore displays leftward and downward displacement. However, the displacement at the wellbore’s end is lesser, given that it is situated at the edge of the geological displacement, while the lateral displacement of the upper geological strata is, generally, more substantial. Consequently, under geological thrust, the wellbore experiences compression and deformation in the lower boundary region of the gas hydrate layer.

Figure 11b,c exhibit the strain and stress distribution within the wellbore, respectively. The maximum strain occurs in the 30 m section at the base of the wellbore, and there is a complex coexistence of tensile and compressive stress distributions with the maximum principal stress distribution. This significantly impacts the structural stability of the wellbore. Hence, as illustrated in Figure 11d, the entire wellbore displays noticeable deformation in the z-y plane.

4.4. Different Pressure Drops on Wellbore Stability

Figure 12 presents the distribution of hydrate saturation observed during a one-year production period of hydrates, utilizing various pressure drops. It is clear that the lower the production pressure, the more hydrates are produced in the same timeframe, and the broader the range of hydrate decomposition within the formation. This implies that the extent of deformation in the formation is greater, and the risk of underwater landslides increases, leading to a heightened instability in the deformation of the wellbore.

Figure 13a depicts the maximum formation subsidence and maximum wellbore deformation obtained under various pressure drops. As illustrated, when the bottom hole pressure is at 6 MPa, the formation experiences a maximum subsidence of 7.54 m. With the increase in bottom hole pressure, the maximum subsidence of the formation within the studied area demonstrates an almost linear decline. When the production pressure is at 10 MPa, the maximum formation subsidence is 4.43 m. Contrarily, the maximum displacement of the formation in the wellbore area escalates with the pressure drop. Considering the total volume of hydrate decomposition during production at various pressures, we can make an initial analysis that the maximum displacement of the wellbore formation is directly proportional to the absolute amount of hydrate decomposition. However, the precise pattern still requires further exploration.

Figure 13b illustrates the comprehensive wellbore deformation while producing hydrates under various bottom hole pressures during subsidence. It is clear that the lower the production pressure, the more pronounced the wellbore deformation and the larger the displacement. It is also noteworthy that, even though the wellbore typically experiences substantial displacements under these conditions, the displacement distribution along the wellbore remains roughly uniform, without significant abrupt changes in displacement. Therefore, under these circumstances, it is fair to consider that the wellbore can still maintain a safe level of production. However, the stability of the wellhead equipment and its secure connection with the wellbore warrant further verification.

Figure 13c portrays the comprehensive wellbore deformation while producing hydrates under different bottom hole pressures during landslides. It is apparent that, regardless of the wellbore’s production pressure, severe deformation will occur in the hydrate section of the wellbore when a landslide takes place, resulting in an overall failure of the development process. Therefore, the key to safe and efficient hydrate development lies in how to avoid landslides or the measures to be taken to prevent landslides during hydrate development.

5. Discussion

This article is based on Abaqus secondary development and establishes a submarine slope formation model containing hydrates and a developed wellbore. It realizes the multi−field coupling of heat-flow-solid and hydrate phase transformation and analyzes the deformation of the wellbore during hydrate decomposition. The risk of wellbore instability at different stages is discussed, and the impact of hydrate production pressure on wellbore stability is evaluated.

Although this study extensively simulates wellbore and slope stability under wellbore extraction conditions, it does not fully capture the entire process of gas and fluid flow out of the wellbore. Additionally, the wellbore structure has been significantly simplified. In future analyses, emphasis should be placed on the wellbore structure and horizontal wells extending up to kilometers in length. Furthermore, full dynamic processes of the slope and wellbore during multiphase flow should be a focal point. These considerations are crucial for the large-scale extraction of hydrates in the future.

6. Conclusions

(1) Minimal hydrate decomposition during drilling results in minor formation deformation near the wellbore, confined within 10 m, ensuring overall stability;

(2) Hydrate decomposition enhances formation settlement and induces wellbore deformation. Despite uniform wellbore displacement, formation settlement and wellhead stability require careful monitoring;

(3) Submarine landslides, triggered by hydrate decomposition, cause severe damage to the formation and wellbore. Complex stress concentration at 30 m depth leads to significant wellbore instability. Advanced geotechnical surveys should be employed to identify stress concentration zones, followed by the implementation of stabilization measures;

(4) Reduced production pressure increases formation settlement and wellbore deformation during hydrate production. Post-landslide, major wellbore deformation occurs in compressed areas, necessitating preventative measures within the hydrate formation wellbore section. Maintaining wellbore and formation stability is crucial for ensuring a safer and more efficient production process.