Discrete Element Modeling of the Effect of Hydrate Distribution Heterogeneity on the Mechanical Behavior of Cemented Hydrate-Bearing Sediments

Abstract

Hydrate distribution heterogeneity is often observed in natural and artificial hydrate-bearing sediments (HBSs). To capture hydrate distribution heterogeneity, a pore-scale method is proposed to model cemented HBSs using the discrete element method (DEM). This method takes into account the quantitative effect of hydrate saturation in a sediment pore on the contact bond parameters surrounding the pore. A series of DEM specimens with different macroscopically and microscopically heterogeneous hydrate distributions are prepared. The mechanical behavior of heterogeneous HBSs is investigated by performing biaxial compression tests with flexible boundaries. The simulation results show that both macroscopic and microscopic hydrate distribution heterogeneity can influence the mechanical properties of HBSs. The shear strength is promoted in both macroscopically and microscopically heterogeneous HBSs. Longitudinally heterogeneous HBSs have a higher secant modulus, while transversely heterogeneous HBSs have a lower secant modulus than homogeneous HBSs. The secant modulus of microscopically heterogeneous HBSs first increases and then decreases with increasing pore hydrate saturation. It is found that the deformation behavior and bond breakage evolution of HBSs depend on hydrate distribution heterogeneity. These findings can provide insights into understanding the mechanical behavior of natural HBSs with heterogeneous hydrate distributions.

1. Introduction

Natural gas hydrates are ice-like crystalline solids in which gas molecules (mainly methane) are trapped in cages of water molecules under low-temperature and high-pressure conditions [1,2]. Owing to the specific conditions of formation, they are widely distributed in the pores of permafrost soils and deep-sea sediments [3,4,5,6,7]. The reserves of natural gas hydrates are estimated to be as much as twice other carbon-based fuels including coal, conventional gas, and petroleum reserves [8]. Thus, natural gas hydrates are considered a promising energy resource and have attracted great interest worldwide.

Only several field trial tests of natural gas hydrate production have been carried out to date. Canada carried out trial tests in Mallik from 2007 to 2008 using a staggered depressurization method [9]. Gas was successfully produced for 6 days at a gas production rate of 2000–3000 m³/d. In 2013, Japan conducted the world’s first offshore production trial in Nankai Trough. The production lasted for 6 days, and a total gas of 1.2 × 10⁵ m³ was obtained [10]. In 2017, China carried out an offshore methane gas hydrate production test in the South China Sea. A total volume of 3.0 × 10⁵ m³ gas was produced for 60 days [11]. In 2020, China conducted another field test using a horizontal well. The gas production increased to 8.6 × 10⁵ m³ in 30 days [12]. These field trials lasted for only a short period, and the dissociated area of hydrate reservoirs is relatively small. The mechanical responses of hydrate reservoirs and sediment layers are unclear after long-term exploitation. One of the most critical challenges is the lack of understanding of the mechanical properties of hydrate-bearing sediments (HBSs). The main methods of exploiting gas hydrates are to dissociate hydrates by changing the pressure and temperature conditions [13,14]. The dissociation of natural gas hydrates can cause degradation of the mechanical properties of HBSs, reducing the stability of the seabed slopes. Potential geohazards and environmental problems such as seabed subsidence, wellbore failure, marine landslides, and methane leakage may be triggered [15,16,17,18,19,20,21]. Owing to the difficulty in conducting in situ tests, experimental and numerical methods are used to study the mechanical properties of HBSs. It has been recognized that many factors can affect the mechanical properties of HBSs, such as hydrate saturation, hydrate morphology, confining pressure, and host sand type [22,23,24,25,26].

Hydrate formation is a complex physical and chemical process that can be divided into two stages: hydrate nucleation and hydrate growth [27]. Hydrate nucleation is a stochastic and heterogeneous process at the molecular level. Many hypotheses have been proposed to explain the hydrate nucleation process, including the labile clusters hypothesis, the local structuring nucleation hypothesis, the interfacial nucleation hypothesis, and the blob hypothesis [28,29,30,31,32]. The process of hydrate nucleation can be classified as homogeneous or heterogeneous. Heterogeneous nucleation occurs when impurities exist, which is the most likely phenomenon. After the nucleation stage, hydrates continue growing and exhibit different morphologies at the macroscopic level. The growth rate of hydrate depends on gas consumption, temperature, etc. As a result, the distribution of hydrates in sediments is complex. As shown in Figure 1, the hydrates in natural HBS cores within fine-grained sediments are shown as veins, lenses, and nodules, while those within coarse-grained sediments distribute in sediment pores and are not directly visible. The hydrate morphologies in sediment pores are categorized as pore-filling, cementing, load-bearing, grain-coating, and patchy. Figure 2 shows the coarse-grained HBS specimens synthesized in a laboratory. The hydrates in sediment pores are not single but a combination of two or more morphologies. Moreover, the hydrate saturations in different pores are different, resulting in an obvious distribution heterogeneity. The physical and mechanical properties of heterogeneous HBSs are less investigated than that of homogeneous HBSs. Liu et al. [33] synthesized hydrate-bearing sediment specimens with heterogeneous and homogeneous hydrate distributions using the sand-water mixing method. The results indicate that the hydrate distribution heterogeneity influences the mechanical properties of hydrate-bearing sediments. Therefore, the distribution heterogeneity of hydrates should be taken into account to accurately characterize the mechanical properties of hydrate-bearing sediments, which is fundamental to predicting the mechanical responses during hydrate exploitation.

Numerical modeling is commonly used to investigate the mechanical properties of HBSs. Continuum methods, such as the finite element method (FEM) and finite difference method (FDM), often consider HBSs a continuous and homogeneous material. However, HBSs are a binary material consisting of sediment particles and filled hydrates. The hydrate distribution heterogeneity cannot be efficiently treated using the continuum methods. As a discontinuous method, the discrete element method (DEM) is an efficient numerical tool that can accurately model granular materials with individual particles. The applications of DEM have been increasingly extended to HBSs, including biaxial/triaxial compression tests, direct shear tests, oedometer tests, bender element tests, hydraulic fracturing, and landslide simulations [18,37,38,39,40,41,42,43]. The DEM can provide pore-scale insight into the mechanical behavior of granular materials. It is difficult to take into account the hydrate distribution heterogeneity for laboratory specimens and numerical models; therefore, most previous studies focused on homogeneous HBSs. Liu et al. [44] studied the mechanical responses of pore-filling type of HBSs with macroscopically heterogeneous hydrate distributions. They found that the heterogeneity of hydrate distribution promotes the shear strength, Young’s modulus, and dilation of HBSs. Ding et al. [45] further investigated the mechanical behavior of pore-filling, load-bearing, and cementing types of HBSs with pore-scale heterogeneity of hydrate distributions. The results indicate that the enhancement of heterogeneous HBSs depends on hydrate morphology. In these DEM models, hydrates are represented by spherical particles that cannot quantitatively describe the hydrate saturation due to their geometric nature. Quantitative approaches are needed to accurately characterize the mechanical behavior of HBSs with heterogeneous hydrate distributions.

This study presents a new pore-scale modeling approach using the DEM to investigate the effect of hydrate distribution heterogeneity on the mechanical behavior of HBSs. Pore hydrate saturation is used to describe the cemented HBSs. The relationship between pore hydrate saturation and associated bond parameters is established according to previous experimental data. A series of microscopically and macroscopically heterogeneous HBS specimens of HBSs are created using the proposed approach. The mechanical behavior of HBSs under biaxial compression tests with flexible boundaries is investigated.

2. DEM Modeling of Cemented Hydrate-Bearing Sediments

2.1. Model Setup

A commercial DEM program, Particle Flow Code 2D (PFC2D), is used to investigate the mechanical behavior of HBS specimens [46]. DEM simulations of biaxial compression tests are performed. In the DEM models of HBSs, two-dimensional rigid circular discs are used to represent sediment particles. The dimensions of the initial HBS specimens are 8.0 mm in width and 16.0 mm in height. The HBS specimen consists of 3402 particles. The initial planar void ratio of the HBS specimen is 0.18. The host sediments modeled in this study are Toyoura sands, which have been widely tested and used in artificial methane hydrate-bearing sediments. The particle sizes of sediments range from 0.15 to 0.25 mm following the particle size distribution of Toyoura sands [47], as shown in Figure 3.

Flexible boundary conditions can capture the reasonable deformation behavior of granular materials in DEM simulations [48,49,50,51]. The HBS specimen with flexible boundaries is prepared as follows. First, the initial HBS specimen is consolidated to a small isotropic stress of 10 kPa under lateral rigid walls. Then, the lateral walls are replaced by membrane boundaries consisting of bonded membrane particles. Each membrane particle is constrained by its adjacent membrane particles, as shown in Figure 4. To maintain constant confining pressures during compression, the force applied to one membrane particle (e.g., Particle O) can be calculated by the following equation:

$${\mathit{F}}_{\mathrm{O}}=\frac{1}{2}{\sigma }_3{L}_{\mathrm{AO}}{\mathit{n}}_{\mathrm{AO}}+\frac{1}{2}{\sigma }_3{L}_{\mathrm{OB}}{\mathit{n}}_{\mathrm{OB}}$$

(1)

where ${\sigma }_3$ is the confining pressure, $L_{\mathrm{AO}}$ and $L_{\mathrm{OB}}$ are the distances between a membrane particle and its adjacent particles, respectively, and ${\mathit{n}}_{\mathrm{AO}}$ and ${\mathit{n}}_{\mathrm{OB}}$ are unit normal vectors of the lines between a membrane particle and its adjacent particles, respectively.

2.2. Contact Models for Hydrate-Bearing Sediments

HBSs are binary materials consisting of sediments and hydrates. For the contacts between sediments without hydrates, the rolling resistance linear contact model is adopted to describe the cohesionless nature and consider the particle shape effect of sediments. There are two methods to take into account hydrates in DEM models. In the first method, hydrates are solid particles distributed in the voids of sediment particles. In the second method, hydrates are virtually considered only in the force-displacement relationship. Although previous studies have demonstrated that the first method can reasonably capture the mechanical behavior of HBSs [24,37,39,52], the high computational cost limits the DEM modeling scale. Spherical and circular particles of hydrates have difficulty accurately describing the various patterns of hydrates owing to their geometric nature. The second method emphasizes the cementing effect of hydrates and is used to investigate the mechanical behavior of cementing and grain-coating types of HBSs [53,54,55]. Ding, Qian, Lu, Li, and Zhang [24] found that cementing hydrates can greatly contribute to the shear strength and stiffness of HBSs. In addition, the most commonly used hydrate formation method in a laboratory, the excess-gas method, usually forms cementing hydrates in host sands [56]. As a result, the second method is applicable to capture the main mechanical behavior of HBSs.

One of the difficult problems for the second method is to determine the relationship between hydrate saturation and the enhancement induced by hydrates. Shen and Jiang [55], Jiang et al. [57], and Jiang et al. [58] proposed a bond contact model considering hydrate saturation for HBSs. However, the effect of hydrates is homogeneously considered for a whole specimen, and the heterogeneous characteristics are not considered. To study the effect of hydrate distribution heterogeneity on the mechanical behavior of HBSs, a new method is proposed to describe the hydrate saturation at the pore scale.

Figure 5 presents the hydrate formation process under excess-gas conditions in a pore based on the micro-CT technique. At a low hydrate saturation, hydrates grow around the interface between gas and water. Then, hydrates continue growing from the edge to the center of the gas area, and the pore volume decreases as hydrate forms. At a high hydrate saturation, hydrates can strongly cement the sediment particles. The cementing behavior induced by the cementing hydrates between sand grains is similar to glass beads cemented with epoxy, as shown in Figure 6. The cementing behavior can be described by the parallel bond contact model [46].

Following the above hydrate formation process in a sediment pore, a conceptual model is proposed to describe the cementing hydrates between sediments. The cementing hydrates are uniformly distributed at the sediment–sediment contacts surrounding the pore, as shown in Figure 7a. The volume of hydrates determines the bond strength and stiffness of the sediment–sediment contact. The bonds are weak at low hydrate saturation, and they are enhanced as hydrate saturation increases, as shown in Figure 7b. In DEM models, the pores are characterized by connecting the contacts between particles, as shown in Figure 8. A pore is formed by a few sediment–sediment contacts, and the hydrates are assumed to cement all contacts. Thus, sediment–sediment contact is affected by two adjacent pores, and the parallel bond parameters are determined by the hydrate saturations of the two adjacent pores, as shown in Figure 9a.

Figure 9b shows a representative model to consider cemented hydrates in one pore. The equivalent radius of the two adjacent particles is calculated by Equation (2). The width of the gap between the two particles is w, and the thickness of the hydrates is h. Then, the volume of the cementing hydrates can be obtained by Equation (5). It is a function of the gap width, hydrate thicknesses, and particle diameter. For a given hydrate saturation in a pore, the volume of cemented hydrates for each contact can be calculated while the thickness of cemented hydrates is unknown. It is difficult to directly obtain an analytical solution from Equation (6). Alternatively, an approximate solution can be obtained by using an iterative method. Thus, the thickness of cemented hydrates can be calculated when the volume of hydrates is given. Figure 10 shows that the approximate solutions agree with the analytical solutions, validating the iterative method.

$$r=\frac{r_1+r_2}{r_1{r}_2}$$

(2)

$$\lambda =\frac{h}{r}$$

(3)

$$\eta =\frac{w}{r}$$

(4)

$$V_H=f\left(\lambda ,\eta ,r\right)=\left\{\begin{array}{c}\lambda r^2-\frac{1}{2}\lambda r^2\sqrt{1-{\lambda }^2}-\frac{1}{2}{r}^2\arcsin \lambda +\frac{1}{2}\eta \lambda r^2\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\lambda <1\\ r^2-\frac{1}{4}\pi r^2+\frac{1}{2}\eta r^2+\left(\lambda -1\right)\left(\frac{1}{2}\eta +1\right)r^2\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\lambda \ge 1\end{array}\right.$$

(5)

$$\lambda =f^{-1}\left(V_H,\eta ,r\right)$$

(6)

where $r_1$ and $r_2$ are the radii of two adjacent sediment particles; $r$ is the average radius of the adjacent sediment particles; $\lambda$ is the hydrate thickness to sediment radius ratio; $\eta$ is the gap to sediment radius ratio.

2.3. Relationship between Hydrate Saturation and Parallel Bond Parameters

In the parallel bond contact model, the normal stiffness, tangential stiffness, tensile strength, and cohesion of bonds are essential parameters. To reduce the bond parameters, the normal to tangential stiffness ratio of parallel bonds is equal to that of sediments, and the tensile strength of bonds is equal to cohesion. These parameters are functions of the hydrate thickness to sediment radius ratio λ, as shown in Equations (7) and (8).

$$\left\{\begin{array}{c}{\overline{k}}_n=f_{\mathrm{pb}\_\mathrm{kn}}\left(\lambda \right)\\ {\overline{k}}_s=\frac{{\overline{k}}_n}{\overline{\kappa }}\end{array}\right.$$

(7)

$$\left\{\begin{array}{c}{\overline{\sigma }}_c=f_{\mathrm{pb}\_\mathrm{ten}}\left(\lambda \right)\\ \overline{c}={\overline{\sigma }}_c\end{array}\right.$$

(8)

where ${\overline{k}}_n$ and ${\overline{k}}_s$ are the normal and tangential stiffnesses of bonds, respectively; ${\overline{\sigma }}_c$ and $\overline{c}$ are the tensile strength and cohesion of bonds, respectively; the subscripts “pb_kn” and “pb_ten” denote the normal stiffness and tensile strength of parallel bonds.

A series of DEM simulations are conducted to determine the relationship between hydrate saturation and parallel bond parameters. According to the experimental conditions [25], homogeneous HBS specimens with 25, 40, and 55% hydrate saturations are created, as shown in Figure 11a. In these specimens, the hydrate saturations in all pores are 25, 40, and 55%. The bond parameters are calibrated by the trial-and-error method. Figure 11b shows the bond strength distributions of sediment contacts induced by cemented hydrates. Then, the relationship between the bond parameters and the hydrate thickness to sediment radius ratio is shown in Figure 12. The functions are obtained by fitting the curves of the parallel bond parameters and hydrate thickness to particle radius ratio, as shown in Equations (9) and (10). The parameters of sediment and membrane particles are summarized in

Table 1. Parameters of sediment and membrane particles in DEM simulations.

Parameters	Host Sediments	Membrane	Unit
Sediment particle density	2650	1000	kg/m3
Sediment particle radii	0.15–0.25	0.004	mm
Normal stiffness	9.0 × 108	9.0 × 107	N/m
Normal to tangential stiffness ratio	1.5	1.5	-
Interparticle friction of sands	0.5	0.0	-
Rolling resistance coefficient	0.5	0.0	-
Tensile strength of bonds	Calculated by Equation (10)	1 × 10100	MPa
Cohesion of bond	Calculated by Equation (10)	1 × 10100	MPa
Damping coefficient	0.7	0.7	-
Parallel bond normal stiffness	Calculated by Equation (9)	9.0 × 1011	N/m
Parallel bond stiffness ratio	1.5	1.5	-
Parallel bond radius multiplier	1.0	-	-

.

$$f_{\mathrm{pb}\_\mathrm{kn}}\left(\lambda \right)=1.2297\times {10}^{11}+1.1462\times {10}^{10}{e}^{6.65563\lambda }\left(\mathrm{Pa}/\mathrm{m}\right)$$

(9)

$$f_{\mathrm{pb}\_\mathrm{ten}}\left(\lambda \right)=1.3665+0.0315e^{5.84519\lambda }\left(\mathrm{MPa}\right)$$

(10)

Figure 13 compares the mechanical responses of DEM models and experimental specimens at different hydrate saturations. The peak deviatoric stresses of the DEM models are similar to those of the experimental specimens at different hydrate saturations. Note that the DEM models show much more prominent strain-softening behavior than the experimental specimens. The reason can be that DEM models use circular particle shapes, while experimental specimens comprise irregular particle shapes. When the shear stress exceeds the shear resistance between bonded sediment particles, the circular particles tend to exfoliate more rapidly than irregular particles since the rotation of irregular particles is restricted. The volumetric strain responses of the DEM models agree with those of the lab specimens. The dilation of HBSs depends on hydrate saturation, and the dilation of HBSs increases with an increase in hydrate saturation for both DEM simulations and experiments.

3. Numerical Simulation Results

3.1. Effect of the Macroscopic Heterogeneity of Hydrates

From a macroscopic perspective, natural HBS specimens contain heterogeneous hydrates such as segregated hydrate lenses in fine-grained sediments and patchy or interconnected hydrates in coarse-grained sediments [62]. In this section, HBS specimens with segregated hydrate lenses are investigated to analyze the mechanical behavior of macroscopically heterogeneous HBSs.

3.1.1. Transversely Heterogeneous Hydrate Distributions

A series of DEM specimens of HBSs with transversely heterogeneous hydrate distributions are prepared. Hydrates are distributed in one, two, or four specific layers in the DEM specimens. The average hydrate saturation of each DEM specimen is 40%. The hydrate layer is located at the top, middle, or bottom of the host sediments when hydrates are only distributed in a layer. In addition, two and four hydrate layers are uniformly distributed in the host sediments. Figure 14a illustrates the DEM specimens with different transversely heterogeneous hydrate distributions. In these specimens, the pore hydrate saturations of the hydrate layers are 100%. Figure 14b shows the bond strength of sediment–sediment contacts induced by cemented hydrates.

The mechanical behaviors of the DEM specimens with transversely heterogeneous hydrate distributions are different from those with homogeneous hydrate distributions, as shown in Figure 15. The initial increments of the deviatoric stress of the transversely heterogeneous DEM specimens are similar to those of the homogeneous DEM specimens. The peak deviatoric stresses for different types of transversely heterogeneous DEM specimens are larger than those of homogeneous DEM specimens. Note that the residual strengths of the transversely heterogeneous DEM specimens are larger than those of homogeneous DEM specimens. The transversely heterogeneous DEM specimens first exhibit larger contraction and then exhibit larger dilation than the homogeneous DEM specimens.

3.1.2. Longitudinally Heterogeneous Hydrate Distributions

Longitudinally heterogeneous DEM specimens are created that are similar to transversely heterogeneous DEM specimens. Hydrates are distributed in one, two, and four perpendicular layers, as shown in Figure 16a. The pore hydrate saturations are 100% in the specified perpendicular layers. The average hydrate saturation for each DEM specimen is 40%. For one perpendicular hydrate layer, the layer can be located at the left, middle, and right sides of the DEM specimens. The bonds between sediments caused by cemented hydrates are shown in Figure 16b.

Figure 17 presents the mechanical responses of the DEM specimens with longitudinally heterogeneous hydrate distributions. The deviatoric stress of longitudinally heterogeneous DEM specimens increases more rapidly than that of homogeneous DEM specimens. The peak deviatoric stresses fluctuate for different types of longitudinally heterogeneous DEM specimens. The DEM specimens with longitudinally heterogeneous hydrate distributions show more pronounced dilation than those with homogeneous hydrate distributions.

3.2. Effect of Microscopic Heterogeneity of Hydrates

From a microscopic (pore-scale) perspective, hydrate formation depends on gas supply, temperature, and pressure, resulting in different hydrate saturations in sediment pores. The hydrate distribution can be heterogeneous at the pore-scale level.

A series of DEM specimens with microscopically heterogeneous hydrate distributions are created to investigate the effect of the pore-scale heterogeneity of hydrates on the mechanical behavior of HBSs. In comparison to the homogeneous DEM specimens, the hydrate saturation in each pore space is assumed to be 50, 60, 70, 80, 90, and 100%. The pores are randomly assigned the given hydrate saturations until the average hydrate saturation reaches 40%, as shown in Figure 18a. The bonds induced by pore-scale hydrates are presented in Figure 18b. Note that the bond number decreases while the bond strength increases with increasing pore hydrate saturation.

Figure 19 presents the mechanical responses of the microscopically heterogeneous DEM specimens. The initial developments of deviatoric stress for different types of microscopically heterogeneous DEM specimens are similar. The peak deviatoric stress increases with an increase in pore hydrate saturation. The microscopically heterogeneous DEM specimens experience larger dilation than the homogeneous DEM specimens. The dilation becomes more pronounced with increasing pore hydrate saturation.

4. Discussion

4.1. Shear Strength and Secant Modulus

Figure 20 shows the shear strength and secant modulus of the heterogeneous and homogeneous DEM specimens. The strength of the longitudinally heterogeneous DEM specimens is slightly larger than that of the transversely heterogeneous model. The average strength of the longitudinally heterogeneous DEM specimens is 9.17 MPa, while that of the transversely heterogeneous DEM specimens is 8.17 MPa. The secant moduli of longitudinally and transversely heterogeneous HBSs are 888 MPa and 370 MPa, respectively. The results suggest that the shear strength is promoted, but the secant modulus is reduced for macroscopically heterogeneous DEM specimens. The reason can be explained as follows. The sediments in the hydrate layers of longitudinally and transversely heterogeneous DEM specimens are bonded, and the bonded sediments can prevent particle rotation during shearing to promote shear strength. The hydrate layers in longitudinally heterogeneous DEM specimens are perpendicular to the loading and they are much stronger than weakly bonded sediments to resist large deformation. Compared with homogeneous DEM specimens, the pure sediments in transversely heterogeneous DEM specimens can undergo large deformation. As a result, the secant modulus of the transversely heterogeneous DEM specimens is smaller than that of the homogeneous DEM specimens.

The shear strength of microscopically heterogeneous DEM specimens is smaller than that of macroscopically heterogeneous DEM specimens, as shown in Figure 20a. However, the secant modulus of microscopically heterogeneous DEM specimens is larger than that of transversely heterogeneous DEM specimens but smaller than that of longitudinal heterogeneous DEM specimens, as shown in Figure 20b.

To quantitatively characterize the heterogeneity of HBSs, the degree of hydrate heterogeneity is defined by the following equation [63]:

$$C=\frac{1}{{\overline{S}}_H}\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(S_{Hi}-{\overline{S}}_H\right)}^2}{n-1}}$$

(11)

where $S_{Hi}$ and ${\overline{S}}_H$ are the pore hydrate saturation and average hydrate saturation of the specimen, respectively.

Figure 21 presents the shear strength and secant modulus of HBSs with microscopically heterogeneous hydrate distributions. The shear strength tends to linearly increase with increasing pore hydrate saturation; however, the secant modulus shows different behavior. The number of pores filled by hydrates decreases with an increase in pore hydrate saturation. The secant modulus of HBSs increases with increasing pore hydrate saturation (less than 47%); however, the secant modulus of HBSs decreases with increasing pore hydrate saturation until the pore hydrate saturation is 100%.

4.2. Deformation Characteristics

Figure 22 shows the deformation of transversely heterogeneous DEM specimens during biaxial compression tests. For one hydrate layer, major deformation is from the pure sediments, while only a longitudinal crack occurs after shearing. Two obvious shear bands occur in the DEM specimens when hydrates are distributed at the top and bottom. For the DEM specimen with a middle hydrate layer, the pure sediments at the top and bottom parts undergo large deformation. There are two shear bands for each part of the pure sediments. When hydrates are distributed in two layers, the thickness of one hydrate layer is approximately half of that in the previous case, resulting in a smaller strength of each hydrate layer. As a result, a large shear band occurs throughout the specimen. It is observed that some parts of the hydrate layers are crushed and become part of the shear band. When hydrates are distributed in four layers, the thickness of each layer is smaller than in previous cases. After shearing, the hydrate layers are crushed and two crossed shear bands are formed across the hydrate layers. Figure 23 presents the deformation of longitudinally heterogeneous DEM specimens after biaxial compression. Apparent shear bands in the DEM specimens are observed after shearing. For the DEM specimens with one perpendicular hydrate layer, a dominant inclined shear band is formed. The hydrate layers are crushed during shearing. The thickness of the shear band in the DEM specimen with a middle hydrate layer is narrower than that with a left or right hydrate layer. For the DEM specimens with two and four perpendicular hydrate layers, all hydrate layers are crushed. Figure 24 shows the deformation of microscopically heterogeneous DEM specimens after biaxial compression tests. Similar deformation is observed in the microscopically heterogeneous DEM specimens, and one remarkable shear band is formed after shearing.

4.3. Bond Breakage Evolution

The evolution of bond breakage of heterogeneous DEM specimens is presented in Figure 25. For macroscopically heterogeneous DEM specimens, longitudinally heterogeneous DEM specimens exhibit more broken bonds than transversely heterogeneous DEM specimens. The broken bond number in transversely heterogeneous DEM specimens depends on the position of the hydrate layers. Only a few bonds are broken during shearing when the hydrate layer is located at the top or bottom of DEM specimens, while more bonds are broken when the hydrate layer is located in the middle of DEM specimens. The broken bond number increases with increasing hydrate layer number since the thickness of the hydrate layer decreases with increasing hydrate layer number. Compared with transversely heterogeneous DEM specimens, the broken bonds increase faster at the initial shearing stage as the longitudinal hydrate layer directly resists loading. For longitudinally heterogeneous specimens, the hydrate layers are parallel to the axial loading direction. The hydrate layers can directly resist the axial loading. The total widths of the hydrate layers are similar; thus, the evolutions of bond breakage in longitudinally heterogeneous DEM specimens are similar even though the distributions of hydrate layers are different. More broken bonds occur in microscopically heterogeneous DEM specimens than in macroscopically heterogeneous DEM specimens. The broken bond evolution depends on the pore hydrate saturation. The broken bond number is likely to increase with decreasing pore hydrate saturation. The reason can be that the contact bonds surrounding a sediment pore are strong when the hydrate saturation of the pore is high.

5. Conclusions

This study investigates the effect of macroscopic and microscopic hydrate distribution heterogeneity on the mechanical behavior of HBSs. The main novelty of this study is that a new pore-scale method is developed to accurately model cemented HBSs using the DEM. The hydrate saturation in each sediment pore is taken into account to strengthen the sediment skeleton. This method is capable of capturing the hydrate distribution heterogeneity at the pore scale. A series of DEM simulations are performed, involving macroscopically and microscopically heterogeneous HBS specimens. The main findings are summarized as follows.

• The parallel bond contact model is applied to describe the cementing effect induced by hydrates considering the hydrate saturation at the pore-scale level. The relationship between the hydrate saturation of a sediment pore and the contact bond parameters surrounding the pore is determined based on experimental data. The results demonstrate that this method can reasonably capture the mechanical behavior of HBSs with homogeneous and heterogeneous hydrate distributions.
• Both macroscopic and microscopic hydrate distribution heterogeneity can influence the shear strength and secant modulus of HBSs. The shear strength is promoted for both macroscopically and microscopically heterogeneous HBSs. Macroscopic hydrate distribution heterogeneity can lead to HBSs with a higher shear strength than microscopic hydrate distribution heterogeneity. The secant modulus of longitudinally heterogeneous HBSs is significantly enhanced, while that of transversely heterogeneous HBSs is reduced. The secant modulus of microscopically heterogeneous HBSs first increases and then decreases with increasing pore hydrate saturation.
• The deformation behavior of HBSs depends on hydrate distribution heterogeneity. For transversely heterogeneous HBSs, pure sediments undergo major deformation for one hydrate layer; however, hydrate layers can be broken when hydrates are distributed in multiple layers. The hydrate layers are broken for all cases of longitudinally heterogeneous HBSs. Microscopically heterogeneous HBSs show similar deformation behavior with one remarkable shear band.
• In future work, other influencing factors, such as sediment type, temperature, and capillary pressure, will be incorporated into the DEM model. Furthermore, the mechanical responses of heterogeneous hydrate reservoirs during exploitation can also be investigated.