Abstract
The dissociation of natural gas hydrate in the seabed during hydrate production can lead to substantial weakening in soil mechanic properties, which in turn results in the deformation and even destruction of the seafloor sediments, thereby posing a potential threat to the safety of the undersea engineering infrastructures. A better understanding of the effect of soil mechanical weakening on seabed deformation will benefit the prediction and/or determination of geological conditions for safer hydrate production. In this work, a three-dimensional model of the submarine geological body was established and engineering properties of hydrate reservoir under different degrees of hydrate dissociation were determined. Such model and engineering parameters were applied for simulation calculation on examining the effect of hydrate dissociation degree on soil subsidence using FLAC3D software. The simulation results show that the subsidence of sediments is centrally distributed within the dissociation area and the subsidence of hydrate-bearing sediment’s top surface is much more severe than that of the overlying layer. Subsidence on the horizontal section above the hydrate reservoir presents an essentially symmetric dish-like subsidence trough, with the greatest subsidence in the centre and approximately constant subsidence beyond the trough range. The radius of the subsidence trough increases with the increase of the degree of dissociation, and the increasing extent in the early dissociation stage (lower dissociation degree) is greater than in the middle and the later stage (higher dissociation degree). As the degree of hydrate dissociation increases, subsidence of both the overlying layer and the HBS increases. Especially when the decomposition exceeds a certain threshold, the seafloor subsidence increases rapidly. The simulation results in this work are of guiding significance for deformation monitoring and infrastructure deformation prevention in practical hydrate production projects.
1. Introduction
Natural gas hydrate is one of the most promising energy sources characterized by widespread distribution, huge reserves, high energy efficiency, and little pollution, thereby bearing a broad development prospect [1–3]. The hydrate is typically distributed in submarine continental shelves, slopes, and permafrost regions. The sediment deformation is to be inevitably occurred due to the dissociation of gas hydrate during production, which is very likely to induce geological hazards such as seabed landslide, posing a potential threat to the stability of submarine engineering infrastructure and the marine ecological environment [4–7]. Therefore, the study on the prediction and understanding of the seabed deformation behavior is of practical significance for the development of hydrate mining.
In this regard, many scholars have performed numerical simulation studies to investigate the effect of hydrate dissociation on the subsidence deformation of seafloor sediments. Kimoto et al. [8] reported the behavior of multiphase materials by using the theory of porous media, and a numerical simulation method was developed to analyze the dissociation process of gas hydrate in hydrate-bearing sediment (HBS). They showed that sediment deformation is concentrated around the dissociation area and the cause of deformation is the generation and seepage of water and gas. Wang et al. [9] used the finite element software ABAQUS to simulate the effect of hydrate dissociation on the deformation and stress distribution of pipelines and seafloor sediments in the case that the gas pipelines are buried in the HBS. They found that the displacement and stress of the pipelines and the adjacent soil layer are closely related to the hydrate dissociation range. In the hydrate production process, different hydrate production modes will directly affect the subsidence deformation of the seabed. By using a coupling model built with FLAC3D, Gong et al. [10] studied the effect of gas hydrate dissociation in the South China Sea on the subsidence of seabed soil under different production modes. They reported that the prolonged heating duration on subsidence scope results in more obvious deformation than heating temperature. Kwon et al. [11] established a two-dimensional model using the finite difference method to investigate how HBS loses its stability in the case where the production well transfers heat to the surrounding hydrate deposits. Their work suggested that excess pore water pressure is generated during the dissociation of the gas hydrate, causing the plastic deformation of deposit in the dissociation zone and the hunch at the seafloor. Lu et al. [12] carried out numerical simulation analysis on the influence of control factors on the seabed deformation and the stress distribution during natural gas dissociation and proposed a simplified critical condition to evaluate the stability of gently inclined stratum upon natural gas hydrate dissociation. Considering the influence of stress on the elastic constant, Zhao and Shang [13] constructed the E-B constitutive model of Duncan Zhang and applied it to the finite element software ABAQUS to simulate the deformation of sediments; the results showed that the deformation of seabed increases nonlinearly as dissociation radius of hydrate increases. Jin et al. [14] used TOUGH and Biot theory to study the effect of production pressure and formation permeability on submarine soil subsidence under depressurization production mode in horizontal wells.
Note that most of the above studies focused on the effects of different factors such as production methods, formation dip angles, and hydrate dissociation range on seabed deformation. In addition, the weakening of soil properties caused by hydrate decomposition is also one of the key factors leading to seabed sediment deformation. To date, various methods of extracting hydrates were reported and represent the state of the art of hydrate mining technologies. These methods include depressurization, thermal stimulation, inhibitor injection, and carbon dioxide replacement [15, 16]. Yet, such production methods are mainly based on the dissociation of natural gas from its original crystalline, icy substance to gaseous phase. As such, the fine equilibrium state between the pristine hydrates and the surrounding soils will be broken during hydrate production [17, 18] because the solid hydrates act as soil skeleton and cementation to maintain the strength of HBS. The dissociation of natural gas hydrates derived from production in the seabed can lead to the weakened physical and mechanical properties of the soil, incurring excessive deformation of overlying layers [19]. Although a nondissociation-based hydrate extraction method called underground solid fluidization was recently reported [16], which involves hinged suction crushing and depressurization of hydrates, it remains in its proof-of-concept stage as only a few cubic meter gas hydrates were produced by this method.
Besides, hydrate production is essentially a dynamic process. As hydrate production proceeds, the dissociation degree of the hydrate reservoir increases continuously. Thus, the engineering properties of hydrate reservoir gradually become weakened and the deformation of the soil layer and its distribution characteristics will also change accordingly. Hence, based on determining the engineering properties under different dissociation degrees, the degree of hydrate dissociation that represents the transient production status can be used to assess and predict the influences of various production stages on the subsidence during hydrate production. As such, this work used variant physical and mechanical parameters of soil under different production stages to study the subsidence evolution of seabed. The regularity and spatial distribution of seabed subsidence are investigated via simulation using the finite difference method, and the effects of the degree of hydrate decomposition on the subsidence distribution characteristics are examined. The results could enrich the knowledge of the seabed deformation during hydrate production and help develop protective measures for the engineering facilities.
2. Numerical Simulation
FLAC3D is a piece of numerical simulation software based on the Lagrangian algorithm and is applicable for the simulation of plastic failure and large-scale deformation [20]. It is adopted to simulate the deformation of seabed in this work.
2.1. Geological Model and Numerical Model
The existing drilling data of the Shenhu sea area in the South China Sea show that the main hydrate layers are typically distributed in the depth of 100∼200 m below the seabed, and their thickness is generally 10∼80 m [21–25]. Taking the stratigraphic conditions of W19 drilling in the Shenhu sea area as a reference, a geological model was established with the soil layer divided into three parts, from top to bottom: the overlying layer was set at 150 m, the HBS, 20 m, and the lower layer, 100 m thickness, respectively. The slope of the seafloor soil layer is assumed to be horizontal, and the depth of water is set as 800 m. A production well is located at the centre of the geological model. A sectional sketch of the geological model is shown in Figure 1.
Considering the practical production conditions of the hydrate production well, the hydrate dissociation range is defined as a cylindrical column with a radius of 20 m around the centroid production well, which is similar to the previous study [26]. According to the Saint-Venant principle, to decrease the boundary effects, the length and width of the numerical model are selected to be 400 m, which is up to 20 times the hydrate dissociation radius. To check whether the boundary constraints of such scale have an impact on the simulation results, a larger horizontal model is established. Through the comparative analysis of the simulation results of the larger horizontal model and the model used in this paper, the boundary constraints imposed in this paper have no significant impact on the simulation results. According to the spherical shape of the hydrate dissociation range, the zone within 20 m around the production well is subdivided more finely to improve the accuracy of numerical calculation with a focus on the soil around the production well.
The production well simulated by the unit pile has an inner diameter of 20 cm and a wall thickness of 2 cm, which is drilled to the top surface of the lower seabed layer. The soil and well pipe have different stiffness, which accordingly result in deformation incompatibility between the well pipe and soil during the actual mining process. Therefore, a contact surface is set between the outer surface of the well and sediments. The boundary constraints of the model are as follows: the top surface of the model is free, the bottom boundary is fixed in the z-direction, and the surrounding boundary is set as the normal displacement constraint boundary, that is, the face perpendicular to the X axis is the X-direction constraint and the face perpendicular to the Y axis is the Y-direction constraint. Because of the large disturbance around the mining well, the mesh around the mining well is subdivided properly to improve the accuracy of numerical simulation. The numerical simulation mesh of each mass is shown in Figure 2.
2.2. Calculation Hypothesis and Simulation Process
The dissociation of hydrates can cause the weakening of soil layers, which in turn leads to deformation of the rock and soil mass above and/or below the hydrate layer. During the extraction of hydrates, the dissociation occurs progressively. This causes the dynamic weakening of mechanical properties of the HBS layer in different production stages, which further results in the seafloor subsidence and its distribution of seafloor sediments changing with time. To explain the effect of hydrate dissociation degree of the subsidence, the method of controlling variables is utilized. The hydrate dissociation is set under four hydrate dissociation extents within the range from 30% to 100%. The calculated deformations of seafloor sediments, in tandem with that under pristine, undisturbed properties (before production), depict an essentially full evolution of the subsidence during hydrate production.
In the numerical simulation, it is assumed that there is a perfect seepage channel in the soil allowing no accumulation of excess pore water pressure and no volume change of the hydrate-bearing layer when gas and liquid seepage during hydrate dissociation. Besides, the overlying layer is assumed to be saturated throughout the production. Also, the production well pressure decreases and the HBS pressure variation is not considered.
Under the pristine status (without disturbing by the mining), there exists a stress field in the soil to balance the self-weight load and hydrostatic pressure. To obtain initial stress field in the model, simulation calculation on the established model was first performed under self-weight load and hydrostatic pressure. Then, the speed and displacement of soil were set to zero, while retaining the calculated initial stress distribution in the model. Then, the model material parameters are modified to weakened parameters under specific hydrate dissociation degree and calculations are conducted until the equilibrium is achieved. The calculations are repeated in the order of under 30%, 60%, 80%, and 100% hydrate dissociation degree, and the resultant subsidence data are acquired corresponding to each dissociation degree.
2.3. Constitutive Model and Parameters of Soil and Structure
Engineering studies show that the Mohr–Coulomb constitutive relation fits well with the mechanical behavior of rock and soil mass [27], and elastic constitutive relation is suitable for artificial material’s mechanics behavior [20] because the force on it is often lower than the ultimate strength. Therefore, Mohr–Coulomb constitutive relation is employed in soil elements, elastic constitutive relation is used in well elements, and column shear constitutive relation is used in contact elements.
According to the geotechnical investigation data in mechanical properties of hydrate sediments in the South China Sea, the solid hydrate is formed in the voids of sediments. Hydrate acts as not only cementation that makes soil particles closely linked together but also the soil skeleton to bear load directly. These effects reinforce mechanical properties of deposits such as elastic modulus, cohesive force, and internal friction angle. When hydrates are decomposed, solid hydrates are translated into the water and natural gas, leading to weakened cementation and destructed soil skeleton. As such, the elastic modulus and cohesive force are significantly reduced, whereas the internal friction angle is slightly reduced. In this paper, the internal friction angle assumed decreases monotonically with the increase of hydrate dissociation degree. Meanwhile, as gas escapes, water fails to invade the unsaturated space left completely, which leads to a slight reduction in soil density. Hence, by synthesizing reported data from the literature [26, 28–31], physical and mechanical parameters of soil layers under different degrees of hydrate dissociation are selected for simulation, which are tabulated in Table 1. Referring to the parameters of similar materials reported in the literature [20], physical and mechanical parameters of well material are shown in Table 2 and the parameters of contact surface are given in Table 3.
3. Results and Discussion
3.1. Spatial Distribution Characteristics of Vertical Displacement
To clarify the spatial distribution characteristics of the vertical displacement field of overlying soil and HBS upon hydrate dissociation, we employ the simulation results under 100% degree of dissociation conditions for analysis. The selected positions in the simulated cubic model and corresponding subsidence results are shown in Figure 3.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 3(a) illustratively gives the selected locations of concerned surfaces/lines in the model cubic that present the typical subsidence simulation results. Figure 3(b) is the subsidence contour of a vertical section across the centroid production well. Obvious subsidence can be observed in the upper part of the HBS and the overlying layer, whereas obvious rebound displacement can be observed at the HBS bottom surface and the underlying layer. The biggest subsidence occurred at the dissociation centre of the HBS top surface. The subsidence decreased significantly from the HBS top surface upward to the overlying layer and from the well axis centre to the surroundings. Likewise, the biggest rebound displacement was located at the HBS bottom surface, and the rebound displacement reduced greatly from the HBS bottom surface down to the lower layer.
Figure 3(c) also presents a similar subsidence configuration of the overlying layer top surface (seafloor): an essentially symmetric dish-like subsidence trough, i.e., the effect sphere with relatively significant subsidence change on the seafloor, is concentrated in a circular area, with the greatest subsidence in the centre and nearly constant subsidence beyond the subsidence trough.
Comparing the subsidence from Figures 3(c) and 3(d), one can rationally infer that the subsidence increases as the sediment buried depth increases. The HBS top surface exhibits the greatest centroid subsidence of 95.15 cm (Figure 3(d)), 4.75% of HBS thickness, while the centre subsidence of the top surface is rather small, around 9.06 cm (Figure 3(c)), just 0.45% of the thickness of HBS (20 m). Beyond the subsidence trough, vertical deformation of the overlying layer and the HBS becomes gradually stable at about 4.72 cm (subsidence) and 0.48 cm (rebound), respectively.
The rebound displacement occurred at the bottom of HBS (Figure 3(e)) is also similar to the subsidence distribution: the centre of the HBS bottom surface shows the maximum rebound displacement, and it decreases significantly as the stress change transmits either beyond the centre or downward to the lower layer (Figure 3(b)).
For a detailed analysis of the subsidence trend, Figure 3(f) gives the curves of subsidence with respect to the distance from the production well in the cubic model surface and HBS top surface. It can be seen that the subsidence curve for either overlying layer or HBS exhibits similar symmetric parabolic shape with maximum subsidence at the production centre and a level-off value to both sides, between which is the concave subsidence trough. The opening of the concave subsidence trough indicates the range where hydrate dissociation has a relatively predominant influence on the subsidence. The radius of such subsidence trough of the top surface is around 125 m, 6.25 times as big as the hydrate production range, whilst that of the HBS top surface is around 30 m, 1.5 times of the hydrate production range, only distributed in the upper range of decomposition area, indicating that the influence range of the subsidence of the upper layer is far greater than that of HBS top surface; the influence area shrinks as the depth increases.
Generally, Figure 3 shows the weakened mechanical properties of HBS because the dissociation of HBS can result in the stress change and subsidence of the overlying layer. The largest subsidence occurred at the central dissociation area. When the stress change transmits upwards, it gradually attenuates, but the lateral dispersion range will expand, leading to the decrease of subsidence amount and expansion of subsidence trough range in the overlying layer. Because of the decrease of the density caused by the decomposition of HBS, the underlying layer is affected by the unloading stress; the bottom of HBS is deformed upward, and the maximum rebound is also located in the centre of dissociation. The rebound is rapidly reduced when it transmits downward, and the rebound only occurs in a very small range. The spatial distribution of settlements follows the law of stress diffusion, indicating the simulation results are reasonable.
The above calculated spatial distribution subsidence suggested specific safeguard requirement for production infrastructure is needed for sustainable deformation durability. As the relative subsidence and subsidence range is quite significant either in HBS or at the production well, the adaptability of the production well and other facilities located in these areas to deformation should be sufficiently considered. Similarly, the construction materials and/or structure of the seabed infrastructure located in the subsidence trough should be elaborately designed with high tolerance to deformation to enhance structural stability.
3.2. Effect of the Hydrated Dissociation Degree on Subsidence
The subsidence of the overlying layer top surface and the HBS top surface under different degrees of hydrate dissociation is presented in Figure 4. Likewise, the subsidence curves for either the overlying layer or HBS exhibit subsidence trough with the greatest value centred in the well and decreased gradually to both sides. Also, the subsidence curve of the HBS top surface is much narrower with a confined opening compared with the overlying layer top surface, indicating the influence range of subsidence for HBS is much smaller than that for the top overlying layer.
(a)
(b)
Figure 4 also shows, for all the degrees of hydrate dissociation, the subsidence of the HBS layer far exceeds that of the overlying layer; the maximum subsidence of the HBS is 6∼10 times larger than the overlying layer. For example, the subsidence for the HBS layer is at the magnitude of tens of centimetres (8.43 cm for 30% dissociation degree and 95.15 cm for 100% dissociation degree), while for the overlying layer at few centimetres (0.71 cm for 30% dissociation degree to 9.07 cm for 100% dissociation degree). In contrast, the influence ranges in HBS (the radius of the concave trough of 22 m for 30% dissociation degree to 28 m for 100% dissociation degree) are rather narrower than those in the overlaying layer (76 m for 30% dissociation degree to 125 m for 100% dissociation degree). Such results indicate that the seafloor subsidence is concentrated in the lower HBS and diverged to a more wide area in the upper overlying layer when the deformation stress transmits upwards from the well.
For more detailed analysis, we extracted the maximum subsidence values and the subsidence trough radius of the top surface of overlying layer from Figure 4(a) and replotted as Figure 5. It can be seen that both the maximum subsidence and the radius of influence range increase monotonically with the increase of the degree of dissociation.
(a)
(b)
In the initial stage of hydrate extraction, e.g., when the dissociation degree is below 30%, the HBS dissociation and softening effect is relatively weak. Coupled with the self-stability of overlying strata and the diffusion of corresponding force, very small subsidence at the top surface of the overlying layer is induced. The maximum subsidence is only 1.26 cm, and the radius of the subsidence trough is 76 m, about 3.8 times the production area radius, which is 20 m.
As the production progressively proceeds, the dissociation degree increases accordingly to a medial range (between 30% and 60%), the HBS dissociation and softening effect is thus intensified, and the subsidence at the top surface of the overlying layer is also increasing. The maximum subsidence grew to ∼3 cm, and the radius of subsidence trough increased to 100 m, 5 times the production area radius.
When the dissociation degree increased from 60% to 80%, the increase in the range of the subsidence trough is slowed down; the radius of subsidence trough increased to 109.5 m. However, the subsidence at the top surface of the overlying layer above the production area increased more rapidly, reaching 7.06 cm.
In the later stage of production, near the full dissociation degree (between 80% and 100%), the increase of subsidence at the overlying layer top surface slowed down obviously. The final maximum subsidence was 9.07 cm, with a 125 m subsidence trough radius.
The above results implicitly indicated that, in the early stage of production, the scope of the subsidence trough will increase rapidly and the amount of subsidence at the overlying layer top surface is relatively small; thus, the geological body and infrastructure on the seafloor will not be that much greatly affected by the deformation.
However, the subsidence intensity and the gradient of subsidence change with distance in the dissociation range and its adjacent upper part are relatively large, and the deformation of the production well and shaft near the top surface of the dissociation range still needs to be monitored. With the development of production to a high extent (60%–80% dissociation degree in this work), the increasing range of subsidence trough will slow down, whereas the amount of subsidence will obviously increase, and uneven subsidence may occur on the top surface of overlying layer. Also, at this stage, the geological body or production infrastructure which is sensitive to deformation on the overlying layer should be included in the monitoring range and corresponding reinforcement measures should be fulfilled in routine production operations to avoid engineering accidents.
4. Evaluation
There are many factors causing seabed deformation due to gas hydrate production, including the weakening of mechanical properties, the change of gas-liquid seepage field, the change of temperature field, and the sand production. This paper focuses on the influence of mechanical properties’ weakening caused by hydrate decomposition on the deformation under different degrees of dissociation.
The above results indicate that the subsidence of sediments mainly concentrates around the dissociation range and the maximum subsidence is located at the dissociation centre, which is essentially consistent with Kimoto et al.’s [8] and Zhao and Shang’s [13] conclusions. Also, the deformation caused by properties’ weakening may reach the maximum of tens of centimetres on the top of HBS at the maximum. Such deformation value is rather significant and must be taken into consideration for safety in practical natural gas production. It is also concluded that there may exist a critical gas hydrated dissociation degree (80% in this paper possibly due to HBS’s large plastic deformation) for the sediment subsidence, beyond which a rapid sediment subsidence would occur.
Considering that the modulus of elasticity is not replaced by the modulus of resilience, the rebound amount may be relatively large. Because of the limitations of numerical simulation software and acquired data, gas-liquid seepage in sediments and change in temperature field are not considered in this work. If these factors are taken into consideration, more comprehensive results may be achieved in better agreement with practical conditions.
5. Conclusions
This paper focused on the simulation of the sediment subsidence caused by the weakening of HBS properties during hydrate production. A piece of finite element difference software, FLAC3D, was used to calculate the subsidence on an established geological model, which comprises hydrate-bearing sediment (HBS), top overlaying layer, and drilling well, under different hydrate dissociation extents that represent progressive mining stages from the undisturbed pristine (no dissociation) to completed status (100% dissociation). The simulation results show that the hydrate dissociation gave rise to significant subsidence in sediments. Such subsidence is centrally concentrated within the dissociation area; and, beyond the dissociation area, the subsidence decreases rapidly and spreads to a wider area, exhibiting a symmetrical relationship of subsidence with respect to the distance from the central well with the parabolic subsidence trough around the drilling well. The subsidence of the top surface of the HBS layer is far more severe than that of the overlying layer. The radius of the subsidence trough (the influenced area by HBS mining) increases with the increase of the dissociation degree, and the increasing speed is greater in the early stage (lower hydrate dissociation degree) than in the middle and the later stage (higher hydrate dissociation degree). In the later stage of production, the radius of subsidence trough can reach more than 6 times that of the production area. As the degree of hydrate dissociation increases, subsidence of both the overlying layer and the HBS increases, especially when the decomposition exceeds a critical degree (60%–80% for the model in this paper), the subsidence increases more rapidly.
The above evolution of subsidence indicates that, in the initial mining stage, the deformation and stability of the infrastructure and deformation-sensitive geological body in the expanding subsidence trough should be carefully monitored. While during the whole mining stage, the critical dissociation degree of accelerated development of deformation rate should be determined by the subsidence development trend to take corresponding treatment measures in time to alleviate the impact of deformation on mining facilities.
Data Availability
The data used to support the findings of this study are included within the article.
Disclosure
Part of the calculation results (the subsidence data under 100% HBS dissociation) was used to analyze the influence of the thickness ratio of overlaying layer to HBS on sediment settlement. Such work was presented as a conference paper in the 2nd International Conference on Advances in Civil Engineering, Energy Resources and Environment Engineering, May 22–24, 2020, Nanning, China, in IOP Conference Series: Earth and Environmental Science. This paper can be found in the following link: https://iopscience.iop.org/article/10.1088/1755-1315/526/1/012135.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under contract No. 41672309 and the National Key Technologies R&D Program of China under contract No. 2017YFC0307700.