# Ground Subsidence Prediction Model and Parameter Analysis for Underground Gas Storage in Horizontal Salt Caverns

Authors: Guimin Zhang; Yuxuan Liu; Tao Wang; Hao Zhang; Zhenshuo Wang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Due to a great demand of natural gas or oil storage in these years, horizontal caverns were proposed to fully use bedded salt formations of China. Under the same geological and operating conditions, the horizontal cavern would shrink more than traditional pear-shaped cavern, which might bring larger ground subsidence and affect the safety of storage facilities. A new prediction model was proposed in this paper for the time-dependent ground subsidence above horizontal caverns. The proposed model considered the impurity of bedded salt formations and simplified the horizontal cavern to an ideal cylinder. The shape of the subsidence trough was determined by the probabilistic integration method, and corresponding calculation formulas for the tilt, curvature, horizontal displacement, and horizontal strain were derived. Based on the assumption that the subsidence volume at the ground was proportional to the reduced volume of horizontal cavern, a formula for the reduced volume over time was established. FLAC3D was introduced to simulate the ground subsidence, and the results show that the proposed prediction model agreed well with the simulation results. Finally, the proposed prediction model was used to analyze the impacts of different stratigraphic parameters and design parameters. The results mainly show that, as the draw angle increases, the subsidence trough becomes deeper and narrower; as the depth of the cavern increases, the maximum subsidence first increases and then decreases, and the subsidence trough gradually becomes round; with the increase of the purity, the subsidence gradually decreases; with the increase of the creep properties and the stress exponential constant, the maximum subsidence first increases rapidly and then slowly approaches the limit; increasing the brine extraction velocity can shorten the cavern construction period and then reduce excessive ground subsidence; the subsidence decreases nonlinearly with the increase of internal pressure; with the increase of the cross section diameter and length of the horizontal cavern, the subsidence presents a significant nonlinear increase. In addition, unlike the traditional pear-shaped cavern, under the same conditions, the ground subsidence above the horizontal cavern according to this newly proposed model is much larger, and the ground subsidence contour line is no longer a standard circle. The findings of this study can help for better understanding of the prediction of ground subsidence above salt caverns and also provide a reference for the design and construction. However, the proposed prediction method is ideal and theoretical and should be further improved by engineering practice in the future.


---

## Body

## 1. Introduction

In China, most of the salt formations were deposited bedded with various interlayers between salt rocks [1, 2]. Comparing bedded salt formations with salt domes, bedded salt formations always have thinner single salt layer, and the grade or purity of salt rock is much lower [1, 2]. Due to a great demand of natural gas or oil storage in these years [1–4], horizontal caverns were proposed to fully use bedded salt formations of China (Figure 1(a)) [5–8]. However, compared with traditional pear-shaped cavern (Figure 1(b)), the shape of horizontal cavern is approximately horizontal cylinder, which is thought to be less stable [7, 8]. Under the same geological and operating conditions, the horizontal cavern would shrink more than pear-shaped cavern [9, 10]. The large volume shrinkage might bring ground subsidence and even lead to the surface collapse [11–16], which seriously affects the safety of the underground storage facilities.Figure 1
Schematic diagram of the salt caverns inferred from SONAR exploration [5]. (a) Horizontal cavern. (b) Pear-shaped cavern.
(a)(b)As a matter of fact, the ground subsidence resulted from the usage of salt formations has been researched by many scholars. Although field observation of the ground subsidence is the basis for establishing an acceptable prediction system, there are still various difficulties. In particular, field observation needs a long time and is a slow process [10, 17–24], and the required financial and material resources are rather large. Therefore, methods such as theoretical research and numerical analysis have been widely used to predict ground subsidence and guide engineering practice, and thus, considerable progress for avoiding or preventing surface subsidence have been made, including the reuse of the abandoned salt cavern for underground energy storage or filling salt cavern with alkali wastes [25, 26].Among various numerical analysis methods, the SALT_SUBSID was developed for the first time to summarize the subsidence characteristics of different mining stages and storage periods [24]. Ehgartner [27] used the finite element method to analyze the ground subsidence of the salt caverns for the petroleum reserve of the United States under different cavern spacings and internal pressures. Hoffman [28] used ABAQUS to study the impact of the buried depth of the storage cavern on the reduced volume and ground subsidence. Neal [29] used the finite element method to predict the ground subsidence and subsidence velocity caused by the cavern shrinkage in salt domes. Li et al. [30] studied the ground subsidence above underground salt cavern under different internal pressures and overlying strata by numerical simulations. The numerical models have outstanding advantage in the prediction accuracy, but usually require more time to build and calculate [31]. Therefore, it is still necessary to further research the theoretical analysis models.In theoretical research studies, Fokker and Visser [22] successfully predicted the ground subsidence caused by salt cavern shrinkage in the Netherlands by a theoretical model. In their model, surface movement matches the influence functions derived from geomechanics, and these influence functions may change over time. Compared with the ground subsidence caused by the point source volcanic eruption, Li et al. [13] tried to predict the ground subsidence induced by gas storage in the salt caverns, but their model did not consider the viscous plasticity of rock salt. Kong et al. [32] assumed that the amount of surface subsidence was linear to the reduced volume of the pear-shaped cavern and used a Gaussian function to fit the subsidence curve, but the reduced volume of the salt cavern should be computed by numerical simulation in advance. A dynamic subsidence model for gas storage in vertical cylinder salt cavern during the whole life cycle was put forward by Wang et al. [16], but their model did not include the impurities in the salt rock, and the vertical cylinder cannot represent the horizontal caverns. Considering that the mechanical properties of bedded salt rocks are different from the almost pure rock salt in the salt domes, a time-dependent prediction model was proposed by Zhang et al. [11] for the ground subsidence of the underground gas storage cavern in impure salt formations. However, the proposed model by Zhang et al. [11] is only suitable for pear-shaped caverns (Figure 1(b)), but not completely suitable for horizontal caverns (Figure 1(a)).This paper proposed a new prediction model for the time-dependent ground subsidence above the horizontal cavern in impure bedded salt formations, considering both the construction period and the storage period. Firstly, the horizontal cavern was simplified to be an ideal cylinder. Then, the shape of the subsidence trough was determined by the probabilistic integration method, and corresponding calculation formulas for the tilt, curvature, horizontal displacement, and horizontal strain were derived. Thirdly, by supposing that the subsidence volume at the ground was proportional to the reduced volume of horizontal cavern, a formula for the reduced volume over time was established. Fourthly, FLAC3D was introduced to verify the proposed prediction model. Finally, the influence of some factors on the ground subsidence was carried out with this model, such as the creep time, the draw angle, the subsidence adjusting coefficient, the cavern depth, the purity, the creep properties of the rock salt, the brine extraction rate, the internal pressure, the cross-sectional diameter, and the length of the horizontal cavern. The work of this paper is helpful to establish the prediction theory of ground subsidence above horizontal salt caverns and also provides a reference for the design and construction.


## 2. Prediction Model of Ground Subsidence above a Horizontal Cavern

### 2.1. Ground Subsidence, Tilt, Curvature, Horizontal Displacement, and Strain

As to the ground subsidence of underground gas storage in salt caverns, the root cause is the volume shrinkage of caverns [9–16]. Due to the creep effect of rock salt [32, 33], the reduced volume of salt caverns changes over time. Assuming that the relationship between the volume of subsidence trough and the reduced volume of salt caverns is linear, the ground subsidence of any location over time could be calculated as following:(1)Wx,y,t=a⋅wx,y⋅VstV∗,where x and y represent the position of a surface point, m, t represents time, month, a is the subsidence adjusting coefficient, 0 ≤ a ≤ 1, which is determined by the parameters of the overlying strata and should be generally calculated on the field monitoring data of subsidence, wx,y indicates the shape function of the ground subsidence induced by a unit reduced volume, m/m3, V∗ indicates the unit volume, 1 m3, and Vs(t) represents the reduced volume of the horizontal cavern, m3.Once the time-dependent ground subsidenceW (x, y, t) is obtained, the surface tilt in φ direction can be calculated:(2)ix,y,t,φ=∂Wx,y,t∂φ=∂Wx,y,t∂xcosφ+∂Wx,y,t∂ysinφ.And, the surface curvature inφ direction can be obtained as follows:(3)Kx,y,t,φ=∂ix,y,t,φ∂φ=∂ix,y,t,φ∂xcosφ+∂ix,y,t,φ∂ysinφ=∂2Wx,y,t∂2xcos2φ+∂2Wx,y,t∂x∂ysin2φ+∂2Wx,y,t∂2ysin2φ.Then, the horizontal displacement inφ direction can be calculated as follows:(4)Ux,y,t,φ=bx+byRx+Ry4ix,y,t,φ=bx+byRx+Ry4∂Wx,y,t∂xcosφ+∂Wx,y,t∂ysinφ,where bx and by are the horizontal movement coefficient, which needs to be inferred from the field observation data, Rx represents the influence distance along the length of subsidence trough, m, and Ry represents the influence distance along the width of subsidence trough, m.In addition, the horizontal strain inφ direction can be calculated as follows:(5)ξx,y,t,φ=bx+byRx+Ry4Kx,y,t,φ=bx+byRx+Ry4∂2Wx,y,t∂2xcos2φ+∂2Wx,y,t∂x∂ysin2φ+∂2Wx,y,t∂2ysin2φ.Above all, the ground subsidence, tilt, curvature, horizontal displacement, and horizontal strain can be calculated according to equations (1)–(5). However, the shape function wx,y and the reduced volume of the horizontal cavern Vs(t) should be obtained first.


### 2.2. Shape Functionwx,y

For a deep buried horizontal cavern (500∼2000 m), the cross-sectional diameter (no more than 60 m) is much shorter, and then, it can be assumed that the shape functionwx,y is not affected by the size of horizontal cavern. Therefore, probability integral method can be introduced to predict the ground subsidence, which was originally justified by Knothe [33] and then widely used all around the world [11, 34–38]. According to the basic principle of the probability integral method, it can be assumed that the probabilities of the ground subsidence along the length and width of the subsidence trough do not affect each other.In the case of sufficient mining, the draw angleβ indicates the angle between the line of the critical deformation value along the main section of the subsidence trough and the boundary point of the goaf and the horizontal line on the outside of the goaf, as shown in Figure 2. As the length L and the width D are the two dimensions of the horizontal salt cavern along the main sections, if the draw angle β is known, the ground influence range along the main sections and the shape function can be calculated as follows [11, 32]:(6)wx,y=1RxRye−πx/Rx2+y/Ry2,Rx=Hcotβ+L2,Ry=Hcotβ+D2,where β represents the draw angle, rad, which is determined by the parameters of the overlying strata and should be generally calculated on the field monitoring data of subsidence, D represents the cross-sectional diameter of horizontal cavern, m, L represents the length of horizontal cavern, m, and H indicates the depth of the cavern center, m.Figure 2
Shape simplification of the horizontal cavern. (a) Along the length (x axis). (b) Along the width (y axis).
(a)(b)Specially, along the length direction of the subsidence trough, the shape function can be justified as follows:(7)wx,0=1RxRye−πx/Rx2.Also, along the width of subsidence trough, the shape function is as follows:(8)w0,y=1RxRye−πy/Ry2.


### 2.3. Reduced VolumeVs(t)

Because rock salt has significant creep characteristics, if the brine pressure or natural gas pressure in the horizonal cavern is smaller than the original lithostatic stress, the cavern volume would shrink over time [39, 40]. In fact, the internal pressure cannot be exactly equal to the lithostatic pressure at the cavern depth [9, 11, 14, 16, 39–41]. This indicates that the convergence of the salt cavern will not stop.Because there are a number of insoluble interlayers between the salt layers and the flow of brine is really difficult to exactly control, the actual shapes of the horizontal caverns (Figure1(a)) are still irregular. In addition, the insoluble sediment might accumulate in a large amount at the cavern bottom. With solution mining by double-well convection, the height of horizontal cavern is usually lower than 60 m, which is much smaller than the buried depth (500 to 2000 m). Thus, the horizontal cavern can be assumed to be a cylinder (Figure 2). Since the underground gas storage has the construction period and operation period in a life cycle, the law of reduced volume in these two periods will be discussed below separately.


#### 2.3.1. Reduced Volume of the Horizontal Cavern in the Construction Period

In order to avoid the fluctuations of the manual operation, the brine extraction rate is usually kept constant in the construction period. If the salt layer is uniform and the solution rate of the cavern is constant, a series of equations for the reduced volume of the horizontal cavern can be calculated as follows [11]:(9)Vt1=AV•EBk11−e−Bk1t1,Vst1=ABVE+AB−1AV•EBk11−e−Bk1t1,A=ρb−ρwα⋅ρs−ρb−ρw,B=α⋅ρsα⋅ρs−ρb−ρw,t1=tn−t0,where Vt1 represents the horizontal cavern volume at time t1, m3, Vst1 represents the reduced cavern volume at time t1, m3, k1 represents the shrinkage rate of the unit volume in the construction period, m3/m3⋅d, VE represents the extracted brine volume, m3, V˙E indicates the brine extraction rate, m3/d, α indicates the proportion of the soluble components in salt rocks, namely, the purity, %, ρw indicates the water density, 1.0 × 103 kg/m3, ρb indicates the brine density, kg/m3, ρs indicates the rock salt density, kg/m3, t0 indicates the start time of the construction period, month, tn represents the current time, month, and t1 represents the characteristic time, month.In particular, both the cavern volumeVt1and the reduced cavern volumeVst1 are equal to zero when the characteristic time t1=0. Assuming that the time tp is when the solution mining is completed, the cavern volume and reduced cavern volume can be written as(10)Vtp=AV•EBk11−e−Bk1tp−t0,Vstp=ABVE+AB−1AV•EBk11−e−Bk1tp−t0,where Vtp indicates the cavern volume at time tp, m3, and Vstp represents the reduced cavern volume at time tp, m3.


#### 2.3.2. Reduced Volume of the Horizontal Cavern in the Operating Period

Assuming that natural gas is immediately stored in the horizontal cavern after the construction, the relationship between the cavern volume and the reduced cavern volume and the characteristic timet2 during operation can also be obtained:(11)Vt2=Vtpe−k2t2,Vst2=Vstp+Vtp1−e−k2t2,t2=tn−tp,where Vt2 indicates the cavern volume at time t2, m3, Vst2 represents the reduced cavern volume at time t2, m3, and k2 represents the shrinkage rate of the unit volume in the operating period of the horizontal cavern, m3/m3⋅d.


#### 2.3.3. Shrinkage Rate of the Unit Volumek1 and k2

If only the steady-state stage is considered in the creep deformation of bedded salt formations, the most widely applied creep constitutive model for rock salt is Norton power exponent model [9, 10]:(12)ε•=Cα⋅e−Q/RmTσdiffσ∗n,for uniaxial behavior,ε=•=Cα⋅e−Q/RmT1n+1∂3J2n+1∂σ=,for3D behaviour,J2=12SijSji,Sij=σij−13σiiδij,where C(α) is a material property of the rock salt, 1/month, which is assumed to be related to the grade α, as shown in equation (13), n indicates the stress exponent constant, σ∗ represents the unit stress, 1 MPa, σdiff represents differential stress, MPa, Q represents the activation energy, kJ/mol, Rm is the gas constant, 8.314 × 10−3 kJ/(mol/K), T represents the absolute temperature of the strata at the cavern center, K, which is related to the depth H, as shown in equation (13), J2 is the second invariant of the deviator stress tensor, MPa2, Sij and Sji represent the deviatoric stress tensor, MPa, σij indicates the stress tensor, MPa, σii represents the sum of diagonal elements of the stress tensor, MPa, and δij is the Kronecker delta symbol.Then, the shrinkage rate per unit volume of the cylindrical cavern can be inferred from the 3D formula of the Norton-power exponent model [9–11, 16]:(13)k=33P0−Pn⋅σ∗n⋅Cα⋅e−Q/RmT,Cα=Cs⋅α+Ci⋅1−α,T=273.15+20+0.02⋅H,where k stands for k1andk2and indicates the shrinkage rate of unit volume of a cylindrical salt cavern, m3/m3, P represents the brine pressure P1 during the building period or the gas pressure P2 during the operation period, MPa, P0 represents the original in situ pressure, MPa, Cs and Ci represent the creep properties of the rock salts and interlayers, 1/month, and H is the depth of the cavern center, m.Finally, substituting equation (6) and equations (9)–(13) into equation (1), the value of ground subsidence can be finally obtained, and then, the tilt, curvature, horizontal displacement, and horizontal strain can also be calculated.


## 2.1. Ground Subsidence, Tilt, Curvature, Horizontal Displacement, and Strain

As to the ground subsidence of underground gas storage in salt caverns, the root cause is the volume shrinkage of caverns [9–16]. Due to the creep effect of rock salt [32, 33], the reduced volume of salt caverns changes over time. Assuming that the relationship between the volume of subsidence trough and the reduced volume of salt caverns is linear, the ground subsidence of any location over time could be calculated as following:(1)Wx,y,t=a⋅wx,y⋅VstV∗,where x and y represent the position of a surface point, m, t represents time, month, a is the subsidence adjusting coefficient, 0 ≤ a ≤ 1, which is determined by the parameters of the overlying strata and should be generally calculated on the field monitoring data of subsidence, wx,y indicates the shape function of the ground subsidence induced by a unit reduced volume, m/m3, V∗ indicates the unit volume, 1 m3, and Vs(t) represents the reduced volume of the horizontal cavern, m3.Once the time-dependent ground subsidenceW (x, y, t) is obtained, the surface tilt in φ direction can be calculated:(2)ix,y,t,φ=∂Wx,y,t∂φ=∂Wx,y,t∂xcosφ+∂Wx,y,t∂ysinφ.And, the surface curvature inφ direction can be obtained as follows:(3)Kx,y,t,φ=∂ix,y,t,φ∂φ=∂ix,y,t,φ∂xcosφ+∂ix,y,t,φ∂ysinφ=∂2Wx,y,t∂2xcos2φ+∂2Wx,y,t∂x∂ysin2φ+∂2Wx,y,t∂2ysin2φ.Then, the horizontal displacement inφ direction can be calculated as follows:(4)Ux,y,t,φ=bx+byRx+Ry4ix,y,t,φ=bx+byRx+Ry4∂Wx,y,t∂xcosφ+∂Wx,y,t∂ysinφ,where bx and by are the horizontal movement coefficient, which needs to be inferred from the field observation data, Rx represents the influence distance along the length of subsidence trough, m, and Ry represents the influence distance along the width of subsidence trough, m.In addition, the horizontal strain inφ direction can be calculated as follows:(5)ξx,y,t,φ=bx+byRx+Ry4Kx,y,t,φ=bx+byRx+Ry4∂2Wx,y,t∂2xcos2φ+∂2Wx,y,t∂x∂ysin2φ+∂2Wx,y,t∂2ysin2φ.Above all, the ground subsidence, tilt, curvature, horizontal displacement, and horizontal strain can be calculated according to equations (1)–(5). However, the shape function wx,y and the reduced volume of the horizontal cavern Vs(t) should be obtained first.


## 2.2. Shape Functionwx,y

For a deep buried horizontal cavern (500∼2000 m), the cross-sectional diameter (no more than 60 m) is much shorter, and then, it can be assumed that the shape functionwx,y is not affected by the size of horizontal cavern. Therefore, probability integral method can be introduced to predict the ground subsidence, which was originally justified by Knothe [33] and then widely used all around the world [11, 34–38]. According to the basic principle of the probability integral method, it can be assumed that the probabilities of the ground subsidence along the length and width of the subsidence trough do not affect each other.In the case of sufficient mining, the draw angleβ indicates the angle between the line of the critical deformation value along the main section of the subsidence trough and the boundary point of the goaf and the horizontal line on the outside of the goaf, as shown in Figure 2. As the length L and the width D are the two dimensions of the horizontal salt cavern along the main sections, if the draw angle β is known, the ground influence range along the main sections and the shape function can be calculated as follows [11, 32]:(6)wx,y=1RxRye−πx/Rx2+y/Ry2,Rx=Hcotβ+L2,Ry=Hcotβ+D2,where β represents the draw angle, rad, which is determined by the parameters of the overlying strata and should be generally calculated on the field monitoring data of subsidence, D represents the cross-sectional diameter of horizontal cavern, m, L represents the length of horizontal cavern, m, and H indicates the depth of the cavern center, m.Figure 2
Shape simplification of the horizontal cavern. (a) Along the length (x axis). (b) Along the width (y axis).
(a)(b)Specially, along the length direction of the subsidence trough, the shape function can be justified as follows:(7)wx,0=1RxRye−πx/Rx2.Also, along the width of subsidence trough, the shape function is as follows:(8)w0,y=1RxRye−πy/Ry2.


## 2.3. Reduced VolumeVs(t)

Because rock salt has significant creep characteristics, if the brine pressure or natural gas pressure in the horizonal cavern is smaller than the original lithostatic stress, the cavern volume would shrink over time [39, 40]. In fact, the internal pressure cannot be exactly equal to the lithostatic pressure at the cavern depth [9, 11, 14, 16, 39–41]. This indicates that the convergence of the salt cavern will not stop.Because there are a number of insoluble interlayers between the salt layers and the flow of brine is really difficult to exactly control, the actual shapes of the horizontal caverns (Figure1(a)) are still irregular. In addition, the insoluble sediment might accumulate in a large amount at the cavern bottom. With solution mining by double-well convection, the height of horizontal cavern is usually lower than 60 m, which is much smaller than the buried depth (500 to 2000 m). Thus, the horizontal cavern can be assumed to be a cylinder (Figure 2). Since the underground gas storage has the construction period and operation period in a life cycle, the law of reduced volume in these two periods will be discussed below separately.


### 2.3.1. Reduced Volume of the Horizontal Cavern in the Construction Period

In order to avoid the fluctuations of the manual operation, the brine extraction rate is usually kept constant in the construction period. If the salt layer is uniform and the solution rate of the cavern is constant, a series of equations for the reduced volume of the horizontal cavern can be calculated as follows [11]:(9)Vt1=AV•EBk11−e−Bk1t1,Vst1=ABVE+AB−1AV•EBk11−e−Bk1t1,A=ρb−ρwα⋅ρs−ρb−ρw,B=α⋅ρsα⋅ρs−ρb−ρw,t1=tn−t0,where Vt1 represents the horizontal cavern volume at time t1, m3, Vst1 represents the reduced cavern volume at time t1, m3, k1 represents the shrinkage rate of the unit volume in the construction period, m3/m3⋅d, VE represents the extracted brine volume, m3, V˙E indicates the brine extraction rate, m3/d, α indicates the proportion of the soluble components in salt rocks, namely, the purity, %, ρw indicates the water density, 1.0 × 103 kg/m3, ρb indicates the brine density, kg/m3, ρs indicates the rock salt density, kg/m3, t0 indicates the start time of the construction period, month, tn represents the current time, month, and t1 represents the characteristic time, month.In particular, both the cavern volumeVt1and the reduced cavern volumeVst1 are equal to zero when the characteristic time t1=0. Assuming that the time tp is when the solution mining is completed, the cavern volume and reduced cavern volume can be written as(10)Vtp=AV•EBk11−e−Bk1tp−t0,Vstp=ABVE+AB−1AV•EBk11−e−Bk1tp−t0,where Vtp indicates the cavern volume at time tp, m3, and Vstp represents the reduced cavern volume at time tp, m3.


### 2.3.2. Reduced Volume of the Horizontal Cavern in the Operating Period

Assuming that natural gas is immediately stored in the horizontal cavern after the construction, the relationship between the cavern volume and the reduced cavern volume and the characteristic timet2 during operation can also be obtained:(11)Vt2=Vtpe−k2t2,Vst2=Vstp+Vtp1−e−k2t2,t2=tn−tp,where Vt2 indicates the cavern volume at time t2, m3, Vst2 represents the reduced cavern volume at time t2, m3, and k2 represents the shrinkage rate of the unit volume in the operating period of the horizontal cavern, m3/m3⋅d.


### 2.3.3. Shrinkage Rate of the Unit Volumek1 and k2

If only the steady-state stage is considered in the creep deformation of bedded salt formations, the most widely applied creep constitutive model for rock salt is Norton power exponent model [9, 10]:(12)ε•=Cα⋅e−Q/RmTσdiffσ∗n,for uniaxial behavior,ε=•=Cα⋅e−Q/RmT1n+1∂3J2n+1∂σ=,for3D behaviour,J2=12SijSji,Sij=σij−13σiiδij,where C(α) is a material property of the rock salt, 1/month, which is assumed to be related to the grade α, as shown in equation (13), n indicates the stress exponent constant, σ∗ represents the unit stress, 1 MPa, σdiff represents differential stress, MPa, Q represents the activation energy, kJ/mol, Rm is the gas constant, 8.314 × 10−3 kJ/(mol/K), T represents the absolute temperature of the strata at the cavern center, K, which is related to the depth H, as shown in equation (13), J2 is the second invariant of the deviator stress tensor, MPa2, Sij and Sji represent the deviatoric stress tensor, MPa, σij indicates the stress tensor, MPa, σii represents the sum of diagonal elements of the stress tensor, MPa, and δij is the Kronecker delta symbol.Then, the shrinkage rate per unit volume of the cylindrical cavern can be inferred from the 3D formula of the Norton-power exponent model [9–11, 16]:(13)k=33P0−Pn⋅σ∗n⋅Cα⋅e−Q/RmT,Cα=Cs⋅α+Ci⋅1−α,T=273.15+20+0.02⋅H,where k stands for k1andk2and indicates the shrinkage rate of unit volume of a cylindrical salt cavern, m3/m3, P represents the brine pressure P1 during the building period or the gas pressure P2 during the operation period, MPa, P0 represents the original in situ pressure, MPa, Cs and Ci represent the creep properties of the rock salts and interlayers, 1/month, and H is the depth of the cavern center, m.Finally, substituting equation (6) and equations (9)–(13) into equation (1), the value of ground subsidence can be finally obtained, and then, the tilt, curvature, horizontal displacement, and horizontal strain can also be calculated.


## 2.3.1. Reduced Volume of the Horizontal Cavern in the Construction Period

In order to avoid the fluctuations of the manual operation, the brine extraction rate is usually kept constant in the construction period. If the salt layer is uniform and the solution rate of the cavern is constant, a series of equations for the reduced volume of the horizontal cavern can be calculated as follows [11]:(9)Vt1=AV•EBk11−e−Bk1t1,Vst1=ABVE+AB−1AV•EBk11−e−Bk1t1,A=ρb−ρwα⋅ρs−ρb−ρw,B=α⋅ρsα⋅ρs−ρb−ρw,t1=tn−t0,where Vt1 represents the horizontal cavern volume at time t1, m3, Vst1 represents the reduced cavern volume at time t1, m3, k1 represents the shrinkage rate of the unit volume in the construction period, m3/m3⋅d, VE represents the extracted brine volume, m3, V˙E indicates the brine extraction rate, m3/d, α indicates the proportion of the soluble components in salt rocks, namely, the purity, %, ρw indicates the water density, 1.0 × 103 kg/m3, ρb indicates the brine density, kg/m3, ρs indicates the rock salt density, kg/m3, t0 indicates the start time of the construction period, month, tn represents the current time, month, and t1 represents the characteristic time, month.In particular, both the cavern volumeVt1and the reduced cavern volumeVst1 are equal to zero when the characteristic time t1=0. Assuming that the time tp is when the solution mining is completed, the cavern volume and reduced cavern volume can be written as(10)Vtp=AV•EBk11−e−Bk1tp−t0,Vstp=ABVE+AB−1AV•EBk11−e−Bk1tp−t0,where Vtp indicates the cavern volume at time tp, m3, and Vstp represents the reduced cavern volume at time tp, m3.


## 2.3.2. Reduced Volume of the Horizontal Cavern in the Operating Period

Assuming that natural gas is immediately stored in the horizontal cavern after the construction, the relationship between the cavern volume and the reduced cavern volume and the characteristic timet2 during operation can also be obtained:(11)Vt2=Vtpe−k2t2,Vst2=Vstp+Vtp1−e−k2t2,t2=tn−tp,where Vt2 indicates the cavern volume at time t2, m3, Vst2 represents the reduced cavern volume at time t2, m3, and k2 represents the shrinkage rate of the unit volume in the operating period of the horizontal cavern, m3/m3⋅d.


## 2.3.3. Shrinkage Rate of the Unit Volumek1 and k2

If only the steady-state stage is considered in the creep deformation of bedded salt formations, the most widely applied creep constitutive model for rock salt is Norton power exponent model [9, 10]:(12)ε•=Cα⋅e−Q/RmTσdiffσ∗n,for uniaxial behavior,ε=•=Cα⋅e−Q/RmT1n+1∂3J2n+1∂σ=,for3D behaviour,J2=12SijSji,Sij=σij−13σiiδij,where C(α) is a material property of the rock salt, 1/month, which is assumed to be related to the grade α, as shown in equation (13), n indicates the stress exponent constant, σ∗ represents the unit stress, 1 MPa, σdiff represents differential stress, MPa, Q represents the activation energy, kJ/mol, Rm is the gas constant, 8.314 × 10−3 kJ/(mol/K), T represents the absolute temperature of the strata at the cavern center, K, which is related to the depth H, as shown in equation (13), J2 is the second invariant of the deviator stress tensor, MPa2, Sij and Sji represent the deviatoric stress tensor, MPa, σij indicates the stress tensor, MPa, σii represents the sum of diagonal elements of the stress tensor, MPa, and δij is the Kronecker delta symbol.Then, the shrinkage rate per unit volume of the cylindrical cavern can be inferred from the 3D formula of the Norton-power exponent model [9–11, 16]:(13)k=33P0−Pn⋅σ∗n⋅Cα⋅e−Q/RmT,Cα=Cs⋅α+Ci⋅1−α,T=273.15+20+0.02⋅H,where k stands for k1andk2and indicates the shrinkage rate of unit volume of a cylindrical salt cavern, m3/m3, P represents the brine pressure P1 during the building period or the gas pressure P2 during the operation period, MPa, P0 represents the original in situ pressure, MPa, Cs and Ci represent the creep properties of the rock salts and interlayers, 1/month, and H is the depth of the cavern center, m.Finally, substituting equation (6) and equations (9)–(13) into equation (1), the value of ground subsidence can be finally obtained, and then, the tilt, curvature, horizontal displacement, and horizontal strain can also be calculated.


## 3. Verification of the Proposed Model and Case Analysis

### 3.1. Basic Parameters of the Horizontal Cavern

MathCAD software was introduced to realize the above equations. Assume that the depth of horizontal cavern center (H) was −1000 m, with a cross-sectional diameter (D) of 40 m and a length (L) of 200 m. In salt formations, it is usually assumed that the natural stress state generated by the overburden weight is isotropic [9, 11, 16]. If the average density of the salt formations and overburden strata ρs was 2.30 × 103 kg/m3, the original in situ pressure P0 would be 23 MPa. If the saturated density of the brine ρb was set to be 1.20 × 103 kg/m3, the brine pressure P1 at the cavern center would be 12 MPa.The values of the basic parameters for the horizontal cavern were referred from research [11]. In detail, the purity α was 0.85, the subsidence adjusting coefficient a was supposed to be 0.4, and the draw angle β was set to be 37°. VE˙ was primarily set to be 6 × 104 m3/month, while the creep property Cs and Ci were 6 × 10−3 month−1 and 2 × 10−5 month−1, respectively. The activity energy Q was set to be 16 kJ/mol. In addition, the solution mining lasted for approximately 45.4 months during the construction period, while the operation period was about 30 years with the gas pressure P2 = 15 MPa. Moreover, the other parameters can be obtained from the above parameters.


### 3.2. Verification of the Proposed Model

As a verification of the proposed model, the software FLAC3D was adopted to carry out stability analysis, while the software ANSYS was used to establish the model in advance [1, 11, 13, 16, 30, 32]. Figure 3 shows a diagram of the geological strata and the horizontal cavern. Because of the symmetry, only a quarter of the horizontal cavern was established in the numerical model. The center of the horizontal cavern was set to be the origin of the entire numerical model, and the x and y coordinate plane is the plane passing through the cavern center. The x axis is along the length direction, the y axis is along the width direction, and the z axis is along the vertical direction. The entire height of the numerical model is 1240 m, with 1000 m above the x and y coordinate plane and −240 m below. The width and length of the numerical model is 2000 m. There are four generalized strata in the entire model, respectively, the surface sedimentary strata (260 m), the upper mudstone (660 m), the bedded salt formations (160 m), and the lower mudstone (160 m).Figure 3
Diagram of the geological strata and the horizontal cavern in FLAC3D.The bottom surface of the numerical model was fixed along thez direction, while the four sides of the model were fixed in the normal. The CPower model was adopted for rock salt, mudstone, and sedimentary stratum, and the required visco-elastic-plastic parameters in FLAC3D are detailed listed in Table 1.Table 1
Visco-elastic-plastic parameters for the numerical simulation.
StrataRock saltMudstoneSedimentary stratumCreep parameterAc (a−1)1.6 × 10−45.2 × 10−75.2 × 10−7Stress index constantn3.83.83.8Young’s modulusE (GPa)18.010.00.5Poisson’s ratioμ0.300.270.35Cohesionc (MPa)1.01.00.5Internal friction angleθ (°)403535Tensile strengthσt (MPa)1.01.00.5Specially, the parameterAc is calculated as follows:(14)Ac=12Cs⋅e−Q/RmT,for rock salt,12Ci⋅e−Q/RmT,for not rock salt.Using this model and FLAC3D, the ground subsidence along thex axis after 30 years of operation was calculated, respectively, and the results can be found in Figure 4. It can be seen that the shapes of the two curves are almost the same. Compared with the calculation result of FLAC3D, the prediction result of the ground subsidence by the proposed formulas is somewhat overestimated from the ground center out to around 1200 m and somewhat underestimated at the model boundary. In detail, the maximum subsidence calculated by the proposed model is 44.79 mm, while it is 43.69 mm by the FLAC3D. The difference between the results by the two methods is only 1.10 mm, which occupies only 2.5% of the result by the proposed model. Therefore, the proposed prediction model is feasible.Figure 4
Ground subsidence alongx axis after 30 years of operations by the proposed model and FLAC3D.As a matter of fact, the horizontal caverns cannot be seen as idealized cylinders. Considering the idealized assumption of the ideal cylinder used in the model and the subsidence equations, the difference in the ground subsidence profiles between the proposed model and the in situ data is inevitable. Therefore, the proposed prediction model is rather ideal and needs to be improved further for being applied in the gas storage engineering.


### 3.3. Case Analysis

Using MATLAB, the subsidence trough and contour lines can be figured out, as shown in Figure5. From Figure 5(a), the largest subsidence happens at the center, while the subsidence becomes smaller as the distance farther away from the center. From Figure 5(b), the shape of each contour line looks like a circle but actually an ellipse.Figure 5
Ground subsidence trough and contour lines above a horizontal salt cavern by MATLAB. (a) Subsidence trough. (b) Subsidence contour.
(a)(b)The proposed model was used to calculate thex-direction subsidence at the year 0, 5, 10, 15, 20, 25, or 30 of the operation periods, and the results are shown in Figure 6. In particular, the operating year 0 refers to the year when the horizontal cavern is completely built. The maximum subsidence is 11.04, 19.34, 26.32, 32.19, 37.14, 41.29, or 44.79 mm, respectively. Figure 7 shows the subsidence tilt above horizontal cavern after 30 years of operation. The maximum tilt appears at about 580 m from the center the ground trough, and its value is about 0.0477 × 10−3 m/m, which is much smaller than the maximum allowable value (2.0 × 10−3 m/m) according to the Chinese Building Standard [11].Figure 6
Subsidence curves during the operating period alongx axis.Figure 7
Ground subsidence tilt above horizontal cavern after 30 years of operation alongx axis.


## 3.1. Basic Parameters of the Horizontal Cavern

MathCAD software was introduced to realize the above equations. Assume that the depth of horizontal cavern center (H) was −1000 m, with a cross-sectional diameter (D) of 40 m and a length (L) of 200 m. In salt formations, it is usually assumed that the natural stress state generated by the overburden weight is isotropic [9, 11, 16]. If the average density of the salt formations and overburden strata ρs was 2.30 × 103 kg/m3, the original in situ pressure P0 would be 23 MPa. If the saturated density of the brine ρb was set to be 1.20 × 103 kg/m3, the brine pressure P1 at the cavern center would be 12 MPa.The values of the basic parameters for the horizontal cavern were referred from research [11]. In detail, the purity α was 0.85, the subsidence adjusting coefficient a was supposed to be 0.4, and the draw angle β was set to be 37°. VE˙ was primarily set to be 6 × 104 m3/month, while the creep property Cs and Ci were 6 × 10−3 month−1 and 2 × 10−5 month−1, respectively. The activity energy Q was set to be 16 kJ/mol. In addition, the solution mining lasted for approximately 45.4 months during the construction period, while the operation period was about 30 years with the gas pressure P2 = 15 MPa. Moreover, the other parameters can be obtained from the above parameters.


## 3.2. Verification of the Proposed Model

As a verification of the proposed model, the software FLAC3D was adopted to carry out stability analysis, while the software ANSYS was used to establish the model in advance [1, 11, 13, 16, 30, 32]. Figure 3 shows a diagram of the geological strata and the horizontal cavern. Because of the symmetry, only a quarter of the horizontal cavern was established in the numerical model. The center of the horizontal cavern was set to be the origin of the entire numerical model, and the x and y coordinate plane is the plane passing through the cavern center. The x axis is along the length direction, the y axis is along the width direction, and the z axis is along the vertical direction. The entire height of the numerical model is 1240 m, with 1000 m above the x and y coordinate plane and −240 m below. The width and length of the numerical model is 2000 m. There are four generalized strata in the entire model, respectively, the surface sedimentary strata (260 m), the upper mudstone (660 m), the bedded salt formations (160 m), and the lower mudstone (160 m).Figure 3
Diagram of the geological strata and the horizontal cavern in FLAC3D.The bottom surface of the numerical model was fixed along thez direction, while the four sides of the model were fixed in the normal. The CPower model was adopted for rock salt, mudstone, and sedimentary stratum, and the required visco-elastic-plastic parameters in FLAC3D are detailed listed in Table 1.Table 1
Visco-elastic-plastic parameters for the numerical simulation.
StrataRock saltMudstoneSedimentary stratumCreep parameterAc (a−1)1.6 × 10−45.2 × 10−75.2 × 10−7Stress index constantn3.83.83.8Young’s modulusE (GPa)18.010.00.5Poisson’s ratioμ0.300.270.35Cohesionc (MPa)1.01.00.5Internal friction angleθ (°)403535Tensile strengthσt (MPa)1.01.00.5Specially, the parameterAc is calculated as follows:(14)Ac=12Cs⋅e−Q/RmT,for rock salt,12Ci⋅e−Q/RmT,for not rock salt.Using this model and FLAC3D, the ground subsidence along thex axis after 30 years of operation was calculated, respectively, and the results can be found in Figure 4. It can be seen that the shapes of the two curves are almost the same. Compared with the calculation result of FLAC3D, the prediction result of the ground subsidence by the proposed formulas is somewhat overestimated from the ground center out to around 1200 m and somewhat underestimated at the model boundary. In detail, the maximum subsidence calculated by the proposed model is 44.79 mm, while it is 43.69 mm by the FLAC3D. The difference between the results by the two methods is only 1.10 mm, which occupies only 2.5% of the result by the proposed model. Therefore, the proposed prediction model is feasible.Figure 4
Ground subsidence alongx axis after 30 years of operations by the proposed model and FLAC3D.As a matter of fact, the horizontal caverns cannot be seen as idealized cylinders. Considering the idealized assumption of the ideal cylinder used in the model and the subsidence equations, the difference in the ground subsidence profiles between the proposed model and the in situ data is inevitable. Therefore, the proposed prediction model is rather ideal and needs to be improved further for being applied in the gas storage engineering.


## 3.3. Case Analysis

Using MATLAB, the subsidence trough and contour lines can be figured out, as shown in Figure5. From Figure 5(a), the largest subsidence happens at the center, while the subsidence becomes smaller as the distance farther away from the center. From Figure 5(b), the shape of each contour line looks like a circle but actually an ellipse.Figure 5
Ground subsidence trough and contour lines above a horizontal salt cavern by MATLAB. (a) Subsidence trough. (b) Subsidence contour.
(a)(b)The proposed model was used to calculate thex-direction subsidence at the year 0, 5, 10, 15, 20, 25, or 30 of the operation periods, and the results are shown in Figure 6. In particular, the operating year 0 refers to the year when the horizontal cavern is completely built. The maximum subsidence is 11.04, 19.34, 26.32, 32.19, 37.14, 41.29, or 44.79 mm, respectively. Figure 7 shows the subsidence tilt above horizontal cavern after 30 years of operation. The maximum tilt appears at about 580 m from the center the ground trough, and its value is about 0.0477 × 10−3 m/m, which is much smaller than the maximum allowable value (2.0 × 10−3 m/m) according to the Chinese Building Standard [11].Figure 6
Subsidence curves during the operating period alongx axis.Figure 7
Ground subsidence tilt above horizontal cavern after 30 years of operation alongx axis.


## 4. Parameter Analysis by the Proposed Prediction Model

According to theoretical analysis in Section3, there are many influence parameters of the ground subsidence above horizontal salt caverns, for example, the draw angle β, the subsidence adjusting coefficient a, the purity α, the creep characteristics n and Cs of rock salt, the cavern depth H, the brine extraction velocity VE˙, the cross-sectional diameter D and length L of the horizontal cavern, and the internal pressure P2. In these parameters, the draw angle β, the subsidence adjusting coefficient a, the purity α, the creep characteristics n and Cs of rock salt, and the cavern depth H belong to the stratigraphic parameters, while the brine extraction velocity VE˙, the cross-sectional diameter D and length L of the horizontal cavern, and the internal pressure P2 are the design parameters. With all other parameters listed above in Section 3, the impacts of the main parameters on the ground subsidence along x axis were analyzed through the proposed prediction model after 30 years of operation.


### 4.1. Impacts of the Stratigraphic Parameters

#### 4.1.1. The Draw Angleβ

The impact of the draw angleβ on the ground subsidence is shown in Figures 8(a) and 9. As the draw angle β varies from 35° to 39°, the subsidence trough has tremendous shape changes (Figure 9). It can be seen that the maximum subsidence gradually increases at the ground, but the impact on the boundary of the model becomes smaller, and the turning point is about 750 m from the ground center. In other words, the subsidence trough becomes deeper and narrower with the bigger of the draw angle β. However, the accurate value of the draw angle β needs to be inferred from the field monitoring data.Figure 8
Influence of the stratigraphic parameters on the ground subsidence.
(a)(b)(c)(d)(e)(f)Figure 9
Impact of draw angleβ on the ground subsidence.


#### 4.1.2. The Subsidence Adjusting Coefficienta

The impact of the subsidence adjusting coefficienta on the ground subsidence is shown in Figure 8(b). It can be seen from the figure that, as the subsidence adjusting coefficient a increases from 0.3 to 0.5, the maximum ground subsidence increases proportionally at the ground, which can also be predicted from equation (1). In view that the subsidence adjusting coefficient a relies fiercely on the characteristic of the overlying strata, its accurate value needs to be inferred from the field monitoring.


#### 4.1.3. The Cavern DepthH

The impact of the cavern depthH on the ground subsidence is shown in Figures 8(c) and 10. As the horizontal cavern buried deeper from −600 m to −1600 m, the area with large subsidence on the ground surface increases gradually. However, the value of the maximum subsidence changes with an inflection, and the inflection depth is approximately −1000 m. The maximum subsidence gradually increases with the depth shallower than −1000 m and decreases with the depth deeper than −1000 m.Figure 10
Impact of cavern depthH on the ground subsidence.With the increasing of the depth of the cavern, the deviatoric stress between the brine pressure or internal pressure and the original in situ stress increases, and the temperature of the stratum surrounding the cavern also increases. The increase of deviator stress and temperature will promote the creep deformation of rock salts, raise the volume shrinkage rate, and finally increase the ground subsidence. However, if the parametersL, β, and D remain constant, the impact area of the ground subsidence would also enlarge with the increase of the depth H. Moreover, under the same amount of ground subsidence volume, the larger the impact area, the smaller the maximum subsidence. Therefore, if the impact area increases slowly at the ground, the maximum subsidence would definitely increase; if the impact area increases rapidly at the ground, the maximum subsidence might decrease. As a result, as the depth increases, the maximum ground subsidence increases first and then decreases, as shown in Figure 8(c).In addition, the buried depth also affects the shape of subsidence trough. Definingλ indicates the ratio of the major axis to the minor axis of the subsidence trough, as follows:(15)λ=RxRy=Hcotβ+L/2Hcotβ+D/2.Assuming that the draw angleβ is 37°, the length L is 200 m, and the cross-sectional diameter D is 40 m, the impact of the cavern depth H on the shape of the ground trough can be figured out, as shown in Figure 11, with the subsidence contour 0.01 m. From the figure, with the increasing of the cavern depth, the ratio λ gradually tends to be 1 (Figure 11(b)), which means that the subsidence trough gradually tends to be a cycle (Figure 11(a)).Figure 11
Impact of cavern depthH on the shape of the subsidence trough with the subsidence contour 0.01 m. (a) Shape of the subsidence trough. (b) Ratio of the major axis to the minor axis.
(a)(b)


#### 4.1.4. The Purityα

By setting the volume of the cavern to be a constant, the effect of purityα on ground subsidence was carried out. The calculating results are shown in Figure 8(d). It can be seen from the figure that, as the purity increases from 0.75 to 0.95, the subsidence gradually decreases from 48.58 mm to 41.61 mm. Although the differences are not too large, the impact of the purity cannot be ignored. The reason is that building storage caverns in high-purity salt formations cannot only significantly reduce sediment and improve the utilization rate of the dissolved cavern but also reduce solution time and corresponding subsidence during the construction period.


#### 4.1.5. The Stress Index Constantn

The impact of the stress index constantn on the ground subsidence is shown in Figure 8(e). Along with the increasing of the stress index constant n from 3.0 to 6.0, the maximum subsidence first increases (from 3.0 to 5.5) and then decreases slightly (from 5.5 to 6.0), and the inflection point appears at about n = 5.5. The reason is that the unit volume shrinkage rate k is not linear to the stress index constant n according to equation (13). Therefore, under the premise of controlling ground subsidence, salt rocks with appropriate exponent constant n should be selected to build the underground gas storage cavern.


#### 4.1.6. The Creep PropertyCs

The impact of the creep propertyCs of the rock salts on the ground subsidence is shown in Figure 8(f). With the increasing of the creep property Cs from 2 × 10−3 month−1 to 20 × 10−3 month−1, the maximum ground subsidence increases rapidly first and then slowly approaches its limit. As a matter of fact, the bigger the creep property Cs, the rapider the creep deformation of the rock salts. If the creep property Cs is large enough, the cavern will shrink to its limit (namely, the cavern will eventually close), which leads to the subsidence towards its limit. In addition, the unit volume shrinkage rate k has a linear relationship with Cs from equation (13), but has a nonlinear relationship with the reduced volume Vs(t) from equations (9)–(11). In view that Vs(t) is linear with the subsidence wx,y,t from equation (1), the creep property of rock salt Cs will be nonlinear with subsidence wx,y,t. In short, the greater the creep property Cs, the greater the ground subsidence, thus the steady-state creep parameter should be given exactly when predicting.Above all, the research studies on the impact of stratigraphic parameters on ground subsidence are conducive to the site selection of horizontal caverns for gas storage. Appropriate selection of bedded salt formations can reduce the subsidence risk during the entire construction and operation period. However, once the location of the cavern is selected, the values of the stratigraphic parameters will be specific and cannot be changed. Moreover, the exact value of the stratigraphic parameters needs to be inferred from the field monitoring data.


### 4.2. Impacts of the Design Parameters

#### 4.2.1. The Brine Extraction VelocityVE˙

According to research [11], if the brine extraction velocity is faster in the construction period, the natural gas could be stored earlier. From an economic point of view, it is appropriate to increase the brine extraction rate. In fact, once the final cross-sectional diameter of the built cavern is determined, the required building time will become shorter with the increasing of the brine extraction velocity. If the prescribed cross-sectional diameter of the cavern is 40 m, assuming the brine extraction velocity is 6.0, 6.5, 7.0, 7.5, or 8.0 × 104 m3/month, correspondingly the building time will be 45.4, 41.6, 38.0, 35.0, or 32.4 months, respectively, by equation (10).Since the brine pressure in the cavern during the construction period is lower than the internal pressure during gas storage, the deviator stress on the surrounding salt is relatively large, and the cavern will inevitably shrink more. Therefore, increasing the speed of brine extraction can shorten the cavern construction period, thereby reducing the excessive shrinkage of the cavern and thus decreasing the excessive ground subsidence, as shown in Figure12(a).Figure 12
Impact of design parameters on the subsidence.
(a)(b)(c)(d)


#### 4.2.2. The Internal Pressure in the Horizontal CavernP2

In order to study the impact of the internal pressure on the ground subsidence, the internal pressureP2 was kept at 9, 11, 13, 15, or 17 MPa as a constant in each case. The calculating results are as shown in Figure 12(b). It can be seen from the figure that the ground subsidence decreases nonlinearly with the increasing of the internal pressure. This reason is that the increasing of the internal pressure reduces the deviatoric stress of the rock salt, which leads to a decrease of the volume shrinkage and ground subsidence.Specially, when the internal pressure is low (9 and 11 MPa), after 30 years of creep, the ground subsidence tends to its maximum, which would happen when the cavern totally closed. When the internal pressure is high (19 and 21 MPa), the ground subsidence tends to be the smallest, which is equivalent to the value after the end of the construction period. As a matter of fact, the total subsidence includes both the subsidence at the construction period and that at the gas storage period. Although changes of the internal gas pressure during operation may change the ground subsidence, the subsidence during the construction period should be constant. Moreover, the internal pressure in the horizontal cavern should be not too high, which may induce tensile fracture in the surrounding rock salts at the top of the cavern and even other parts. Therefore, the internal pressure should be raised appropriately.


#### 4.2.3. The Cross-Sectional DiameterD

As it is known to all, it is much better to obtain a larger horizontal cavern for natural gas storage. If the length of the cavern is constant (200 m), the cross-sectional diameter of the horizontal cavern will determine the whole volume. Assuming the brine extraction velocity is kept 6.0 × 104 m3/month as a constant, if the cross-sectional diameters of the horizontal caverns are 20, 30, 40, 50, or 60 m, respectively, the cavern volume Vtp will be 6.28, 14.13, 25.12, 39.25, or 56.52 × 104 m3 by equation (10), and then, the building time will be, respectively, 9.6, 23.1, 45.4, 83.7, or 165.8 months. It can be seen that, as the diameter of the cross-section increases, the required time to build the cavern will become longer. In particular, it takes nearly 14 years to build a horizontal cavern with a diameter of 60 m and a length of 200 m. This is obviously uneconomical for the storage company, and thus, a reasonable choice of the cavern diameter is required.The reason why it needs such a long time is that the shrinkage of the cavern gradually increases with the increasing of the cavern volume. In fact, it should dissolve much more rock salt formations to form a horizontal cavern with diameter 60 m and length 200 m. Under a longer construction period, the salt cavern will shrink more because the brine pressure during the construction period is much lower than pressure during gas storage. Therefore, when the cross-sectional diameter increases proportionally, the required construction time will increase nonlinearly.Figure12(c) shows the impact of different cross-sectional diameters on the ground subsidence. From the figure, the ground subsidence increases obviously with the increasing of the cross-sectional diameter, which is disadvantageous to promote a larger volume capacity of the storage. Therefore, in the shape design period, an appropriate cross-sectional diameter should be selected.


#### 4.2.4. The Length of the Horizontal CavernL

If the cross-sectional diameterD of the horizontal cavern is kept 40 m as a constant, the length L will determine the volume of the salt cavern. Assuming the brine extraction velocity is kept constant 6.0 × 104 m3/month, if the length of the cavern is set to be 100, 150, 200, 250, or 300 m respectively, correspondingly the cavern volume Vtp will be 12.56, 18.84, 25.12, 31.40, or 37.68 × 104 m3 by equation (10), and then, the building time will be 20.3, 32.0, 45.4, 60.8, or 78.5 months, respectively. It can also be seen that the required building time will become longer with the increasing of the length.Figure12(d) shows the impact of different lengths on the ground subsidence. From the figure, the ground subsidence increases obviously with the increasing of the cavern length. Thus, an appropriate length should be selected for the economic and safety considerations. In addition, according to equation (15), as the length increases, the ratio of the major axis to the minor axis of the subsidence trough (λ) becomes larger, which indicates that the changing of the length will also change the roundness of the subsidence contour.


## 4.1. Impacts of the Stratigraphic Parameters

### 4.1.1. The Draw Angleβ

The impact of the draw angleβ on the ground subsidence is shown in Figures 8(a) and 9. As the draw angle β varies from 35° to 39°, the subsidence trough has tremendous shape changes (Figure 9). It can be seen that the maximum subsidence gradually increases at the ground, but the impact on the boundary of the model becomes smaller, and the turning point is about 750 m from the ground center. In other words, the subsidence trough becomes deeper and narrower with the bigger of the draw angle β. However, the accurate value of the draw angle β needs to be inferred from the field monitoring data.Figure 8
Influence of the stratigraphic parameters on the ground subsidence.
(a)(b)(c)(d)(e)(f)Figure 9
Impact of draw angleβ on the ground subsidence.


### 4.1.2. The Subsidence Adjusting Coefficienta

The impact of the subsidence adjusting coefficienta on the ground subsidence is shown in Figure 8(b). It can be seen from the figure that, as the subsidence adjusting coefficient a increases from 0.3 to 0.5, the maximum ground subsidence increases proportionally at the ground, which can also be predicted from equation (1). In view that the subsidence adjusting coefficient a relies fiercely on the characteristic of the overlying strata, its accurate value needs to be inferred from the field monitoring.


### 4.1.3. The Cavern DepthH

The impact of the cavern depthH on the ground subsidence is shown in Figures 8(c) and 10. As the horizontal cavern buried deeper from −600 m to −1600 m, the area with large subsidence on the ground surface increases gradually. However, the value of the maximum subsidence changes with an inflection, and the inflection depth is approximately −1000 m. The maximum subsidence gradually increases with the depth shallower than −1000 m and decreases with the depth deeper than −1000 m.Figure 10
Impact of cavern depthH on the ground subsidence.With the increasing of the depth of the cavern, the deviatoric stress between the brine pressure or internal pressure and the original in situ stress increases, and the temperature of the stratum surrounding the cavern also increases. The increase of deviator stress and temperature will promote the creep deformation of rock salts, raise the volume shrinkage rate, and finally increase the ground subsidence. However, if the parametersL, β, and D remain constant, the impact area of the ground subsidence would also enlarge with the increase of the depth H. Moreover, under the same amount of ground subsidence volume, the larger the impact area, the smaller the maximum subsidence. Therefore, if the impact area increases slowly at the ground, the maximum subsidence would definitely increase; if the impact area increases rapidly at the ground, the maximum subsidence might decrease. As a result, as the depth increases, the maximum ground subsidence increases first and then decreases, as shown in Figure 8(c).In addition, the buried depth also affects the shape of subsidence trough. Definingλ indicates the ratio of the major axis to the minor axis of the subsidence trough, as follows:(15)λ=RxRy=Hcotβ+L/2Hcotβ+D/2.Assuming that the draw angleβ is 37°, the length L is 200 m, and the cross-sectional diameter D is 40 m, the impact of the cavern depth H on the shape of the ground trough can be figured out, as shown in Figure 11, with the subsidence contour 0.01 m. From the figure, with the increasing of the cavern depth, the ratio λ gradually tends to be 1 (Figure 11(b)), which means that the subsidence trough gradually tends to be a cycle (Figure 11(a)).Figure 11
Impact of cavern depthH on the shape of the subsidence trough with the subsidence contour 0.01 m. (a) Shape of the subsidence trough. (b) Ratio of the major axis to the minor axis.
(a)(b)


### 4.1.4. The Purityα

By setting the volume of the cavern to be a constant, the effect of purityα on ground subsidence was carried out. The calculating results are shown in Figure 8(d). It can be seen from the figure that, as the purity increases from 0.75 to 0.95, the subsidence gradually decreases from 48.58 mm to 41.61 mm. Although the differences are not too large, the impact of the purity cannot be ignored. The reason is that building storage caverns in high-purity salt formations cannot only significantly reduce sediment and improve the utilization rate of the dissolved cavern but also reduce solution time and corresponding subsidence during the construction period.


### 4.1.5. The Stress Index Constantn

The impact of the stress index constantn on the ground subsidence is shown in Figure 8(e). Along with the increasing of the stress index constant n from 3.0 to 6.0, the maximum subsidence first increases (from 3.0 to 5.5) and then decreases slightly (from 5.5 to 6.0), and the inflection point appears at about n = 5.5. The reason is that the unit volume shrinkage rate k is not linear to the stress index constant n according to equation (13). Therefore, under the premise of controlling ground subsidence, salt rocks with appropriate exponent constant n should be selected to build the underground gas storage cavern.


### 4.1.6. The Creep PropertyCs

The impact of the creep propertyCs of the rock salts on the ground subsidence is shown in Figure 8(f). With the increasing of the creep property Cs from 2 × 10−3 month−1 to 20 × 10−3 month−1, the maximum ground subsidence increases rapidly first and then slowly approaches its limit. As a matter of fact, the bigger the creep property Cs, the rapider the creep deformation of the rock salts. If the creep property Cs is large enough, the cavern will shrink to its limit (namely, the cavern will eventually close), which leads to the subsidence towards its limit. In addition, the unit volume shrinkage rate k has a linear relationship with Cs from equation (13), but has a nonlinear relationship with the reduced volume Vs(t) from equations (9)–(11). In view that Vs(t) is linear with the subsidence wx,y,t from equation (1), the creep property of rock salt Cs will be nonlinear with subsidence wx,y,t. In short, the greater the creep property Cs, the greater the ground subsidence, thus the steady-state creep parameter should be given exactly when predicting.Above all, the research studies on the impact of stratigraphic parameters on ground subsidence are conducive to the site selection of horizontal caverns for gas storage. Appropriate selection of bedded salt formations can reduce the subsidence risk during the entire construction and operation period. However, once the location of the cavern is selected, the values of the stratigraphic parameters will be specific and cannot be changed. Moreover, the exact value of the stratigraphic parameters needs to be inferred from the field monitoring data.


## 4.1.1. The Draw Angleβ

The impact of the draw angleβ on the ground subsidence is shown in Figures 8(a) and 9. As the draw angle β varies from 35° to 39°, the subsidence trough has tremendous shape changes (Figure 9). It can be seen that the maximum subsidence gradually increases at the ground, but the impact on the boundary of the model becomes smaller, and the turning point is about 750 m from the ground center. In other words, the subsidence trough becomes deeper and narrower with the bigger of the draw angle β. However, the accurate value of the draw angle β needs to be inferred from the field monitoring data.Figure 8
Influence of the stratigraphic parameters on the ground subsidence.
(a)(b)(c)(d)(e)(f)Figure 9
Impact of draw angleβ on the ground subsidence.


## 4.1.2. The Subsidence Adjusting Coefficienta

The impact of the subsidence adjusting coefficienta on the ground subsidence is shown in Figure 8(b). It can be seen from the figure that, as the subsidence adjusting coefficient a increases from 0.3 to 0.5, the maximum ground subsidence increases proportionally at the ground, which can also be predicted from equation (1). In view that the subsidence adjusting coefficient a relies fiercely on the characteristic of the overlying strata, its accurate value needs to be inferred from the field monitoring.


## 4.1.3. The Cavern DepthH

The impact of the cavern depthH on the ground subsidence is shown in Figures 8(c) and 10. As the horizontal cavern buried deeper from −600 m to −1600 m, the area with large subsidence on the ground surface increases gradually. However, the value of the maximum subsidence changes with an inflection, and the inflection depth is approximately −1000 m. The maximum subsidence gradually increases with the depth shallower than −1000 m and decreases with the depth deeper than −1000 m.Figure 10
Impact of cavern depthH on the ground subsidence.With the increasing of the depth of the cavern, the deviatoric stress between the brine pressure or internal pressure and the original in situ stress increases, and the temperature of the stratum surrounding the cavern also increases. The increase of deviator stress and temperature will promote the creep deformation of rock salts, raise the volume shrinkage rate, and finally increase the ground subsidence. However, if the parametersL, β, and D remain constant, the impact area of the ground subsidence would also enlarge with the increase of the depth H. Moreover, under the same amount of ground subsidence volume, the larger the impact area, the smaller the maximum subsidence. Therefore, if the impact area increases slowly at the ground, the maximum subsidence would definitely increase; if the impact area increases rapidly at the ground, the maximum subsidence might decrease. As a result, as the depth increases, the maximum ground subsidence increases first and then decreases, as shown in Figure 8(c).In addition, the buried depth also affects the shape of subsidence trough. Definingλ indicates the ratio of the major axis to the minor axis of the subsidence trough, as follows:(15)λ=RxRy=Hcotβ+L/2Hcotβ+D/2.Assuming that the draw angleβ is 37°, the length L is 200 m, and the cross-sectional diameter D is 40 m, the impact of the cavern depth H on the shape of the ground trough can be figured out, as shown in Figure 11, with the subsidence contour 0.01 m. From the figure, with the increasing of the cavern depth, the ratio λ gradually tends to be 1 (Figure 11(b)), which means that the subsidence trough gradually tends to be a cycle (Figure 11(a)).Figure 11
Impact of cavern depthH on the shape of the subsidence trough with the subsidence contour 0.01 m. (a) Shape of the subsidence trough. (b) Ratio of the major axis to the minor axis.
(a)(b)


## 4.1.4. The Purityα

By setting the volume of the cavern to be a constant, the effect of purityα on ground subsidence was carried out. The calculating results are shown in Figure 8(d). It can be seen from the figure that, as the purity increases from 0.75 to 0.95, the subsidence gradually decreases from 48.58 mm to 41.61 mm. Although the differences are not too large, the impact of the purity cannot be ignored. The reason is that building storage caverns in high-purity salt formations cannot only significantly reduce sediment and improve the utilization rate of the dissolved cavern but also reduce solution time and corresponding subsidence during the construction period.


## 4.1.5. The Stress Index Constantn

The impact of the stress index constantn on the ground subsidence is shown in Figure 8(e). Along with the increasing of the stress index constant n from 3.0 to 6.0, the maximum subsidence first increases (from 3.0 to 5.5) and then decreases slightly (from 5.5 to 6.0), and the inflection point appears at about n = 5.5. The reason is that the unit volume shrinkage rate k is not linear to the stress index constant n according to equation (13). Therefore, under the premise of controlling ground subsidence, salt rocks with appropriate exponent constant n should be selected to build the underground gas storage cavern.


## 4.1.6. The Creep PropertyCs

The impact of the creep propertyCs of the rock salts on the ground subsidence is shown in Figure 8(f). With the increasing of the creep property Cs from 2 × 10−3 month−1 to 20 × 10−3 month−1, the maximum ground subsidence increases rapidly first and then slowly approaches its limit. As a matter of fact, the bigger the creep property Cs, the rapider the creep deformation of the rock salts. If the creep property Cs is large enough, the cavern will shrink to its limit (namely, the cavern will eventually close), which leads to the subsidence towards its limit. In addition, the unit volume shrinkage rate k has a linear relationship with Cs from equation (13), but has a nonlinear relationship with the reduced volume Vs(t) from equations (9)–(11). In view that Vs(t) is linear with the subsidence wx,y,t from equation (1), the creep property of rock salt Cs will be nonlinear with subsidence wx,y,t. In short, the greater the creep property Cs, the greater the ground subsidence, thus the steady-state creep parameter should be given exactly when predicting.Above all, the research studies on the impact of stratigraphic parameters on ground subsidence are conducive to the site selection of horizontal caverns for gas storage. Appropriate selection of bedded salt formations can reduce the subsidence risk during the entire construction and operation period. However, once the location of the cavern is selected, the values of the stratigraphic parameters will be specific and cannot be changed. Moreover, the exact value of the stratigraphic parameters needs to be inferred from the field monitoring data.


## 4.2. Impacts of the Design Parameters

### 4.2.1. The Brine Extraction VelocityVE˙

According to research [11], if the brine extraction velocity is faster in the construction period, the natural gas could be stored earlier. From an economic point of view, it is appropriate to increase the brine extraction rate. In fact, once the final cross-sectional diameter of the built cavern is determined, the required building time will become shorter with the increasing of the brine extraction velocity. If the prescribed cross-sectional diameter of the cavern is 40 m, assuming the brine extraction velocity is 6.0, 6.5, 7.0, 7.5, or 8.0 × 104 m3/month, correspondingly the building time will be 45.4, 41.6, 38.0, 35.0, or 32.4 months, respectively, by equation (10).Since the brine pressure in the cavern during the construction period is lower than the internal pressure during gas storage, the deviator stress on the surrounding salt is relatively large, and the cavern will inevitably shrink more. Therefore, increasing the speed of brine extraction can shorten the cavern construction period, thereby reducing the excessive shrinkage of the cavern and thus decreasing the excessive ground subsidence, as shown in Figure12(a).Figure 12
Impact of design parameters on the subsidence.
(a)(b)(c)(d)


### 4.2.2. The Internal Pressure in the Horizontal CavernP2

In order to study the impact of the internal pressure on the ground subsidence, the internal pressureP2 was kept at 9, 11, 13, 15, or 17 MPa as a constant in each case. The calculating results are as shown in Figure 12(b). It can be seen from the figure that the ground subsidence decreases nonlinearly with the increasing of the internal pressure. This reason is that the increasing of the internal pressure reduces the deviatoric stress of the rock salt, which leads to a decrease of the volume shrinkage and ground subsidence.Specially, when the internal pressure is low (9 and 11 MPa), after 30 years of creep, the ground subsidence tends to its maximum, which would happen when the cavern totally closed. When the internal pressure is high (19 and 21 MPa), the ground subsidence tends to be the smallest, which is equivalent to the value after the end of the construction period. As a matter of fact, the total subsidence includes both the subsidence at the construction period and that at the gas storage period. Although changes of the internal gas pressure during operation may change the ground subsidence, the subsidence during the construction period should be constant. Moreover, the internal pressure in the horizontal cavern should be not too high, which may induce tensile fracture in the surrounding rock salts at the top of the cavern and even other parts. Therefore, the internal pressure should be raised appropriately.


### 4.2.3. The Cross-Sectional DiameterD

As it is known to all, it is much better to obtain a larger horizontal cavern for natural gas storage. If the length of the cavern is constant (200 m), the cross-sectional diameter of the horizontal cavern will determine the whole volume. Assuming the brine extraction velocity is kept 6.0 × 104 m3/month as a constant, if the cross-sectional diameters of the horizontal caverns are 20, 30, 40, 50, or 60 m, respectively, the cavern volume Vtp will be 6.28, 14.13, 25.12, 39.25, or 56.52 × 104 m3 by equation (10), and then, the building time will be, respectively, 9.6, 23.1, 45.4, 83.7, or 165.8 months. It can be seen that, as the diameter of the cross-section increases, the required time to build the cavern will become longer. In particular, it takes nearly 14 years to build a horizontal cavern with a diameter of 60 m and a length of 200 m. This is obviously uneconomical for the storage company, and thus, a reasonable choice of the cavern diameter is required.The reason why it needs such a long time is that the shrinkage of the cavern gradually increases with the increasing of the cavern volume. In fact, it should dissolve much more rock salt formations to form a horizontal cavern with diameter 60 m and length 200 m. Under a longer construction period, the salt cavern will shrink more because the brine pressure during the construction period is much lower than pressure during gas storage. Therefore, when the cross-sectional diameter increases proportionally, the required construction time will increase nonlinearly.Figure12(c) shows the impact of different cross-sectional diameters on the ground subsidence. From the figure, the ground subsidence increases obviously with the increasing of the cross-sectional diameter, which is disadvantageous to promote a larger volume capacity of the storage. Therefore, in the shape design period, an appropriate cross-sectional diameter should be selected.


### 4.2.4. The Length of the Horizontal CavernL

If the cross-sectional diameterD of the horizontal cavern is kept 40 m as a constant, the length L will determine the volume of the salt cavern. Assuming the brine extraction velocity is kept constant 6.0 × 104 m3/month, if the length of the cavern is set to be 100, 150, 200, 250, or 300 m respectively, correspondingly the cavern volume Vtp will be 12.56, 18.84, 25.12, 31.40, or 37.68 × 104 m3 by equation (10), and then, the building time will be 20.3, 32.0, 45.4, 60.8, or 78.5 months, respectively. It can also be seen that the required building time will become longer with the increasing of the length.Figure12(d) shows the impact of different lengths on the ground subsidence. From the figure, the ground subsidence increases obviously with the increasing of the cavern length. Thus, an appropriate length should be selected for the economic and safety considerations. In addition, according to equation (15), as the length increases, the ratio of the major axis to the minor axis of the subsidence trough (λ) becomes larger, which indicates that the changing of the length will also change the roundness of the subsidence contour.


## 4.2.1. The Brine Extraction VelocityVE˙

According to research [11], if the brine extraction velocity is faster in the construction period, the natural gas could be stored earlier. From an economic point of view, it is appropriate to increase the brine extraction rate. In fact, once the final cross-sectional diameter of the built cavern is determined, the required building time will become shorter with the increasing of the brine extraction velocity. If the prescribed cross-sectional diameter of the cavern is 40 m, assuming the brine extraction velocity is 6.0, 6.5, 7.0, 7.5, or 8.0 × 104 m3/month, correspondingly the building time will be 45.4, 41.6, 38.0, 35.0, or 32.4 months, respectively, by equation (10).Since the brine pressure in the cavern during the construction period is lower than the internal pressure during gas storage, the deviator stress on the surrounding salt is relatively large, and the cavern will inevitably shrink more. Therefore, increasing the speed of brine extraction can shorten the cavern construction period, thereby reducing the excessive shrinkage of the cavern and thus decreasing the excessive ground subsidence, as shown in Figure12(a).Figure 12
Impact of design parameters on the subsidence.
(a)(b)(c)(d)


## 4.2.2. The Internal Pressure in the Horizontal CavernP2

In order to study the impact of the internal pressure on the ground subsidence, the internal pressureP2 was kept at 9, 11, 13, 15, or 17 MPa as a constant in each case. The calculating results are as shown in Figure 12(b). It can be seen from the figure that the ground subsidence decreases nonlinearly with the increasing of the internal pressure. This reason is that the increasing of the internal pressure reduces the deviatoric stress of the rock salt, which leads to a decrease of the volume shrinkage and ground subsidence.Specially, when the internal pressure is low (9 and 11 MPa), after 30 years of creep, the ground subsidence tends to its maximum, which would happen when the cavern totally closed. When the internal pressure is high (19 and 21 MPa), the ground subsidence tends to be the smallest, which is equivalent to the value after the end of the construction period. As a matter of fact, the total subsidence includes both the subsidence at the construction period and that at the gas storage period. Although changes of the internal gas pressure during operation may change the ground subsidence, the subsidence during the construction period should be constant. Moreover, the internal pressure in the horizontal cavern should be not too high, which may induce tensile fracture in the surrounding rock salts at the top of the cavern and even other parts. Therefore, the internal pressure should be raised appropriately.


## 4.2.3. The Cross-Sectional DiameterD

As it is known to all, it is much better to obtain a larger horizontal cavern for natural gas storage. If the length of the cavern is constant (200 m), the cross-sectional diameter of the horizontal cavern will determine the whole volume. Assuming the brine extraction velocity is kept 6.0 × 104 m3/month as a constant, if the cross-sectional diameters of the horizontal caverns are 20, 30, 40, 50, or 60 m, respectively, the cavern volume Vtp will be 6.28, 14.13, 25.12, 39.25, or 56.52 × 104 m3 by equation (10), and then, the building time will be, respectively, 9.6, 23.1, 45.4, 83.7, or 165.8 months. It can be seen that, as the diameter of the cross-section increases, the required time to build the cavern will become longer. In particular, it takes nearly 14 years to build a horizontal cavern with a diameter of 60 m and a length of 200 m. This is obviously uneconomical for the storage company, and thus, a reasonable choice of the cavern diameter is required.The reason why it needs such a long time is that the shrinkage of the cavern gradually increases with the increasing of the cavern volume. In fact, it should dissolve much more rock salt formations to form a horizontal cavern with diameter 60 m and length 200 m. Under a longer construction period, the salt cavern will shrink more because the brine pressure during the construction period is much lower than pressure during gas storage. Therefore, when the cross-sectional diameter increases proportionally, the required construction time will increase nonlinearly.Figure12(c) shows the impact of different cross-sectional diameters on the ground subsidence. From the figure, the ground subsidence increases obviously with the increasing of the cross-sectional diameter, which is disadvantageous to promote a larger volume capacity of the storage. Therefore, in the shape design period, an appropriate cross-sectional diameter should be selected.


## 4.2.4. The Length of the Horizontal CavernL

If the cross-sectional diameterD of the horizontal cavern is kept 40 m as a constant, the length L will determine the volume of the salt cavern. Assuming the brine extraction velocity is kept constant 6.0 × 104 m3/month, if the length of the cavern is set to be 100, 150, 200, 250, or 300 m respectively, correspondingly the cavern volume Vtp will be 12.56, 18.84, 25.12, 31.40, or 37.68 × 104 m3 by equation (10), and then, the building time will be 20.3, 32.0, 45.4, 60.8, or 78.5 months, respectively. It can also be seen that the required building time will become longer with the increasing of the length.Figure12(d) shows the impact of different lengths on the ground subsidence. From the figure, the ground subsidence increases obviously with the increasing of the cavern length. Thus, an appropriate length should be selected for the economic and safety considerations. In addition, according to equation (15), as the length increases, the ratio of the major axis to the minor axis of the subsidence trough (λ) becomes larger, which indicates that the changing of the length will also change the roundness of the subsidence contour.


## 5. Discussion on the Advantages and Disadvantages of the Proposed Model

### 5.1. Comparison with the Traditional Pear-Shaped Cavern

Compared with pear-shaped salt cavern, horizontal cavern has its own characteristics. In terms of shape, horizontal cavern is more like a horizontal cylinder. According to equation (13) in this paper and equation (15) in the research [11], the unit volume shrinkage rate k is bigger than that of pear-shaped cavern under the same condition. Once the cavern reduction transmits to the surface, the ground subsidence above the horizontal cavern will be higher than that of the pear-shaped cavern.Generally, the volume of horizontal cavern is larger than that of pear-shaped cavern. If the brine extraction velocity is the same, it needs more time to construct a horizontal cavern than a pear-shaped cavern. Specially, if the brine extraction velocity is 6.0 × 104 m3/month, the required construction time will be about 165.8 months for a 60 m diameter horizontal cavern with length of 120 m. Consuming so long time is obviously uneconomical for storage companies. Therefore, the cross-sectional diameter and length of the horizontal cavern needs to be selected reasonably.Different from pear-shaped cavern, the contour of the ground subsidence above horizontal cavern is no longer a standard circle. According to equation (15), the shallower the horizontal cavern buried, the more elliptical the contour of the subsidence is. As the cavern buried deeper, the ratio λ of the major axis to the minor axis of the subsidence trough gets closer to 1, which indicates that the contour is getting closer to a circle. In addition, changing the length and cross-sectional diameter of the horizontal cavern will also change the roundness of the subsidence contour.As to the ground subsidence above horizontal cavern, the impact on the surface structures is different from that of the pear-shaped cavern, especially the calculation methods of the surface tilt, curvature, horizontal displacement, and horizontal strain, namely, the calculation formulas need to be adjusted appropriately.


### 5.2. Applicability of the Proposed Prediction Model

After all, by comparing to FLAC3D, the prediction results by the proposed model are acceptable. However, the real horizontal salt caverns are not idealized cylinder. The difference between the ideal cylinder and the actual horizontal cavern in the ground subsidence profile is inevitable. In fact, the shape of the impact function even changes over time [22]. Hence, this proposed model in this article is ideal and should be improved for practitioners to apply.The proposed prediction method is only suitable for the horizontal caverns with thicker rock salts at the roof and the floor. In this model, the shape of the impact function is thought to be not relying on the dimensions of the horizontal cavern; thus, the above proposed model is not applicable for horizontal caverns whose size is close to the cavern depth. Moreover, it is currently difficult to consider the impact of different stratigraphic layers. Since the ground subsidence above the horizontal cavern is very complex and time dependent, the above model should be verified by more field subsidence data and then be corrected accordingly.


## 5.1. Comparison with the Traditional Pear-Shaped Cavern

Compared with pear-shaped salt cavern, horizontal cavern has its own characteristics. In terms of shape, horizontal cavern is more like a horizontal cylinder. According to equation (13) in this paper and equation (15) in the research [11], the unit volume shrinkage rate k is bigger than that of pear-shaped cavern under the same condition. Once the cavern reduction transmits to the surface, the ground subsidence above the horizontal cavern will be higher than that of the pear-shaped cavern.Generally, the volume of horizontal cavern is larger than that of pear-shaped cavern. If the brine extraction velocity is the same, it needs more time to construct a horizontal cavern than a pear-shaped cavern. Specially, if the brine extraction velocity is 6.0 × 104 m3/month, the required construction time will be about 165.8 months for a 60 m diameter horizontal cavern with length of 120 m. Consuming so long time is obviously uneconomical for storage companies. Therefore, the cross-sectional diameter and length of the horizontal cavern needs to be selected reasonably.Different from pear-shaped cavern, the contour of the ground subsidence above horizontal cavern is no longer a standard circle. According to equation (15), the shallower the horizontal cavern buried, the more elliptical the contour of the subsidence is. As the cavern buried deeper, the ratio λ of the major axis to the minor axis of the subsidence trough gets closer to 1, which indicates that the contour is getting closer to a circle. In addition, changing the length and cross-sectional diameter of the horizontal cavern will also change the roundness of the subsidence contour.As to the ground subsidence above horizontal cavern, the impact on the surface structures is different from that of the pear-shaped cavern, especially the calculation methods of the surface tilt, curvature, horizontal displacement, and horizontal strain, namely, the calculation formulas need to be adjusted appropriately.


## 5.2. Applicability of the Proposed Prediction Model

After all, by comparing to FLAC3D, the prediction results by the proposed model are acceptable. However, the real horizontal salt caverns are not idealized cylinder. The difference between the ideal cylinder and the actual horizontal cavern in the ground subsidence profile is inevitable. In fact, the shape of the impact function even changes over time [22]. Hence, this proposed model in this article is ideal and should be improved for practitioners to apply.The proposed prediction method is only suitable for the horizontal caverns with thicker rock salts at the roof and the floor. In this model, the shape of the impact function is thought to be not relying on the dimensions of the horizontal cavern; thus, the above proposed model is not applicable for horizontal caverns whose size is close to the cavern depth. Moreover, it is currently difficult to consider the impact of different stratigraphic layers. Since the ground subsidence above the horizontal cavern is very complex and time dependent, the above model should be verified by more field subsidence data and then be corrected accordingly.


## 6. Summary and Conclusions

Due to a large demand of oil or gas storage in these years, horizontal caverns were proposed to make full use of bedded salt formations in China. In this paper, a new prediction model was put forward for time-dependent ground subsidence of the horizontal caverns. Using this proposed model, the impact of some main factors on the ground subsidence was calculated and analyzed, and the results show that(1)
Among all of the stratigraphic parameters, the draw angle has a great impact on the shape of the subsidence trough. In particular, the surface subsidence trough becomes deeper and narrower as the draw angle increases. The subsidence adjusting coefficient is direct proportional to the maximum subsidence. With the increasing of the cavern depth, the maximum subsidence first increases and then decreases, and the shape of the ground subsidence trough gradually becomes rounded. The amount of subsidence decreases gradually with the increasing of the purity of rock salts. With the increasing of the creep propertyCs and the stress index constant n, the maximum subsidence first increases rapidly and then slowly approaches its limit. The research studies on the impact of stratigraphic parameters on ground subsidence are helpful to the site selection of horizontal caverns.(2)
In terms of design parameters, increasing the brine extraction velocity can shorten the construction period of the cavern, thereby reducing the excessive shrinkage of the cavern volume and excessive ground subsidence. The amount of ground subsidence decreases nonlinearly with the increasing of the internal pressure. However, the internal pressure should be appropriately increased to avoid failure or damage at the cavern top. As the cross-sectional diameter and the length of the cavern increase, the construction time of the horizontal cavern will be extended nonlinearly, and accordingly, the subsidence will increase significantly. Although a larger volume capacity of the storage is better in economic aspects, the proper cross-sectional diameter and cavern length should be selected. The research conclusions on the design parameters are conducive to put forward the reasonable construction and operation plan.(3)
Compared with pear-shaped salt cavern, the unit volume shrinkage ratek of the horizontal cavern is bigger under the same condition, and thus, the volume of horizontal cavern is larger. Once the cavern reduction transmits to the surface, the ground subsidence above the horizontal cavern will be much higher than that of the pear-shaped cavern. Different from pear-shaped cavern, the contour of the ground subsidence above horizontal cavern is no longer a standard circle. The deeper the horizontal cavern buried, the rounder the contour of the subsidence will be.Above all, the proposed prediction method is only suitable for horizontal caverns with thicker rock salts at the roof and the floor. In this model, the shape of the impact function is thought to be not relying on the dimensions of the horizontal cavern; thus, the above-proposed model is not applicable for horizontal caverns whose size is close to the cavern depth. Moreover, it is currently difficult to consider the impact of different stratigraphic layers. Since the ground subsidence above the horizontal cavern is especially complex and time dependent, the above model should be verified by more field subsidence data and then be corrected accordingly.



---

*Source: 9504289-2021-09-27.xml*