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

: This research explores the geomechanical challenges associated with gas hydrate extraction in submarine slope zones, a setting posing a high risk of signiﬁcant geological calamities. We investigate slope and wellbore deformations driven by hydrate decomposition within a subsea environment. Utilizing Abaqus, a ﬂuid-solid-thermal multi-ﬁeld 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, inﬂicting 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.


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%, porefilling 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 2 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.

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 hy-drate decomposition experimental results and numerical simulation results and analyzing benchmark cases.

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.

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): 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 −1 to 58 kJ•mol −1 .
The phase equilibrium conditions for methane hydrate are shown in the following equation [33,34]: 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: where g h is the rate of decomposition of the natural gas hydrate, kg•m −3 •s −1 , g g is the rate of generation of natural gas, kg•m −3 •s −1 , g w is the rate of water generation, kg•m −3 •s −1 , 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 reac is the kinetic constant of hydrate dissociation, mol•m −2 •s −1 •Pa −1 , 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 2 •m −3 , Γ 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 Γ 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 mD, 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 3 kJ•m −3 •s −1 , and h r is the reaction heat of the hydrate decomposition per unit amount of substance, 53 kJ•mol −1 .

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]: where: K is the layer volume modulus of the natural gas hydrates, Pa, K sh1 is the layer volume modulus of the natural gas hydrates when saturation is 1, Pa, K sh0 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 sh1 is the layer shear modulus of the natural gas hydrates when saturation is 1, Pa, G sh0 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 sh1 is the internal cohesion of the natural gas hydrates layer when saturation is 1, Pa, and c sh0 is the internal cohesion of the natural gas hydrates layer when saturation is 0, Pa.

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: where σ * is the effective stress, n t is the proportion of wet, σ denotes the skeleton effective stress, and 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: where σ and f 1 represent the total stress matrix and the body-force vector, respectively, t is the surface traction vector in unit area, δv and δε 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 ρ 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]: 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.

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.
1, Pa, and  sh0 is the internal cohesion of the natural gas hydrates layer when saturation is 0, Pa.

Analysis of Porous Media in Abaqus
The effective stress principal is adopted in Abaqus to express the basic stress deter mining rules for porous media, shown as: where  * is the effective stress,  is the proportion of wet,  denotes the skeleton effec tive stress, and ̅ 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 us ing the virtual work principal: where  and   represent the total stress matrix and the body-force vector, respectively  is the surface traction vector in unit area,  and  represent the virtual velocity field and the matrix of virtual strain rate, respectively,  is the liquid saturation,  is the liquid volume ratio to pore volume, and  is the density of liquid; (2) Determining the constitutive behavior of both the fluid and skeleton structur within the porous medium [8]; (3) Ensuring the continuity of the wetting liquid phase throughout the porous me dium, as described [19]: where  and  separately denote the volumes of unbound liquid and bound liquid and  is the proportion of unbound liquid in a unit volume.

Flowchart of the Methodology
Figure 1 presents a flowchart outlining the methodology.The yellow section illus trates how hydrate saturation is updated within the computational domain, a crucial as pect of the simulation.The green box denotes the primary theoretical framework em ployed in this method.The blue box highlights the macroscopic objects, loads, and ulti mate research objectives.

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.
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.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 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.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.

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.

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.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 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 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.

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.

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.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.

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

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 −3 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 −4 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.
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 −4 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.

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

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.
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.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 Figures 6c and 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.
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.

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.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.

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.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.

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.

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 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.
along the wellbore remains roughly uniform, without significant abrupt changes in dis placement.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 equip ment and its secure connection with the wellbore warrant further verification.Figure 13c portrays the comprehensive wellbore deformation while producing hy drates under different bottom hole pressures during landslides. is apparent that, regard less 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.

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 in stability 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 well bore extraction conditions, it does not fully capture the entire process of gas and fluid flow 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.

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

Figure 2 .
Figure 2. Gupta et al.'s 2D natural gas hydrates experiment and simulation conditions: (a) pressure drop and gas production experimental conditions; and (b) numerical simulation boundary conditions [22].

Figure 2 .
Figure 2. Gupta et al.'s 2D natural gas hydrates experiment and simulation conditions: (a) pressure drop and gas production experimental conditions; and (b) numerical simulation boundary conditions [22].

Figure 3 .
Figure 3.Comparison of simulation results between this model and other simulation experimental results: (a) cumulative gas production; and (b) total volumetric strain.

Figure 3 .
Figure 3.Comparison of simulation results between this model and other simulation experimental results: (a) cumulative gas production; and (b) total volumetric strain.

Figure 4 .
Figure 4. Seismic reflection images of the South China Sea submarine slope.

Figure 4 .
Figure 4. Seismic reflection images of the South China Sea submarine slope.

Figure 5 .
Figure 5. Model grid division.The red area represents the wellbore section, the blue area represents the hydrate layer, and the gray area represents the formation.

Figure 5 .
Figure 5. Model grid division.The red area represents the wellbore section, the blue area represents the hydrate layer, and the gray area represents the formation.

Figure 10 .
Figure 10.(a) Distribution of equivalent plastic strains within the geological layers and (b) distribution of displacement of the geological layers during the landslide.

Figure 10 .
Figure 10.(a) Distribution of equivalent plastic strains within the geological layers and (b) distribution of displacement of the geological layers during the landslide.

Figure 11 .
Figure 11.Distribution of (a) displacement, (b) maximum principal strain, (c) maximum principal stress, and (d) displacement contour lines of the wellbore.

Figure 11 .
Figure 11.Distribution of (a) displacement, (b) maximum principal strain, (c) maximum principal stress, and (d) displacement contour lines of the wellbore.

Figure 12 .
Figure 12.Hydrate saturation distribution cloud map when using different pressure drops for 1 year of production.

Figure 12 . 21 Figure 13 .
Figure 12.Hydrate saturation distribution cloud map when using different pressure drops for 1 year of production.J. Mar.Sci.Eng.2023, 11, x FOR PEER REVIEW 18 of 21

Figure 13 .
Figure 13.(a) Maximum subsidence and displacement in wellbore after 1 year production with var ying pressure drops, (b) overall wellbore deformation under formation subsidence, and (c) overal wellbore deformation under formation subsidence and landslides.

Figure 13 .
Figure 13.(a) Maximum subsidence and displacement in wellbore after 1 year production with varying pressure drops, (b) overall wellbore deformation under formation subsidence, and (c) overall wellbore deformation under formation subsidence and landslides.

Table 1 .
Parameters used in the experiment.

Table 1 .
Parameters used in the experiment.

Table 2 .
Parameter for the numerical model.