Spatiotemporal destabilization modes of upper continental slopes undergoing hydrate dissociation

Journal Pre-proof Spatiotemporal destabilization modes of upper continental slopes undergoing hydrate dissociation

Fang Liu, Lin Tan, Giovanni Crosta, Yu Huang PII:

S0013-7952(19)30635-0

DOI:

https://doi.org/10.1016/j.enggeo.2019.105286

Reference:

ENGEO 105286

To appear in:

Engineering Geology

Received date:

6 April 2019

Revised date:

3 September 2019

Accepted date:

4 September 2019

Please cite this article as: F. Liu, L. Tan, G. Crosta, et al., Spatiotemporal destabilization modes of upper continental slopes undergoing hydrate dissociation, Engineering Geology (2019), https://doi.org/10.1016/j.enggeo.2019.105286

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier.

Journal Pre-proof Spatiotemporal destabilization modes of upper continental slopes undergoing hydrate dissociation

First author: Fang Liu (corresponding author) Affiliation: State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, China; Key Laboratory of Geotechnical and Underground Engineering (Tongji University), Ministry of Education, China Address: 1239 Siping Road, Shanghai 200092, China

f

Tel: 0086-21-65983142

oo

Email: [email protected]

pr

Second author:Lin Tan

e-

Affiliation: State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, China; Key Laboratory of Geotechnical and Underground Engineering (Tongji

Pr

University), Ministry of Education, China

Address: 1239 Siping Road, Shanghai 200092, China

rn

al

Email: [email protected]

Third author: Giovanni Crosta

Jo u

Affiliation: Department of Earth and Environmental Sciences, University Milano-Bicocca, Italy

Address: Piazza della Scienza, 4 - 20126 Milan, Italy Email: [email protected]

Fourth author: Yu Huang Affiliation: Key Laboratory of Geotechnical and Underground Engineering (Tongji Un iversity), Ministry of Education, China Address: 1239 Siping Road, Shanghai 200092, China Email: [email protected]

-1-

Journal Pre-proof

Abstract: Gas hydrates dissociation could induce or trigger submarine landslides, especially in upper continental slopes where hydrates are vulnerable to natural and artificial perturbations. This work investigates destabilization mechanisms of an upper continental slope undergoing hydrate dissociation and identifies spatiotemporal failure modes influenced by characteristics of the overburden

above

the

hydrate-bearing

layer

(i.e.

the

hydrate

reservoir).

A

oo

f

Thermo-Hydro-Chemical coupled numerical model of transient pore pressure induced by hydrate

pr

dissociation is coupled with the limit equilibrium slope analysis method to study the spatiotemporal evolution of the potential sliding plane and to calculate the corresponding factor of

e-

safety. The results suggest that overpressure generated by the liberated fluid from hydrate

Pr

dissociation is the primary reason for instability in a gentle marine slope. The study identif ies

al

three sliding modes, namely co-melting non-interface sliding, co-melting interface sliding, and

rn

post-melting non-interface sliding, depending on the overburden’s characteristics, including

Jo u

overburden thickness, permeability, and cohesion. Co-melting non-interface sliding takes place during hydrate dissociation if the hydrate reservoir underlies a thin, pervious and low-cohesion overburden cover. For less permeable and more cohesive overburdens, the potential sliding plane is deeper and co-melting interface sliding could be triggered due to overpressure developed at the reservoir-overburden interface. If the hydrate reservoir is covered by a thick, low-permeability and slightly cohesive overburden, post-melting non-interface sliding could occur after the hydrates are completely dissociated. This failure is delayed, because the gas/water trapped at the interface during hydrate dissociation is insufficient to trigger instability due to very high overburden stresses. However, as the gas migrates upwards over time and encounters a weak zone in the

-2-

Journal Pre-proof

overburden deposits, failure could happen within the overburden deposits even after hydrate dissociation stops. The findings help to improve our fundamental understanding about the destabilization mechanism and failure modes of the continental slopes undergoing hydrate dissociation, and to delineate the vulnerable configurations of the slopes.

Keywords: gas hydrate; submarine landslide; upper continental slope; THC coupled analysis;

Jo u

rn

al

Pr

e-

pr

oo

f

pore pressure

-3-

Journal Pre-proof

1 Introduction Gas hydrates are ice-like crystalline compounds with natural gas (predominantly methane) molecules trapped in water lattices. The continental margins at water depths of greater than 300 m at high latitude or 500 m in temperate latitude sequester a large quantity of methane hydrates (MHs) at low temperature and high pressure. MHs have increasingly attracted worldwide research attention not only because of their economical value as a tremendous source of energy but also due

oo

f

to their potential to cause catastrophic geo-hazards. Hydrate dissociation could prime or directly trigger massive submarine landslides, posing enormous threat to offshore infrastructures and

pr

coastal cities (Maslin et al., 2010; Paull et al., 2011). For instance, the Storegga Slide is believed to have been triggered by hydrate dissociation (Paull et al., 2007). MHs are susceptible to change

e-

in temperature and pressure and may dissociate under natural or artificial perturbations, such as

Pr

sea level drop (Maslin et al., 1998), ocean warming (Skarke et al., 2014), drilling activities, and oil/gas production (Moridis et al., 2018). As MHs dissociate, the host sediments will be severely

al

weakened due to loss of hydrate cementation (Hyodo et al., 2014) and pore-pressure buildup attributed to fluid volume expansion (Xu and Germanovich, 2006) and re-consolidation

rn

(Handwerger et al., 2017). Thus, stability of continental slopes could be severely undermined by

Jo u

extensive hydrate dissociation. In particular, upper continental slopes with water depth ranging from 300 to 500 m are considered the most vulnerable to environmental changes (Mestdagh et al., 2017; Ruppel and Kessler, 2017), in that the base of the stable hydrate zone (BSHZ) reaches the seafloor and a slight variation in temperature or pressure could induce hydrate dissociation. For instance, Phrampus and Hornbach (2012) reported that MHs along a broad swathe of the North American margin are being rapidly destabilized by recent changes in intermediate-depth ocean temperature associated with the Gulf Stream. Hydrate-associated submarine landslides have been investigated in the past two decades with physical modelling (Zhang et al., 2016), analytical models (Nixon and Grozic, 2007; Xiong and Zhang, 2012), and numerical simulations (Jiang et al., 2015; Kwon and Cho, 2012; Sultan et al., 2004). A consensus among these studies is that overpressure must be involved to destabilize a 4

Journal Pre-proof gentle submarine slope undergoing hydrate dissociation (Dugan and Flemings, 2000; Elger et al., 2018; Talling et al., 2014). However, the pore pressure models incorporated in most analyses of hydrate-induced slope instability remain simple. For instance, Nixon & Grozic (2007) quantified the excess pore pressure assuming undrained conditions for a pre-specified amount of dissociation, and proposed a conservative approach for assessing seafloor instability by neglecting pressure diffusion during hydrate dissociation. Xiong and Zhang (2012) took a similar approach but considered the effect of methane solubility. These models do not consider the temporal evolution of pore-pressure development associated with hydrate dissociation, so they cannot predict the

oo

f

commencement of slope destabilization. To more accurately capture pore-pressure buildup associated with hydrate dissociation, Kwon and Cho (2012) created a sequentially coupled model,

pr

which decoupled hydrate dissociation (causing pressure buildup) from consolidation (resulting in

e-

pressure dissipation) in each time step. Due to the lack of accurate models of pore pressure, the spatiotemporal destabilization modes predicted by such slope analyses are not necessarily

Pr

accurate.

A holistic model for slope stability under hydrate dissociation should properly predict the rate of

al

pore pressure development according to hydrates’ thermodynamic behavior and also the

rn

simultaneous dissipation through faults or other pathways in the deposits. This pressure evolution

Jo u

arises from the coupling of chemical reaction, mass and heat transfer, and mechanical deformation. A number of numerical simulators (Kurihara et al., 2009; Rutqvist, 2017; White and McGrail, 2006) have been developed with rigorous implementations of the pore pressure models. These tools have been used to simulate hydrate production ( e.g., Myshakin et al., 2018) and the resultant wellbore instability (e.g., Qiu et al., 2015). However, the applications of such high-fidelity pore pressure models to submarine landslides associated with hydrate dissociation are sparse. This paper investigates the destabilization mechanisms of upper continental slopes undergoing hydrate dissociation using a model rooted in the limit equilibrium analysis with the capacity to capture pore pressure’s generation and dissipation. We build a one-dimensional slope model by correlating the weakening of cohesion to hydrate dissociation and computing transient pore pressure using a thermo-hydro-chemical (THC) coupled analysis. With this model, we investigate 5

Journal Pre-proof the role of overburden deposits above the hydrate-bearing layer (i.e., hydrate reservoir) and identify the key properties of the deposits controlling the slope failure patterns in time and space. Among a number of possible causes of hydrate dissociation, bottom water warming (Skarke et al., 2014) is chosen as an example here, although the proposed method can be extended to other mechanisms inducing hydrate dissociation and the associated seafloor instability. The main contributions of this study include: (1) the development of a slope stability model considering the transient pore water pressure induced by hydrate dissociation, and (2) the recognition of different spatiotemporal destabilization modes of continental slopes undergoing hydrate dissociation with

oo

f

the aid of the developed slope model. This paper is organized as follows. After presenting the framework of limit equilibrium analysis with transient pore pressure, Section 3 formulates and

pr

numerically solves the evolving pore pressure by considering the THC coupled process es of

e-

hydrate dissociation. Section 4 configures a computational model for upper continental slopes undergoing hydrate dissociation due to bottom water warming. Section 5 presents the results with

Pr

an emphasis on the spatiotemporal destabilization modes inferred from a parametric study. Section

2

al

6 highlights the implication of this study and discusses the limits of the presented model.

Limit equilibrium analysis with transient pore pressure

rn

To maintain s implicity, we adopted the infinite slope model for a submarine slope with a buried

Jo u

hydrate reservoir. This infinite-slope approximation has been widely used particularly for submarine slides with the lateral-longitudinal dimensions much greater than the thickness (Kwon and Cho, 2012; Lizárraga et al., 2016; Sultan et al., 2004). This approximation implies that: (1) a planar potential sliding surface running parallel to the seafloor, (2) heat transfer and fluid flow perpendicular to the seafloor, and (3) constant properties of all formations along the direction parallel to the seafloor. Figure 1 illustrates the submarine slope with slope angle  at water depth Hw. The hydrate reservoir with thickness hh runs parallel to the seafloor below an overburden layer with thickness ho . We assume that the process of hydrate dissociation causes a potential slip (denoted with dash line in Figure 1(a)) at vertical depth H from the seafloor. Within the framework of the limit equilibrium method, the safety factor of the slope can be defined for the given failure plane as: 6

Journal Pre-proof Fs 

 f c '  'tan  '  d d

(1)

where d and f are the shear stress and the shear strength in the effective stress on the failure plane, respectively. The shear strength includes the cohesive strength (characterized by the effective cohesion c' ) and the frictional strength depending on the effective normal stress  ' on the failure plane and the effective internal friction angle '. The stress components on the failure plane (i.e. d and  ' ) can be formulated by analyzing the forces on a representative unstable mass (denoted with a white box in Figure 1(a) and detailed in Figure1 (b)), to which a rigid body is the analogue.

f

Given the configuration shown in Figure 1 and by following the formulation based on force

oo

equilibrium given in (Brunsden and Prior, 1984; Chen and Zhang, 2014), the stresses on the failure plane can be written as

pr

 d   ' H sin  cos 

(3)

e-

 '   ' H cos2   ue   '0  ue

(2)

Pr

where  ' is the buoyant unit weight of the soil; ue and  '0 are the excess pore pressure and the initial effective overburden stress on the failure plane. Therefore, Eq. (1) can be re-written as: 𝑐′

𝑢

tan𝜑′

𝐹𝑠 = 𝛾′𝐻sin𝛽cos𝛽 + (1 − 𝜎 𝑒′ ) tan𝛽

(4)

al

0

rn

The term ue is computed from the difference between the total pore pressure (u) and the hydrostatic water pressure (us ) determined from the average water depth of the slope. Note that us

Jo u

increases with increasing water depth. This variation is small at a limited length of a very gentle slope, and therefore neglected in this study. Following Bishop’s principle (Bishop, 1959) for multi-phase flow (e.g., unsaturated soils), the total pore pressure is calculated from: u

Sg Sw uw  ug Sw +S g Sw +S g

(5)

where Sw and Sg are the water and gas saturation (defined as the volume fraction of the pore space occupied by water or gas), respectively; and uw and ug are the pore water pressure and the pore gas pressure, respectively. Equation (4) suggests that the pressure ratio r = ue / '0 has a vital effect on the slope stability. This ratio equals zero under hydrostatic conditions and the slope remains statically stable; the ratio could be one (i.e. liquefaction condition) or even greater due to overpressure in poorly-permeable 7

Journal Pre-proof deposits, which favors slope instability. Equation (4) can be re-written as Fs  Fsc  Fs

(6)

where Fsc and F s  correspond to the first and the second term of the right-hand side of Eq. (4), representing the partial factor of safety contributed by the cohesion and the effective overburden stress at the failure plane, respectively. F sc is negatively correlated to the thickness of the overburden layer; the overburden thickness affects F s  in a less straightforward manner, since it favors the excess pore pressure buildup (yielding a reduction in F s ) and increases the overburden

rn

al

Pr

e-

pr

oo

f

stress (leading to an increased F s ).

Figure 1. A schematic illustration of an infinite gentle submarine slope undergoing hydrate

Jo u

dissociation. The unstable mass is analogy to a rigid block with one unit out-of-plane contact area, denoted with a white box in (a), and the force system of the block is shown in (b).

Equation (6) indicates that a marine slope undergoing hydrate dissociation can be destabilized because of two reasons. First, hydrate dissociation weakens the host deposits and thus reduce the partial factor F sc. As proved by extensive experiments (Masui et al., 2005; Waite et al., 2009; Yoneda et al., 2017), the effective cohesion c' of the host deposits decreases during hydrate dissociation due to loss of cementation, although the effective internal frictional angle ' of the deposits is nearly unaffected. It is noteworthy that the reduction in F sc alone is not sufficient to destabilize a gentle submarine slope where the slope angle is much less than the internal friction angle of the sediments that form the slope. The second reason of destabilization is the pressure 8

Journal Pre-proof buildup causing a dramatic drop in the partial factor F s . Pressure buildup due to hydrate dissociation has been attributed to two mechanisms: 1) gas and water release from solid hydrates, which causes a fluid volume expansion impeded in relatively impervious medium and thereby results in a rapid pressure buildup; 2) the pressure taken by the solid hydrates is transferred to pore fluid following hydrate dissociation (Handwerger et al., 2017), causing an instant increase in pore pressure and a rapid consolidation of the host deposits. In this study we consider the former mechanism, noting that the latter is mostly encountered in high-saturation reservoirs where the host deposits are originally unconsolidated.

f

Pore pressure buildup by hydrate dissociation

oo

3

The evolution of pore pressure during hydrate dissociation arises from an interplay of hydrate

pr

decomposition, multi-phase flow, and heat transfer. Here we formulate the transient pore pressure

e-

by considering these coupling processes.

3.1 Governing equations and numerical implementation

Pr

As illustrated in Figure 2, a mass of MH-bearing sediment (MHBS) can be viewed as a porous system composed of a soil skeleton and three species (i.e. water, methane and hydrate) hosted in

al

the pores of the skeleton. For simplification, we assume that: (1) the temperature is the same

rn

across all phases at the same point at any time; (2) the soil skeleton is poroelastic; (3) hydrates are immobile and incompressible solid in the pores; and (4) the three pore species are each present in

Jo u

an individual phase by ignoring presence of dissolved methane, water vapor and ice. The first three assumptions are widely accepted (Gamwo and Liu, 2010; Kimoto et al., 2010; Moridis and Pruess, 2014); the last assumption enables us to concentrate on the primary effect of hydrate dissociation without being trapped into an over-complicated numerical implementation.

Figure 2. A schematic representation of MHBS as (a) a multiphase, multispecies porous media and (b) the volume fractions of different species/phases. 9

Journal Pre-proof The governing equations of the heat and flow through this idealized system can be formulated based on mass and energy conservation equations. A complete derivation under three-dimensional conditions can be found elsewhere (Chen et al., 2016). Here we present the formulation under one-dimensional (1D) conditions as shown in Figure 1, under which fluid and heat transport occurs in the direction perpendicular to the seafloor (i.e., along z-axis shown in the figure). The mass conservation of the three pore species yields: q md   S   =   +  , t z t

  w, g ,h

(7)

f

where  is the porosity; the subscript  denotes the species in the pores (w: water in aqueous phase;

oo

g: methane in gas; h: hydrate in solid); S, ρ, q are the saturation degree, the density, and the mass flux (i.e., the rate of mass flow per unit area), respectively; and md / t is the mass rate of

pr

the released species  per unit volume of the sediment due to hydrate dissociation. The density of

e-

water is a function of temperature and pressure (Moridis and Pruess, 2014). The change of porosity is connected with the pore pressure by a compressibility coefficient :

Pr

   u

(8)

al

The chemical equation for dissociation (endothermic process) and formation (exothermic process) of MH is written as: CH4 ∙ ηH2 O→ ←

rn

endothermic

CH4 + ηH2 O

(9)

Jo u

exothermic

where η is the hydrate number. Here we use η = 5.75, assuming complete hydration in structure I (Sloan and Koh, 2008). This suggests that forming (or melting) a certain mole of hydrates consumes (or releases) an equal mole of methane; that is, given χ as the reaction rate of hydrates (defined as the mole of hydrates formed or dissociated per unit time), the consuming (or releasing) molar rate of methane and water during hydrate formation (or dissociation) is χ and ηχ, respectively. Here we define consuming as negative and releas ing as positive. Thus, the mass rates of different pore species during formation or dissociation can be computed from: mwd  M w; t

mgd t

 M g ;

mhd  M h  t

(10)

where Mw, Mg and Mh are the molar mass of water, methane and hydrate, respectively. The reaction rate χ can be estimated based on Kim-Bishnoi’s model (Kim et al., 1987): 10

Journal Pre-proof  Ea  1 3 2 3  N0 N h ueq  u g   RT 

  K0 exp  

(11)

where K0 is the intrinsic dissociation constant; E a is the activation energy of hydrate dissociation; R = 8.314 J/mol/K is the universal gas constant; N0 and Nh are the initial and transient moles of hydrates, respectively; and ueq is the equilibrium pressure for stable MHs at temperature the relationship given in Table 1. By introducing Darcy's law, Eq. (7) can be re-written as   w, g , h

(12)

f

 md    k k  u  S   =   eff r      g cos   +  , t z    z   t

oo

where k eff is the effective permeability of the soil depending on the hydrate saturation (see Table 1); k rα and μα are the relative permeability and viscosity of species , respectively. Note that krh is zero

pr

since hydrate is assumed to be immobile, while k rw and krg depend on water and gas saturations

e-

according to several available models (Moridis and Pruess, 2014). Besides mass conservation, energy conservation provides an additional governing equation:

Pr

  T H d h     q c T   1     s cs    Sa  c    z z  w, g t   w, g , h   t

(13)

where cs and s are the specific heat capacity and density of the soil grains, respectively; h is the

(14)

rn

H d  H t

al

heat flux; and 𝜕𝐻 𝑑 ⁄𝜕𝑡 is the dissociation heat rate computed through

Jo u

where ΔH is the latent heat, which can be estimated through Kamath’s equation (1984):

B T A  56552, B  16.8137; 273.15  T  298.15    A  27329, B  50.0647; 248.151  T  273.15 H  A 

(15)

By introducing Fourier's law, Eq. (13) is re-written as   T  T      1     s cs    S  c     t  z z     w, g , h    keff kr   u  H d      g cos   c T      z   w, g    z t  

(16)

where λΘ is the composite thermal conductivity of the system, which can be estimated through the hydrate and water saturation (Reagan and Moridis, 2008) according to an earlier model (Somerton, 11

Journal Pre-proof et al., 1974a; Somerton, et al., 1974b). Table 1. Model and parameters adopted in the heat and flow analysis 1) Composite thermal conductivity model (Moridis et al., 2005)

 





Sh + Sw  sw  sd   sd

sd — Dry thermal conductivity, sd =1.0 W/m/K sw —Wet thermal conductivity, sw =3.3 (W/m/K) 2) Capillary pressure model (Van Genuchten, 1980) 1

 1 

Sirw — Irreducible water saturation, S irw = 19% P 0 — Gas entry pressure value, P 0 = 2000 Pa  — Pore size distribution index,  

f

1/ 

oo

Pc  P0  S *  Sr  Sirw S*  1  Sirw

pr

3) Model for permeability of hydrate-bearing medium (Phillips, 1991)

e-

keff  k (1  Sh )n k — Intrinsic permeability

Pr

n —Permeability reduction exponent, n = 3

4) Relative permeability model (modified version of (Stone, 2007)) n

  w, g

al

 S  Sir  kr     ,  1  Sir 

n— The exponent, n = 4

rn

Sirg —Irreducible gas saturation, S irg = 2% Sirw —Irreducible water saturation, S irw =20%

ln ueq          ln ueq        

Jo u

5) The phase equilibrium relationship (Moridis and Pruess, 2014)

 1.94138504464560  105 3.31018213397926  103 T 2.25540264493806  101T 2 7.67559117787059  102 T 3

，T  273.2 K

1.30465829788791 104 T 4 8.86065316687571 108 T 5  4.38921173434628  101 7.76302133739303  101T 7.27291427030502  103 T 2 3.85413985900724  105 T 3

，T  273.2 K

1.03669656828834  107 T 4 1.09882180475307  1010 T 5

12

Journal Pre-proof In summary, Eqs. (12) and (16) form a set of governing equations of the system (consisting of four individual equations) that involve six unknown variables: pore gas pressure ug, pore water pressure uw, temperature T, water saturation Sw, gas saturation Sg and hydrate saturation Sh . These equations are closed by the following auxiliary equations: S w  S g  Sh  1

(17)

ug  uw  uc

(18)

where uc is the capillary pressure to be determined by van Genuchten’s function (1980). Table 1

oo

f

summarizes the major relations or models used for calculating the necessary parameters in this study.

pr

No analytical solution is available for the governing equations due to nonlinearity of many relations in use (e.g., permeability and capillary pressure). We achieved a numerical solution using

e-

a MATLAB code by linearizing the governing equations with the implicit finite-difference method.

Pr

The central and forward difference approximations were used for the spatial and time derivatives, respectively. Then, the linearized difference equation set was solved with the Newton-Raphson

al

iteration method. We followed the implementation given by Gamwo and Liu (2010), and the

3.2 Validation

rn

validity of our implementation is illustrated through two cases in the next section.

Jo u

To validate our code for the transient pore pressure associated with hydrate dissociation, we replicated two simulations (Tang et al., 2007; Reagan and Moridis, 2008) performed with TOUGH+HYDRAT E, a simulator for nonisothermal multiphase flow under conditions typical of common natural methane hydrate deposits (Moridis, 2014). TOUGH+HYDRATE is among the most prevailing codes for hydrate-related simulations nowadays and has received many applications (Moridis et al., 2018; Myshakin et al., 2018). This replication is a necessary step but not the ultimate goal of our study. The validated code is for achieving the transient pore pressure required for a slope stability analysis, which has not been studied in the literature.

3.2.1

1D dissociation test

We validated our code by simulating a dissociation test conducted in a cylindrical, stainles s steel pressure vessel with an internal radius of 1.9 cm and a length of 50 cm (Tang et al., 2007). The 13

Journal Pre-proof vessel cell was jacketed with an adiabatic, impervious layer, and filled with a sand (grain size: 300 ~ 450 µm, porosity: 30.8%, intrinsic permeability: 3 × 10-13 m2 ) containing methane hydrates and water in 21.83% and 29.61% of the void volume, respectively. Initially, the cell was at a pressure of 3.535 MPa and a temperature of 1.54 ºC. To trigger hydrate dissociation, the outlet valve was opened and maintained at a constant pressure of 0.93 MPa and a temperature of 1.54 ºC. The gas released through the outlet valve was collected. Tang et al. (2007) also reported their simulation data obtained with TOUGH+HYDRAT E. We adopt the same relationships (e.g., relative permeability model, hydrate dissociation model) and

oo

f

corresponding parameters as reported. To be consistent with the reported simulation, the reaction

 Ea   f A As ueq  u g   RT 

(19)

e-

  K 0 exp  

pr

rate of hydrate dissociation was computed from the following rather than Eq. (11):

where A s and f A are the hydrate reaction area and the adjustment factor, respectively. More detailed

Pr

information about the model parameters can be read from Tang et al. (2007). Figure 3(a) illustrates the geometric and boundary conditions of the numerical model of the 1D

al

dissociation test. In the experiment, the 1D flow assumption is invalid near the outlet valve of the

rn

cell due to boundary effect. However, for simplification, we assume 1D flow in the entire domain and use an intrinsic permeability (i.e., the so-called interface permeability k inter ) in the boundary

Jo u

element different from that in the other elements in order to resemble the 3D effect of the outlet. Therefore, the formulation in Section 3.1 remains valid for the boundary element as well. Since k inter represents the overall flowability near the outlet, it is difficult to measure. Instead, we performed a sensitivity study on k inter. Figure 3(b) compares our simulation with the test and the reported simulation. With k inter = 3×10-16 m2, our model is consistent with the test data and the reported simulation.

14

Cumulative volume of gas at STP(cm3)

Journal Pre-proof

1x104

3x10-15 kinter= 3x10-16 m2

8x103

6x103

4x103 Tang et al. 2007 (test data) Tang et al. 2007 (simulation) This study

2x103

0

0

10

20

30

40

50

Elapsed time (min)

(a)

(b)

3.2.2

pr

oo

f

Figure 3. Validation of a 1D dissociation test: (a) the numerical model setup; (b) comparison between experimental and numerical results.

An overburden-free hydrate reservoir subject to ocean warming

e-

We re-visited a case study simulated by Reagan and Moridis (2008) with TOUGH+ HYDRAT E for

Pr

investigating possible methane release from hydrate dissociation. This case represents a typical scenario of warm and shallow hydrate reservoirs encountered in upper continental slopes of the

al

Gulf of Mexico. Below a water depth of 570 m, a hydrate reservoir extends vertically from the seafloor to a depth of 16 m, and the initial hydrate saturation is 3%. The intrinsic permeability of

rn

the host sediment is 10-15 m2 (representative of a silty reservoir). Initially, the temperature is 6 ºC

Jo u

at the seafloor and linearly increases with depth at a geothermal gradient of 26 ºC/km inferred from the reported profile of initial temperature. Hypothetically, hydrates thermally dissociate because of a temperature rise at the seafloor at an annual rate of 0.03 ºC as attributed to a climate change scenario. We replicated the boundary/initial conditions and adopted the models (see Table 2) associated with hydrate dissociation as reported in Reagan and Moridis (2008). Different from Reagan and Moridis (2008), we assume insoluble methane and a physical absence of inhibitors (e.g. salts) to maintain simplic ity in our model. The effect of methane solubility is not significant and will be discussed in section 6.3. In the natural oceanic hydrate reservoir, the presence of salts dissolved in water impacts the equilibrium conditions of hydrates. However, the salinity dependence can be indirectly considered by justifying the temperature in the phase 15

Journal Pre-proof

equation of hydrate, given that an increase of salinity results in a dec rease of the equilibrium temperature at an ambient pressure (Moridis, 2014). To tackle the effect of inhibitors, we shifted the phase equation of hydrate (provided in Table 1) by a temperature decrease ∆𝑇𝐷 calculated from (Moridis, 2014): TD =TD , r

  ln 1  Y  ln 1  YAi

(20)

i A, r

f

𝑖 where 𝑌𝐴𝑖 is the inhibitor mole fraction in the aqueous phase for a given salinity; 𝑌𝐴,𝑟 is the

oo

reference mole fraction of the inhibitor corresponding to a temperature decrease ∆𝑇𝐷,𝑟 . Here we

pr

𝑖 = 0.01335 and ∆𝑇 = 2 ºC. Given salinity 0.035 in mass fraction in this study, the used 𝑌𝐴,𝑟 𝐷,𝑟

e-

phase equation is shifted to a lower temperature by ∆𝑇𝐷 = 1.65 o C. Figure 4 shows that our results are in good agreement with the reported ones. Both studies suggest

Pr

that hydrate dissociation starts at the top and base of the hydrate reservoir. The base of the reservoir locates at the lowermost border of the stable hydrate zone, and therefore a slight

al

downward heat flux is sufficient to trigger dissociation at that locality. Nevertheless, dissociation

rn

is more prominent at the top of the reservoir, resulting in a quicker downward propagation of dissociation front from the seafloor. The hydrate reservoir is depleted before 100 yrs in our

Jo u

simulation and the reported one, since the profiles of hydrate saturation obtained from both simulations at 100 yrs are vertical overlapping the axis Sh =0. We also noted that the dissociation front propagation is slightly delayed in our model at 10 yrs. This discrepancy is consistent with our assumption of insoluble methane (that is, all the released methane is in the gaseous phase only). This assumption leads to overestimated pore pressure in the dissociated zone, implied in Figure 4b by comparing the gas saturation profiles. The gas saturation at the dissociated zones obtained from our model is larger than the reported one. Due to higher pore pressure, more heat input is required for further dissociation and thus dissociation is slower in our model.

16

Journal Pre-proof Table 2. Parameters of the slope model Layer

Hydrate reservoir

Parameters Intrinsic permeability (m2 )

Value 1×10-15

Thickness (m) Initial hydrate saturation Sh0 (%)

16 30

Initial gas saturation (%) Internal frictional angle (º)

0 35 Varying 250～400

Cohesion (kPa) Intrinsic permeability (m2 ) Thickness (m)

Overburden layer

Varied, 10-19 ~10-15 Varied, 1～100 Varied, 5～200

Cohesion (kPa) Internal friction angle (º)

35 2

-15

Internal friction angle (º) Cohesion (kPa)

35 20

pr

layer

oo

Underburden

10 344

f

Intrinsic permeability (m ) Thickness (m)

e-

Slope angle (º) Salinity

28

Hydrate number Initial porosity (%) Pore compressibility (MPa) Unit weight of water (kN/m3 )

5.75 30 -8 10 10.37

Saturated unit weight of deposits (kN/m 3 ) Soil grain density (kg/m3 )

16.5 2650

Hydrate density (kg/m3 ) Specific heat capacity of soil grain ( J/(kg·ºC))

9200 1000

Specific heat capacity of gas ( J/(kg·ºC)) Specific heat capacity of water ( J/(kg·ºC))

2160 4200

Specific heat capacity of hydrate ( J/(kg·ºC))

2080

Viscosity of gas (Pa•s) Viscosity of water (Pa•s)

10 10-3

al

Pr

Geothermal gradient (ºC/km)

Jo u

rn

All formations

3 0.035

-5

17

Journal Pre-proof

0

0

Depth (mbsf)

5

Reagan & Moridis, 2008 10 yr 100 yr This study 10 yr 100 yr

5 Hydrate reservoir

10

10

15

15 Underburden

20

0

1

2

(a) 3

(b) 4

20

0

Sh (%)

2

4

6

8 10

f

Sg (%)

oo

Figure 4. Revisit of an overburden-free hydrate reservoir under shallow water: (a) the profiles of

pr

the hydrate saturation, and (b) the profiles of the gas saturation.

The excess pore pressure induced by hydrate dissociation is in the order of several kilopascals (the

e-

plot is not presented). This pressure is insufficient to destabilize a gentle cohesive slope common

Pr

in submarine settings. To examine a wider range of vulnerable upper continental slopes, we extend this case study towards less stable scenarios in the subsequent sections to probe the spatiotemporal

Model setup

rn

4

al

modes of slope destabilization attributed to hydrate dissociation.

We modify the preceding case study towards more vulnerable conditions by considering an

Jo u

overburden above the hydrate reservoir with Sh = 30%. This hydrate saturation is selected corresponding to the average value in submarine sediments (Fujii et al., 2015; Wang et al., 2014). o

We assume a constant slope angle 3 to represent a typical continental slope (Haflidason et al., 2004; He et al., 2014; Hühnerbach and Masson, 2004). The ocean warming scenario with an annual rate of 0.03 ºC was maintained, although a long-term continuous temperature rise may not be realistic. With this hypothetical triggering mechanism, w e concentrate on a theoretical study of possible destabilization modes of slopes constrained with a range of typical conditions rather than a detailed field-specific study. Other triggering mechanisms are discussed in Section 6.2. Tables 1 and 2 summarizes the models and parameters adopted in the study, respectively. Based on a thorough review on mechanical properties of MHBS, we chose the strength parameters obtained by Miyazaki et al. (2011) from drained triaxial compression tests on synthetic MHBS specimens 18

Journal Pre-proof that resembled the conditions of Nankai Trough of Japan. As suggested, the change in hydrate saturation has a marginal effect on the internal friction angle of the host soil but a considerable effect on the cohesion. Accordingly, we assumed that the internal frictional angle remains 35º, while the cohesion of the hydrate reservoir linearly decreases from 400 kPa to 250 kPa as the hydrate saturation reduces from 30% to zero. Although more sophisticated relationships have been empirically established to link MHBS’s cohesion to hydrate saturation (Liu et al., 2017; Masui et al., 2005; Salehabadi et al., 2008; Yoneda et al., 2015), we used a linear one since a universal relationship is still lacking.

oo

f

To understand the role of the overburden in affecting the destabilization modes of the slope, we conducted a sensitivity study by varying the thickness ho, cohesion co and permeability k o of the

pr

overburden based on a baseline case with ho = 20 m, co = 50 kPa, and ko = 10-17 m2. The thickness

e-

ho ranges from 0 to 100 m to encompass a typical range inferred for submarine lands lides (Chaytor et al., 2009), the cohesion co ranges from 5 to 200 kPa to cover low-cohesion to cohesive

Pr

overburden, and the permeability ranges from 10-19 m2 to 10-15 m2 to represent a silty or clayey

5

Results

al

cover of reservoirs being considered optimal sites for production (Li et al., 2013; Sun et al., 2015).

rn

5.1 Propagation of dissociation front and induced pore pressure Figures 5a-b illustrate the evolving profiles of temperature and hydrate saturation obtained from

Jo u

the baseline case. The results obtained from the other cases exhibit a similar trend in general. As the heat is transferred from the seafloor downward through the marine deposits to the hydrate reservoir, dissociation starts at the top of the reservoir after 18-year rise of temperature and the dissociation front moves downward over time. Although not significant, a local dissociation is also found at the base of the reservoir, which, lining up with the bottom of the hydrate stable zone, is sensible to variation in temperature and pressure. Dissociation terminates in 276 yrs after all hydrates melt. Since hydrate dissociation is an endothermic process, a reduction in temperature is evident in the proximity of the dissociation front at 276 yrs as shown Figure 5a. Figures 5c-f present the evolving gas saturation, induced capillary pressure, excess pressure, and pressure ratio. As shown in Figure 5c, liberated gas accumulates at the top and base of the 19

Journal Pre-proof reservoir in the early stage of temperature rising, producing two peaks of gas saturation in the profile. As the upper dissociation zone propagates downwards, the gas saturation increases over the entire reservoir and reaches the maximum as the reservoir completely melts (in 276 yrs); meanwhile, the light gas migrates into the overburden layer, resulting in an increase in the gas saturation (and a decrease in the water saturation) of the overburden layer. After dissociation terminates due to reservoir depletion, the gas saturation decreases throughout the entire depth as the gas is released from the seafloor. As shown in Figure 5d, capillary pressure appears as gas liberates and migrates in the deposits producing multi-phase flow. Accordingly, as shown in Figure

oo

f

5e, the excess pore pressure builds up forming an upward-moving peak during hydrate dissociation. Meanwhile, the excess water pressure diffuses upwards and downwards. However,

pore pressure is noticeable after reservoir depletion.

Pr

20

30

al

30

40

10 15 T (ºC)

20

0

20

0

60

80

0

0

0

(d)

40

20

30

40 0

10

20

30

5

10

15

20

uc (kPa) 0 yr 200 yr

(g)

20

30

40 0

Sh (%)

400

10

Depth (mbsf)

Depth (mbsf)

30

200

ue (kPa)

10

20

30

Sg (%)

(b)

10

20

40 40

Jo u

5

(e)

10

rn

20

40

Depth (mbsf)

(c)

10

Depth (mbsf)

Depth (mbsf)

10

0

e-

0

(a)

Depth (mbsf)

0

pr

pressure diffusion only dominates after dissociation completely terminates, so reduction of excess

18 yr (diss. start) 276 yr (diss. end)

-1 0

1

2

3

4

5

r 70 yr 400 yr

122 yr (sliding start) 500 yr

Figure 5. (color online) Evolving profiles of (a) temperature, (b) hydrate saturation, (c) gas saturation, (d) capillary pressure, (e) excess pore pressure, and (f) pressure ratio obtained from the baseline simulation. The initial profile of hydrate saturation (0 yr) is covered by the one of 18 yr, so it is not shown in (b). 20

Journal Pre-proof 5.2 Evolution of slope stability Figure 6 presents the evolution of the partial factors of safety (F s and F sc) and the overall factor of safety (F s ) obtained from the baseline case. As shown in Figure 6a, F s  is initially constant along depth because of a constant internal friction angle in all formations. With pressure build-up due to hydrates dissociation, F s drops and forms a valley (i.e., the minimum value) in the overburden layer, which propagates upwards over time due to an upward-migrating peak of excess pore pressure (see Figure 5e). After dissociation ceases, Fs recovers as the excess pore pressure dissipates through the seafloor. As shown in Figure 6b, the initial profile of F sc decreases with

oo

f

depth in each layer; F sc is higher in the reservoir than that in the overburden due to higher cohesion attributed to hydrate cementation. As hydrates dissociate, F sc declines in the dissociated

pr

zone due to reduced cohesion, and F sc remains unaffected in other formations. We note that the

e-

reduction in F sc alone will not destabilize a gentle marine slope, and a declining F s is a necessity for slope instability. Figure 6c shows the evolving profiles of Fs , from which the magnitude and

Pr

the position of the minimum F s (i.e. F s,min ) can be recognized to indicate the safety factor and the depth of the most dangerous plane (i.e., the potential sliding plane), respectively. Note that the

al

factor of safety becomes negative between 122 to 400 years, because the pore pressure becomes so

rn

great that it significantly exceeds the normal stress on the sliding surface. A negative factor of safety should be interpreted as slope failure, although the magnitude of the negative factor is

Jo u

non-physical since the limit equilibrium slope analysis model becomes invalid. Figure 7 presents Fs,min and the depth of the potential sliding plane over time obtained from the baseline simulation. The magnitude of Fs,min decreases due to pressure build-up and de-cementation in the sediments as hydrates dissociate. It drops to below 1 right above the interface between the reservoir and the overburden after 122-year rising in temperature before the reservoir completely melts. This indicates occurrence of a slide (denoted as point b). For easy interpretation, we define the timing (the depth) of sliding onset as the timing (the depth) when (where) F s,min first reaches one. Intuitively, we expect different scenarios in terms of timing and position of sliding onset under various properties of the overburden. This is confirmed through a parametric study to be presented subsequently.

21

Journal Pre-proof 0

(a)

0

0

(b)

0 yr

(c)

18 yr (diss. start) 70 yr

Depth (mbsf)

10

10

10

122 yr (sliding start) 200 yr 276 yr (diss. end)

20

20

20

30

30

30

40 40 -45 -30 -15 0 15 30 0 Fs

25

50

75

400 yr 500 yr

40 -25

0

25 Fs

50

75

f

Fsc

oo

Figure 6. (color online) Evolving profiles of (a) stress-associated partial factor of safety, (b)

pr

0

e-

16

a

Factor of safety Depth

12

a: Diss. starts (t=18 yr) b: Slide onset (t=122 yr)

Pr

8

Overburden layer

Hydrate reservoir

0

30

rn

1 0

60

20 Post-failure

al

4

10

b

90

120

30

Depth of potential sliding surface (mbsf)

Minimum safety factor Fs,min

cohesion-associated partial factor of safety, and (c) the overall factor of safety obtained from the baseline simulation.

150

Jo u

Elapsed time (yr)

Figure 7. Evolution of the safety factor and the potential sliding plane during dissociation.

5.3 Co- and post-melting sliding In general, the timing of sliding onset is governed by a competition between buildup and dissipation of excess pore pressure, which is impacted by the thickness and the permeability of the overburden. Besides, the cohesion of the overburden plays a secondary role in slope stability. The timing of sliding onset under different configurations of the overburden were investigated through a parametric study on the 1D slope stability model, and Figure 8 presents the results. As shown in Figure 8, the timing of sliding onset is visualized as a three-dimensional color-coded surface, which is developed from polynomial fitting based on a set of data points obtained from the 22

Journal Pre-proof parametric study. For easy interpretation, the surface is projected on a two-dimensional space formed by the axes of thickness and cohesion. For comparison, the timing of dissociation onset and termination (denoted as the wire-frame surfaces) is superimposed in the same figure. When -16

the overburden is relatively permeable with k o =10

2

m (see Figure 8a), slope instability occurs in

cases with a low-cohesion and relatively thin (less than 50 m thick) overburden. As the overburden becomes less permeable (see Figure 8b-c), slope failure is possible in cases with overburden sediments of higher cohesion up to 200 kPa, and however hindered in an overburden thicker than 60 m. It is noteworthy that the color-coded surface lies in between the two wire-frame surfaces 2

m , indicating that slope instability takes place during the process of

f

-18

oo

when k o is no more than 10

hydrate dissociation. This type of instability is termed as co-melting sliding for short in this paper.

pr

As shown in Figure 8d, the color-coded surface extrudes above the upper wire-frame surface when

e-

the permeability of the overburden sediments is reduced further to 10 -19 m2. This indicates that a slide occurs after dissociation completely terminates. Such a slide, termed as post-melting sliding

Pr

for short, is possible if the overburden is relatively thick (thicker than 50 m) with cohesion less than 100 kPa, beyond which the slope otherwise becomes stable (denoted as the gray area in the

al

projected maps). Under these scenarios, the released gas is trapped at the base of very tight

rn

overburden as hydrates dissociate. Although the pressure is high at the base of the overburden, sliding is hindered by the very thick overburden that provides high overburden stress. However,

Jo u

the released gas keeps ascending slowly over time, forming relatively high excess pressure in the upper zone of the overburden where the stress is less and thus a post-melting sliding is triggered. The black dashed line in the surface projection delineates the transition between co- and post-melting sliding. We interpret the line as a narrow zone of transition rather than a clear cut border because of fitting errors. Note that the application of different fitting methods will affect the smoothness of the surfaces (which is beyond the scope of this study) but will not change the main conclusion of this study on three distinct spatiotemporal destabilization modes of the slope.

23

e-

pr

oo

f

Journal Pre-proof

Pr

Figure 8. (color online) The timing of sliding onset under various conditions of the overburden layer with different permeabilities: (a) k o =10-16 m2; (b) k o =10-17 m2 ; (c) k o =10-18 m2; (d) k o =10-19 m2. The base is the projection of the surface of the timing. The gray areas in the base indicate that the slope

rn

al

remains stable.

5.4 Depth of the potential sliding surface

Jo u

Figure 9a presents the depth of the potential sliding plane in a parametric space forming by overburden thickness and cohesion. The light blue surface represents the position of the seafloor. Figures 9b to 9e project the surfaces of sliding depth shown in Figure 9a into two-dimensional maps under different values of the permeability. The slide conditionally occurs either adjacent to or high above the interface between the overburden and the reservoir depending on permeability, cohesion and thickness of the overburden. For a relatively permeable overburden (see Figure 9b with ko = 10

-16

2

m ), the excess pore pressure diffuses upwards rapidly, and sliding occurs high

above the interface within a relatively thin and low-cohesion overburden. In contrast, in less permeable cases (see Figure 9c-d with k o ranging from 10

-17

to 10

-18

2

m ), the sliding plane is

likely localized near the interface, since pressure diffusion is much longer than that of pressure buildup ruled by hydrate dissociation, the liberated gas and water is trapped at the interface, 24

Journal Pre-proof generating high pressure that exceeds the overburden stress provided by a relatively thin overburden. In cases with overburden sediments of low permeability with k o = 10-19 m2 and low cohesion (see Figure 9e), we note a transition from the interface sliding to the non-interface sliding when the overburden layer is thicker, since the sliding pattern turns to post-melting s liding

Jo u

rn

al

Pr

e-

pr

oo

f

from co-melting sliding.

Figure 9. (color online) The depth of the potential sliding plane under various conditions of the overburden layer. (a) plots four three-dimensional surfaces under different values of permeability to represent the sliding depth (i.e., the distance from the interface between the hydrate reservoir and the overburden) at different overburden thickness and cohesion. (b-e) show the projected surfaces under a given permeability.

6 Implications and discussion 6.1 Destabilization modes Informed by the preceding parametric study, we identify three possible distinct spatiotemporal 25

Journal Pre-proof modes of destabilization of upper continental slopes undergoing hydrates dissociation under different characteristics of the overburden, namely, (1) co-melting non-interface sliding, (2) co-melting interface sliding, and (3) post-melding non-interface sliding. The examples of different modes are marked as point A, C, D in Figures 8 and 9, respectively. Figure 10 presents the evolution of the safety factor and the position of the potential sliding plane corresponding to the three examples. For comparison, a stable example without sliding (i.e., B) is also included.

B

Diss. starts

Diss. ends

Slide onset

10 D

5

A

C

Fs = 1

0 200

400

600

800

1000

Pr

Elapsed time (yr)

e-

0

50

(b) 40

al

D 30

rn

20

A

10

B

C 0 0

Jo u

Height of potential sliding surface from reservoir top (m)

oo

15

f

A (Co-melting non-interface sliding) C (Co-melting interface sliding) D (Post-melting non-interfaces sliding) B (Stable)

(a)

pr

Minimum factor of safety Fs

20

200

400

Hydrate reservoir

600

800

1000

Elapsed time (yr)

Figure 10. The evolution of the safety factor and the potential sliding surface under different destabilization modes. The conditions of case A-D are shown in Figures 8 and 9.

Inspired by the concept of pipe structures formed by migrating overpressured fluid (Elger et al. 2018), we illustrate the three recognized scenarios in two dimensions as shown in Figure 11. In general, co-melting s liding occurs in shallow reservoirs. As shown in Figure 11a-b, in case of a shallow reservoir covered by relatively pervious deposits, the liberated fluid finds pathways and migrates upward so that the overpressure in the reservoir is insufficient to trigger interface sliding. However, the slope could fail retrogressively above the interface, as the migrating fluid encounters 26

Journal Pre-proof an upper zone of loose and low-cohesion deposits where a slight overpressure is adequate to cause instability. If the cover is very tight so that extremely-high pressure is generated below the base of the cover as illustrated in Figure 11c-d, interface sliding could be possible. As shown in Figure 11e-f, in case of deep reservoirs covered by tight overburden deposits, although a large volume of fluid is trapped in the reservoir forming a high-pressure zone, the slope could remain stable till the reservoir is depleted. However, the overpressured fluid continues to escape, although slowly, through the tight cover, and slope failure can be triggered if the upward fluid encounters a shallow weak zone in the overburden. Such failure is often delayed by slow fluid migration and finally

Jo u

rn

al

Pr

e-

pr

oo

f

takes place long after reservoir depletion.

Figure 11. A schematic illustration of three modes of sliding (modified from Elger et al. 2018): co-melting non-interface sliding (a to b), co-melting interface sliding (c to d), and post-melting non-interface sliding (e to f). 27

Journal Pre-proof 6.2 Effect of internal friction angle of the overburden We investigate the effect of the internal friction angle ( ) of the overburden by repeating the slope



with

=

20º.

e-

pr

oo

f

analyses

Pr

Figure 12 compares the results developed from different values of  for a reservoir with low permeability of the overburden (ko = 10-19 m2 ). The gray areas indicate the stable conditions, i.e.

al

the stable zone. The dashed line, corresponding to the pore ratio r = 1, divides the unstable zone

rn

into two parts: compressive shear failure (the shaded area with r < 1) and tensile shear failure

Jo u

(outside the shaded area with r > 1). The variation of the internal friction angle does not eliminate or add spatiotemporal destabilization modes of the slope, while it alters the conditions of the transition (denoted by the bold solid lines) among different modes. Another observation is that reduction in the internal friction angle widens the unstable zone of compressive shear failure, since Mohr-Coulomb failure criterion is adopted in the model for both tension and compression regimes.

28

f

Journal Pre-proof

oo

Figure 12. The effect of the internal friction angle (  on the spatiotemporal destabilization modes -19 2 given k o = 10 m : (a)  = 35º (a), (b)  = 20º. The gray areas indicate stable conditions, and the

Pr

6.3 Effect of dissolution of methane

e-

pr

shaded areas roughly indicate compressive shear failure with r < 1. The bold solid lines delineate the transitions among three destabilization modes: co-melting non-interface sliding (mode I), co-melting interface sliding (mode II), and post-melting non-interface sliding (mode III).

To maintain s implicity, the proposed pore-pressure model does not account for methane solubility.

al

This approximation overestimates the volume of free gas in the system and thus the magnitude of pore-pressure buildup induced by hydrate dissociation. Given gas pressure 6.5 MPa and

rn

temperature 7o C, the methane solubility is about 2.0 g/kg. Ignoring this methane solubility, we

Jo u

overestimate the excess pore pressure by less than 2%, which is considered an acceptable error.

6.4 Limits of infinite slope model There are several limits stemmed from the assumption of the infinite slope, which however can be extended to two- or three-dimensional slope stability analysis models. Constrained by this assumption, hydrates are assumed to be located at a constant depth below seabed and constant saturation along the failure plane. However, gas hydrates can form both in nodules and pockets as well as in continuous layers. For instance, gas hydrates are found to be vein or nodular in some areas of the South China Sea (Zhang et al., 2014); therefore, the content of gas hydrates would possibly vary in the vertical direction and along the failure plane. The infinite slope model cannot tackle the influence of such variation in hydrate saturation. As in situ detection methods develop, the accuracy of determining the spatial distribution of hydrates will be improved 29

Journal Pre-proof and further research will be required to account for the effect of the variation of hydrates on the slope stability in two or three dimensions. Mechanisms of hydrate dissociation may be of natural or artificial origins. Ocean warming is considered in this study the main factor in large-scale dissociation over large areas. Such extensive dissociation of hydrate and the resulting massive slope failure can be approximated with the simplified infinite slope model. In contrast, man-made perturbation mechanisms (e.g., heat transfer during oil or gas production) cause hydrate dissociation in the vicinity of the extraction structure, resulting in local, or small-scale, dissociation that may be detrimental to the engineering

oo

f

projects (e.g., foundation piles, pipelines, extraction wells and offshore infrastructures). Therefore, two- or three-dimensional s lope stability models become fundamental to investigate the harmful

Conclusions

e-

7

pr

impact of man-made mechanisms causing local disturbances to hydrate stability.

This paper explores the destabilization mechanism of an upper continental slope with a hydrate

Pr

reservoir undergoing dissociation. Within the framework of the limit equilibrium analysis, we developed a slope model by computing the transient pore pressure from a high-fidelity procedure

al

for THC coupled analyses of hydrate dissociation. With this model, we capture the spatiotemporal

rn

destabilization modes depending on the characteristics of the overburden above the reservoir. The timing and depth of slope failure vary in a parametric space formed by the thickness, the

Jo u

permeability, and the cohesion of the overburden deposits. Three distinct sliding modes are identified under different combinations of the parameters under consideration. In general, co-melting sliding is possible in a shallow reservoir, otherwise post-melting sliding becomes feasible. In cases of shallow reservoirs (less than about 50 m thick in our simulations), a transition is noted from non-interface sliding to interface sliding as the permeability of the overburden deposits decreases. In cases of deep reservoirs (thicker than about 50 m in our simulations), slope failure could occur after reservoir is depleted and post-melting migration of fluid encounters an upper zone of loose and low-cohesion deposits in the overburden layer where a non-interface sliding is triggered. These findings are helpful to delineate the vulnerable configurations of the continental slopes hosting hydrates. 30

Journal Pre-proof

Acknowledgement The work was supported by the Chinese National Natural Science Foundation (grant numbers 41572267, 51639008, and 41877241) and the Program of Shanghai Academic Research Leader (grant No. 17XD1403700). The collaboration between Tongji University (China) and the University of Milano-Biccoca (Italy) was financially supported by 111 Project (grant number B14017) and International Cooperative Training Program for Postgraduates with Interdisciplinary and Innovative Academic Talents (grant number 2018XKJC-015). This study was also partly

f

funded by Project MIUR – Dipartimenti di Eccellenza 2018–2022.

oo

Reference

Bishop, A.W., 1959. The principle of effective stress. Tek. Ukebl. 39, 859–863.

pr

Brunsden, D., Prior, D.B., 1984. Slope instability.

Chaytor, J.D., Brink, U.S. ten, Solow, A.R., Andrews, B.D., 2009. Size distribu tion of submarine

e-

landslides along the US Atlantic margin. Mar. Geol. 264, 16–27. https://doi.org/10.1016/j.margeo.2008.08.007

Pr

Chen, H.X., Zhang, L.M., 2014. A physically-based distributed cell model for predicting regional rainfall-induced shallow slope failures. Eng. Geol. 176, 79–92. https://doi.org/10.1016/j.enggeo.2014.04.011

Chen, L., Yamada, H., Kanda, Y., Lacaille, G., Shoji, E., Okajima, J., Ko miya, A., Maruyama, S., 2016.

al

Numerical analysis of core-scale methane hydrate dissociation dynamics and multiphase flow in porous media. Chem. Eng. Sci. 153, 221–235. https://doi.org/10.1016/j.ces.2016.07.035

rn

Dugan, B., Flemings, P.B., 2000. Overpressure and fluid flow in the New Jersey continental slope : Implications for slope failure and cold seeps. Science (80-. ). 289, 288–291.

Jo u

https://doi.org/10.1126/science.289.5477.288 Elger, J., Berndt, C., Rüpke, L., Krastel, S., Gross, F., Geissler, W.H., 2018. Submarine slope failures due to pipe structure formation. Nat. Commun. 9, 715. https://doi.org/10.1038/s41467-018-03176-1 Fujii, T., Suzuki, K., Takayama, T., Tamaki, M., Komatsu, Y., Konno, Y., Yoneda, J., Yamamoto, K., Nagao, J., 2015. Geological setting and characterization of a methane hydrate reservoir distributed at the first offshore production test site on the Daini-Atsumi Knoll in the eastern Nankai Trough, Japan. Mar. Pet. Geol. 66, 310–322. https://doi.org/10.1016/j.marpetgeo.2015.02.037 Gamwo, I.K., Liu, Y., 2010. Mathematical modeling and numerical simulation of methane production in a hydrate reservoir. Ind. Eng. Chem. Res. 49, 5231–5245. https://doi.org/10.1021/ie901452v Haflidason, H., Sejrup, H.P., Nygård, A., Mienert, J., Bryn, P., Lien, R., Forsberg, C.F., Berg, K., Masson, D., 2004. The Storegga Slide: Architecture, geometry and slide develop ment. Mar. Geol. 213, 201–234. https://doi.org/10.1016/j.margeo.2004.10.007 Handwerger, A.L., Rempel, A.W., Skarbek, R.M., 2017. Submarine landslides triggered by destabilization of high-saturation hydrate anomalies. Geochemistry Geophys. Geosystems 18, 31

Journal Pre-proof 2429–2445. https://doi.org/10.1002/2016GC006706 He, Y., Zhong, G., Wang, L., Kuang, Z., 2014. Characteristics and occurrence of submarine canyon-associated landslides in the middle of the northern continental slope , South China Sea. Mar. Pet. Geol. 57, 546–560. https://doi.org/10.1016/j.marpetgeo.2014.07.003 Hühnerbach, V., Masson, D.G., 2004. Landslides in the North Atlantic and its adjacent seas: An analysis of their morphology, setting and behaviour. Mar. Geol. 213, 343–362. https://doi.org/10.1016/j.margeo.2004.10.013 Hyodo, M., Li, Y., Yoneda, J., Nakata, Y., Yoshimoto, N., Nishimura, A., 2014. Effects of dissociation on the shear strength and deformation behavior of methane hydrate-bearing sediments. Mar. Pet. Geol. 51, 52–62. https://doi.org/10.1016/j.marpetgeo.2013.11.015 Jiang, M., Sun, C., Crosta, G.B., Zhang, W., 2015. A study of submarine steep slope failures triggered 190, 1–16. https://doi.org/10.1016/j.enggeo.2015.02.007

f

by thermal dissociation of methane hydrates using a coupled CFD-DEM approach. Eng. Geol.

oo

Kamath, V.A., 1984. Study of heat transfer characteristics during dissociation of gas hydrates in porous media. University of Pittsburgh, Pittsburgh, PA.

pr

Kim, H.C., Bishnoi, P.R., Heidemann, R.A., Rizvi, S.S.H., 1987. Kinetics of methane hydrate decomposition. Chem. Eng. Sci. 42, 1645–1653. https://doi.org/10.1016/0009-2509(87)80169-0 Kimoto, S., Oka, F., Fushita, T., 2010. A chemo-thermo-mechanically coupled analysis of ground

e-

deformation induced by gas hydrate dissociation. Int. J. Mech. Sci. 52, 365–376. https://doi.org/10.1016/j.ijmecsci.2009.10.008

Pr

Kurihara, M., Sato, A., Ouchi, H., Narita, H., Masuda, Y., Saeki, T., Fujii, T., 2009. Prediction of gas productivity fro m Eastern Nankai Trough methane hydrate reservoirs. SPE Reserv. Eval. Eng . 12, 477–499.

al

Kwon, T., Cho, G., 2012. Submarine slope failure primed and triggered by bottom water warming in oceanic hydrate-bearing deposits. Energies 5, 2849–2873. https://doi.org/10.3390/en5082849

rn

Li, G., Li, X. Sen, Zhang, K., Li, B., Zhang, Y., 2013. Effects of impermeable boundaries on gas production from hydrate accumulations in the shenhu area of the South China sea. Energies 6,

Jo u

4078–4096. https://doi.org/10.3390/en6084078 Liu, Z., Wei, H., Peng, L., Wei, C., Ning, F., 2017. An easy and efficient way to evaluate mechanical properties of gas hydrate-bearing sediments: the direct shear test. J. Pet. Sci. Eng. 149, 56–64. https://doi.org/10.1016/j.petrol.2016.09.040 Lizárraga, J.J., Buscarnera, G., Frattini, P., Crosta, G.B., 2016. Spatially distribu ted modeling of landslide triggering: An approach based on principles of unsaturated soil plasticity. Landslides Eng. Slopes. Exp. Theory Pract. 2, 1287–1294. Maslin, M., Mikkelsen, N., Vilela, C., Haq, B., 1998. Sea-level- and gas-hydrate-controlled catastrophic sediment failures of the Amazon Fan. Geology 26, 1107–1110. https://doi.org/10.1130/0091-7613(1998)026<1107:SLA GHC>2.3.CO;2 Maslin, M., Owen, M., Betts, R., Day, S., Jones, T.D., Ridgwell, A., 2010. Gas hydrates: Past and future geohazard? Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 368, 2369–2393. https://doi.org/10.1098/rsta.2010.0065 Masui, A., Haneda, H., Ogata, Y., Aoki, K., 2005. Effects of methane hydrate formation on shear strength of synthetic methane hydrate sediments, in: Proceedings o f the Fiftheenth (2005) International Offshore and Polar Engineering Conference. Seoul, Korea, pp. 364–369. 32

Journal Pre-proof Mestdagh, T., Poort, J., De Batist, M., 2017. The sensitivity of gas hydrate reservoirs to climate change: perspectives from a new combined model for permafrost-related and marine settings. Earth-Science Rev. 169, 104–131. https://doi.org/10.1016/j.earscirev.2017.04.013 Miyazaki, K., Masui, A., Sakamoto, Y., Aoki, K., Tenma, N., Yamaguchi, T., 2011. Triaxial compressive properties of artificial methane-hydrate-bearing sediment. J. Geophys. Res. Solid Earth 116, 1–11. https://doi.org/10.1029/2010JB008049 Moridis, G.J., 2014. User’s manual for the HYDRATE v1.5 option of TOUGH+ v1.5: A code for the simulation of system behavior in hydrate-bearing geologic media. Moridis, G.J., Pruess, K., 2014. User’s manual of the TOUGH + CORE code v1.5: a general-purpose simulator of non-isothermal flow and transport through porous and fractured media. Moridis, G.J., Queiruga, A.F., Reagan, M.T., 2018. Geomechanical Stability and Overall System Behavior of Sloping Oceanic Accumulations of Hydrates Responding to Dissociation Stimuli, in:

f

Offshore Technology Conference Asia. Kuala Lumpur, Malaysia.

oo

https://doi.org/10.4043/28496-MS

Moridis, G.J., Seol, Y., Kneafsey, T.J., 2005. Studies of reaction kinetics of methane hydrate

pr

dissocation in porous media, Lawrence Berkeley National Laboratory. Lawrence Berkeley National Laboratory.

Myshakin, E.M., Seol, Y., Lin, J.S., Uchida, S., Collett, T.S., Boswell, R., 2018. Numerical simulations

e-

of depressurization-induced gas production from an interbedded turbidite gas hydrate-bearing sedimentary section in the offshore India: Site NGHP-02-16 (Area-B). Mar. Pet. Geol.

Pr

https://doi.org/10.1016/j.marpetgeo.2018.10.047

Nixon, M.F., Grozic, J.L.H., 2007. Submarine slope failure due to gas hydrate dissociation: a preliminary quantification. Can. Geotech. J. 44, 314–325. https://doi.org/10.1139/t06-121

al

Paull, C.K., Ussler, I., Holbrook, W.S., 2007. Assessing methane release from the coloss al Storegga submarine landslide. Geophys. Res. Lett. 34, 1–5. https://doi.org/10.1029/2006GL028331

rn

Paull, C.K., Ussler, W., Dillon, W.P., 2011. Potential Role of Gas Hydrate Decomposition in Generating Submarine Slope Failures, in: Natural Gas Hydrate. pp. 149–156.

Jo u

https://doi.org/10.1007/978-94-011-4387-5_ 12 Phrampus, B.J., Hornbach, M.J., 2012. Recent changes to the Gulf Stream causing widespread gas hydrate destabilization. Nature 490, 527–530. https://doi.org/10.1038/nature11528 Qiu, K., Yamamoto, K., Birchwood, R., Chen, Y., 2015. Well integrity evaluation for methane hydrate production in the deepwater Nankai Trough. SPE Drill. Complet. 30, 52–67. Reagan, M.T., Moridis, G.J., 2008. Dynamic response of oceanic hydrate deposits to ocean temperature change. J. Geophys. Res. 113, C12023. https://doi.org/10.1029/2008JC004938 Ruppel, C.D., Kessler, J.D., 2017. The interaction of climate change and methane hydrates. Rev. Geophys. 55, 126–168. https://doi.org/10.1002/2016RG000534 Rutqvist, J., 2017. An overview of TOUGH-based geomechanics models. Comput. Geosci. 108, 56–63. https://doi.org/10.1016/j.cageo.2016.09.007 Salehabadi, M., Jin, M., Yang, J., Haghighi, H., Ahmed, R., Tohidi, B., 2008. Finite element modelling of casing in gas hydrate bearing sediments. Eur. Conf. Exhib. https://doi.org/https://doi.org/10.2118/113819-M S Skarke, A., Ruppel, C., Kodis, M., Brothers, D., Lobecker, E., 2014. Widespread methane leakage from the sea floor on the northern US Atlantic margin. Nat. Geosci. 7, 657–661. 33

Journal Pre-proof https://doi.org/10.1038/ngeo2232 Sloan, E.D., Koh, C., 2008. Clathrate hydrates of nautral gases, 3rd Editio. ed. Taylor and Francis, Inc., Boca Raton, FL. https://doi.org/https://doi.org/10.1201/9781420008494 Somerton, W.H., El-Shaarani, A.H., Mobarak, S.M., 1974a. High temperature behavior of rocks associated with geothermal type reservoirs, in: SPE California Regional Meeting. Society of Petroleum Engineers, San Francisco, California. https://doi.org/10.2118/4897-MS Somerton, W.H., Keese, J.A., Chu, S.L., 1974b. Thermal behavior of unconsolidated oil sands. Soc. Pet. Eng. J. 14, 513–521. https://doi.org/10.2118/4506-PA Stone, H.L., 2007. Probability model for estimating three-phase relative permeability. J. Pet. Technol. 22, 214–218. https://doi.org/10.2118/2116-pa Sultan, N., Cochonat, P., Foucher, J., Mienert, J., 2004. Effect of gas hydrates melting on seafloor slope instability. Mar. Geol. 213, 379–401. https://doi.org/10.1016/j.margeo.2004.10.015

f

Sun, J., Ning, F., Li, S., Zhang, K., Liu, T., Zhang, L., Jiang, G., Wu, N., 2015. Numerical simulation

oo

of gas production from hydrate-bearing sediments in the Shenhu area by depressurising: The effect of burden permeability. J. Unconv. Oil Gas Resour. 12, 23–33.

pr

https://doi.org/10.1016/j.juogr.2015.08.003

Talling, P., Clare, M., Urlaub, M., Pope, E., Hunt, J., Watt, S., 2014. Large submarine landslides on https://doi.org/10.5670/oceanog.2014.38

e-

continental slopes: geohazards, methane Release, and climate change. Oceanography 27, 32–45. Tang, L., Li, X., Feng, Z., Li, G., Fan, S., 2007. Control mechanisms for gas hydrate production by

Pr

depressurization in different scale hydrate reservoirs. Energy and Fuels 21, 227–233. Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.

al

Waite, W.F., Santamarina, J.C., Cortes, D.D., Dugan, B., Espinoza, D.N., Germaine, J., Jang, J., Jung, J.W., Kneafsey, T.J., Shin, H., Soga, K., Winters, W.J., Yun, T.-S., 2009. Physical properties of

rn

hydrate-bearing sediments. Rev. Geophys. 47, RG4003. https://doi.org/10.1029/2008RG000279.Table

Jo u

Wang, X., Collett, T.S., Lee, M.W., Yang, S., Guo, Y., Wu, S., 2014. Geological controls on the occurrence of gas hydrate from core, downhole log, and seismic data in the Shenhu area, South China Sea. Mar. Geol. 357, 272–292. https://doi.org/10.1016/j.margeo.2014.09.040 White, M.D., McGrail, B.P., 2006. STOMP-HYD: a new numerical simulator for analysis of methane hydrate production from geologic formations., in: Proceedings of 2nd International Sympo- Sium on Gas Hydrate Technology. Daejeon, Korea. Xiong, Z., Zhang, J., 2012. Effect of dissociation of gas hydrate on the stability of submarine slope, in: 31st ASME International Conference on Ocean, Offshore and Arctic Engineering. Rio de Janeiro, Brazil. Xu, W., Germanovich, L.N., 2006. Excess pore pressure resulting from methane hydrate dissociation in marine sediments: A theoretical approach. J. Geophys. Res. Solid Earth 111, B01104. https://doi.org/10.1029/2004JB003600 Yoneda, J., Masui, A., Konno, Y., Jin, Y., Egawa, K., Kida, M., Ito, T., Nagao, J., Tenma, N., 2015. Mechanical properties of hydrate-bearing turbidite reservoir in the first gas production test site of the Eastern Nankai Trough. Mar. Pet. Geol. 66, 471–486. https://doi.org/10.1016/j.marpetgeo.2015.02.029 34

Journal Pre-proof Yoneda, J., Masui, A., Konno, Y., Jin, Y., Kida, M., Katagiri, J., Nagao, J., Tenma, N., 2017. Pressure-core-based reservoir characterization for geomechanics: Insights from gas hydrate drilling during 2012–2013 at the eastern Nankai Trough. Mar. Pet. Geol. 86, 1–16. https://doi.org/10.1016/j.marpetgeo.2017.05.024 Zhang, G., Yang, S., Ming, Z., Liang, J., Jingan, L., Holland, M., Schultheiss, P., 2014. GMGS2 expedition investigates rich and complex gas hydrate environment in the South China Sea. Fire Ice 14, 1–5. Zhang, X.H., Lu, X.B., Chen, X.D., Zhang, L.M., Shi, Y.H., 2016. Mechanism of soil stratum instability induced by hydrate dissociation. Ocean Eng. 122, 74–83.

Jo u

rn

al

Pr

e-

pr

oo

f

https://doi.org/10.1016/j.oceaneng.2016.06.015

35

Journal Pre-proof

Highlights: A slope stability model is developed for slopes undergoing hydrate dissociation.



The transient pore pressure is computed from a THC coupled analysis.



Three failure modes depending on overburden characteristics are identified.

Jo u

rn

al

Pr

e-

pr

oo

f



36