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

: We present a fully coupled thermo-hydro-mechanical formulation for the simulation of sediment deformation, ﬂuid and heat transport and ﬂuid/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 ﬁrst 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 inﬂuence on sediment deformation.


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.

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

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.

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. Please note that the formulation considers the existence of unfrozen water below freezing temperature due to capillary effects. 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: and 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.

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.

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.

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, θ m l ) is the product of the mass fraction of that component (ω m l ) 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: where j gr s is the flux of solid and can be expressed in terms of the solid velocity ∂u ∂t as: Applying the chain rule for all the derivatives, Equation (3) can be rewritten as: 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: where the term ∇ ∂u ∂t can be expressed in the form: and ε x + ε y + ε z are Cartesian strains.

•
Mass balance of methane: Methane component is present in liquid, gas and hydrate phases, and its total mass balance is expressed as: where θ m l , θ m g and θ m h are the volumetric mass of methane in the liquid, gas and hydrate phases, respectively. For a given temperature, pressure and salinity the term θ m l can be obtained according to [42,43], θ m g is equal to one as the gas phase is considered mono-component and θ m h = 1 1+n h assuming n h as 6.176.
The mass flux terms j m l and j m g are 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: The mass flux term j m h denote the relative motion of methane in the hydrate phase with respect to the solid phase as a result of the medium deformation: The term f m is 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: where θ w l , θ w i and θ w h are the volumetric mass of water in the liquid, ice and hydrate phases, respectively. The term θ w l depends on the concentration of salt and methane dissolved in the liquid phase, θ w i is equal to one and θ w h = 1 − θ m h . The mass flux term of water in liquid (j w l ) is computed similarly as in Equation (9), while the flux terms in hydrate (j w h ) and ice phases (j w i ) are computed similarly as in Equation (11) for S h and S i , respectively .
The f w is 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: where θ st l is 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 (j st l ) 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 f st is 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: 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: with Please note that E m l and the dispersive terms (i E l and i E g ) 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: 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 f Q is the external supply of energy per unit volume. Table 4. Specific energy expression and representative values of specific and latent heat for each phase of the system.

Phase
Specific Energy Specific and Latent Heat

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: 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, P l ; gas pressure, P g ; temperature, T and salinity, S) and the dependent variables of the system (listed in Tables 5 and 6). These equations allow capturing the coupling among the various physical phenomena considered in the system and close the mathematical formulation. Table 5. Constitutive equations for MHBS system and dependent variables computed using each constitutive law. Please note that values for the parameters used in the TH simulation are also given.

Constitutive Law Equation Dependent Variable
Advective fluid flow Advective fluid flow, q l and q g Kozeny's model Gas density [42]

Stress-strain behavior
Effective stress [47] σ = σ − α B P p I Effective stress tensor, σ α B = 1 Hydrate-CASM [41] Subloading yield function Stress tensor, σ 0 † λ c is not computed in terms of S h or S i . This simplification valid for low saturations of ice [48]; ‡ m and ‡ P 0 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 e ah , κ h and p 0h recover the hydrate-free parameters e, κ and p 0 when (S h + S i ) = 0; The use of bold symbols denote vectors and tensors. Table 6. Equilibrium restrictions for MHBS system. The dependent variables computed using each of the equilibrium restrictions are also included.

Equilibrium Restriction Equation Dependent Variable
Hydrate phase change Hydrate phase boundary [43] ln(P eq ) = −1.644866 × Hydrate kinetic rate [52,53] R h (T, P g ) = φ p S h A h K d P eq (T) − P g − K f P g − P eq (T) Hydrate mass change, dm h

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).  [65,66,74] models. The simulations are plotted against the experimental data from [58] for (a) pore-filling and (b) cementing hydrate. Adapted from [41].

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., P g > P eq , 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 P g = 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 (ω m l ) 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., P g < P eq , 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.

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.

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.

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 (S h = 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 (S h = 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.  Figure 7. Simulated evolution of the (a,b) deviatoric stress, (c,d) potential porosity, (e,f) temperature, (g,h) liquid pressure and (i,j) hydrate and (j,k) gas saturation of the host and hydrate-bearing specimens subjected to thermal stimulation and depressurization during shear. Dotted lines represent the behavior of the host sediment specimens while solid lines refer to hydrate-bearing specimens. Please note that variations in hydrate and gas saturation before and after dissociation respectively, are due to the volumetric deformation of the specimen.
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).

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 S h 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.

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.

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

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: where dm h is 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: whereε is the sediment volumetric strain rate (sign convention as in continuum mechanics, negative compression):ε Replacing Equation (A3), Equation (A2) can be rearranged as: The first term in Equation (A4) expresses the effect of hydrate formation/dissociation reaction in the computation of dS h , 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 S h 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: where the term dm i ρ i φ p (1−S h ) corresponds to changes in ice saturation due to freezing and melting (computed using the freezing characteristic function proposed by [54]), and the term ∂(S h ) ∂t does only consider changes in hydrate saturation due to hydrate density variations.