Geomechanical Properties of Deep-Sea Pore-Filled Methane Hydrate-Bearing Soils at Critical State Using DEM Analysis

Abstract

The recognition of the geomechanical properties of methane hydrate-bearing soil (MHBS) is crucial to exploring energy resources. The paper presents the mechanical properties of a pore-filled MHBS at a critical state using the distinct element method (DEM). The pore-filled MHBS was simulated as cemented MH agglomerates to fill the soil pores at varying levels of methane hydration (MH) saturation. A group of triaxial compression (TC) tests were conducted, subjecting MHBS samples to varying effective confining pressures (ECPs). The mechanical behaviors of a pore-filled MHBS were analyzed, as it experienced significant strains leading to a critical state. The findings reveal that the proposed DEM successfully captures the qualitative geomechanical properties of MHBS. As MH saturation increases, the shear strength of MHBS generally rises. Moreover, higher ECPs result in increased shear strength and volumetric contraction. The peak shear strength of MHBS increases with rising MH saturation, while the residual deviator stress remains mainly unchanged at a critical state. There is a good correlation between fabric changes of the MHBS with variations in principal stresses and principal strains. With increasing axial strain, the coordination number (CN) and mechanical coordination number (MCN) increase to peak values as the values of MH saturation and ECPs increase, and reach a stable value at a larger axial strain.

1. Introduction

Methane hydrate-bearing soil (MHBS) deposits are found abundantly in continental margins and permafrost regions, which contain MH within the pore spaces. MH is a promising resource to alleviate the current energy shortage ([1,2,3,4,5]). Therefore, there is a growing focus on surveying the physical, thermodynamic, chemical, and mechanical properties of MHBS, as well as developing technological methods for identifying and exploiting MHBS deposits ([6,7]). However, it is important to consider the potential hazards associated with safe mining, such as submarine landslides and tsunamis ([8,9,10,11,12,13,14,15]). Therefore, the geomechanical properties of MHBS are crucial to understanding from the view of geotechnique ([16,17,18,19]).

As a type of geotechnical material, when it comes to MHBS, our understanding of its mechanical properties mostly relies on conventional triaxial compression tests ([20,21,22]). Undeniably, triaxial tests are widely used and effective for investigating the macroscopic mechanical characteristics of geomaterials. However, as our understanding of geotechnical materials deepens, it has become evident that the microstructure within the material significantly influences its macroscopic mechanical properties ([23,24,25,26]). For certain specific macroscopic phenomena, such as shear bands and anisotropy in soil, it is necessary to investigate the microstructural aspects of specimens to obtain a more profound and reasonable explanation.

Nevertheless, due to limitations in experimental equipment, it is often challenging to directly observe mechanical responses at the particle level within the soil during conventional triaxial tests. Compared to traditional geotechnical materials, MHBS specimens possess more complex material compositions, typically consisting of soil particles, water, gas, and methane hydrates. These factors exert multiple influences on the mechanical properties, particularly the content, distribution of hydrates, and hydrate morphology within MHBS ([27,28,29,30,31,32]). Therefore, relying solely on conventional compression tests may not sufficiently illuminate the geomechanical behaviors of MHBS.

The difficulty in obtaining undisturbed MHBS samples and the requirement of specialized low-temperature and high-pressure testing apparatus in indoor triaxial tests significantly hinder the systematic investigation of MHBS mechanical properties. Additionally, their mechanical characteristics are relevant to the attributes and statuses of their constituent components, the interactions between solid hydrates and soil particles playing a prominent role in mechanical behavior. Thus, understanding the mechanical properties of MHBS necessitates considering the effects of different hydrate forms on the soil framework and the resulting diverse mechanical behaviors.

Consequently, the distinct element method (DEM) has garnered attention from researchers ([29,33,34,35,36,37,38,39,40]). The DEM allows for generating MHBS specimens with different microstructural distributions from a microscopic perspective, enabling numerical simulations to explore the macro-/micro-mechanical responses of particle assemblies under varying loading situations, thereby strengthening the understanding of the mechanical performances of MHBS. The DEM is instrumental in enhancing our comprehension of the geomechanical behavior of MHBS with different MH saturations in the critical state, thereby facilitating the establishment of a suitable constitutive model.

The sand-based MHBS specimens with a pore-filling type of MH, which are composed of a mixture of sand particles and hydrate particles, represent a special granular material with distinct discontinuities. A reproducible numerical specimen preparation method is firstly introduced. A group of drained triaxial compression (TC) tests were conducted to investigate the mechanical and deformation responses of a pore-filled MHBS, subjecting the MHBS samples to varying effective confining pressures (ECPs). The mechanical and deformation properties of the pore-filled MHBS were analyzed when the specimen was experienced to a significant strain, leading to a critical state. The qualitative properties of the pore-filling type of MHBS are then evaluated in different respects.

2. Materials and Methods

A cube with a side length of 5.44 mm was generated, which contained 10,000 soil particles, using a particle size distribution (PSD) illustrated in Figure 1a.

A novel numerical fabrication method for preparing the pore-filled MHBS was introduced in a previous study ([18]). The main idea of the approach to generate a pore-filled MHBS is as follows: Firstly, the MH cluster is regarded as a whole cluster entity formed by the agglomeration of granular bonding. This implies that the hydrate deposits are modeled as assemblies consisting of numerous small particles that are interconnected through bonding interactions. Next, these hydrate clusters are randomly filled within the pores of soil particles, creating MHBS samples with the corresponding MH saturation. This technique aims to simulate and understand the formation process of porous filling hydrate deposits. Here is just a brief description of this method. Detailed information can be referred to in the previous reference ([18]).

MH saturation is calculated by the following equation:

$$S_{MH}=\frac{V_{MH}}{V_V}\times 100\%$$

(1)

where $V_{MH}$ is the MH’s volume of MH and $V_V$ is the total pore volume containing MH.

Four MHBS samples were generated with MH saturations ($S_{mh}$ = 0%, 5%, 10%, and 20%), in which the numbers of MH clusters (N = 0, 540, 1080, and 2160) were filled in the soil skeleton. One MH cluster with a void ratio $e_a$ = 0.82 and a radius $r$ = 0.024 m was produced which contained 60 MH particles cemented by the parallel bond model (PBM) ([41,42]). The parameters of the PBM used were as follows: the normal bond strength was ${\overline{\sigma }}_{\mathrm{c}}=11.37 \mathrm{MPa}$, and the shear parallel bond strength was ${\overline{\tau }}_{\mathrm{c}}=11.37 \mathrm{MPa}$. The normal parallel bond stiffness was to be ${\overline{k}}^{\mathrm{n}}$ = 9.244 × 10¹² $\mathrm{N}\cdot {\mathrm{m}}^{-3}$. The shear parallel bond stiffness was to be ${\overline{k}}^{\mathrm{s}}$ = 6.163 × 10¹² $\mathrm{N}\cdot {\mathrm{m}}^{-3}$. The ratio of parallel bond radius to ball radius was $\overline{\lambda }=1$. The schematic diagram after sample preparation with a specified MH saturation is shown in Figure 1b, in which MH particles are highlighted in red. Three ECPs (${\sigma }_3$ = 2, 4, and 6 MPa) were chosen and were constant with the stress boundary to conduct the simulation of drained triaxial compression tests. The major principal stress was loaded at a compression strain rate of 1%/min.

The deviatoric shear stress, $q$, can be obtained as follows:

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

(2)

The volumetric strain, ${\epsilon }_v$, is calculated as follows ([43]):

$${\epsilon }_v=\frac{△V}{V}=(1+{\epsilon }_1)(1+{\epsilon }_2)(1+{\epsilon }_3)-1$$

(3)

where ${\epsilon }_1$, ${\epsilon }_2$, and ${\epsilon }_3$ are the major, intermediate, and minor principal strain, respectively.

The void ratio is calculated as follows:

$$e=\frac{V_v}{V_s}=\frac{V-V_s}{V_s}=\frac{V}{V_s}-1$$

(4)

where $V_v$ is the volume of the voids, $V_s$ is the volume of the soil, and V is the total volume of the MHBS.

The parameters used were listed in

Table 1. Parameters of granules used in the DEM.

Parameters	Value	Parameters	Value
Soil density/(g/cm3)	2.65	MH density/(g/cm3)	0.9
Wall normal stiffness/(N/m)	1.5 × 105	Wall tangential stiffness/(N/m)	1.0 × 105
$\mathrm{Soil} \mathrm{normal} \mathrm{contact} \mathrm{stiffness}, k_n^p$/(N/m)	1.5 × 105	$\mathrm{Soil} \tan \mathrm{gential} \mathrm{contact} \mathrm{stiffness}, k_s^p$/(N/m)	1.0 × 105
$\mathrm{MH} \mathrm{normal} \mathrm{contact} \mathrm{stiffness}, k_n^p$/(N/m)	1.5 × 104	$\mathrm{MH} \tan \mathrm{gential} \mathrm{contact} \mathrm{stiffness}, k_s^p$/(N/m)	1.0 × 104
$\mathrm{Soil} \mathrm{interparticle} \mathrm{coefficient} \mathrm{of} \mathrm{friction}, {\mu }^p$	0.5	$\mathrm{MH} \mathrm{interparticle} \mathrm{coefficient} \mathrm{of} \mathrm{friction}, {\mu }^p$	0.04

as follows.

3. Results

The contact model ([32]) was chosen to simulate the mixed type of MHBS, in which part of the MHBS is in the bonding state. If readers need detailed information, the reference ([32]) can be referred to. Here is just a brief introduction to the basic mechanical properties of the model. As shown in Figure 1, the mechanical behaviors include tension, compression, shearing, bending, and twisting modes of loading. The solid lines represent the mechanical responses with cementation, while the dashed lines represent mechanical responses without cementation. These actions lead to two failure mechanisms: one involving the breakage and fracture of the hydrate itself, and the other involving the cementation failure of the hydrate particles with the skeleton particles.

3.1. Macro-Mechanical Properties of the Pore-Filling Type of MHBS

3.1.1. Stress–Strain Response

Four MH saturations ($S_{mh}$ = 0%, 5%, 10%, and 20%) and three ECPs (${\sigma }_3$ = 2, 4, and 6 MPa) were selected to study the macro-/microscopic mechanical properties of the pore-filling type of MHBS. Figure 2a–d presents the stress–strain relationship (${\sigma }_1-{\sigma }_3~{\epsilon }_1$) of the pore-filled MHBS under different MH saturations and effective confining pressures. At a given MH saturation level, the deviatoric stress of the MHBS increases as the ECP increases. Increasing the ECP enhances the contact forces between soil particles, resulting in a denser soil structure and increased shear strength. As the ECP increases, the MHBS sample reaches its peak strength, and with the increasing axial strain, the deviatoric strain gradually stabilizes, eventually reaching a critical state.

Figure 3a–c presents the deviatoric stress–axial strain curves of the MHBS obtained from TC tests with different MH saturations at the same ECP. It can be observed that the variation in MH saturation had a minimal influence on the initial stiffness of the MHBS sample. The peak deviatoric stress strengths of the pore-filled MHBS increased with an increase in MH saturation.

Figure 4a–c shows the relationship between stress ratio and axial strain at different MH saturations ($S_{mh}$ = 0%, 5%, 10%, and 20%) under different ECPs (${\sigma }_3$ =2.0 MPa, 4.0 MPa, and 6.0 MPa). Based on the simulated observations, the peak stress ratio of the MHBS increases gradually with increasing MH saturation. The formation of hydrates in the soil enhances its shear strength. The presence of hydrates reduces the contact area between soil particles, thereby increasing the cohesion of the soil and consequently its strength.

Figure 5 presents the peak strength envelopes on the plane for the pore-filled MHBS samples at different MH saturations. With an increase in hydrate saturation, the slope of the peak strength envelope of the MHBS samples tends to increase. Based on the stress–strain relationships obtained from TC tests on the MHBS samples, the internal friction angle and cohesion of the MHBS samples at different MH saturations could be determined, as shown in Figure 6. The internal friction angle of the pore-filled MHBS samples approximately increased linearly with an increase in MH saturation. However, the trend of cohesion variation for the pore-filling type of MHBS samples with increasing MH saturation is not clear. This is because the hydrate clusters act as a whole to fill the pores of the soil skeleton, where the soil skeleton primarily carries the load. This is consistent with the characteristics of the pore-filled MHBS observed by the reference ([28]), indicating that the contribution of hydrates to the mechanical properties of the pore-filled MHBS samples was mainly in terms of frictional characteristics rather than cementation characteristics.

3.1.2. Volumetric–Axial Strain Response

Figure 7a–e presents the volumetric strain–axial strain curves of MHBS samples with different MH saturations ($S_{mh}$ = 0%, 5%, 10%, and 20%) under different ECPs (${\sigma }_3$ = 2, 4, and 6 MPa). It can be observed that in Figure 7a–d, at the same MH saturation, the volumetric shrinkage of the MHBS samples becomes more pronounced with increasing ECP. This volume change behavior is qualitatively consistent with previous results ([28]). The sediment endures initially volumetric contraction then volumetric dilatancy. The volumetric strain shows less dilation with increasing ECP at a given axial strain. A higher ECP will strengthen the breakage of MH bonds, then the particles move into the pore space during deformation, decreasing the MHBS dilation.

As shown in Figure 7e, volumetric strain shows marked dilation after an early shear contraction. Under the same axial strain, the phenomenon of shear dilation becomes more pronounced with increasing MH saturation. The presence of MH thus enhances dilation. As the MH saturation increases, this results in an increase in the contact area between soil particles and an increase in the interparticle contact forces. Additionally, the presence of hydrate clusters provides resistance to shear deformation. Therefore, as the hydrate saturation increases, the soil becomes more resistant to deformation under shear stress, leading to a more pronounced shear dilation phenomenon.

3.1.3. Critical State

Figure 8 shows the relations between the deviatoric stress and the mean principal stress of MHBS. It can be seen that as the MH saturation increases, the slope of the critical state line of the MHBS specimen also increases. By fitting the data from Figure 8, the variation of the critical-state stress ratio of the MHBS specimen with MH saturation can be plotted, as illustrated in Figure 9. It can be observed that as the MH saturations increase, the critical-state stress ratio of the MHBS also increases.

The critical-state shear stress ratio is expressed in a nonlinear relationship with MH saturation, as shown in the following equation:

$$M(S_{mh})=M+{\alpha }_1{S}_{mh}^{{\alpha }_2}$$

(5)

where $M(S_{mh})$ is the critical stress ratio of the MHBS at a given MH saturation, $M$ is the critical-state stress ratio of the MHBS sample indicating zero hydrate saturation, ${\alpha }_1$ and ${\alpha }_2$ are the model parameters, and ${\alpha }_1{S}_{mh}^{{\alpha }_2}$ is the parameter that considers the presence of MH in the MHBS and the influence of the mechanical properties of the critical state of the soil sample. $M=0.7571$, ${\alpha }_1=0.0022$, and ${\alpha }_2=1.0$ can be obtained based on Equation (5).

Figure 10 shows the critical state lines of the MHBS specimens on a plane of $e-lnp$ at different MH saturations ($S_{mh}$ = 0%, 5%, 10%, and 20%). It can be observed that the critical state lines of the MHBS specimens slightly shift upward with increasing MH saturation. The following equation was used for fitting:

$$e_{cr}(S_{mh})=e_{cr0}(S_{mh})-{\lambda }_{c}ln\left(\frac{p}{p_a}\right)$$

(6)

where $p_a$ is the standard atmospheric pressure, 1.01 × 10⁵ Pa, and ${\lambda }_c=0.087$.

When the average principal stress is taken at the standard atmospheric pressure, the reference pore ratio at the critical state under different MH saturations can be obtained, as shown in Figure 11. The relations between the critical-state void ratio and MH saturation are fitted using the following formula:

$$e_{cr0}(S_{mh})=e_{cr00}+e_1{S}_{mh}^{e_2}$$

(7)

where $e_{cr00}$, $e_1$, and $e_2$ represent the model parameters, taking $e_{cr00}=0.9598$, $e_1=0.0003$, and $e_2=1.0$, respectively. $e_1{S}_{mh}^{e_2}$ considers the influence of MH presence on the critical-state mechanical properties of the MHBS specimens.

3.1.4. Dilatancy

The dilatancy index is stated as follows ([44,45]):

$$d=d{\epsilon }_v/d{\epsilon }_s$$

(8)

where $d{\epsilon }_v$ is the change in volumetric strain and $d{\epsilon }_s$ is the increment in the equivalent deviatoric strain.

Figure 12a–d represents the stress ratio and dilatancy of the pore-filled MHBS samples with MH saturation under different ECPs. The stress ratio increases nonlinearly with increasing dilatancy, and then the growth rate gradually decreases. Under the same hydrate saturation, it increases with the increase in ECP. When the peak intensity is reached, the stress ratio is reduced and still maintains a certain growth. Subsequently, the stress ratio is reduced and ultimately fluctuates in the critical state.

3.2. Micro-Mechanical Properties of the Pore-Filled MHBS

3.2.1. Fabric Component

A unit normal contact, $F_{ij}$, is defined as the normal to the tangent plane of the contact. When the diagonal elements, i.e., $F_{11}$, $F_{22}$, and $F_{33}$, are equal, this indicates isotropy in the spatial distribution of contact directions. A higher value indicates a greater tendency for the contact directions to be distributed in that specific direction. The strong normal contact force, $F_{{}_{}}^s$, is stated as those forces equal to and larger than the average on every particle, with the respective contacts indexed as $n_{}^s$, isolating the response of the strong and weak force network. The unit normal vector for the strong contacts is calculated as follows ([40,45]):

$$F_{{}_{ij}}^s=\frac{1}{N_c}{\displaystyle \sum _{s=1}^{N_{{}_c}^s}{n}_{{}_i}^s{n}_{{}_j}^s}, i=j=1,d$$

(9)

where $N_c^s$ is the number of the strong contacts and $n_i^s$ is the component of the unit vector $n^s$ at a strong contact.

The mean strong contact is used to calculate using the following equation:

$$F_{mean}^s=\frac{F_{11}^s+F_{22}^s+F_{33}^s}{3}$$

(10)

Then, the anisotropic enhancement coefficient, ${\Delta }_c^s$, is used as follows:

$${\Delta }_c^s=\sqrt{\frac{3}{2}[{(F_{11}^s-F_{mean}^s)}^2+{(F_{22}^s-F_{mean}^s)}^2+{(F_{33}^s-F_{mean}^s)}^2]}$$

(11)

Figure 13 shows three principal values of fabric components ($F_{11}$, $F_{22}$, and $F_{33}$), strong fabric components ($F_{11}^s$, $F_{22}^s$ and $F_{33}^s$), and ${\Delta }_c^s$ of the pore-filled MHBS plotted against axial strain under ${\sigma }_3$ = 2.0 MPa and $S_{mh}=5\%$ with the growth of axial strain.

.

In Figure 13a, the major principal value, $F_{11}$, of the contact fabric tensor initially increases with the growth of axial strain. As the axial strain increases, $F_{11}$ slightly decreases and gradually approaches a stable value. On the other hand, the intermediate principal value, $F_{22}$, and the minor principal value, $F_{33}$, of the contact fabric tensor decrease initially with the increase in axial strain and then slightly increase towards a stable state. For the strong contact fabric tensor illustrated in Figure 13b, the trend of the major principal value with axial strain development is qualitatively similar to that of the regular contact fabric tensor. However, the variation of $F_{11}^s$ is more pronounced, and the stable value of the major principal value in the strong contact fabric tensor is significantly higher than that of the regular contact fabric tensor. The values of $F_{22}^s$ and $F_{33}^s$ decrease initially with the increase in axial strain and then increase slightly with further axial strain development, approaching a stable state. Figure 13c illustrates the variation of the anisotropic enhancement coefficient, ${\Delta }_c^s$, for the strong contact fabric with axial strain. The value of ${\Delta }_c^s$ initially increases to a peak value with the increase in axial strain and then decreases gradually with further axial strain, indicating a softening behavior in the anisotropic enhancement coefficient. The observed pattern of the anisotropic enhancement coefficient aligns with the increases in deviatoric stress and axial strain in the MHBS specimen, illustrated in Figure 3, demonstrating a strain-softening characteristic.

3.2.2. Contact Rose Diagram

Figure 14, Figure 15 and Figure 16 show the distribution of contact normal vectors and strong contact normal vectors projected on the x–y, x–z, and y–z planes with ${\sigma }_3$ = 2.0 MPa and $S_{mh}=5\%$ under different axial strains (${\epsilon }_1=2.669\%$, ${\epsilon }_1=6.647\%$, and ${\epsilon }_1=12.615\%$). In the x–y plane, the contact point fabric of the specimen is relatively uniformly distributed. However, in the x–z and y–z planes, the proportion of horizontally oriented strong contact fabric decreases, while the proportion of vertically oriented contact points significantly increases. The strong contact fabric exhibits a higher proportion of vertical contact points in the x–z and y–z planes, while the proportion of horizontal contact points decreases. The major principal value of the strong contact fabric tensor in the MHBS specimen is significantly greater than that of the regular contact fabric tensor.

3.2.3. Coordination and Mechanical Coordination Number

The coordination number (CN) of a particle represents the number of contacts it has with surrounding particles. The average coordination number is usually calculated as follows ([46]):

$$Z=\frac{2C}{N}$$

(12)

where C is the number of contacts and N is the number of particles.

During the shearing process of granular material specimens, there may exist particles that have only one contact or no contacts. These contacts do not contribute to the load transfer within the particle system and do not contribute to the resistance against external forces. Therefore, a mechanical coordination number (MCN) is defined as:

$$Z_m=\frac{2C-N_1}{N-N_0-N_1}$$

(13)

where $N_1$ and $N_0$ are the numbers of particles with only one or no contacts, respectively.

Figure 17a,b shows the CN and MCN of the soil skeleton in the pore-filling type of MHBS with $S_{mh}=5\%$ under different confining pressures. As the axial strain increases, the CN and MCN of the soil skeleton initially increase to a peak value, then decrease and gradually reach a stable state. The values of $Z$ and $Z_m$ increase as the value of ECPs increases. Figure 18a,b shows the CN and MCN of the soil skeleton in the pore-filling type of MHBS under different MH saturations. The CN and MCN of the soil skeleton initially increase to a peak value, then decrease and gradually lead to a stable state. The values of $Z$ and $Z_m$ increase with an increase in the value of MH saturations.

4. Discussion

Gas hydrates have different occurrence modes, i.e., pore-filling type, supporting matrix, cementation, coating, etc. Thus, different occurrence modes will have different influences on the mechanical and deformation properties of MHBS. This paper only studies one of the modes, i.e., pore-filled MHBS. Other modes and the existence of mixtures of different modes still need further research. This is a work of basic research, and its future application lies in the formulation of a constitutive theory of MHBS and large-scale boundary issues.

5. Conclusions

A series of TC tests was carried out using the DEM under different ECPs and MH saturations. The macro-/micro-mechanical properties of the pore-filling type of MHBS were analyzed by linking the evolution of micro-mechanical parameters to the macro-scale behavior of the granular material. The main conclusions are as follows:

(1) There is a good correlation between fabric changes of the MHBS with the variations in principal stresses and principal strains.

(2) The slope of the critical lines of the pore-filled MHBS increases as the MH saturation increases.

(3) The stress ratio increases nonlinearly with increasing dilatancy, and then the growth rate gradually decreases.

(4) The deviator stress causes anisotropy as the proportion of vertical contact points tends to increase and produce a vertical deflection.

(5) The CN and MCN of MHBS increase as the ECPs and MH saturations increase.