Influence of hydrate exploitation on stability of submarine slopes

The decomposition of natural gas hydrate will reduce the cementation effect of hydrate and produce ultra-static pore pressure, which will change the mechanical characteristics of the reservoir. Eventually, a series of geological disasters could be triggered, of which the submarine landslide is a typical example. In order to analyze the stability of hydrate-bearing submarine slopes and to explore the internal relationship between hydrate decomposition and submarine landslides, a “two-step reduction method” was described in this paper. This method was based on a strength reduction approach, which can be used to assess the effects of the initial geostress balance and hydrate decomposition on substrate strength reduction. This method was used to reveal the essence of hydrate decomposition, and then, a joint operation mode of multi-well was proposed. The internal relationship between hydrate decomposition and submarine landslides was analyzed in detail. And the development process and mechanism of submarine landslide were deeply discussed. The results showed that hydrate decomposition is a dynamic process of stress release and displacement, where the “stress inhomogeneity” distributed along the slope is transformed into “displacement inhomogeneity.” We concluded that hydrate decomposition could trigger a submarine landslide, especially along a sliding surface. The formation of submarine landslide is a gradual development process and presents the dual characteristics of time and space.


Introduction
Natural gas hydrate (NGH) is a non-stoichiometric caged crystalline substance formed by natural gas and water molecules under high pressure and low temperature (Englezos 1993). It is widely distributed in global offshore continental shelf environments and the onshore permafrost zones . Due to its huge resource reserves, NGH is considered an alternative energy source of great potential for the twenty-first century (Thakur 2010). The exploitation of hydrates has been successfully tested in several countries, in order to assess the technical feasibility of hydrate exploitation. Russia first exploited NGH in Messoyakha field, Siberia, in the last century (Makogon et al. 2005). At the end of 2011, the accumulative gas obtained by hydrate decomposition reached 5.4 × 10 9 m 3 (Li et al. 2016). In 2007, Canada adopted the depressurization method for NGH production, but the production only lasted for 60 h due to sand production (Chand and Minshull 2004;Song et al. 2019b). In 2012, the USA successfully exploited NGH in Alaska permafrost region via CO 2 injection (Boswell et al. 2012;Schicks et al. 2018). In 2013, Japan achieved the successful trial production of oceanic NGH for the first time, and the total gas production reached 1.2 × 10 4 m 3 (Kimura et al. 2019). Since then, China successfully tested and produced NGH in the northern South China Sea in 2017 and 2020 Sun et al. 2021;Zhu et al. 2021), establishing world records in the process for total gas production and average daily gas production. They demonstrated that muddy silt hydrate resources could be exploited in a safe and controlled manner. Figure 1 shows the second NGH production test in China.
Production methods for hydrate include heating, depressurization, chemical reagent injection, carbon dioxide replacement or a combination of methods (Koh et al. 2002;  Moridis et al. 2007;Wang et al. 2014). Although the methods are different, the essence of hydrate exploitation is hydrate decomposition (Tupsakhare et al. 2017;. However, when the hydrate decomposes, the contribution made by hydrate as an agent of cementation also disappears and excess pore pressure is generated. The strength of the hydrate and hence surrounding strata decreases, which can trigger a submarine landslide or submarine turbidity current and cause wellbore instability and platform overturning (Sultan et al. 2004;Rutqvist et al. 2009;Li et al. 2014;Liu et al. 2017;Li and Han 2021;Zhang et al. 2021). In addition, one of these geological events may also trigger another, initiating what has been called a "disaster chain" (Lu et al. 2019). With the continuous expansion of the scale of an event, long-term and large-scale geological disasters can be induced. A submarine landslide, which is associated with hydrate decomposition, is a typical example (Mountjoy et al. 2014).
Terrestrial slope stability has been intensively studied, with the limit equilibrium method and numerical analysis methods widely used to analyze slope stability (Zaniboni et al. 2016;Luo et al. 2018;Li and Yang 2019;Wang et al. 2019;Liu et al. 2021). At the same time, a series of detailed research methods have been developed based on different considerations. Examples of the latter include the simplified Bishop method (Bishop 1955;Cheng and Yip 2007), Janbu method (Tang et al. 2020) and the upper bound method of limit analysis based on plastic mechanics . In terms of numerical analysis, the strength reduction method is well known and has been widely used (Bolla and Paronuzzi 2020;Wang et al. 2020). In contrast to terrestrial slopes, submarine slopes have a more complex geological environment, such as the influence of overlying seawater pressure, continuous seepage, biological disturbance and submarine volcanism. Meanwhile, they still have some similarities, the most obvious of which is the similar slope forms and the stability risks they all face. As a result, research methods similar to those developed for terrestrial slopes can be used for analysis of the submarine slopes. Tan et al. (2021) analyzed the stability of submarine slopes during the exploitation of gas hydrates, within the framework of the limit equilibrium method. Leynaud et al. (2004) used the limit equilibrium and finite element methods to study the static and dynamic (seismic) stabilities of the dipping seabed. Song and co-workers (Song et al. 2019a, b) investigated the stability of hydrate-bearing submarine slopes (HBSS) at different gradients using the strength reduction method. Kong et al. (2018) carried out multi-factor sensitivity analysis of submarine slopes based on strength reduction method. Zheng et al. (2019) combined limit analysis with strength reduction to discuss the influence of wave pressure on the stability of submarine sensitive clay slopes. Chen et al. (2020) established the finite element model of wave-seabed interaction and evaluated the stability of artificial submarine slopes based on strength reduction method.
A relatively limited amount of research has been carried out on the stability of submarine slopes affected by hydrate decomposition. There are few systematic researches on the internal relationship between hydrate decomposition and submarine landslides. In addition, although the stability of HBSS has been explored by some scholars (Kwon and Cho 2012;Jiang et al. 2015;Elger et al. 2018;Kong et al. 2018;Song et al. 2019a), most of them made qualitative assessments based on the obtained safety factor, while the influence of hydrate decomposition before a landslide was ignored. Some studies have not even considered the influence of initial geostress balance on subsequent hydrate decomposition and submarine landslide. In order to investigate the influence of hydrate decomposition on submarine slope and analyze the development process of submarine landslide after hydrate decomposition, we propose a "two-step reduction method" which can consider the initial stress balance, hydrate decomposition and strength reduction. The influence of hydrate decomposition on the stability of submarine slopes was comprehensively studied using this method. In this study, the potential effects of hydrate decomposition were revealed and reasonable engineering solutions proposed. The processes and mechanisms associated with submarine landslides were analyzed in detail.
2 Two-step reduction method

Strength reduction method
The strength reduction method was first proposed by Zienkiewicz (1975). As a method to evaluate slope stability evaluation, it has been widely used by many scholars due to its simplicity and the ease with which it can be understood. The reduced shear strength parameter can be expressed as: where c and φ are the shear strength that soil can provide; c m and φ m are the shear strength needed to maintain balance or actually exerted by soil; and F r is the strength reduction factor.
In the finite element analysis, different strength reduction factors (F r ) are assumed first and then calculated and analyzed according to the reduced strength parameters. During the calculation process, the value of F r increases continuously until critical failure is reached. The strength reduction factor (F r ) at critical failure is then the safety factor of slope stability.

Two-step reduction method
When the strength reduction method is used, the initial geostress should be balanced first and then for the strength to be reduced. For most research objects, such as slopes, hills and mines, they have reached a stable state under the long-term geological action. Therefore, it is reasonable to think that their mechanical parameters remain unchanged. When carrying out strength reduction, it is only necessary to balance the initial geostress and then carry out strength reduction. However, for some geological objects (slopes, hills, mines, etc.) that have been affected by recent disturbances (e.g., earthquake, blasting, etc.), their mechanical parameters would have changed to an unknown extent. The following issues therefore need to be taken into consideration in relation to disturbances: (1) For an unspecified time after a disturbance, the new geostress field would not have stabilized. If new mechanical parameters were given to geological attributes in this post-disturbance transitional stage of stress field evolution, it would obviously be inconsistent with the actual situation when the stress field is in a stable state. This scenario (transitional stage) could not be used to model the impact of hydrate decomposition on slope stability. (2) Sometimes, we need to analyze the disturbance process or disturbance influence, which occurs after the balance of initial geostress.
(1) c m = c∕F r m = arctan(tan ∕F r ) 1 3 (3) The change of model parameters caused by geological disturbance may be local or global, thereby making the simulation process and interpretation of any results difficult.
If the limitations listed above are taken into account and combined with the strength reduction implementation of ABAQUS (Sun et al. 2016;Schneider-Muntau et al. 2018;Song et al. 2019a), a simple and effective "two-step reduction method" can be used to assess slope stability. The main steps are as follows: (1) Firstly, the field variable is added, and then, the mechanical parameters that need initial balance are set to correspond to the initial value of the field variable (the initial value of the field variable should be less than 1). (2) According to the actual situation, the mechanical parameters affected by the disturbance are given to the whole or part of the strata of the geological object and set to the corresponding field variable value of "1." (3) Strength reduction is carried out on the basis of the mechanical parameters corresponding to the field variable value "1," and the mechanical parameters corresponding to each field variable value are given. (4) Establish the first analysis step (initial balance) and adopt the mechanical parameters corresponding to the initial values of field variables. (5) Establish the second analysis step and reduce the strength to the field variable value "1." (6) Establish the third analysis step and reduce the strength to the final value of the field variable.
When the influence of different degrees of disturbance is analyzed, or the rock stratum in a local area is disturbed, it is only necessary to modify the mechanical parameters corresponding to the variable value "1" of the disturbed rock stratum field. The two-step reduction method can be expressed by formula (2): where c and φ are the shear strengths that undisturbed soil can provide, while c 1 and φ 1 are the shear strengths that the disturbed soil can provide. F(c) and F(φ) represent the relationship between undisturbed and disturbed soil shear strength. The values for F(c) and F(φ) are determined by the actual situation or research needs. c m and φ m are the shear strength needed to maintain balance, or the shear force exerted by the soil. F r is the strength reduction factor.
To sum up, the "two-step reduction method" can be classified as a type of "variable strength stability analysis" that can be used to assess the potential impact of disturbances on slope sediments (e.g., earthquake, blasting, hydrate decomposition, preloading). It is especially suitable for comparative analysis of overall stability, which is accompanied by local strength changes. The main advantages of the "two-step reduction method" are: (1) Avoidance of repeated geostress balance and mechanical parameter assignment. (2) (2) It can be used to assess slope stability for different degrees of disturbance and different disturbance forms (linear and nonlinear). (3) It can be used to assess the effect small-scale changes in the geostress field on slope stability. (4) It can be used for a comprehensive analysis of initial stress balance (the first analysis step), the effect of disturbance (the second analysis step) and strength reduction (the third analysis step) through a single ".cae" file containing multiple steps for analysis.

Slope instability criteria
Generally, there are three criteria for slope instability (Masson et al. 2011;Song et al. 2019a): (1) Formation of a continuous plastic zone (plastic zone penetration).
(2) Sudden change of displacement at the top of a slope.
In the process of strength reduction, the three criteria show a progressive relationship in time. First, the plastic zone is formed at the bottom of the slope and develops upward, then the plastic zone penetrates, and finally the displacement of the top of the slope changes suddenly. With a further reduction of strength, the plastic zone expands continuously, accompanied by increased displacement and ultimately large-scale instability of the slope. This is represented by non-convergence of the calculation. As long as the consistent criterion is adopted for an actual calculation, the safety factor obtained can be used for qualitative and quantitative analysis. Considering that the first two criteria need to be combined with cloud images for artificial judgment, this paper has used the interruption of calculation due to non-convergence as the consistent criterion for judging slope instability.

Applicability of the two-step reduction method
The focus of this paper was the stability of HBSS. The paper investigated the influence of hydrate decomposition on the stability of submarine slopes exposed to seawater pressure (depth) and gravity. The decomposition of NGH will decrease the mechanical strength of the sediments comprising the hydrate reservoir. Therefore, the key problem to be considered is the change of mechanical properties of hydrate reservoir when the decomposition degree of hydrate is different. When the strength reduction method is used to analyze the slope after hydrate decomposition, the cardinal number "1" (the overall mechanical parameters of the macroscopic submarine slope) corresponding to different degrees of decomposition is also different. In addition, the hydrate reservoirs with different degrees of decomposition in the natural environment are in the same state before decomposition (all in equilibrium state under overlying seawater pressure and self-weight stress). Therefore, it is impossible to simply give the mechanical parameters after hydrate decomposition to the reservoir and then carry out initial balance and strength reduction. By using the approach described above, the stability of HBSS can be solved simply and effectively by adopting a "two-step reduction method."

Site description
The Shenhu Sea area is located in the Baiyun Sag in the South China Sea (as shown in Fig. 2). In 2017 and 2020, China tried to produce hydrates in this region and established world records in the process for production duration of 60 days, total gas production of 8.6 × 10 5 m 3 and daily average gas production of 2.87 × 10 4 m 3 Sun et al. 2021;Zhu et al. 2021). The average slope angle in the area is 3.3-3.8°, while the maximum slope angle can reach 30°. According to the seismic reflection profile and the location of the BSR (bottom simulating reflector), the dimensions of the hydrate reservoir such as depth of burial and thickness could be determined. The depth of seawater in this area can exceed 1000 m. According to the lithology logs for the Shenhu drill hole, the distance from the seafloor to the top of the hydrate-bearing sediment ranges from 50 to 350 m. The thickness of the hydrate-bearing strata ranged from 10 to 120 m.

Dimensions of the numerical model
Based on the actual geological conditions, referring to the relevant literature and making appropriate optimization, the dimensions of model are finally determined. The numerical model had a total length of 2590 m, a slope length of 1600 m and horizontal submarine planes with lengths of 500 m above the slope and 500 m beyond the toe of the slope (Fig. 3). The slope angle was set at 6°, with the minimum height set at 700 m (Fig. 3). The trend (slope) of the hydrate reservoir was assumed to be parallel to the surface of the submarine slope (Fig. 3). The reservoir was assumed to be 60 m thick and the depth of burial 200 m. In addition, the model was divided into overlying strata (pink strata in Fig. 3

Constitutive model
The Mohr-Coulomb model is widely applicable in many fields such as geotechnical engineering. It is not only simple in form but also can meet the needs of engineering. In this paper, a "two-step reduction method" based on strength reduction was used for simulation analysis (as described in Sect. 2). Because the reduced mechanical parameters only relate to cohesion c and internal friction angle φ (as shown in Formula 2), and considering the high uncertainty of coupling analysis. It is assumed that water and gas can be quickly discharged during hydrate decomposition. The decomposition process of sandy reservoir is simulated without considering the influence of pore pressure. This is also supported by other scholars (Kong et al., 2018). Therefore, the Mohr-Coulomb model selected in this paper is completely reasonable and meets the needs. Other secondary contents such as integration method are not the focus of this paper. They are not discussed here, and the default options of ABAQUS can be adopted.

Model parameters
The mechanical properties of hydrate reservoir are closely related to hydrate saturation. There is a significant difference in strength between the reservoir and the surrounding strata. Triaxial tests and hydrate decomposition tests with different hydrate saturations (Miyazaki et al. 2011;Hyodo et al. 2013;Luo et al. 2020) show that prior to decomposition, the strength of the hydrate reservoir was stronger than that of the surrounding strata due to the cementation and filling of hydrate. And the higher the hydrate saturation, the higher the strength of hydrate-bearing sediments. Therefore, we have reasons to believe that with the increase of decomposition degree, the strength of hydrate-bearing sediments gradually decreases. This assumption is also a common practice for hydrate simulation (Kong et al. 2018;Song et al. 2019a, b). Based on the experiments of studying the mechanical properties of hydrate-bearing sediments and related simulation researches (Zhang et al. 2010;Ghiassian and Grozic 2013;Shi et al. 2015Shi et al. , 2019Kong et al. 2018;Choi et al. 2020), the mechanical parameters of the strata were determined after trial calculation, with the hydrate saturation about 25%. The model is a two-dimensional plane model, with a quadrilateral plane strain element. The mesh of hydrate reservoir and the area connected to it were refined. The model size and dimensions are shown in Fig. 3. The mechanical parameters of the strata are presented in Table 1.

Boundary conditions
The boundary conditions should correspond to the simulated geological conditions. The boundary conditions of the model in this paper are as follows: Horizontal displacement is limited at the left and right boundaries of the model, with displacement in all directions limited at the bottom boundary of the model. The load on the model included the overlying seawater pressure (depth) and its own gravity, with a simulated water depth of 1000 m. The division of strata and boundary conditions of the model are shown in Fig. 4.

Solution process
As mentioned earlier, the multiple analysis steps based on "two-step reduction method" was adopted to calculate in this paper. Different analysis steps have corresponding research objectives and show a progressive relationship, which will be further introduced here. The whole solution process was divided into three sub-solution processes: initial geostress balance, hydrate decomposition and landslide process (safety factor calculation). First, balance the initial geostress. The purpose of this analysis step is to restore the original in situ stress state of HBSS under long-term geological action. This is also a prerequisite for subsequent analysis of influence of hydrate decomposition. Then, the influence of hydrate decomposition was analyzed (the first analysis step of strength reduction). In this analysis step, the influence of different decomposition degrees on HBSS can be obtained via providing mechanical parameters corresponding to different decomposition degrees. Finally, the analysis of development process of submarine landslide was carried out (the second analysis step of strength reduction). It should be noted that the strength reduction of the whole model is required in this calculation process, while only the strength reduction of hydrate reservoir was carried out in the previous analysis step. In addition, the last calculation process requires the model to automatically stop due to non-convergence, so as to obtain the safety factor of HBSS.

Influence of hydrate decomposition
Hydrate decomposition can be attributed to many potential factors, including earthquakes, volcanic eruptions, climate change and a fall in sea level (Max 2003;Huseynov and Guliyev 2004;Mau et al. 2007;Lapham et al. 2008;Dzyuba and Zektser 2013). These events cause significant changes in temperature and pressure in sediments, which can subsequently cause large-scale decomposition of the hydrate (Lu et al. 2019). At the same time, considering the model size and the effective verification of the "two-step reduction method," the hydrate decomposition in this paper is assumed to occur evenly throughout the whole the reservoir.

Displacement change caused by hydrate decomposition
The displacement nephogram with different degrees of decomposition is shown in Fig. 5. The illustrations in Fig. 5b-d are unified to highlight the displacement caused by the increase in the degree of decomposition. The automatically generated legend in Fig. 5a can be used to highlight the distribution of the magnitude of displacement in different areas of the nephogram. It can be seen from Fig. 5 that the decomposition of hydrate causes a significant change in displacement. Interestingly, with the hydrate reservoir as the boundary, the upper displacement is wedge-shaped and the lower displacement is arc-shaped. Moreover, the peak displacement of the upper strata is concentrated near the top of the slope, while the displacement decreases closer to the foot of the slope. With the increase of degree of decomposition, the wedge-shaped area spreads from the top to the foot of the slope and the overall displacement increases continuously. The wedge-shaped distribution of displacement and the characteristics of large displacement at the top of the slope and smaller displacement at the foot of slope are closely related to the slope gradient. In addition, the "arc-shaped" distribution of the lower strata is very similar to the "arc-shaped" sliding surface of the landslide. Therefore, the formation of submarine landslides in hydrate-bearing sediments is directly related to the displacement changes caused by hydrate decomposition.
In order to obtain insight into additional changes in displacement caused by hydrate decomposition, two survey lines S1 (sloped) and S2 (vertical) are drawn in Fig. 5d. S1 is in the hydrate reservoir direction, while S2 is a vertical section through the reservoir. Figure 6 shows the displacement changes under different degrees of decomposition.
It can be clearly observed that the displacement curves for different degrees of decomposition show similar change rules. The slope displacement presents an inclined distribution with larger displacement at the top of slope and smaller displacement at the foot of slope. The peak displacement occurs at a certain distance from the top of the slope, not at the top of the slope. With an increase in the degree of hydrate decomposition, the overall displacement of the survey line increases and the peak value rises continuously. At the same time, it can be noted that the displacement of the area outside the hydrate reservoir in the horizontal direction changes little.
The displacement change in the vertical direction can be obtained from S2. In the vertical direction, the change in displacement caused by hydrate decomposition is in two different states, with the hydrate reservoir as the boundary. The displacement of strata above the hydrate reservoir is greatly affected by hydrate decomposition, with the extent of displacement positively related to the degree of hydrate decomposition. The displacement peak appears at the top of the model and decreases with an increase in the depth of burial. With the decrease of displacement, there is an obvious steep drop in the extent of displacement within the hydrate reservoir. The displacement of strata below the hydrate reservoir changes little and hence not by the same order of magnitude as that of the upper part, although it still shows a decreasing trend from top to bottom.
The distribution of displacement vectors with degree of decomposition set at 50% is shown in Fig. 6c. The largest displacements are concentrated in the area above the hydrate reservoir, with the displacement at the top of the slope larger than that at the foot of the slope. In terms of the direction of displacement, the change in displacement caused by hydrate decomposition is dominated by vertical displacement, which leads to noticeable subsidence of the surface of the seabed. At the same time, due to the inclination of the slope, it will also produce horizontal displacement to the right.
To sum up, the range of influence of hydrate decomposition is limited. In the horizontal direction, the extent of influence is limited to the width of the hydrate reservoir. In the vertical direction, the extent of influence is largely limited to the hydrate reservoir and strata above it. Other areas in the model are less affected by hydrate decomposition.

Stress change and strain caused by hydrate decomposition
Hydrate decomposition is a dynamic process, with displacement changing in response to the redistribution of stress. The change in stress is not as obvious in the nephogram as the extent of displacement caused by hydrate decomposition (as shown in Fig. 7a).
The size of the hydrate reservoir is relatively small in comparison with the model as a whole model. The change of local stress caused by hydrate decomposition is much smaller than that of the whole model, so the change is not obvious in the nephogram a Displacement change of S1  1 3 as a whole. The stress analysis of the S2 line provided a better means to understand the nature of changes in the reservoir. The stress variation curve of S2 is shown in Fig. 7b. It can be seen that stress shows strong regularity. In the lateral direction, the hydrate reservoir area is greatly affected by hydrate decomposition, while the area outside the hydrate reservoir is much less affected. It is precisely because of the different degree of influence that significant stress and abrupt change occurs at the boundary between the hydrate reservoir and surrounding strata. In the hydrate reservoir area, with the increase in the degree of decomposition, the stress decreases continuously, which is opposite to the change law of displacement. At the same time, it can be noted that when the degree of decomposition of hydrate is low (25%), the stress is characterized by a large value at the top of the slope and a small value at the foot of the slope. That is, there is a "non-uniform" distribution of stress in terms of magnitude along the slope. The stress peak appears at a certain distance from the top of the slope, which coincides with the displacement peak. With a further increase in the degree of hydrate decomposition, the "non-uniform" distribution of stress gradually disappears and the stress curve gradually tends to be flat. Therefore, hydrate decomposition has the effect of "dilution and homogenization" on the "non-uniform" distribution of slope stress.
Elastic strain is not the main cause of material deformation and failure, while the tensile strength of geotechnical materials is very small. Inelastic shear strain was therefore a Stress nephogram with degree of decomposition (50%) b Stress variation of S2  selected for analysis in this study. At the same time, the true principal strain is used as a reference for comparative purposes. The distribution of strain with a degree of decomposition set at 50% is shown in Fig. 8. Figure 8a is the nephogram of inelastic shear strain, while Fig. 8b represents the distribution of the true principal strain vector. As shown in Fig. 8a, inelastic shear strain was significant in the hydrate reservoir area, while continuous strain below the reservoir can be regarded as the potential slip surface. The peak value of inelastic shear strain appears at both ends of the hydrate reservoir, that is, at the transition zone between the hydrate reservoir and the non-hydrate strata. Combined with the sudden change of stress in Fig. 7b, it can be seen that the strata are subjected strong shear stress at the boundary of the hydrate reservoir and surrounding strata. As a result, shear failure could easily occur in this area. According to Fig. 8b, it can be clearly observed that the peak value of principal strain appears at the left end point of the hydrate reservoir, which is far greater than the principal strain value at other positions. In addition, the end point is connected with the principal strain at the top boundary of the model. Comparing Fig. 8a with Fig. 8b, the inelastic shear strain is highly consistent with the true strain distribution. This high level of concordance indicates that inelastic shear strain is probably the dominant type of strain caused by hydrate decomposition. Stated differently, inelastic shear strain can better characterize the failure of strata than other strain types. As the decomposition of the hydrate becomes more pronounced, the strain (deformation) in different regions of the model increases simultaneously. This synchronous increase is likely to cause continuous deformation and shear failure in local areas. Finally, a through deformation zone (plastic failure zone) is formed, which further forms a continuous sliding surface and induces a large-scale submarine landslide.

Essence of hydrate decomposition
Through the comparative analysis of displacement and stress, it can be seen that the decomposition of hydrate is a dynamic process of stress release and displacement increase, with the magnitude of displacement increase related to the degree of stress release. Due to the inclined nature of a submarine slope, the stress distribution along the slope after geostress balance is non-uniform. The stress values are greater at the top of a slope and smaller at the foot of a slope. In the process of hydrate decomposition, stress is released and its "non-uniform" distribution along the slope disappears. However, the displacement distribution is "non-uniform" again and gradually decreases from the top to the foot of the slope, which corresponds to the released stress. Therefore, the main effect of hydrate decomposition is to transform the "stress inhomogeneity" distributed along the slope into "displacement inhomogeneity," with the displacement positively correlated with the stress difference before and after decomposition.
The patterns and processes discussed above can be explained in terms of momentum. Formula (3) can be used to describe the forces acting on a microelement of the hydrate reservoir during the decomposition of hydrate as follows: where F is the average combined force on the micro element during hydrate decomposition; Δt is the time of hydrate decomposition; dm is the mass of the micro element; and Δv is the variation of average velocity of the micro element during hydrate decomposition.
It can be seen from Fig. 5 that the hydrate decomposition process is accompanied by an increase in displacement, that is, the combined force changes. The dm of each microelement is almost equal (the mass of discharged water and gas is very small relative to the mass of the reservoir), the Δt is the same, and Δv is determined by F. According to Fig. 7b, before hydrate decomposition, the stress along the slope is larger at the top and smaller at the foot, but after hydrate decomposition, the stress along the slope is almost equal. Therefore, in this scenario, the F on the microelement from the foot of the slope to the top of the slope gradually increases. That is, Δv gradually increases from the foot of the slope to the top of the slope, and the closer to the top of the slope, the greater the displacement in the same period of time.

Suggestions for hydrate exploitation
According to the foregoing analysis, the range of displacement and changes in stress caused by hydrate decomposition are limited. In the horizontal direction, the extent of influence is the width of hydrate reservoir. In the vertical direction, it extent of influence is largely confined to the hydrate reservoir and strata above the hydrate reservoir. Especially near the top of the slope and in the transition area between the hydrate reservoir and non-hydrate strata, there will be a large change in displacement or shear stress. At present, a vertical hole is mainly used for gas hydrate test production (Yu et al. 2019b). Therefore, for low-risk gas hydrate production in the actual submarine slope area, the hole should be located at the foot of the slope, or far away from the slope area, instead of near the top of the slope.
In addition, NGH exploitation may need to be carried out in a submarine slope with poor stability, which will lead to a submarine landslide. At this time, the joint exploitation mode of "vertical shaft + horizontal well" in combination with multi-branch (3) F ⋅ Δt = dm ⋅ Δv holes could be considered, which can significantly mitigate the risk of slope failure. In this multi-well operational mode, horizontal wells should be set in the lower area of the hydrate reservoir, while the vertical wells should be set away from the slope area. This multi-well joint operation mode is also the development trend for future hydrate exploitation (Wilson et al. 2011;Yu et al. 2019a). The multi-well joint operational mode is shown in Fig. 9.

Landslide development process
By adopting the "two-step reduction method," after analyzing the effect of hydrate decomposition, it is convenient to reduce the overall strength of a submarine slope through a second reduction analysis step and then to obtain the development process of a submarine landslide. Figure 10 shows the development process for the plastic zone with different reduction factors, when the degree of decomposition is 75%.
As shown in Fig. 10, the development process of submarine landslide presents a strong regularity and shows an obvious progressive relationship in stages, and hydrate decomposition can induce submarine landslides. First of all, when hydrate decomposition is finished, the plastic damage caused by hydrate decomposition is mainly manifested at the boundary between two ends of the hydrate reservoir (the boundary between the hydrate reservoir and non-hydrate reservoir), and the left boundary area is connected with the horizontal submarine plane at the top of the model (as shown in Fig. 10a). With the weakening of the overall strength of the submarine slope, the plastic zone develops further. An obvious plastic failure zone extends downward from the junction at the left end of the hydrate reservoir (as shown in Fig. 10b). If Fig. 10b is compared with Fig. 8, it can be seen that the development of the plastic zone is highly consistent with the shear failure and principal strain distribution caused by hydrate decomposition, especially the potential slip surface, with which it basically coincides. This indicates that hydrate decomposition can induce a submarine landslide. In other words, the submarine landslide after hydrate decomposition is the further development of shear failure based on the influence of hydrate decomposition.
With a further increase of the reduction factor, the plastic zone at the right boundary of the hydrate reservoir connects with the right submarine plane at the top of the model (as shown in Fig. 10c). On this basis, the reduction factor continues to increase and the lower potential slip surface continues to expand, which connects with the plastic zone at the right end of the reservoir (as shown in Fig. 10d Fig. 9 Schematic diagram of a multi-well joint operation 1 3 overall connectivity is formed from the top to the foot of the slope and a landslide takes place. After that, the landslide continued to develop with the increase of reduction factor. When the reduction factor is 1.866, the maximum moving distance of the landslide Fig. 10 Landslide process on a submarine slope (decomposition by 75%) reaches 174 m. It can be seen that a large-scale submarine landslide has occurred at this time. When the model is stopped due to non-convergence, the final safety factor is 1.875. To sum up, the processes associated with submarine landslides, under the influence of hydrate decomposition, can be divided into the following stages: (1) First, hydrate decomposition leads to plastic failure at the interface between the hydrate reservoir and non-hydrate reservoir. The plastic zone at the junction near the top of the slope runs through the top of the slope.
(2) Second, the plastic zone at the junction near the top of the slope extends deep into the strata and forms a potential slip surface. (3) Third, the plastic zone near the junction of the foot of the slope extends to the top and connects with the horizontal submarine plane at the foot of the slope. (4) Then, the potential slip surface is connected with the plastic zone at the foot of the slope, forming a plastic zone from the top to the foot of the slope. (5) Finally, the penetrating plastic zone develops continuously and a large-scale submarine landslide takes place.
Through the above analysis, it can be known that hydrate decomposition will cause shear failure of the reservoir and adjacent strata. The subsequent submarine landslide is the ultimate result of further development and intensification of previous shear failure caused by hydrate decomposition. The mechanism of submarine landslide under the influence of hydrate decomposition can be summarized as follows: hydrate decomposition causes shear failure in local areas, especially at the boundary between reservoir and non-reservoir at both ends of the reservoir (based on the model in this paper). The existing shear failure is further developed, and the local failure areas are gradually connected with each other to form a through failure area with a large span. Then, the through failure zone develops into potential slip surface, which develops toward the top and foot of the slope. Eventually, a submarine landslide is induced.

Safety factor
The safety factors of submarine slopes under different degrees of decomposition were counted. At the same time, after hydrate decomposition (F r = 1), the peak value of the plastic zone before overall strength reduction of the submarine slope was also counted. The variation of the safety factor values and peak values for the plastic zone under different degrees of decomposition is shown in Fig. 11. According to Fig. 11, with an increase in the degree of hydrate decomposition, the safety factor gradually decreases, while the peak values of the plastic zone increase monotonically with an increase in the degree of decomposition. According to the relationship between plastic failure and safety factor, we can explain the relationship between the degree of hydrate decomposition and submarine slope stability from different perspectives. On the one hand, the increase in hydrate decomposition will inevitably lead to a decrease in hydrate reservoir strength. According to the principle of strength reduction, the decrease of cohesion and internal friction angle of the hydrate reservoir will lead to a decrease of "overall strength" of the macroscopic submarine slope. On the other hand, according to the analysis above, the decomposition of the hydrate will cause significant changes in stress and displacement, with the geostress redistributed as a consequence. Especially in the boundary area between hydrate reservoir and non-hydrate reservoir, it is accompanied by significant stress concentration and shear strain. This dual action of stress and displacement will have a profound impact on the submarine slope. In other words, plastic failure (landslide) will be induced by hydrate decomposition, or what could be called a "memory effect and historical effect" of strata stress. This comprehensive influence increased with an increase in the degree of hydrate decomposition.
The safety factor of hydrate changed from 1.883 to 1.862 during the transition from undecomposed to completely decomposed. This showed that the influence of hydrate decomposition on the overall stability of a submarine slope is relatively small, which is consistent with the conclusion that the range of influence of hydrate decomposition is limited. It should be pointed out that the influence of hydrate decomposition on submarine slope is controlled by the relative size of the hydrate reservoir and submarine slope. That is, the influence of hydrate decomposition is relative.
The displacement nephogram and vector distribution within the plastic zone just passing through from the top to the foot of the slope are shown in Fig. 12 (decomposition by 100%). In Fig. 12a, we can see that in the initial stage when the plastic zone runs from the top to the foot of the slope, a large range of strata move along the potential slip surface as a whole. At this time, the overall form of the landslide has taken shape. The displacement nephogram of the area above the hydrate reservoir is similar to that in Fig. 5, which still shows a wedge-shaped distribution of displacement near the top of the slope, rather than a circular arc-shaped distribution similar to the potential slip surface. This indicates that at this time, the potential slip surface is formed in the deep stratum (below the hydrate layer) and the landslide begins. However, the shallow stratum (above the hydrate reservoir) has not yet formed a continuous slip surface. At the same time, the whole is distributed in circular arc shape, while the upper part is distributed in a wedge shape, which can also be regarded as a precursor to the formation of a large-scale submarine landslide.
As shown in Fig. 12b, when the safety factor changes from only 1.851 to 1.852, the maximum displacement of the slope increases from 21.54 to 87.2 m. The displacement nephogram is arc-shaped, and a large-scale submarine landslide is formed. It can be seen Peak of plastic zone (F r =1)

Fig. 11
Variation curve of safety factor and peak value of plastic zone that a submarine landslide presents dual characteristics in time and space. In time, the change is extremely fast, while in space, it is accompanied by a large-scale change from the top to the foot of the slope. This is also the most fundamental reason why submarine landslides cannot be effectively monitored in the natural environment. In Fig. 12, the maximum displacement appears at the top of the sliding surface, with the top of the slope moving downward, while the foot of the slope moves upward under the action of the landslide. In addition, the displacement in this process is mainly horizontal. It should be pointed out that after the landslide begins to form (the plastic zone runs through), the large-scale displacement is mainly dominated by the strata sliding forming the continuous sliding surface and is no longer controlled by the strata strength.

Conclusions
(1) The proposed "two-step reduction method" can be used to consider the effects of the initial stress balance, hydrate decomposition and strength reduction on submarine slope stability. The success of the method meant the tedious steps of repeated initial balancing and parameter assignment could be avoided.
(2) Hydrate decomposition is a dynamic process of stress release that results in an increase in slope displacement. The main effect of hydrate decomposition is transformation of the "stress inhomogeneity" distributed along the slope into "displacement inhomogeneity." The displacement is positively correlated with the stress difference before and after decomposition. (3) When the stability of submarine slopes is low and hydrate exploitation takes place, submarine landslides can be triggered. The joint exploitation mode of "vertical well + horizontal well" combined with multi-branch holes should be considered to mitigate the risk. Horizontal wells should be set in the lower reaches of the hydrate reservoir, while vertical wells should be set away from the slope. (4) A submarine landslide, especially the formation of a sliding surface, will be induced by hydrate decomposition. (5) The final formation of submarine landslide is a gradual development process, while the submarine landslide is an event with dual characteristics of sudden change in time and large scale in space.