A Geomechanical Model for Gas Hydrate Bearing Sediments Incorporating High Dilatancy, Temperature, and Rate Effects

Abstract

The geomechanical behavior of methane hydrate bearing sediments (MHBS) is influenced by many factors, including temperature, fluid pressure, hydrate saturation, stress level, and strain rate. The paper presents a visco-elastoplastic constitutive model for MHBS based on an elastoplastic model that incorporates the effect of hydrate saturation, stress history, and hydrate morphology on hydrate sediment response. The upgraded model is able to account for additional critical features of MHBS behavior, such as, high-dilatancy, temperature, and rate effects. The main components and the mathematical formulation of the new constitutive model are described in detail. The upgraded model is validated using published triaxial tests involving MHBS. The model agrees overly well with the experimental observations and is able to capture the main features associated with the behavior of MHBS.

1. Introduction

Methane hydrates are crystalline structures consisting of methane trapped in frozen-water molecules cages (e.g., [1,2]). They are stable at high pressure (P) and low temperature (T), and are generally found in seabed sediments and permafrost settings in the hydrate (P-T) stability zone. Methane hydrates bearing sediments (MHBS) are stable soils, however if changes in pressure and/or temperature (and/or chemical conditions) shift the hydrate outside the (P-T) stability zone, the methane-hydrate dissociates releasing gas-methane and liquid-water (e.g., [3,4,5,6,7]). According to Ruppel and Kessler [8], the global amount of methane carbon buried in MHBS is estimated to be close to 2000 gigatons.

The vast amount of methane hydrate hidden in both Arctic permafrost and marine continental margins worldwide has the potential to satisfy present and future global energy demand, provided efficient and economical production strategies are developed (e.g., Makogon et al. [9]). However, this is a challenging task owing to the complex behavior exhibited by methane hydrates upon dissociation (e.g., [6,10]). For example, the large volume changes related to hydrate dissociation (i.e., 1 m³ of hydrate produces during dissociation approximately 164 m³ of free gas and around 800 L of water) triggers significant increments in the fluid pressures, resulting in effective stress reduction, with the associated sediment deformations and changes in both soil porosity and permeability. Furthermore, because hydrate dissociation is a strong endothermic reaction, the reduction of the sediment temperature may freeze the pore water, blocking the sediment permeability. Moreover, methane production is generally accompanied with sand migration, which impact on both borehole mechanical stability and fluids flow. Sand production is considered as one of the major limitations for the commercial exploitation of gas hydrate (e.g., [11,12]). It is then apparent that the profound perturbations in the sediment condition and the multiphysics nature of this problem, with the strong couplings between the different thermo-hydro-mechanical and chemical (THMC) phenomena that control the sediment behavior require the development of advanced and robust models.

The interest in studying the behavior of MHBS does not only limit to energy production, but it also impacts on the stability of wellbores, offshore platforms, pipelines, and others subsea structures associated with hydrocarbon production. In addition, hydrate dissociation could trigger submarine landslides, which may affect underwater infrastructure. Furthermore, uncontrolled release of methane from hydrates will severely contribute to greenhouse effects and ocean warming (e.g., [8,13]). For example, after hydrate dissociation from submarine MHBS, the gas methane is not directly released to the atmosphere, but it is transformed (after biochemical processes in the ocean) into carbon dioxide, resulting in oxygen consumption and alterations in the chemical composition of the sea water (e.g., [14]). Therefore, the uncontrolled release of methane could severely harm our environment because the amount of carbon dioxide produced and released in this process is considerable.

One possible solution to safely produce energy from MHBS is via the exchange of methane by carbon dioxide in the hydrate cage (i.e., CO₂-CH₄ exchange). It has been proven in the laboratory that this replacement is energetically favorable and effective (e.g., [15,16,17]). The CO₂-CH₄ hydrate exchange will solve two problems simultaneously, i.e., stability of the hydrate bearing sediment during methane production, and capture of carbon dioxide in the sediment. However, the practical implementation of this technique in the field still requires further research and in-situ testing.

From the discussion above, we can conclude that a good understanding and modelling of MHBS is a critical component to address the challenges and opportunities associated with this type of sediment. The geomechanical behavior of MHBS is particularly challenging. Experimental studies in this area have primarily focused on how the strength, stiffness, and others mechanical properties of MHBS are affected by the following factors: hydrate saturation (e.g., [6,18,19,20]), effective mean stress (e.g., [21]), confining pressure (e.g., [22]), temperature (e.g., [23,24,25,26,27,28,29,30]), pore pressure e.g., [31]), strain rate (e.g., [32,33,34,35,36,37,38,39,40,41]), drainage conditions (e.g., [42]), and sediment skeleton (e.g., [20,43]).

Several THMC formulations and computer codes have been proposed to model the behavior of MHBS (e.g., [44,45,46,47,48,49,50,51,52,53]). In particular, several constitutive models have been proposed in the last few years to simulate the mechanical behavior of MHBS. The initial efforts in this area assumed and elastic behavior of the material. For example, the Duncan and Chang [54] model was proposed to formulate nonlinear elastic models for MHBS (e.g., [23,55,56,57,58]). An extended Mohr-Coulomb model (MCM) able to account for the influence of hydrate dissociation on strength (which is considered via a linear decrease in the cohesion) was proposed by Freij-Ayoub et al. [29] to evaluate wellbore stability. Extended MCMs that include the influence of hydrate saturation on both material strength and stiffness were also adopted to simulate the behavior of MHBS (e.g., [46,59,60,61]). More recent developments used critical state concepts ([62]) to model the response of MHBS. This type of approaches generally incorporates the effect of hydrate saturation, hydrate morphology, and other key MHBS factors into the mathematical formulation (e.g., [63,64,65,66,67,68,69,70,71,72]). The strains partition concepts proposed by Pinyol et al. [73] to model claystones was adopted by Sanchez et al. [68] to develop a model for MHBS able to evaluate the relative contribution of both sediment-skeleton and methane-hydrate to the mechanical behavior of the material. As for the effect of temperature on MHBS behavior, Yu et al. [28] propose a temperature dependent nonlinear elastic model and consider the effect of temperature on the initial modulus of the sediment. The upgraded Duncan-Chang model developed by Song et al. [23] adopts two coefficients to account for the effect of temperature and hydrate dissociation on the mechanical behavior of MHBS. The time-dependent behavior of hydrate bearing sediment have been mainly simulated using viscoelastic concepts (e.g., [37]). The THM formulations for MHBS proposed in Kimoto et al. [38] and Akaki et al. [39] adopted a mechanical model based on viscoplastic concepts, however the mechanical law was not validated. In Deusner et al. [40], the rate-dependent behavior of MHBS was explained considering a mechanism based on the kinematic rearrangement of the material fabric.

In this paper we present the extension of an existing model to deal with features of MHBS that have received less attention in the constitutive modeling of this type soil, such as, the influence of temperature, high dilatancy, and strain-rate effects. The starting point is the elastoplastic mechanical for MHBS proposed by Gai and Sanchez [67] which is able to account for the impact of methane-hydrates saturation on pre-consolidation stress, material stiffeners, and strength. The model was able to satisfactorily capture the overall MHBS experimental behavior observed in triaxial tests under different conditions. However, this model was not formulated to consider temperature and strain-rate effects. In this work we extend the original model (i.e., Gai and Sanchez [67]) to account for these two effects, and we also propose and enhanced hardening law to capture better the high dilatant behavior typically observed in MBHS. In the following sections we present first the mathematical formulation of the upgraded model, and we then evaluate the model performance when compared against published experimental results.

2. Methodology

We focused our efforts on the development of an advance geomechanical model for MHBS. The constitutive model is formulated in the general framework of viscoplasticity and adopts Perzyna’s theory [74] to extend an inviscid elastoplastic model to handle rate effects in soils. The proposed model is then validated based on available data in the literature. In this section, we present first some basic aspects associated with MHBS, we then briefly introduce the model proposed by Gai and Sanchez [67], afterward we extend it to handle soils exhibiting high dilatancy, non-isothermal conditions, and rate effects.

Methane hydrate saturation (S_H: ratio between the volume of hydrates and the volume of voids) has a strong influence on the mechanical response of MHBS. It has been shown that stiffness, peak deviatoric stress, and dilation of the MHBS specimens increase with S_H (e.g., [18,75]). However, the impact of methane-hydrates on sediment behavior not only depend on the amount of hydrates, but also on its morphology. Gas hydrates are generally found in three main form types in the sediment structure (Figure 1): (a) cementation, (b) pore-filling, and (c) supporting matrix (e.g., [1,76]). In the first pore-habit type, the methane-hydrates are mainly present at the contact between the grains and they act as a bonding material. In the pore filling form, the hydrates generally tend to grow freely in the pore space, without bridging particles together. Hydrates in the supporting matrix form are part of the solid skeleton. For a similar hydrate concentration, the cementing type hydrate morphology provides the maximum stiffness, strength and dilatancy (e.g., [18]).

In the following Section 2.1 we briefly present Gai and Sanchez [67] model and then its upgrade to handle high-dilatancy, temperature, and time-dependent behavior in Section 2.2, Section 2.3 and 2.4, respectively. In Section 3.1, Section 3.2 and Section 3.3 we present the model validation using published experimental data from triaxial tests.

2.1. Basic Model for MHBS—Brief Introduction

The starting point is the mechanical constitutive model for methane hydrate-bearing soils presented in Gai and Sanchez [67]. The model adopts a hierarchical single surface (HISS) proposed by Desai et al. [77] and Desai [78], and incorporates some key ingredients initially proposed by Uchida et al. [79] to deal with particular features of MHBS, namely: sub-loading concepts (e.g., [80,81]); cementing effects associated with the presence of hydrates; and bonding damage during shearing. The shape of the continuous HISS yield surface can be adapted to the particular conditions of the soil under investigation depending on the selected parameters. The HISS-MHBS model [67] is able to account for cementing effects, hydrate morphologies, and for the effect of confinement. The performance of the HISS-MHBS model was very satisfactory when compared against experimental data from triaxial tests based on natural core specimens and synthesized samples involving different hydrate saturations, hydrate morphology, and confinements [67]. However, the model underpredicted the volume expansion observed in tests involving hydrate bearing soils exhibiting high-dilatancy behavior. The model was not developed originally to handle the effect of both temperature and loading-rate on MHBS response. In the next sections we present the upgrade of the model to handle these relevant features of MHBS behavior.

The yield function of the HISS-MHBS model incorporating the strength enhancement effects (related to the presence of methane hydrate) and sub-loading concepts can be written as [67]:

$$F=\frac{a}{M^2}{q}^2-9\gamma p^{\prime 2}+9\gamma p^{\prime n}{\left[R\left(p_c+p_d\right)\right]}^{2-n}$$

(1)

where the constants a and $\gamma$ are related to the shape of the yield surface; n is a parameter associated with the transition from compressive to dilative volume change; p’ and q are the mean effective and deviatoric stresses, respectively; M is the slope of critical line in the q − p’ space; p_c is the effective pre-consolidation mean stress (which control the size of the elastic domain), and p_d controls the increase of the sediment strength associated with the presence of hydrates.

The evolution variable R (with 0 < R ≤ 1) is related to the sub-loading yield surface, which is introduced into the mathematical formulation to model: (i) irrecoverable strains that may develop when the stress state is inside the yield surface (aspect that cannot be modeled with a standard elasto-plastic model), and (ii) a smooth transition between elastic and plastic states. Figure 2 illustrates the three yield surfaces we consider in this model.

The hardening law is isotropic and depends on the plastic volumetric strains (ε_v^p):

$$\frac{d{p}_c}{p_c}=\frac{\left(1+e\right)}{\lambda -\kappa } d{\epsilon }_v^p$$

(2)

where e is void ratio; κ and λ are the slopes for the elastic and plastic isotropic paths in the e − (ln)p’ plane; respectively; and $d{\epsilon }_v^p$ is the volumetric plastic strain

The other equations related to model formulation are introduced in the Appendix A. In the Appendix we also explain how the basic model accounts for the effect of both S_H and hydrate morphology, more details can be found in Gai and Sanchez [67].

2.2. Model Upgrade to Consider MHBS Exhibiting High Dilatancy

Owing to the additional bonding provide by hydrates and the complex interactions between hydrates and soil skeleton, hydrate bearing soils generally exhibits high dilation. This feature of MBHS behavior is challenging to model (i.e., [67,82]). To overcome this issue the upgraded model incorporates an enhanced strain hardening law that depends on both, shear and volumetric plastic strains. The new hardening law extends Equation (2):

$$d{p}_c=\frac{1+e}{\lambda +\kappa }{p}_{c}d{ϵ}_v^p+D_s\frac{1+e}{\lambda +\kappa }{p}_{c}d{ϵ}_d^p$$

(3)

where $d{ϵ}_d^p$ is the deviatoric plastic strains; and $D_s$ is a parameter proposed by Nova [83] to control the effect of the deviatoric strain on hardening. This parameter varies between 0 and 1. The validation of this model is presented in Section 3.1.

2.3. Model Upgrade to Incorporate the Effect of Temperature on MHBS Behavior

The study of thermal effects on soil behavior has attracted growing attention in the last few years, particularly because of the increasing number of non-isothermal problems in energy geotechnics ([84]), such as, energy geostructures (e.g., [85,86]), and high-level nuclear waste disposal (e.g., [87,88,89]). It has been observed that temperature affects soil behavior in different manners, e.g., thermal volumetric strains depend on the stress history, the initial elastic modulus increases with the increase of temperature; and the preconsolidation pressure decreases with the temperature increase. It also appears that the friction angle at critical state and the normally consolidated line are independent of temperature. Campanella and Mitchell [90] proposed the first conceptual framework to consider the effect of temperature on soil behavior. Afterward, Hueckel and Baldi [91] adopted plasticity theory to model the thermomechanical behavior of saturated soils. Then, these ideas were extended by Gens [92] to model the behavior of unsaturated soils under non-isothermal conditions. Thermo-plastic models for soils were improved in subsequent research (e.g., [85,93]).

As discussed in the introduction, most of the analyses associated with the effect of temperature on MHBS have been based on thermal-dependent nonlinear elastic models (e.g., [23,28]). However, several laboratory investigations (e.g., [6,20,31,94,95,96]) focused on the behavior of MHBS in the hydrate-stability zone indicate that the thermal effects impact on their mechanical response beyond the elastic domain. In this section we extend the constitutive model presented in Section 2.1 to account for the effect of temperature in the MHBS hydrate stability zones. The extension of the HISS-MHBS model to account for thermal effects considers a dependence of the p_c (Equation (1)) on temperature, as suggested in Laloui and Cekerevac [93] for other type of soils:

$$p_c=p_{c0}\left(1-r_{T}log\frac{T}{T_0}\right)$$

(4)

where T₀ is the reference temperature, p_c0 is the preconsolidation mean stress at T₀, and r_T is a model parameter that considers the effect of temperature on the preconsolidation pressure. Figure 3 illustrates the effect of temperature on the MHBS yield. The validation of this model is presented in Section 3.2.

2.4. Model Upgrade to Consider Rate-Dependent Effects on the Behavior of MHBS

The time-dependent behavior of soils has been investigated for a long time. Particular emphasis has been placed in the understanding of creep phenomena in soils (e.g., [97,98,99]). The impact of strain-rate on shear behavior of soils has also been intensively studied (e.g., [100,101,102,103]). Rate effects in soils have been generally modeled with success using viscoplastic models (e.g., [104,105,106]). The interest of rate effects on the behavior of MHBS has increased in the last few years. Triaxial tests at different strain rates based on both artificial samples and natural specimens have been conducted by Miyazaki et al. [32,33,34,35,36] and Yoneda et al. [41], respectively. As for modeling, Miyazaki et al. [37] proposed a nonlinear viscoelastic constitutive equation and validated against tests conducted at different strain-rates on artificial methane-hydrate-bearing sand specimens. Kimoto et al. [38] and Akaki et al. [39] implemented a viscoplastic model in a coupled THM formulation to solve boundary value problems involving MHBS. However, the mechanical model was not validated in these works.

In this research we adopt the well-known Perzyna’s overstress theory [74] to extend the model presented in the previous sections to account for rate-dependent effects in MHBS. Perzyna’s concepts allow upgrading in a single manner elastoplastic constitutive models to include viscoplastic behavior. Perzyna’s approach has been successfully applied to develop several stress-strain rate-dependent models for soils (e.g., [107,108,109,110,111,112]). In this theory, stress states outside the yield surface are allowed and the distance between the yield surface (i.e., from the previous time step) and the predictor stress (associated with the new time/strain increment) is a measured of the rate of the viscoplastic strains.

Perzyna’s overstress concept is typically expressed through the following equation (e.g., [74,113]):

$$〈\Phi \left(F\right)〉=\left\{\begin{array}{c}\Phi \left(F\right) if F>0\\ 0 if F\le 0\end{array}\right.$$

(5)

where $〈∎〉$ is the Macaulay brackets, $\Phi$ is a scalar function (called sometimes the flux function [113]) that grows monotonically with F and defines the magnitude of the plastic strain rate.

The elastoplastic yield surface from the previous time-step is generally called the ‘static yield surface’ (F_S), and the homothetic yield surface passing through the current (predictor) stress state during yielding (i.e., outside F_S) is often called the ‘dynamic yield surface’ (F_D) (e.g., Hinchberger and Rowe [111]). In our model F_S and F_D are given by Equation (1), depending on whether we use the stresses and internal variables at the beginning of the time step (identified with the subscript S), or the stresses and internal variables associated with the predicted stresses (identified with the subscript D), respectively.

The direction of the viscoplastic strains is given by the gradient of the plastic potential function at the (current) predicted stress point (G_D) calculated with the following equation (assuming associated plasticity):

$$G_D=F_D=\frac{a}{M^2}{q}_D^2-9\gamma \left\{{\left(p_D^{\prime }\right)}^2-{\left(p_D^{\prime }\right)}^n{\left[R\left(p_c+p_d\right)\right]}^{2-n}\right\}$$

(6)

The visco-plastic strain rate is given by the following equation:

$$\frac{d{\epsilon }_D^{vp}}{dt}={\gamma }_f〈\Phi \left(F\right)〉\frac{\partial G_D}{\partial {\sigma }_D}$$

(7)

where ${\gamma }_{\mathrm{f}}$ is a model parameter called fluidity.

As for the flow function, we adopted the following expression, Desai and Zhang [114]:

$$\Phi \left(F\right)={\left(\frac{f}{f_0}\right)}^{n_f}$$

(8)

where f is the over-stress index; f₀ a reference value (i.e., such that the expression is non-dimensional); and n_f is a model parameter. Note that in this initial model we do not propose any dependence of the flow function on hydrate concentration, but it can be incorporated if needed. The over-stress index is calculated using the internal (static) variables at ‘S’ and the predicted (dynamic) stresses ‘D’, as follows:

$$f=\frac{a}{M^2}{q}_D^2-9\gamma \left\{{\left(p_D^{\prime }\right)}^2-{\left(p_D^{\prime }\right)}^n{\left[R_S\left(p_{cS}+p_{dS}\right)\right]}^{2-n}\right\}$$

(9)

and the viscoplastic strains are obtained after:

$$d{\epsilon }_D^{vp}=\frac{d{\epsilon }_D^{vp}}{dt}dt$$

(10)

The static and dynamic yield surfaces are updated during yielding. Figure 4 shows the yield surfaces adopted in this work. The validation of this model is presented in Section 3.3.

3. Results and Discussions

In the following sections we present the application and validation of the upgraded model based on published experimental results.

3.1. Model Application Involving MHBS Exhibiting Large Dilation

The first set of tests selected to validate the proposed model corresponds to the triaxial tests reported in Miyazaki et al. [82]. These tests were simulated in previous works using a model for MHBS based on the modified Cam-Clay model (e.g., Uchida et al. [64]) and the basic MHBS-HISS model, Gai and Sanchez [67]. These models were able to capture qualitatively well the overall behavior of the MHBS samples, but the high-dilatancy exhibited during shearing. The hydrate sediments were prepared in the lab using Toyoura sand at different hydrates saturations. The samples were then sheared at different confining stresses ($\sigma$₃). Figure 5, Figure 6 and Figure 7 present the experimental results (symbols) of the tests sheared at $\sigma$₃ = 1 MPa (S_H = 0%; S_H = 34%; and S_H = 41%), $\sigma$₃ = 2 MPa (S_H = 0%; S_H = 31%; and S_H = 43%), and $\sigma$₃ = 3 MPa (S_H = 0%, S_H = 27%, and S_H = 42%), respectively. These figures also present the modeling results obtained with the original MHBS-HISS model [67] (dash lines, identified as M_B), and the upgraded model considering Equation (3) for the strain-hardening law (solid lines, identified as M_U). The model parameters are listed in

Table A1. Parameters adopted for Model Validation in Figure 5, Figure 6 and Figure 7.

Properties	Upgraded Model MU	Basic Model MB
M	1.28	1.28
λ	0.2	0.2
κ	0.004	0.004
n	1.3	1.3
a	3	3
γ	−1/9	−1/5
α	35	20
β	1	1
μ	7	7
υ	0.1	0.1
Ds	0.8	0.0

. Note that p_c = 7.7 MPa, 12.1 MPa and 16.1 MPa were adopted for the cases associated with $\sigma$₃ = 1 MPa, 2 MPa, and 3 MPa, respectively.

In terms of deviatoric stress (i.e., Figure 5a, Figure 6a and Figure 7a), both constitutive models are able to reproduce quite satisfactorily the main trends observed in the tests, namely: increase of MHBS strength and stiffens with S_H, a more marked post-peak softening behavior with the increase of S_H, and increase of MHBS strength with the increase of confinement. The main difference between the two models is apparent when comparing the volumetric behavior (Figure 5b, Figure 6b and Figure 7b). The new model reproduces noticeably much better the large dilation observed in these tests, particularly for the MHBS samples with high S_H. The new model also captures well the decrease in soil dilation with the confinement increase.

The second case analyzed in this section corresponds to the triaxial tests reported in Hyodo et al. [6]. The synthetic MHBS samples were based on Toyoura sand and prepared at three hydrate saturations: S_H = 0%, S_H = 39%, and S_H = 54%. The samples were sheared at $\sigma$₃ = 3 MPa. We modeled these tests using both, the original MHBS-HISS model and the upgrade version proposed in this work.

Table A2. Parameters adopted for Model Validation in Figure 8.

Properties	Upgraded Model MU	Basic Model MB
M	1.28	1.28
λ	0.24	0.24
κ	0.005	0.005
pc (MPa)	9.6	9.6
n	0.95	0.95
a	3	3
γ	−1/10	−1/10
α	8	6
β	0.7	0.7
μ	4.5	4.5
υ	0.35	0.35
Ds	0.1	0.0

lists the model parameters. As in the previous case, the effect of hydrate saturation on MHBS behavior is characterized by a marked influence of S_H on material stiffens, strength (Figure 8a), and volumetric response (Figure 8b). It can be observed that the proposed model is able to reproduce more accurately the large dilation observed in these experiments and overcomes the shortcomings experienced by previous approaches to model this critical feature of MHBS behavior. These results suggest that the proposed model will be able to properly predict the sediment volume changes (with the associated settlements and ground subsidence) in engineering problems involving MHBS.

3.2. Model Application Considering the Effect of Temperature on MHBS Behavior

In this section we adopted the thermo-plastic model for MHBS introduced in Section 2.3 to simulate three cases reported in the literature. All the triaxial tests were conducted inside the hydrate stability zone at three constant temperatures. The first series of triaxial tests analyzed in this section were performed on synthetic MHBS samples prepared from Toyoura sand by Hyodo et al. [31]. The triaxial tests were conducted at two temperatures: T = 1 °C and T = 5 °C. The parameters adopted to model these tests are listed in

Table A3. Parameters adopted for simulation in Figure 9.

Parameter	Value
M	1.25
λ	0.25
κ	0.015
pc (MPa)	3
n	1
a	2
γ	−1/9
α	20
β	0.5
μ	5
υ	0.15

. Figure 9a shows that the upgraded model (lines) agreed very well with the observed experimental behavior, reproducing the increase of strength with temperature (although, the model slightly underpredicted it at the lower temperature), as well as the post-peak behavior observed in the experiments (symbols). The model is also capable of properly capturing the effect of temperature on the volumetric behavior of MHBS including the increase of dilation with the decrease of temperature (Figure 9b).

The second series of triaxial tests analyzed in this section correspond to the three synthetic MHBS samples based on kaolin clay reported in Li et al. [96]. The three tests were carried out at the same confining pressures ($\sigma$₃ = 5 MPa), same hydrate saturation (S_H = 30%), same loading rate, and they were shared at three constant temperatures, T = −5 °C, T = −10 °C, T = −20 °C.

Table A4. Parameters adopted for simulation in Figure 10.

Parameter	Value
M	1.2
λ	0.5
κ	0.02
pc (MPa)	5
n	1
a	6
γ	−1/9
α	15
β	1.2
μ	1
υ	0.15

lists the adopted parameters. The updated model captures well both the increase of strength and the change in material stiffness with temperature (Figure 10).

The last set of experiments analyzed in this section corresponds to the triaxial tests conducted by Hyodo et al. [6] on two MHBS Toyoura sand samples that were prepared at S_H = 53.1% and S_H = 51.6%, and sheared at T = 1 °C and at T = 10 °C, respectively. The two samples were tested at the same effective confining pressure ($\sigma$₃′ = 3 MPa). Figure 11a,b present the experimental results together with the modeling outputs. Figure 11a shows that the model predicts very well, the impact of temperature on soil stiffness, strength, and post-peak behavior. The volumetric behavior of the MHBS (Figure 11b) at different temperatures is also well capture by the upgraded model. The adopted parameters are listed in

Table A5. Parameters adopted for simulation in Figure 11.

Parameter	Value
M	1.25
λ	0.25
κ	0.015
pc (MPa)	3
n	1
a	3
γ	−1/9
α	20
β	1
μ	1
υ	0.15

.

3.3. Model Application Considering the Effect of Loading Rate on MHBS Behavior

To validate the capability of the upgraded model (Section 2.4) to deal with time-rate effects in MHBS we adopted two types of tests, namely, (a) drained triaxial tests conducted at different (constant) loading rates (LR), reported in Miyazaki et al. [33], and (b) drained triaxial experiments conducted at two loading rates during shearing, reported in Miyazaki et al. [33] and Yoneda et al. [41].

As for the set of tests at constant LR, the artificial MHBS samples were based on Toyoura sand and intended to have a target hydrate saturation S_H~40%. Given the difficulties to form hydrate in experimental conditions and to determine its saturation, in the examined samples S_H varies between 35~45%. Here we adopt a unique set of parameters (

Table A6. Model Parameters adopted for Case 1 in Figure 12 and Case 2 in Figure 13.

Parameter	Value
M	1.22
λ	0.2
κ	0.017
pc (MPa)	9
n	1
a	5
γ	−1/9
α	20
β	1
μ	20
υ	0.2
γf	5 × 10−6
nf	1

) to model all the experiments. Three loading rates were investigated, LR = 0.1%/min, LR = 0.05%/min, and LR = 0.01%/min. Three samples were tested per each experiment (i.e., per each S_H and LR). Figure 12a–c presents the comparisons between the experimental (symbols) and modeling (lines) results for the tests conducted at LR = 0.1%/min, LR = 0.05%/min, and LR = 0.01%/min, respectively.

The proposed model agrees well with the overall observed behavior involving different loading rates and hydrates saturations. Particularly, it properly captures the increase in strength with the increase of LR observed in MHBS samples. It also reproduces well the softening behavior after the peak in the MHBS samples, particularly at the two lower loading rates. There are some slight differences between experiments and simulations in some tests, but the overall performance of the model is very satisfactory, particularly considering that despite the variations in the S_H of the tested specimens, we adopted a unique set of parameters to model all the experiments.

In the first set of tests at varying LR during shearing analyzed in this section the artificial MHBS samples were also based on Toyoura sand, as reported in Miyazaki et al. [33]. The three tests involve samples prepared at two hydrate saturations, namely, S_H = 35% and S_H = 48%. The tests considered different loading rates and different axial strain intervals ($\Delta e$). This type of experiment generally contemplates the following test protocol: (i) preparation of the hydrate-sand or sand (S_H = 0%) specimen by the water saturation method in the test apparatus, (ii) application of the loading at a low loading rate (LR_L) until the axial strain reached a target value (e.g., 0.25%); (iii) modification of the loading rate to the high loading rate (LR_H) until the target axial strain interval (i.e., until e.g., another 0.25% is reached). This sequence is repeated successively until the total axial strain is reached (i.e., until shear failure). The three tests were conducted at the same LR_H = 0.1%/min but considering different low loading rates and axial strain intervals. In test (I) LR_L was 0.05%/min, $\Delta e$ = 0.5%, and S_H = 48% (Figure 13a); in test (II) LR_L was 0.05%/min, $\Delta e$ = 0.25%, and S_H = 48% (Figure 13b); and in test (III) LR_L was 0.01%/min, $\Delta e$ = 0.25%, and S_H = 35% (Figure 13c).

We modeled these tests following the same test protocol explained above. We adopted a unique set of parameters to simulate all the experiments, which are listed in Table A6. Figure 13a–c presents the tests (blue lines) and the numerical outputs in red lines. The two stress-strain curve envelops related to LR_H (curve 2) and LR_L (curve 1) are plotted using green dash lines in Figure 13a.

Test (III) (i.e., Figure 13c) exhibited a slightly lower strength that tests (I) (i.e., Figure 13a) and (II) (i.e., Figure 13b). This response was expected because in test (III) both S_H and LR_L were smaller than in tests (I) and (II). Comparing Figure 13b, c, one can conclude that (for the same $\Delta e$), the larger the difference between the LR_H and LR_L, the larger the difference between the deviatoric stresses associated with these loading rates. Comparing Figure 13a, b, we can see that (for the same loading rates and hydrate saturation), the change in $\Delta e$ does not significantly impact on the deviatoric stresses associated with the LR_H and LR_L. The proposed model is able to capture qualitatively well the main features of MHBS behavior discussed above. The results obtained with the viscoplastic model are also satisfactory in quantitative terms.

The final set of experiments analyzed in this section corresponds to the two drained compression triaxial tests conducted by Yoneda et al. [41] at varying LR involving natural samples from different boreholes (Borehole 16B-4P and Borehole 17C-9P). Both tests adopted the same axial strain interval ($\Delta e$ = 1%), but in one of them LR_L was 0.01%/min, LR_H = 0.05%/min, and S_H~80% (Figure 14a); and in the other test LR_L was 0.001%/min, LR_H = 0.01%/min, and S_H~65% (Figure 14b).

In this case it is also apparent that the higher strength corresponds to the test conducted at the higher loading rate and on the specimen with the higher hydrate saturation (Figure 14a). A tendency to dilate was also observed in the first test that corresponds to the specimen with the higher S_H (Figure 14b). The parameters adopted to model these two experiments are listed in

Table A7. Model parameters adopted for Case 3 in Figure 14.

Parameter	Borehole 16B-4P	Borehole 17C-9P
M	1.2	1.2
λ	0.28	0.28
κ	0.025	0.025
pc (MPa)	10	6.8
n	1	1
a	5	4.3
γ	−1/12	−1/12
α	10	10
β	10	10
μ	5	5
υ	0.2	0.2
γf	2 × 10−6	2 × 10−6
nf	1	1

. In addition, in this case, it can be observed that the model captures well the global behavior observed in these experiments, reproducing very satisfactorily both the stress developed at the different loading rates and the volumetric behavior displayed by these specimens during shearing. These are critical features of MHBS behavior that the proposed approach is able to properly reproduce, enabling a more reliable modeling of engineering problems involving this type of soil.

4. Conclusions

We upgraded a mechanical constitutive model for MHBS to account for features of this type of soil that are generally not well simulated by standard models. We selected an existing model that considers, amongst others factors, the effect of hydrate saturation on material stiffness, strength and pre-consolidation stress. We incorporated into the formulation three main new components: an enhanced hardening law (to capture the large dilatancy behavior observed in MHBS); a modified evolution law for the pre-consolidation pressure (to account for the effect of temperature on the sediment mechanical properties); and Perzyna’s concepts to include rate effect on the response of MHBS. We then validated the model using published data associated with the testing of methane hydrate bearing soils under different conditions. The proposed approach agreed well with the overall response observed in the selected experiments, demonstrating that the new model is able to successfully account for key features of MHBS behavior, such as, effect hydrate saturation and morphology, large dilatancy, thermal effects, and rates effects. These are key features of material behavior to properly model different engineering problems involving MHBS behavior, amongst others, wellbore geomechanical stability and seabed subsidence during methane production from hydrate bearing sediments reservoirs; and submarine landslides. Furthermore, as MHBS behavior engages strongly coupled THMC phenomena, where the mechanical behavior controls the changes in material porosity (with its direct impact on hydraulic permeability, fluids storage, thermal conductivity, and diffusion) a geomechanical constitutive able to predict the actual material behavior is instrumental to conduct reliable numerical simulations.