Surface Deformations Caused by the Convergence of Large Underground Gas Storage Facilities

: The article presents a method of forecasting the deformation of the land surface over large ﬁelds of underground gas storage facilities located in salt caverns. The solution allows for taking into account many parameters characterising the operation of underground gas storage facilities, such as cavern processes (leaching, enlargement, operational, etc.), their depth, distribution, diameter, shape, and many others. The advantage of the applied method over other available options is the possibility of using it for large ﬁelds of caverns while keeping the calculations simple. The effectiveness of the method has been proven for predicted surface subsidence for the EPE ﬁeld with 114 underground caverns. The hypothesis was compared with the measurement outcomes.


Introduction
In recent decades, there has been an increase in the demand for natural gas in the world [1]. This has been noticed by the European Commission Directorate-General for Energy, which defined energy security as the ability of the gas transmission system to ensure continuous and reliable gas supplies to consumers on economic terms and the ability to face interruptions in natural gas supplies. This definition indicates that the actions of the EU state bodies should aim both at ensuring operational security (current at the time), tactical security (planned and predicted for a specific season) and strategic security (planned and predicted in subsequent years) [2].
Considering the fluctuations of natural gas consumption in an annual period, such as the increase in the autumn-winter period and subsequent decrease in the spring-summer period (Figure 1), operators have been forced to create systems to accumulate gas surpluses. This is partly because energy storage can play a pivotal role in future energy systems compatible with a carbon-neutral and environmentally friendly society.
In general, there are four types of underground gas storage facilities contained in: depleted reservoirs, aquifers, mines and salt caverns. Each type has its own physical characteristics such as; retention capability, porosity, permeability, and economic issues including site preparation and maintenance costs, deliverability rates, and cycling capability [4,5]. The most common natural gas storage type globally is depleted reservoirs because their greatest advantages are their wide availability and existing underground and surface infrastructure (existing wells, gathering systems, and pipeline connections) [6][7][8]. However, it is the gas storage facilities built in salt dome formations (salt caverns) that are currently building popularity. These storage facilities are characterised by significantly higher gas withdrawal rates than other types of UGS and may serve as peak gas storage facilities (withdrawal and injection cycles). Currently, 9% of underground gas storage facilities worldwide are located in salt caverns. Table 1 presents the number of salt caverns fields in Europe.  These fields consist of a variety of caverns. For example, Solvay has 114 caverns in the EPE field (in Gronau-Epe-Germany), 31 production, 5 crude oil, 77 natural gas and 1 helium (as of 09/2019). However, the unquestionable advantage of gas (and oil) storage in salt caverns, i.e., their availability to cover short-term, very large gas deficits, has also become a problem for them. UGS operators in caverns, given their benefits, intensified the work of the caverns in terms of market instruments by creating a network of commercial storage of gas and oil. These activities led to frequent annual changes in pressure inside the caverns. The related injection, unloading and additional leaching caused a significant disruption of the cavern convergence process, including their significant acceleration [9].
The problem of the progressive accelerated convergence of salt caverns is currently very important for UGS operators. The correct method of forecasting the convergence and the related change in the state of rock mass deformation allows for the prevention of damage to the operating pipelines, damage to the surface transmission network, and the determining of the required period of leaching processes-i.e., restoring or increasing the storage volume, or specifying the periods of required echometric or geodetic measurements [10]. For example, for the mentioned EPE cavern field, the maximum measured subsidence of the land surface in 2013 slightly exceeded 700 mm [11].
The presented article shows a modern and effective method of predicting the deformation of the rock mass and the land surface above the used large cavern field. Presented calculation results were made at the end of 2019. Then, the forecasted deformation indicators were confirmed by geodesic measurements carried out by the cavern field operator (SWG Solvey Epe, Ahaus, Germany). The effectiveness of the forecast depends on the details of the information provided (including, among others, cavern dimensions, production phases, convergence rate), which allows for the advance determining of the potential environmental hazards associated with possible damage to the operating pipelines and leakage of the medium. This knowledge permits for preventive work to be carried out and thus the safety of the transmission of energy medium.

Methods for Determining Changes in Rock Mass Deformation for Salt Caverns
Many scientists have dealt with the influence of underground gas storage facilities located in salt caverns on the stress and deformation of the rock mass, its tightness and surface area [6,[12][13][14][15]. This is related to the growing interest in this method of storing energy factors. Depending on the challenges faced by scientists, they used various tools, most suitable for obtaining satisfactory results. The literature contains many works in which the authors used numerical models to describe the behaviour of a single salt cavern. These analyses were mainly related to the problem of heterogeneity in the rock structure of the underground storage and the assessment of their tightness, among others [16]. The authors used various physical models for the description of the rock mass and the processes occurring in it; for example, in the work of [17] the authors have used Drucker-Prager constitutive law with equivalent Mohr-Coulomb parameters, including a tension cut-off with the WIPP creep model to simulate plastic yielding and creep. These calculations allowed for the determination of the deformations in the area of a single salt cavern with mudstone interbeds. A similar approach was presented by [18] in their work to describe the behaviour of a single horizontal cavern. They used the creep model with Mohr-Coulomb criterion. Wang et al. [13] carried out numerical analyses for the two adjacent caverns and their mutual influence on stability. Numerical analyses for determining the stability of the pillar separating the two adjacent caverns were also conducted [19] however, numerical methods were not the only engineering tool for assessing the behaviour of underground gas storage facilities [20]. The growing interest in neural networks encouraged scientists to use them also for UGS analyses [21]. Despite the popularity and significant advantages of the methods presented above, analytical methods have also found their place in the description of the changing state of stress and strain in the areas of underground salt caverns [22]. In the work of [23], the authors presented interesting analytical and numerical solutions for the hexagonal system of seven underground caverns. The authors analysed the optimal distribution of these caverns and their influence on the surrounding rock mass. Polański [24] carried out analyses of changes in the strain-stressing state of the caverns during thermal change, which is observed during gas injection and cavern emptying processes. The numerical analysis carried out for an exemplary cavern modelled by Norton's rheological model [15,[24][25][26] showed that the thermal change implemented in the solutions affects the state of the cavern in a small range and should be taken into account when analysing the area up to 50m from the cavern axis in detail (for caverns with a diameter of 30-35 m).
In previous work [27], the authors presented analytical solutions enabling the calculation of the volume of a single cavern, its capacity and convergence. Analytical methods are used especially when conducting analyses for large areas of cavern fields and when it comes to determining potential deformations of the ground surface [28,29]. The advantage of analytical solutions is the simplicity of building a model of the entire cavern field in the salt dome and taking into account their mutual relations and nature. However, this solution also has disadvantages. They are associated with one-parameter describing the nature of the rock mass. This assumption causes an inability to analyse rock mass stratification.
This deformation significantly affects the safety of pipelines and gas pipelines [30,31]. Examples of damage can be found in the works of, among others [32][33][34][35].

Method of Surface Deformation Prediction for a Large Field of Salt Caverns
As shown in the previous section, numerical modelling of cavern fields usually involves several caverns. This limitation is related to the number of elements and the calculation of equilibrium equations that a computational machine must perform in order to obtain a result. Therefore, in order to determine the deformation state of the rock mass and terrain surface for the whole cavern field covering more than 100 salt caverns, the authors used their own analytical solutions ( Figure 2). They were based on the modified geometrical-integral condition of Knothe theory [36]. The modification was based on the basic assumptions presented by Sroka and Schober [37]. The size of the deposit element was determined in a way, so that the trough formed by the extraction of this element was practically matching the fundamental solution in the form of a Gaussian function and it allowed for, among other things, precise consideration of the three-dimensional geometry of the selected deposit, i.e., the variation in its thickness and depth and of other necessary information such as the distribution of pressure in the case of fluid deposit exploitation. That is why, based on presented research, simulations, and parameter identifications, a new methodology (modification) has been developed. The theoretical solution of modification was presented in the conference [29], where: and: S(r, z, t)-subsidence of surface point at moment t, situated in the ground level z at a distance of r from the cavern axis, s max (z, t)-maximum subsidence at moment t, H r -distance between the terrain surface and the depth of cavern roof, H f -distance between the terrain surface and the depth of cavern floor, R r -parameter of the horizontal influence scale, so-called radius of main influences, calculated from the cavern roof, R f -radius of main influences, calculated from the cavern floor, h = H f -H r -cavern's height, β-so-called angle of main influences (Knothe 1953), λ-coefficient of horizontal displacement, c-coefficient/characteristics of cavern roof, n-surface factor of the range of main influences of exploitation in the rock mass, z-height on the terrain surface, z r -depth of cavern roof occurrence, z f -depth of cavern floor occurrence.
But it needed verifications. The relationship between the maximum subsidence S max (t) and the subsidence trough volume M(t) can be presented as: where: R r (z)-radius of main influences, calculated at any z point in the strata from the cavern roof, and R f (z)-radius of main influences, calculated at any z point in the strata from the cavern floor.
In the case of a cavern field, the calculation of any surface point subsidence is carried out while assuming a linear superposition, i.e., the summation of subsidence from individual caverns. The subsidence trough volume M(t) depends on the volumetric convergence K(t), the delaying effect of the overlying rock mass, and possible volume losses connected with the deformation of the overlying rock mass. The convergence in time may be described analytically using a logarithmic function or an exponential function [38].
Assuming an incremental leaching model, Sroka and Schober [37] ascertained that the subsidence trough volume on the surface may be described by the following formula: The analyses of the subsidence results carried out in situ indicated that the a factor was practically equal to one, i.e., no volume losses occured in the rock mass [39].
The S max (t) value also depends on the shape of the cavern (e.g., is it a cylinder, a sphere, a cone) and on the geometrical convergence model [39,40]. The results of the comparative calculations carried out lead to ascertaining that the maximum subsidence value may be determined using the formula (6) with a very good approximation: where: K(t−∆t)-volumetric convergence of cavern over time t−∆t, ∆t = 1 f -delay time caused by the delaying action of the overlying rock mass. The basic component of the computational model consists of the correct definition of the initial condition in the form of volumetric convergence of a single cavern [41].
The history of the operation of single caverns indicates that the leaching or enlargement phases may change multiple times [42]; also, the same pertains to the stored medium (petroleum or gas). In this connection, the multiphase cavern use model was adopted for the description of the convergence course in time, introducing the phases of leaching, standby (preparedness), operation and enlargement (Figure 3). It is obvious that each phase leads to different convergence values, and thus to different values of the ξ parameter. Examples of average values of the relative convergence rate ξ depending on the production phase for the EPE cavern field are reported in Table 2 [11]. When adopting the exponential model of the volumetric convergence process, the convergence increase for the i-th phase of use in this phase ( Figure 3) may be calculated using the formula: Thus, the cavern volume at the end of the i-th phase of use equals: For the leaching phase, assuming the linear model of the leaching process, the convergence value may be calculated acc. to the Formula (9) for any time point during the leaching phase (with the assumption that t 0 < t < t 1 and T = t 1 −t 0 ): where: T-leaching time, ξ 0 -relative volumetric convergence rate. According to various literature sources [43], the ξ 0 values are in the range of: For a definite calculation of time t, convergence values for the operational phases should be added up between the beginning of the leaching t 0 and the t−∆t time.
The presented solution allows for the calculations to be carried out for every phase of cavern operation; among others, the following values may be calculated: convergence and free volume, subsidence, tilt, curvature, horizontal displacement, and deformation. The calculations may be carried out for any moment in time, not only for the final period of a given phase.
This solution has been successfully applied to many large cavern fields in Germany. In the next chapter, the authors present a sample case study.

Prediction of Land Surface Deformation Along with Its Verification for an Exemplary Area of an Underground Gas Storage Facility
The Epe underground gas storage facilities are built in the form of caverns leached in the salt dome. Its size is estimated at a cumulative capacity volume of 42 million m 3 (114 salt caverns up to 2019) (Figure 4). Using the presented algorithm, based on the data provided by the cavern operator SGW (Salzgewinnungsgesellschaft Westfalen mbH & Co. KG, Ahaus, Germany) the authors made forecast calculations of the future land surface deformation for the period to 12.30.2019.
All 114 caverns were included in the calculation and their geometries and the course of the phases of use in time. The data pertains to the time period starting in 1972. Exemplary data attributed to the cavern no. S_004, characterising its individual production phases, are presented in Table 3. The calculations were carried out for the following values of characteristic parameters of the EPE cavern field: angle of main influence-β = 34 • , subsidence coefficient-a = 1.0, and delay time ∆t = 4 months. As a result of the calculations, the forecasted values of deformation indicators for the period of 12.30.2019 were obtained. The calculated maximum subsidence value for an exemplary field of 114 caverns equaled 81 cm, while the measured value amounted to 79 cm. In Figure 5, the course of the levelling point subsidence no. 211024 located above the cavern field in the period of 1975-2019 is presented as the comparison of theoretical course and geodetic measurements. Geodetic measurements have been carried out since 1973. Approximately 850 points were levelled each year. The levelling campaign was over 130 km. In addition, 10 fixed points were located outside the subsidence trough.  In Figure 7 results of theoretical calculations of horizontal displacement in mm over the EPE cavern field-as per 12.30.2019 are shown with three measuring points (u1, u2 and u3) where horizontal displacements were measured. In Table 4, comparison of calculated and measured data for 3 example points for horizontal displacements for three different time periods (2006-2019, 2003-2016 and 2000-2016) with their differences in cm are presented. The data shown in Figures 5-7 and Table 4 confirm that the mathematical model proposed in this paper is able to describe the process of subsidence of the ground surface over a large cavern field in a salt rock mass relatively well and sufficiently accurately.

Conclusions
Recently the prediction of UGS facilities inducted surface problem is becoming increasingly important every year. Changes in the use of underground gas storage and frequent fluctuations of the pressure inside the caverns cause stress-strain imbalance. This creates an increase in the velocity of convergence of surrounding salt rock mass. The consequence of this behaviour is the failure of salt caverns or damage in pipelines (operating and transport).
The application of the presented method allows for the determination of future deformations of the ground surface, taking into account the processes in the cavern (leaching, enlargement, operational, etc.). The results of these calculations can be easily applied to analysing changes in the stress of underground gas storage infrastructure objects such as pipelines or other construction objects.
Despite necessary geometrical and physical idealisations related to the cavern geometry, course of convergence, and phases of operation, comparative computations performed for the EPE cavern field with its 114 caverns, have fully confirmed the value and the integrity of the presented solution.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to company's policy.

Conflicts of Interest:
The authors declare no conflict of interest.