A State-Dependent Constitutive Model for Gas Hydrate-Bearing Sediments Considering Cementing E ﬀ ect

: This paper presents a state-dependent constitutive model for gas hydrate-bearing sediments (GHBS), considering the cementing e ﬀ ect for simulating the stress–strain behavior of GHBS. In this work, to consider the inﬂuence of hydrate on matrix samples in theory, some representative GHBS laboratory tests were analyzed, and it was found that GHBS has obvious state-related characteristics. At the same time, it was found that GHBS has high bonding strength. In order to describe these characteristics of GHBS, the cementation strength related to hydrate saturation is introduced in the framework of a sand state correlation model. In addition, in order to accurately reﬂect the inﬂuence of cementation on the hardening law of GHBS, the degradation rate of cementation strength is introduced, and the mixed hardening theory is adopted to establish the constitutive model. The model presented in this paper reproduces the experimental results of Masui et al. and Miyazaki et al., and the prediction performance of the model is satisfactory, which proves the rationality of this work.


Introduction
With the continuous development and progress of exploration technology and equipment, the exploitation of oil and gas resources has gradually shifted from land and shallow sea to deep sea [1]. There are abundant resources in the ocean, especially natural gas hydrate, which can reach twice the total amount of traditional fossil fuels [2]. Gas hydrates are considered as one of the most promising strategic resources for commercial development in the future [3].
Gas hydrates are caged crystalline compounds formed by water and natural gas under certain temperature and pressure conditions. They are widely distributed in seabed strata and have the advantages of being clean, with large reserves and high heat value [3]. Some countries have carried out on-site exploitation of hydrate, and Mallik in Canada was the first place to carry out gas hydrate test mining [4]. Gas hydrate test mining was carried out three times, in 1998, 2002, and 2007-2008, respectively, and the production response characteristics of gas hydrate-bearing sediments (GHBS) were directly identified for the first time. However, due to the limitations of the production test cycle, this study is not enough to fully evaluate the resource potential of gas hydrate. In 2012, American researchers conducted the world's first hydrate production test project combining the carbon dioxide replacement method and depressurization method on the northern slope of Alaska, and it confirmed the feasibility of developing hydrate by using the replacement method on the technical level [5]. Since 2000, Japan has carried out the world's first sea area hydrate test production in Nankai on the stress-strain behavior of GHBS. Masui et al. [22] prepared two types of GHBS using percolated methane gas in a specimen of Toyoura sand mixed with powder ice and excessive water and later carried out a series of triaxial tests with different effective confining pressures. Yun et al. [23], Miyazaki et al. [24], Ghiassian et al. [25], Liu et al. [26], Madhusudhan et al. [27] and Dong et al. [28] conducted experimental studies on GHBS.
On the basis of existing research results, Soga et al. [29] have described the mechanical properties of GHBS in detail, but because the mechanical properties of GHBS are the basis of this study, in order to facilitate understanding, this paper briefly introduces the main mechanical properties of hydrate sediments: (1) GHBS may show strain hardening or softening during shearing, and its characteristics depend on S h , compactness and stress level. (2) For the deformation characteristics of GHBS, shear shrinkage and dilatancy may occur, which are mainly related to internal states, such as hydrate occurrence mode and S h . (3) The cementation of GHBS has an important influence on its characteristics of strength and deformation. With the decline of cementation, it often shows a more obvious phenomenon of softening and dilatancy.
To sum up, the strength characteristics of GHBS not only depend on the void ratio and stress level, making it similar to sand, but are also related to S h and cementation, thus producing essential differences in the mechanical properties of the ordinary soil and sand. Therefore, GHBS is also a kind of state-dependent material. By constructing appropriate expressions of the state parameter, the state-dependent characteristics of GHBS can be described.

Yield Function
Since GHBS exhibits characteristics similar to sand, the yield function of the sand model is modified. A large number of studies on sand have shown that the deformation caused by consolidation is very small, while the deformation of sand is mainly caused by stress ratio [30][31][32]. Poorooshasb et al. [33] assumed that the yield trajectory of sand is a series of straight lines with a constant stress ratio. In Figure 1, the schematic diagram of the yield surface shows the shape of the yield surface, which is represented by Equation (1).
where q*, and p* are the mean normal stress and generalized shear stress of sand, respectively, and M* is the stress ratio of critical state.  [27] and Dong et al. [28] conducted experimental studies on GHBS.
On the basis of existing research results, Soga et al. [29] have described the mechanical properties of GHBS in detail, but because the mechanical properties of GHBS are the basis of this study, in order to facilitate understanding, this paper briefly introduces the main mechanical properties of hydrate sediments: (1) GHBS may show strain hardening or softening during shearing, and its characteristics depend on Sh, compactness and stress level. To sum up, the strength characteristics of GHBS not only depend on the void ratio and stress level, making it similar to sand, but are also related to Sh and cementation, thus producing essential differences in the mechanical properties of the ordinary soil and sand. Therefore, GHBS is also a kind of state-dependent material. By constructing appropriate expressions of the state parameter, the statedependent characteristics of GHBS can be described.

Yield Function
Since GHBS exhibits characteristics similar to sand, the yield function of the sand model is modified. A large number of studies on sand have shown that the deformation caused by consolidation is very small, while the deformation of sand is mainly caused by stress ratio [30][31][32]. Poorooshasb et al. [33] assumed that the yield trajectory of sand is a series of straight lines with a constant stress ratio. In Figure 1, the schematic diagram of the yield surface shows the shape of the yield surface, which is represented by Equation (1).
where q*, and p* are the mean normal stress and generalized shear stress of sand, respectively, and M* is the stress ratio of critical state. Considering that the cementation of GHBS has an important influence on the yield strength, the straight line represented by Equation (1) is translated to the left to establish the yield function of GHBS, as shown by Equation (2), whereas the shape of the yield surface is shown in Figure 1. If the attenuation of cementation is considered, the subsequent yield surface of GHBS will simultaneously move to the right with the loading process. Considering that the cementation of GHBS has an important influence on the yield strength, the straight line represented by Equation (1) is translated to the left to establish the yield function of GHBS, as shown by Equation (2), whereas the shape of the yield surface is shown in Figure 1. If the attenuation of cementation is considered, the subsequent yield surface of GHBS will simultaneously move to the right with the loading process.
where q and p are the mean normal stress and generalized shear stress of the GHBS, respectively, p c is the cementation of GHBS, and M is the critical state stress ratio, which is equal to M*. When the stress ratio reaches M, GHBS enters the state of destruction. The stress ratio η is given by Equation (3).

Cementation
Referring to Sanchez's and Uchida's method for determining the cementation of GHBS [14,34], it is assumed that the initial value of cementation p c0 satisfies the exponential relationship with S h (see Equation (4)).
where α and β are the model parameters, while all other terms have the same meaning as defined previously. Compared with ordinary sand, shear behavior not only causes rotation, slip and redistribution of particles in the material, but it also causes hydrate breakage and debonding between the hydrate and sand particles [35]. With the development of strain, cementation will decline. This difference directly affects the macroscopic mechanical behavior of GHBS. In research on the evolutionary pattern of cementation of GHBS, several authors [14,18,34] assume that it is a function of plastic shear strain (ε p s ) and have obtained ideal prediction effect. In this paper, the same assumption is made, and it is found that ε p s cannot be recovered. The evolutionary pattern of cementation of GHBS is expressed using Equation (5). dp c = −a c p c0 dε p s where a c indicates the attenuation rate of the cementation of GHBS.

State Parameters and Dilatancy
The dilatancy of GHBS depends on the internal state and stress level, which can be reflected by the yield function. Therefore, it is necessary to construct a state parameter expression that reflects the internal state and then to clarify the dilatancy of GHBS by describing the current state.
If e is the current void ratio of the material, and e c is the critical void ratio, then the state parameter ϕ can be expressed as the difference [36] and the ratio [37] between e and e c . Yao et al. [32] use the void ratio e η corresponding to the mean normal stress p under the condition of equal stress ratio to replace e c and construct the state parameter ϕ considering the stress ratio. In this paper, the state parameter ϕ is expressed using Equation (6).
According to the basic theory of soil mechanics, the relationship between the current void ratio e of GHBS, the initial void ratio e 0 and the volume strain ε v is given by Equation (7).
In order to obtain the expression of the state parameter ϕ of the GHBS, it is necessary to establish the expression of the initial void ratio e 0 and the critical void ratio e c . The hydrate has the effects of filling and cementation on GHBS. Furthermore, for the GHBS mainly filled with the filling effect, e 0 may decrease with the increase of S h ; however, for the GHBS-containing structure, the larger the S h , the smaller is the e 0 . The influence of the content and occurrence mode of hydrate on e 0 is not completely clear at present. Additionally, the samples of GHBS prepared indoors are affected by methane gas scouring, which may lead to uneven samples. Therefore, it is difficult to establish the mathematical relationship between e 0 and S h . In order to reduce the model parameters, this paper assumes that e 0 and S h of GHBS satisfy a simple inverse proportional relationship, as given by Equation (8).
where e * 0 is the initial void ratio of the host sand, and k is the model parameter that reflects the influence of hydrate occurrence mode and formation conditions on the initial void ratio of GHBS.
Yoneda et al. [38] verified that the deformation process affects the occurrence conditions of hydrate and leads to hydrate decomposition, whereas the consolidation line of GHBS gradually approaches the consolidation line of the host sand, due to which the GHBS should have the same critical state line (CSL) as the host sand. The CSL of sand can be expressed using Equation (9).
where e c0 is the void ratio of CSL when p = 1kPa, λ is the slope of CSL in e-lnp space, ξ is the material parameter with the generally accepted value of 0.7 for sand, and p a is the atmospheric pressure (101kPa) introduced for dimensional unification. Using the expression of dilatancy equation proposed by Wan and Guo [37], the state parameter ϕ, as given by Equation (6), is introduced into the dilatancy equation, which is given by Equation (10).
where m is the model parameter, which can be obtained from the results of the triaxial test. It can be seen that when m = 0, Equation (10) is reduced to the dilatancy equation of the Cambridge model.

Hardening Modulus
The plastic strain of GHBS can be divided into two parts. One part is related to the expansion of the yield surface. For this part, η is used to reflect the isotropic hardening of GHBS with the evolution of plastic strain. The other part results from the movement of the yield surface. Cementation reflects the motion hardening of GHBS with the evolution of plastic shear strain. Both of them constitute the mixed hardening of GHBS.
The hardening caused by the stress ratio is the same as the hardening modulus expression of sand, while the isotropic hardening modulus A 1 is given by Equation (11).
where h is the model parameter and G is the shear modulus of GHBS. For sand, the shear modulus has the following empirical correlation (Equation (12)).
where G* is the shear modulus of sand and G * 0 is the model parameter. Yoneda et al. [38] have given the empirical correlation of the shear modulus of GHBS under 1MPa effective confining pressure, which is given by Equation (13).
Based on Equations (12) and (13), the shear modulus of GHBS is rewritten and given by Equation (14).
where χ is the material parameter of GHBS and reflects the degree of influence of the hydrate on the hardening modulus. The cementation of GHBS also causes material hardening during shearing, which has the characteristics of kinematic hardening. In this paper, the kinematic hardening modulus is represented by Equation (15).
Therefore, the hardening modulus of GHBS can be expressed using Equation (16).
It can be seen that, when S h = 0, and p c0 = 0, Equation (16) returns to the hardening modulus expression of sand.

Model Parameters
The model established in this paper contains thirteen parameters and can be divided into two categories, as represented by the parameters given in Tables 1 and 2. Among them, sand material parameters ξ, λ, e c0 and µ can be obtained through the isotropic compression test of the host sand, while M, h, m and G * 0 are calibrated according to the results of the triaxial test. The hydrate-related parameters k, α, β and a c are obtained by fitting the triaxial test results of different S h values, whereas χ is calculated using Equation (14).

Model Prediction
Masui et al. [22] took Toyoura sand as the matrix material with an initial void ratio of 0.65, prepared GHBS samples with different S h and conducted triaxial drained tests under the effective confining pressure of 1MPa. The test data of Masui are predicted using the model proposed in the current paper, and the parameters used are presented in Tables 1 and 2. Figure 2 shows the comparison between experimental and model results for different S h values in terms of deviatoric stress. It can be seen that the predicted results of the model are consistent with the test data, and the model can satisfactorily reflect the strength and strain softening characteristics of GHBS. With the increase in S h , the stiffness and strength of GHBS gradually increase, and the phenomenon of strain-softening becomes more obvious. The prediction curve of the volume strain of GHBS is shown in Figure 3. Compared to the experimental data, it can be seen that the model can satisfactorily describe the effect of hydrates on dilatancy. When the strain is small, GHBS mainly shows shear shrinkage. With the increase in axial strain, the volume strain gradually increases. Furthermore, higher the value of S h , more obvious is the dilatancy.  In order to test the applicability of the model under different effective confining pressures, it is compared with the test data of Miyazaki [24]. Miyazaki et al. also used Toyoura sand and methane to prepare samples with approximately equal Sh values. The initial void ratio of host sand samples is 0.6, while the effective confining pressures are varied consecutively through values of 0.5MPa, 1MPa, 2MPa and 3MPa. During the tests, the axial strain εa and radial strain εr of GHBS are monitored.
As the test conditions of Miyazaki et al. [24] and Masui et al. [22] are not the same, the model parameters are changed. Table 3 presents the parameters used to predict the test data of Miyazaki et al. The results predicted by the model are compared with the experimental data. It can be seen from Figure 4 that, for the same Sh, the increase in effective confining pressure makes the GHBS change from strain-softening to strain-hardening, because the constraint effect of effective confining pressure improves the rigidity and strength of GHBS. Meanwhile, the rate of decline of cementation of GHBS decreases. Figure 5 shows the comparison between the relationship curve of predicted radial strain and the experimental data. The predicted curve is in good agreement with the experimental data, indicating that the model proposed in the current paper has a very good applicability. . Figure 2. Comparison of the predicted stress-strain curves with the experimental results [22]. In order to test the applicability of the model under different effective confining pressures, it is compared with the test data of Miyazaki [24]. Miyazaki et al. also used Toyoura sand and methane to prepare samples with approximately equal Sh values. The initial void ratio of host sand samples is 0.6, while the effective confining pressures are varied consecutively through values of 0.5MPa, 1MPa, 2MPa and 3MPa. During the tests, the axial strain εa and radial strain εr of GHBS are monitored.
As the test conditions of Miyazaki et al. [24] and Masui et al. [22] are not the same, the model parameters are changed. Table 3 presents the parameters used to predict the test data of Miyazaki et al. The results predicted by the model are compared with the experimental data. It can be seen from Figure 4 that, for the same Sh, the increase in effective confining pressure makes the GHBS change from strain-softening to strain-hardening, because the constraint effect of effective confining pressure improves the rigidity and strength of GHBS. Meanwhile, the rate of decline of cementation of GHBS decreases. Figure 5 shows the comparison between the relationship curve of predicted radial strain and the experimental data. The predicted curve is in good agreement with the experimental data, indicating that the model proposed in the current paper has a very good applicability.  In order to test the applicability of the model under different effective confining pressures, it is compared with the test data of Miyazaki [24]. Miyazaki et al. also used Toyoura sand and methane to prepare samples with approximately equal S h values. The initial void ratio of host sand samples is 0.6, while the effective confining pressures are varied consecutively through values of 0.5MPa, 1MPa, 2MPa and 3MPa. During the tests, the axial strain ε a and radial strain ε r of GHBS are monitored.
As the test conditions of Miyazaki et al. [24] and Masui et al. [22] are not the same, the model parameters are changed. Table 3 presents the parameters used to predict the test data of Miyazaki et al. The results predicted by the model are compared with the experimental data. It can be seen from Figure 4 that, for the same S h , the increase in effective confining pressure makes the GHBS change from strain-softening to strain-hardening, because the constraint effect of effective confining pressure improves the rigidity and strength of GHBS. Meanwhile, the rate of decline of cementation of GHBS decreases. Figure 5 shows the comparison between the relationship curve of predicted radial strain and the experimental data. The predicted curve is in good agreement with the experimental data, indicating that the model proposed in the current paper has a very good applicability.

Conclusions
In this study, an elastoplastic model is proposed based on the framework of the sand statedependent constitutive model, in which GHBS is considered as a bonded material. The main conclusions are as follows: (1) By summarizing the test rules of GHBS, it is found that GHBS has certain cohesive strength and obvious state-related characteristics. In this paper, GHBS is regarded as a special cementing material, and the sand state correlation model is selected as the basic framework.

Conflicts of Interests:
The authors declare no conflict of interest.

Conclusions
In this study, an elastoplastic model is proposed based on the framework of the sand statedependent constitutive model, in which GHBS is considered as a bonded material. The main conclusions are as follows: (1) By summarizing the test rules of GHBS, it is found that GHBS has certain cohesive strength and obvious state-related characteristics. In this paper, GHBS is regarded as a special cementing material, and the sand state correlation model is selected as the basic framework. (2) By constructing appropriate state parameter expressions, introducing the cementation and its evolutionary pattern and adopting the mixed hardening rule, the state-dependent constitutive model of GHBS is established.

Conflicts of Interests:
The authors declare no conflict of interest. Figure 5. Comparison of the predicted stress-strain curves with the experimental results [24].

Conclusions
In this study, an elastoplastic model is proposed based on the framework of the sand state-dependent constitutive model, in which GHBS is considered as a bonded material. The main conclusions are as follows: (1) By summarizing the test rules of GHBS, it is found that GHBS has certain cohesive strength and obvious state-related characteristics. In this paper, GHBS is regarded as a special cementing material, and the sand state correlation model is selected as the basic framework. (2) By constructing appropriate state parameter expressions, introducing the cementation and its evolutionary pattern and adopting the mixed hardening rule, the state-dependent constitutive model of GHBS is established.

Conflicts of Interest:
The authors declare no conflict of interest.