Effect of pillar width on the stability of the salt cavern field for energy storage

The assessment of optimal pillar width was based on numerical modelling. The numerical simulations aimed to find the optimal pillar width in given geological conditions. The pillars between solution mined caverns are irregular in shape, and this makes geomechanical analysis difficult (Zhang et al., 2021). To avoid the influence of the cavern shape on the results of analysis, a shape of cavern was based on cylinder.

Previous study (Cyran & Kowalski, 2021) demonstrated that sharp edges in caverns lead to stress concentration and reduced stability. To address this issue, the cavern edges were smoothed by incorporating a quadrant with a radius of 10 m in the corners of a vertical cross-section of the cylindrical cavern shape (Fig. 1). Consequently, the original cylindrical shape of the cavern was transformed into a cylinder with slightly smoothed edges. This adjustment aimed to mitigate the stresses associated with sharp edges, which have an impact on pillar stability. The dimensions of examined caverns were tailored to the geological conditions of the Mechelinki salt deposit. The analysed shape is based on the assumed cavern height (H) of 120.0 m and a maximum diameter (D) of 80.0 m (Fig. 1). The thickness of the roof pillar was set to 15 m, while the bottom pillar had a thickness of 5 m.

Figure 1:

To precisely asses the pillar width in a cavern filed, all the caverns in the filed should be modelled. In the present paper, the layout of the cavern field was based on concentric rings to ensure an equal distance (pillar width) between adjacent caverns and optimize the utilization of the salt deposit. The cavern field consisted of three rings (circles) with varying numbers of caverns. In the numerical analysis, three variations of cavern layouts associated with circles, as shown in Fig. 2, were considered. The first variant of the cavern layout comprised seven caverns (one central and six surrounding caverns). In the second variant, the cavern layout encompassed circles 1 and 2, totalling 19 caverns, including 7 from circle 1 and 12 from circle 2. The third variant of the cavern layout incorporated three circles, resulting in a total of 37 caverns, including 7 from circle 1, 12 from circle 2, and 18 from circle 3. Additionally, three variants of pillar width were evaluated: 1D (80.0 m), 2D (160.0 m), and 3D (240.0 m), respectively. A 60° wedge section represents caverns from each circle, including the cavern in the centre.

Figure 2:

The geomechanical model of the cavern field was built based on the geological structure of the Mechelinki salt deposit. The salt beds belong to the Zechstein salt formation (cyclothem PZ1) and are situated 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. The salt beds are underlain and overlain by anhydrite layers: the lower anhydrite (A1d) and the upper anhydrite (A1g). The thicknesses of these anhydrite layers range from about 1.5 to 20.0 m. The salt beds are underlain by Kupferschiefer (T1), Zechstein Limestone (Ca1), and Rotliegend and Silurian sediments. The PZ2 and PZ3 sediments overlay the salt beds, along with Triassic, Jurassic, Cretaceous, and Cenozoic strata (Czapowski et al., 2007).

The numerical simulations were conducted for cavern fields consisting of 7, 19, and 37 individual caverns (Fig. 2). Boundary conditions were applied in the form fixed displacement in the normal direction on the side and bottom planes. The geometry of the model took the shape of a wedge (Fig. 2) with the base being an equilateral triangle, with sides measuring 2500 m and a height of 1500 m. The direct roof and the bottom of the rock salt beds consisted of anhydrite layers with a thickness of 30 m at the roof and 20 m at the bottom. Above and below the anhydrite layers, rock mass was located (Fig. 3). The upper layer with a thickness of 60 m was made up of soils. The roof of the rock salt beds was located at 975 m b.g.l. (Fig. 3). The mesh consisted of approximately 3.8 million elements, predominantly hexahedral. To enhance the efficiency of numerical simulations and ensure an accurate representation of the pillar width, the element size at the sidewalls of the pillars was smaller than that of the others (approximately 4 m). Outside the pillars, the sidewall elements were larger (around 50 m). The initial value of the hydrostatic stress varied with depth, ranging from 0 MPa at the top to 36.00 MPa at the bottom of the 3D model.

Figure 3:

The rock mass mechanical behaviour was simulated with the use of the mechanical parameters and the constitutive models of the rock salt and the nonsalt rocks. Two different constitutive models were chosen to accurately project the viscoelastic plastic behaviour of rock salt and the elastic–plastic response of nonsalt rocks. The two-component Norton Power Law with a Mohr–Coulomb plasticity criterion was applied to describe the plastic yielding and creep of rock salt. However, the Mohr–Coulomb elastic–plastic model was used to simulate the mechanical behaviour of the surrounding rocks and the anhydrite.

The mechanical parameters of rock salt and nonsalt rocks were determined based on laboratory tests and validated in previous studies (Cała et al., 2018; Cyran et al., 2021). The parameters presented in Table 1 describe four materials: rock salt, anhydrite, rock mass, and soil. These mechanical parameters were applied to the numerical calculations. In addition, creep parameters for rock salt were determined in laboratory tests and calculated based on the Norton Power law. The creep parameters for rock salt used in the numerical simulations are n = 5.0 and A = 1.08 · 10⁻⁴⁵ Pa^−5.0 s⁻¹.

The evaluation of the pillar width was based on numerical modelling. The 3D geomechanical model of the cavern field was introduced into FLAC3D software, which is a tool for solving geomechanical and geotechnical problems, including rheological phenomena.

The relationship between pillar width and stability was assessed with the use of the following criteria (stability factors): deformations including the total displacements, vertical and horizontal displacements, von Mises stress (vMS), strength/stress ratio (SSR), and safety factor (SF). The total displacements were analysed at the sidewalls of the pillars because they reflected the changes in the pillar shape and volume caused by the salt creep. The vertical displacements show an overburden impact on the pillar. The horizontal displacements show the influence of the volume changes in salt caverns on the pillar, causing the pillar to shrink or expand. The vMS is the equivalent effective stress that is related to the creep rate within the primary and secondary creep stages. The von Mises stress (vMS) is expressed by equation (1): (1) σvm=[(σ1−σ2)2+(σ2−σ3)2+(σ3−σ1)2]2 {\sigma _{vm}} = \sqrt {{{[{{({\sigma _1} - {\sigma _2})}^2} + {{({\sigma _2} - {\sigma _3})}^2} + {{({\sigma _3} - {\sigma _1})}^2}]} \over 2}} where σ_vm is a von Mises stress, σ₁ is a major principal effective stress, σ₂ is an intermediate principal effective stress, and σ₃ is a minor principal effective stress. The vMs indicates zones of large stress concentration which are important due to difference in mechanical behaviour between rock salt and other rocks such as anhydrite overlying and underlying rock salt. Rock salt exhibits rheological behaviour, but in the case of anhydrite, rheology can be ignored.

The strength/stress ratio (SSR) is used to indicate the dangerous areas in the pillars and surroundings. The strength–stress ratio (SSR) is a local indicator of a current stress state and its proximity to failure. SSR is described by equation (2): (2) SSR=|σ′1−σ3σ1−σ3|≤10 SSR = \left| {{{{{\sigma '}_1} - {\sigma _3}} \over {{\sigma _1} - {\sigma _3}}}} \right| \le 10 where σ′₁ is the minimum effective principle stress in the failure state, σ₁ is the current minimum effective principle stress, and σ₃ is the current maximum effective principle stress.

An SSR equal to 1.0 implies a material failure, and furthermore, an SSR of 2 indicates that the material reaches 50% of its strength.

The safety factor (SF) is based on a dilatancy damage criterion for rock salt that can be expressed by equation (3): (3) SF=b⋅I1I2 SF = {{b \cdot {I_1}} \over {\sqrt {{I_2}} }} where b is a material constant, I₁ is the first invariant of the stress tensor, and I₂ is the second invariant of the deviatoric stress tensor. When SF < 1, the shear stresses in the salt (J₂) are large relative to the mean stress (I₁), and dilatant behaviour is predicted. SF indicates material damage when it is below 1.0 (DeVries et al., 2005; DeVries, 2006; Sobolik et al., 2006; Sobolik & Ehgartner, 2006; Van Sambeek et al., 1993). This damage threshold indicates the onset of damage. It is important to note that a very short-term occurrence of a SF less than 1.0 does not necessarily result in immediate salt fracturing. The greater concern would be a value less than 1.0 over a period of several weeks or months, indicating that the accumulation of damage would cause fracture generation (Sobolik, 2016).

In addition, the changes in the pillar volume reflected by the cavern volume shrinkage (convergence) were evaluated. The cavern volume shrinkage was calculated as the ratio of the volume loss to the initial cavern volume.

The numerical analysis was carried out for a period of 12 years. Water leaching was employed for excavating the salt caverns. Upon completion of excavation, brine filled the entire salt cavern for a duration of 2 years subjecting the internal surface of the salt cavern to hydraulic pressure. Subsequently, brine was discharged from the caverns by natural gas, initiating cavern operation. During operation, load conditions were simulated by changes in cavern pressure. Ten operation cycles were simulated with cavern pressure oscillating between a minimum of 4 MPa and a maximum of 17.5 MPa. For each operation cycle, the assessment criteria mentioned above were analysed for 4 operation periods: at the beginning of the gas injection period, at the end of the gas injection period, at the beginning of the gas withdrawal period, and at the end of the gas withdrawal period. The numerical simulations assumed simultaneous and instantaneous leaching of all caverns in the applied software. Consequently, the displacements during the leaching period were minimal compared to those during the operation period, and their role in long-term cavern stability could be neglected. To mitigate the adverse effect of this assumption on the results, in-situ stress redistribution was calculated initially and used as the initial conditions for subsequent calculations, including creep.

Displacements in the pillars between adjacent caverns are primarily caused by the creep behaviour of rock salt. The analysis focused on the total, vertical, and horizontal displacements of the rock mass within the pillars and surrounding the caverns. As a result of rock salt creep, the volume of the caverns decreases over time, leading to changes in the pillar width. These changes are visible in the displacements maps.

The distribution of total, vertical, and horizontal displacements during the last withdrawal period is presented in Figs. 4, 5, and 6, respectively. In the maps depicting vertical displacements, upward displacements are represented as positive values, while downward displacements are shown as negative values (Fig. 5 A–I). The vertical displacements demonstrate the influence of the overburden load on the pillars. Regarding horizontal displacements, both positive and negative values are observed (Fig. 6A–I). The sign of the displacement indicates its direction along the X axis. Positive horizontal displacements are directed outward from the centre of the cavern field, while negative horizontal displacements are directed inward toward the centre of the cavern field.

The maximum displacements, indicated in red, are predominantly located in the contour zone surrounding the caverns. These values gradually increase as the width of the pillar between two adjacent caverns decreases and the number of caverns in the ring increases (Fig. 4 A–I).

The results indicate that when only caverns in variant 1 of the cavern layout are analysed and the pillar width is 1D, 2D, or 3D, the maximum displacements are 114.99 cm, 109.34 cm, and 106.57 cm, respectively (Fig. 4A, D, G). However, when caverns in variant 3 of cavern layout are taken into account, the maximum displacements in the pillar are 155.05 cm, 157.37 cm, and 139.50 cm, respectively (Fig. 4C, F, I). Increasing the distance between caverns leads to changes in the distribution and magnitude of displacements in the pillar between adjacent caverns and the area above the caverns. Notably, the contour zone around the caverns in variant 3 of cavern layout exhibits significant differences among the pillar widths of 1D, 2D, and 3D (Fig. 4C, F, I).

The maximum displacements of 139.5 cm and 157.37 cm are observed around caverns 1 and 2 (inner caverns) with a pillar width of 3D and 2D (Fig. 4F, I), respectively. However, for a pillar width of 1D, the maximum displacements are only found around cavern 4 (outer cavern) (Fig. 4C). It should be noted that increasing the pillar width to 2D results in an increase in the maximum displacement values for caverns in ring nos 2 and 3. However, for a pillar width of 1D, large displacements ranging from 90.0 to 100.0 cm are observed in a larger area (shown in green) (Fig. 4C). This area includes the upper part of the pillars between caverns 1, 2, and 3, as well as a significant zone above these caverns. Increasing the pillar width to 2D has a positive impact, as the displacements in the upper part of the pillars and above the caverns decrease to 60.0–70.0 cm between caverns 1 and 2 and 50.0–60.0 cm between caverns 2 and 3. Further increasing the pillar width to 3D results in additional decreases in displacements to 30.0–40.0 cm (shown in blue) in the upper part of the pillars and above the caverns.

Increasing the number of caverns in the ring affects the total displacements and their distribution. The large displacements, reaching 80.0–90.0 cm (Fig. 4B) between adjacent caverns 1 and 2, are observed in the upper part of the pillars with a width of 1D in ring no. 2. These displacements increase to 100.0–110.0 cm in ring no. 3 (Fig. 4C). Similar trends are observed for a pillar width of 2D. The displacements in the upper part of the pillars between adjacent caverns 1 and 2 range from 40.0 to 50.0 cm in ring no. 2 (Fig. 4E) and increase to 60.0–70.0 cm in ring no. 3 (Fig. 4F).

Vertical displacements in the pillar decrease as the pillar width increases. However, an increase in the number of caverns in the ring leads to higher values and extensions of vertical displacements. The largest downward displacements are mainly located above the cavern roof (blue) and below the bottom (orange), transferring to the pillars (Fig. 5A–I). The vertical displacements in the upper parts of the pillars and caverns are predominantly downward, while those in the lower parts are mainly upward. On average, the vertical displacements of the pillars are small, but locally they can be larger. For a pillar width of 1D, vertical displacements range from −70 cm to −90 cm in the upper part of the pillars and from 30 cm to 50 cm in the lower part (Fig. 5B, C). The load is mainly applied to the inner caverns (1 and 2) in ring nos 2 and 3. Increasing the pillar width to 2D decreases the magnitude and extension of vertical displacements (Fig. 5D–F), but the area of downward displacements (−50 cm to −60 cm) in the upper part of the pillar between the inner caverns (1 and 2) remains. Enlarging the pillar width to 3D further reduces this area.

The horizontal displacements in the caverns and pillars contribute to the shrinkage of cavern volume and the tension in the pillars (Fig. 6A–I). The distribution of horizontal displacements shows that they are the largest in the contour zones around the caverns and gradually decrease towards the centre of the pillar. However, this zone is wide when the pillar width is 1D and extends over most of the pillar (Fig. 6A–I). For a pillar width of 1D, there is a significant difference between the values of positive and negative horizontal displacements. This difference decreases as the pillar width increases and reaches almost the same value for a pillar width of 3D. In the cross-section through ring no. 1, the maximum values of horizontal displacements are −114.92 cm and 63.35 cm for a pillar width of 1D (Fig. 6A), whereas they are −106.38 cm and 103.90 cm for a pillar width of 3D (Fig. 6C). This large difference in horizontal displacements for a pillar width of 1D causes unevenly distributed tension in the pillars. The largest horizontal displacements, directed inside the centre of the cavern field, are distributed in the contour zone of all caverns when the pillar width is 2D and 3D (Fig. 6D–I), but they are distributed in the contour zone of the outer caverns when the pillar width is 1D (Fig. 6A–C). An increase in the number of caverns in the field leads to an increase in both positive and negative horizontal displacement values for all analysed pillar widths (Fig. 6A–I). Interestingly, there is an increase in the value of horizontal displacements for a pillar width of 2D in both rings 2 and 3 (Fig. 6E, F) compared to a pillar width of 1D and 3D.

The results indicate that the total and horizontal displacements in the pillars are largest in the contour zone around the caverns and gradually decrease towards the centre of the pillar (Figs. 4, 6). The vertical displacements are largest in the upper and lower parts of the pillars (Fig. 5). The value and range of total and vertical displacements decrease with an increase in the pillar width (Figs. 4–5), but the value of horizontal displacements increases with the pillar width (Fig. 6). In general, the values of total, vertical, and horizontal displacements increase when the number of caverns increases (Figs. 4–6). The largest vertical displacements are associated with the pillars between adjacent caverns located in the inner part of the cavern field (Fig. 5). The largest total displacements are found in the inner part of the cavern field consisting of 19 or 37 caverns when the pillar width is 2D and 3D (Fig. 4E, F, H, I). However, when the cavern field consists of only 7 caverns, the largest total displacements occur at the edge of the cavern field (Fig. 4D, G). Similarly, when the pillar width is 1D, the largest total and horizontal displacements are located at the edge of the cavern field (Fig. 4 B, C). The increase in the value of total and horizontal displacements for ring nos 2 and 3, and a pillar width of 2D, indicates that this pillar width is insufficient to provide stability. Therefore, the allowable pillar width for a cavern field consisting of 19 and 37 caverns should be between 2D and 3D.

The maps in Figure 7A–I present the values and distribution of von Mises stress (vMS) during the last withdrawal period. The maximum vMS value (red colour) is associated with the area below the cavern bottom and ranges from 21.63 to 21.13 MPa. The extent of these areas remains the same regardless of the pillar width and number of caverns. The range of vMS values between 14.0 and 16.0 MPa is located in the area above the roof of the cavern. The extent of these areas is more influenced by the number of caverns rather than the pillar width.

The overall vMS values in the pillars are relatively low, mostly ranging between 8.0 and 10.0 MPa (Fig. 7 A–I). However, in certain areas of the pillars, the vMS ranges from 10.0 to 12.0 MPa when considering ring nos 2 and 3 (Fig. 7 B, C, E, F, H, I). Additionally, as the pillar width increases, the vMS decreases. For example, in variant 3 of cavern layout, the vMS reaches 8.0–10.0 MPa and 6.0–8.0 MPa when the pillar width is 2D and 3D, respectively.

The results in Figure 8 show that areas of the pillars with slightly higher vMS values (10.0–12.0 MPa) are associated with the contours surrounding the caverns. The distribution of these areas varies for different pillar widths and the number of caverns in the ring. When the pillar width is 1D, the vMS of 10.0–12.0 MPa is associated with the contours surrounding the outer caverns (Fig. 7 A–C). The same distribution of vMS values is observed for in variant 1 of cavern layout and pillar widths 2D and 3D (Fig. 7 D, G). However, for pillar widths 2D and 3D and ring nos 2 and 3 (Fig. 7 E, F, H, I), this value (10.0–12.0 MPa) of vMS is found in the contours surrounding all caverns. Moreover, the results show an increasing proportion of areas where the vMS value is 6.0–8.0 MPa with an increase in the pillar width.

It is worth noting that there is a difference in vMS between the beginning and end of the last withdrawal period. Although the maximum vMS value is slightly higher (21.33–21.63 MPa) at the end of the withdrawal period compared to the beginning (21.13–21.33 MPa), the distribution of vMS values differs between these two periods (Figs. 7, 8). The vMS ranges from 10.0 to 12.0 MPa in the contour zone around the caverns (Fig. 8A–D). This zone is larger at the beginning of the withdrawal period (Fig. 8B, D) and occurs around all caverns in cross-sections through ring nos 2 and 3. Furthermore, the vMS values of 8.0–10.0 MPa in the contour zone around all caverns are larger and occupy a larger area of the pillars at the beginning of the withdrawal period compared to the end (Fig. 8A–D).

The above results indicate that the vMS values in the pillars are generally low compared to the area below the cavern bottom (Figs. 7, 8). The largest value of vMS (10.0–12.0 MPa) in the pillars is related to the contour zone around the caverns and gradually decreases to the centre of the pillar where their value is the lowest (6.0–8.0 MPa). The overall vMS magnitude decrease when pillar width increases and increase with the number of caverns (Fig. 7A–I). Based on the above results, the pillar width 2D is safe for the long-term stability of the cavern field consisting of 7, 19, and 37 caverns (Figs. 7, 8).

During the analysed time period of 12.5 years, there were variations in the strength/stress ratio (SSR) in the pillars between adjacent caverns. The lowest SSR values (1.73–1.87, dark blue colour) were observed in the area below the caverns (Fig. 10A–I). The extent of these areas is primarily influenced by the pillar width and decreases as the pillar width increases. Additionally, the SSR ranges from 3 to 3.5 over the roof of the cavern when the pillar width is 1D (Fig. 10 A–C) and 3.5–4.0 when the pillar width is 2D and 3D (Fig. 9 D–I). These SSR values and their distribution do not significantly affect the overall stability of the caverns.

In general, the SSR in the pillars between caverns is high, with most areas having values above 8.0 (Fig. 9A–I). In the contour zone around the caverns, the SSR value decreases to 4.5–5.5 when the pillar width is 2D and 3D (Fig. 9 D–I). When the pillar width is 1D, the SSR value of 7.0–8.0 covers most of the pillar area (Fig. 9 A–C).

The results presented in this section indicate that the SSR in the pillar is generally high compared to the area below the caverns (Fig. 9A–I). The SSR value is primarily influenced by the pillar width, where a larger pillar width corresponds to a higher SSR value. The lowest SSR values (4.5–5.5) in the pillars are associated with the contour zone around the caverns, and the SSR gradually increases towards the centre of the pillar (Fig. 9A–I). Based on these results, the SSR in the pillar, considering all analysed pillar widths and the number of caverns, is sufficiently high to ensure cavern stability. However, for overall cavern stability, including the area above the roof and below the bottom, the minimum allowable pillar width is 2D for all analysed cavern fields.

The SF (safety factor) value and its distribution exhibit variations for the analysed pillar width and number of caverns in the cross-section. The lowest SF values range from 3.09 to 3.32 and are locally indicated in the contour zone around the caverns (Fig. 10A–I). These SF values satisfy safety requirements for cavern stability, indicating no risk of dilatancy damage in the pillars during the analysed time period of 12.5 years.

Overall, the SF in the pillars is high, with most values exceeding 5.0 and some falling below 4.0 (Fig. 10A–I). In the contour zone around the caverns, the SF reaches 3.5–4.5 and gradually increases towards the centre of the pillar, where its value is the highest (Fig. 10A–I). These SF values are observed in areas above the roof and below the bottom of the caverns for all analysed variants. Some changes in the SF value and distribution are mainly influenced by the pillar width. For a pillar width of 1D, most areas of the pillar exhibit SF values above 5.0–6.0 (Fig. 10G–I). However, when the pillar width is 2D and 3D, the SF of 5.0–6.0 occurs only in the contour zone around the caverns (Fig. 10D–I). The distribution of the lowest SF is affected by both the pillar width and the number of caverns in the cross-section. When the pillar width is 1D, the lowest SF values (3.09–3.32) occur in individual elements (Fig. 10A–C). The lowest SF values of 3.17–3.28 are locally indicated around the middle part of the caverns for a pillar width of 2D and 3D in ring nos 2 and 3 (Fig. 11D–I).

Based on the above results, it is evident that the overall SF in the pillar is sufficiently high to meet safety requirements. The value and distribution of SF are primarily influenced by the pillar width (Fig. 10A–I). The SF in the pillars increases with the pillar width. The lowest SF values (3.0–3.5) in the pillars are associated with the contour zone around the caverns and gradually increase towards the centre of the pillar (Fig. 10A–I).

The results presented in Figure 11 demonstrate the changes in the cavern volume caused by injection and withdrawal cycles, resulting in volumetric shrinkage or convergence. The convergence is significantly higher during the withdrawal cycle compared to that during the injection cycle. These changes in the cavern volume also affect the pillar volume. Figure 11 illustrates that the impact of pillar width on convergence increases with the number of caverns. The influence of pillar width on convergence is relatively low for ring no. 1 (Fig. 11A) due to lower convergence values ranging from 4.5% (1D) to 5.3% (2D and 3D). For a cavern field consisting of 19 caverns (Fig. 11B, ring no. 2), the convergence ranges from 4.9% (1D) to 6.0% (2D). The largest influence of pillar width on volumetric shrinkage is observed for pillar ring no. 3 (Fig. 11C, 37 caverns), with convergence values of 5.6%, 6.2%, and 6.8% for pillar widths of 1D, 2D, and 3D, respectively. The highest volume shrinkage (6.80%) is estimated for a pillar width of 2D and a cavern field consisting of 37 caverns (ring no. 3). Similarly, for ring no. 2 (19 caverns), the pillar width 2D (Fig. 11B) showed the largest convergence (6.0%). The convergence is the lowest for pillar width 1D (Fig. 11A–C) and was estimated between 4.50% and 5.70%. These results align with the values and distribution of total and horizontal displacements (Figs. 4A–I, 6A–I).