A State-Dependent Elasto-Plastic Model for Hydrate-Bearing Cemented Sand Considering Damage and Cementation Effects

Abstract

In order to optimize the efficiency and safety of gas hydrate extraction, it is essential to develop a credible constitutive model for sands containing hydrates. A model incorporating both cementation and damage was constructed to describe the behavior of hydrate-bearing cemented sand. This model is based on the critical state theory and builds upon previous studies. The damage factor D_s is incorporated to consider soil degradation and the reduction in hydrate cementation, as described by plastic shear strain. A computer program was developed to simulate the mechanisms of cementation and damage evolution, as well as the stress-strain curves of hydrate-bearing cemented sand. The results indicate that the model replicates the mechanical behavior of soil cementation and soil deterioration caused by impairment well. By comparing the theoretical curves with the experimental data, the compliance of the model was calculated to be more than 90 percent. The new state-dependent elasto-plastic constitutive model based on cementation and damage of hydrate-bearing cemented sand could provide vital guidance for the construction of deep-buried tunnels, extraction of hydrocarbon compounds, and development of resources.

1. Introduction

Natural gas hydrates are seen as a novel and ideal energy source with environmental value in the twenty-first century because of their wide distribution, vast geologic reserves, high energy density, and minimal environmental pollution [1,2]. Natural gas hydrates are crystals of water and natural gas, primarily methane, and the detected reserves are widely distributed [3,4] in marine sediments and permafrost on Qingdao Plateau along the Sichuan-Tibet Railway [5]. The wetland structure in the gas hydrate-rich region of the Tibetan Plateau has also been anthropogenically damaged by the gas hydrate exploration activities that have recently taken place in China’s onshore area [6,7,8,9]. The microstructure and elastic properties of the rocks [10,11,12,13,14] in sedimentary strata are altered by the presence of gas hydrates. Natural gas hydrates decompose during extraction, weakening the mechanical qualities of the strata that contain them. This causes deformation and destroys the reservoirs, which in turn causes major geological disasters like landslides, stratigraphic deformation, and collapse.

The effect of hydrates on the mechanical properties of soils is primarily based on both filling and cementation, which results in energy soils exhibiting properties similar to those of dense sandy or cemented soils [15,16,17]. According to research, hydrates are primarily distributed in four ways in natural energy soils: pore filling, cementation, wrapping, and skeletal support [18,19,20]. According to Brown et al. [21], the breakdown of hydrates weakens the cementation of soil particles, increases porosity, and results in underconsolidated or loose sand soil. The development of a plausible constitutive model of hydrate-bearing cemented sand is of great value and significance for optimizing the efficiency and safety of natural gas hydrate exploitation.

To study the mechanical properties of hydrate-containing cemented sand, numerous researchers have produced hydrate-containing soils indoors [22,23,24,25,26]. By putting methane gas into sand with varying water contents and applying temperature and pressure conditions, the effects of various saturations, confining pressures and temperatures on the shear strength of energy soils were examined [27]. The gas saturation method was used to produce test specimens. As gas hydrate dissociates, the shear strength gradually decreases from 30% to 90% of the initial level [28], the shear zone thickens and the soil displays a brittle damage pattern [17]. Soil with more dispersed distribution of hydrates yields later, has higher peak strengths and slower development of cementation damage [29].

Experts and academics both domestically and internationally have conducted a number of studies on the constitutive model of hydrate-containing cemented soils, which may be broadly divided into three groups based on the experimental results: The nonlinear elastic constitutive model is the first kind. Using the least-squares method, Miyazaki et al. [22] determined the tangential Poisson’s ratio related to the matrix material and defined the tangential modulus expression related to hydrate saturation (S_h). The effective circumferential pressure and hydrate saturation were taken into account when establishing the nonlinear elasticity principle model. In order to improve the model’s ability to represent the stress-strain relationship of naturally occurring soft clay soil, Wang et al. [30] introduced the concept of soil damage to modify the Duncan-Chang model. The second type is a damage constitutive model for hydrate-bearing cemented soils based on geotechnical damage theory [31,32,33,34]. The use of effective hydrate saturation was proposed to characterize the effect of hydrate storage mode [35,36]. Weibull distribution and residual strength correction coefficients were introduced into the damage statistical constitutive model [37]. Nevertheless, these models do not take into account the impact of gas hydrate cementation or alternative storage methods, and obtaining the Weibull parameters might be challenging. The third category is the elastic-plastic constitutive model. By considering the filling effect, cementation effect and other forms of fugacity of hydrates, the existing models such as the Cam-Clay model [24] and the uniform hardening model [38] are expanded accordingly. Uchida et al. [26] took into account the shear expansion and cohesion of the soil, introduced the concept of mechanical hydrate saturation related to the skeleton structure of the sediments, and transformed the yield surface equation of the modified Cam-Clay model to model the critical state of hydrate-bearing soils. Considering that hydrates have no effect on the residual strength of soils, parameters related to hydrate saturation were introduced to correct the hardening law [25]. Zhu et al. [39] introduced cementation stresses s in the uniform hardening model of structural soils under the basic framework of uniform hardening theory. Two parameters $\beta$ (controlling the degree of strength of the cemented structure) and the parameter $\omega$ (controlling the degree of breakage of the cemented structure) related to the structural nature of the cementation were introduced into the elliptic-parabolic double yield surface model [40]. The influence of cementation was considered by translating the yield surface and the critical state surface based on the double hardening parameter model [41]. Li et al. connected the void ratio and the stress ratio to derive an equation of state-dependent dilatancy [42,43]. The critical state model of super-consolidated structural soils was revised in terms of an expression for the law of evolution of the structural strength p to represent the process of cemented structural damage of the soil [44].

A thorough knowledge of the damage and deformation features of sand affected by hydrates and weakened by hydrate cementation remains elusive after years of intensive research on hydrate-bearing cemented sand. The evolution law of the cementation stress gradual exertion process is not described by the majority of previously built models; instead, they explain the evolution law of the soil cementation stress decay process. Furthermore, the damage factor was not taken into account in the previously proposed hydrate-bearing cemented sand, making it impossible to adequately depict the deterioration of the cementation and the soil body damage. For hydrate-bearing cemented sand, we thus require a new state-dependent critical state model that takes cementation and damage into account.

An elastic-plastic model of hydrate-bearing cemented sand that takes cementation and damage into account was developed using the critical state theory, building on earlier research. To account for soil damage and the weakening of hydrate cementation, the damage factor D_s, which is defined by plastic shear strain ${\epsilon }_{\mathrm{d}}^{\mathrm{p}}$ is added. The damage is then incorporated into the cementation strength function P_h to describe the evolution law of the cementation stress progressive exertion process. By altering the shear dilatancy equations, yield surface equations, and hardening parameters of the constitutive model, a computer program was created to simulate the mechanisms of cementation and damage evolution as well as the stress-strain curves of hydrate-bearing cemented sand under various hydrate saturations S_h, effective perimeter pressures CSP, and initial porosities e₀. The results demonstrate the effectiveness of the model in reproducing the constitutive mechanical behavior of impaired soil cementation and soil deterioration with the same set of parameters when compared to the experimental data of earlier researchers. The conclusions can offer crucial directions for building deep-buried tunnel construction, hydrocarbon compound mining, and resource development.

2. Model Formulation

Soil was first treated as a work-hardening plastic material by Drucker et al. in 1957 [45]. Roscoe et al. found that the final portion of all loading paths lies on a unique surface and that the paths end with a unique critical porosity ratio line. The presentation of the Cam-Clay model and the modified Cam-Clay model marked the formalization of critical state theory [46,47]. And for the first time, strength theory and deformation theory (compression, dilatancy) were integrated into one framework. Roscoe et al. found that the deformation development of soil under external loading ends at a particular point regardless of the initial state and stress path. Soil in the large deformation stage of shear test tends toward the final critical condition, that is, the volume and total stress are unchanged, while the shear strain continues to develop and flow.

$$\stackrel{.}{\mathsf{\sigma }}=0,{\dot{\epsilon }}_v=0,{\dot{\epsilon }}_q\ne 0$$

(1)

Since the solid particles are incompressible, the change in the volume of the soil specimens depends entirely on the change in the voids, i.e., the physical significance of the description of the soil specimens volume v and porosity e is equivalent. The effective stress in the critical state is constant, it is equivalent to the ratio of q and p* being constant.

$$e=e_c={\widehat{e}}_c(p)$$

(2)

$$\eta ={(q/p^{*})}_c=M_c$$

(3)

The drainage conditions and the state of the soil have no bearing on the requirements for attaining the critical state. Dimaggio and Sandler proposed the cap model based on Drucker et al.’s study and the Cambridge model [48,49]. A straightforward and useful method for easily obtaining the cap model coefficients from traditional geomechanical test data was presented by Huang and Chen (1990) [50].

When the soil is inside a yield surface, the model predicts an elastic behavior; after the yield surface is reached, a plastic behavior starts. The p*-q space is used for establishing the model. Deviator stress and mean effective stress are described as

$$p*=\frac{\left({{\sigma }_1}^{*}+2{{\sigma }_3}^{*}\right)}{3}$$

(4)

$$q={\sigma }_1^{*}-{\sigma }_3^{*}$$

(5)

where ${{\sigma }_1}^{*}$ and ${{\sigma }_3}^{*}$ are the major and minor primary stresses.

For $p*$ and q, the calculated conjugate strain rates are $d{\epsilon }_v$ and $d{\epsilon }_d$.

$$d{\epsilon }_v=d{\epsilon }_1+2d{\epsilon }_2$$

(6)

$$d{\epsilon }_d=\frac{2\left(d{\epsilon }_1-d{\epsilon }_3\right)}{3}$$

(7)

$$d\epsilon =d{\epsilon }_v+d{\epsilon }_d$$

(8)

$$d{\epsilon }_d=d{\epsilon }_d^e+d{\epsilon }_d^p$$

(9)

$$d{\epsilon }_v=d{\epsilon }_v^e+d{\epsilon }_v^p$$

(10)

where $d{\epsilon }_d^e$ and $d{\epsilon }_d^p$ are the elastic and plastic shear strain increments, and $d{\epsilon }_v^e$ and $d{\epsilon }_v^p$ are the elastic and plastic volumetric strain increments.

2.1. Implication of Damage

Mining leads to increased strain in soils, as well as performance degradation from hydrate cementation. Damage theory is applied to represent the degradation of hydrate cementation and the damage to cemented sand. It is hypothesized that the equation between the plastic shear strain ${\epsilon }_{\mathrm{d}}^{\mathrm{p}}$ and the damage variable D_s of hydrate cementation and soils is as follows [51]:

$$D_s=1-{\mathrm{e}}^{s{\epsilon }_{\mathrm{d}}^{\mathrm{p}}}$$

(11)

where s is the damage index, reflecting the damage rate of hydrate-bearing cemented sand.

2.2. Bonding and Debonding Regulations of Hydrate-Bearing Cemented Sand

The hydrates can increase the strength and cohesive strength of soils, giving them properties similar to structural and cemented sand [52]. In the critical state, the bonding is gradually weakened until it is completely destroyed, and only the friction strength plays a role [53,54].

The expression for the cementation factor P_h0 is referenced [26]:

$$p_{\mathrm{h}0}=c{(S_h)}^g$$

(12)

Due to the very limited volumetric strain, the bonding damage D_s is introduced into the bonding strength function P_h, which describes the evolution law of the gradual degradation of the bond stress [24]:

$$\begin{array}{l}d\ln p_{\mathrm{h}}=-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p\\ \ln p_{\mathrm{h}}=-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p+C\\ p_{\mathrm{h}}=\mathrm{C}\cdot {\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\\ p_{\mathrm{h}0}=c{(S_h)}^g\\ p_{\mathrm{h}}={\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\cdot c{(S_h)}^g\end{array}$$

(13)

where v_c, c and g are model parameters. Introducing the damage factor D_s yields

$${p^{*}}_{\mathrm{h}}=p_{\mathrm{h}}\cdot (1-D_s{)=\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\cdot c{(S_h)}^g\cdot (1-D_s)$$

(14)

where ${p^{*}}_{\mathrm{h}}$ is the bonding strength function considering damage.

2.3. Elasticity Stage

The definition of the hydrate content in cemented sand is as listed below:

$$S_h=\frac{V_h}{V_h+V_0}$$

(15)

where V₀ and V_h stand for the pore space and hydrate volumes.

Based on the concept of elastic theory:

$$d{\epsilon }_{\mathrm{d}}^e=\frac{dq}{3G}$$

(16)

$$d{\epsilon }_v^e=\frac{d{p}^{*}}{K}$$

(17)

$$K=\frac{2(1+\mu )}{3(1-2\mu )}G$$

(18)

where G and K are the elastic shear and bulk parameters, $\mu$ is the Poisson’s ratio. Applying the hydrate-bearing cemented sand damage index D_s results in

$$d{\epsilon }_{v{D}_s}^e=\frac{d{p}^{*}}{K(1-D_s)}$$

(19)

$$d{\epsilon }_{\mathrm{d}{D}_{\mathrm{s}}}^e=\frac{dq}{3G(1-D_s)}$$

(20)

where $d{\epsilon }_{v{D}_s}^e$, $d{\epsilon }_{\mathrm{d}{D}_{\mathrm{s}}}^e$ are the increments of elastic volumetric and shear strain considering damage D_s.

The shear stiffness resulting from the presence of hydrate and the sand skeleton’s shear modulus are added together to form G. It is presumed that the increase in shear stiffness $G_h\left(S_h\right)$ increases linearly with hydrate saturation [26].

$$G_h\left(S_h\right)=G_h+\tau S_h$$

(21)

where G_h(S_h) is the shear stiffness index for hydrate-bearing cemented sand, G_h is the shear stiffness index for sand particles, and $\tau$ is a model parameter.

In Figure 1, $\tau$ represents the gradient of the simulation line of experimental data from Clayton et al. [55]. Equation (21) illustrates the cumulative effect of the hydrate hole-filling and cementation processes on the hydrate-bearing cemented sand stiffness. It is impossible to ascertain the individual contributions of each mechanism to the rise in stiffness, though.

For sand particles, G can be represented using the empirical formula that follows [24]:

$$G=\frac{G_{\mathrm{p}}\cdot \sqrt{p^{*}{p}_0}\cdot {(2.97-e)}^2}{1+e}$$

(22)

where P₀, G_p indicates the atmospheric pressure and shear stiffness parameter.

2.4. Yield Function Considering the Effects of Cementation and Damage

In this paper, the cementation factor ${p^{*}}_{\mathrm{h}}$ and damage factor D_S are introduced to characterize the weakening of hydrate cementation and sand damage. Prior to grain breaking, a steady stress-ratio loading only causes comparatively minor plastic volumetric strains [56]. The yield function for hydrate-bearing cemented sand could be expressed as

$$f=q-{\eta }^{*}\left(p^{*}+{\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\cdot c{(S_h)}^g\cdot (1-D_s)\right)$$

(23)

where ${\eta }^{*}$ is the stress ratio at yielding. Figure 2 shows the yield criterion for hydrate-bearing cemented sand in the p*-q space.

According to the theory of plasticity [57], the plastic consistency requirement, which is stated as follows, can be used to determine a loading index L

$$L=\frac{1}{K_p}\left(\frac{\partial f}{\partial p^{*}}d{p}^{*}+\frac{\partial f}{\partial q}dq+\frac{\partial f}{\partial {p^{*}}_{\mathrm{h}}}d{p_{\mathrm{h}}}^{*}\right)=\frac{1}{K_p}\left(dq-{\eta }^{*}d{p}^{*}-{\eta }^{*}{\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\cdot c{(S_h)}^g\cdot (1-D_s)\right)$$

(24)

where K_p is the parameter of plastic hardening.

The $d{\epsilon }_v^p$ and $d{\epsilon }_{\mathrm{d}}^p$ can then be expressed as [42]:

$$d{\epsilon }_{\mathrm{d}}^p=L=\frac{1}{K_p}\left(dq-{\eta }^{*}d{p}^{*}-{\eta }^{*}{\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\cdot c{(S_h)}^g\cdot (1-D_s)\right)$$

(25)

$$d{\epsilon }_v^p=LD=\frac{D}{K_p}\left(dq-{\eta }^{*}d{p}^{*}-{\eta }^{*}{\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\cdot c{(S_h)}^g\cdot (1-D_s)\right)$$

(26)

where D is the dilatancy equation.

The damage factor D_s is introduced and at the same time the effective stress is changed. Based on the assumption of equivalent strain, it can be obtained:

$$d{\epsilon }_v=d{\epsilon }_{v{D}_s}^e+d{\epsilon }_{v{D}_s}^p=\frac{d{p}^{*}}{K(1-D_s)}+\frac{D}{K_p(1-D_s)}\left(dq-{\eta }^{*}d{p}^{*}-{\eta }^{*}\cdot {\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\cdot c{(S_h)}^g\cdot (1-D_s)\right)$$

(27)

$$d{\epsilon }_q=d{\epsilon }_{d{D}_s}^e+d{\epsilon }_{d{D}_s}^p=\frac{dq}{3G(1-D_s)}+\frac{1}{K_p(1-D_s)}\left(dq-{\eta }^{*}d{p}^{*}-{\eta }^{*}\cdot {\mathrm{e}}^{-v_{\mathrm{c}}d{\epsilon }_{\mathrm{d}}^p}\cdot c{(S_h)}^g\cdot (1-D_s)\right)$$

(28)

2.5. Critical State

Hyodo et al. [58,59] examined the influence of density and hydrate saturation on the mechanical characteristics of the hydrate-bearing cemented sand. Additionally, it was presumable that each drained triaxial test ended with the approximate achievement of the critical state.

The influence of hydrate saturation on the critical state lines is depicted in Figure 3. It is evident that the average effective stress and the void ratio have a roughly linear relationship. As value of S_h increases, the critical state lines moves upward, indicating that the hydrate-bearing cemented sand exhibits a tendency to dilate.

The impact of S_h on the critical state lines is shown in Figure 4. The data from the experiment [59] matches well to a number of straight lines, the slope of which rises as one increases the value of S_h.

The change in stress ratios with hydrate saturation M(S_h) at the critical state is shown in Figure 5. M(S_h) rises in a nonlinear manner as hydrate saturation rises. Because there is more contact between the hydrate and the sand particles, more hydrate will add to sand-hydrate frictional resistance, the critical stress ratios increase more quickly as value of S_h approaches 60%.

It appears that the critical stress ratio M(S_h) and critical state lines for hydrate-bearing cemented sand change with the hydrate saturation based on the earlier studies [24].

$$e_{\varpi }\left(S_h\right)=e_{\varpi }+\varsigma S_h^{\kappa }$$

(29)

$$e_c\left(S_h\right)=e_{\varpi }\left(S_h\right)-{\gamma }_c\ln \left(\frac{p^{*}}{p_0}\right)$$

(30)

$$M\left(S_h\right)=M+\chi S_h^{\delta }$$

(31)

where M and M(S_h) are the critical stress ratio for sand particles and the hydrate-bearing cemented sand, and $e_{\varpi }$, ${\gamma }_c$, $\varsigma$, $\kappa$, $\chi$ and $\delta$ are model critical state parameters.

2.6. State-Dependent Dilatancy Function

Tests have also shown that the dilatancy of the hydrate-bearing cemented sand is dependent upon bonding strength, hydration saturation, and stress ratio [60,61,62]. A general formulation for the dilatancy was put forward by Li and Dafalias [43].

$$D=D(\eta , e, R, G)$$

(32)

where R is an internal state that is changeable beyond the void ratio, G stands for the intrinsic material constant, and $\eta$ is the stress ratio.

There are two conditions that the generic dilatancy statement must fulfill. First of all, in a critical state, the dilatancy must be zero. Second, there is a chance to arrive at a “phase transformation state”, where e ≠ e_c, $\eta$ ≠ M and D = 0.

$$D\left(R,G,e\ne e_c,\eta \ne M\right)=0$$

(33)

$$D\left(R,G,e=e_c,\eta =M\right)=0$$

(34)

In order to account for the hydrate accumulation behavior in the dilatancy of sand, the following modified dilatancy expression is suggested (24):

$$D=D\left({\eta }^{*},e,{p^{*}}_{\mathrm{h}},R,G\right)$$

(35)

where p_h stand for the impacts on dilatancy of the cementing of hydrates.

Been and Jefferies [63] introduced a state factor $\varsigma =e-e_c$ to quantify the difference between the present situation and the critical state condition of sand particles by measuring the density. Li and Dafalias [43] introduced the state factor $\varsigma$ to express the state-dependent for the dilatancy equation.

$$D=\frac{n_0}{M}\left(M{e}^{b\varsigma }-\eta \right)$$

(36)

where b and n₀ are two model parameters.

In order to fulfill the conditions of Equations (33) and (34), the dilatancy equation is modified [24].

$$D=\frac{n_0}{M\left(S_h\right)}\left[M\left(S_h\right)e^{b\varsigma \left(S_h\right)}{n}_c-{\eta }^{*}\right]$$

(37)

$$\varsigma \left(S_h\right)=e-e_{\varpi }\left(S_h\right)$$

(38)

where $\varsigma \left(S_h\right)$ is the state index for the hydrate-bearing cemented sand. n_c is an equation of the bonding strength ${p^{*}}_{\mathrm{h}}$.

Considering the hydrate cementation weakening and soil damage as reflected in the bonding strength ${p^{*}}_{\mathrm{h}}$. The damage factor D_s is introduced to propose the formulation for n_c. n_c is employed to explain the cementing process on dilatancy.

$$n_c=1-s\cdot D_{\mathrm{c}}$$

(39)

2.7. Hardening Rule

A significantly revised constitutive equation is provided for the plastic index K_p, as suggested by Li and Dafalias [43]. For a wide variety of confining pressures and void ratios, this kind of equation could effectively reflect the strain-hardening and strain-softening properties of hydrate-bearing cemented sand [43,64]:

$$K_p=\frac{lG{e}^{j\varsigma \left(S_h\right)}}{{\eta }^{*}}\left[M\left(S_h\right)e^{-j\varsigma \left(S_h\right)}-{\eta }^{*}\right]$$

(40)

where j is a hardening rule index and l is a factor of the plasticity.

$$l=l_1-l_2{e}_0$$

(41)

where l₁ and l₂ are model parameters. The plastic modulus is employing l to represent the impact of density.

3. Calibration of the State-Dependent Elasto-Plastic Model Parameters

There are nineteen parameters in the model. Ten of these are substrate sand mechanical properties parameters, which are listed in

Table 1. Model parameters of the substrate sand.

Symbol	Physical Meaning of Parameters	Value (Dimensionless)
$\mu$	Poisson’s ratio	0.05
Gp	The shear modulus factor	300
${\gamma }_c$	Model parameter, the slope of lines of critical state	0.127
M	Critical stress ratio for substrate sand	1.16
$e_{\varpi }$	Model parameter for void ratio of substrate sand	1.119
n0	Model parameter for dilatancy	1.62
b	Model parameter for dilatancy	1.31
l1	The plastic modulus parameter	1.96
l2	The plastic modulus parameter	2.23
j	Model hardening rule parameter	0.78

and can be calibrated using the techniques suggested by Xiao et al. [54] and Li and Dafalias [43]. These parameters were derived by data fitting or formulae. Moreover, the calibrated values of these parameters are ratios, which are dimensionless, and are used to model and accurately predict the stress-strain relationships of hydrate-bearing cemented sand.

We thank the reviewer for this criticism. In Figure 1, $\tau$ represents the gradient of the simulation line of experimental data from Clayton et al. [55]. In Figure 3, ${\gamma }_c$ denotes the slopes of void ratio e and ln(p*/p₀). The values can be found by fitting the experimental data [59]. The parameter $e_{\varpi }$ of Equation (30) could be calibrated by fitting the test data for the critical stress ratio and the critical state line in the e-p* plane [43]. In Figure 5, M denotes the critical state stress ratio, which can be obtained from the slope of the deviatoric stress q and mean effective stress p* curve fitted to the experimental data [59]. The parameter n₀ could be calibrated based on the ${\epsilon }_v-{\epsilon }_q$ curves. Based on the experimental q-${\epsilon }_q$ curves, a suitable value of $\mu$ was firstly chosen and then the parameter G_p can be calculated. Once G_p has been determined, l, and therefore l₁ and l₂, could be calibrated on Equation (41) [43]. For hydrate-bearing cemented sand, nine more parameters are required (see

Table 2. Extra model parameters for hydrate-bearing cemented sand.

Symbol	Physical Meaning of Parameters	Value (Dimensionless)
$\tau$	An elastic modulus parameter denoting the gradient of the simulation line in Figure 1	215.7
$\varsigma$	The critical state model parameters	0.58
$\kappa$	The critical state model parameters	1.99
$\chi$	The critical state model parameters	0.36
$\delta$	The critical state model parameters	1.38
g	Bonding parameter	1.59
c	Bonding parameter	0.67
vc	Debonding model parameter density on the plastic modulus	7.89
s	Damage factor	−0.10

); these are calibrated using the methodology described in the Shen et al. publication [24].

The damage index s represents the rate of damage to hydrate-bearing cemented sand. Since it is not possible to determine the damage index s directly through testing, this study examines the effects of various values of s on the simulation curves of hydrate-bearing cemented sand using the results of deviatoric stress q—axial strain ${\epsilon }_1$ curves that were predicted with an efficient restriction pressure of 5 MPa and a hydrate saturation degree of S_h = 24.2%.

In the simulation process, the s value range is between 0 and −0.1. The s values are closer to the test results when 0.1. The s value can be adjusted according to the impact of different external factors on the rate of degradation of hydrates [51]. An analysis of Figure 6 demonstrates that the higher the value of s, the higher the peak intensity of the deviatoric stress-axial strain curve. The lower the value of s, the more evident the cementation’s weakening effect is, resulting in a smaller peak value overall. When the axial strain is small, the impact of the model parameter s on the q-${\epsilon }_1$ curve is small, and with an increase in axial strain, the impact on the q-${\epsilon }_1$ curve begins to appear.

4. Model Validation and Analysis

Utilizing the host material Toyoura sand for the drainage triaxial test, cylindrical frozen specimens with a diameter of 30 mm and a height of 50 mm were created by combining the volume of sand corresponding to the target density with the amount of water in the range of 0–60% saturation of the target hydrate [59]. Triaxial compression tests were conducted using initial porosities of around 40% and 45% at effective confining pressures of 1 MPa, 3 MPa, and 5 MPa and temperatures of 1 °C, 5 °C, and 10 °C at a shear rate of 0.1%/min.

The deviatoric stress and volumetric strain actual and theoretical curves for substrate sand and cemented sand with 53.7% hydrate saturation at various effective confining pressures are displayed in Figure 7 and Figure 8. Deviatoric stress q-${\epsilon }_1$ and columetric strain ε_v-${\epsilon }_1$ relationships predicted by the state-dependent elasto-plastic model broadly match the experimental data [59], where a significant increase in perimeter pressure led to a corresponding increase in the stiffness and strength of cemented sand notably, and the peak intensity of the curve increases with increasing perimeter pressure. Because of the presence of bond strength and an increase in the critical stress ratio, the peak strength of hydrate-bearing cemented sand is higher than that of substrate sand. It is theorized that there is a small discrepancy in the model calculation when the perimeter pressure is lower since the predicted curve peak strength is smaller than the test curve. The stress-strain curve of the cemented sand develops towards strain hardening as the perimeter pressure increases, as evidenced by the continual saturation of hydrate and the disappearance of the softening tendency with the gradient of the perimeter pressure. The inclusion of the damage component in the model, which accounts for both the damage and the breakdown of hydrate, is likely the reason why the peak intensity of the simulated curve is marginally lower than that of the experimental curve.

For substrate sand with the same initial porosity and hydrate-bearing cemented sand with 24.2%, 35.1%, and 53.1% hydrate saturation, the deviatoric stress and volumetric strain simulation results are compared in Figure 9 and Figure 10. The pattern of strain softening with increasing hydrate-bearing cemented sand stiffness and maximum strength with increasing hydrate saturation can be adequately captured by the model since it accounts for both cementation factors and damage factors. As the effective restricting pressure decreases, the dilatancy becomes more noticeable. The curve’s peak strength appears at a very similar axial strain of about 2.3%. The maximum strength of the curve occurs earlier at hydrate saturations of 35.1% and 53.1% than at hydrate saturations of 24.2% and 0 at an axial strain of 2.1%, and the curve changes more steeply at the inflection point. This suggests that when the hydrate content rises, the hydrate and the yield deformation of soil break down sooner and more quickly. The soil can withstand more damage than gas hydrates, and its strength can be strengthened by increasing saturation.

The general trend of deformation of substrate sand and hydrate-bearing cemented sand with 52% saturation for various initial porosities is shown in Figure 11 and Figure 12. In Figure 8, as the hydrate saturation increases, the deviation of the theoretical simulation curves from the experimental curves of the volumetric strain ${\epsilon }_v$- axial strain ${\epsilon }_1$ relationship becomes more pronounced. This behavior could be explained by the fact that strain softening weakens and the strength of hydrate-bearing cemented sand decreases dramatically with increasing hydrate saturation. This suggests that while dense specimens show slight volume expansion behavior and strain softening, loose specimens exhibit compression behavior and strain hardening. There is a significant difference between the experimental data for a restricting pressure of 3 MPa and the predicted deviatoric stress q and volumetric strain ε_v, in contrast to the curves for an effective confinement pressure of 5 MPa shown by Figure 9 and Figure 10. The model introduces cementation factor, and the stress-strain relationship behaves more clearly with increasing hydrate saturation relative to the conversion of the experimental data from strain hardening to softening. It is demonstrated that the model for hydrate-bearing cemented sand with cementation damage, which considers the state dependency, is valid and can accurately simulate both the mechanical characteristics of cemented sand damage and the weakening of cementation.

5. Conclusions

(1)

A nonlinear state-dependent elasto-plastic constitutive model based on cementation and damage in hydrate-bearing cemented sand is constructed for the phenomena of natural gas hydrate cementation weakening and soil damage. This is achieved by introducing a quantitative description of the hydrate’s cementation and decomposition mechanism, taking into account the damage factor D_s and the cementation damage strength factor P_h. The model captures the deviatoric stress q, volumetric strain ε_v, and axial stress of hydrate-bearing cemented sands with acceptable accuracy at different hydrate saturation, effective confining pressure, and void ratio levels, when compared to data from drainage triaxial tests by Hyodo et al. [59] and other scholars, suggesting that the model is valid in representing the intricate mechanical behavior of compromised cementation and soil particles degradation in cemented sand.

(2)

The state parameters represent the potential strength of cemented soil, and further introduce the state-dependent dilatancy equations. Compared to the previous model, a set of model parameters can be used to reproduce the mechanical properties of hydrate-bearing cemented sand with varying effective perimeter pressures, hydrate saturations, and void ratios. The strength of hydrate-bearing cemented sand against deviatoric stress increased with increasing hydrate saturation and surrounding pressure, and the peak strength was higher than that of substrate sand. The dilatancy phenomenon becomes more pronounced with decreasing effective confining pressure. Additionally, this agrees with experimental results of Hyodo et al. [59]. For different porosities, the loose specimens show strain hardening and compression behavior, and the dense specimens show slight strain softening and volume expansion behavior, with minor deviations from the model calculations at lower peripheral pressures.

(3)

The simulated curves have slightly lower peak intensities compared to the experimental curves due to the inclusion of damage factors in the model that consider soil deterioration and gas hydrate decomposition. As the saturation of the gas hydrate increases, both the soil yield and the gas hydrate become more susceptible to distortion and breakdown. Unlike previous hydrate models, the damage index s controls the magnitude of the peak intensity of the stress-strain curve by measuring the rate of damage to the sand. The gas hydrate decomposition rate can be modified based on the impact of different external elements, which indicates the mechanical damage characteristics of cemented sand containing cemented hydrate under intricate circumstances.

(4)

This study exclusively focuses on the cementing impact of hydrate-bearing cemented sand, disregarding other mechanical consequences of hydrate, such as the reduction in sediment porosity. Consequently, the model utilized in this study requires improvement. The next step will consist of examining the legal consequences of hydrate cementation damage and developing a constitutive model that incorporates hydrate compaction.