Thermo-Hydro-Mechanical Coupled Modeling of Methane Hydrate-Bearing Sediments: Formulation and Application

Abstract

We present a fully coupled thermo-hydro-mechanical formulation for the simulation of sediment deformation, fluid and heat transport and fluid/solid phase transformations occurring in methane hydrate geological systems. We reformulate the governing equations of energy and mass balance of the Code_Bright simulator to incorporate hydrate as a new pore phase. The formulation also integrates the constitutive model Hydrate-CASM to capture the effect of hydrate saturation in the mechanical response of the sediment. The thermo-hydraulic capabilities of the formulation are validated against the results from a series of state-of-the-art simulators involved in the first international gas hydrate code comparison study developed by the NETL-USGS. The coupling with the mechanical formulation is investigated by modeling synthetic dissociation tests and validated by reproducing published experimental data from triaxial tests performed in hydrate-bearing sands dissociated via depressurization. Our results show that the formulation captures the dominant mass and heat transfer phenomena occurring during hydrate dissociation and reproduces the stress release and volumetric deformation associated with this process. They also show that the hydrate production method has a strong influence on sediment deformation.

1. Introduction

Methane hydrates hold vast amounts of methane gas below the sea floor and in permafrost regions [1]. If economically producible, hydrate-sourced methane could supply 10% of the global methane consumption in the coming decades [2,3,4] and bring a new level of energy self-sufficiency to countries that lack conventional reserves (e.g., Japan, India, and South Korea) [5]. Current techniques for methane production from hydrate include thermal stimulation, depressurization, and inhibitor injection [6,7], among which, depressurization is deemed the most mature approach. These techniques perturb the thermodynamic and chemical conditions of the reservoir to destabilize the hydrate phase and force its dissociation into water and gas. The release of mobile phases during dissociation and the loss of the solid hydrate phase has a significant impact in the mechanical (e.g., strength, stiffness, volumetric behavior) and hydraulic (e.g., retention curve, hydraulic conductivity) properties of the sediment and may destabilize it (e.g., [8,9]).

Over the last 20 years, several numerical models have been proposed to simulate the behavior of gas hydrate reservoirs. Initially, modeling efforts focused on evaluating the productivity of methane extraction from in situ sediments. As a result, numerous coupled thermo-hydraulic (TH) models were developed (e.g., [10,11,12,13,14,15,16,17]), which improved considerably the understanding of the TH behavior of MHBS. However, these models generally disregarded the mechanical response of the sediment or considered it with simple approaches. Numerous experimental studies have investigated the impact of hydrate dissociation on the mechanical properties of MHBS. These studies show that the sediment tends to weaken during hydrate dissociation due to a reduction of the sediment stiffness and strength with decreasing hydrate saturation (e.g., [18,19,20,21]). Additionally, they report significant changes in porosity, permeability and stress state of the sediment during dissociation (e.g., [22,23]). Field observations also note the impact of hydrate dissociation in the large-strain response of MHBS (e.g., [15,24,25,26]). Borehole instabilities, gas blowouts and seafloor subsidence have been widely reported during conventional hydrocarbon exploitation nearby hydrate reservoirs (e.g., [8,9,27]). Moreover, large underwater landslides have been related to hydrate dissociation in continental margins (e.g., [25,26,28]).

Hence, geomechanics of MHBS has recently become a key factor in reservoir simulation to ensure well-bore stability in production scenarios and mitigate hydrate-related geohazards [29]. Even though knowledge of the mechanical behavior of MHBS is still limited, current modeling efforts have focused on developing suitable constitutive laws for capturing it. Table 1 lists some of the most notable thermo-hydro-mechanical (THM) formulations that incorporate advanced constitutive models for MHBS. The models listed in Table 1 consider from linear elastic to elastoplastic and visco-plastic formulations in which different dependencies of the yield function with the hydrate saturation are assumed.

Here, we develop a new fully coupled THM formulation to simulate the mechanical behavior of MHBS, particularly focusing on hydrate dissociation scenarios. The formulation builds on the established finite element simulator Code_Bright [37], which has been highly validated at solving multiphase mass and heat transport problems in the geological media (e.g., [38,39,40]) and recently applied to examine the THM behavior of MHBS [36]. We reformulate Code_Bright’s governing equations of energy and mass balance in terms of the potential porosity of the sediment (i.e., space between the mineral grains) and the remaining available porosity after the formation of hydrate and/or ice. This volumetric distinction allows isolation of the effects of mechanical deformation from the effects of hydrate and ice fluid/solid phase transformations, both affecting the hydraulic and mechanical properties of the porous medium. We also implement the elastoplastic constitutive model Hydrate-CASM [41] to capture the effect of hydrate saturation in the sediment stress-strain response, and integrate the Peng-Robinson equation of state (EoS) [42] and the thermodynamic equations proposed by [43] to compute methane gas density and solubility, respectively.

The TH capabilities of our formulation are validated against the results from a series of state-of-the-art simulators involved in the NETL-USGS first international gas hydrate code comparison study [44]. Finally, the mechanical coupling is investigated by modeling synthetic dissociation tests and validated by reproducing published experimental data from triaxial tests performed in hydrate-bearing sands dissociated via depressurization [21].

2. Methodology This section sets the theoretical framework of the problem, develops the mathematical formulation proposed to simulate the THM behavior of MHBS and describes the strategy adopted for the numerical solution. 2.1. Components, Phases and Partial Saturations

The multicomponent and multiphase approach followed by [37] is adopted here to describe the MHBS system (Figure 1). The porous medium is considered to be composed by four mass components; mineral grains (gr), methane (m), water (w), and salt (st), that can be partitioned among five possible phases; solid (s), hydrate (h), ice (i), gas (g) and liquid (l). The incompressible mineral grains form the non-reactive solid continuum that provides the skeletal structure to the porous medium. Within the pores, hydrate and ice can grow as solid and immobile phases. Finally, fluid flow and storage are restricted to the available space between the mineral grains and hydrate and ice phases.

Similar to [45], the pore-space of the sediment is divided into two porosities (Table 2) to isolate the effects of mechanical deformation from those related to hydrate and ice phase change on both the hydraulic and mechanical properties of the porous medium.

The potential porosity (_ϕp) is the volume of pores between the mineral grains of the sediment. This porosity is strictly related to the volumetric deformation of the sediment and is the maximum available space for the formation of hydrate and ice phases. The available porosity (_ϕa ) is the volume of pores existing between solid phases (i.e., part of the potential porosity non-occupied by hydrate or ice). This porosity determines the available space for the fluids to flow or be stored, governs the sediment permeability and evolves according to variations in potential porosity and hydrate and ice phase transformations. Each phase saturation is obtained as the ratio between the its volume and the corresponding volume of voids in which it is partitioned (Table 2). Thus, the following volumetric restrictions apply:

_Sl+_Sg=1

and

_{ϕa ϕp}=(1−_Si)(1−_Sh)

Equation (1) allows the preservation of the basic structure of the classical two-phase flow models for the gas-water system and Equation (2) provides the coupling between the gas-water system and the solid pore-filling phases.

2.2. Multiphysical Coupled System

The behavior of the multiple phases considered within the pores involves four physical processes: mechanical deformation, fluid flow, thermal flow and phase change reactions (Figure 2). Each of these processes interacts with each other generating a multiphysical coupled system. Accurate representation of the interactions between these processes in the mathematical formulation is vital for adding certainty to the simulation of the mechanical behavior of MHBS.

2.3. Governing Equations

The governing equations are divided here into two main groups; (i) balance equations and (ii) constitutive equations and equilibrium constraints. All the parameters used in the formulation are listed and defined in Table 3.

(i) Balance equations

Mass balance equations

The mass balance equations are established following a compositional approach. In these equations the volumetric mass of a component in a phase (e.g., methane in the liquid phase,_θlm) is the product of the mass fraction of that component (_ωlm) and the bulk density of the phase (_ρl).

• Mass balance of mineral grains:

  The mineral grains coincide with the permanent solid phase and define the skeletal structure of the porous medium. The mass conservation of this component can be written as:

  ∂∂t_ρs(1−_ϕp)+∇(_jsgr)=0

  where_jsgris the flux of solid and can be expressed in terms of the solid velocity (∂u∂t) as:

  _jsgr=_ρs(1−_ϕp)∂u∂t

  Applying the chain rule for all the derivatives, Equation (3) can be rewritten as:

  (1−_ϕp)∂_ρs∂t+(1−_ϕp)(∂u∂t∇_ρs)−_ρs∂_ϕp∂t−_ρs(∂u∂t∇_ϕp)+_ρs(1−_ϕp)∇∂u∂t=0

  Neglecting the gradients of density and porosity convected by the solid phase and under the assumption of small strain, Equation (5) can be rewritten as:

  ∂_ϕp∂t=1_ρs(1−_ϕp)∂_ρs∂t+(1−_ϕp)∇∂u∂t

  where the term∇∂u∂tcan be expressed in the form:

  ∇∂u∂t=^∂2 _ux∂x∂t+^∂2 _uy∂y∂t+^∂2 _uz∂z∂t=∂∂t(_εx+_εy+_εz)=∂_εv∂t

  and_εx+_εy+_εzare Cartesian strains.

• Mass balance of methane:

  Methane component is present in liquid, gas and hydrate phases, and its total mass balance is expressed as:

  ∂∂t_{θlm Sl}+_{θgm Sgϕa}+_{θhm Sh ϕp}+∇(_jlm+_jgm+_jhm)=^fm

  where_θlm,_θgmand_θhmare the volumetric mass of methane in the liquid, gas and hydrate phases, respectively. For a given temperature, pressure and salinity the term_θlm can be obtained according to [42,43],_θgmis equal to one as the gas phase is considered mono-component and_θhm=11+_nhassuming_nhas 6.176.

  The mass flux terms_jlmand_jgmare the relative motion of methane in the liquid and gas phases, respectively, with respect to the solid phase. These terms are obtained as the sum of advective and diffusive flux terms as follows:

  _jlm=_ilm+_{θlm ql}+_{ϕa Sl θlm}∂u∂t

  _jgm=_igm+_{θgm qg}+_{ϕa Sg θgm}∂u∂t

  The mass flux term_jhmdenote the relative motion of methane in the hydrate phase with respect to the solid phase as a result of the medium deformation:

  _jhm=_{ϕp Sh θhm}∂u∂t

  The term^fmis the external sink/source of methane per unit volume. Please note that because of the compositional approach adopted in the formulation, this term do not include methane mass changes from hydrate kinetics.

• Mass balance of water:

  Water is present in liquid, ice and hydrate phases but it is neglected as vapour in the gas phase. Thus, the water total mass balance of water is expressed as:

  ∂∂t_{θlw Sl ϕa}+_{θiw Si}(1−_Sh)+_{θhw Shϕp}+∇(_jlw+_jiw+_jhw)=^fw

  where_θlw,_θiwand_θhware the volumetric mass of water in the liquid, ice and hydrate phases, respectively. The term_θlwdepends on the concentration of salt and methane dissolved in the liquid phase,_θiwis equal to one and_θhw=1−_θhm.

  The mass flux term of water in liquid (_jlw ) is computed similarly as in Equation (9), while the flux terms in hydrate (_jhw) and ice phases (_jiw ) are computed similarly as in Equation (11) for_Shand_Si, respectively.

  The^fwis an external sink/source of water per unit volume and do not include water mass changes from hydrate kinetics.

• Mass balance of salt:

  Salt is present as a dissolved component in the liquid phase, and it is not allowed to precipitate as a solid phase. Its concentration modifies the liquid density, influences the solubility of methane and can inhibit hydrate stability conditions. The total mass balance of salt is expressed as:

  ∂∂t_{θlst Sl ϕa}+∇_jlst=^fst

  where_θlstis the volumetric mass of salt in the liquid phase. Changes in salt concentration due to “freshening” of the pore water can be used as a tracer of ongoing hydrate dissociation and/or ice melting.

  The flux term of salt in the liquid phase (_jlst ) is the relative motion of salt in the liquid phase with respect to the solid phase and is computed similarly as in Equation (9). Finally, the term^fstis the external sink/source term of salt.

Energy balance equation:

The equation for internal energy balance in the porous medium is established according to [37] considering the internal energy in each phase of the system, so that:

∂∂t(_{El ρl Sl}+_{Eg ρg Sg})_ϕa+_{Ei ρi Si}(1−_Sh)+_{Eh ρh Shϕp}+_{Es ρs}(1−_ϕp)+∇(_ic+_jEl +_jEg +_jEi +_jEh +_jEs )=^fQ

where the energy flux of fluid phases include energy dispersion, advective transport by fluid mass relative to the solid phase and that related to the soil velocity:

_jEl =_iEl +_{jlm Elm}+_{jlw Elw}+_{jlst Elst}+_{El ρl ϕa Sl}∂u∂t

_jEg =_iEg +_{jgm Eg}+_{Eg ρg ϕa Sg}∂u∂t

with

_{El ρl}=_{θlm Elm}+_{θlw Elw}+_{θlst Elst}

Please note that_Elmand the dispersive terms (_iEl and_iEg ) have been neglected in our simulations.

The energy flux of solid phases only depend on the advective transport related to the soil velocity, so that:

_jEh =_{Eh ρh ϕp Sh}∂u∂t

_jEi =_{Ei ρi ϕp}(1−_Sh)_Si∂u∂t

_jEs =_{Es ρs}(1−_ϕp)∂u∂t

For which exothermic/endothermic changes caused by hydrate formation/dissociation and ice formation/melting are considered through the latent heat of each phase transformation (Table 4) and the term^fQis the external supply of energy per unit volume.

Momentum balance equation:

The momentum balance equation for the bulk porous medium reduces to the equilibrium of total stresses if the inertial terms are neglected:

∇σ+b=0

Equation (21) is formulated based on the infinitesimal strain theory, uses an additive decomposition of the total strain tensor and assumes a Biot’s coefficient_αB = 1 for the fluid-solid coupling (Table 5).

The constitutive equations for MHBS allow rewriting Equation (21) in terms of the solid velocities, gas and fluid pressure and temperature.

(ii) Constitutive equations and equilibrium constraints

A set of constitutive laws and equilibrium restrictions are used in the formulation to establish the link between the primary (i.e., solid displacements,u; liquid pressure,_Pl; gas pressure,_Pg ; temperature, T and salinity, S) and the dependent variables of the system (listed in Table 5 and Table 6). These equations allow capturing the coupling among the various physical phenomena considered in the system and close the mathematical formulation.

Stress - strain behavior of MHBS

Several experimental studies show that MHBS have a greater stiffness, strength, dilation, and softening behaviors with increasing hydrate saturation (e.g., [55,56,57,58,59,60,61]). The increase in sediment strength due to hydrate has been widely attributed to a physical bonding between the hydrate crystals and the sediment grains. This bonding has been modeled using different strategies; (i) addition of a cohesion constituent in the failure criteria [60,62,63,64], (ii) enlargment of the yield surface by cohesion and dilation [65,66], (iii) partition of the stress between hydrate and matrix in a bonding damage framework [65,67], (iv) attribution of physical bonding properties in discrete element methods (DEM) [68] and (v) expansion of the failure envelope in a spatially mobilized plane (SMP) model [69]. Refs. [70,71] propose that kinematics might govern the increase in strength observed in hydrate-bearing sands. Alternatively, [41] propose that the greater strength and dilatancy observed in MHBS can be explained by densification and stiffening of the host sediment due to pore invasion by hydrate.

The constitutive model Hydrate-CASM [41] is integrated in our formulation to predict the mechanical response of MHBS. This elastoplastic model extends the formulation of the CASM model [72] by implementing the subloading surface model [73] and the densification mechanism. The densification is a novel mechanism for MHBS that attributes hydrate-related changes in the void ratio, the swelling line slope and the volumetric yield stress of the sediment to stress-strain changes (Table 5). The performance of the Hydrate-CASM model has been tested over a range of hydrate saturations and confining stress using experimental data from triaxial tests performed in synthetic hydrate-bearing sands and validated against the outputs from other notable constitutive models for MHBS that simulated the same experimental data (Figure 3).

Hydrate and ice phase boundaries

The dominant parameters governing the stability field of methane hydrates in marine sediments are pressure, temperature and salinity. Our formulation employs the empirical solution proposed by [43] to define the P-T conditions that describe the equilibrium of methane hydrates in saline water (Figure 4). Above the hydrate phase equilibrium curve (i.e.,_Pg>_Peq , shaded area in Figure 4), hydrate can form when the concentration of methane dissolved in the liquid phase is at or above saturation value [75]. Our formulation computes methane solubility at dissociation pressure as a function of salinity and temperature following the empirical proposal of [43] (Table 6). This function is combined with the IUPAC recommendation at_Pg = 0.1 MPa [51] to extend the calculation of methane solubility in water at any other given pressure (Table 6). The resulting estimation of methane concentration in water (_ωlm ) is key to model the upward transport of gas by means of diffusion and advection. It can also be used to predict the resulting shape of sulfate profiles in pore water (e.g., [76,77]) and the volume and distribution of hydrate formations in marine sediments (e.g., [78]).

At P-T conditions below the hydrate phase equilibrium curve (i.e.,_Pg<_Peq , white area in Figure 4), hydrate dissociates into its forming components. Our formulation assumes hydrate phase change to be governed by first order kinetics as suggested by [52,53], in which the rate of phase change is proportional to the product of the hydrate surface area and the driving force (i.e., pressure difference between the phase equilibrium and gas pore pressure) (Table 6).

Water-to-ice transformation may also take place in MHBS during temperature dropping if hydrate formation is restricted by insufficient methane gas, or due to endothermic cooling at fast depressurization. Ice changes are computed using the freezing characteristic function proposed by [54] (Table 6). This function is derived from the thermodynamic equilibrium between liquid and ice phases in pores in frozen soils and employs van Genuchten’s (1980) model to relate the degree of unfrozen water saturation to the thermodynamic properties of the sediment. As stated in Table 5 our formulation consider changes in ice saturation to have the same mechanical effects that hydrate phase change.

2.4. Numerical Solution Strategy

The PDE’s (Partial Differential Equations) system presented is solved numerically following the procedure established in [37]. The Galerkin finite element method (FEM) is adopted for the spatial discretization, while finite differences are used for the temporal discretization via an implicit scheme (backward Euler method) that incorporates an automatic sub-stepping procedure based on error control [79]. Also, a mass conservative approach is adopted by directly discretizing the storage terms [80,81], and the mechanical problem is solved using an implicit stress point algorithm (SPA) [82].

Once the formulation is discretized, this can be represented as a system of algebraic equations that are linearized via the Newton-Raphson scheme. Then, LU decomposition and back-substitution (non-symmetric matrix) or conjugate squared gradients are used to solve the system of equations simultaneously (monolithic solution) in a fully coupled manner. 3. Results and Discussion 3.1. Thermo-Hydraulic Validation

In this section, we reproduce the benchmark problem P2 proposed in the NETL-USGS first international gas hydrate code comparison study. The initial problem conditions, model parameters and further specifications can be found in https://www.netl.doe.gov/node/7285. Problem P2 analyzes a 20 m one-dimensional horizontal closed domain (no flow boundary conditions) initialized with gradients in liquid and gas pressure, and temperature that yield aqueous saturated conditions on half of the domain and aqueous unsaturated conditions on the other half (Figure 5a). We have simulated the problem in two-dimensions and have used a non-uniform mesh, with a higher concentration of elements in the contact between the two domains, to represent better the problem solution.

The initial conditions of the problem P2 are designed to trigger complete dissociation of the hydrate phase via the thermal capacitance of the right half of the domain. From these conditions, the simulation proceeds to reach the equilibrium in temperature and pressure along the domain. During the transition, the hydrate begins to dissociate at the contact between both domains and the dissociation front (vertical dotted line Figure 5e) advances progressively to the left. After 10 days of simulation hydrate dissociation occurs simultaneously with hydrate reformation. As shown in Figure 5e, hydrate saturation at 100 and 1.000 days of simulation increases with respect to the previous value as we move away from the dissociation front towards the left side. The low temperature in this side of the domain in combination with the endothermic nature of the dissociation reaction and the increase in pore pressure due to the migration of free gas allows reaching the P-T conditions for hydrate reformation.

The good matching observed between our results and those from the simulators involved in the code comparison study (i.e., HydResSim, MH-21, CMG STARS, STOMP-HYD, TOUGH-FX and Univ-Houston) allows validating the capacity of our formulation at capturing the dominant mass and heat transfer phenomena during gas hydrate dissociation. 3.2. THM Modelling of Synthetic Experimental Tests

In this section, various synthetic triaxial tests are simulated to investigate the stress-strain and volumetric deformation characteristics of MHBS during hydrate dissociation under strain controlled shear. For the simulations, we assume two-dimensional axisymmetric jacketed sand specimens of 10 cm × 5 cm, of which a radial plane of 10 cm × 2.5 cm is used as computational domain (shaded area in Figure 6). The right-hand side of the domain has a boundary condition of constant total stress and the bottom is immobile. All the simulations are conducted under constant compression rate of 0.012%/min with drained conditions only at the top of the sample. The particular conditions employed in each simulation are listed in Table 7.

We label as Hs_Temp and Hs_Dep the two simulations used to characterize the mechanical response of the host sediment specimen (_Sh=0%) subjected to thermal stimulation (8 → 11.3 °C) and depressurization (6 → 4.6 MPa). These simulations provide a reference behavior that allows distinguishing the mechanical effects of temperature and pressure variations from those related to hydrate dissociation in MHBS. For the dissociation tests, MHBS specimens (_Sh≠0%) are sheared up to an axial strain of approximately 16% (Mhbs_Temp1, Mhbs_Dep) and 40% (Mhbs_Temp2). Then, depressurization and/or thermal stimulation are imposed over a period of 100 seconds and the final liquid pressure (4.6 MPa) and temperature (11.3 °C) are kept constant thereafter.

Figure 7 shows the simulated evolution of the deviatoric stress, potential porosity, temperature, liquid pressure and hydrate and gas saturation during the complete triaxial loading history for each simulation.

The host specimen Hs_Temp, sheared under constant pore pressure and subjected to thermal stimulation, is characterized by a peak in the deviatoric stress followed by softening (Figure 7a) and dilatant response (Figure 7c). The specimens with hydrate (Mhbs_Temp1 and Mhbs_Temp2) show a greater peak in the deviatoric stress and a more dilatant response than the corresponding host sediment before dissociation is simulated. Then, dissociation via thermal stimulation leads to a drop in the deviatoric stress (points 1 and 2 in Figure 7a) as well as to a reduction in the dilatant response (point 1 in Figure 7c) or a compressive behavior of the dissociated sediment (point 2 in Figure 7c).

The mechanical behavior of the host specimen Hs_Dep subjected to depressurization is highly influenced by the increase in effective confining stress associated with it. Depressurization leads the host sediment to shift from a softening to a hardening response and from a dilatant to a compressive volumetric behavior (point 1 in Figure 7b,d). For the hydrate-bearing specimen Mhbs_Dep, the effects on the effective stress induced by depressurization adds to the mechanical effect of hydrate dissociation. Upon dissociation specimen Mhbs_Dep starts contracting under shear (point 2 in Figure 7d) rather than reducing its dilatant behavior, as observed in the equivalent specimen subjected to thermal stimulation (point 1 in Figure 7c).

These results capture two important features of the mechanical response of MHBS during hydrate dissociation under strain controlled shear. Hydrate dissociation leads to a drop of the deviatoric stress which is associated with the shrinking of the yield surface due to sediment mechanical weakening [41]. In consequence, the sediment shows a less dilatant behavior than the corresponding host sediment or even a contractant response under shear. Our results also show that the volumetric deformation of the sediment after dissociation strongly depends on the hydrate production method. Point 2 in Figure 7d clearly shows that dissociation via depressurization leads the specimen contraction while dissociation by thermal stimulation does not (point 1 in Figure 7c).

3.3. THM Modelling of the Experimental Tests of Li et al. (2018)

The multistage triaxial shear test conducted by [21] is used in this section to examine the capabilities of our formulation at capturing the effect of hydrate dissociation via depressurization in the mechanical response of synthetic MHBS. Figure 8 summarizes the mechanical properties and the initial and boundary conditions considered in the simulation. The modeling is approached here in two steps, firstly the mechanical behavior of the specimen with a constant hydrate saturation of 38.5% is calibrated based on the experimental deviatoric stress-axial strain and volumetric strain-axial strain relationships obtained before depressurization is applied (black solid lines in Figure 9a,b). Once the mechanical properties are calibrated, hydrate dissociation is simulated at approximately 2.5% of axial strain by depressurizing the specimen from a liquid pressure of 8 MPa to 3 MPa. The depressurization is applied as a boundary condition on top of the specimen (blue line in Figure 8) while this is subjected to a constant triaxial compression rate of 0.2%/min. Temperature is imposed on the top, base and right-hand side of the computational domain to simulate the circulation of fluid at a constant temperature around the specimen.

Results in Figure 9 show that once depressurization is induced, our formulation captures a subtle drop in the deviatoric stress (point 1 in Figure 9a) as well as an increase in volumetric strain (points 2 to 3 in Figure 9b). However, the model overestimates the maximum deviatoric stress after dissociation and slightly underestimates the volumetric strain at the end of the triaxial. For comparison purposes, the evolution of hydrate saturation with the axial strain is plotted removing the change in_Sh due to sediment deformation (i.e., changes in hydrate saturation are related only to the chemical reaction). Figure 9d shows that the model captures the rate of hydrate dissociation.

4. Conclusions Hydrate dissociation has a significant effect on the THM properties of MHBS and may cause its mechanical destabilization. Our formulation captures the mass and heat transfer phenomena during hydrate dissociation in the porous medium. Particularly, the TH performance shows that the migration of the dissociation front caused by thermal stimulation and depressurization and the volume of hydrate that reforms due to gas migration are closely reproduced. Regarding the coupling with the mechanical formulation, our synthetic results show that dissociation leads to a drop in the deviatoric stress and to a less dilatant or even contracting behavior of the dissociated MHBS in comparison to the host sediment. Particularly, it is observed that sediment may contract in cases when the deviatoric stress drops below the failure envelope of the host sediment. Our synthetic results also suggest that the deformation of MHBS strongly depends on the hydrate production method. They show that specimens subjected to hydrate dissociation via depressurization are more likely to show contractant behavior under shear than those subjected to thermal stimulation because of the increase in effective stress. We finally illustrate that our formulation successfully reproduces the main THM features observed in dissociation tests performed in real specimens.

Model Reference	Mechanical Approach
Kimoto et al. (2010) [30]	Viscoplasticity withSh dependency
Rutqvist (2011) [31]	Mohr-Coloumb elastoplasticity withSh dependency
Kim et al. (2012) [32]	Mohr-Coloumb elastoplasticity withSh dependency
Klar et al. (2013) [33]	Mohr-Coloumb elastoplasticity withSh dependency
Gupta et al. (2016) [34]	Poroelasticity withSh dependency
Sun et al. (2018) [35]	Thermodynamics-based elastoplastic model withSh dependency
Sánchez et al. (2018) [36]	Elastoplasticty withShdependency + Damage model

Bulk volume	Vt=Vs+Vh+Vi+Vl+Vg
Potential void space	Vp=Vt−Vs
Potential porosity	ϕp=Vp/Vt
Available void space	Va=Vp−Vh−Vi
Available porosity	ϕa=Va/Vt
Hydrate saturation	Sh=Vh/Vp
Ice saturation	Si=Vi/(Vp−Vh)
Gas saturation	Sg=Vg/Va
Liquid saturation	Sl=Vl/Va

Roman Symbols
a(T)	Peng-Robinson EoS parameter	M	Slope of critical state line in p’-q space
Aβ	Viscosity parameter	Mβ	Molecular mass of phaseβ
A(T)	Viscosity thermal function	n	Stress-state coefficient: yield surface shape parameter
Ah	Hydrate surface area	nh	Hydration number
b	Body forces	p′	Mean effective stress
b	Peng-Robinson EoS constant	P0	Van Genuchten parameter
Bβ	Viscosity parameter of phaseβ	p0	Isotropic yield stress of hydrate-free sediment
cβα	Specific heat of componentαin phaseβ	p0h	Isotropic yield stress of MHBS
dmβ	Mass change of phaseβ	Peq	Hydrate phase equilibrium pressure
Dβα	Diffusion coefficient of componentαin phaseβ	Pβ	Pressure of phaseβ
e	Void ratio of hydrate-free sediment	Pβ0	Reference pressure for phaseβ
eah	Available void ratio of the MHBS	Pc	Critical pressure
eh	Hydrate void ratio	Pp	Pore pressure
Eβ	Energy of phaseβ	q	Deviatoric stress
Eβα	Energy of componentαin phaseβ	qβ	Advective fuid flow
f	Hydrate-CASM yield function	r	Yield surface spacing ratio
fα	External mass supply of componentα	R	Subloading ratio
fQ	Internal/external energy supply	Rg	Regnault constant
fw	Peng-Robinson EoS function	Rh	Rate of hydrate mass change
g	Gravity forces	S	Salinity
iβα	Diffusive flux of componentαin phaseβ	Sβ	Saturation of phaseβ
ic	Conductive heat flow	Se	Effective liquid saturation
iEβ	Energy dispersivity in phaseβ	Sls	Maximum liquid saturation
I	Identity matrix	Srl	Residual liquid saturation
jβα	Mass flux of componentαin phaseβ	t	Time
jEβ	Advective flux of energy of phaseβ	T0	Reference temperature
k0	Intrinsic permeability of hydrate-free sediment	T	Temperature
k	Intrinsic permeability of MHBS	Tc	Critical temperature
krβ	Relative permeability of phaseβ	Teq	Hydrate phase equilibrium temperature
KCH4	Solubility of methane in water	Tr	Reduced temperature
Kd	Hydrate dissociation constant	u	Displacement
Kf	Hydrate formation constant	u	Subloading parameter controlling the
Lh	Latent heat of hydrate dissociation		plastic deformations before yielding
Li	Latent heat of ice melting	v	Molar volume
m	Van Genuchten parameter	V	Unitary volume
Greek Symbols
αB	Biot’s coefficient	∇	Differential operator = (∂∂x;∂∂y;∂∂z)
βt	Thermal expansion coefficient of the liquid phase	ν	Poisson’s ratio
∂	Partial derivative	ωβα	Mass fraction of componentαin phaseβ
|εp|	Norm of the incremental plastic strain	μβ	Viscosity of phaseβ
εv	Volumetric strain	ϕp	Potential porosity
ε(x,y,z)	Cartesian strains	ϕa	Available porosity
γ	Solute variation	ρβ	Density of phaseβ
κ	Slope of swelling line of hydrate-free sediment	ρβ0	Reference density of phaseβ
κh	Slope of swelling line of MHBS	σ	Cauchy total stress (compression positive)
κrf	Swelling line slope reduction factor	σ′	Effective stress (compression positive)
λ	Slope of normal compression line	σc	Confining stress (compression positive)
λc	Composite thermal conductivity	θβα	Volumetric mass of componentαin phaseβ
λdry	Dry thermal conductivity	τ	Tortuosity coefficient
λsat	Liquid saturated thermal conductivity	ϑ	Peng-Robinson acentric factor

Phase	Specific Energy	Specific and Latent Heat
Gas (g)	Eg=cg(T−T0)	cg = 2500 J (Kg K)−1
Hydrate (h)	Eh=Lh+ch(T−T0)	ch = 2108 J (Kg K)−1; †Lh = 3.39 e5 J Kg−1
Ice (i)	Ei=Li+ci(T−T0)	ci = 3144 J (Kg K)−1; †Li = 3.34 e5 J Kg−1
Liquid (l)	Elw=clw(T−T0)	clw = 4184 J (Kg K)−1
	Elst=1.42e5+clst(T−T0)	clst = 2200 J (Kg K)−1
Solid (s)	Es=cs(T−T0)	cs = 874 J (Kg K)−1

†_Lhand_Li are positive for hydrate dissociation and ice melting respectively.

Constitutive Law	Equation	Dependent Variable
Advective fluid flow
Darcy’s law	qβ=−kkrβ μβ(∇Pβ−ρβg)β=l,g	Advective fluid flow,qlandqg
Kozeny’s model	k=k0ϕa3 (1−ϕa)2(1−ϕp)2 ϕp3
	k0=1×10−13 m2
	krl =Se1−(1−Se1/m)m2
	krg =1−krl
	m‡=0.645
	μβ=Aβexp(Bβ273.15+T)
	Al=2.1×10−12MPa s;Bl= 1808.5 K
	Ag=1.48×10−12MPa s;Bg = 119.4 K
Retention curve [46]	Se=Sl−SrlSls−Srl=1+Pg−PlP011−m−m	Saturation of mobile phases,SlandSg
	P0‡=0.075 MPa
Liquid density [37]	ρl=ρl0 1+βt(Pl−Pl0 )+A(T)+γωlst	Liquid density,ρl
	A(T)=(T+288.9414)(T−3.9863)250,8929.2(T+68.12963)
	ρl0 = 1002.6 kg/m3
	βt=4.5×10−4 MPa−1
	Pl0 =0.1MPa
	γ=0.6923
	ωlst=0
Gas density [42]	Pg=RgTv−b−a(T)v(v+b)+b(v−b)	Gas density,ρg
	a(T)=0.45724Rg2 Tc2/Pc [1+fw(1−Tr0.5)]2
	fw=0.37464+1.54226ϑ−0.26992ϑ2
	ϑ=0.0015
	b=0.0778RTc/Pc
	Rg=8.3144598J/mol K
	Pc=4.60MPa
	Tr=T/Tc,Tc=190.4K
	ρg=Mg PgRgT
	Mg=0.016042 Kg/mol
Non-advective fluid flow
Fick’s law	iβα=ϕaτρl Sl DlαI∇ωlαα=m,st	Diffusive flux,ilmandilst
	Dlm=5.9×10−6exp(273.15+T)2.3 Pg
	Dlst=1.1×10−4exp−24530Rg(273.15+T)
	τ=1
Fourier’s law	ic=−λc∇T	Conductive heat flow,ic
	†λc=λsatSl+λdry(1−Sl)
	λsat=2W/mK,λdry=2.18 W/mK
Stress-strain behavior
Effective stress [47]	σ′=σ−αB PpI	Effective stress tensor,σ′
	αB=1
Hydrate-CASM [41]	Subloading yield function	Stress tensor,σ
	f=(qMp′)n+1ln(r)ln(p′Rp0h′)
	dR=−ulnR|dεp|
	Densification mechanism
	eah=e−eh;eh=e(Sh+Si)
	κh=κκrf
	κrf= 0 if(Sh+Si)=0
	κrf=3(Sh+(1−Sh)Si)2−2.68(Sh+(Sh+(1−Sh)Si)+0.9934
	if0<(Sh+(1−Sh)Si)≤0.42
	p0h′=ee(Sh+Si)λ−κh p′0λ−κλ−κh

†_λcis not computed in terms of_Shor_Si . This simplification valid for low saturations of ice [48]; ‡ m and ‡_P0 values correspond to those given in the benchmark problem analyzed in Section 3.1. [44]. The m value is included within the range of values reported in the literature for gas hydrate numerical simulation studies (e.g., [49,50]); Please note that_eah,_κhand_p0h′recover the hydrate-free parameterse,κand_p0′when(_Sh+_Si) = 0; The use of bold symbols denote vectors and tensors.

Equilibrium Restriction	Equation	Dependent Variable
Hydrate phase change
Hydrate phase boundary [43]	ln(Peq)=−1.644866×103−0.1374178T+5.4979866×104/T+2.64118188×102ln(T)+S(1.1178266×104+7.67420344T−4.515213×10−3 T2−2.04872879×105/T−2.17246046×103ln(T))+S2(1.70484431×102+0.118594073T−7.0581304×10−5 T2−3.0979619×103/T−33.2031996ln(T))	Equilibrium pressure,Peq
Methane solubility	ln(KCH4)Pg =ln(KCH4)Peq ln(KCH4)0.1MPaln(KCH4)Teq(Pg)	Dissolved methane concentration,ωlm
[43]	ln(KCH4)Peq =−2.5640213×105−1.644805×102T+9.1089042×10−2 T2+4.90352929×106/T+4.93009113×104ln(T)+S[−5.1628513×102−0.33622376T+1.8819904×10−4 T2+9.76525718×103/T+99.523354ln(T)]
[51]	ln(KCH4)0.1MPa=1×10−9−417.5053+599.8626(100/T)+380.3636ln(T/100)−62.0764T/100+S[−0.06423+0.03498(T/100)−0.0052732(T/100)2]
	ln(KCH4)Teq(Pg)=1×10−9−417.5053+599.8626(100/Teq(Pg))+380.3636ln(Teq(Pg)/100)−62.0764Teq(Pg)/100+S[−0.06423+0.03498(Teq(Pg)/100)−0.0052732(Teq(Pg)/100)2
	ωlm=1.604×10−2 KCH4Pg Kg/Kg
Hydrate kinetic rate [52,53]	Rh(T,Pg)=ϕp Sh AhKd〈Peq(T)−Pg〉−Kf〈Pg−Peq(T)〉	Hydrate mass change,dmh
	Ah=0.375μ/m
	Kd=124×103exp(−9400/T(K))mol/m2Pas
	Kf=0.5875×10−11mol/m2Pas
	†dmh=Mh Rh(T,Pg)
	Mh=0.018016nh+0.016042Kg/mol
Ice phase change
Freezing characteristic function [54]	1−† Si=1+1−(1−ρi/ρl)Pl−ρi Liln(T/273.15)P011−m)−m	Ice saturation,Si
	ρi=0.91·ρl

† See Appendix A for the derivation of hydrate and ice mass conservation equations.

Test Conditions
Test Name	Temperature (°C)	Confining Pressure (MPa)	Liquid Pressure (MPa)	Hydrate Saturation (%)	Remarks
Hs_Temp	8 → 11.3	8	6	0	Thermal stimulation atεa=16%
Mhbs_Temp1	8 → 11.3	8	6	20	Thermal stimulation atεa=16%
Mhbs_Temp2	8 → 11.3	8	6	20	Thermal stimulation atεa=40%
Hs_Dep	8	8	6 → 4.6	0	Depressurization atεa=16%
Mhbs_Dep	8	8	6 → 4.6	20	Depressurization atεa=16%

Author Contributions

Formal analysis, M.D.L.F., J.V. and H.M.-M.; Methodology, M.D.L.F., J.V. and H.M.-M.; Simulations, M.D.L.F.; Supervision, J.V. and H.M.-M.; Writing original draft, M.D.L.F.; Writing, review and editing, M.D.L.F., J.V. and H.M.-M.

Funding

Marin-Moreno was supported by the United Kingdom Natural Environment Research Council Standard Grant (NE/J020753/1).

Acknowledgments

We thank two anonymous reviewers for their thoughtful comments that improved the quality of the manuscript. The data used in this paper are published in the cited sources and access to our simulations results can be obtained by contacting the lead author.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

MHBS	Methane hydrate-bearing sediments
NETL	National Energy Technology Laboratory
P	Pressure
T	Temperature
TH	Thermo-hydraulic
THM	Thermo-hydro-mechanical
USGS	United States Geological Survey

Appendix A. Development of Hydrate and Ice Mass Balance Equations

Mass balance of hydrate

The mass conservation of hydrate considers mass exchange due to phase transformations, sediment volumetric deformation and porosity variations due to changes in the solid phase density, so that:

∂∂t_{ρh ϕp Sh}+∇_{ρh ϕp Sh}∂u∂t=d_mh

whered_mhis the mass change of hydrate caused by kinetics.

Neglecting the gradients of density and porosity convected by the solid phase and under the assumption of small strain, Equation (A1) can be rewritten as:

∂(_Sh)∂t=d_{mhρh ϕp}+_Shε˙-_{Sh ϕp}∂(_ϕp)∂t-_{Sh ρh}∂(_ρh)∂t

whereε˙is the sediment volumetric strain rate (sign convention as in continuum mechanics, negative compression):

ε˙=-1(1-_ϕp)∂(_ϕp)∂t+(1-_ϕp)_ρs∂(_ρs)∂t

Replacing Equation (A3), Equation (A2) can be rearranged as:

∂(_Sh)∂t=d_{mhρh ϕp}-_{Sh ϕp}ε˙-_{Sh ρh}∂(_ρh)∂t

The first term in Equation (A4) expresses the effect of hydrate formation/dissociation reaction in the computation ofd_Sh, the second the effect of the volumetric deformation of the sediment and the third one those associated with changes hydrate density variations with temperature. Since our formulation defines_Sh independently of the volume of ice and assumes that hydrate and ice phase transformations do not occur simultaneously (see volumetric definitions in Table 2), Equation (A4) disregards derivatives with respect to ice saturation.

Mass balance of ice

As with Equation (A4) variations in ice saturation can be written as:

∂(_Si)∂t=d_{miρi ϕp}(1-_Sh)-_{Si ϕp}ε˙-_{Si ρi}D(_ρi)Dt+_Si(1-_Sh)∂(_Sh)∂t

where the termd_{miρi ϕp}(1-_Sh) corresponds to changes in ice saturation due to freezing and melting (computed using the freezing characteristic function proposed by [54]), and the term∂(_Sh)∂tdoes only consider changes in hydrate saturation due to hydrate density variations.