Submarine Slope Failure Primed and Triggered by Bottom Water Warming in Oceanic Hydrate-Bearing Deposits

Abstract

Many submarine slope failures in hydrate-bearing sedimentary deposits might be directly triggered, or at least primed, by gas hydrate dissociation. It has been reported that during the past 55 years (1955–2010) the 0–2000 m layer of oceans worldwide has been warmed by 0.09 °C because of global warming. This raises the following scientific concern: if warming of the bottom water of deep oceans continues, it would dissociate natural gas hydrates and could eventually trigger massive slope failures. The present study explored the submarine slope instability of oceanic gas hydrate-bearing deposits subjected to bottom water warming. One-dimensional coupled thermal-hydraulic-mechanical (T-H-M) finite difference analyses were performed to capture the underlying physical processes initiated by bottom water warming, which includes thermal conduction through sediments, thermal dissociation of gas hydrates, excess pore pressure generation, pressure diffusion, and hydrate dissociation against depressurization. The temperature rise at the seafloor due to bottom water warming is found to create an excess pore pressure that is sufficiently large to reduce the stability of a slope in some cases. Parametric study results suggest that a slope becomes more susceptible to failure with increases in thermal diffusivity and hydrate saturation and decreases in pressure diffusivity, gas saturation, and water depth. Bottom water warming can be further explored to gain a better understanding of the past methane hydrate destabilization events on Earth, assuming that more reliable geological data is available.

1. Introduction

Submarine slope failure is a fundamental shear failure process that occurs in submarine sediment formations [1,2,3]. Increase in shear stress or decrease in shear strength can lead to shear failure in sediments, which in turn leads to slope failure. The mechanisms that potentially cause shear failure include (a) removal of lateral supports, (b) increase in weight of overlying sediment, (c) generation of excess pore pressure (overpressure), (d) upward fluid flow, and (e) earthquakes [3]. Here, one mechanism, e.g., generation of overpressure, may simply prime a slope for failure before another one, e.g., earthquake, actually triggers slope failure [3].

1.1. Gas Hydrate Dissociation as a Trigger or Primer for Slope Failure

Gas hydrates are ice-like solid compounds in which gas molecules are trapped within cages of hydrogen-bonded lattices of water molecules. Natural gas hydrates, mostly methane hydrate, are found to occur in continental margins where the pressure and temperature thermodynamically satisfy the stability conditions for hydrate formation. However, it has been reported that gas hydrates may dissociate as a result of: (a) sea-level changes, e.g., Amazon Fan [4] and the western Mediterranean Sea [5]; (b) plate movement due to tectonic activities, e.g., Tokai thrust in the Eastern Nankai wedge [6]; (c) earthquake events, e.g., north-east Sakhalin margin [7]; (d) temperature or salinity change in seawater due to abnormal thermohaline circulation, e.g., Storegga slide [8]; and (e) frictional heat from the subduction zone [9]. Some evidence indicates that submarine slope failures that occurred in hydrate-bearing sedimentary basins might have been directly triggered by gas hydrate dissociation, or at least primed by the same after which some other mechanisms triggered failures (Table 1).

Methane hydrate dissociates to produce methane-saturated water and water-saturated methane gas. Owing to the huge quantity of methane gas compressed in the hydrogen bond cages of solid hydrate crystals, there is a pronounced pressure-dependent fluid expansion during dissociation, resulting in the generation of excess pore pressure (overpressure) in sediments [10,11,12,13,14,15]. The excess pore pressure reduces the effective stress of sediments and the shear strength, which results in shear failure. Because a number of oceanic gas hydrate deposits are classified as low-permeability fine-grained sediments, e.g., Ulleung basin sediments as reported in Kwon et al. [16], hydrate deposits in the Gulf of Mexico as reported in Francisca et al. [17], and Krishna–Godavari basin in India as reported in Yun et al. [18], the generation of excess pore pressure due to hydrate dissociation is known to be the most relevant process responsible for destabilizing hydrate-bearing sediments [10,11,12,13,14,15]. For example, a 1 °C increase in the temperature of hydrate-bearing sediments results in the pore fluid pressure increasing by the order of several megapascals under the no mass flux condition [15,19,20]. Therefore, any thermal change can stimulate hydrate dissociation and thus mechanically destabilize sediments [21,22,23,24].

Table 1. Case review of hydrate-related landslides.

Area	Features	Possible Primer or Triggers	References
Storegga Slide	Gradient = 10–20° Maximum thickness = 200 m Area = 34,000 km2 Volume of loss = 5500 km3	- Earthquake - Methane hydrate dissociation due to bottom water warming	Bugge et al. [25] Canals et al. [26] Vogt and Jung [8]
Buford Sea Margin (Alaska)	Width = 500 km Length = 40–50 km Thickness = 100–150 m	- Methane hydrate dissociation due to sea-level lowering	Kayen and Lee [10]
Blake Ridge	Subsidence by seafloor collapse Area = 100 km2 Volume of loss = 13 km3	- Overpressurized free gas beneath gas hydrate zone	Dillon et al. [27] Holbrook [28]
Northern East Sakhalin Margin	Gradient = 2.5° Length = 70 km Volume of mass wasting = 660 km3	- Shallow earthquakes - Methane hydrate dissociation due to sea-level lowering	Wong et al. [7]
Tokai Thrust in Eastern Nankai Wedge	Gradient = 2°–8° Failure layer concurs with BSR	- Methane hydrate dissociation due to plate uplift caused by subduction of paleo-Zenisu ridge	Cochonat et al. [6]
Amazon Fan	Gradient = less than 1° Area = 15,000 km2 Thickness = 200 m	- Methane hydrate dissociation due to sea-level lowering - Overburdening by Amazon River sediment discharge during deglacial time	Maslin et al. [4]
Balearic Basin of the Western Mediterranean Sea	Volume of mass wasting = 500 km3	- Methane hydrate dissociation due to sea-level lowering	Rothwell et al. [5]

Table 1. Case review of hydrate-related landslides.

Area	Features	Possible Primer or Triggers	References
Storegga Slide	Gradient = 10–20° Maximum thickness = 200 m Area = 34,000 km2 Volume of loss = 5500 km3	- Earthquake - Methane hydrate dissociation due to bottom water warming	Bugge et al. [25] Canals et al. [26] Vogt and Jung [8]
Buford Sea Margin (Alaska)	Width = 500 km Length = 40–50 km Thickness = 100–150 m	- Methane hydrate dissociation due to sea-level lowering	Kayen and Lee [10]
Blake Ridge	Subsidence by seafloor collapse Area = 100 km2 Volume of loss = 13 km3	- Overpressurized free gas beneath gas hydrate zone	Dillon et al. [27] Holbrook [28]
Northern East Sakhalin Margin	Gradient = 2.5° Length = 70 km Volume of mass wasting = 660 km3	- Shallow earthquakes - Methane hydrate dissociation due to sea-level lowering	Wong et al. [7]
Tokai Thrust in Eastern Nankai Wedge	Gradient = 2°–8° Failure layer concurs with BSR	- Methane hydrate dissociation due to plate uplift caused by subduction of paleo-Zenisu ridge	Cochonat et al. [6]
Amazon Fan	Gradient = less than 1° Area = 15,000 km2 Thickness = 200 m	- Methane hydrate dissociation due to sea-level lowering - Overburdening by Amazon River sediment discharge during deglacial time	Maslin et al. [4]
Balearic Basin of the Western Mediterranean Sea	Volume of mass wasting = 500 km3	- Methane hydrate dissociation due to sea-level lowering	Rothwell et al. [5]

Despite previous attempts to link hydrate dissociation to submarine slope failures (see examples in Table 1), the underlying processes for submarine slope failure—from geological changes (e.g., sea-level lowering [2,4,5,10], bottom water warming [8], and plate tectonism [9]) through hydrate dissociation and excess pore pressure generation to shear failure—remain poorly understood. Furthermore, with global warming currently leading to an increase in the temperature of the ocean-bottom water [29], referred to as bottom water warming, the potential of bottom water warming to cause hydrate dissociation and eventually to trigger slope failures has not been investigated in adequate detail.

1.2. Bottom Water Warming—Past History and Current Trend

The carbon isotope ratio (¹³C/¹²C ratio or δ¹³C) is a useful indicator to identify the source of sedimentary organic matter, specifically, whether the carbon is from marine plants or continental plants. In particular, any δ¹³C anomaly in past sedimentary formations can be interpreted as evidence of abrupt and massive input of carbon to the ocean or atmosphere from an external reservoir that is rich in either ¹²C or ¹³C [30]. The extraordinary decrease in δ¹³C, as reported by Kennett and Stott [31], coincides with the warming of the deep sea by approximately 5–7 °C ~55 million years ago at the Paleocene/Eocene transition [30,31,32,33]. This event is known as the late Paleocene thermal maximum (LPTM) or the Paleocene/Eocene thermal maximum (PETM). The anomalous decrease in δ¹³C in marine sediment deposits at the end of the Paleocene can only be explained by the massive injection of carbon into the ocean and atmosphere, possibly in the form of methane from hydrates buried in oceanic sediments [32,33,34]. The first trigger that destabilized methane hydrate has variously been considered to be heat transport from high latitudes to low latitudes due to abnormal thermohaline circulation [35], mantle-derived intrusion [36], pressure decrease due to lowering of sea levels, uplift due to plate tectonics [9], volcanism [37], or catastrophic slope failure in hydrate deposits due to seismic activity [34,38]. Irrespective of the causative factor, after methane was released to the atmosphere, it accelerated global warming and further thermally destabilized methane hydrate through conduction and advection processes [34]. Consequently, it triggered another series of slumps in the hydrate-bearing deposits and added more carbon to the ocean and atmosphere, aggravating the positive feedback. This has significant implications for the anthropogenic global warming that is being witnessed currently, because it suggests that hydrate decomposition could occur in oceanic sediments.

With anthropogenic activities already having a significant and rapid effect on the global climate, the huge volume of methane hydrate in marine sediments is a potential hazard [39]. The methane hydrate in the permafrost region of the Arctic has been reported to have already started to decompose [39,40,41,42,43,44,45]. The data of the near-surface air temperature for the past 140 years reveal a striking warming trend since the 1990s. The increase in greenhouse gases during the last 50 years has elevated the near-surface air temperature by approximately 0.5 °C. Furthermore, several climate models predict a temperature increase of ~2–6 °C by 2100 [46]. Although the warming of the ocean is delayed relative to that of air because of the large heat capacity of water [34], several studies have recently reported that the temperature of oceans has increased too [25,44,47,48,49,50]. A study showed that the ocean water up to a depth of 300 m warmed by 0.31 °C from the mid-1950s to the mid-1990s, i.e., 0.007 °C per year [47]. Another showed that ocean water up to a depth of 3000 m warmed by 0.0037 °C from 1955 to 1998, i.e., 0.00086 °C per year [29]. The updated estimate of the warming records of oceans as presented by Levitus et al. [48] indicated that the temperature of oceans worldwide for the 0–700 m layer increased by 0.18 °C from 1955 to 2010, i.e., 0.00321 °C per year, whereas that of the 0–2000 m layer warmed by 0.09 °C over the same period, i.e., 0.00161 °C per year. The oceans at high latitudes and in the Arctic regions are more sensitive to climate change than others. For instance, the sea waters above the West Spitsbergen continental margin, west of Norway’s Svalbard archipelago, have warmed by 1 °C over the last 30 years, possibly because of global warming [43]. Similarly, the Weddell Sea bottom water (1500–3500 m) was reported to have warmed by approximately 0.01 °C per year from 1970 to 1990 [50].

If bottom water warming continues in the deep oceans, it would cause the dissociation of natural gas hydrates, which in turn could trigger a slope failure in hydrate-bearing deposits. Thus, the submarine slope instability induced in hydrate-bearing deposits by bottom water warming was explored to gain a better insight into methane hydrate destabilization caused by global warming. Here, the generation of excess pore pressure due to hydrate dissociation was considered as a primer or a trigger for slope failure. In addition, analytical models for calculating the excess pore pressure caused by hydrate dissociation were presented. A hypothetical scenario, in which the temperature of the bottom water increases by 1 °C over 10,000 years and remains constant over the following 10,000 years, was examined for slope instability in deep oceanic environments. An in-situ hydrate deposit was selected as a model site, and a coupled thermal-hydraulic-mechanical (T-H-M) finite difference analysis that captured thermal conduction through sediments, thermal dissociation of gas hydrate, generation of excess pore pressure, consequential pressure diffusion, and hydrate dissociation against depressurization was performed. The effects of key parameters on the generation of excess pore pressure were explored via a parametric study.

2. Gas Hydrate Dissociation in Sediments

2.1. Thermal Dissociation of Gas Hydrate

Hydrate dissociation starts when the pressure-temperature (PT) state reaches the equilibrium boundary. In reality, the hydrate stability curve should be expressed as a broad band that reflects different boundaries for different pore sizes (see for example the works by Anderson et al. [51] and Kwon et al. [52]) rather than a single phase boundary. A single phase boundary causes instantaneous dissociation as temperature increases; in contrast, a broad margin across the PT boundary is expected to cause gradual dissociation. This is due to hydrates in smaller pores dissociating at lower temperatures than those in larger pores [51,52]. However, for simplicity, the present study assumes pure methane as a hydrate-forming gas, pure pore water, and no pore size effect on hydrate stability.

During thermal stimulation, hydrate dissociates into gas-saturated water and water-vapor-saturated gas. The released fluids raise the pore fluid pressure to the equilibrium pressure corresponding to a given temperature. The increase in pressure along the PT boundary hinders further dissociation (i.e., self-preservation; [15]). Therefore, partial dissociation during thermal stimulation is characterized by a pressure-temperature evolution along the phase boundary until all the hydrate dissociates. Due to the pronounced fluid expansion of hydrate, thermal stimulation of hydrate-bearing sediments under an undrained condition can cause a large increase in the fluid pressure (of the order of several megapascals) when volume expansion is restricted by the high skeletal stiffness B_sk of the medium. This study used the analytical model proposed by Kwon et al. [15] to calculate the pore fluid pressure and the volume fraction of methane hydrate in a unit volume of sediments, as shown in Figure 1. This study can also be referred to for the details of the analytical model for the thermal dissociation of gas hydrates. Figure 1 shows the pressure evolution and changes in hydrate saturation during thermal dissociation of methane hydrate under isochoric condition (i.e., constant volume condition with no mass flux). As seen in Figure 1, a temperature increase of 1 °C results in a pressure increase of the order of several megapascals, showing self-preservation behavior along the phase equilibrium boundary. Higher hydrate saturation results in longer self-preservation behavior and leads to higher overpressure upon complete dissociation.

Three phases—water (with dissolved gas), free gas (with water vapor), and mineral—remain in sediments when the gas hydrate is completely dissociated. Under constrained volume expansion, the thermal expansion of the phases causes an increase in the pore fluid pressure of sediments. This study disregards the thermal expansion of the phases because it causes only a minute increase in fluid pressure (e.g., less than 100 kPa with a temperature increase of 1 °C) in comparison with hydrate dissociation.

Figure 1. Thermal stimulation of methane hydrate-bearing sediments under a constant volume condition: (a) pressure evolution and (b) hydrate saturation remaining in pore spaces. Note that the initial gas saturation (S_g) is 0.1, and the thermodynamic phase boundary is plotted using the following equation given by Sloan [53]: P [kPa] = exp(38.908 − 8533.80/T [K]).

Figure 1. Thermal stimulation of methane hydrate-bearing sediments under a constant volume condition: (a) pressure evolution and (b) hydrate saturation remaining in pore spaces. Note that the initial gas saturation (S_g) is 0.1, and the thermodynamic phase boundary is plotted using the following equation given by Sloan [53]: P [kPa] = exp(38.908 − 8533.80/T [K]).

2.2. Hydrate Dissociation by Pressure Diffusion

The excess pore fluid pressure produced by hydrate dissociation subsequently generates a pressure gradient within the sediment. Overpressurized pore fluids (water and free gas) are then dissipated to lower-pressure zones. The assumption of no mass flux that was previously made for calculating the excess pore pressure is relaxed at this stage. This pressure diffusion process in the sediments can be described by using Terzaghi’s consolidation theory [54]. Considering a one-dimensional sediment column, the dissipation of the excess pore pressure can be calculated with the time scale and spatial scale by using Equation 1 [54]:

(1)

where u is the excess pore pressure; t, the time; z, the depth of the sediment column; and C_V, the coefficient of consolidation (or pressure diffusivity).

The elevated pore fluid pressure, previously calculated by using the thermal dissociation model, is applied as an initial condition for the pressure diffusion analysis. It is assumed in the flow analysis that sediments are water-saturated and have a single-phase flow of water, although there is a change in the volume fractions of phases in the pore space during hydrate dissociation (e.g., transient flow characteristics). It is presumed that water is the dominant phase that flows over gas because the scope of the present study is confined to low hydrate saturation and fine-grained sediments. When the hydrate saturation is greater than 40%, and consequently, the free gas released by hydrate dissociation percolates through pore throats, multiphase flow analyses are required [19,20].

Because of pressure diffusion, the pressure in the hydrate-bearing system becomes lower than the equilibrium pressure. When the gas hydrate stability is disrupted via depressurization, gas hydrates start dissociating. Thus, as long as hydrate remains, the system will reestablish the equilibrium state by dissociating hydrate [15,20]. The pressure will increase to the equilibrium pressure for a given temperature because of the hydrate dissociation during pressure diffusion. When the system is depressurized to a certain PT condition (P₀, T₀), the pressure P₁ (=P₀ + ΔP) increases and the hydrate-bearing sediments are stabilized again (P₁, T₀) by gas hydrate dissociation (Figure 2a). During the rebuilding of the pressure, the volume fraction of hydrate melting (S_hm) in an isochoric and isothermal condition can be derived by Equation (2):

(2)

where ρ_w and ρ_h are the density of water and hydrate, respectively; M_w and M_h, the molar mass of water and hydrate, respectively; and S_g and S_h, the volume fraction of gas and hydrate in the pore space, respectively. χ is the stoichiometric number of hydration; R, the gas constant; P, the pressure; and T, the absolute temperature. The gas solubility k_H [mol·m⁻³·atm⁻¹] in water at temperature T is calculated by Henry’s law [Equation (3)]:

(3)

where ΔH_solution [J·mol⁻¹] is the enthalpy of the solution. k⁰_H,298.15K [mol·L⁻¹·atm⁻¹] is Henry’s law constant, which represents the solubility under standard conditions. Methane has very low solubility even in pure water, and its solubility decreases further in saline water. Herein, because of the lack of data about the solubility of methane in saline water, the solubility in pure water was used throughout this study (i.e., ΔH_solution = 13,000 J·mol⁻¹, k⁰_H,298.15K = 1.4 × 10⁻³ mol·m⁻³·atm⁻¹ from [55]).

As the initial saturation of free methane gas (S_g) increases, the volumetric quantity of methane hydrate (S_hm) that is required to dissociate to re-establish the equilibrium increases. When S_g is 10% and the initial pressure is 7.3 MPa, S_hm of 10% is adequate for rebuilding the pressure by approximately 10 MPa, as shown in Figure 2(a). The pressure recovery by the dissociation of methane hydrate (S_hm) can be calculated for different initial free gas saturations (S_g) in pore spaces, as shown in Figure 2(b). In the following analysis, the pressure diffusion process and the consequential hydrate dissociation are evaluated by considering the amount of dissociated hydrate and remaining hydrate in the dissociation layer.

Figure 2. Dissociation of methane hydrate against depressurization: (a) pressure generation by hydrate dissociation and (b) the volume fraction of hydrate dissociated to reestablish the equilibrium pressure. Note that the thermodynamic phase boundary is plotted using the following equation given by Sloan [53]: P [kPa] = exp (38.908 − 8533.80/T [K]).

Figure 2. Dissociation of methane hydrate against depressurization: (a) pressure generation by hydrate dissociation and (b) the volume fraction of hydrate dissociated to reestablish the equilibrium pressure. Note that the thermodynamic phase boundary is plotted using the following equation given by Sloan [53]: P [kPa] = exp (38.908 − 8533.80/T [K]).

3. Numerical Analysis—Effect of Bottom Water Warming

3.1. Scope of Numerical Analysis

Destabilization of marine gas hydrate-bearing sediments subjected to thermal loading was explored through a one-dimensional numerical analysis. The increase in greenhouse gases over the last 50 years (1955–1998) has elevated the temperature of deep oceans up to a depth of 3000 m by 0.0037 °C, or by 0.00086 °C per year [29]. If this trend continues, resulting in the bottom water of deep oceans being warmed by 1 °C over 10⁴ years and kept constant for the following 10⁴ years, the following question arises: would it be large enough to dissociate natural gas hydrates and eventually trigger slope failure in hydrate-bearing deposits? To examine this problem, a coupled T-H-M algorithm that can capture the underlying physical processes, including thermal conduction through sediments, thermal dissociation of gas hydrate, generation of excess pore pressure, consequential pressure diffusion, and hydrate dissociation against depressurization, was developed based on one-dimensional finite difference modeling. In particular, thermal dissociation of hydrate and dissociation against pressure depressurization were considered by incorporating the analytical models described above.

Gas hydrates are thermodynamically most sensitive at the base of the hydrate occurrence zone (hereafter, BHOZ: it can also be detected as the bottom simulating reflector, BSR), where the pressure and temperature conditions become unstable more easily because of thermal or pressure change, than in the region above the BHOZ. Thus, the BHOZ was considered as a potential hydrate dissociation layer. The depth of the BHOZ varies with water depth but is mostly within several hundred meters below the seafloor.

3.2. Algorithm for Sequentially Coupled T-H-M Analysis

Thermal conduction, hydrate dissociation, pressure diffusion response were sequentially solved by updating the parameters for the next step in the numerical simulation. An algorithm for the sequentially coupled T-H-M analysis is described here:

(1) Thermal Conduction Analysis. Heat flows into the hydrate-bearing sediments from the warmed bottom water at the seafloor, and it can be transported by convection (or advection) and conduction. In marine sediments, conduction is the primary mechanism, and not convection, because of the low permeability of gas hydrate-bearing sediments. In the numerical modeling, the thermal diffusivity of the media was considered to be constant. The temperature increase at the region of the BHOZ (450 mbsf) was calculated through a one-dimensional finite difference analysis using Fourier’s heat conduction law [Equation (4)]:

(4)

where T is the temperature; κ, the thermal diffusivity of the medium (m²·s⁻¹), which indicates how a thermal front moves through the medium; k_t, the thermal conductivity (W m⁻¹·K⁻¹); C_p, the specific heat (J·kg⁻¹·K⁻¹); and ρ_sed, the density of the sediment.

(2) Hydrate Dissociation Analysis. The excess pore fluid pressure induced by thermal change and hydrate dissociation at the bottom of a gas hydrate occurrence region at 450 mbsf, which was one grid point in the mesh, was calculated using the aforementioned analytical model. If any hydrate existed (S_h > 0) and the temperature exceeded the hydrate stability boundary (T > T_eq), the pore pressure P of the grid point increased to reestablish a new equilibrium pressure P_eq The remaining hydrate fraction S_h was calculated with consideration of the volumetric strains and pressure increase. If the conditions for hydrate dissociation (i.e., hydrate presence and thermodynamical stability) were not satisfied, the pore pressure and hydrate fraction remained unchanged.

(3) Pressure Diffusion Analysis. A pressure diffusion analysis was performed using Equation 1 in order to obtain the distribution of pore fluid pressure after the pressure diffusion in the medium. The excess pore fluid pressure (ΔP = P_eq − P₀ calculated from the hydrate dissociation analysis) was used as an initial condition for the pressure diffusion analysis. Then, if the pressure calculated after the pressure diffusion was less than the hydrate equilibrium pressure at a given temperature (P < P_eq), hydrate dissociation by depressurization (pressure diffusion) was calculated by using Equation (2).

(4) Calculation of Excess Pore Pressure and Factor of Safety (FS). The final excess pore pressure (u_ex) was calculated by subtracting the initial pore fluid pressure from the current pore fluid pressure after the pressure diffusion analysis. The factor of safety (FS) was then estimated using the excess pore pressure and Equation 5. To calculate FS, it was assumed that the hydrate-bearing reservoir was an infinite slope with the BHOZ parallel to the seafloor, implying that the potential failure plane was also parallel to the seafloor and excess water pressure was developed at the failure zone [10]. The FS on an infinite hydrate-bearing slope was derived based on force equilibrium, as shown in Equation (5) [10]:

$$FS=\frac{{\text{F}}_{\text{resisting}}}{{\text{F}}_{\text{driving}}}=\frac{\left[\left({\rho }_{sat}-{\rho }_w\right)gz\cos \beta -\frac{u_{ex}}{\cos \beta }\right]\tan {\varphi }_{cs}}{\left({\rho }_{sat}-{\rho }_w\right)gz\sin \beta }$$

(5)

where ρ_sat and ρ_w are the density of seawater-saturated sediment and the density of seawater, respectively. ϕ_cs is the friction angle of sediments at a critical state; β, the slope angle; u_ex; the excess pore pressure induced by hydrate dissociation; z, the depth of the failure plane; and g, the gravitational acceleration. After updating changes in properties including hydrate saturation, water saturation, gas saturation, pressure and temperature, and the time increment, this procedure was repeated for a time period of interest.

3.3. Geological Conditions and Input Properties

The Blake Ridge 997 site (ODP Leg 164) was selected as a model site for hydrate stability and hydrate occurrence [56]. The geological conditions and sediment properties for the thermal, hydraulic, and mechanical analyses are listed in Table 2. The hydrate stability zone was determined on the basis of the water column depth (i.e., 2750 m) and the geothermal conditions of formation (i.e., seafloor temperature of 3 °C and geothermal gradient of 36.9 °C/km; see Table 2). It was considered that methane hydrate occurred from the seafloor to the BHOZ (450 mbsf) with a pore saturation of 10%, as shown in Figure 3. It was assumed that the hydrate-forming gas was pure methane and that the pore water was an aqueous solution containing 3.5 wt % of NaCl. A free gas zone has been observed to exist beneath the BHOZ, and the free gas occupies ~10% of the pore spaces [57]. Free gas with a pore saturation of 10% was assumed to be present below 450 mbsf. When a free gas zone exists beneath the hydrate occurrence zone, the sediments at the bottom of the hydrate layer may overlap with the free gas zone and could be considered to be in a three-phase equilibrium (i.e., hydrate, water, free gas), as in the case of Blake Ridge [56,57,58,59]. Thus, the BHOZ located at 450 mbsf was assumed to contain both methane hydrate (pore saturation of 10%) and free methane gas (pore saturation of 10%), establishing a three-phase equilibrium. The three-phase equilibrium was simply considered as a thermodynamic border irrespective of whether the methane hydrate was stable or unstable. In contrast, the hydrates above the BHOZ, which were in two-phase equilibrium, would hardly undergo decomposition by a temperature increase of 1 °C.

The mechanical and thermal properties of the hydrate-bearing sediments used here were obtained from the published data of natural hydrate-bearing sediments (e.g., Paull et al. [56]; Yun et al. [60]; Mesri [61,62]; Henninges et al. [63]; Santamarina and Ruppel [64]; Cortes et al. [65]; Waite et al. [66]). This study considered uniform fine-grained sediment as a host medium for simplification. The sediment was assumed to have a thermal diffusivity of 1 × 10⁻⁶ m²·s⁻¹, which was within the range of the published data [63,64,65,66,67]. Assuming homogeneous distribution of hydrate in sediments, the presence of hydrate affects permeability by reducing the pore size and/or blocking pore. The permeability in hydrate-bearing sediments is presumed to be relatively lower than that in gassy sediments below the BHOZ [68,69,70]. Due to a lack of direct measured data on the coefficient of consolidation or permeability of hydrate-bearing sediments, an arbitrary value of low-permeability clay sediments (C_V = 10⁻⁸ m²/s) was used whereas for the gassy sediments below the BHOZ, the general value of silty-clay sediments (C_V = 10⁻⁷ m²/s) was chosen. A typical value for the friction angle of 30° was taken in this study. The contribution of the cohesion intercept c' to the shear strength (τ = c' + σ'tan ϕ_cs) is assumed to be insignificant relative to the one by the effective stress (σ') and internal friction angle (ϕ_cs).

Figure 3. Geological conditions of the reconstructed model site. Note that the BHOZ represents the base of the hydrate occurrence zone and the methane hydrate stability curve is plotted based on the following equation given by Peltzer and Brewer [71]: 1/T [K] = 3.83 × 10⁻³ − 4.09 × 10⁻⁴ log₁₀P [MPa] + 8.64 × 10⁻⁵ (log₁₀P [MPa])².

Figure 3. Geological conditions of the reconstructed model site. Note that the BHOZ represents the base of the hydrate occurrence zone and the methane hydrate stability curve is plotted based on the following equation given by Peltzer and Brewer [71]: 1/T [K] = 3.83 × 10⁻³ − 4.09 × 10⁻⁴ log₁₀P [MPa] + 8.64 × 10⁻⁵ (log₁₀P [MPa])².

Table 2. Geological conditions and properties of hydrate-bearing sediments used in this study.

Property	Values
Hydrostatic pressure at the seafloor [MPa] a,b	27.8
Temperature at the seafloor [°C] a,c	3.3
Geothermal gradient [°C·km−1] a,c	36.9
Depth of the base of the hydrate occurrence zone (BHOZ) [mbsf] a,b	450
Pore volume fraction of methane hydrate in hydrate layer [%] a,b	10
Pore volume fraction of methane gas in free gas layer [%] a,d	10
Thermal diffusivity of hydrate-bearing sediments [m2·s−1] e	10−6
Coefficient of consolidation of hydrate-bearing sediments [m2·s−1]	10−8
Coefficient of consolidation of gassy sediments [m2·s−1]	10−7
Mass density of water-saturated sediments [g·cm−3]	1.8
Friction angle of sediments at critical state [°]	30

Data from: ^a Paull et al. [56]; ^b Henry et al. [72]; ^c Ruppel [73]; ^d Dickens et al. [57]; ^e Briau and Chaouch [74].

Table 2. Geological conditions and properties of hydrate-bearing sediments used in this study.

Property	Values
Hydrostatic pressure at the seafloor [MPa] a,b	27.8
Temperature at the seafloor [°C] a,c	3.3
Geothermal gradient [°C·km−1] a,c	36.9
Depth of the base of the hydrate occurrence zone (BHOZ) [mbsf] a,b	450
Pore volume fraction of methane hydrate in hydrate layer [%] a,b	10
Pore volume fraction of methane gas in free gas layer [%] a,d	10
Thermal diffusivity of hydrate-bearing sediments [m2·s−1] e	10−6
Coefficient of consolidation of hydrate-bearing sediments [m2·s−1]	10−8
Coefficient of consolidation of gassy sediments [m2·s−1]	10−7
Mass density of water-saturated sediments [g·cm−3]	1.8
Friction angle of sediments at critical state [°]	30

Data from: ^a Paull et al. [56]; ^b Henry et al. [72]; ^c Ruppel [73]; ^d Dickens et al. [57]; ^e Briau and Chaouch [74].

A 1000-m-thick sediment deposit was modeled. The seafloor, which acted as an isothermal boundary (heat source) with a temperature rise of 0.0001 °C per year, was placed at the top boundary. Homogeneous pore-filling methane hydrate was assumed to exist from the seafloor to 450 mbsf with a pore saturation of 10%. The time period covered in each simulation was 20,000 years with a finite time increment of 1 year. During the numerical simulation period, the geothermal temperature, excessive pore fluid pressure, hydrate volume fraction, and FS at the BHOZ (450 mbsf) were monitored.

4. Analysis Results

4.1. Thermal Responses at the Base of Hydrate Occurrence Zone (BHOZ)

When the temperature of the bottom water increased by 1 °C over 10⁴ years and stayed constant over the following 10⁴ years, the thermal response at a depth of 450 mbsf was calculated via a one-dimensional finite difference analysis, as shown in Figure 4. It was assumed that heat conduction was the primary mechanism for heat transfer in the sediment and that thermal conductivity and diffusivity were kept constant. Because of the thermal diffusion from the seafloor, the temperature of the BHOZ at 450 mbsf increased gradually, showing a time delay of more than 5000 years. For shallow hydrate deposits, in which the BHOZs were located at shallower depths of 100–200 m, e.g., the Cascadia continental margin, Gulf of Mexico, and the Arctic regions in permafrost environments, the temperature increase at the BHOZ is expected to occur more rapidly.

Figure 4. Thermal response at the BHOZ due to the increase in the seafloor temperature.

Figure 4. Thermal response at the BHOZ due to the increase in the seafloor temperature.

The thermal diffusivity of hydrate-bearing sediments is a complex feature of sediment mineralogy, porosity, and pore fluid composition (i.e., water, gas hydrate, and free gas). Any physical or chemical reaction caused by a temperature change changes the thermal properties of hydrate-bearing sediments. For instance, an increase in the temperature of hydrate-bearing sediments during thermal diffusion results in gradual dissolution of hydrate, which leads to a change in the volumetric composition of the pore fluids. This is because the solubility of methane in pore water increases with increasing temperature in the presence of hydrate at two-phase equilibrium [75,76,77]. In addition, a temperature increase that is large enough to dissociate hydrate would significantly alter the thermal properties of the media. However, in the numerical modeling, it was considered that the thermal conductivity and thermal diffusivity of the media remained constant during the heat transfer process.

Water and methane hydrate have similar thermal conductivity values (0.56 W·m⁻¹·K⁻¹ for water at 0 °C and ~0.6 W·m⁻¹·K⁻¹ for methane hydrate; [78]). Thus, the hydrate presence does not increase the thermal conductivity of the sediments significantly (e.g., there is an increase of less than 1% for a hydrate saturation of 10%) [78]. By using an effective medium model, it has been shown that an increase in methane hydrate saturation from 0% to 10% leads to a decrease of ~5% in the specific heat and an increase of ~5% in the thermal diffusivity [78]. In contrast, hydrate dissociation releases free gas, which is likely to cause a decrease in the density, specific heat, and thermal conductivity of sediment. However, the effect of free gas generation on the thermal diffusivity of hydrate-bearing sediments is not well understood, although it is presumed that reduction of density and specific heat may compensate the thermal conductivity decrease [i.e., κ = k_t/(C_p·ρ_sed)] during free gas generation. When hydrate saturation is fairly high, e.g., more than 10%, the free gas saturation resulting from the complete dissociation of the hydrate would be significant, and it would alter the thermal diffusivity of a medium. Because hydrate saturation was low (i.e., 10%) and thermal propagation was a fairly slow process (i.e., over 10,000 years), however, it was justifiable to assume a constant thermal conductivity (and thermal diffusivity) and disregard the heat of hydrate dissociation in this study. However, the assumption that hydrate dissociation would not significantly alter the thermal diffusion rate is a reasonable approximation only under limited conditions of low hydrate saturation less than 10%.

4.2. Generation of Excess Pore Pressure and Pressure Diffusion

Using the thermal response at the BHOZ as an isothermal input (Figure 4), the excess pore pressure caused by the elevated temperature was calculated by using the analytical model presented by Kwon et al. [15]. Figure 5 shows the changes in excess pore pressure and hydrate saturation remaining and dissociated with time when no mass flux and no volume expansion are allowed. The maximum pressure was ~6.5 MPa and the volume fraction of hydrate dissociated for re-establishing equilibrium was ~0.07 (Figure 5).

On the other hand, the generation of excess pore pressure also creates flows of pore fluids to surrounding regions, dissipating the overpressure, which is referred to as pressure diffusion. The pressure diffusion process continues until all the excess pore pressure is dissipated. However, as long as methane hydrate remains in the system, gas hydrate dissociates against depressurization during the pressure diffusion, recovering the equilibrium pressure. Equation (2) previously shown in Section 2 can capture this physical process of hydrate dissociation against depressurization during pressure diffusion. Therefore, two different physical processes of hydrate dissociation—(1) thermal dissociation of hydrate during temperature increase and (2) hydrate dissociation against depressurization due to pressure diffusion (depressurization)—need to be considered. By performing the sequential analysis that coupled thermal diffusion, thermal dissociation of hydrate, pressure diffusion, and hydrate dissociation against depressurization, the excess pore pressure and the volume fraction of gas hydrate remaining at the BHOZ were calculated, as shown in Figure 6. As gas hydrate dissociated and excess pore pressure was accumulated, the excess pore pressure at the dissociation zone reached the maximum value (~3.26 MPa) after 8827 years. Subsequently, as all the hydrate was completely dissociated [Figure 6(b)], no further hydrate dissociation occurred, and the excess pore pressure at the BHOZ started to decline because of pressure diffusion [Figure 6(a)].

Figure 5. Isochoric destabilization of methane hydrate at BHOZ: (a) excess pore pressure developed and (b) volume fractions of remaining and dissociated hydrates corresponding to temperature increase.

Figure 5. Isochoric destabilization of methane hydrate at BHOZ: (a) excess pore pressure developed and (b) volume fractions of remaining and dissociated hydrates corresponding to temperature increase.

Figure 6. (a) Excess pore pressure developed, (b) volume fractions of remaining and dissociated gas hydrates and volume fractions of dissociated gas hydrates against temperature increase and pressure diffusion at the BHOZ.

Figure 6. (a) Excess pore pressure developed, (b) volume fractions of remaining and dissociated gas hydrates and volume fractions of dissociated gas hydrates against temperature increase and pressure diffusion at the BHOZ.

Time-lapse one-dimensional spatial distributions of the excess pore pressure over the hydrate occurrence zone and free gas zone are shown in Figure 7. As the dissociation occurred at the BHOZ, the excess pore pressure was diffused centered on the BHOZ at 450 mbsf. A downward flow to the free gas zone existing beneath the BHOZ was also considered. The pressure continued increasing because of hydrate dissociation [Figure 6(a)], whereas the overpressurized fluids in the dissociation layer flowed upward and downward [Figure 7(a)]. Because hydrate-bearing sediments had lower permeability than did underlying gassy sediments, more overpressurized pore fluid flowed downward, resulting in an overpressurized free gas zone. The diffusion calculation of this analysis suggests a possible mechanism for the occurrence of the overpressurized free gas zone, where sediment failure may occur if the free gas formation is critically pressurized (e.g., Hornbach et al. [79]). While the time scales for thermal diffusion and pressure diffusion were of the order of hundreds of years in this study, the time scales for hydrate dissociation were of the order of hours. Thus, the rates of thermal diffusion and pressure diffusion governed the rate of the excess pore pressure generation rather than the hydrate dissociation kinetics.

Figure 7. One-dimensional spatial distribution of excess pore pressure due to hydrate dissociation.

Figure 7. One-dimensional spatial distribution of excess pore pressure due to hydrate dissociation.

4.3. Factor of Safety (FS)

The stability of an infinite hydrate-bearing slope undergoing gas hydrate dissociation was evaluated using the factor of safety (FS). FS is the ratio of the resisting shear strength to the driving shear stress along an assumed planar failure surface. Thus, when FS < 1, the slope will fail, and when FS > 1, the slope is stable. The shear strength of the sediment at a critical state (τ) can be determined using the Mohr-Coulomb relationship τ = c' + σ'_n tanϕ_cs, where c' is the sediment cohesion (assumed to be insignificant in this study); σ'_n, the effective stress normal to the plane of shear failure; and ϕ_cs, the friction angle of the sediment at a critical state. The effective stress σ'_n is closely related to the pore pressure. As the temperature increase causes hydrate dissociation in sediments and elevates the pore pressure, it reduces the effective stress σ'_n and shear strength τ, and in turn, the resisting force [i.e., F_resisting, numerator in Equation (5)]. The pore pressure can increase until the downward component of the gravitational force [i.e., driving force F_driving, denominator in Equation (5)] becomes greater than the resisting force. At this point, the FS becomes less than one, and failure occurs.

The FS for an infinite hydrate-bearing slope undergoing hydrate dissociation was calculated using the excess pore pressure distributions (u_ex) induced by hydrate dissociation (Figure 7) and Equation (5). It was assumed that methane hydrate occurred in an infinite slope and that the BHOZ of the hydrate deposit was parallel to the seafloor. This implies that a potential failure plane would be parallel to the seafloor [10]. The FS of the hydrate-bearing slope with an inclination of 3° was calculated based on the parameters listed in Table 2, as shown in Figure 8. The FS became less than unity after 8827 years when the maximum excess pore pressure was generated. This suggests that the excess pore pressure generated by hydrate dissociation could be large enough to trigger a slope failure when the inclination is 3°. The probability of a slope failure increased with the slope inclination, as shown in Figure 9.

Figure 8. Factor of safety of an infinite hydrate-bearing slope having an inclination of 3°.

Figure 8. Factor of safety of an infinite hydrate-bearing slope having an inclination of 3°.

Figure 9. Factor of safety of an infinite hydrate-bearing slope with varying inclinations.

Figure 9. Factor of safety of an infinite hydrate-bearing slope with varying inclinations.

The hydrate dissociation mechanism proposed here could be both a primer and a trigger for slope failure. In some cases, the excess pore pressure generated by hydrate dissociation could be significant enough to trigger a slope failure with no other external trigger. When the sediment is less permeable, hydrate dissociation against pressure diffusion will be less. Thus, the maximum excess pore pressure will be higher than the value calculated in Figure 6. In addition, when the hydrate saturation is higher, self-preservation persists for a longer duration, and thus, the excess pore pressure will build further. In all cases, the risk of failure increases. The effect of the relevant parameters on the slope stability is discussed in the following section.

5. Discussion and Implications

5.1. Parametric Study

With the simulation results shown in Section 4 as a reference case (REF), nine cases in which parameters including thermal diffusivity, pressure diffusivity (coefficient of consolidation), and hydrate and gas saturations (Table 3) were varied were explored. The results are presented in terms of excess pore pressure.

Table 3. Summary of Simulated Cases.

Case	Thermal diffusivity [m−2·s−1]	Coefficient of consolidation above BHOZ [m−2·s−1]	Hydrate saturation [%]	Gas Saturation [%]
REF	10−6	10−8	10	10
TH1	10−5	10−8	10	10
TH2	5 × 10−7	10−8	10	10
TH3	10−7	10−8	10	10
HP1	10−6	10−9	10	10
HP2	10−6	10−7	10	10
HP3	10−6	10−6	10	10
SH1	10−6	10−8	5	10
SH2	10−6	10−8	20	10
SH3	10−6	10−8	10	30

Table 3. Summary of Simulated Cases.

Case	Thermal diffusivity [m−2·s−1]	Coefficient of consolidation above BHOZ [m−2·s−1]	Hydrate saturation [%]	Gas Saturation [%]
REF	10−6	10−8	10	10
TH1	10−5	10−8	10	10
TH2	5 × 10−7	10−8	10	10
TH3	10−7	10−8	10	10
HP1	10−6	10−9	10	10
HP2	10−6	10−7	10	10
HP3	10−6	10−6	10	10
SH1	10−6	10−8	5	10
SH2	10−6	10−8	20	10
SH3	10−6	10−8	10	30

(1) Thermal Diffusivity. The rate of thermal front propagation is proportional to the thermal diffusivity of a medium. In addition to the REF case, three cases with different thermal diffusivities (10⁻⁵ m²·s⁻¹ for TH1, 10⁻⁶ m²·s⁻¹ for REF, 5 × 10⁻⁶ m²·s⁻¹ for TH2, and 10⁻⁷ m²·s⁻¹ for TH3) were simulated. Figure 10a and Figure 10b show the thermal response and the generation of the excess pore pressure at the BHOZ for various levels of thermal diffusivity. Higher thermal diffusivity appeared to cause a faster temperature response (larger temperature increase for a given time), which resulted in the faster generation of a larger excess pore pressure. These results suggest that increasing the thermal diffusivity increases the susceptibility of a slope for failure because of the faster thermal responses to bottom water warming and larger excess pore pressure.

(2) Pressure Diffusivity. The pressure diffusion rate, which is determined by the sediment permeability, plays an important role in determining the generation of excess pore pressure. Figure 10c shows the excess pore pressures in relation to various levels of the pressure diffusivity (coefficient of consolidation) of hydrate-bearing sediments (i.e., above the BHOZ). These levels were expressed in order of pressure diffusivity as follows: 10⁻⁹ m²·s⁻¹ for HP1 (stiff clays), 10⁻⁸ m²·s⁻¹ for REF (soft clays), 10⁻⁷ m²·s⁻¹ for HP2 (silt-to-clay sediments), and 10⁻⁶ m²·s⁻¹ for HP3 (fine sands). While the largest excess pore pressure was found for HP1, the dissociation of methane hydrate generated minimal excess pore pressure when the pressure diffusivity was 10⁻⁶ m²·s⁻¹, which corresponds to a typical value for fine sands. Higher pressure diffusivity (higher permeability) allowed the fluid to flow more freely, resulting in faster dissipation of excess pore pressure. As a result, the hydrate dissociated faster to maintain the equilibrium pressure. Consequently, a slope was found to be less susceptible to failure when the grain size coarsened and the sediment was more permeable.

Figure 10. Results of the cases studied: (a) thermal responses at the BHOZ for the cases of TH1, TH3, and TH3; (b) excess pore pressure at the BHOZ for the cases of TH1, TH2, and TH3; (c) excess pore pressure at the BHOZ for the cases of HP1, HP2, and HP3; (d) excess pore pressure at the BHOZ for the cases of SH1, SH2, and SH3.

Figure 10. Results of the cases studied: (a) thermal responses at the BHOZ for the cases of TH1, TH3, and TH3; (b) excess pore pressure at the BHOZ for the cases of TH1, TH2, and TH3; (c) excess pore pressure at the BHOZ for the cases of HP1, HP2, and HP3; (d) excess pore pressure at the BHOZ for the cases of SH1, SH2, and SH3.

(3) Hydrate and Free Gas Saturations. Hydrate and gas saturations are critical to the sediment instability. As hydrate saturation increased, the excess pore fluid pressure was maintained longer at the equilibrium pressure during dissociation, despite the pressure diffusion process. Thus, higher hydrate saturations yielded higher pressure and a larger overpressurized zone. Consequently, this phenomenon increased the risk of failure, as shown in Figure 10(d). In contrast, the presence of free gas in the pores increased stability by acting as a cushion against the volume expansion of pore fluids during hydrate dissociation (see SH3 in Figure 10d; S_g = 30%).

5.2. Thermal Destabilization of Hydrate-Bearing Sediments

A small increase of even 1 °C in the bottom water temperature in the deep ocean could lead to slope failure in gas hydrate-bearing deposits in some cases, even though the heat transfer to the BHOZ requires time. Considering that heat diffusion by bottom water warming starts from the seafloor, a failure would be primed from the upper part of a continental margin, where the BHOZ occurs at a shallower depth. This would cause catastrophic or progressive slope failure. When the thermal diffusivity of sediments was 10⁻⁶ m/s², the time delay from initiation of bottom water warming to slope failure was estimated to be of the order of 10³–10⁴ years. When the thermal diffusivity was 10⁻⁷ m/s², the time delay would be of the order of 10⁴–10⁵ years. This implies that the bottom water warming during the LPTM, which lasted approximately 10⁴ years [30,31,32,33,34], might have been significant enough to dissociate hydrate and to prime or trigger slope failure in oceanic hydrate-bearing sediments. The hypothesis that major slides at the Norwegian-Barents margin during the Holocene era were primed or triggered by hydrate dissociation caused by bottom water warming and not by sea-level lowering [8] can be supported by the analysis results. Massive releases of methane hydrate into the ocean and atmosphere accompanying negative excursions of carbon isotope (¹³C) have occurred in the past (Late Paleocene, 55 ma; Early Cretaceous, 117 ma; Early Toarcian, 183 ma; [38]). If bottom water warming stimulates a naturally occurring hydrate-bearing slope, the area of the landslide (volume of detachment) can be determined, and consequently the amount of released methane can be derived. Similarly, integration of the temperature records of bottom water, thermal diffusivity of sediments, isotope records, and geophysical data of a submarine slope will allow more reliable interpretation of historical events.

5.3. Shallow Methane Hydrate Deposits

In addition, shallow methane hydrate deposits, in which the BHOZ is located at shallow depth, e.g., in the Artic regions, Gulf of Mexico, and Cascadia Margin, were explored. Under the circumstances where water depth was less than 1000 m and the BHOZ occurred at a shallow depth of less than 200 mbsf, thermal diffusion was faster, and thus, excess pore pressure was generated more rapidly. This increased the susceptibility of a slope to failure. It was found that gas hydrates remained during pressure diffusion over 2 × 10⁴ years in shallow hydrate deposits; therefore, excess pore pressure generation was more dependent on the heat diffusion rate than on the pressure diffusion rate. Permeable sediments can still be liquefied, and shallow debris flows can occur, because the upward fluid flow may compensate the effective stress in the near seafloor, as pointed out by Xu and Germanovich [12].

5.4. Implications for Sea-Level Lowering

Besides bottom water warming, sea-level lowering is commonly hypothesized as a trigger that destabilizes methane hydrate and causes slope instabilities [2,4,5,7,10]. To date, an assumption that sea-level change is proportionally equivalent to the change in hydrostatic pressure (i.e., pressure = ρ_w·g·h) is overly simplified, ignoring thermal expansion or shrinkage. Sea-level change results from two main mechanisms: (a) volume change due to seawater density change in response to temperature and salinity variations and (b) mass change due to exchange of water with the atmosphere and continents, including glaciers and ice sheets [80]. Sea-level rise due to global warming in the 21st century is expected to be mainly due to the thermal expansion of ocean water; in fact, it has been reported that 70% of the rise is due to thermal expansion [42,80]. While correction for salinity variation and thermal expansion is necessary in order to accurately access the hydrostatic pressure change corresponding to the eustatic sea-level change, no such correction has been made to date.

The maximum sea-level variation in the past is thought to be ~100–200 m based on the literature [81,82,83,84]. Thus, when considering the change induced in seawater density by salinity or heat, the actual pressure drop during sea-level change is presumed to be less than 1 MPa. Thus, the maximum excess pore pressure resulting from gas hydrate dissociation would be of the order of several hundred kilopascals. This implies that sea-level lowering can prime or trigger slope failure under limited conditions with high hydrate saturation, low permeability, or shallow methane hydrate deposits.

6. Conclusions

The submarine slope failure of hydrate-bearing deposits induced by bottom water warming, which has been accelerating since the 1950s because of global warming, was explored in the presented study. The analysis method and results presented here advance the understanding of hydrate-related slope instability problems. A coupled T-H-M finite difference analysis was performed to capture the underlying physical processes including thermal conduction through sediments, thermal dissociation of gas hydrate, excess pore pressure generation, pressure diffusion, and hydrate dissociation against depressurization. In particular, the analytical formulations for calculating excess pore pressure and hydrate quantity consumed during dissociation were incorporated in the coupled analyses.

A hypothetical scenario, in which the temperature of the bottom water increased by 1 °C over 10,000 years and remained constant over the following 10,000 years, was examined for slope failure in the deep oceanic environment. The BHOZ was considered as a potential shear plane of slope failure and the initial weakened zone before catastrophic shear band growth. A temperature rise of 1 °C over 10,000 years at the seafloor was observed to generate huge excess pore pressure at the BHOZ, which was large enough to reduce the stability of a slope. The case study results show that a slope becomes more susceptible to failure with increases in thermal diffusivity and hydrate saturation and decreases in pressure diffusivity, gas saturation, and water depth. The scenario of bottom water warming, as shown in the present study, can be further explored to gain a better understanding of the past methane hydrate destabilization events on earth.