Simulated geomechanical responses to marine methane hydrate recovery using horizontal wells in the Shenhu area, South China Sea

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Marine and Petroleum Geology journal homepage: www.elsevier.com/locate/marpetgeo

Research paper

Simulated geomechanical responses to marine methane hydrate recovery using horizontal wells in the Shenhu area, South China Sea Guangrong Jina, Hongwu Leib, Tianfu Xua,∗, Xin Xina, Yilong Yuana, Yingli Xiaa, Jinge Juoc a

Key Laboratory of Groundwater Resources and Environment, Ministry of Education, Jilin University, Changchun 130021, China State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, China Academy of Sciences, Wuhan, Hubei 430071, China c School of Civil and Environmental Engineering, Georgia Institute of Technology, 790 Atlantic Dr NW, Atlanta, GA 30313, USA b

A R T I C L E I N F O

A B S T R A C T

Keywords: Natural gas hydrate THM model Geomechanical response Depressurization Seafloor subsidence South China sea

The recovery of natural gas hydrate (NGH) from marine sediments faces challenging due to not only the gas productivity itself but also the possible geohazards, such as seafloor subsidence, submarine landslide, wellbore instability and possible sand production. The coupled thermal-hydrodynamic-mechanical (THM) processes during NGH recovery are generally complex, and numerical simulation tools are needed to assess related geomechanical responses. We have extended the Biot consolidation model in our previous simulator TOUGH2Biot, and incorporated into the existing TOUGH + hydrate code, resulting in a THM simulator for the NGH recovery. The THM simulator is used to assess the geomechanical responses to gas recovery from an unconfined hydratebearing sediment (HBS) in the Shenhu area, South China Sea. We investigated depressurization using constant bottom-hole pressure through a horizontal well. Results show that methane production quickly reaches a stabilized state and the water rate increases linearly. The drawdown of pore pressure around the well controls the increase in effective stress. Subsidence becomes significant after depressurization due to the quickly propagation of pore pressure. Depressurization in early stage could contribute to more than half of the total subsidence. The decreasing production pressure leads to an increase in the methane production rate but deterioration in subsidence. A decrease in intrinsic permeability of overlying and underlying layer is undesirable due to its decrease in methane production rate and the worse of subsidence. A balance between gas productivity and related geomechanical response must be achieved. The methods and preliminary results presented in this study could help us to understand the geomechanical behaviors during NGH recovery and to design trial production schemes under similar conditions.

1. Introduction Natural gas hydrate (NGH) is a crystalline solid compound formed by gas molecules (such as methane, ethane, propane, and carbon dioxide) that occupy the cage structures of water molecules (Sloan and Koh, 2007; Liu et al., 2015a, 2015b). NGH has been found in either permafrost regions or deep ocean sediments (Max et al., 2005; Max and Johnson, 2014) where the ambient conditions of high pressure and low temperature are existed. As the large reserve of trapped natural gas, NGH is considered as a potential unconventional fossil fuel resource (Collett, 2002), and its recovery attracts a wide attentions from academic institutions, governments and oil companies (Moridis et al., 2009; 2011a). Several methods, such as depressurization, thermal stimulation, inhibitor injection, are employed to destabilize hydrates by altering the conditions of pressure and/or temperature for NGH, and methane gas is

∗

then recovered through a production well (Jin et al., 2016; Zhang et al., 2010; Bhade and Phirani, 2015). Recent studies are conducted to reveal the mechanism, productivity and efficiency of these methods from laboratory-scale experiments to field trials (Kneafsey et al., 2007; Marinakis et al., 2015; McGrail et al., 2004; Yamamoto et al., 2014; Sun et al., 2014; Li et al., 2016; Moridis et al., 2011a). In 2013, the first offshore field trial was conducted in the Nankai Trough, Japan. Depressurization was technologically validated as feasible method and a total of 12000 m3 of methane gas was recovered during the 6-day of production (Yamamoto et al., 2014). In general, depressurization by reducing the pressure in NGH formation is regarded as the most economic method (Moridis et al., 2011a; Zhang et al., 2010). Thermal stimulation is more suitable as an assisted approach, due to the costly heating, to improve production performance (Zhang et al., 2010; Moridis and Kowalsky, 2005; Reagan et al., 2014). The effect of inhibitor weakens because of the dilution of seawater. Gas exchange

Corresponding author. E-mail address: [email protected] (T. Xu).

https://doi.org/10.1016/j.marpetgeo.2017.11.007 Received 31 December 2016; Received in revised form 28 July 2017; Accepted 6 November 2017 0264-8172/ © 2017 Published by Elsevier Ltd.

Please cite this article as: Jin, G., Marine and Petroleum Geology (2017), https://doi.org/10.1016/j.marpetgeo.2017.11.007

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

method is recently proposed by using CO2 and/or N2 to replace methane molecules in the cage structure (McGrail et al., 2004; White et al., 2011). However, the exchange rate highly depends on permeability. In addition to the feasibility and efficiency of production, NGH recovery also arouses the concern of mechanical response due to the decreased pore pressure in hydrate-bearing sediment (HBS) (Rutqvist and Moridis, 2008; Gupta et al., 2015). Furthermore, the solid hydrate dissociates into mobile water and methane gas after hydrate dissociation, and the sediment strength and stiffness of the HBS decrease accordingly. This may causes seafloor subsidence, deterioration of wellbore stability and possible sand production (Uchida et al., 2016). Therefore, recovery of subsurface NGH involves thermal (T), hydraulic (H), mechanical (M) and chemical (especially when gas exchange) processes. Several simulators of TH processes, such as TOUGH + HYDRATE, STOMP-HYD, MH21-HYDRES and HydrateResSim (Moridis et al., 2008; White et al., 2011; Gupta et al., 2015), have been developed to assess the productivity of various methods. However, geomechanical responses were not considered in these numerical models. This study will focus on the coupled THM processes without considering the chemical gas exchange method. To improve understanding of the geomechanical responses during NGH recovery, several simulators have been developed. TOUGH + HYDRATE code was incorporated with the commercial geomechanical code FLAC3D by Rutqvist and Moridis (2008). Klar et al. (2010) used the FLAC to assess the effect of hydrates on the stress-strain behavior. Zhou et al. (2014) investigated deformation during the Nankai production trial. Recently, a chemo-thermo-hydro-mechanical simulator using elasto-viscoplastic model was proposed to predict ground stability (Kimoto et al., 2010). Gupta et al. (2015) developed a hydro-geomechanical hydrate simulator by coupling fluid flow with linear elasticity. Stress and strain distributions were calculated for gas production from horizontal and vertical wells using the depressurization method (Rutqvist et al., 2012; Klar et al., 2010; Gupta et al., 2015). The effective stress increases in response to the depletion of pore pressure. Stress around well are driven by pressure depletion and this may reduce the integrity and stability of wellbore. Moreover, increased shear stress may cause shear failure in formation and grain detachment, which leads to sand production (Rutqvist et al., 2012; Uchida et al., 2016). Recently, a 2D thermo-hydro-mechanical bonded contact model was proposed to study the changes of macro-scale and micro-scale mechanical behaviors by the distinct element method code, PFC2D (Jiang et al., 2016). The number of coupled THM simulations during hydrate dissociation has increased in recent years. However, there is still a need to further explore the geochemical response. A powerful simulator is the basis to improve our understanding of geomechanical behaviors in a complex geological setting and/or when using various production strategies. Furthermore, the geomechanical studies for the Shenhu area of South China Sea are rarely reported although some tri-axial tests have been investigated (Zhang et al., 2015; Sun et al., 2017). This study takes module Biot of our previous THM simulator TOUGH2Biot (Lei et al., 2015) to characterize the mechanical behavior of HBS. The geomechanical module Biot is integrated into the existing TOUGH + hydrate code (Moridis et al., 2008), resulting in the specific THM simulator TOUGH + hydrate + Biot (simply called hydrateBiot) for the NGH recovery. The geomechanical response under depressurization through a horizontal well in Shenhu area is then investigated. The relationship between the geomechanical response and the gas productivity under different production pressure and permeability of bounded layer are discussed.

Table 1 Governing equations of fluid and heat flow in TOUGH + hydrate. Description

Equation

Mass and energy conservation

d ∫ M κdV = ∫ F κ•ndΓ + ∫ qκdV dt Vn Γn Vn M κ = ∑β = A, G, I , H φSβ ρβ Xβκ , κ =

Mass accumulation Mass flux (aqueous phase) Mass flux (Gas phase) Energy accumulation Heat flux Reaction heat of hydrate dissociation

FκA = −k

krA ρA κ XA (∇PA μA

− ρA g),

w, m, i, h

κ = w, m, i

b

slippage krG ρG κ κ κ ⎞ FG = −k ⎛1 + XG (∇PG − ρG g) + JG , PG ⎠ μG ⎝ θ M = (1 − φ) ρR CR T + ∑β = A, G, H , I φSβ ρβ Uβ + Qdiss

κ = w, m

F θ = −λ∇T + ∑β = A, G hβ Fβ Δ (φρH SH ΔUH ) for equilibrium dissociation Qdiss = ⎧ ⎨ for kinetic dissociation ⎩QH ΔUH

developed by the Lawrence Berkeley National Laboratory, provides a reliable and open source base to simulate the thermo-hydrological processes during the NGH recovery. TOUGH + hydrate can simulate the transport of four components (water, CH4, hydrate, water-soluble inhibitors such as salts or alcohols) among four phases (gas phase, liquid phase, ice phase and hydrate phase), and also the non-isothermal dissociation of hydrate and heat flow. The governing equations of mass and energy balances used in TOUGH + hydrate are summarized in Table 1 (see Nomenclature for definitions of all symbols used). To simulate the mechanical process, the Biot module of TOUGH2Biot (Lei et al., 2015) is extended to characterize the geomechanical response associated to the NGH recovery. Although the triaxial tests show that the stress-strain relationship of a hydrate-bearing sample is not expressed as elastic behavior (Zhang et al., 2015; Masui et al., 2008; Miyazaki et al., 2010, 2011), but the HBS can be regarded as elastic material as long as the range of application is sufficiently limited to small-strain cases far away from the critical state (Gupta et al., 2015). Based on the principle of effective stress (The effective stress, σ ′, is the difference between the total stress, σ , and the pore pressure, Pa , as σ ′ = σ − ξPa and the ξ is the Biot' coefficient) and the assumption of linear elasticity, an extended Biot consolidation model with displacement as the primary unknown variables was formulated as shown in Table 2 (Lei et al., 2015). The effect of temperature change on stress is also considered in the model. The strength increases linearly with hydrate saturation and also confining pressure, but the confining pressure has a minor effect. Therefore, the bulk modulus, shear modulus and cohesion are simply expressed as (Rutqvist et al., 2012):

K = (KSH1 − KSH0) × SH + KSH0

G = (GSH1 − GSH0) × SH + GSH0

c = (cSH1 − cSH0) × SH + cSH0 Where KSH0 and KSH1 are the bulk modulus without hydrate and when the hydrate saturation is 1, respectively. GSH0 and GSH1 are the shear modulus without hydrate and when the hydrate saturation is 1, respectively. cSH0 and cSH1 are the cohesion without hydrate and when the hydrate saturation is 1, respectively. 2.2. Coupling between TH and M

2. Modeling approaches The coupled THM processes can be decoupled into the fluid flow and heat flow models (TH) and the mechanical model. The hydrateBiot inherits the fully coupled TH processes from TOUGH + hydrate. Stress and strain is obtained by solving the extended Biot mechanical

2.1. Governing equations for THM processes The TOUGH + hydrate code (Moridis et al., 2008), which was 2

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Table 2 3D extended Biot mechanical model. Description displacement

Stress and strain

Governing equations

− G∇2 wx −

G ∂ 1 − 2υ ∂x

− G∇2 w y −

G ∂ 1 − 2υ ∂y

− G∇2 wz −

G ∂ 1 − 2υ ∂z

( σ ′ = 2G ( σ ′ = 2G ( z

τyz

∂w x ∂x

+

∂w x ∂x

+

∂w x ∂x

+

) ) )

∂w y ∂y ∂w y ∂y ∂w y ∂y

+

∂wz ∂z

+

∂wz ∂z

+

∂wz ∂z

+ εx + 3βT ΔT ⎫ ⎪ ⎪ + εy + 3βT ΔT ⎪ ⎬ υ ε + εz + 3βT ΔT ⎪ 1 − 2υ v ⎪ = Gγyz , τzx = Gγzx , τxy = Gγxy ⎪ ⎭

σx′ = 2G y

υ ε 1 − 2υ v

( ( (

υ ε 1 − 2υ v

equations. The TH processes are first calculated first and then sequentially the mechanical model is called in a single time step. Because the HBS is a mixture containing hydrate, sediment skeleton, aqueous phase and gas phase, they jointly bear the stress. Similar to Rutqvist et al. (2012) and Lei et al. (2015), the average pore pressure is calculated according to Pavg = SA PA + SG PG . Changes in effective stress give feedback to the hydraulic parameters through porosity and permeability as ′ ) and k = k 0 exp[b × (φ / φ0 − 1)] (the φ = φr + (φ0 − φr )exp(a × σM ′ = (σx′ + σy′ + σz′)/3), respectively. average effective stress is σM

)+ )+ )+

∂Pa ∂x

+ 3βT K

∂T ∂x

∂Pa ∂y

+ 3βT K

∂T ∂y

∂Pa ∂z

+ 3βT K

∂T ∂z

εx = − εy = − εz = −

∂w x , γyz ∂x ∂w y ∂y

=−

, γzx = −

∂wz , γxy ∂z

=−

=0 ⎫ ⎪ ⎪ =0 ⎬ ⎪ = γsat ⎪ ⎭

(

∂w y ∂z

+

∂wz ∂y

)⎫⎪

∂w x ∂z

)⎬

∂w y

) ⎪⎭

(

∂wz ∂x

+

(

∂w x ∂y

+

∂x

⎪ ⎪

initialization of mechanics includes reading the geometry, mechanical properties and initial and boundary conditions. The mechanical calculations are sequentially executed after TH is solved. The porosity and permeability is updated after stress calculation. 2.3.2. Simulator verification The THM code testing with a 1D analytical solution has been reported in Lei et al. (2015). The case of the gas recovery from HBS in the Gulf of Mexico (Rutqvist et al., 2012) was chosen to test the accuracy of hydrateBiot. The case was first studied by Rutqvist et al. (2012) using the link code of TOUGH + hydrate and FLAC3D (hereafter referred to as hydrateFLAC). The modeled domain using the hydrateBiot simulator involves the structure with a unit length (= 1 m) in y-axis and a formation thickness of 548.25 m (z-axis) and a lateral extension of 500 m (Fig. 2) in x-axis, which represents a system of parallel horizontal wells with 1000 m well spacing. A bottom hole pressure of 2.7 MPa is maintained at the well. The top of model is at a water depth of 2750 m, which is the seafloor depth. The HBS is about 18 m in thickness and is confined by shale. The main physical parameters are summarized in Table 3. The initial pore pressure distribution is hydrostatic, and temperature increases with depth according to geothermal gradient. Hydrate saturation is 0.7 throughout the HBS. The initial stress is assumed to be isotropic and to increases with depth below the seafloor with a bulk

2.3. Simulator development and verification 2.3.1. Simulator development The TH model is spatially discretized using an integral finite difference approach and the fully implicit method is used for time discretization. Further details are shown in Moridis et al. (2008). The mechanical model is implemented by the Galerkin finite element method as reported in Lei et al. (2015). The mechanical modular is developed using Fortran 90/95 and is integrated into the open source TOUGH + hydrate code, forming the improved hydrateBiot simulator. The modular structure of hydrateBiot is shown as Fig. 1. The

Fig. 1. Modular structure of hydrateBiot (the dotted line denotes modification).

Fig. 2. Model setup of the verification case (from Rutqvist et al., 2012).

3

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Table 3 Physical and mechanical parameters used in the verification case. Parameter

Value

Parameter

Value

Initial pressure at the base of HBS Initial temperature at the base of HBS Porosity of HBS and shale

33 MPa

95 MPa

0.7 0.1 μD

Bulk modulus at SH = 0 Bulk modulus at SH = 1 Shear modulus at SH = 0 Shear modulus at SH = 1 Cohesion at SH = 0 Cohesion at SH = 1

0.03

Friction angle

30°

Thermal expansion coefficient Biot's coefficient (−)

1 × 10−5 °C

Intrinsic permeability of HBS Hydrate saturation Intrinsic permeability of shale Water salinity (mass fraction)

21 °C 0.3 1 mD

575 MPa 87 MPa 525 MPa 0.5 MPa 2.0 MPa

−1

1

Fig. 4. Comparison of subsidence evolution at seafloor (at x = 0 m, y = 0.5 m and z = 0 m).

density of 2600 kg/m3. The model sides are no flow and heat exchange boundary because of symmetry. The top and bottom of model are set to be the constant temperature and pressure boundaries. The model sides and bottom are fixed in normal direction, and the top surface is freely moving. The geomechanical parameters in Table 3, are derived from experimental results by Masui et al. (2008) and Miyazaki et al. (2010, 2011). The certain mechanical properties (bulk and shear modulus and cohesion) increase linearly with hydrate saturation. The Biot' coefficient of the heavily cemented sedimentary rocks is less than one (Alam et al., 2010). However, it is close to one for the lightly cemented sediments (e.g. loose grainy sediments, sands or marine sediments). Thus the Biot' coefficient of 1 is used in this study. Pore pressure, temperature and hydrate saturation at point A (x = 0.5 m, z = 471 m) and B (x = 10 m, z = 471 m) were picked out (Fig. 3). Fig. 3 shows very similar variation between hydrateBiot and hydrateFLAC. Pore pressure drops quickly by pumping water at the production well and the pressure driving force causes hydrate dissociation. The decrease of hydrate saturation is slower than temperature because of the endothermic nature of hydrate dissociation. After the hydrate dissociates completely, temperature rises due to the heat convection of fluid flow and heat exchange with the overlying shale. As shown in Fig. 4, the vertical subsidence at seafloor increases with continuing production. The seafloor subsidences by the two simulators show similar trend and fit well. Although additional subsidence could be obtained from the Mohr-Coulomb stress-strain relationship, the elastic relationship could provide us a rough estimate. The elastic-

plastic relationship will be developed in subsequent version of hydrateBiot. 3. Problem setup 3.1. Site description The Shenhu area in the north slope of South China Sea is the first exploration place of the NGH reservoir in China (Fig. 5). It is located in Pearl River Mouth Basin between the Xisha Trough and Dongsha Islands (Wu et al., 2011; Li et al., 2011). The drilling campaign by China Geological Survey in 2007 discovered the sediments to be rich in NGH in the SH2, SH3 and SH7 drill holes. The thickness of the hydrate reservoir is estimated to be 10–44 m, which overlies the base of the gas hydrate stability zone (BGHSZ) (Wu et al., 2011). According to the lithology logs in the SH2 drill hole, the hydrate reservoir sits at about 185 m below sea floor whit a thickness of 44 m. The sea floor is about 1235 m below the sea level. The hydrate disseminates in the sediments that is primarily composed of clays (Zhang et al., 2010), with a porosity from 0.33 to 0.48 and a hydrate saturation from 25% to 48% (Li et al., 2011; Zhang et al., 2010; Wu et al., 2011; Huang et al., 2015). The hydrate reservoir is overlaid and underlaid by permeable strata, which have the same lithology as that in the hydrate reservoir but lack hydrate. The gases contained in the hydrate in the

Fig. 3. Comparison of the evolution of (a) pore pressure, (b) temperature and (c) hydrate saturation at the observed point A and B.

4

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Fig. 5. Location of the study area in South China Sea and the drilling sites (Wu et al., 2011).

3.3. Initial and boundary conditions

Shenhu area primarily consist of methane (96.10–99.91%) with minor ethane and propane presence (Wu et al., 2011).

The initial pore pressure at the base of HBS is about 14.97 MPa, and pore pressure increases hydrostatically with depth. Initial temperature at the base of HBS is 14.87 °C, and temperature increases with depth with a geothermal gradient of 0.047 °C/m. The constant temperature and fluid pressure are assumed on the top and bottom of the model, because the overlying and underlying layers are sufficient for the thermal and flow process in HBS (Jin et al., 2016). No mass and heat flux boundaries are applied at x = −10 km and x = 10 km as the domain laterally extends enough to insignificantly affect modeling results. It is assumed that the initial stress distribution is isotropic because of the lack of field data. Initial stress increases with depth and follows the stress gradient of approximate 22.6 kPa/m. The zero displacement condition is assumed at the direction normal to the model sides and bottom. The surface is freely moving.

3.2. Conceptual model According to the logs in SH2 site, the conceptual model for numerical simulation consists of three horizontal layers including the HBS with a thickness of 44 m and overlying and underlying layers with a thickness of 185 m, respectively (Fig. 6a). These three layers are considered as homogeneous and isotropic media. One horizontal well is employed to recover the gas. The horizontal well is preferred because it allows the depressurization over a wider range in the unconfined HBS (Zhang et al., 2010; Moridis et al., 2011b; Jin et al., 2016). The horizontal well locates at a depth of 208 m, which is near the middle part of HBS. The modeled domain laterally extends to 10 km in x direction (Fig. 6). To reduce the computational resources consumption for this 3D geological model, the pressure and temperature drop in the wellbores of the horizontal well are neglected. This approach simplifies the modeled domain into two-dimensional domain (x-z plane) with a unit width (1 m) in y direction. The size of modeled domain is 20 km × 414 m, which is discretized into 10296 cells (132 × 78) for simulation. The grid is horizontally refined near the well with the minimum interval of 0.1 m. The grid size increases logarithmically with the distance to the well. Vertical grid size in HBS is 1.0 m and it gradually increases in both overlying and underlying layers, and reaches to the maximum size of 50 m at the top and bottom of the model.

3.4. Multiphase flow and thermal parameters Related hydraulic and thermal parameters are summarized in Table 4. The porosity in the sediment is 0.38, and the intrinsic permeability is 10 milidracy (mD) (Zhang et al., 2010; Su et al., 2012). It is assumed that only methane gas molecules are contained in hydrate. Hydrate saturation (SH) in the HBS is set based on the core analysis (Wu et al., 2011). NGH is heavily accumulated in the central part of the HBS, and SH decreases towards both shallow and deep parts (Fig. 7). The aqueous saturation of HBS can be inferred by 1-SH. The layers overlying and underlying the HBS are saturated by water with aqueous saturation 5

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Fig. 6. (a) 3D view of the position of hydratebearing sediment and the position of horizontal well, and (b) cross-sectional view of the well position.

of 1.0. The relative aqueous and gas permeability varies with the aqueous and gas saturation, respectively, which can be inferred from:

Table 4 Thermal-physical parameters of the formation. Parameter Porosity Hydrate density Hydrate thermal conductivity Hydrate specific heat Rock grain density Rock grain specific heat Wet thermal conductivity Dry thermal conductivity

Value 0.38

920.0 kg/m

3

0.50 W/m/K

2.1 kJ/kg/K

2600 kg/m3 1.0 kJ/kg/°C 3.0 W/m/K

Parameter

Value

Intrinsic permeability Irreducible aqueous saturation (SirA ) Index for aqueous phase (nA ) Irreducible gas saturation (SirG ) Index for gas phase (nG ) Index for pore structure (m ) The entry capillary pressure (P0 )

10 mD 0.30

krA = (SA∗ )nA , SA∗ = (SA − SirA)/(1 − SirA) krG = (SG∗)nG , SG∗ = (SG − SirG )/(1 − SirG )

5.0

where krA and krG is the relative permeability of aqueous phase and gas phase, respectively, SirA is irreducible aqueous saturation, nA is index for aqueous phase, SirG is irreducible gas saturation and nG is index for gas phase. The capillary pressure in the model is a function of saturation, which can be inferred from (van Genuchten, 1980):

0.03 3.5 0.45 1.0 × 105 Pa

1.0 W/m/K

Pcap = −P0 ([SA∗ ]−1/ m − 1)1 − m , SA∗ = (SA − SirA)/(1 − SirA) Where SirA is irreducible aqueous saturation and m is index for pore structure, P0 is the entry capillary pressure.

3.5. Mechanical parameters The tri-axial test data (Fig. 8) show that the strength of the artificial samples containing NGH (Zhang et al., 2015), which is similar to that of HBS from South China Sea, is lower than the sand sample from either the artificial or natural core sample (Masui et al., 2008). The strength increases with hydrate saturation and confining pressure and the effect of confining pressure on strength is ignored. The secant modulus can be calculated from test data and it could be served as the elastic modulus. The mechanical properties (the bulk modulus and shear modulus) are derived as listed in Table 5. The artificial samples here could provide us a rough limitation of the mechanical properties because of the lack of test data on natural core sample in South China Sea. The friction is assumed to be independent of hydrate saturation and equal to 5° as results of Zhang et al. (2015).

Fig. 7. The hydrate saturation profile inferred from core analysis (Wu et al., 2011).

6

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

1.12 × 104 m3 and 1.86 × 103 m3, respectively. 4.2. Evolution of fluid and physical parameters Fig. 10 shows the evolution of the distribution of pore pressure drawdown induced by depressurization. The largest pore pressure drawdown is −6.9 MPa at the production well, which responds to the initial pore pressure of 14.9 MPa at the depth of production well. The pore pressure drawdown near the production well is round in shape. This pattern of pore pressure drawdown closely follows that of a typical point sink in a homogeneous and isotropic formation. The pore pressure decreases by about −1.5 MPa in the range of 20 m from well. Thus, pore pressure sharply decreases from −6.9 MPa to −1.5 MPa within the short distance of 20 m from well. The pore pressure then decrease slowly when measured farther away from well. The pressure driving force consequently causes the roundness of the dissociation front of the hydrate around the point sink (Fig. 11a), and its endothermic nature induces the decrease in temperature (Fig. 12a). The hydrates at the top and bottom of the HBS also dissociate (Fig. 11a), and the temperature decreases accordingly (Fig. 12a). However, the gaseous methane accumulate more at the bottom of the HBS than that at the top because of the high-temperature of the initial water at the bottom of the HBS, and water flow from the underlying layer (Fig. 13). As depressurization proceeds, the range of pore pressure drawdown expands slowly. The range of the dissociation front of the hydrate expands similarly but the propagation of the hydrate dissociation front lags behind that of the pore pressure drawdown. Furthermore, hydrates at the top and bottom of the HBS where hydrate content is initially low is fully dissociated earlier than that in middle part (Fig. 11). Thus the freeing up pore space is conductive to fluid flow. The flow of cold water in the overlying layer and warm water in underlying layer flow toward the well results in the vertical profile of temperature change, and they support the dissociation of hydrate (Fig. 12). The presence of hydrate decreases the effective permeability, therefore, water flows directly from the overlying and underlying layers, or preferentially flows along the interface of the hydrate dissociation front (Fig. 11).

Fig. 8. Strength of HBS samples similar to that from South China Sea (Zhang et al., 2015) and the test date of sand samples (σc is confining pressure).

Table 5 Mechanical parameters derived. Parameter

Value

Bulk modulus at SH = 0 Bulk modulus at SH = 1 Shear modulus at SH = 0 Shear modulus at SH = 1 Friction angle Cohesion at SH = 0 Cohesion at SH = 1 Thermal expansion coefficient Biot's coefficient (−)

2.0 MPa 32 MPa 1.5 MPa 29 MPa 5° 0.1 MPa 2.0 MPa 1 × 10−5 °C 1

−1

3.6. Recovery scenario 4.3. Stress and displacement evolution As an economic method, depressurization has been tested in the marine HBS in the Nankai Trough of Japan. Depressurization through the horizontal well is more useful to increase gas recovery because of the increasing contact area between the well and the HBS (Zhang et al., 2010; Moridis et al., 2011b). Therefore, depressurization is employed to recover gas. The production pressure is 8 MPa, which is considered as a mild depressurization, and is maintained at a constant level during the production. Since gas production is affected by not only the physical properties (such as porosity and permeability) but also operational factors (such as production pressure, well spacing and perforation intervals) (Jin et al., 2016; Sun et al., 2016), the effects of production pressure and permeability of the bounded layers on the geomechanical response are also investigated.

The stress change induced by depressurization after 2 year is shown in Fig. 14. The effective stress increases in response to depressurization. The largest change in stress is located at the depth of the production well (Fig. 14) because of the largest drawdown of pore pressure. The maximum and minimum principal stress changes for the total and effective stress are approximately equal in value, which is the result of the assumption of the homogeneous medium, the initially isotropic stress distribution and the presence of permeable layers. Sediment strength decreases due to the hydrate dissociation of HBS (Fig. 11) and this may cause the slight difference in changes in the maximum and minimum principal stress (Fig. 14). Two positions at distances of 0.6 m and 9.3 m from the production well (but the same depth) are respectively selected to discuss the evolution of effective principal stress (Fig. 15a). When depressurization begins, the effective principal stress first increases sharply, and then reaches a relative stable state because of the quick propagation of pore pressure. With the dissociation of hydrate at top and bottom of HBS (Fig. 11), effective permeability mainly depending on hydrate saturation increases, and water in overlying and underlying layer flows into HBS in an easy way. The domain of pore pressure drawdown decreases accordingly and pore pressure drawdown is weakened. Therefore, the effective principal stress at 0.6 m flattens out after 50 days. The effective principal stress at 9.3 m from well is lower than that at 0.6 m because of the pore pressure loss along the distance away from the well. The shear stress, which is proportional to the difference between the maximum and minimum effective principal stresses, increases because of the lateral limitation of contraction. When comparing the evolution

4. Results and discussion 4.1. Gas and water production As shown in Fig. 9a, there is significant methane release from hydrate in the initial 73 days. The highest release rate (QR) reaches 40 m3/ day, however, the released methane cannot be totally extracted from the production well because of the capillary pressure. Methane extracted from well (QP) is low than QR over the 73 days. Due to quick propagation and low storability, pore pressure around the production well decreases and then stabilizes. As a result, QR decreases sharply in the initial 73 days, then both QR and QP decrease slowly. The water discharge rate (QW) increases because of the presence of unconfined layers. The cumulative methane and water production after 2 year is 7

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Fig. 9. Evolution of (a) gas production rate (QP), release rate (QR) and water production rate (QW) and (b) cumulative gas volume (VP) and water volume (VW) using depressurization method. The unit of the volumetric rate is measured at standard temperature and pressure.

Fig. 11. Evolution of hydrate saturation (SH) distributions at different times.

et al., 2012). This may suggest that the prediction of HBS failure heavily depends on sediment texture, and additional tri-axial test data are need to support its reliability. The corresponding vertical displacement at the seafloor is shown in Fig. 16. The seafloor level drops sharply in response to the initial dramatic drop in pore pressure at the well, and then decreases slowly. The estimated seafloor subsidence in the first 50 days reaches 1.52 m. It is more than half of the subsidence that occurs during the 2 years of production (Fig. 16a). The seafloor subsides about 1.8 m and 2.5 m after 0.5 year and 2 year (Fig. 16b), respectively. Due to the low strength of the unconsolidated sediment, the range of seafloor subsidence is relatively large in horizontal direction. This may implies that a minor pressure decrease in production well would cause subsidence at seafloor. However, since the mechanical behavior containing hydrate is complex, related tri-axial experiments for sediment strength and constitutive relationships need to be further investigated and developed, respectively, to improve our understanding of the geomechanical responses.

Fig. 10. Evolution of pore pressure drawdowns (△P) at different times (the negative value indicate the decrease in pore pressure relative to background condition).

of maximum and minimum effective principal stress to the Mohr-Coulomb (MC) failure criterion, the effective principal stress path at 0.6 m from well is very close to the MC failure surface estimated from tri-axial test data using silt clay (Zhang et al., 2015), which is similar to the HBS sample from South China Sea. This implies that yield failure may quickly occurs in vicinity the well, and sand grains may be detached and flow toward the production well. The effective principal stress path at 9.3 m from the well also reaches the MC failure surface estimated from Zhang et al. (2015) after 2 year. Furthermore, the effective principal stress path is different from the MC failure surface estimated from a sand rich formation, which is similar to samples from Nankai Trough (Masui et al., 2008; Miyazaki et al., 2011; Yoneda et al., 2015; Rutqvist

4.4. Effect of production pressure The geomechanical responses to production pressure (PW), a major parameter that affects gas productivity, are investigated here. The production pressure of 10 MPa corresponds to a mild pore pressure drawdown of about −5 MPa at the well. This pressure causes a low QP of 10 m3/day (Fig. 17a). Despite of the low QP, the seafloor subsides 8

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Fig. 14. Vertical profile of changes in total and effective principal stress after 2 years (along wellbore). The vertical profile is extracted along the intersecting line of the plane x = 100 m and y = 0.5 m.

range of pore pressure drawdown and effective stress (Fig. 18). Consequently, worse seafloor subsidence is obtained. As shown in Fig. 17a, the QP under different PW decrease with time and tend to be consistent at later period. Although a high methane production rate is obtained in the early period, the subsidence of 5 MPa in the initial 50 days could reaches the subsidence of 8 MPa after 730 days. This quick subsidence may induce wellbore instability and tilt of the production platform. Therefore, the gas productivity and related geochemical response should be balanced during NGH recovery. Fig. 12. Evolution of temperature changes (△T) at different times (the negative value indicate temperature decrease).

4.5. Effect of intrinsic permeability of bounded layers Because intrinsic permeability affects fluid flow and the propagation of pore pressure, the intrinsic permeability of the overlying and underling layers (OL and UL, respectively) are simultaneously reduced to investigate the effect of intrinsic permeability of the bonded layers on QP and seafloor subsidence. QP decreases with intrinsic permeability of bonded layers and greater (worse) seafloor subsidence is obtained (Fig. 19). The decreasing intrinsic permeability of the bonded layers implies the decrease in penetrability within the OL and UL. The water supply from surrounding formation (i.e. OL and UL) decreases with the decreasing intrinsic permeability of OL and UL. Although the effective permeability of HBS is initially low because of the presence of hydrate, the water are gradually tend to flow from the HBS. Water from HBS increases and a larger domain (in HBS) for water flow is needed to maintain the constant pressure at production well. As a consequence, the worse mobility for fluid flow within the OL and UL and the low storability of water contributes the expanding range of the pore pressure drawdown (Figs. 10a and 20a, b). The wider range of pore pressure drawdown in the HBS indicates the increasing pressure driving force in the HBS for hydrate dissociation. Unfortunately, this increasing pressure driving force does not improve the methane production (Fig. 19a). This result is primarily caused by the fact that its endothermic nature restricts the hydrate dissociation and the low level of water from bounded layers. The largest change in temperature is limited in the middle part of the HBS (Fig. 20 c and d). For an intrinsic permeability of 5 mD, cold water in the OL and warm water in the UL could flows into the HBS (Fig. 20c). This effect is supported by the change in temperature in vertical direction, and it partly provide heat for hydrate dissociation. For an intrinsic permeability of 1 mD, however, water flows difficultly into the HBS and the range of temperature reduction is

Fig. 13. Evolution of gas saturation (SG) distributions at different times.

sharply after the initially dramatic drop of pore pressure at the well and then quickly subsides to a very low rate (Fig. 17b). This is caused by the quick propagation of pore pressure out the production well and the limited influence range under these low depressurization conditions. When decreasing the production pressure, the methane production rate increases but the seafloor subsidence increases accordingly. This is caused by the fact that a lower production pressure at the well causes not only the increase in effective stress but also the larger influence 9

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Fig. 15. (a) Effective principal stress evolution and corresponding (b) effective principal stress path at the positions of 0.6 and 9.3 m horizontally away from well. σ1′ and σ3′ are the maximum and minimum effective principal stress, respectively. The Mohr-Coulomb Failure Criterion is the critical line to determine the failure of sediment.

Fig. 16. (a) Temporal evolution of vertical displacement at the seafloor (at x = 100 m, y = 0.5 m and z = 0 m) and corresponding (b) spatial distribution (along the intersecting line of the plane y = 0.5 m and z = 0 m) at different recording time.

homogeneous and isotropic formation and isotropic initial stress distribution. This results in the fact that the seafloor subsidence worsens with the decrease in the intrinsic permeability of bounded layer (Fig. 19b). Moreover, the seafloor greatly subsides in early stage of production and the final subsidence is extremely high for an intrinsic permeability of 1 mD. This mainly caused by the larger domain of pore pressure drawdown. The subsidence mainly depends on the strength of sediment, the reduced degree and the domain of pore pressure, and the

largely limited in the HBS (Fig. 20d). Heat for hydrate dissociation may come from heat conduction from surrounding environment, and thus hydrate dissociation is restrained (Fig. 20f). Therefore, heat convection by mobile water is significant for better performance of the methane production process. The increasing range of pore pressure drawdown causes the increase in effective stress (Fig. 21). Likewise, the maximum and minimum effective principal stress are roughly equal due to the assumption of the

Fig. 17. Evolution of (a) gas production rate (QP), and (b) vertical displacement at the seafloor (at x = 100 m, y = 0.5 m and z = 0 m) under different production pressure (PW).

10

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

to geomechanical analysis for gas production from hydrate. 2. For depressurization in a horizontal well system, the gas production quickly reach a stable state. The increase in effective stress is controlled by pore pressure drawdown around the well. The seafloor quickly subsides after depressurization due to the quick propagation of the pore pressure. The seafloor subsidence in early stages contributes more than half of total subsidence. 3. The seafloor subsidence worsens (increases) with decrease in production pressure although corresponding gas production increases. The seafloor subsidence in the initial period for a low production pressure may reach the total subsidence of a high production pressure. A decrease in the intrinsic permeability of bounded layers is undesirable due to the decrease in the methane production rate and the wider range of pore pressure drawdown. The corresponding increase in effective stress results worse (increased) subsidence. Therefore, the gas productivity and related geochemical response should be balanced. Although preliminary geomechanical analysis for the Shenhu area is obtained, the elastic-plastic relationship for stress-strain is needed further development to analyzing the complex geomechanical behavior of the HBS formations. The range of problems concerning the processes in gas production from HBS is very broad. The present modeling results are specific to the conditions and parameters considered. The “numerical experiments” do give a detailed understanding of the dynamic evolution, and provide useful insight into multiphase fluid flow, heat transfer and geomechanical responses.

Fig. 18. Vertical profile of changes in the maximum and minimum effective principal stress ( Δσ1′ and Δσ3′, respectively) after 2 years. The vertical profile is extracted along the intersecting line of the plane x = 100 m and y = 0.5 m. The PW is production pressure.

duration of production. For the HBS is similar to the low hydrate content in the top and bottom of HBS but high in the middle part, therefore, decrease in intrinsic permeability is seemingly unfavourable for methane production and controlling subsidence.

Acknowledgments This work is jointly supported by the National Key Research and Development Program of China (2017YFC0307304), the National Natural Science Foundation of China (41602255), and China Geological Survey Project (121201005000151216).

5. Conclusions This study present an coupled thermal-hydrodynamic-mechanical simulator hydrateBiot, which is obtained by incorporating Biot module of TOUGH2Biot into the existing TOUGH + hydrate code, for simulating geomechanical processes during NGH recovery. It is then applied to evaluate the geomechanical response for methane production in an unconfined hydrate-bearing sediment in the Shenhu area, South China Sea. The following conclusions can be drawn from this study:

Nomenclature

M Fκ q V Γ t ϕ Sβ ρβ , ρR

1. A sequential coupling scheme is adopted between thermal-hydrodynamic process and mechanical process. The effect of hydrate content on sediment strength is considered in the model. The verification show that hydrateBiot is a reliable tool that can be applied

mass or energy accumulation [kg/m3 or J/m3] flux of mass or energy of component κ [kg/(m2· s)] sink/source volume [m3] surface area [m2] time [s] porosity [-] saturation of phase β [-] density of phase β or rock grain [kg/m3]

Fig. 19. Evolution of (a) gas production rate (QP), and (b) vertical displacement at the seafloor (at x = 100 m, y = 0.5 m and z = 0 m) under different intrinsic permeability of overlying layer (OL) and underlying layer (UL).

11

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Fig. 20. The distribution of (a) and (b) pore pressure drawdown, (c) and (d) temperature change, (e) and (f) hydrate saturation at 730 days for the intrinsic permeability of 5 mD and 1 mD corresponding to the overlying and underlying layer (OL and UL).

CR specific heat of rock grain [J/(kg·°C)] internal energy of phase β [J/kg] Uβ λ average thermal conductivity of grid [W/(K·m)] T temperature [°C] specific enthalpy of phase β [J/kg] hβ mass change of hydrate component under kinetic dissociation QH ΔUH specific enthalpy of hydrate dissociation/formation Superscripts κ component, κ = w, i, g is water, salt and gas, respectively Subscripts β

G K ν wl Pa βT σl′ εl τl′ γl γsat c φ ϕ0 ϕr k0 a, b

Fig. 21. Vertical profile of changes in the maximum and minimum effective principal stress ( Δσ1′ and Δσ3′, respectively) under different intrinsic permeability of overlying layer and underlying layer (OL and UL). The vertical profile is extracted along the intersecting line of the plane x = 100 m and y = 0.5 m.

X βκ k krβ μβ Pβ g bslippage J κβ

mass fraction of componentκ in phase β [-] permeability [m2] relative permeability of phaseβ [-] viscosity of phaseβ [Pa·s] pressure of phase β [Pa] gravitational acceleration vector Klinkenberg b factor accounting for gas slippage effects mass diffusion of componentκ in phase β [kg/(m2·s)]

Phase, β = A, G, H , I is aqueous, gas, hydrae and ice phase, respectively shear modulus [Pa] bulk modulus [Pa] Poisson ratio [-] displacement, l = x , y, z [m] Averaged pressure [Pa] Thermal expansion [1/K] effective stress l = x , y, z [Pa] normal strain l = x , y, z [-] shear stress l = xy, yz , zx [Pa] shear strain l = xy, yz , zx [-] Rock weight with saturated fluid kg/(m2·s2) cohesion [Pa] internal friction angle [°] porosity at zero stress [-] residual porosity at “infinite” stress [-] permeability at zero stress [m2] experimental coefficient for change of ϕ andk

References Alam, M.M., Borre, M.K., Fabricius, I.L., Hedegaard, K., Rogen, B., Hossain, Z., Krogsboll, A., 2010. Biot's coefficient as an indicator of strength and porosity reduction:

12

Marine and Petroleum Geology xxx (xxxx) xxx–xxx

G. Jin et al.

Moridis, G.J., Collett, T.S., Boswell, R., Kurihara, M., Reagan, M.T., Koh, C., Sloan, E.D., 2009. Toward production from gas hydrates: current status, assessment of resources, and simulation-based evaluation of technology and potential. SPE Reserv. Eval. Eng. 12 (5), 745–771. Moridis, G.J., Kowalsky, M.B., Pruess, K., 2008. TOUGH + HYDRATE V1. 0 User's Manual: a Code for the Simulation of System Behavior in Hydrate-bearing Geologic Media. Report LBNL-0149E. Lawrence Berkeley National Laboratory, Berkeley, CA, USA. Moridis, G.J., Kowalsky, M., 2005. Gas production from unconfined Class 2 hydrate accumulations in the oceanic subsurface. Chapter 7. In: Max, M., Johnson, A.H., Dillon, W.P., Collett, T. (Eds.), Economic Geology of Natural Gas Hydrates. Kluwer Academic/Plenum Publishers, pp. 249–266. Reagan, M.T., Moridis, G.J., Johnson, J.N., Pan, L., Freeman, C.M., Boyle, K.L., Keen, N.D., Husebo, J., 2014. Field-scale simulation of production from oceanic gas hydrate deposits. Transp. Porous Media 108 (1), 151–169. Rutqvist, J., Moridis, G.J., 2008. Development of a numerical simulator for analyzing the geomechanical performance of hydrate-bearing sediments. In: Proceedings of the 42nd U.S. Rock Mechanics Symposium. San Francisco American Rock Mechanics Association, ARMA Paper No. 139. Rutqvist, J., Moridis, G.J., Grover, T., Silpngarmlert, S., Collett, T.S., Holdich, S.A., 2012. Coupled multiphase fluid flow and wellbore stability analysis associated with gas production from oceanic hydrate-bearing sediments. J. Pet. Sci. Eng. 92–93, 65–81. Sloan, E.D., Koh, C.A., 2007. Clathrate Hydrates of Natural Gases, third ed. CRC Press, Taylor and Francis. Su, Z., Moridis, G.J., Zhang, K., Wu, N., 2012. A huff-and-puff production of gas hydrate deposits in Shenhu area of South China Sea through a vertical well. J. Pet. Sci. Eng. 86, 54–61. Sun, J., Ning, F., Zhang, L., Liu, T., Peng, L., Liu, Z., Li, C., Jiang, G., 2016. Numerical simulation on gas production from hydrate reservoir at the 1st offshore test site in the eastern Nankai Trough. J. Nat. Gas. Sci. Eng. 30, 64–76. Sun, J., Zhang, L., Ning, F., Lei, H., Liu, T., Hu, G., Lu, H., Lu, J., Liu, C., Jiang, G., Liang, J., Wu, N., 2017. Production potential and stability of hydrate-bearing sediments 1 at the site GMGS3-W19 in the South China Sea: a preliminary feasibility study. Mar. Pet. Geol. 86, 447–473. Sun, Y., Li, B., Guo, W., Lü, X., Zhang, Y., Li, K., Wang, P., Jin, G., Qu, L., 2014. Comparative analysis of the production trial and numerical simulations of gas production from multilayer hydrate deposits in the Qilian Mountain permafrost. J. Nat. Gas. Sci. Eng. 21, 456–466. Uchida, S., Klar, A., Yamamoto, K., 2016. Sand production model in gas hydrate-bearing sediments. Int. J. Rock Mech. Min. 86, 303–316. van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. White, M.D., Wurstner, S.K., McGrail, B.P., 2011. Numerical studies of methane production from Class 1 gas hydrate accumulations enhanced with carbon dioxide injection. Mar. Pet. Geol. 28 (2), 546–560. Wu, N., Zhang, H., Yang, S., Zhang, G., Liang, J., Lu, J.A., Su, X., Schultheiss, P., Holland, M., Zhu, Y., 2011. Gas hydrate system of Shenhu area, northern South China Sea: geochemical results. J. Geol. Res. http://dx.doi.org/10.1155/2011/370298. Article ID 370298. Yamamoto, K., Terao, Y., Fujii, T., Ikawa, T., Seki, M., Matsuzawa, M., Kanno, T., 2014. Operational overview of the first offshore production test of methane hydrates in the Eastern Nankai Trough. Offshore Technol. Conf. http://dx.doi.org/10.4043/ 25243-MS. Yoneda, J., Masui, A., Konno, Y., Jin, Y., Egawa, K., Kida, M., Ito, T., Nagao, J., Tenma, N., 2015. Mechanical behavior of hydrate-bearing pressure-core sediments visualized under triaxial compression. Mar. Pet. Geol. 66, 451–459. Zhang, K., Moridis, G.J., Wu, N., Li, X., Reagan, M.T., 2010. Evaluation of Alternative Horizontal Well Designs for Gas Production from Hydrate Deposits in the Shenhu Area, South China Sea. International Oil and Gas Conference and Exhibition in China. SPE, Beijing, China, pp. 8–10 131151. Zhang, X.H., Lu, X.B., Shi, Y.H., Xia, Z., 2015. Study on the mechanical properties of hydrate-bearing silty clay. Mar. Pet. Geol. 67, 72–80. Zhou, M., Soga, K., Xu, E., Uchida, S., Yamamoto, K., 2014. Numerical study on eastern Nankai Trough gas hydrate production test. Offshore Technol. Conf. http://dx.doi. org/10.4043/25169-MS.

calcareous sediments from Kerguelen Plateau. J. Petrol. Sci. Eng. 70 (3), 282–297. Bhade, P., Phirani, J., 2015. Gas production from layered methane hydrate reservoirs. Energy 82, 686–696. Collett, T.S., 2002. Energy resource potential of natural gas hydrates. AAPG Bull. 86 (11), 1971–1992. Gupta, S., Helmig, R., Wohlmuth, B., 2015. Non-isothermal, multi-phase, multi-component flows through deformable methane hydrate reservoirs. Comput. Geosci. 19, 1063–1088. Huang, L., Su, Z., Wu, N.Y., 2015. Evaluation on the gas production potential of different lithological hydrate accumulations in marine environment. Energy 91, 782–798. Jiang, M., Fu, C., Cui, L., Shen, Z., Zhu, F., 2016. DEM simulations of methane hydrate exploitation by thermal recovery and depressurization methods. Comput. Geotech. 80, 410–426. Jin, G., Xu, T., Xin, X., Wei, M., Liu, C., 2016. Numerical evaluation of the methane production from unconfined gas hydrate-bearing sediment by thermal stimulation and depressurization in Shenhu area, South China Sea. J. Nat. Gas. Sci. Eng. 33, 497–508. Kimoto, S., Oka, F., Fushita, T., 2010. A chemo–thermo–mechanically coupled analysis of ground deformation induced by gas hydrate dissociation. Int. J. Mech. Sci. 52, 365–376. Klar, A., Soga, K., Ng, M.Y.A., 2010. Coupled deformation-flow analysis for methane hydrate extraction. Geotechnique 60 (10), 765–776. Kneafsey, T.J., Tomutsa, L., Moridis, G.J., Seol, Y., Freifeld, B.M., Taylor, C.E., Gupta, A., 2007. Methane hydrate formation and dissociation in a partially saturated core-scale sand sample. J. Pet. Sci. Eng. 56 (1), 108–126. Lei, H., Xu, T., Jin, G., 2015. TOUGH2Biot - a simulator for coupled thermal- hydrodynamic-mechanical processes in subsurface flow systems. Comput. Geosci. 77, 8–19. Li, G., Moridis, G.J., Zhang, K., Li, X.S., 2011. The use of huff and puff method in a single horizontal well in gas production from marine gas hydrate deposits in the Shenhu Area of South China Sea. J. Pet. Sci. Eng. 77 (1), 49–68. Li, X.S., Xu, C.G., Zhang, Y., Ruan, X.K., Li, G., Wang, Y., 2016. Investigation into gas production from natural gas hydrate: a review. Appl. Energy 172, 286–322. Liu, C., Meng, Q., He, X., Li, C., Ye, Y., Lu, Z., Zhu, Y., Li, Y., Liang, J., 2015a. Comparison of the characteristics for natural gas hydrate recovered from marine and terrestrial areas in China. J. Geochem. Explor 152, 67–74. Liu, C., Meng, Q., He, X., Li, C., Ye, Y., Zhang, G., Liang, J., 2015b. Characterization of natural gas hydrate recovered from Pearl River Mouth Basin in South China sea. Mar. Pet. Geol. 61, 14–21. Marinakis, D., Varotsis, N., Perissoratis, C., 2015. Gas hydrate dissociation affecting the permeability and consolidation behaviour of deep sea host sediment. J. Nat. Gas. Sci. Eng. 23, 55–62. Masui, A., Miyazaki, K., Haneda, H., Ogata, Y., Aoki, K., 2008. Mechanical characteristics of natural and artificial gas hydrate bearing sediments. In: Proceedings of the 6th International Conference on Gas Hydrates (ICGH2008), Vancouver, British Columbia, Canada, July 6–10, 2008. Max, M.D., Johnson, A.H., 2014. Hydrate petroleum system approach to natural gas hydrate exploration. Pet. Geosci. 20 (2), 187–199. Max, M.D., Johnson, A.H., Dillon, W.P., 2005. Economic Geology of Natural Gas Hydrate (Vol. 9). Springer Science & Business Media. McGrail, B.P., Zhu, T., Hunter, R.B., White, M.D., Patil, S.L., Kulkarni, A.S., 2004. A new method for enhanced production of gas hydrates with CO2. Gas hydrates energy Resour. potential Assoc. Geol. hazards 12–16. 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. 116, B06102. Miyazaki, K., Masui, A., Tenma, N., Ogata, Y., Aoki, K., Yamaguchi, T., Sakamoto, Y., 2010. Study on mechanical behavior for methane hydrate sediment based on constant strain-rate test and unloading-reloading test under triaxial compression. Int. J. Offshore Polar Eng. 20, 61–67. Moridis, G.J., Collett, T.S., Pooladi-Darvish, M., Hancock, S., Santamarina, C., Boswell, R., Kneafsey, T., Rutqvist, J., Kowalsky, M., Reagan, M.T., Sloan, E.D., Sum, A.K., Koh, C., 2011a. Challenges, uncertainties, and issues facing gas production from gas-hydrate deposits. SPE Reserv. Eval. Eng. 14 (1), 76–112. Moridis, G.J., Reagan, M.T., Boyle, K.L., Zhang, K., 2011b. Evaluation of the gas production potential of some particularly challenging types of oceanic hydrate deposits. Transp. Porous Media 90 (1), 269–299.

13