Numerical Simulation on Sand Production Based on Laboratory Gas Hydrate Production Experiment

Abstract

Gas from natural gas hydrate (NGH) is priced competitively with gas prices. Most marine NGH is stored in low cementing strata, which easily cause sand production problems, restricting the commercial production and environmental safety of NGH’s development. This study applied a numerical simulation on sand production in hydrate-bearing sediments’ (HBS) exploitation. The numerical simulation on sand production was carried out for different productions of laboratory NGH exploitation. The results show radial strain appeared to be deformed away from the wellbore and show radial displacement close to the wellbore during mining. Due to the overburden stress condition, the boundary condition wall was a displace less rigid body. The radial displacement was greatly affected by depressurization, which showed the displacement to the wellbore and sanding. The radial strain was dominant by the shear shrinkage phenomenon in the mechanical model, while the reservoir’s radial displacement was away from the wellbore instead. The balance between the fluid driving force of production rates towards the wellbore and radial displacement drawing away from the wellbore is significant to sand production in HBS. The dominant forces of sanding were different mechanical and hydraulic combinations in three periods of GH production.

1. Introduction

Depressurization-driven gas production from coarse-dominant and fines-controlled sediments of HBS may cause sand production [1], which causes environmental pollution and production equipment damage [2]. During the laboratory sand production of HBS production, it is difficult to observe the internal process of sanding [3] and monitor the spatiotemporal evolution process [4]. Therefore, thermal-hydraulic-mechanical-chemical process (THMC) numerical simulation research can simulate the sand production process in the hydrate laboratory, coupling the temperature field, flow field, and stress field of the hydrate reservoir together to analyze the sand production mechanism [5,6]. Effectively extending the laboratory time–space scale sand production experiment to the on-site time–space scale helps the field mining for long-time production [7]. Therefore, some scholars have carried out numerical simulation research on sand production during natural gas hydrate mining.

TOUGH+HYDRATE+FLAC3D were applied to simulate hydrate reservoirs, and it was found that the shear failure of hydrate reservoirs during depressurization mining will promote formation settlement and is a great risk to sand production [3,8]. A thermal-hydraulic-mechanical sand production theoretical model (THMC) was established for hydrate formations. Researchers analyzed the migration process and location of sand during the hydrate trial production in Japan in 2013, and they concluded that the uneven distribution of reservoir stress will cause the shear deformation of the reservoir; that is, the uneven distribution of stress in the low saturation hydrate reservoir will lead to the shear deformation of the high hydrate saturation reservoir [5,9,10,11]. A TOUGH+HYDRATE model was established for the hydrate depressurization mining process. Parameters such as temperature, pressure, and mass were imported into FLAC3D to calculate stress and other parameters, and then the stress and other parameters were imported into the sand production model of PFC3D. The results show that an increase in bottom hole pressure will increase seabed subsidence and increase sand production, but it does not affect normal short-term test production. For long-term test production operations, the relationship between productivity, reservoir stability, and sand production need to be balanced. For the South China Sea low-permeability hydrate reservoir, under a high bottom hole pressure difference, most of the sand production was due to the formation of squeezing shear damage, and the sand created by gas and water migration accounted for a relatively small proportion [12]. For the enlightenment of sand production control in hydrate wells and the key issues that need to be solved, according to the 2017 test wells, the precision design of gravel packing and sand control was made, and the layered sand control was proposed [13,14]. Yan et al. found that increasing the pressure drop and the production pressure difference would lead to the concentration of rock stress around the well and lead to sand production. They believed that the initial sand production had nothing to do with the initial hydrate saturation, while the later sand production had an impact on the initial high hydrate saturation storage. The hydrate saturation of the layer was sensitive to pressure decrease [15]. Yu et al. proposed a sand control method for hydrating high argillaceous silt reservoirs through numerical simulation, evaluated the impact of different sand control methods on gas production efficiency and sand production, and proposed appropriate sand control design standards that could effectively improve the hydrate’s performance [16]. Wan et al. concluded that the greater the permeability and depressurization of the South China Sea hydrate reservoir, the greater the amount of subsidence and the faster the subsidence rate would be; they believed that subsidence mainly occurred in the early stages of mining, and that the effective stress increased [17]. It is the main cause of the subsidence of the reservoir. Jin et al. conducted a numerical simulation study on the exploitation of South China Sea hydrate through horizontal wells. It is believed that the settlement in the early stage of exploitation is very likely to account for more than half of the total settlement; when considering the uneven vertical distribution of reservoir hydrate, the reservoir’s lower permeability of the upper and lower strata contributes little to the production of hydrate extraction, but it makes the subsidence worse [18]. Zhang et al. found that sand production in gas hydrate reservoirs is mainly caused by the dissociation of gas hydrate, leading to the fallout of the cemented fines or the eroding of host rock by the fluid’s flow [19].

However, the expensive cost of field tests and the black box of hydrate experiments points to the numerical simulation of sand production for hydrate exploitation as the best option. To study the sand production mechanism during hydrate exploitation, a sand production numerical simulation was carried out in this paper.

2. Sand Production Model in NGH Exploitation

2.1. Establishment

In this paper, the sand production simulation was carried out jointly using TOUGH+HYDRATE and Abaqus (Figure 1). The coupling method is similar to the coupling between TOUGH+HYDRATE and FLAC. In this case, the coupling was unidirectional. The pressure and hydrate saturation were calculated by TOUGH+HYDRATE, which were imported into the commercial mechanics software Abaqus. In ABAQUS, the mechanical response to pore pressure changes was calculated, and mechanical parameters such as stress and sediment deposition were obtained. The Fortran program was applied to program the mechanical properties of the HBS during exploitation, which improved the accuracy of the stress–strain model. Because of the unidirectional coupling, the effect of porosity changes on seepage after Abaqus calculations was not considered.

2.2. Basic Assumptions

The No. 5 test from previous sand production tests was simulated in this paper [20]. The experiment was simplified to a two-dimensional axisymmetric model, where the simulated area was the sediment portion of the sediment reactor. Because of the axisymmetry, the simulated area was an rz profile (Figure 2). The whole model was set up according to the experimental conditions. The model’s geometry was as follows: r, 0.16 mm~76 mm; z, 0~100 mm; mesh profile, 13 meshes equally divided vertically and 21 meshes divided radially, 21 × 13 = 273 meshes in total.

The initial pore pressure of the model was 11 MPa, and the sides and top were watertight boundaries where no material exchange could occur with the outside boundaries (

Table 1. Parameters used in the model.

Elastic Modulus (SH = 50%)	450 MPa	Porosity Ratio	0.667
Elastic modulus (SH = 0)	100 MPa	Permeability	10 × 10−13 mD
Poisson’s ratio	0.15	Hydrate saturation	46%
Dry density	1.5 g/cm3

). The left side was a variable pore pressure boundary condition, where the pore pressure values were obtained by the interpolation of the measured pore pressure. The left side of the model corresponded to the location of the production well, and the production well pressure was reduced in two stages: 0–2.5 min from 11 MPa to 3.4 MPa; the second stage of production pressure was reduced slowly from 3.4 MPa to 0.1 MPa at atm. The top side of the model had a stress boundary of 12 MPa, and the right side and bottom side were set to a zero displacement boundary due to the constraints of the reactor wall. The mechanical boundary conditions on the left side were simplified to a free displacement boundary. The initial pore ratio in the model was 0.67 based on the porosity conversion, the permeability was 10 × 10⁻¹³ mD, and the hydrate saturation was 46%. The initial moment of the model simulation corresponded to the moment when the sediment was consolidated at 12 MPa at the top and ready to open the valve at the bottom to start mining.

2.3. Verification of Discrete Element Model

After the model was built and meshed and the boundary conditions were set, the stress equilibrium of the system was first established, and the stress equilibrium obtained in this model is shown in Figure 3. In the initial state of the model, the vertical displacement at the top was larger than that at the bottom, and the vertical displacement obtained was in the order of 10⁻⁶ mm (Figure 3a), which was much smaller than the experimental settlement of 10⁻³ mm.

The radial displacements of the model were on the order of 10⁻⁶ mm (Figure 3b). The strains after stress equilibration were extremely small. The radial strains (Figure 3c) and axial strains (Figure 3d) were loaded by the overlying stresses, which varied more regularly and both showed an order of 10⁻⁷ mm. The pore pressure was 11 MPa (Figure 3e), which was consistent with the experimental tests. The vertical and radial displacements of the model met the ABAQUS’s requirements for initial stress balance, and this state could be used as the initial state for hydrate extraction.

3. Results and Discussion

3.1. Numerical Simulation Results of Sand Production in the Process of Gas Hydrate Exploitation

As shown in Figure 4, the pore pressure P of the numerical simulation was consistent with the pore pressure variation pattern during the experimental mining process, which was consistent with the two main stages of the model assumptions. The axial displacement of the numerical simulation was related to the pressure. However, the subsidence caused by changes in permeability and fluid properties during the extraction process could not be characterized yet, which affected the fitting results. The axial displacement U2 at the beginning of the pressure drop was larger, but the gradually slow subsidence during the extraction process could not be characterized. Finally, the total subsidence was still similar. Therefore, the effect of settlement during gas hydrate exploitation was influenced not only by changes in stress but also by other factors that need to be explored. As shown in Figure 5 and Figure 6, the stress variation, pore pressure, radial strain, and axial strain during the buckling process were simulated based on the intrinsic model of the gas hydrate exploitation process. The total reduction in radial displacement was considered to be the volume loss of the reservoir, i.e., the volume of sand released. Both the axial strain and axial displacement were expressed as the downward movement by overburden stress. It is noteworthy that the radial strain was deformed to the right during the extraction process, away from the wellbore, while the radial displacement was displaced to the left during the extraction process, close to the wellbore. We assumed that under the overburden stress, the right wall was a non-displaced rigid body and that the radial strain was manifested as the shear contraction phenomenon in the South China Sea gas hydrate mechanics model, but the radial displacement of the reservoir was characterized as displacement towards the wellbore.

3.2. Spatial Distribution of Pressure in the Process of Gas Hydrate Exploitation

Figure 7 shows the spatial distribution of pressure during the extraction process. The original formation pressure showed a more uniform axial distribution, while the pressure changed to a more uniform radial distribution during the opening minutes of the well dewatering. The pore pressure decreased from an initial 11 MPa to 8.64 MPa, 6.04 MPa, and 3.56 MPa at 50 s, 100 s, and 150 s, respectively. The pressure drop was large, but the radial pressure difference within the reservoir was 0.2 MPa. Thus, during the first stage of production, the pressure near the wellbore was low, and the pressure away from the wellbore was high, showing a radially uniform pressure gradient. The distribution of the pressure gradient was radially uniform. However, as production progressed, the pressure dropped to 2.945 MPa, 2.4342 MPa, 1.4067 MPa, and 0.3789 MPa at 50 min, 100 min, 200 min, and 300 min, respectively, and the pressure difference within the reservoir decreased to 0.12 MPa, showing a non-uniform pressure gradient distribution. In the radial direction, the pressure was higher away from the wellbore; in the axial direction, the pressure gradient distribution showed the influence by the overlying stress. Thus, in the second and third stages of production, the pressure distribution in the reservoir was mainly influenced by the gas pressure, while the pressure differential was mainly due to the different permeability of the reservoir. During these two stages of production, the combined effect of overburden stress loading and the gravity of the overlying sand, the lower sediment of the reservoir was gradually compacted, resulting in lower sediment porosity and the permeability being lower than that of the upper sediment of the reservoir, further resulting in higher pressure in the lower reservoir than the upper. On this basis, the radial depressurization gradient was gradually deflected by compaction, showing the lowest pressure in the upper reservoir near the borehole wall and the highest pressure in the lower reservoir away from the borehole wall. No further change in the pressure gradient was observed during the continuous pressure drop.

In summary, the pressure gradient in the first stage of production was strongly influenced by the dewatering pressure drop, while in the gas production stage, the pressure gradient was influenced not only by the gas production pressure drop gradient but also by the overlying stress loading, resulting in a deflection of the pressure gradient.

3.3. Spatial Distribution of Radial Strain in the Process of Natural Gas Hydrate Exploitation

Figure 8 shows the spatial distribution of radial strains during the extraction process. The radial strains in the original formation showed a relatively uniform distribution, with a strain of −2 × 10⁻¹⁰. The whole reservoir formation showed increasing strain away from the wellbore, and formation near the wellbore showed increasing strain to the right. The strain near the wellbore at 50 s, 100 s, and 150 s varied from 0.003, 0.006, and 0.0085, respectively, away from the wellbore to 0.0006, 0.001, and 0.0015, respectively, at the boundary, and the strain difference between the wellbore and the boundary increased from 0.0024, 0.005, and 0.007, respectively, indicating that the strain was gradually concentrated in the interior of the reservoir. Thus, during the first stage of production, the rapid depressurization due to drainage and the shearing of the reservoir due to the pressure drop resulted in a more pronounced change in the radial strain. It is worth noting that the smaller points of radial strain in the initial formation had an effect on the radial strain distribution during the first stage of mining, where the radial strain at the upper part of the boundary was smaller, resulting in a deflection of the radial strain in the upper part, which gradually became smaller as mining progressed. As mining progressed, the radial strains at 50 min, 100 min, 200 min, and 300 min were 0.0045, 0.004, 0.004, and 0.0028, respectively, showing a more uniform strain gradient distribution with less overall variation. Therefore, during the second and third stages of production, the radial strain in the reservoir was low during the process of gas production and pressure reduction.

In summary, the radial strains showed a shearing phenomenon, with the radial strains in the first stage of production being more influenced by the drainage depressurization and initial strain, while the radial strains in the reservoir varied less during the gas production stage.

3.4. Spatial Distribution of Axial Strain during Natural Gas Hydrate Exploitation

Figure 9 shows the spatial distribution of the axial strain during the extraction process. The axial strain in the original formation showed a relatively uniform distribution with a strain of −8 × 10⁻⁸, and the whole axial strain showed a downward strain under the overlying stress. In the first few minutes of well drainage and depressurization, the axial strain was subjected to the force field change by depressurization. The axial strain showed a deflection towards the near-borehole wall. While this deflection was the largest in the upper part of the near-borehole wall and the smallest in the upper part of the boundary, the whole axial strain showed an increasing downward strain. Thus, during the first stage of production, rapid depressurization was caused by the drainage, and the axial strain of the reservoir was deflected due to the shear force generated by the pressure drop. The deflection of the axial strain was gradient. As mining progressed, axial strains were −0.0154, −0.0163, −0.01829, and −0.02021 at 50 min, 100 min, 200 min, and 300 min, respectively, with small overall variations. Although the variation in axial strain was small, the strain gradient was very unevenly distributed, with the smallest axial strain at the top of the boundary deposit and the largest axial strain near the top and bottom of the wellbore. This was mainly due to the overburden stress at the bottom and the compaction of the upper sediments, resulting in a larger axial strain at the bottom. The top of the near-wellbore was subject to greater axial strain variation due to the combined effects of suction from the pressure drop and shear due to pressure drop and overburden stress. The top of the boundary, conversely, was less affected by the pressure drop, although it was affected by the overburden stress. Therefore, during the second and third stages of production, the axial strain in the reservoir was smaller but not uniformly distributed during the gas production and pressure drop.

In summary, the axial strain exhibited a downward strain, and although the overall axial strain was small, its axial strain gradient was unevenly distributed.

3.5. Spatial Distribution of Radial Displacement in the Process of Natural Gas Hydrate Exploitation

Figure 10 shows the spatial distribution of the radial displacement during the extraction process. The radial displacement of the original formation showed a relatively uniform distribution, with a strain of −0.07 mm. The radial displacement near the wellbore at 50 s, 100 s, and 150 s was −0.05 mm, −0.09 mm, and −0.15 mm away from the wellbore, respectively, while the boundary was a gradual increase in radial displacement from −0.01 mm to −0.03 mm. Therefore, during the first stage of production, the rapid depressurization caused by drainage resulted in a significant change in radial displacement and deflection of −0.15 mm (more than 2%) of the sand body near the borehole wall, with some of the sand body possibly entering the wellbore and exiting the sand.

As mining progressed, the radial displacements near the borehole wall were −0.15 mm, −0.013 mm, −0.011 mm, and −0.07 mm at 50 min, 100 min, 200 min, and 300 min, respectively, and the radial displacement at the boundary gradually increased by −0.01 mm, showing a more uniform stress gradient distribution. Therefore, during the second and third stages of production, the radial displacement of the reservoir was sheared by the radial strain during the process of gas production and decompression, showing a significant reduction of radial displacement.

In summary, the radial displacement showed the trends towards the wellbore increased and then decreased. In the first stage of production, the radial displacement was influenced by the drainage and depressurization, presenting a significant displacement towards the borehole wall and sanding, while in the second production stage with gas production, the reservoir radial displacement was influenced by the radial strain shearing, and its displacement decreased. Therefore, in the first stage of production, the radial displacement was mainly influenced by the combined effect of fluid and overburden stress, which may lead to partial sand production. In the second and third stages of production, the radial shear of the overburden stress was influenced by the shear of the sand body, and some sand was still produced despite partial shrinkage.

3.6. Spatial Distribution of Axial Displacement in the Process of Gas Hydrate Exploitation

Figure 11 shows the spatial distribution of the axial displacement of the pressure during the extraction process. The axial displacement of the original formation showed a relatively uniform distribution, with a displacement of −5 × 10⁻⁶ mm. The whole axial displacement showed a downward displacement under the overlying stress. At 50 s, 100 s, and 150 s, the top axial displacement was −0.4 mm, −0.9 mm, and −0.1 mm, respectively, with the axial displacement differential between the layers and the axial displacement near the borehole wall being slightly larger. Thus, during the first phase of production, the axial displacement of the reservoir resulted in the settlement of the reservoir due to pressure drop. As production progressed, the top axial displacement settled slowly at −1.5 mm, −1.6 mm, −1.8 mm, and −1.9 mm at 50 min, 100 min, 200 min, and 300 min, respectively, with a little overall variation. The axial displacement gradient was evenly distributed, but the difference between layers was large and increasing, with the largest axial displacement at the top of the sediment and the smallest axial displacement at the bottom. This was mainly due to the overburden stress on the bottom sediment and the compaction of the upper sediment, resulting in a smaller axial displacement of the bottom sediment. The top sediment, conversely, was subject to overburden stress and sand exit gap, which further compacted the top sediment, resulting in a larger axial displacement. Therefore, in the second and third stages of production, the axial displacement of the reservoir decreased steadily and slowly during the process of gas production and pressure reduction.

In summary, the axial displacement showed a downward displacement, and although the entire axial displacement gradient was relatively uniform, the axial displacement difference was large.

4. Discussion

Most gas hydrate resources are stored in unconsolidated or weakly consolidated sediments [21]. Sand production is an unavoidable problem during hydrate exploitation, compared with that of the rock reservoirs of gas and oil wells. The volumetric deformation of hydrate-bearing sediments and the driving force of fluid during hydrate exploitation are the key factors of sand production [22,23,24,25,26,27], while the volumetric deformation and sand production of HBS is also a key problem of CO₂ leakage from CO₂ hydrate sequestration [28,29].

From previous studies on the mechanical response of HBS [30,31,32,33,34,35,36,37], the typical sample deformation and stress-strain curves of HBS in tri-axial are shown in Figure 12 and Figure 13. The shear shrinkage (Figure 12c) and dilatancy (Figure 12d) are two typical deformations in the rock and soil tri-axial tests. Most gas and oil wells show shear dilatancy in rock; the mechanical failure of rock means sand production, which is widely applied to principles of sanding in conventional petroleum wells. However, the unconsolidated or weakly consolidated HBS during depressurization shows shear shrinkage, which could show the different sanding behaviour and sand prediction model. Because of the weak strength of HBS, the mechanical failure model for sanding is difficult for determining the sand production of HBS.

The sand production behaviour of HBS is different in the three GH production periods [20,38,39,40]. For the wellbore in the gas hydrate reservoir (Figure 14a), the HBS showed the shear dilatancy and moved into the wellbore (sand production) in the first period of GH production by gas/water fluid (Figure 14b). In the first period, the sanding of HBS was contributed by the HBS dilatancy force (mechanical field) and the gas fluid driving force (hydraulic field) together. In the second period of GH production by high flow rates of gas fluid carrying water (Figure 14c), the HBS showed the shear shrinkage and moved away from the wellbore, meaning that the HBS shrinkage drawing force (mechanical field) resisted the gas fluid driving force (hydraulic field). In the third period of GH production by the low-flow rate gas fluid but high-pressure differential (high effective stress), the HBS showed the shear dilatancy and moved into the wellbore again (Figure 14d), while fluid driving force may have been low due to mud-cake formation by the HBS filtration properties and wall building [40,41]. Above, the sanding behaviours were dominated by different mechanical (M)–hydraulic (H) components, which were the M+H, M-H, and M in the first, second, and third periods of GH production, respectively.

5. Conclusions

The sand production model of gas production from HBS was built, and the sand production mechanism of gas production from HBS was discussed.

(1) The pressure gradient was dominated by drainage and depressurization; meanwhile, the radial strain was controlled by drainage step-down and the initial strain variables effect in the first stage of production. In the second production stage, the pressure gradient was deflected by the depressurization gradient and overburden stress, where the radial strain of the reservoir changed little.

(2) The axial strain presented a downward strain trend. However, the entire axial strain was small, and the axial strain gradient was the non-uniform arrangement.

(3) The radial displacement towards the borehole showed a trend of increase–decrease. The radial displacement was dominated by the drainage step-down, which was expressed as obvious reservoir displacement towards the borehole and sanding in the first stage of production. The radial displacement of the reservoir was affected by the radial strain shear shrinkage, and thus reservoir displacement decreased instead in the second production stage.

(4) The axial displacement showed a trend of downward displacement. However, the entire axial displacement gradient was uniform, and the axial displacement difference was large.

(5) Different mechanical (M)–hydraulic (H) components dominated key factors of sand production in the first, second, and third periods of GH exploitation.