Abstract
Deep-sea hydrate has great commercial exploitation value as a new type of energy, due to huge reserves, wide distribution, cleanliness, and lack of pollution. Accurately, prediction of the mechanical properties of hydrate reservoirs is a key issue for safe and efficient exploitation of deep-sea hydrate. Although there have been some experimental and numerical simulation studies on the borehole stability of the hydrate layer, the influence of temperature and flow on the decomposition of reservoir hydrate is still not well understood. There have been few pure mechanical studies on the stress and strain state of the hydrate formation around the well, and it is impossible to intuitively understand the influence of the wellbore on the original stress state of the hydrate formation. This paper therefore uses a discrete element method to establish a deep-water shallow hydrate reservoir borehole stability model and compares the discrete element numerical model with an elastoplastic analytical model of borehole stability to verify the reliability of the numerical model. A simulation study on the influence of factors such as reservoir depth and hydrate saturation on wellbore stability is carried out. The simulation results effectively present the constitutive characteristics of strain softening of hydrate sediments. According to the different mechanical characteristics, the near-well zone can be divided into a plastic strain softening zone, a plastic strain hardening zone, and an elastic zone. Reservoir depth and hydrate saturation are found to change the stress state near the well. The greater the depth and the lower the hydrate saturation, the greater the borehole shrinkage. The diameter of the optimal horizontal well in the goaf is in the range from 0.6 to 1.2 m.
1. Introduction
Natural gas hydrate (NGH) is a clathrate crystalline compound formed by water molecules enveloping methane molecules in the form of hydrogen bonds. It exists widely in deep-water shallow strata and permafrost zones [1]. One cubic meter of NGH can release 164 cubic meters of methane gas and 0.87 cubic meters of water [2]. According to Kvenvolden’s statistics, the proven carbon reserves of NGH in the world are more than twice the carbon reserves of conventional fossil fuel energy sources (coal, oil, and natural gas) [3]. Therefore, NGH is currently one of the most promising energy types. With the advancement of deep-water oil and gas resource development, how to efficiently, economically, and safely exploit deep-water shallow hydrates has attracted increasing attention.
Borehole instability is one of the main problems faced in the drilling and production of gas hydrate formations. Key factors affecting borehole stability are the petrophysical and mechanical properties of the hydrate formation and the stability of the hydrate in the formation. The stability of NGH is a core factor. Since the stability of hydrate depends on certain temperature and pressure conditions and the exploitation of hydrate will inevitably disturb or even destroy the original temperature, pressure, and stress state of the sedimentary layer, study of the borehole stability of the hydrate layer is a key to effective exploitation. Due to the influence of gas hydrate characteristics and geological environment, this kind of drilling will face more complex problems than general oil and gas drilling. Before drilling, the formation is in pressure balance state. But after drilling, the original support force is lost, and then, the cementing strength becomes weak, which will also lead to the instability of borehole wall. In particular, under the influence of temperature and pressure, the decomposition of the supporting solid hydrate will cause the borehole to collapse. The increasing of water content will further threaten stability of the hydrate.
There has been a series of studies investigating the borehole stability of the hydrate layer. Thermal and poroelastic effects are preferably considered to estimate wellbore stability [4]. Birchwood et al. proposed an elastic-plastic wellbore stability prediction model based on the Mohr-Coulomb criterion, which took into account the effect of temperature on the thermodynamic state of the hydrate layer [5]. By studying the influence of the mud circulation rate and sedimentary salinity on the decomposition of hydrate, they found that the mud circulation rate was the most critical factor in keeping the hydrate stable. Freij-Ayoub et al. established a numerical model of hydrate borehole stability coupled with thermodynamic effects. The simulation results showed that when the drilling fluid temperature was 5°C higher than the reservoir, the yield zone around the borehole expanded by 32% [6]. Rutqvist and Moridis used numerical simulations to study the influence of deep heat flow intrusion, gas production, and the weight of mining equipment on the mechanical state of the hydrate layer [7]. Their results showed that heat flow intrusion and increasing gas production weakened the mechanical stability of the hydrate layer. Mining equipment above the seabed can increase reservoir pressure and help maintain the mechanical stability of the reservoir. Gao and Gray analyzed the wellbore stability through a coupled geomechanics and reservoir simulator [8]. Khabibullin et al. proposed a one-dimensional semianalytical model to describe the transfer of heat and fluid and coupled it into a numerical model for temperature field prediction around the well [9]. They found that the amount of hydrate decomposition depended on the initial reservoir characteristics and the bottom hole temperature and pressure conditions [9]. Kou et al. established a mathematical model of wellbore stability in hydrate formations considering factors such as heat conduction and hydrate decomposition [10, 11]. The simulation results revealed the influence of drilling fluid temperature on the decomposition of hydrate and indicated that the thermal decomposition of hydrate led to deterioration of the mechanical environment of the reservoir. Rutqvist et al. established a numerical model of coupled multiphase flow to study the effect of decompression mining methods on the mechanical states of two different sedimentary layers [12]. They found that the gas production rate and bottom hole pressure drop determined the pressure state of the entire reservoir and also changed the mechanical state of the near-well zone. Chong et al. studied the impact of horizontal well mining on hydrate production with experimental equipment and found that horizontal well mining could increase gas production while reducing water production [13]. Yu et al. studied the process of hydrate decomposition and gas intrusion into reservoirs through experiments and numerical simulations [14]. The experimental results showed that drilling fluid temperature, hydrate saturation, and reservoir pressure were the main factors leading to hydrate decomposition and gas invasion. The numerical simulation results showed that the use of low-temperature drilling fluid and low circulating rate was beneficial for controlling gas intrusion into the reservoir. Sun et al. conducted a numerical simulation analysis of borehole stability based on the characteristics of the hydrate layer in the Shenhu area of the South China Sea [15]. The results showed that the thermal effect of drilling and the high salinity of the drilling fluid caused the release of free gas in the reservoir, which led to an increase of the pore pressure in the near-wellbore zone. Controlling the salinity of the drilling fluid could effectively control the generation of free gas and prevent borehole wall instability.
Although there have been some experimental and numerical simulation studies on the borehole stability of the hydrate layer, the influence of temperature and flow on the decomposition of reservoir hydrate is still not well understood. There have been few pure mechanical studies on the stress and strain state of the hydrate formation around the well, and it is impossible to intuitively understand the influence of the wellbore on the original stress state of the hydrate formation. Considering the effect of stress on rock strength is better to understand the stability of the borehole [16]. Connected to collapse and fracture gradients is the stress around the wellbore [17]. The drilling process disrupted the stress balance of the formation rock and caused the redistribution of stress around the borehole. The stability of the formation will be expressed through stress. When the stress is unbalanced, the formation will be destroyed. Determining the hydrate formation stress is the prerequisite for studying the stability of the borehole wall. The research results are always uncertain because of the complex working conditions and the dynamic changes of underground stress, which makes it difficult to solve the problem of borehole stability.
This paper therefore uses discrete element method to establish a numerical model of deep-water shallow hydrate reservoir borehole with actual stratigraphic environment based on strain softening characteristics to capture the stress and strain states of the borehole, which are used to describe the stability. For verifying the reliability of the numerical model, a comparison of borehole stability in the discrete element numerical model and an elastoplastic analytical model was conducted. It also simulates the influence of factors such as reservoir depth and hydrate saturation on wellbore stability in order to provide more theoretical support for the drilling and production of deep-water shallow hydrates.
2. Borehole Stability Model for Deep-Water Shallow Hydrate Reservoirs
2.1. Analytical Model of Borehole Stability Based on Strain Softening Characteristics
Timoshenko and Goodier deduced elastic mechanics expressions of radial and circumferential stress of a cylindrical vessel with uniformly distributed internal and external forces [18]. Yu further considered the special case when the outer boundary of the cylinder tended to infinity, and the stress expression at this time was the elastic solution of the horizontal isotropic formation stress distribution [19]. Therefore, equations (1) and (2) give the elastic zone stress solution of the borehole wall stability model in the cylindrical coordinate system, and equation (3) gives the radial displacement expression: where is the horizontal stress, is the radial stress at the elastoplastic boundary, is the radial coordinate of the elastoplastic boundary, is the radial coordinate of a point, is the elastic modulus, is Poisson’s ratio, is the radial stress, is the circumferential stress, and is the radial displacement.
Chen and Abousleiman took the strain softening model as the constitutive model of the formation and gave a plastic model of horizontal isotropic formation stress distribution [20]. Based on the Drucker-Prager criterion, a strain softening constitutive model in the form of internal friction angle was given as: where is the initial internal friction angle, is the strain, and and are the parameters related to the constitutive model.
Based on the research of Chen and Abousleiman, the plastic zone stress differential expressions of the hydrate formation borehole stability model are:
According to previous research of Chen and Abousleiman [21, 22], the boundary conditions of the elastoplastic interface are: where is the ratio of radial displacement to coordinates, (; ) and are the transition parameters, is the volumetric strain, and is the initial average effective stress. Thus, equations (1) and (2) are the stress expressions in the elastic zone of the formation around the well, equations (5)–(7) are the stress differential expressions in the plastic zone, and the boundary conditions of equations (8)–(11) are added together to form the wellbore stability model of hydrate formation.
2.2. Numerical Model of Borehole Stability Based on the Discrete Element Method
This paper uses the discrete element numerical simulation method to simulate the stability of a horizontal borehole during the mining process. The discrete element method was first proposed by Cundall and Strack to study the mechanical behavior of rocks and soils [23]. In discrete element simulation, the model body is composed of particles. After setting microscopic parameters (such as stiffness and bonding strength) for particles and contact points and applying external forces, the movement and collision between particles can produce different mechanical responses of the model body. When the mechanical response conforms to the real mechanical behavior of the material, the microscopic parameters are applicable to the material.
Figure 1 shows the discrete element model of the borehole stability. The yellow particles are sand, and the blue particles are hydrates existing in the pores of the sediment. The model is a cuboid space surrounded by six invisible walls with a length and width of 3 m and a depth of 1 m. The “particle expansion method” is adopted to successively generate sand particles and hydrate particles [24] to ensure that all particles can be stably generated within the specified area and conform to the pore characteristics of real sediments. Subsequently, the servo mechanism is adopted to control the movement of the wall [25] so that the model reaches the consolidation state in three directions with the horizontal maximum and minimum ground stress and vertical ground stress.

To simulate the influence of the borehole on the hydrate formation, a through-hole model is selected, and all particles in the range are deleted. By establishing a cylindrical wall that fits the inner wall of the borehole, the servo mechanism is used to apply mud pressure on the inner wall of the borehole. When the inner wall pressure of the wellbore reaches a stable state, the borehole deformation and the stress state around the wellbore reflect the model’s prediction of the real situation. Since the discrete element model simulates the effective stress state, the ground stress and mud pressure of the model are the differences between the actual stress and the pore pressure.
To simulate the composition of sediments, the established model takes sand particles with a porosity of 0.5 and a radius of 4 cm as the sediment skeleton, hydrate particles with a radius of 1.5 cm are generated between the pores of the sediment skeleton, and the hydrate saturation in sediments can be changed by setting the number of hydrate particles. The microscopic parameters of model particles and contact are based on the discrete element numerical simulation study of hydrate sediments [26, 27]. Due to the large discrete element model of borehole stability, to ensure the calculation speed of the model, the sand particles and hydrate particles are also amplified. However, this amplification process is directly related to the setting of the particle microscopic parameters, so the particle size changes will not cause changes in the mechanical response of the model.
2.3. Verification of Numerical and Analytical Models
For deep-water shallow hydrate reservoir, the in situ stresses are decided by the overburden stress. So two horizontal stresses can be set as equal. The model in Figure 1 can be modified as Figure 2. Since the analytical model is based on the assumption of horizontal isotropic stress, to compare and verify the calculation results of the numerical model and the analytical model, it is necessary to establish a cylindrical discrete element model of the hydrate reservoir (Figure 2) and apply the same horizontal and vertical ground stress as the analytical model so that the discrete element model and the analytical model have the same ground stress conditions. The model dimension adapted is large enough to eliminate the boundary effect. During calibration, the horizontal ground stress was 3 MPa, the vertical stress was 5 MPa, and the shaft wall support stress was 1 MPa.

In addition to ensuring that the calibrated models have the same ground stress conditions, it should also be ensured that the constitutive characteristics and mechanical responses of the hydrate formation are the same. The constitutive model in the form of internal friction angle described by equation (4) determines the constitutive characteristics and mechanical response of the reservoir in the analytical model. To fit the strain corresponding to the peak value tan, the expression of strain in equation (4) is transformed to obtain equation (12). By changing , , , , and D in formula (12), a combination of parameters with the same mechanical response as the discrete element model can be fitted. Figure 3 shows the relationship between tan and strain . Table 1 lists the parameter combinations of equation (12) obtained by fitting. To facilitate comparison, when the analytical model is used to study the influence of different saturations, the corresponding parameter combinations in Table 1 are also used. where is the initial internal friction angle, is the strain, and , , , and are the parameters related to the constitutive model.

Figure 4 compares the calculation results of the analytical model and the discrete element model, where is the deformed borehole radius in the analytical model and is the distance between a point and the borehole center. It can be seen from the figure that the stress distribution states around the well according to the analytical model and the discrete element model are very similar. Since the analytical model is based on the continuity assumption of elastoplastic mechanics, the calculation results of the analytical model are more continuous, while the discrete element model is composed of discrete element granular elements, so the measured stress values are more volatile, but the trend of the two changes is consistent.

3. Parametric Analysis of Wellbore Stability
In 2007, China’s NGH drilling project successfully drilled a physical sample of NGH in the Shenhu area in the northern South China Sea. The hydrate layer was 18-34 m thick, with a saturation of 20%-43%, and the methane content in the released gas was higher than 99% [28]. Several trial mining works led by the Guangzhou Marine Geological Survey in the Shenhu waters have proved that there are abundant hydrate resources in the Shenhu seabed [29, 30], and the Shenhu waters of the South China Sea have gradually become an important area for hydrate exploration and mining research in China.
With the advancement of hydrate exploration and research, hydrate reservoir mining methods and concepts are constantly being updated. A deep-sea shallow hydrate reservoir has the characteristics of large reserves and poor cementation and is a weakly consolidated or unconsolidated nondiagenetic hydrate reservoir. The exploitation of this type of hydrate reservoir can be accompanied by environmental pollution and geological disasters such as massive releases of methane gas and submarine landslides [31]. The Marine Hydrate Development and Research Team proposed for the first time a new method of “solid fluidization mining of submarine non-diagenesis hydrate deposits.” This had core advantages such as low pollution, low secondary disasters, and no damage to the lower porous reservoir hydrate [32]. The present paper assumes that the solid fluidization method is used to mine the deep-sea shallow hydrate, and it is necessary to study two types of wells in the hydrate formation drilling: drilling vertical wells and producing horizontal mining wells (see Figure 5) [33]. This paper takes a shallow hydrate reservoir in Shenhu sea area as the research object and adopts an analytical model and discrete element method to study drilling a vertical well and a mining horizontal well, respectively.

Chongyuan et al. studied the in situ stress state of a seabed formation in the northern South China Sea and found that the ratio of the maximum horizontal principal stress to the vertical stress in the northern South China Sea was about 0.76, and the ratio of the maximum to minimum horizontal principal stress ranged from 1.07 to 1.18 [34]. Sun et al. studied the downhole mining environment of the hydrate test well SH2 in the Shenhu sea area of the South China Sea and gave the pore pressure and mud pressure of the SH2 well at different depths [35]. Since the analytical model and discrete element model used in this article are drainage models, the actual stress conditions used are effective stresses. Combined with the abovementioned studies, Table 2 lists the ground stress, mud pressure, pore pressure, three-dimensional effective ground stress and borehole support stress (effective mud pressure), and other parameter values.
Since the borehole stability analytical model is based on the assumption of horizontal isotropic stress, when using this analytical model to study the wellbore stability of a vertical well, the horizontal stress is the maximum horizontal stress in the table. It can be seen from the table that the gap between the maximum level and the minimum ground stress is relatively small, so the use of this approximation method will not cause too much deviation in the analytical model calculation results, and the results still reflect the stress state of the hydrate formation around the well.
3.1. Borehole Stability Analysis in the Vertical Well
3.1.1. Effect of Depth
Figure 6 shows the relationship between borehole support stress and at different depths, where and are the initial borehole radius and the deformed borehole radius in the analytical model. The larger the ratio, the smaller the deformed borehole radius. In Figure 6, the 200 meters below sea floor (mbsf), 400 mbsf, and 600 mbsf borehole support stresses all decrease with increase of , and all have experienced a process of deceleration of decline. When is small, the strain near the borehole is small, and the stratum deforms elastically only. The formations of 200 mbsf, 400 mbsf, and 600 mbsf appear to have a plastic deformation at , 1.0217, and 1.037, respectively. At this time, the speed of the borehole support stress curve decreases significantly.

Figures 7–9 show the relationship between radial stress, circumferential stress, and vertical stress around the well and at different depths, respectively. In the analytical model, represents the distance between a point and the borehole center; the larger the ratio, the further away the point is from the borehole. Figures 7–9 show that the greater the depth, the greater the three-dimensional stress. The radial stress increases with increase of , and the circumferential stress and the vertical stress decrease with increase of , but this trend has nothing to do with the depth. In the figures, the back of the elastoplastic interface is the elastic deformation zone, the area between the elastoplastic interface and the strain softening-hardening interface is the strain hardening zone, and the strain softening zone is before the strain softening-hardening interface. As the depth increases, the ratios of the elastoplastic interface and the strain softening interface around the well gradually decrease. Table 3 lists the positions of the elastic-plastic interface and the strain softening-hardening interface corresponding to different depths.



3.1.2. Effect of Hydrate Saturation
Figure 10 shows the relationship between borehole support stress and under different hydrate saturations. When is small, only elastic deformation occurs in the formation around the well. When is greater than 1.0037, a plastic deformation zone appears around the well, and the stress curve of the borehole support turns and starts to slow down. Figure 10 shows that the hydrate saturation has little effect on the stress of the borehole support.

Figure 11 shows the distribution of radial stress, circumferential stress, and vertical stress around the well under different hydrate saturations. Figure 11 shows that the saturation change has a small effect on the three-dimensional stress state. Table 4 lists the positions of the elastic-plastic interface and the strain softening-hardening interface corresponding to different saturations. It can be seen from the table that the ratio of the elastic-plastic interface under the three saturations remains unchanged, but the ratio of the strain softening-hardening interface changes slightly, which increases with increasing saturation.

3.2. Effect of Mining Depths
Figures 12–14 show the relationships between borehole radial strain, borehole diameter reduction, and time at different initial borehole diameters at depths of 200 mbsf, 400 mbsf, and 600 mbsf. As the number of time increases, the radial strain of the borehole and the shrinkage of the borehole diameter gradually tend to remain unchanged, indicating that the stress state around the borehole reaches equilibrium after deformation and the borehole shape no longer changes.



Figure 15 shows the effect of depth on the maximum radial strain and hole diameter reduction. The maximum radial strain in the figure is the value of borehole radial strain tending to be constant in Figures 12–14. It can be seen that as the initial borehole diameter increases, the absolute reduction in borehole diameter also increases, but the radial strain of the borehole first increases and then decreases. The initial borehole diameter is in the range from 0.8 to 1.2 m, and the radial strain is large, indicating that the relative deformation of the borehole is small at this time. Therefore, the initial borehole diameter is in the range from 0.6 to 1.2 m, which can be used as a relatively optimal diameter for production boreholes.

4. Conclusion
This paper has used a discrete element method to establish a borehole stability model of a deep-water shallow hydrate reservoir. The discrete element numerical model has been compared with an elastic-plastic analytical model of borehole stability based on strain softening to verify the reliability of the numerical model. The analytical model and the discrete element model were used to study drilling vertical wells and production horizontal wells, respectively, and the influence of factors such as depth and saturation on the stress state and borehole strain around the well were analyzed. The results show that the near-wellbore zone can be divided into a plastic strain softening zone, a plastic strain hardening zone, and an elastic zone according to the different mechanical characteristics. Reservoir depth and hydrate saturation can change the stress state near the well. The greater the depth and the lower the hydrate saturation, the greater the borehole shrinkage. The diameter of the optimal horizontal well in the goaf is in the range from 0.6 to 1.2 m, the actual hydrate reservoir may have strong heterogeneity, and the optimal production hole diameter may be slightly smaller than this range. The calculated results should be verified by experiment. However, the experiment has not been carried out owing to difficultly constructing the borehole of hydrate in laboratory. The next step is to carry out a field experiment for the stability study of hydrate borehole.
Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge support from the National Natural Science Foundation of China (No. U19A2097), the National Key R&D Program of China (No. 2018YFC0310200), and the Key Laboratory Construction of Robot Engineering and Intelligent Manufacturing (Grant no. 19SXHZ0033).