Shape Modelling and Volume Optimisation of Salt Caverns for Energy Storage

Salt caverns are an attractive solution to the growing energy demand in view of their large storage capacity, safety of storage operation and long operation time. The designing process of salt caverns is still considered a complex issue despite progress in geotechnical, construction and exploration methods. Finding the optimal shape and dimensions of a salt cavern in given geological conditions is a difficult engineering problem in view of safety and stability requirements. In this paper, the stability of typical cavern shapes (cylindrical, enlarged top, and enlarged bottom), with each of the three variants differing by their diameter, was evaluated against the stability factors of the geological conditions of the bedded salt deposit. Moreover, the analysed shapes were examined in terms of edges. The three-step smoothing of sharp edges was performed, and its impact on the cavern’s stability performance was studied. Moreover, the analysis aimed to find the optimal cavern shape and volume in the implemented geological conditions. The evaluation was based on the following criteria: the displacement, effective strain, von Mises stress, strength/stress ratio and safety factor. The results of this evaluation can be useful in the design of an optimal cavern shape and volume and for planning new cavern fields for storing natural gas, compressed air or hydrogen in the bedded salt deposits.


Introduction
Salt caverns are the best method for the storage of natural gas due to their large capacity, safety and long operation time [1,2]. Most previous studies concerning underground gas storage facilities have been limited to pure and homogeneous salt formations in salt domes. However, in the last two decades, increasing attention has been paid to gas storage in bedded salt formations [3][4][5][6]. In recent years, salt caverns have been considered for large-scale hydrogen and compressed air storage [7][8][9][10].
The pressure limits of salt caverns in bedded salt deposits were studied by Bruno and Dusseault [11]. These authors described cavern deformation and bedding plane slips for a variety of cavern configurations. The optimal operating pressure of salt caverns in a bedded salt formation with high insoluble content was analysed by Zhang et al. [12] and Wang et al. [13]. The cavern roof stability in a bedded rock salt formation was studied by Bruno [14], Bruno et al. [15] and deVries et al. [16]. The influence of interlayers on storage operations was analysed by Cosenza et al. [17] and Consenza and Ghoreychi [18]. Weakness planes and hydraulic fracturing in bedded salt formation were reported by Minkley and Muhlbauer [19]. The deformation characteristics of mudstone interbeds and their impact on cavern stability was described by Yang et al. [20], Li et al. [5], Wang et al. [3] and Yu and Liu [4]. The influence of anhydrite interbeds on cavern stability was studied by Cała et al. [21]. Horizontal caverns in bedded salt deposits were studied by Wei et al. [22] and Jie et al. [23].

Experimental Verification of Cavern Shapes and Dimensions
The evaluation of cavern shapes was based on numerical modelling and divided into two parts: the first part included the evaluation of three standard cavern shapes: cylindrical, enlarged top and enlarged bottom. The shape of enlarged caverns was based on a cylinder whose diameter in the roof of the bottom area was increased. The dimensions of examined caverns were adapted to the geological conditions of the Mechelinki salt deposit. Due to the fact that, in numerical simulations, sharp edges are prone to stress concentration [35,36], the "smoothing" of sharp edges was performed in the second part. The "smoothing" was performed in three steps. In each step, a quadrant characterised by a different radius (10-30 m) was inscribed into the corners of a vertical cross section through the analysed cavern shapes ( Figure 1). Consequently, the cavern shape evolved from a cylinder to a more complex shape. In the first step, the cavern shapes based on a cylinder with sharp edges were examined. In the next three steps, the edges were smoothed, and the radiuses of the quadrant were 10 (Step 1), 20 (Step 2), and 30 m (Step 3). Thus, four cases were considered: caverns with sharp edges (se) and smoothing steps 1 (ss1), 2 (ss2), and 3 (ss3). Finally, the analysed cavern shape consisted of a cylinder in the middle part and two hemispheres in the upper and lower part. As a result of smoothing, the maximal initial diameter that caverns had in the first step (Table 1) decreased, as did the cavern volume. The volume of the cavern from smoothing to Step 3 was about 84% of the cavern volume with sharp edges.
Three shapes (cylindrical, enlarged top and enlarged bottom) varying in their maximum initial diameter with sharp edges and three "smoothing" steps for each shape (28 caverns in total) were analysed. The maximum initial diameter of the cylindrical shape was 60 m, and the height was 120 m. This cylindrical shape was validated as a base shape for further simulations. The base diameter and the height of the cavern were chosen based on the salt bed thickness and the most optimal proportions (H/D = 2/1). The base cavern diameter was enlarged in a roof or bottom of 60 m, 30 m and 15 m. As a result, caverns with an enlarged top or enlarged bottom and an initial diameter of 120 m, 90 m and 75 m were created ( Table 1). The cavern height was constant for all simulated cavern shapes in consideration of the thickness of the rock salt layers in the Mechelinki salt deposit. The roof and bottom pillar thicknesses were set to 15 and 5 m, respectively. Therefore, the cavern diameter and edges were the key factors investigated. The volume of each simulated cavern was calculated (Table 1). All analysed shapes were regular in comparison to the complex shapes of real caverns or underground workings [21,31]. These simulated shapes can form a basis for future cavern design.

Experimental Verification of Cavern Shapes and Dimensions
The evaluation of cavern shapes was based on numerical modelling and divided into two parts: the first part included the evaluation of three standard cavern shapes: cylindrical, enlarged top and enlarged bottom. The shape of enlarged caverns was based on a cylinder whose diameter in the roof of the bottom area was increased. The dimensions of examined caverns were adapted to the geological conditions of the Mechelinki salt deposit. Due to the fact that, in numerical simulations, sharp edges are prone to stress concentration [35,36], the "smoothing" of sharp edges was performed in the second part. The "smoothing" was performed in three steps. In each step, a quadrant characterised by a different radius (10-30 m) was inscribed into the corners of a vertical cross section through the analysed cavern shapes ( Figure 1). Consequently, the cavern shape evolved from a cylinder to a more complex shape. In the first step, the cavern shapes based on a cylinder with sharp edges were examined. In the next three steps, the edges were smoothed, and the radiuses of the quadrant were 10 (Step 1), 20 (Step 2), and 30 m (Step 3). Thus, four cases were considered: caverns with sharp edges (se) and smoothing steps 1 (ss1), 2 (ss2), and 3 (ss3). Finally, the analysed cavern shape consisted of a cylinder in the middle part and two hemispheres in the upper and lower part. As a result of smoothing, the maximal initial diameter that caverns had in the first step (Table 1) decreased, as did the cavern volume. The volume of the cavern from smoothing to Step 3 was about 84% of the cavern volume with sharp edges.

Constitutive Models and Mechanical Parameters
The structural stability of caverns in the bedded salt deposit depends on the strength and deformation characteristics of the salt and nonsalt beds surrounding and overlying the avern [16].
The geomechanical model was built based on the geological characteristics of the Mechelinki salt deposit in order to study the response of the rock mass around the cavern under the bedded rock salt conditions. The salt beds are part of the Zechstein salt formation (cyclothem PZ1). In the geomechanical model, the rock salt beds were located at a depth of 960-980 m below ground level (b.g.l.). The thickness of the salt layers ranged from 160 to 190 m ( Figure 2). The salt beds are underlain and overlain by anhydrite layers (the lower anhydrite being A1d and the main anhydrite being A1g) with thicknesses ranging from about 1.5 to 20.0 m. There are two anhydrite interlayers in the rock salt beds with a thickness of up to 0.5 m [37]. These interlayers were not considered in the simulations described in this paper because their stability performance would obscure the impact of shape and smoothing on the cavern's stability. Moreover, the influences of anhydrite interlayers on the cavern's stability were analysed in a previous paper [21]. The salt beds are underlain by Kupferschiefer (T1), Zechstein Limestone (Ca1) and Rotliegend and Silurian sediments and covered by the sediments of the PZ2 and PZ3 cyclothems as well as Triassic, Jurassic, Cretaceous, and Kenozoic strata [37].
The mechanical behaviour of the rock mass surrounding the examined caverns was simulated with the use of constitutive models and the mechanical parameters of the rock salt and the nonsalt rock. The Mohr-Coulomb elastic-plastic model was applied to describe the mechanical behaviour of the surrounding rocks and the anhydrite; furthermore, a twocomponent Norton Power Law with a Mohr-Coulomb plasticity criterion was used to simulate plastic yielding and creep. These models were chosen to accurately simulate the elastic-plastic response of anhydrite and the viscoelastic plastic behaviour of rock salt.
The mechanical parameters of rock salt and nonsalt rock used in the numerical simulations are given in Table 2. The presented parameters were determined based on laboratory tests and validated in previous studies [21]. Two types of materials were distinguished, as required for the numerical modelling: rock salt and anhydrite. The parameters applied to the numerical calculations were based on the results of laboratory tests (Table 2). In addition, creep parameters for rock salt (determined in laboratory tests and calculated based on the Norton Power law) were n = 5.0 and A = 1.08·10 −45 Pa −5.0 s −1 .
Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 27 Figure 2. The geological structure of the Mechelinki salt deposit in the area of the underground store [21].
The mechanical behaviour of the rock mass surrounding the examined caverns was simulated with the use of constitutive models and the mechanical parameters of the rock salt and the nonsalt rock. The Mohr-Coulomb elastic-plastic model was applied to describe the mechanical behaviour of the surrounding rocks and the anhydrite; furthermore, a two-component Norton Power Law with a Mohr-Coulomb plasticity criterion was used to simulate plastic yielding and creep. These models were chosen to accurately simulate the elastic-plastic response of anhydrite and the viscoelastic plastic behaviour of rock salt.
The mechanical parameters of rock salt and nonsalt rock used in the numerical simulations are given in Table 2. The presented parameters were determined based on laboratory tests and validated in previous studies [21]. Two types of materials were distinguished, as required for the numerical modelling: rock salt and anhydrite. The parameters applied to the numerical calculations were based on the results of laboratory tests ( Table 2). In addition, creep parameters for rock salt (determined in laboratory tests and calculated based on the Norton Power law) were n = 5.0 and A = 1.08•10 −45 Pa −5.0 s −1 .

Numerical Model and Boundary Conditions
The numerical simulations were performed for a set of individual caverns (cylindrical with a maximum initial diameter of 60 m and an enlarged bottom and top with maximum initial diameters of 75 m, 90 m and 120 m (Figure 3)). Boundary conditions were adopted in the form of stress conditions to obtain solutions that tended towards being safe. The model was a cylinder with a diameter of 250 m and a height of 320 m. The direct roof and the bottom of the rock salt beds consisted of anhydrite layers with a thickness of 35 m at the roof and 80 m at the bottom. The mesh contained about 1.65 million tetrahedral elements. To ensure that the cavern shapes were accurately represented, the element size at the cavern sidewalls was smaller than others (about 3 m). Outside the cavern, the sidewall elements were larger than others (about 20 m). The initial value of the hydrostatic stress changed within a depth from 22.32 MPa at the top to 30.00 MPa at the bottom of the 3D model.

Numerical Model and Boundary Conditions
The numerical simulations were performed for a set of individual caverns (cylindrical with a maximum initial diameter of 60 m and an enlarged bottom and top with maximum initial diameters of 75 m, 90 m and 120 m (Figure 3)). Boundary conditions were adopted in the form of stress conditions to obtain solutions that tended towards being safe. The model was a cylinder with a diameter of 250 m and a height of 320 m. The direct roof and the bottom of the rock salt beds consisted of anhydrite layers with a thickness of 35 m at the roof and 80 m at the bottom. The mesh contained about 1.65 million tetrahedral elements. To ensure that the cavern shapes were accurately represented, the element size at the cavern sidewalls was smaller than others (about 3 m). Outside the cavern, the sidewall elements were larger than others (about 20 m). The initial value of the hydrostatic stress changed within a depth from 22.32 MPa at the top to 30.00 MPa at the bottom of the 3D model.

Method of Analysis and Assessment Criteria
The evaluation of the cavern shapes and dimensions was based on numerical modelling. The 3D geomechanical model was introduced into FLAC 3D -software for solving geomechanical and geotechnical problems, including rheological phenomena.

Method of Analysis and Assessment Criteria
The evaluation of the cavern shapes and dimensions was based on numerical modelling. The 3D geomechanical model was introduced into FLAC 3D -software for solving geomechanical and geotechnical problems, including rheological phenomena.
The effect of the cavern shape and dimensions on cavern stability and volumetric closure was examined by the analysis of the following criteria: the displacements, effective strains (ESs), von Mises stress (vMS), strength/stress ratio (SSR) and safety factor (SF). The same criteria were used in the optimisation of the cavern shape and volume. Displacements were analysed at the cavern sidewalls as well as at the bottom and the roof of the caverns. The displacements also illustrated the decrease in cavern volume with time caused by the salt creep (volumetric closure). The vMS is the equivalent effective stress that is related to the creep rate within the primary and the secondary creep stages. ES Equation (1) is a scalar quantity that represents the magnitude of the plastic strain tensor. The ES is expressed by where ε ef is the effective strain and ε dev is the deviatoric strain tensor. The effective strain represents damage associated with plastic deformation. The ES of the rock mass surrounding the cavern should be less than 3% during the period of operation, which is usually 30 years [29,34]. The strength/stress ratio (SSR) is used to indicate the dangerous areas in the sidewalls of the caverns and in its surroundings. An SSR of 1.0 implies a material failure, and an SSR of 2 indicates that the material reaches 50% of its strength. The safety factor (SF) is a dilatancy damage criterion Equation (2) for rock salt that can be expressed by where b is a material constant, I 1 is the first invariant of the stress tensor and I 2 is the second invariant of the deviatoric stress tensor. According to von Sambeek et al. [38], de Vries et al. [39,40] and Sobolik and Ehgartner [26], an SF < 1 indicates failure, and an SF < 1.5 shows local damage. In addition, the volumetric closure (volume convergence) of each cavern shape was evaluated. The volumetric closure was calculated as the ratio of the volume loss to the initial cavern volume. The accepted convergence rate was below 1% per year [41].
The numerical analysis was performed for a period of 9.5 years and consisted of three phases: leaching, 4 years when the caverns were filled with brine, and five operation cycles in which the cavern pressure oscillated between the minimum of 4 MPa to the maximum of 17.5 MPa. For each operation cycle, the aforementioned criteria were analysed in two periods: at the end of the gas injection period, and at the end of the gas withdrawal period. The cavern leaching was simulated to be excavated instantly in the applied software. This assumption resulted from the fact that the displacements in the leaching period were very small in comparison to the operation period. Consequently, their role in the long-term cavern stability was neglected. To eliminate the adverse effect of this assumption on the results, the in situ stress redistribution was calculated first and used as the initial conditions for further calculations including creep.

Results
The cavern shapes and dimensions were validated based on the factors described above. The results of the stability analysis are presented in the form of maps in  To facilitate reading, the description of caverns includes their shape, initial maximum diameter and smoothing step. The maps show, with the help of a colour scale, the values of the analysed factors: the displacements, ES, vMS, SSR and SF in the relevant time-i.e., the initial period (after leaching), the first operation cycle (the first gas withdrawal and injection) and the last operation cycle (the last gas withdrawal and injection). Different colours represent the values of each factor; for instance, the dark blue colour represents a value of 3 for the SSR in the maps. The "smoothing" steps are marked as "s1", "s2" and "s3", and the sharp edges are marked as "se" in the text below. The figures in the text illustrate the most significant changes in the value and distribution of simulated factors. The full results of the analysis are presented in five supplements (marked from A to E) that are available in the Supplementary Materials.

Displacements
In the initial period, the largest displacements (23.8 cm) were indicated for the roof (red area) of the caverns with an enlarged top and diameter of 120 m ( Figure 4A). The gradual "smoothing" of the sharp edges caused a reduction of the largest displacements in this area. The decrease in the cavern diameter corresponded to a decrease in the largest displacement value and area ( Figure 4E-G). The largest displacements in the walls (14.0-16.0 cm) were found in the cavern with an enlarged bottom, sharp edges and an initial diameter of 120 m ( Figure 4E). In the caverns with maximum initial diameters of 90 and 75 m, the area of the largest displacements was smaller and concentrated in the middle part of the caverns ( Figure 4F,G).  At the end of the first withdrawal period, the caverns with an enlarged bottom, a maximum initial diameter of 120 m, sharp edges and ss1 showed the largest displacements (138.44 cm) in the middle part and the surroundings ( Figure 5A,B). The value of displacements reached 120.0-130 cm in the caverns with an enlarged top, maximum initial diameter of 120 m, sharp edges and ss1 for the displacements located over the roof of these caverns ( Figure 5E,F). In the middle part of the caverns with an enlarged top, se and ss1 as well as with an enlarged bottom, ss2 and ss3, the largest displacements decreased to 100.0-110.0 cm ( Figure 5C-F). A further decrease to 90.0-100.0 cm was found in caverns with an enlarged top, ss2 and ss3 ( Figure 5G,H). The largest displacements in the caverns with an initial diameter of 90 and 75 m and sharp edges were 70.0 and 50.0 cm, respectively, and were concentrated in the middle part of these caverns ( Figure 5I,J). The caverns characterised by ss3 ( Figure 5K,L) showed a lower value in the largest displacements from 70.0 cm (diameter: 90 m) to 60.0 cm (diameter: 75 m). At the end of the first withdrawal period, the caverns with an enlarged bottom, a maximum initial diameter of 120 m, sharp edges and ss1 showed the largest displacements (138.44 cm) in the middle part and the surroundings ( Figure 5A,B). The value of displacements reached 120.0-130 cm in the caverns with an enlarged top, maximum initial diameter of 120 m, sharp edges and ss1 for the displacements located over the roof of these caverns ( Figure 5E At the end of the first injection period, the distribution of the largest displacements within the walls of the caverns was similar to the withdrawal period, but their value was lower ( Figures 5 and 6). The largest displacements were registered in the walls of the cavern with an enlarged bottom (128.30 cm) and over the top of the cavern with an enlarged top (110.0 cm), both of which were characterised by an initial diameter of 120 m and se ( Figure 6A  At the end of the first injection period, the distribution of the largest displacements within the walls of the caverns was similar to the withdrawal period, but their value was lower (Figures 5 and 6). The largest displacements were registered in the walls of the cavern with an enlarged bottom (128.30 cm) and over the top of the cavern with an enlarged top (110.0 cm), both of which were characterised by an initial diameter of 120 m and se ( Figure 6A  At the end of the last injection period, the tendency towards displacement distribution related to the "smoothing" of sharp edges was similar to the first period. However, this tendency was more visible in the caverns with an enlarged top than in the caverns At the end of the last injection period, the tendency towards displacement distribution related to the "smoothing" of sharp edges was similar to the first period. However, this tendency was more visible in the caverns with an enlarged top than in the caverns with an enlarged bottom ( Figure 8). The value of the largest displacements was lower than in the withdrawal period. The largest displacements in the caverns with se and a maximum initial diameter of 120 m ( Figure 8A,B) reached 359.86 cm (enlarged bottom) and 350.00 cm (enlarged top). In the caverns with se and a maximum initial diameter of 75 m, the value of the displacements ranged from 125.0 to 175.0 cm ( Figure 8E,F). The smoothing to ss3 caused a decrease of the displacement value to 75.0-125.0 cm and a major change in their distribution (Figure 8C,D,G,H).
The performance of the cylindrical caverns in displacement simulations indicated that the "smoothing" of sharp edges had a less significant effect on displacement value and distribution than in caverns with enlarged shapes (Figures 5-8). Moreover, the displacements in the cylindrical caverns were the lowest of all analysed shapes. The largest value of displacements changed from 10.0 cm in the initial period to 40.0-60.0 cm in the first cycle ( Figure 9A-D) and to 75.0-100.0 cm in the last cycle ( Figure 9E-H). There were visible differences in the displacement value and distribution between withdrawal and injection periods ( Figure 9).
The largest value and extension of displacements in all analysed caverns occurred as the gas pressure decreased (withdrawal period). An increase in gas pressure (injection period) caused a decrease in displacement value and distribution. The largest displacements, found in caverns with an enlarged bottom and a diameter of 120 m, were very high and covered a wide area in the side walls and surroundings ( Figure 5A At the end of the last injection period, the tendency towards displacement distribution related to the "smoothing" of sharp edges was similar to the first period. However, this tendency was more visible in the caverns with an enlarged top than in the caverns with an enlarged bottom ( Figure 8). The value of the largest displacements was lower than in the withdrawal period. The largest displacements in the caverns with se and a maximum initial diameter of 120 m ( Figure 8A   that the "smoothing" of sharp edges had a less significant effect on displacement value and distribution than in caverns with enlarged shapes (Figures 5-8). Moreover, the displacements in the cylindrical caverns were the lowest of all analysed shapes. The largest value of displacements changed from 10.0 cm in the initial period to 40.0-60.0 cm in the first cycle ( Figure 9A-D) and to 75.0-100.0 cm in the last cycle ( Figure 9E-H). There were visible differences in the displacement value and distribution between withdrawal and injection periods (Figure 9). The largest value and extension of displacements in all analysed caverns occurred as the gas pressure decreased (withdrawal period). An increase in gas pressure (injection period) caused a decrease in displacement value and distribution. The largest displacements, found in caverns with an enlarged bottom and a diameter of 120 m, were very high and covered a wide area in the side walls and surroundings ( Figures 5A,B

Effective Strains (ESs)
In the initial period, the effective strains (ESs) showed a value from 0.6% to 0.4% in the bottom of caverns with an enlarged bottom and at the roof of caverns with an enlarged top ( Figure 10A-D). In the five analysed withdrawal and injection periods, the distribution of ES in the cavern walls and surroundings was similar. In the withdrawal periods, the ESs were higher in comparison to the injection periods ( Figure 10). The largest value of ES (red area) was indicated in caverns with a diameter of 120 m and sharp edges ( Figure  10E,G,I,K). The area of these displacements was located in the bottom of caverns with an enlarged bottom and at the roof of caverns with an enlarged top. The highest value of ES

Effective Strains (ESs)
In the initial period, the effective strains (ESs) showed a value from 0.6% to 0.4% in the bottom of caverns with an enlarged bottom and at the roof of caverns with an enlarged top ( Figure 10A-D). In the five analysed withdrawal and injection periods, the distribution of ES in the cavern walls and surroundings was similar. In the withdrawal periods, the ESs were higher in comparison to the injection periods ( Figure 10). The largest value of ES (red area) was indicated in caverns with a diameter of 120 m and sharp edges ( Figure 10E,G,I,K). The area of these displacements was located in the bottom of caverns with an enlarged bottom and at the roof of caverns with an enlarged top. The highest value of ES in these caverns reached 7.0-10.0% (red and orange area locally) at the end of the first withdrawal period ( Figure 10E,G) and 9.5% at the end of the first gas injection period ( Figure 10M,O). A major increase was found at the end of the last withdrawal period (to 27.0%) ( Figure 10I,K) and at the end of the first gas injection period (to 26.0%) ( Figure 10O,X). However, these values of ES occurred locally (red area). The "smoothing" of the sharp edges to ss3 resulted in a decrease in ES ( Figure 10F,H,J,L,N,P,R,W). Moreover, in the last injection and withdrawal periods, the value and distribution of the ES were the same for the cylindrical caverns and the caverns with a maximum initial diameter of 75 m and where ss2 and ss3 were below 2.5% ( Figure 11).
The difference in the value and distribution of ES between the injection and withdrawal period was less visible for caverns with a diameter of 75 m (Figure 11). In the last cycle of both withdrawal and injection periods, the ES reached 5.0-7.0% (light blue area) at the bottom and the roof of caverns with a diameter of 75 m and se ( Figure 11D,G,M,P). In these groups, the distribution of the ES within the walls and surroundings of caverns with an enlarged bottom and top was the same. In the caverns with a maximum initial diameter of 120 m, the ESs were larger, and there were differences in their distribution between the enlarged-top and enlarged-bottom caverns ( Figure 11).
The cylindrical caverns showed better performance than the enlarged caverns ( Figure 11). The largest ES reached 2.5% in the last withdrawal and injection period ( Figure 11A-C,J-L). However, in caverns with a diameter of 75 m and with ss2 and ss3, the ES value and distribution were the same as in the cylindrical caverns ( Figure 11).
The ESs in caverns with a maximum initial diameter of 120 m were too large and involved a danger of plastic damage ( Figure 10). ES values (lower than 3.0%) that involved no potential danger of plastic damage were registered in caverns with a maximum initial diameter of 75 m ( Figure 11). The performed simulations indicated that the value and distribution of ESs were mainly affected by the cavern diameter. However, the "smoothing" of sharp edges influenced the value of ESs more than their distribution. in these caverns reached 7.0-10.0% (red and orange area locally) at the end of the first withdrawal period (Figure 10E,G) and 9.5% at the end of the first gas injection period ( Figure 10M,O). A major increase was found at the end of the last withdrawal period (to 27.0%) ( Figure 10I,K) and at the end of the first gas injection period (to 26.0%) ( Figure 10O,X)). However, these values of ES occurred locally (red area). The "smoothing" of the sharp edges to ss3 resulted in a decrease in ES ( Figure 10F,H,J,L,N,P,R,W). Moreover, in the last injection and withdrawal periods, the value and distribution of the ES were the same for the cylindrical caverns and the caverns with a maximum initial diameter of 75 m and where ss2 and ss3 were below 2.5% ( Figure 11).  The difference in the value and distribution of ES between the injection and withdrawal period was less visible for caverns with a diameter of 75 m ( Figure 11). In the last cycle of both withdrawal and injection periods, the ES reached 5.0-7.0% (light blue area) at the bottom and the roof of caverns with a diameter of 75 m and se ( Figure 11D,G,M,P). In these groups, the distribution of the ES within the walls and surroundings of caverns with an enlarged bottom and top was the same. In the caverns with a maximum initial diameter of 120 m, the ESs were larger, and there were differences in their distribution between the enlarged-top and enlarged-bottom caverns ( Figure 11).
The cylindrical caverns showed better performance than the enlarged caverns ( Figure 11). The largest ES reached 2.5% in the last withdrawal and injection period ( Figure 11A-C,J-L). However, in caverns with a diameter of 75 m and with ss2 and ss3, the ES value and distribution were the same as in the cylindrical caverns ( Figure 11).
The ESs in caverns with a maximum initial diameter of 120 m were too large and involved a danger of plastic damage ( Figure 10). ES values (lower than 3.0%) that involved no potential danger of plastic damage were registered in caverns with a maximum initial diameter of 75 m ( Figure 11). The performed simulations indicated that the value and distribution of ESs were mainly affected by the cavern diameter. However, the "smoothing" of sharp edges influenced the value of ESs more than their distribution.

Von Mises Stress
The value and distribution of vMS were equally related to changes in the cavern diameter and the "smoothing" of the sharp edges. The vMS reached the largest value (31.05 MPa-red area) in the initial period ( Figure 12A-D) and decreased after 4 years to 10.50 MPa in the first injection period. The largest value occurred in the bottom of the enlarged bottom caverns and the roof of the enlarged top caverns with a maximum initial diameter of 120 m, se and ss1. In the caverns with the same parameters but a maximum initial diameter of 90 m, the largest vMS occurred only locally ( Figure 12E,H)

Von Mises Stress
The value and distribution of vMS were equally related to changes in the cavern diameter and the "smoothing" of the sharp edges. The vMS reached the largest value (31.05 MPa-red area) in the initial period ( Figure 12A-D) and decreased after 4 years to 10.50 MPa in the first injection period. The largest value occurred in the bottom of the enlarged bottom caverns and the roof of the enlarged top caverns with a maximum initial diameter of 120 m, se and ss1. In the caverns with the same parameters but a maximum initial diameter of 90 m, the largest vMS occurred only locally ( Figure 12E,H).   (Figure 13A,E). However, the vMS value and area covered by its highe value were larger in the caverns with an enlarged bottom than in the enlarged-top cav erns, regardless of smoothing step (Figure 13). At the end of the withdrawal period, th distribution and values of vMS within the walls of caverns with ss3 and a maximum initi diameter of 75 and 90 m were similar (Figure 13C-D,G-H). At the end of the injectio period, the differences between these caverns were more clearly outlined ( Figure 13I-P) The value of the vMS was two times larger at the end of the withdrawal period than at the end of the injection period; e.g., in the caverns with a maximum initial diameter of 120 m and an enlarged bottom, it was 22.20 MPa and 11.90 MPa, respectively ( Figure 13A,I). vMS concentrations were found under the bottom of all analysed caverns with sharp edges (Figure 13A,E). However, the vMS value and area covered by its highest value were larger in the caverns with an enlarged bottom than in the enlarged-top caverns, regardless of smoothing step (Figure 13). At the end of the withdrawal period, the distribution and values of vMS within the walls of caverns with ss3 and a maximum initial diameter of 75 and 90 m were similar (Figure 13C-D,G-H). At the end of the injection period, the differences between these caverns were more clearly outlined ( Figure 13I-P).
The vMS showed a larger extension and value below the bottom of the enlargedbottom caverns. Generally, the "smoothing" process performed for all analysed caverns had a positive impact on the vMS (Figures 12 and 13). The results of the vMS simulations for the caverns with maximum initial diameters of 90 m and 75 m were positive. The value of the vMS was two times larger at the end of the withdrawal period than at the end of the injection period; e.g., in the caverns with a maximum initial diameter of 120 m and an enlarged bottom, it was 22.20 MPa and 11.90 MPa, respectively ( Figure 13A,I). vMS concentrations were found under the bottom of all analysed caverns with sharp edges (Figure 13A,E). However, the vMS value and area covered by its highest value were larger in the caverns with an enlarged bottom than in the enlarged-top caverns, regardless of smoothing step (Figure 13). At the end of the withdrawal period, the distribution and values of vMS within the walls of caverns with ss3 and a maximum initial diameter of 75 and 90 m were similar (Figure 13C-D,G-H). At the end of the injection period, the differences between these caverns were more clearly outlined ( Figure 13I-P).

Strength/Stress Ratio (SSR)
In the analysed time (9.5 years), there was a major difference in the value and distribution of the strength/stress ratio (SSR) between withdrawal and injection periods. The SSR in the initial period ranged from 4.0 to 6.0 (green and light blue area) at the walls of the caverns. It locally decreased to 3.0-3.5 (dark blue area) in the caverns with an initial maximum initial diameter of 120 m, precisely at the roof of the enlarged-top caverns (from se to ss2, Figure 14 A-C) and in the bottom of the enlarged-bottom caverns with se and ss1 ( Figure 14E,F).
In the first withdrawal period, a low SSR (1.80-2.5) was indicated in an area below the bottom of the enlarged bottom caverns with a maximum initial diameter of 120 m from se to ss2 ( Figure 15A-C). In this group, but with a maximum initial diameter of 90 m, this low SSR (1.80-2.5) occurred locally in the caverns with se and ss1; however, the caverns with a maximum initial diameter of 75 m were only characterised by the se (Figure 15I-L). The SSR increased to 3.0-4.0 below the bottom of the caverns with an enlarged top and a maximum initial diameter of 120 m ( Figure 15E-H). At the walls of these caverns, the SSR ranged from 4.0 to 5.0.
for the caverns with maximum initial diameters of 90 m and 75 m were positive.

Strength/Stress Ratio (SSR)
In the analysed time (9.5 years), there was a major difference in the value and distribution of the strength/stress ratio (SSR) between withdrawal and injection periods. The SSR in the initial period ranged from 4.0 to 6.0 (green and light blue area) at the walls of the caverns. It locally decreased to 3.0-3.5 (dark blue area) in the caverns with an initial maximum initial diameter of 120 m, precisely at the roof of the enlarged-top caverns (from se to ss2, Figure 14 A-C) and in the bottom of the enlarged-bottom caverns with se and ss1 ( Figure 14E,F). In the first withdrawal period, a low SSR (1.80-2.5) was indicated in an area below the bottom of the enlarged bottom caverns with a maximum initial diameter of 120 m from se to ss2 ( Figure 15A-C). In this group, but with a maximum initial diameter of 90 m, this low SSR (1.80-2.5) occurred locally in the caverns with se and ss1; however, the caverns with a maximum initial diameter of 75 m were only characterised by the se (Figure 15I-L). The SSR increased to 3.0-4.0 below the bottom of the caverns with an enlarged top and a maximum initial diameter of 120 m ( Figure 15E-H). At the walls of these caverns, the SSR ranged from 4.0 to 5.0. In the first injection period, the lowest SSR reached 6.0-7.0 in the area below the bottom of the caverns with an enlarged bottom and a maximum initial diameter of 120 m ( Figure 16). Results of numerical simulations for the last injection and withdrawal period In the first injection period, the lowest SSR reached 6.0-7.0 in the area below the bottom of the caverns with an enlarged bottom and a maximum initial diameter of 120 m ( Figure 16). Results of numerical simulations for the last injection and withdrawal period were similar to those of the first cycle. There were only some differences in the SSR value, which slightly decreased to 1.75-2.0 ( Figure 16).
The SSR seemed to be more influenced by caverns' dimensions and shape than by the "smoothing" of the sharp edges (Figures 14-17). This is illustrated by the performance of the cylindrical caverns compared to the caverns with a maximum initial diameter of 120 m (Figures 15 and 17) in the withdrawal periods. In the cylindrical caverns, the SSR values of 2.5-3.5 occurred locally below the bottom and in the middle part of the cavern with se and s1 ( Figure 18A,B). An SSR in a range of 1.75-2.5 covered the area below the bottom of all simulated caverns with an enlarged bottom and a diameter of 120 m, regardless of the "smoothing" of their edges ( Figure 15A,B and Figure 17C,F). The caverns with an enlarged top showed better performance than the enlarged-bottom caverns ( Figure 15E-H). The lowest SSR registered in these caverns was 3.0-4.0 in the withdrawal periods. In the first injection period, the lowest SSR reached 6.0-7.0 in the area below the bottom of the caverns with an enlarged bottom and a maximum initial diameter of 120 m ( Figure 16). Results of numerical simulations for the last injection and withdrawal period were similar to those of the first cycle. There were only some differences in the SSR value, which slightly decreased to 1.75-2.0 ( Figure 16). The SSR seemed to be more influenced by caverns' dimensions and shape than by the "smoothing" of the sharp edges (Figures 14-17). This is illustrated by the performance of the cylindrical caverns compared to the caverns with a maximum initial diameter of 120 m (Figures 15 and 17) in the withdrawal periods. In the cylindrical caverns, the SSR values of 2.5-3.5 occurred locally below the bottom and in the middle part of the cavern with se and s1 ( Figure 18A,B). An SSR in a range of 1.75-2.5 covered the area below the bottom of all simulated caverns with an enlarged bottom and a diameter of 120 m, regardless of the "smoothing" of their edges ( Figures 15A,B and 17C,F). The caverns with an enlarged top showed better performance than the enlarged-bottom caverns ( Figure 15E-H). The lowest SSR registered in these caverns was 3.0-4.0 in the withdrawal periods.

Safety Factor (SF)
The safety factor (SF) reflects the differences between the withdrawal and injection periods, as with other validated factors. In the initial period, the lowest safety factor was 4.12-5.0. The lowest SF of 2.5 (dark blue area) occurred locally in the withdrawal periods. There were slight variations in the value and distribution of the SF between the first and the last cycle. In the first injection period, an SF of 6.4-7.0 (dark blue area) was found at the roof of the enlarged and cylindrical caverns ( Figure 19). This range of SF occurred locally but was most intense at the roof of the enlarged-top caverns with ss1 ( Figure 19B-D). It was also found locally in the caverns with an enlarged bottom, ss1 and a maximum initial diameter of 75 and 90 m ( Figure 19E,F) but was not visible in the cavern with a maximum initial diameter of 120 m ( Figure 19G). The distribution of the SF changed in the last injection period ( Figure 20) and decreased slightly to 6.22. An SF in the range of 6.22 and 7.0 occurred in the roof of cylindrical caverns with ss1 and ss2, an enlarged top and a maximum initial diameter of 75 and 90 ( Figure 20A-H). It also occurred locally in caverns with an enlarged bottom and a diameter of 75 and 90 m ( Figure 20I-K).

Safety Factor (SF)
The safety factor (SF) reflects the differences between the withdrawal and injection periods, as with other validated factors. In the initial period, the lowest safety factor was 4.12-5.0. The lowest SF of 2.5 (dark blue area) occurred locally in the withdrawal periods. There were slight variations in the value and distribution of the SF between the first and the last cycle. In the first injection period, an SF of 6.4-7.0 (dark blue area) was found at the roof of the enlarged and cylindrical caverns ( Figure 19). This range of SF occurred locally but was most intense at the roof of the enlarged-top caverns with ss1 ( Figure 19B   In the first withdrawal period, the SF was lower and reached 2.5-3.0 within the middle part of caverns with ss3 and any diameter ( Figure 21A-E). The same SF value was found at the walls of caverns with a maximum initial diameter of 75 and 90 m and sharp edges ( Figure 21F-H). An SF of 2.48-2.50 occurred locally at the bottom of the enlarged-bottom caverns and at the roof of the enlarged-top caverns, both with a diameter of 120 m and with se ( Figure 21I,J). In the cylindrical caverns, this SF occurred locally only in the middle part of the cavern ( Figure 21A,E). Compared to the first withdrawal period, the changes in the SF value and distribution in the last withdrawal period were negligible ( Figure 21).
The performed evaluation indicated that the value of the SF was higher than 2.5 over the entire analysed time (Figures 19-22). These results show no danger of dilatant failure for all simulated caverns. However, the caverns with an enlarged bottom showed better results than caverns with an enlarged top. The SF was more affected by the shape and dimensions of the evaluated caverns but less influenced by the smoothing of the sharp edges. In the first withdrawal period, the SF was lower and reached 2.5-3.0 within the middle part of caverns with ss3 and any diameter ( Figure 21A-E). The same SF value was found at the walls of caverns with a maximum initial diameter of 75 and 90 m and sharp edges ( Figure 21F-H). An SF of 2.48-2.50 occurred locally at the bottom of the enlargedbottom caverns and at the roof of the enlarged-top caverns, both with a diameter of 120 m and with se ( Figure 21I,J). In the cylindrical caverns, this SF occurred locally only in the middle part of the cavern ( Figure 21A,E). Compared to the first withdrawal period, the changes in the SF value and distribution in the last withdrawal period were negligible ( Figure 21).  The performed evaluation indicated that the value of the SF was higher than 2.5 over the entire analysed time (Figures 19-22). These results show no danger of dilatant failure for all simulated caverns. However, the caverns with an enlarged bottom showed better results than caverns with an enlarged top. The SF was more affected by the shape and dimensions of the evaluated caverns but less influenced by the smoothing of the sharp edges.   The performed evaluation indicated that the value of the SF was higher than 2.5 over the entire analysed time (Figures 19-22). These results show no danger of dilatant failure for all simulated caverns. However, the caverns with an enlarged bottom showed better results than caverns with an enlarged top. The SF was more affected by the shape and dimensions of the evaluated caverns but less influenced by the smoothing of the sharp edges.

Volumetric Closure
At the end of the simulated operation time (9.5 years), the highest volume convergence (13.22%) was estimated for the cavern with se, a maximum initial diameter of 120 m and an enlarged bottom ( Figure 22C). The volume convergence was slightly lower for the same cavern with an enlarged top (12.33%, Figure 23B). The reduction in cavern diameter to 90 m resulted in a decrease to 8.02% (enlarged bottom) and 7.79% (enlarged top) ( Figure 23D,E). A smaller convergence was found for the enlarged-bottom and enlargedtop caverns with a maximum initial diameter of 75 m, with a value of 6.69% and 6.62%, respectively ( Figure 23F,H). The convergence of the cylindrical cavern with se was 5.81%. The "smoothing" of sharp edges had an impact on convergence, but this impact was stronger in the caverns with a maximum initial diameter of 120 and 90 m. The convergence in the cavern with an enlarged top, ss3 and a maximum initial diameter 120 m was 8.50%, but with a maximum initial diameter of 75 m, it decreased to 5.61%. Results for the enlarged bottom caverns were similar: the convergence of caverns with a diameter of 120 m was 8.83% and that of caverns with a diameter of 75 m was 5.65% ( Figure 23). The cylindrical shape was less affected by the "smoothing" of the sharp edges. Consequently, the convergence in cylindrical caverns ranged from 5.81% for se to 5.13% in ss3 (Figure 23).

Volumetric Closure
At the end of the simulated operation time (9.5 years), the highest volume conve gence (13.22%) was estimated for the cavern with se, a maximum initial diameter of 1 m and an enlarged bottom ( Figure 22C). The volume convergence was slightly lower f the same cavern with an enlarged top (12.33%, Figure 23B). The reduction in cavern d ameter to 90 m resulted in a decrease to 8.02% (enlarged bottom) and 7.79% (enlarged to ( Figure 23D,E). A smaller convergence was found for the enlarged-bottom and enlarge top caverns with a maximum initial diameter of 75 m, with a value of 6.69% and 6.62 respectively ( Figure 23F,H). The convergence of the cylindrical cavern with se was 5.81 The "smoothing" of sharp edges had an impact on convergence, but this impact w stronger in the caverns with a maximum initial diameter of 120 and 90 m. The convergen in the cavern with an enlarged top, ss3 and a maximum initial diameter 120 m was 8.50 but with a maximum initial diameter of 75 m, it decreased to 5.61%. Results for the e larged bottom caverns were similar: the convergence of caverns with a diameter of 120 was 8.83% and that of caverns with a diameter of 75 m was 5.65% ( Figure 23). The cyl drical shape was less affected by the "smoothing" of the sharp edges. Consequently, t convergence in cylindrical caverns ranged from 5.81% for se to 5.13% in ss3 ( Figure 23) The volumetric closure in simulated caverns is related to their dimensions and the operating pressure. The increase of gas pressure in the injection periods had a positive effect on the volume convergence because its value decreased. Otherwise, convergence increased in the withdrawal period with a decrease in gas pressure. The shape of caverns had a minor impact on the convergence, but the enlarged-top caverns showed a slightly lower convergence than the enlarged-bottom caverns. The "smoothing" of the sharp edges positively affected the volumetric closure, as the volume of cavern decreased in each smoothing step. An acceptable convergence (lower than 1% per year) was registered in the caverns with a maximum initial diameter of 90 m and 75 m.

Discussion
From the above analysis, it was found that caverns with an enlarged top outperformed caverns with an enlarged bottom when evaluated in terms of stability factors such as displacements, ES, vMS and SSR as well as volumetric closure. These results agree with the findings of Sobolik and Ehgartner [26] and Wang et al. [29]. The worse performance of caverns with an enlarged bottom was caused by the larger exposure of the cavern surface area to higher pressure. This pressure was caused by the difference between the in situ hydrostatic pressure and the operating gas pressure in the deeper part of the salt beds. However, the SF evaluation showed better results for caverns with an enlarged bottom because of the lower deviatoric stress in the cavern roof. As explained by Sobolik and Ehgartner [26], in caverns with an enlarged top, the angle of the slope of the cavern walls directs the pressure to hold up the cavern roof, thereby decreasing downward vertical displacement. This decrease in deformation is compensated by higher deviatoric stresses in the cavern roof, resulting in the possibility of dilatant damage and fracturing. Moreover, in the analysed caverns, the higher thickness of the roof pillar compared to the bottom pillar contributed to the better stability performance of the enlarged-top caverns.
Referring to the results of earlier works [2,25,32,35], cylindrical caverns showed the best performance of the analysed shapes.
The results of the numerical simulation indicated that the cavern dimensions have a great impact on the stability factors. The caverns with a maximum initial diameter of 120 m showed the worst geomechanical performance of all analysed caverns. On the contrary, most of the evaluated factors, except for displacements and SSR, showed similar results in cylindrical caverns with a maximum initial diameter of 75 m. Wang et al. [33] concluded that the dimensions of caverns have a greater effect on cavern deformation than the cavern shape. However, in the analysed cases, the performance of caverns with a maximum initial diameter of 75 m was almost the same, and the caverns with an enlarged top presented better results than caverns with an enlarged bottom if the maximum initial diameter was 90 m or 120 m.
The "smoothing" of the sharp edges performed in the three steps had a significant effect on the cavern performance in all analyses. Shapes with sharp edges had the worst results in the evaluation against stability factors. The assumption that sharp edges are prone to the stress concentration in the rock mass surrounding the caverns that caused the formation of damage zones [35,36] has been confirmed. The "smoothing" steps caused the roof and the bottom of the analysed caverns to be less flat and, as a consequence, better able to withstand deformations and applied load. Yang et al. [42] stated that caverns with a plain bottom should be avoided in design and construction. A flat roof of caverns was mentioned by Berest et al. [43], Coleman et al. [44], Reveille et al. [45], Wang et al. [2] and Osnes et al. [46] as a key factor in incidents of casing and roof stability loss. Moreover, a "smoothing" of sharp edges results in the formation of ellipsoidal caverns that are considered optimal for maintaining stability and resistance to volume convergence [3,29,47].
For the analysed shapes, the "smoothing" effect was stronger in the enlarged caverns than the cylindrical caverns. The geomechanical performance of the cylindrical shape was less affected by "smoothing" than the enlarged cavern. Moreover, the "smoothing" of the sharp edges contributed to the shrinkage in the volume of analysed caverns, resulting in lower susceptibility to volumetric closure.
In addition, all shapes and dimensions evaluated in terms of stability factors had a positive relationship with operating pressure. A high operating pressure (injection periods) improved the stability and performance of the analysed factors. This is because a high gas pressure can balance more load acting on the rock mass around the caverns; consequently, the creep deformation and volume shrinkage are reduced [1,43].
In this paper, the numerical simulations were performed with a k of 1 (the ratio of the horizontal stress to the vertical stress). However, some studies [48,49] that considered different k values showed an influence of stress regime on the relaxation zone and yielding zones around underground openings. An evaluation of different k values applied to numerical simulations of salt caverns and their shapes, especially those with an enlarged top and bottom, would be an interesting issue for future research.

Conclusions
In this paper, 28 caverns differing in shape and diameter were evaluated in terms of stability factors such as displacement, effective strain, von Mises stress, the safety factor and the strength/stress factor as well as volumetric closure. The stability performance of these caverns was compared, and the optimal cavern shape and volume in the geological conditions of the bedded salt deposit (the Mechelinki salt deposit) was found. Moreover, the performed evaluation aimed to determine the impact of sharp and smoothed edges on stability factors analysed in numerical simulations.
Numerical modelling indicated that the analysed stability factors are sensitive to changes in the initial cavern diameter and shape. Moreover, simulated shapes with sharp edges showed worse performance than shapes with smoothed edges. The caverns with a maximum initial diameter of 75 m (D/H ratio 5/8) with an enlarged top and bottom showed positive results. However, the caverns with sharp edges and ss1 showed a local increase in displacement in the middle part of the caverns and in the effective strain in the cavern's roof and bottom. The smoothing of the sharp edges in ss2 and ss3 contributed to the elimination of these areas. An overall comparison between the enlarged-top and enlarged-bottom caverns tended to favour enlarged-top caverns. This shape also showed promising results for caverns with a maximum initial diameter of 90 m, but only with ss3.
The cylindrical caverns showed the best results of all analysed shapes, but their volume was about 22% smaller than the enlarged caverns with a diameter of 75 m. Therefore, in the simulated geological conditions, the most optimal shape is an enlarged top with a maximum initial diameter of 75 m. The enlarged-bottom cavern with a maximum initial diameter of 75 m can also be considered on the grounds of the positive evaluation results. Moreover, the described results indicate that a possible local increase in the cavern diameter resulting from the leaching and petrological diversity of rock salt, causes no danger to cavern stability.
The stability performance of all analysed caverns showed a dependence on operating pressure. There was a large difference in the value of each stability factor between the injection and withdrawal period. The lowering of the operating pressure in the withdrawal periods had a negative impact on the evaluated factors. This impact increased in each operation cycle and was most visible in the displacement and effective strain simulation results. In future research, it would be interesting to consider the various values of k (the ratio between horizontal stress to vertical stress) and its impact on the cavern stability factors. The results presented in this paper can be used for the design process of new cavern fields for storing natural gas, compressed air or hydrogen. The described methodology can be adapted for the estimation of cavern storage capacity in bedded salt deposits. The results of the performed evaluation are also applicable to the design of optimal cavern shapes and dimensions. Moreover, the assessment of risks related to the cavern operation cycles and the planning of injection and withdrawal periods can benefit from the results described in this paper.

Supplementary Materials:
The following are available online at https://www.mdpi.com/2076-341 7/11/1/423/s1. SUPPLEMENT A, full results of the displacements analysis. SUPPLEMENT B, full results of the effective strain analysis. SUPPLEMENT C, full results of the von Mises stress analysis. SUPPLEMENT D, full results of the strength/stress ratio analysis. SUPPLEMENT E, full results of the safety factor analysis.