Critical drawdown pressure prediction for sanding production of underground gas storage in a depleted reservoir in China

Sanding production of the underground gas storage (UGS) in the depleted reservoir will cause the decline of effective storage capacity, and the erosion of well string and auxiliary equipment, and other engineering problems. A deep understanding of the sanding production contributes to maintaining the injection and production capacity of gas storage and the long‐term safety of the UGS. In this paper, core samples of five main sedimentary microfacies were drilled from the UGS in middle China. The mechanical properties with different conditions of the moisture content in multi‐cycling injection–production process were tested. The critical drawdown pressure (CDP) calculation methods for the UGS were established based on the stress analysis near the wellbore and the Mogi–Coulomb criterion. Adopting the finite element method simulation on the geological model of the target UGS, the geo‐stress of the formation was calculated by inversed analysis and validated based on the Kaiser effect of rock. Then, the CDP of two target wells of the UGS were investigated. The effects of the water moisture and the cycling times of the injection–production process on the CDP were also analyzed. The geo‐stress distribution shows that the horizontal maximum principal stress ranges from 53 to 68 MPa, and the horizontal minimum principal stress is 46–58 MPa in the target reservoir section, and the vertical minimum principal stress is 69–77 MPa in the target reservoir. The calculated CDP changed with the types of the sedimentary microfacies in the range of [6.27, 8.77] MPa, in which the CDP of the mixed mud and sand flat sedimentary facies was the lowest. The weakness of the CDP with the increasing of the moisture content and the cycling times were analyzed quantitively.

6][7] UGS can be classified as: depleted oil and gas reservoir, aquifer, salt cavern, and abandoned pit.The UGS built in depleted oil and gas reservoir is the most dominant type and occupies more than 70% of the total number of UGS.In the multi-cycling injection-production process, the sanding production with the high-speed natural gas in the depleted oil and gas reservoir will erose the strings, valves, and other equipment.Many UGS in the depleted oil and gas reservoir have faced the engineering challenges of sanding production, such as Wen 96, Wen 23, Hutubi in China, 8 Redfield, Vincent, 9 Hillsboro 10 in America, Hajdúszoboszló, 11 Haidach 12 in Europe.Thus, a deep understanding of the sanding production in the depleted oil and gas reservoir contributes to maintaining the injection and production capacity of gas storage and the long-term safety of the UGS.
Sand production is a common but serious problem in traditional oil and gas fields, and of the weakly consolidated and unconsolidated sandstones.Many scientific efforts have been devoted to the sanding mechanism and sanding prediction modeling, through theoretical analysis, experiments and numerical simulations. 13,14The sand production in engineering applications is a complex process, influenced by various parameters, including elastic modulus, pressure gradient near the wellbore, mechanical strength of the reservoir rock mass, variations of the geo-stress, fluid pressure and flow rate, and so forth. 15,16From the perspective of perforation, the risk of sand production also depends on the direction of well deviation and orientation or perforation direction. 17Bradford et al. 18 established a semi-analytical elastoplastic wellbore stability model by assuming the isotropic and linear elastic behavior of rock and considering the effect of pore pressure.The constitutive model described the wellbore rock failure by the effects of the geo-stress and pore pressure under ideal conditions.On this basis, Stavropoulou et al. 19 assumed that both rock elasticity and strength (cohesion) depended on porosity and coupled porosity by sand erosion with wellbore rock failure.The material becomes weaker as porosity increases, and the reservoir produced sand when the rock starts to break.Willson et al. 20 proposed a prediction method for sand production rate by introducing load factor (intensity-normalized nearwell formation pressure) and Reynolds number (a function of permeability, viscosity, density and flow rate), and derived the empirical relationship among these parameters.However, only single sanding mechanism was considered in these studies, while the interaction of influencing factors such as perforation size, frequency, and direction, particle size, and shape, depth, pressure, permeability and fluid type, was neglected. 21any scholars have also conducted experiments on sanding production in the lab.Tronvoll 22 conducted large scale simulation experiments on sand production in the lab and found that fluid flow had little effect on cavity failure.They found the geo-stress and fluid flow rate were the key parameters affecting sand production, and the amount of sand production with increasing time.But the sand production rate of the weakly compacted sandstone has little relationship with the time factor. 23,24Ewy et al. 25 considered the effects of the geo-stress on the wellbore failure and predicted the CDP of the wellbore coupling the linear elasticity and linear poroelasticity theory with the modified Lade failure criterion in 2001.Some other scholars have considered the effect of water cuts on sand production.Bruno et al. 26 conducted triaxial compression tests on weakly cemented sandstones of different moisture content, and revealed sand production increases with the moisture content through the visualized sands pack model in the lab.Bianco et al. 27 analyzed the influence of the saturation of the wetting phase on the behavior, shape and stability of sand arches in two-phase flow experiments.It indicated that the bond strength of single-phase saturated sandy soil is not enough to support a stable arch, and cause a large amount of sand production.However, the increase of wetting phase saturation in a small range will not affect the stability of the arch, and the instability of the arch would occur after exceeding the critical value.
In addition, numerical simulations on the sanding predictions based on various assumptions, constitutive laws, failure criteria and numerical procedures have also been developed recently.Papamichos et al. 28 proposed a numerical model coupling the pore mechanics and erosion behavior to predict the external stress at the onset of sand production.The results showed a quadratic polynomial relationship between the sand production rate and the external stress and fluid flow rate, and the sand production rate is constant with time under constant external stress and flow rate.Nouri et al. 29 established a new numerical model based on bilinear strain-hardening/ softening Mohr-Coulomb model to simulate strain softening of the material and effective stress relaxation at the cavity surface.The model could be used to analyze borehole stability over time and predict the sanding CDP and rate in real time.Detournay 30 used FLAC to reproduce the slit pattern of cavity evolution observed in laboratory settings.They found the fracture mechanism was a combination of volume collapse (compaction zone formation) and hydrodynamic transport of failed material, and the sand production volume monotonically increased with simulation time.Volonté et al. 31 established a three-dimensional (3D) coupled model using finite element modeling to simulate fluid flow and rock deformation, and quantified the spatial distribution and damage severity of the rock around the perforation, adopting several optimal geometries and meshes describing the wellbore, perforation tunnel, and surrounding strata.Kim et al. 32 showed a 3D numerical model quantitatively describing the sand failure, erosion, and production.It indicated that sand production is firstly caused by the failure of wellbore or perforation, followed by the erosion process.The model achieved a balance between the mechanical friction for holding sand grains into the formation and hydrodynamic forces that cause separation of sand grains.Azadbakht et al. 33 combined finite difference method and finite element method (FEM) to study the rock fatigue and the sand bond degradation due to cyclic pressure variations, while capturing the effects of key parameters in sand production process.Cai et al. 34 established a dual-scale FEM-FEM model based on homogenization to simulate the compaction process of viscoplastic particles.They used micro scale FEM to model the 2D particle structure and macro scale finite element to analyze the homogenized model.The model solved problems caused by various types of granular media.Climent et al. 35 applied a coupled CFD-discrete element method (DEM) model to predict the sanding production under the drying and hydrostatic conditions.For the dry sandstones, a plastic zone was formed around the perforation where bond fracture occurred and the particle stress was reduced, and caused sand production.While in the hydrostatic conditions, drag force generated by fluid flow, slowed the movement of particles towards the cavity, and thus prevented sand production.Ghassemi et al. 36 coupled the DEM and the lattice Boltzmann method (LBM) to investigate the destabilization mechanism of the sand arch around the perforation and the sand erosion process.The established LBM-DEM coupled model was validated by laboratory experiments and was used to evaluate the effect of different parameters (e.g., flow rate, surrounding stress) on the velocity and volume of sands production.Wu et al. 37 simulated typical borehole progressive failure of high-porosity sandstone by gradually decreasing the internal support stress or increasing the far-field stress using FEM/DEM method.It was found that the failure form around the perforation may change from shear failure to mixed failure, and then to compaction failure with the rising of the average stress, while the rising of the cohesive strength led to the change from compaction failure to mixed failure or mixed failure to shear failure.There have been a large number of relevant studies on sand production in the oil and gas reservoir, including theoretical analysis, experiments and numerical simulations.After the long-term development, the free sand of the depleted oil and gas reservoir can be negligible.
However, sand production has been reported in the UGS of the depleted oil and gas reservoir, such as Wen 96, Wen 23, Hutubi in China, Redfield, Vincent, Hillsboro in America, Hajdúszoboszló, Haidach in Europe.9][40][41][42] These factors were not taken into account in sand production studies on the development of the reservoir, which could not be fully applied to UGS in depleted oil and gas reservoirs.Song et al. accurately predicted CDP based on inversion results of numerical simulation of the geo-stress in the target reservoir, combined with experiments and logging data.It showed that the depletion of reservoir pressure and the increase of water saturation will also lead to the decrease of CDP. 8 Thus, to guarantee the safe production, accurate prediction of the critical drawdown pressure (CDP) of UGS is necessary for predicting the sand production problems.
In this paper, the main sedimentary microfacies of the target formation in the UGS in China were classified, and the corresponding cores were drilled from two target wells.The mechanical properties of the core samples with the different conditions of the moisture content in multi-cycling injection-production process were tested in the laboratory, including the Young's modulus, Poisson's ratio, cohesive force, frictional angle, and so on.The CDP calculation methods for the UGS in the depleted gas reservoir were established based on the stress analysis near the wellbore district and the failure criterion.Adopting the FEM simulation on the geological model of the target UGS, the geo-stress of the formation was calculated by inversed analysis and validated based on the Kaiser effect of rock.Then, the CDP of two target wells of the UGS along the vertical direction were investigated, as well as the effects of the water moisture and the cycling times of the injection-production process on the CDP.UGS IN DEPLETED GAS RESERVOIR

| Stress analysis on wellbore of UGS
For UGS in depleted gas reservoir, the free sand and weakly consolidated sand of reservoir have been brought out with the gas extraction process and are basically ignored due to long-term extraction, as shown in Figure 1A,B.Thus, most of the sand production from UGS in depleted gas reservoir is caused by the failure of rock mass, which mostly starts at the well wall and gradually proceeds into the formation, as shown in Figure 1C.The rock failure at the well wall may cause cavity of the formation and weaken the strength of the rock mass, which would worsen the sands production.
Thus, the stress state of the rock mass near the well is one of the key factors for the sanding production of UGS.In this section, the stress state of the wellbore crosssection in the UGS was analyzed, and the following assumptions were adopted.
pressure at the bottom of the well is P w .Assuming the angle between the wellbore and the normal direction is θ, and the angle between the wellbore and the maximum principal stress is β.A cylindrical coordinate system (R, α, z) was established along the injection-extraction well, as shown in Figure 2A.
As shown in Figure 2B, a micro-element is taken from the rock around the well wall for the stress balance analysis.The stress field of this micro-element is , then the stress tensors in the converted coordinate system can be expressed as 42 : (1) The radial, tangential, and axial stress in the cylindrical coordinate system of the wellbore surrounding are as follows 42 : (2) (3) At this point, the shear stresses around the well are as follows 42 : F I G U R E 2 Force analysis, coordinate transformation and micro-element system decomposition of rocks around the wellbore of an underground gas storage in depleted gas reservoir.(A) Converted coordinate system, (B) micro-element force analysis.
where P w is the flow pressure in the well bottom, MPa; R is the radial distance, m; R w is the radius of wellbore, m; When R = R w , the geo-stress on the rock at the well wall is σ Combining Equations ( 2)-( 8), the three-way effective principal stresses in the rock of the well wall are as follows: The maximum principal stress σ 1 , the intermediate principal stress σ 2 , and the minimum principal stress σ 3 around the well are as follows: The Mogi-Coulomb criterion is based on the improved Mohr-Coulomb criterion by adding the effect of the three-way principal stress on the strength of the rock, which can be expressed as follows 43 : where τ otc is the octahedral shear stress, MPa; σ m,2 is the average value of normal stress, MPa; F is the failure index of matrix rock mass of well wall; A and B are the material constants related to the cohesion and the angle of internal friction; ε is the Biot coeffcient.When F < 0, the rock of well wall will fail leading to sand production; F > 0, the rock of well wall remains stable; thus, the CDP calculated at F = 0 is the CDP of the sand production of the UGS.

| MECHANICAL PROPERTIES TEST ON THE ROCK SAMPLES FROM THE TARGET UGS
The core samples for the mechanical test were drilled from two target wells (named #1 and #2).There are mainly five sedimentary microfacies in the target UGS of the depleted gas reservoir, including the sand flat (SF) sedimentary facies, bar sandstones (BAS) sedimentary facies, flood-channel (FC) sedimentary facies, mixed mud and sand flat (MMSF) sedimentary facies, and beach sandstone (BES) sedimentary facies.The mechanical properties of the core samples with the different conditions of the moisture content in multi-cycling injection-production process were tested in the laboratory using the MTS815 system (as shown in Figure 3), including the Young's modulus, Poisson's ratio, cohesive force, internal friction angle, and so on.

| Mechanical test of core samples of different sedimentary microfacies types
Fifty-seven sandstone core samples of five sedimentary microfacies of the target wells (#1 and #2) were drilled and used for laboratory tests, including weight, volume, density, porosity, permeability, and so forth.Triaxial compression tests were conducted under the confining pressure of 10 and 70 MPa, and the stress-strain curve of one sandstone sample is shown in Figure 4A.Meanwhile, the cohesive force (C) and the internal friction angle (φ) were obtained by fitting the confining pressure and maximum principal stress of the mechanical tests, as shown in Figure 4B.
The results showed that the triaxial compressive strength of the tested core samples ranged from 82. | 4293

| Mechanical properties of rock samples with different moisture content
The core samples of different sedimentary microfacies were vacuumed for 12 h and then saturated with water under the 10 MPa of confining pressure condition for 24 h by using the LY-BH-50 vacuum pressurized saturation device (as in Figure 5A).First, all the core samples were saturated to 100%, then, dried to the moisture of 10%, 30%, 50%, 70%, respectively, as shown in Figure 5.
The mechanical test of the rock samples for different moisture content under the confining pressure of 63 MPa were conducted.The results showed that the mechanical properties of the rock will be weakened after moisture content, and the strength of rock degraded more obviously with the increase of moisture content.Considering that there is no grain broken into the finer rock fragmentation, the internal friction angle φ is found to be almost constant for the same sedimentary microfacies with different moisture content.The fitted curves between the cohesive force and the moisture content were shown in Figure 6

| Mechanical properties of rock samples after multi-cycling loading
To simulate the multi-cycling injection-production process of UGS, the core samples from the target Well #1 were selected for cyclic loading tests with 50, 100, and 200 times of cycling, respectively, under 63 MPa of confining pressure condition.The core samples failed at the end of the cycle.The results showed significant degradation of the rock strength with the increase of the cycling times.Considering that there is no grain broken into the finer rock fragmentation, the internal friction angle φ is found to be almost constant for the same sedimentary microfacies after multi-cycling loading.The fitted relationship between the cohesive force and the number of cycles was shown in  The in situ mechanical properties of rock were estimated by coupling the laboratory test results and well-logging data including the compressional sonic travel time (DTC), shear sonic travel time (DTS), shale volume (VSH), density logging data (DEN) and the true formation resistivity (RT).The empirical formulas were adopted to estimate the Poisson Ratio, Young's modulus, cohesion, and internal friction angle using DTC, DTS, VSH, and DEN, as follows 44 : where μ d is the dynamic Poisson's ratio; E d is the dynamic Young's modulus, MPa; Δt p is the compressional sonic travel time, μs/ft; Δt s is the shear sonic travel time, μs/ft; ρ b is the density of rock, g/cm 3 ; V sh is the shale volume.The mechanical properties of different sedimentary microfacies in the target Well #1 and #2 were obtained by inversed analysis based on the well-logging data.The geo-stress field profiles of the target formation were established based on the geo-stress inversion results, as shown in Figure 9.Then, these mechanical properties of rock and the geo-stress for the target wells by inversed analysis were adopted for the CDP calculation by substituting the established CDP calculation Equations ( 1)-( 12), (17), and (18), as shown in Figure 10.For the target Well #  SF, and MMSF, respectively.The results showed that the CDP in the injection wellbore area is closely related to the strength of the rock sedimentary microfacies itself and the state of the geo-stress.Considering the similarity of the layers, the CDP of different sedimentary microfacies depends on the strength of the rock sedimentary microfacies.The strength of MMSF and FC is lower, so its CDP is correspondingly smaller.compaction of rock particles, leading to a reduction in rock strength, which in turn affects the CDP, especially for those sedimentary microfacies with high mud content, such as FC and MMSF.

| CDP result for multi-cycling loading
The curves of the CDP versus the depth under the conditions of 50 cycling times were presented in.increasing of cycling times.The reason is that the damage caused by alternating loads to rock with lower strength may be more serious.The weakening of rock mechanical properties due to rock fatigue damage is an important reason for sanding in the UGS reservoir, and the increase of cycling times will lead to the gradual accumulation of rock fatigue damage until rock failure and sanding occurs.

| CONCLUSIONS
In this paper, a total of five main sedimentary microfacies of the target formation in the UGS in China were classified, and the corresponding cores were drilled from the target Well #1 and Well #2.The mechanical properties of 57 core samples with the different conditions of the moisture content in multicycling injection-production process were tested in the laboratory, including Young's modulus, Poisson's ratio, cohesive force, frictional angle, and so on.
The CDP calculation methods for the UGS in the depleted gas reservoir were established based on the stress analysis near the wellbore district and the Mogi-Coulomb criterion.Adopting the FEM simulation on the geological model of the target UGS, the geo-stress of the formation was calculated by inversed analysis and validated based on the Kaiser effect of rock.Then, the CDP of two target wells of the UGS along the vertical direction were investigated, as well as the effects of the water moisture and the cycling times of the injection-production process on the CDP.
The following conclusions can be achieved: 1.The geo-stress distribution showed that the horizontal maximum principal stress ranges from 53 to 68 MPa, and the horizontal minimum principal stress is 46-58 MPa in the target reservoir section, and the vertical minimum principal stress is 69-77 MPa in the target reservoir section.2. The calculated average CDP of five sedimentary microfacies showed that: for the target Well #1, the BES is 8.77 MPa, the BAS is 8.85 MPa, the FC is 8.16 MPa, the SF is 8.54 MPa, and the MMSF is 6.79 MPa; for the target Well #2, the BES is 8.48 MPa, the BAS is 8.07 MPa, the FC is 7.7 MPa, the SF is 8.1 MPa, and the MMSF is 6.27 MPa. 3. The weakness of the CDP with the increasing of the moisture content and the cycling times were analyzed quantitively.The decreasing rates varied for different sedimentary microfacies with the increasing of the moisture content and cycling times.The results showed that the CDP of MMSF and FC was more sensitive to moisture content and cycling times.

F I G U R E 4
Triaxial compression test on the core samples of beach sandstone in Well #1.(A) Stress-strain curve of one core sample under 70 MPa confining pressure.(B) The fitted Mohr circle and the calculation of the internal friction angle and cohesion strength.
. The C of BES decreased from 16.59 to 12.89 MPa; the C of BAS decreased from 18.97 to 16.11 MPa, the C of FC decreased from 18.15 to 13.43 MPa, the C of SF decreased from 21.15 to 17.55 MPa, and the C of SF decreased from 14.24 to 11.51 MPa.

Figure 7 .
The C of BES decreased from 16.59 to 14.88 MPa; the C of BAS decreased from 18.97 to 17.45 MPa, the C of FC decreased from 18.15 to 16.87, the C of SF decreased from 21.15 to 19.96, and the C of MMSF decreased from 14.24 to 12.93 MPa.

F 1 |
I G U R E 5 LY-BH-50 vacuum pressurized saturation device and photos of core samples after saturation.F I G U R E 6 Cohesive force of five sedimentary microfacies for different moisture content.BAS, bar sandstone; BES, beach sandstone; FC, flood-channel; MMSF, mixed mud and sand flat; SF, sand flat.SANDING CDP 4.Analysis of mechanical properties of rock mass in injection-production wells

4. 2 | 4 . 3 |
Analysis of the geo-stress distributionThe spatial distribution of the geo-stress fields was calculated by adopting the FEM simulation on the geological model in target UGS.The simulation was conducted as a two-way coupled iteration adopting Petrel and FracMan software, and the Kozeny-Carman model was chosen for the calculation of the stress sensitive permeability model.The simulation results were validated by historical matching of the fluids production data in the operating process of both the gas reservoir and UGS.As shown in Figure8, the distribution of the geostress was heterogeneous.The horizontal maximum principal stress ranges from 53 to 68 MPa, and the horizontal minimum principal stress is 46-58 MPa in the target reservoir section, and the vertical minimum principal stress is 69-77 MPa in the target reservoir section.Adopting the Sample #2 from 2880.1 m as the target well, the horizontal two-way geo-stress are 55.5 and 47.0 MPa measured by the Kaiser effect of rock, while the simulation results are 55.6 and 49.1 MPa, and the prediction accuracy is larger than 95%.Calculation of CDP4.3.1 | CDP result of the target UGS wells

F I G U R E 7
Cohesive force of five sedimentary microfacies for different cycling times.BAS, bar sandstone; BES, beach sandstone; FC, flood-channel; MMSF, mixed mud and sand flat; SF, sand flat.

4. 3 . 2 |
CDP results for different moisture contentThe curves of the CDP versus the depth under the conditions of 10% moisture content was presented in Figure11.For the Well #1 when the moisture content is 10%, the average value of CDP decreases to 8.40, 8.28, 7.49, 7.99, and 5.82 MPa for BES, BAS, FC, SF, and MMSF, F I G U R E 8 Numerical simulation results of geo-stress distribution.(A) Horizontal maximum principal stress, (B) Horizontal minimum principal stress.F I G U R E 9 Well logging and geomechanical properties of the target wells.(A) Well #1.(B) Well #2.respectively.For the Well #2: the average value of CDP decreases to 8.09, 7.71, 7.13, 7.56, and 5.36 MPa for BES, BAS, FC, SF, and MMSF, respectively.The results indicated that, with the increasing of the moisture content, the CDP of different sedimentary microfacies showed different degrees of decreasing, which are significant for FC, SF, and MMSF, as shown in Figure12.The reason lies in that the water cut destroys the cohesion and F I G U R E 10 Critical drawdown pressure of different sedimentary microfacies in the target wells.(A) Well #1.(B) Well #2.F I G U R E 11 Critical drawdown pressure for 10% moisture content of the target wells.(A) Well #1.(B) Well #2.

Figure 13 .
Figure 13.For the Well #1, after 50 cycling times, the average value of CDP decreases to 8.52, 8.53, 7.91, 8.29, and 6.30 MPa for BES, BAS, FC, SF, and MMSF, respectively.For the Well #2: the average value of CDP decreases to 8.09, 7.61, 7.44, 7.81, and 5.87 MPa for BES, BAS, FC, SF, and MMSF, respectively.The results indicated that, with the increasing of cycling times, the CDP of different sedimentary microfacies showed different degrees of decreasing, as shown in Figure 14.The CDP of MMSF decreases more with material constants related to the cohesion and the angle of internal friction μ d dynamic Poisson's ratio E d dynamic Young's modulus (MPa) F failure index of matrix rock mass of well wall P p pore pressure (MPa) P w bottom hole flow pressure (MPa) R radial distance (m) R w radius of wellbore (m) S H maximum horizontal principal stress (MPa)F I G U R E 14 Weakening curve of the critical drawdown pressure for different cycling times.(A) Well #1.(B) Well #2.BAS, bar sandstone; BES, beach sandstone; FC, flood-channel; MMSF, mixed mud and sand flat; SF, sand flat.