Discrete Element Simulation Analysis of Damage and Failure of Hydrate-Bearing Sediments

China  LTDS-CNR UMR  , Ecole Centrale de Lyon, Écully  , France Abstract It is of great significance to study the damage and failure law of hydrate-bearing sediments for for exploration and development, as well as for warning secondary disasters such as tsunami and earthquake. The discrete element modeling and simulation method has the advantages of low cost, strong repeatability and accurate response to the microstructure of samples, therefore , the discrete element method is used to simulate and analyze the damage and failure law of hydrate-bearing sediments in this paper. First, a triaxial undrained shear teat model of hydrate-bearing sediments is established. by discrete element simulation software; Then, the effects of different influencing factors on the fracture characteristics of hydrate hydrate-bearing sediments is studied; Finally, the effects of different factors on the initiation stress and damage stress of hydrate sediments are analyzed, and the damage law of hydrate-bearing sediments is obtained. The results show that: (  ) The cementation of hydrate particles is greater than its bearing effect, thus the lithology of sediments changes from loose sand body to brittle rock with the increase of hydrate saturation, resulting in the corresponding change of fracture characteristics from loose sand body failure to brittle rock failure. (  ) With the increase of hydrate saturation, the initiation stress of sediment shows a step-by-step increase law, and the damage stress is positively correlated with the shear modulus. (  ) The heterogeneity of hydrate distribution is closely related to the failure mode of sediments. With the increase of hydrate distribution heterogeneity, the initiation stress displays an increase-decrease-increase pattern while the damage stress shows a law of increasing and then decreasing, which has the high sensitivity to heterogeneity of hydrate that natural gas changed the effective volume the strength. The properties were carried out through numerical and the model was verified by triaxial shear test. his study concluded that the mechanical strength of sediments did not change significantly when the saturation was lower than  .   . High pressure and low temperature triaxial test equipment has been very mature and has been applied to the mechanical properties of hydrate. T. S. Yun studied the relationship between particle size, confining pressure and saturation through triaxial test, their research shows that when the hydrate concentration is low, the mechanical properties of the sediment are only controlled by the stiffness and contact strength of the sediment particles themselves. When the concentration exceeds a certain value (  ), the sediment particles effectively increase the bonding between soil particles and fill the pore space, thus controlling the structural strength and deformation trend of the sediment  . The results of B. N. Madhusudhan. shows that the strength of hydrate-bearing sediments is affected by the specific surface area and particle shape of hydration cementation. The strength of lower specific surface area and coarser particles changes significantly during hydrate dissociation, and changes the deformation mode of sediment  . Shi Shen obtained that the strength and cohesion of the samples increase with the increase of pore pressure and effective confining pressure the consolidation drainage axial test  . Jeong Hoon Choi conducted multi-stage drained triaxial shear test on hydrate-bearing sediments samples prepared in the laboratory. The test shows that sediment strength is positively correlated with confining pressure and saturation, while dilatancy angle is negatively correlated with confining pressure and saturation. Compared with single-stage test, it is considered that multi-stage test make efficient use of core samples  . Miyazaki studied the effect of strain rate on strength of artificial hydrate samples  ,  . Hyodo analyzed the effect of hydrate the mechanical properties of sediments under triaxial shear Clayton conducted a triaxial test of in-situ hydrate by analyzing


Introduction
Due to the depletion of resources and many environmental problems caused by traditional fossil energy, it has been difficult to meet the development needs, but the development of non fossil energy can not keep up with the pace of economic development, according to statistics, the proportion of non fossil energy in the total energy consumption of the world's major economies is still no more than   . Natural Gas Hydrate are known as flammable ice, it is an ice like crystalline material formed by natural gas and water under high pressure and low temperature. Natural gas hydrate is distributed in deep-sea or land permafrost, with rich reserves - , it only generates a small amount of carbon dioxide and water after combustion, the pollution is far less than that of coal and oil, and the reserves are huge. Therefore, it is internationally recognized as the alternative energy of oil, and it has gradually attracted extensive attention from the international community  , Fig.  shows the distribution of hydrate in the whole region. The occurrence conditions of hydrate are extremely harsh, which makes it easy to decompose in the mining process or under the influence of other inducing factors, so as to change the geomechanical conditions of reservoir, lead to reservoir damage and damage, and lead to landslide, earthquake, tsunami and other disasters - .Therefore, it is of great significance to study the effects of different influencing factors on the damage and failure of hydrate-bearing sediments. Because it is difficult to maintain the stable temperature and pressure conditions of hydrate after mining, the sedimentary and mechanical properties of hydrate sediments are mostly tested by physical and laboratory tests. Liu Huaishan realized the accurate prediction of natural gas hydrate reservoir by synthesizing pseudo acoustic curve, and studied the variation laws of resistivity, acoustic velocity, density and gamma with saturation  . W. J. winters studied the effects of different types, different formation methods and pore pressure on the mechanical properties of sediments through acoustic properties  . Gao W. Hu studied the sound velocity characteristics of hydrate-bearing sediments with different saturation, and found that the sound velocity changes little when the saturation is less than . When the saturation exceeds , the hydrate has an obvious consolidation effect on the sediments, but below this value, it is only filled in the pore structure or attached to the sediment particles, which has little effect on the sediment strength  , M. De la fuente introduced a densification mechanism and considered that the invasion of natural gas hydrate changed the effective pore volume of sediments, which affecting the structural strength. The mechanical properties were carried out through numerical simulation, and the model was verified by triaxial shear test. his study concluded that the mechanical strength of sediments did not change significantly when the saturation was lower than .  . High pressure and low temperature triaxial test equipment has been very mature and has been applied to the mechanical properties of hydrate. T. S. Yun studied the relationship between particle size, confining pressure and saturation through triaxial test, their research shows that when the hydrate concentration is low, the mechanical properties of the sediment are only controlled by the stiffness and contact strength of the sediment particles themselves. When the concentration exceeds a certain value (), the sediment particles effectively increase the bonding between soil particles and fill the pore space, thus controlling the structural strength and deformation trend of the sediment  . The experimental results of B. N. Madhusudhan. shows that the strength of hydrate-bearing sediments is affected by the specific surface area and particle shape of hydration cementation. The strength of lower specific surface area and coarser particles changes significantly during hydrate dissociation, and changes the deformation mode of sediment  . Shi Shen obtained that the strength and cohesion of the samples increase with the increase of pore pressure and effective confining pressure through the consolidation drainage axial test  . Jeong Hoon Choi conducted multi-stage drained triaxial shear test on hydrate-bearing sediments samples prepared in the laboratory. The test shows that sediment strength is positively correlated with confining pressure and saturation, while dilatancy angle is negatively correlated with confining pressure and saturation. Compared with single-stage test, it is considered that multi-stage test can make more efficient use of core samples  . Miyazaki studied the effect of strain rate on the strength of artificial hydrate samples ,  . Hyodo analyzed the effect of hydrate dissociation on the mechanical properties of sediments under triaxial shear conditions,and Clayton conducted a dynamic triaxial test of in-situ hydrate sediments by nuclear magnetic resonance  . By analyzing the effects of different fine particle content and density on the mechanical properties of methane hydrate-bearing sediments, Wu Yang concluded that the increase of fine particle content changed the morphology and distribution mode of hydrate between soil particles, affected the cementation strength of natural gas hydrate, and then controlled the mechanical properties of sediments  . In the above research, different methods are used to test the macro mechanical properties of synthetic hydrate-bearing sediments samples from different angles. However, due to the limitation of test conditions, it is difficult to ensure that the artificial samples have the saturation and distribution law similar to that of hydrate-bearing sediments under natural conditions, and it is impossible to study the damage and failure process of hydrate-bearing sediments from a micro perspective.
Discrete element method (DEM) can explain the macro mechanical properties of materials from the micro scale, and has been widely used in the fields of rock and soil mechanics. Ma et al. Studied the mechanical properties of coal, roof and floor composite rock mass (RCF) with different coal thickness by DEM numerical simulation method, and established the damage constitutive model of RCF  . Yang Shengqi et al established the single hole disc splitting model with and without cracks by particle flow code (PFC), studied the strength and deformation law of the disc, and analyzed the evolution characteristics of micro cracks from the meso mechanism ,  . Many scholars have introduced DEM into the study of mechanical properties of natural gas hydrate-bearing sediments. Helgerud et al. simulated three different occurrence forms and mechanical properties of hydrate in sediments by particle flow program - . Wenrong Feng et al. established the model of deep-water and shallow gas hydrate reservoir, studied the borehole stability of horizontal and vertical wells at different depths, and compared it with the analytical solution to verify the reliability of DEM method  . Salinity is an important factor affecting the stability of sediments, Jiang extended the constitutive model of hydrate particle cementation, and quantitatively analyzed the relationship between salinity and mechanical properties of hydrate sediments by DEM  . Earthquakes and tsunamis affect the stability of hydrate-bearing sediments, under cyclic loading, once the hydrate dissociates, it may lead to submarine landslide, Wang et al. Conducted several groups of equal stress drainage cyclic loading biaxial tests through DEM, and analyzed the effect of cyclic loading on the mechanical properties of sediments from a micro perspective  . Jiang studied the effects of temperature and water pressure on the mechanical properties and strain of hydrate-bearing sediments by using the thermal-water-mechanical contact model in DEM ,  . In his further research, he coupled computational fluid dynamics and DEM to study the effect of earthquake on hydrate bearing seabed landslide ,  . Hydraulic fracturing can effectively transform hydrate reservoirs, Yao used particle flow simulation software to test the influence of water pressure and saturation on fracture propagation law  . The above numerical simulation shows the successful application of DEM in different directions of mechanical properties of hydrate-bearing sediments, but there is still a lack of understanding of the process from cracking to damage of hydrate-bearing sediments affected by different influencing factors. Crack initiation stress and damage stress are important characteristics of rock strength, and they are also important nodes in different stages of rock damage and failure process. The study on the initiation stress and damage stress of hydrate-bearing sediments is of great significance for understanding its failure mechanism and progressive failure process and early warning of secondary disasters caused by hydrate decomposition.
The existence of hydrate makes the sediments have the characteristics of brittle rocks, and the brittleness increases with the increase of saturation, and the corresponding failure mode changes  . The occurrence characteristics of hydrate lead to its uneven distribution in sediments. according to previous studies - , the heterogeneity of hydrate distribution is an important factor affecting the mechanical properties of sediments  . In this paper, we established the triaxial undrained shear test model of hydrate-bearing sediments by using three-dimensional particle flow discrete element numerical simulation software (PFCD).and then designed different working conditions to anslyze the effects of hydrate saturation and hydrate distribution heterogeneity on sediment failure mode, initiation stress and damage stress.

Numerical simulation process
Particle flow discrete element method.
Cundall and potyondy studied the theory of granular media, and developed the particle flow dispersion element method on this basis , . The particle flow dispersion element method simplifies the medium into particles and simulates the motion and interaction of medium particles. Particles are represented as spheres, in the particle flow discrete element simulation software, which can be stacked with each other to form blocks  . The load is realized by the displacement of the wall element. The interaction between particles and between particles and walls is realized through contact. Different contact models determine the form and size of the action. The successful operation of the program is based on the following assumptions  . () The particle element is a rigid body without deformation; () The contact between particles is flexible, and a certain amount of overlap is allowed, but the amount of overlap is relatively small relative to the particle size, and the amount of overlap is related to the force; () The contact between media can be bonded; () Spherical particles can be overlapped and combined into different condensation units.
This simulation program is calculated by the display difference method. In a sufficiently small time and short step, the particle motion only affects the adjacent particles. Only the units in contact with the particle affect the force acting on the particle. The calculation process includes two aspects  : () The motion of particles is described by Newton's second law, and the velocity and acceleration of particles are determined according to the resultant force and torque; () The contact force between particles is calculated according to the force displacement law. The contact is determined by the normal and tangential stiffness and corresponding displacement, as shown in Fig   g s is the normal contact interval. k n and k s are normal and tangential. A is the indirect contact surface of particles stiffness respectively. M t is the torque.
J is the polar moment of inertia between the keys. s is the tangential displacement.
M b is the contact bending moment. I is the moment of inertia.  b is the rotation angle of the bonding contact surface. The maximum normal contact force and maximum shear contact force on the contact surface can be calculated from the contact force and contact moment (including bending moment and torque). The calculation formula is as follows  :  In this paper, the triaxial test model space of three-dimensional cylinder is established，the diameter of the cylinder sample is mm and the height is mm. In this model space, the boundary effect can be ignored. Initially, each wall property is set to rigid frictionless. Refer to hydrate characteristics in Shenhu sea area of the South China Sea  , Sediment particles with porosity of . and particle size of .-.mm are generated according to particle gradation, and  particles are finally generated. As shown in Fig. , the particle size distribution curve of sediment particles is generated in the simulation. Set the particle friction coefficient to . and the damping coefficient to ., and then control the servo mechanism to consolidate the sediment particles under the confining pressure of MPa. Generally, there are two methods of hydrate generation. One is to randomly produce gas hydrate immediately after sediment generation and run under zero constraint condition to reach equilibrium, and then to the confining pressure of MPa to solidify particles. The other is to set the confining pressure to MPa after the formation of sediment to consolidate the sediment, then randomly generate hydrate between sediment pores by using the stepwise expansion method, and then consolidate the sediment and hydrate by maintaining the servo pressure of MPa. Brugada simulated the two preparation methods by DEM, and found that the bias stress peaks of the samples generated by the two preparation methods were basically the same  . However, at the same saturation, the initial stiffness of the first preparation method is greater than that of the second. The hydrate produced by the first preparation method mainly acts in the early stage of strain. In fact, the transformation process of sediments from brittleness to plasticity is greatly affected by hydrate. At the same time, considering that hydrate is formed after methane gas migrates into the framework pores of deep-sea shallow soil with initial bearing capacity through appropriate channels,thus ,the second method is used in this paper.  . The small blue particles in the sample are hydrate particles, the number is , and the particle size is .mm, It can be seen from the figure that the area with large particle displacement in the damaged sample corresponds to the fracture concentration area. Therefore, the failure mode of the sample can be well evaluated through the fracture and particle displacement nephogram. The bond properties between hydrate-hydrate and hydrate-sediment are set as parallel bond model to realize the cementation of hydrate, and the bond properties between sediment and sediment are set as linear contact model, Considering the irregularity of sediment particle surface, the characteristic of anti-rolling resistance is added between sediment-sediment and fractured parallel bonding bond, and the particle attribute is set to inheritable.
The density values of sediments and hydrates in the model refer to the hydrate and sediment characteristics in Shenhu sea area of the South China Sea , . The micro parameters of the model were adjusted by trial and error method to make the macro parameters of the model close to the mechanical parameters in the laboratory, so as to calibrate the parameters of particle mechanical properties and particle contact properties, as shown in table  and table . After the assignment of particle and contact properties, start the triaxial loading test, maintain the servo confining pressure of MPa on the side wall, turn off the servo of the upper and lower walls, and set the normal phase stiffness of the side wall to . times the average stiffness to simulate the soft rubber sleeve in the laboratory, then set the loading speed to mm / min.

Model parameter verification.
In order to verify the reliability of the model, the sediment particle size was set to .-.mm and the confining pressure was MPa, and then the load was carried out. The results obtained were compared with the previous experimental and simulation results ,, , as shown in Fig. . The sediment particle size gradation generated by this simulation is basically consistent with the sample particle size gradation prepared in the test. The trend of stress-strain curve after the compaction stage is more consistent with the test, showing the post peak stage more completely. The deviation of the curve may be due to the limitation of model size. The reliability of the model is verified by comparison with previous studies. Simulate the loading process after model parameter verification: Firstly, fix the lower wall of the sample and apply a falling speed of mm / s to the lower wall until the post peak strength reaches  of the peak stress after the sample is damaged. Secondly, the shear stress-strain relationship, volume strain, number of cracks, number of acoustic emission events and failure mode during loading are recorded. Then change the hydrate saturation and hydrate distribution, and repeat the above two steps. Finally, the damage law of hydrate bearing sediments is analyzed according to the factors recorded during loading under different hydrate saturation and hydrate distribution conditions.

Results and Discussion
Determination of initiation stress and damage stress.
Scholars have done a lot of research on the determination of rock initiation stress and damage stress, such as direct observation method based on CT, SEM and other observation equipment ,  、acoustic emission method ,  、method based on crack and volume strain ,  , however, a unified standard has not yet been formed. Li Cunbao  conducted theoretical and Experimental Research on the initiation stress and damage stress of Shale Based on fracture mechanics. He believed that it was reliable to determine the initiation stress and damage stress by crack and volume strain, and established the corresponding model，On this basis, Yin Dawei  used PFCD program to simulate and analyze the method of crack and volumetric strain, and obtained that taking the axial stress corresponding to the rapid growth point of crack as the damage stress is more consistent with the test results. In this paper, the initiation stress and damage stress of sediment samples are determined by monitoring microcracks. Generally, rock compression deformation can be divided into five stages. Under the action of loading stress, the internal natural microcrack is closed, which corresponds to the first stage of rock deformation, i.e. compaction stage; The second stage of rock deformation is linear elastic stage, in which the stress-strain relationship is linear, and the deformation recovers after unloading; When the loading is continued, the microcrack begins to sprout and expand stably in the rock, the rock appears plastic deformation, the curve deviates from the linearity, and the starting point of this stage corresponds to the crack initiation stress( i ); When the micro cracks in the rock propagate and penetrate one step, the number of cracks shows an unstable and rapid growth trend, and the rock failure can not be recovered. The starting point of this stage corresponds to the damage stress( d ), and the peak stress corresponding to termination point ; The fifth stage is the With the increase of hydrate saturation, the lithology characteristics of sediments change, the brittleness characteristics become more and more obvious, and the characteristics of corresponding deformation stages change accordingly.
In Fig. , the deviatoric stress-strain characteristics, the number of microcracks and the number of acoustic emission events during loading are shown under the condition of hydrate saturation of  -, then stop loading when the post peak stress accounts for  of the peak stress.On the whole, With the increase of hydrate saturation, the shear yield of sediments occurs between  - of axial strain, and the stress required for further deformation decreases, that is, the strain softening characteristic is more obvious. When the hydrate saturation is low, the number of cracks and acoustic emission events are less, with the reason that the sediments show the characteristics of loose sand body. In addition, with the increase of hydrate saturation, the brittle characteristics of sediments increase, and the number of corresponding microcracks and acoustic emission events increase.
As shown in Fig. , When the hydrate saturation is , the first and second stages of sample deformation experience a long time, the first microcrack occurs when the axial strain reaches .. When the hydrate saturation is , the first and second stages of sample deformation experience a shorter time, and the axial strain is .. However, peak stress and total axial strain of the sample have not changed significantly at the point of crack growth rapidly , indicating that when the water core saturation is less than , the hydrate particles only reduce the porosity of the sample by filling, the bonding between sediment particles is still linear contact, and the failure of the sample is mainly in the form of loose sand body. The time interval between crack initiation and rapid growth increases first and then decreases when the hydrate saturation is  -, this is mainly because with the increase of hydrate saturation, hydrate particles are no longer just filling and pores, but show the role of cemented sediment particles, which increases the strength of the sample. But, due to the low hydrate concentration and the small number of cemented bonds, the strain softening phenomenon is obviously enhanced at this stage. When the hydrate saturation exceeds , with the increase of saturation, the time of the first and second stages of sample deformation increases, the time interval between crack initiation and growth rapidly increases, the amplitude of peak stress increases, and the time after peak decreases. In Fig.  , with the increase of hydrate saturation, the peak value of acoustic emission events gradually moves towards the direction of increasing strain. Acoustic emission events mainly occur in the post peak stage, indicating that with the increase of hydrate saturation, the brittleness characteristics of sediment samples become more and more obvious. The heterogeneity of hydrate distribution has an obvious influence on the mechanical properties of sediments, so different hydrate particle distribution models are designed to study its influence law. According to the research of Liu  et al., the influence of the longitudinal distribution heterogeneity of hydrate is significantly greater than the transverse distribution heterogeneity. Therefore, only the longitudinal distribution heterogeneity is studied. Similarly, it is set to stop loading when the post peak strength is  of the peak stress. Firstly, we assumed that when the hydrate saturation is , the hydrate particles are distributed in the bottom , bottom  and bottom  of the sediment sample, and the comparison is made with the case of uniform distribution. The strain curve, the number of microcracks and the number of acoustic emission events are shown in of ., showing a trend of first decreasing, then increasing and decreasing finally. As the hydrate particle distribution tends to be uniform, the time interval between crack initiation and growth rapidly increases first and then decreases, the axial strain decreases, while the peak stress increases first and then decreases.
Typical failure law of hydrate-bearing sediments. Different hydrate saturations and hydrate particle distributions affect the failure mode of sediment samples, and the changes of initiation stress and damage stress are closely related to the failure mode  . PFCD program can monitor the particle displacement, crack generation and crack surface propagation and penetration during the loading process of the model, which provides convenience for us to understand the progressive damage and failure of the sample from the microscopic point of view. Fig  and Fig  show the failure modes of samples with hydrate saturation of  - and heterogeneous particle distribution of hydrate saturation of , respectively.
It can be seen from Fig.  that, there are fewer cracks in the sample, and the number of microcracks increases with the increase of hydrate saturation when the hydrate saturation is low. With the increase of hydrate saturation, the particle displacement of the damaged sample generally shows a decreasing trend, and the particle displacement at both ends of the sample is greater than that in the middle.
When the hydrate saturation is less than , the microcracks are evenly distributed in the sample, and there is no obvious fracture surface after the sample is damaged.
When the hydrate saturation exceeds , the damaged sample produces an obvious fracture surface, and with the increase of saturation, the fracture surface becomes more obvious, mainly one main fracture surface. From the sample with saturation of  -, it can be seen that the included angle between the fracture surface and the loading direction is about  °. From the microcrack profile, the distribution of cracks in the sample is basically consistent with the morphology of the fracture surface. The number of tensile cracks produced by the sample is significantly more than that of shear cracks, indicating that the failure of the sample is mainly tensile failure, which is more obvious when the hydrate saturation is low. At the same time, it can be seen that the shear cracks are mainly distributed on the main fracture surface, indicating that the failure of the sample is caused by the penetration of cracks caused by a large number of shear cracks. in which the first row is the crack surface after sample failure, the second row is the particle displacement nephogram, and the third row is the microcrack profile.
Fig.  shows that the particle displacement increases gradually from the hydrate particles distributed in the bottom 3 to the bottom 3. When the hydrate particles are evenly distributed, the particle displacement decreases. At the same time, the particle displacement in the area without hydrate distribution is greater than that in the area with hydrate histribution. When the hydrate particles are distributed in the bottom 3, the microcracks of the sample are mainly distributed on the regional interface. As the hydrate particle distribution tends to be uniform, the microcracks expand in the direction of hydrate particle distribution. From the hydrate particles distributed in 3 to 3 of the bottom, it can be seen that the microcrack distribution of the sample is more and more concentrated, the fracture surface is gradually clear, and the main fracture surface appears, and the included angle between the main fracture surface and the loading direction is about °. Compared with the uniform distribution and distribution of hydrate particles in the bottom 3, when hydrate particles are distributed in the bottom 3, the density of hydrate particles is large, and the number of bonding bonds between hydrate particles and sediment particles is large, which increases the strength of the sample, and the brittleness characteristics are more obvious. Therefore, the fracture surface generated is clear and concentrated. Through the above analysis, it can be seen that the heterogeneity of hydrate distribution has a significant impact on the failure mode of the sample. The failure form of the area where hydrate particles are not distributed in the sample is the failure of loose sand body. The intrusion of hydrate particles makes the sample have brittle characteristics, and the greater the density of hydrate particles, the more obvious the brittle characteristics are.

Effects of different saturation on initiation stress and damage stress of sediments.
Through the analysis of the above two sections, the determination methods of initiation stress and damage stress under different hydrate saturation and distribution conditions are determined, and the micro damage and failure law of the sample is more clearly understood by analyzing the typical failure mode of the sample. The saturation of hydrate significantly affects the mechanical characteristics of sediments. In Fig.  (a), the sample volume shows the phenomenon of early shrinkage and later expansion. When the hydrate saturation is less than , the increase of saturation increases the volume strain, because the hydrate particle invasion increases the friction between sediment particles and increases the peak strength and post peak residual strength of the sample. Corresponding to Fig.   (a), the hydrate saturation increases from  to , and the ratio of initiation stress to peak stress decreases from . to ., which may be caused by the increase of hydrate particles reducing the anti rolling resistance between sediment Volume strain ratio (%)

Hydrate Saturation (%)
Elastic modulus of shear (b) particles and the weak cementation of hydrate particles. When the hydrate saturation is more than  but less than , the volumetric strain decreases with the increase of hydrate saturation, and the initiation stress increases and increases rapidly. When the hydrate saturation is , the volumetric strain increases, and the rate of initiation stress increases decreases. This may be because the cementation of hydrate particles gradually increases and the bearing capacity gradually decreases with the increase of hydrate saturation. As can be seen from Figure  (a), when hydrate saturation exceeds , with the increase of saturation, the overall initiation stress shows a step increase trend, and the ratio of initiation stress to peak stress shows a periodic change rule of increase-decrease-increase-decrease, and the ratio is between .-..This may be related to the number of cemented bonds formed between hydrate particles and sediment particles, that is, the microscopic structure of hydrate particles is inconsistent with the size distribution of sediment particles and the size of specific surface area, resulting in periodic changes of cemented bonds formed between them.  shows that the shear modulus of the sample gradually increases, and the shear modulus increases from Mp to Mp, with an increase range of , indicating that the strength of the sediment sample increases significantly with the increase of hydrate saturation. As can be seen from Fig.  (b), with the increase of hydrate saturation, the damage stress of the sample gradually increases, and the ratio of damage stress to peak stress also increases with the ratio between .-.. At the same time, it can be concluded that there is a positive correlation between damage stress and shear modulus. The stress state at the damage point corresponds to the unstable propagation of microcracks in the sample. From this stage, irreversible plastic damage occurs. The failure mode of hydrate-bearing sediments is closely related to the initiation stress and damage stress. In order to deeply understand the influence of hydrate distribution on the damage and failure law of sediments, it is necessary to analyze the influence of hydrate distribution heterogeneity on the initiation stress and damage stress of sediments. It can be seen from Fig.  (a) that when the axial strain is less than ., the influence of different hydrate distribution on the volumetric strain is not obvious. When the axial strain exceeds ., the volumetric strain of hydrate distributed in the bottom 3s is greater than that of hydrate evenly distributed, and the volumetric strain of hydrate distributed in the bottom 3 is greater than that of hydrate distributed in the bottom 2. When the specimen is damaged, the difference between the maximum volumetric strain and the minimum volumetric strain is . Correspondingly, in Fig.  (c), the ratio of initiation stress to peak stress is between .-.. With the increase of hydrate distribution heterogeneity, the initiation stress shows the law of increase decrease increase, that is, the volume strain has a good corresponding relationship with the initiation stress. The reason for the above phenomenon may be that with the change of hydrate distribution heterogeneity, the volume proportion, shear strength and failure mode of hydrate containing area and hydrate free area are different; With the increase of the heterogeneity of hydrate distribution, the shear modulus of the sample decreases gradually, indicating that the more uniform the hydrate distribution is, the stronger the resistance to failure of the sample is( Fig.(b)). According to Fig.  (d), the damage stress increases first and then decreases, and the ratio of damage stress to peak stress increases first, then decreases and increases finally, that is, the damage stress has no obvious corresponding relationship with peak stress and shear modulus, but is related to specific failure mode. The ratio of damage stress to peak stress is .-., which indicates that the damage stress of the sample is sensitive to the distribution of hydrate.

Conclusions
In this paper, the triaxial undrained shear test of hydrate-bearing sediments under different hydrate saturations and hydrate distributions by numerical simulation is established firstly. Then, the effects of hydrate saturation and distribution heterogeneity on the initiation stress, damage stress and failure mode of the sample are studied by analyzing the shear stress-strain relationship, the number of tensile and shear cracks, the law of crack propagation and the number of acoustic emission events respectively. On this basis, the damage law of hydrate containing sediments is obtained. The main conclusions are as follows.
()The triaxial undrained shear test of sediment samples under different saturations and distributions conditions of hydrate is simulated. Through the analysis of the characteristics of partial stress and strain, crack number and acoustic emission time in different deformation stages, it is found that the number of cracks has an obvious corresponding relationship with the deformation stage of the sample. Therefore, the method for determining crack initiation stress and damage stress in sediment samples is determined by monitoring the number of cracks.
()Hydrate particles have cementation on sediment soil particles. With the increase of hydrate saturation, the number of cementation bonds in sediments gradually increases, the strength of samples increases, and the lithologic characteristics of samples change from loose sand body to brittle rock; When the hydrate concentration is lower than , the hydrate particles mainly play the role of bearing, but when the hydrate concentration exceeds , the hydrate particles mainly play the role of cementation; The more homogeneous the distribution of hydrate particles, the smaller the density of hydrate particles, the weaker the bearing effect of hydrate, the stronger the cementation of hydrate particles, and the stronger the sample strength, indicating that the cementation of hydrate particles has a greater impact on the mechanical strength of sediments than the bearing effect of hydrate particles filling pores.
()With the increase of hydrate saturation, the number of microcracks generated during the failure process of the sample gradually increases, and the number of tensile cracks is significantly greater than that of shear cracks, indicating that the sample is mainly tensile failure; More obvious main crack occurs with the higher the hydrate saturation and it is roughly at an angle of  ° with the loading direction. In addition, the shear cracks are mainly concentrated near the main crack, showing that the final failure of the sample is due to the initiation, propagation and penetration of a large number of shear cracks; It is concluded that the failure in the hydrate area of the sample is the form of loose sand body failure, and the hydrate area is the form of brittle rock failure by the hydrate distribution heterogeneity test.
()When the hydrate saturation exceeds , the initiation stress of the sample is positively correlated with the hydrate saturation, In addition, the increase of saturation results in a step-by-step increasing trend of initiation stress and a periodic trend of increase-decrease-increase-decrease of the ratio between .-. of initiation stress to peak stress; With the increase of hydrate saturation, there is a good positive correlation between damage stress and shear modulus, and the ratio of damage stress to peak stress ranges from . to ..
()The initiation stress shows the law of increase decrease increase when hydrate distribution heterogeneity increase, and the volumetric strain has a good corresponding relationship with the initiation stress; With the increase of hydrate heterogeneity, the shear modulus decreases gradually, while the damage stress increases first and then decreases. The ratio of damage stress to peak stress is .-., which indicates that the damage stress of the sample is sensitive to the distribution of hydrate.

Data availability
The data used to support the findings of this study are available from the corresponding author upon request.