Geomechanical simulation of energy storage in salt formations

A promising option for storing large-scale quantities of green gases (e.g., hydrogen) is in subsurface rock salt caverns. The mechanical performance of salt caverns utilized for long-term subsurface energy storage plays a significant role in long-term stability and serviceability. However, rock salt undergoes non-linear creep deformation due to long-term loading caused by subsurface storage. Salt caverns have complex geometries and the geological domain surrounding salt caverns has a vast amount of material heterogeneity. To safely store gases in caverns, a thorough analysis of the geological domain becomes crucial. To date, few studies have attempted to analyze the influence of geometrical and material heterogeneity on the state of stress in salt caverns subjected to long-term loading. In this work, we present a rigorous and systematic modeling study to quantify the impact of heterogeneity on the deformation of salt caverns and quantify the state of stress around the caverns. A 2D finite element simulator was developed to consistently account for the non-linear creep deformation and also to model tertiary creep. The computational scheme was benchmarked with the already existing experimental study. The impact of cyclic loading on the cavern was studied considering maximum and minimum pressure that depends on lithostatic pressure. The influence of geometric heterogeneity such as irregularly-shaped caverns and material heterogeneity, which involves different elastic and creep properties of the different materials in the geological domain, is rigorously studied and quantified. Moreover, multi-cavern simulations are conducted to investigate the influence of a cavern on the adjacent caverns. An elaborate sensitivity analysis of parameters involved with creep and damage constitutive laws is performed to understand the influence of creep and damage on deformation and stress evolution around the salt cavern configurations. The simulator developed in this work is publicly available at https://gitlab.tudelft.nl/ADMIRE_Public/Salt_Cavern.

are similar to pure elastic deformations and can be significant in the secondary phase. The tertiary creep phase, seen when mean stresses are low enough to allow the formation of microscopic cracks, is characterized by microcrack damage evolution, hence creep acceleration, and ultimately brittle rupture of the material 21 . Creep in salt occurs at rates that are significant on short and long engineering time scales at temperatures in the range of 20-200 • C 22 and stresses as low as 0.2 MPa 22,23 . The deformation mechanics are governed by dislocation motion within the crystalline grains and by a range of grain boundary processes. As explained in the literature 24 , two main deformation mechanisms have been investigated by employing laboratory tests and micro-structural analyses. The first is the dislocation creep mechanism, which in the steady case is characterized by a power-law dependence of creep strain rate on deviatoric stress and an Arrhenius dependence on temperature. This type of behavior is grain size insensitive, and generally shows stress exponents in the range of 3.5-5.5 and activation energies around 60 kJ/mol 10,13,25 . When small (natural) quantities of brine are present, creep strains are over 5-10 %, and mean normal stresses are high enough to suppress micro-cracking, the process is accompanied by a growth of new crystals during deformation, which slightly enhances creep rate, but without changing the overall flow law significantly 24,[26][27][28] . Dislocation creep is favored towards higher stresses and temperatures through its highly nonlinear power-law stress dependence and Arrhenius temperature dependence. The second creep mechanism is solution-precipitation creep or pressure solution. It is a linear viscous creep mechanism that involves stress-driven dissolution-precipitation transfer of salt around water-bearing grain boundaries and is favored in fine-grained materials at low stresses, and temperatures 11,24 . Various formulations used in the literature to model creep and dilatancy are presented in Table 1. They are the power-law, Hou/Lux, Munson and Dawson (MD), and lastly, Hoek Brown models. For further details on these constitutive models, refer to the references in the Table 1. Note Table 1 does not include models that explicitly account for pressure solution. It should be mentioned that there are many other formulations, including the Hampel/Schulze 14 one used here, and several that do explicitly account for pressure solution 24,[29][30][31] .
The present work uses the Carter 10 constitutive model constants to describe the material behavior. The model is based on the power-law creep relationship derived by quantifying the creep processes observed in the laboratory tests on natural salt at differential stresses > 5 MPa. Laws for natural salt fitted above about 5 MPa will not contain pressure solution effects as it only becomes dominant in natural salt at lower stresses 10,13,14,26 . We do not include the pressure solution creep and recrystallization effects mentioned above, in the present analysis, due to the uncertainties in the ranges of stress and strain where these processes operate (see 27,28,32 ) and because we anticipate that hydrogen penetration into cavern walls will suppress these effects of water, due to desiccation. In salt bodies with a large lateral extent, pressure solution may become critical in controlling the far-field (low stress) behavior. Damage continuum mechanics is also incorporated in this work. Several cyclic loading experiments were conducted to study the damage characteristics, and fatigue of rock salt [33][34][35][36][37] . The formulations employed previously to study damage evolution of rock salt are presented in Table 2. In this work, the damage evolution process during the tertiary creep stage introduced by the Kachanov law is incorporated by introducing the dual damage variable D with the assumption that creep rate depends on the current damage state and the stress state of the system, according to the literature 38 . The damage law incorporated increases the effective stress which has been caused due to the dilatancy of rocksalt occurring from the damage of existing microcracks or the creation of new microcracks. It can cause a change in the permeability of rock salt from 10 −21 m 2 to 10 −17 m 2 , 39 . In this work, only geomechanical modeling of rock salt is performed without considering any change in permeability. Table 1. Constitutive relations presented in the literature to express creep strain rate and creep strain for Rock salt.
In the past years, several numerical studies [44][45][46][47][48][49][50][51][52][53] were conducted to describe geo-mechanical rock salt deformation. More recently, multi-scale finite element schemes (MSFEM) where also developed for reservoir rock mechanics [54][55][56][57][58] , for both linear and non-linear elastic deformations. Single-scale finite element models 59,60 were also developed to simulate salt caverns deformation under cyclic storage loading for homogeneous 2D caverns with simplified geometries using structured grids. However, only a few recent researchers have employed lab-scale constitutive models to reservoir scale analysis. One group of studies studied the creep behavior of Brazilian salt rocks using a multi-mechanism deformation creep model for symmetrical cavern shapes using commercial software ABAQUS 61,62 . They conduct axis-symmetrical simulations with different elastic properties for heterogeneities such as shale to confirm the cavern tightness and integrity for energy storage. However, caverns with complex geometries and different governing creep mechanisms for heterogeneities such as carnallite were not considered in this work. Another group of studies 63 conducted a stability analysis of symmetrical caverns subjected to cyclic loading for different salt mines using Perzyna's visco-plastic model in code bright finite element solver. Each salt mine was considered with different elastic and visco-plastic properties. However, individual mines have a different lithological composition comprising of different minerals and they can undergo plastic deformation with distinct governing equations due to varying crystal structure. To address the above challenges from the literature and distinctly account the heterogeneity in the lithology around salt caverns, we present an open-source python-based computational FE framework that addresses the influence of complex cavern shapes, varying elastic and creep properties, and suitable creep governing equations in reservoir scale on the state of stress and deformations around the caverns.
The heterogeneous composition of the geological domain around salt caverns will affect the deformation of salt caverns. Elastic and creep properties of these rock formations vary along with the depth of the geological domain. Interlayers in salt formations are difficult to dissolve in water and pose many challenges in the process of designing caverns and during their operation [64][65][66][67] . Caverns within bedded salt formations are not so stable as those created within salt domes 3 . Heterogeneity of the rock salt may affect the solution mining process, which in turn will also affect the shape of the constructed cavern 6 . In field scale cases, salt properties are heterogeneous due to impurities such as anhydrite or potassium/magnesium salts, shale rocks, bischofite 7,19,68 . Heterogeneity in the geological domain surrounding salt caverns involves many insoluble interlayers (anhydrites, potash salts, shale, gypsum, mudstone, etc.) 35,65,69 . Few experimental researchers have attempted to study the impact of heterogeneity in the composition of rock salt 19,65,68,70 . However, very few researchers have tried to study the impact of heterogeneity by considering elastic and creep properties of interlayers in the domain 19 . Wang et al. 64 showed the effect of heterogeneity by considering only elastic properties and by considering equivalent elastic modulus. Because of different crystal lattices and compositions, every material behaves differently under stress subjected to longer timescales. Hence, it becomes critical to incorporate creep properties of the heterogeneity. In this work, we study two approaches to study the impact of heterogeneity. The first method employs only elastic properties of the heterogeneity. The second method uses different elastic properties and creep constitutive relations to model creep aptly subjected for extended periods. One more aspect that becomes critical is studying the effect of the irregularly shaped cavern on the state of stress and deformation. This work attempts to study this effect allowing us to capture real field scenarios in the computational domain.
The computational framework employs an implicit simulation model on non-uniform fully-unstructured triangular mesh and constant strain triangular (CST) elements. The model is further expanded by incorporating Eulerian strains and considering the stored product (hydrogen) density and corresponding hydrostatic pressure. The rest of the paper is structured as follows. Firstly, creep and damage constitutive laws governing equations are presented. Numerical methodology is presented to account for creep consistently. Next, a comparison of numerical and experimental results is shown. Followed by, the influence of creep under monotonic and cyclic loading on salt caverns is discussed. Then the impact of complex, realistic cavern geometries and material heterogeneity inside the geological domain with different material properties and governing equations are elaborated. Then, the evolution of rock salt damage is studied, and a detailed sensitivity analysis of all the parameters in the chosen constitutive laws is presented. Lastly, the impact of multi-cavern simulations on stress is studied.

Governing equations
In this work, conservation of momentum is employed to solve for drained solid with constant pore pressure, where the σ and f b are the stress and body force terms. Note that the inertial terms are not included as for their minimal impacts. Through the theory of the linear elasticity, only the elastic part of the total strain is proportional to the stress field, i.e., here, C is the 4th rank elasticity or stiffness tensor 56 . For 2D isotropic homogeneous discrete elements, one can express it through the Lame's constants and µ as where δ ij is the Dirac delta.
When time dependent plastic flow (creep) deformation is considered, the elastic deformation is only part of the total strain 8 . In this case, the total strain can be expressed as where ∇ s = 1 2 (∇u + ∇u T ) 54 . Note that u T is the transpose of the displacement vector. In addition, in this work, derivatives of displacements are taken with respect to the deformed discrete geometries. Finally, the governing equation for displacement is found as For Eq. (7) to be well-posed, a closure term for the creep strain is necessary. It is elaborated in the next section.
Creep formulation. In this work, creep strain rate is expressed as Norton-Bailey power law using the parameters used for Halites by Carter law 10 as presented in Table 1. Using Carter law and under multiaxial stress condition the creep strain is expressed as 41,71 , The above formulation incorporates temperature dependency as expressed by the Arrhenius law. here s is deviatoric part of the stress tensor, and Q, R and T denote the activation energy, Boltzmann's constant and temperature, respectively 8,13,27 . In this work it is also assumed that this flow rule is still valid for small volume changes that can accommodate damage law.
Tertiary creep. In this work, damage is only considered in the tertiary stage of creep. To study tertiary phase of the creep which might involve initiation and propagation of microcracks leading to rupture it is important to study the damage mechanics of the rock salt. Damage mechanics of rock salt is still an active field of research both on the experimental and numerical studies.
In this work, the Kachanov law 19,38,40 is employed to incorporate damage in studying rock salt. The creep damage rate is a function of stress σ and the current damage state D. The constitutive creep equation is therefore stated as where the damage state variable is expressed through the evolution equation as Here, D * is the critical value of the damage, at which the given material breaks, and Ḋ is the damage evolution rate, expressed as Finally, the creep strain rate is give ny 19,40 Variables a, B, n, r in the Eqs. (11) and (12) above represent the material dependent constants. It is also important to note that in the case of D = 0 , Eq. (12) represents the well known power-law constitutive equation. More generally, Eq. (12) can be written as, where σ d is the dilatancy intensified stress given by σ d = σ /(1 − D).
The material parameters to incorporate damage strain used in 59 cannot be used here due to different time scales. The time period used to study failure criteria of bedded rock salt in 59 is less than 10 days. However, since the focus of our work is mainly for long term energy storage (months/years) the time scale is off by a magnitude of 10-100. The value of the parameters B is the main parameter that is considered to vary in our work. Due to lack of available experimental data in large timescales, the value of the parameters B is assumed to be 100 times the magnitude in the paper 19 mainly due to the difference in timescales.

Numerical methodology
Equation (7) together with Eq. (12) forms a well-posed system for nonlinear time-dependent deformation vector u = (u, v) of salt rock with elastic and inelastic deformation. Here, u and v stand for displacement in x and y directions, respectively.
The discrete system in space is found by using finite-element method (FEM), which can be stated as (4) ε = ε el + ε cr , ε cr =ε cr (σ , D),  54 . Hereafter, since the solution is obtained on a single mesh resolution, the super-index h is avoided for simplicity of the descriptions. In this work, 2D triangular mesh is used to employ finite element methodology. The deformation is stored in the nodes of the triangle and the stress is computed in the centre. Figure 3a shows the mesh element with deformation stored at the corner nodes. The shape functions in the local coordinates system ( ξ , η ) are written as, The shape functions are derived from the ratio of area of the triangular element opposite to the node i to total area. The areas can be seen in Fig. 3a. Using these shape functions in Eq. (14) allows us to compute deformation and the same shape functions can be used to compute global coordinates of an interior point. The strain is computed using horizontal and vertical deformation (u, v), To compute differentials of displacements, shape functions and chain rule is employed considering local coordinates as shown 8 .
Stress is further computed from Eq. (2). When heterogeneity is involved, each element will have different Lame parameters because of different elastic properties. Using the principle of virtual displacements 8 and Gauss divergence theorem, the final discrete equation for an element with unit thickness is obtained from Eq. (7) is, The above area integrals are approximated using Gauss quadrature rule by transforming it from global to local coordinate system. Here, D e = div T N e , the applied stress on the boundary along the normal is t (Neumann bc) which is integrated over area Ŵ , the volumetric force is f , integrated over an elemental volume e . Here, fictitious forces F cr = e D T Cε cr dA are computed from creep strain. When the time dependent creep strain is included, two main time-integration approaches can be pursued: implicit (Euler backward) and explicit (Euler forward). In this work, both approaches are developed to allow for stability and efficiency analyses in a more comprehensive approach. In explicit formulation, the creep strain in continuum form is is obtained by integrating its rate functions, i.e., In discrete form, it can be simply stated as The above formula states that creep at new time step (n + 1) can be obtained by the already-available information at time step n. The explicit approach is summarised in Algorithm 1. Here, the stress term in the creep model is obtained based on the new time step (n + 1) . In this case, the Newton-Raphson iterative approach is needed to solve for the displacement iteratively. More precisely, the nonlinear residual reads The objective is to find u n+1 , which also allows one to find ε n+1 and σ n+1 . The linear form of this equation is then obtained using which is then solved iteratively until convergence is achieved, i.e. R n+1 = 0 . Here, the Jacobian J = ∂R ∂u i is computed based on the elastic part of the residual equation. The creep term is lagged by one iteration, i.e., and as such does not contribute in the Jacobian matrix. This allows for considering the creep term in the righthand-side of the iterative linearised system, i.e., The iterations i are repeated until the residual norm falls below the prescribed threshold, i.e., ||R i || 2 < ǫ r . This is called the convergence state. Algorithm 2 provides an overview of the implicit time-integration scheme for non-linear displacement modelling.

Results
This chapter presents numerical results of a series of 2D test cases beginning with the benchmarking of the developed simulator, and then to quantify the impact of nonlinear time-dependent creep physics. For this study, 8 test cases are considered: (1) linear elastic deformation under constant load, (2) consistency verification of linear elastic model, (3) nonlinear (creep) deformation under constant load, (4) nonlinear (creep) deformation under cyclic load, (5) complex heterogeneous model with creep deformation, (6) irregular cavern shape model with and without heterogeneity, (7) real field geometrical data, and (8) tertiary creep and material failure. The input parameters, used for the first 6 test cases, are summarized in Table 3. The constants chosen here are obtained from lab-scale experimental data. These constants are used as a base test case for the remaining test cases. However, these constants can vary when energy is stored in the reservoir scale. Then accordingly, these constants have to be scaled considering the sensitivity of the parameter, loading conditions, and local heterogeneities beyond the scope of this work.
Defining the load for the boundary conditions. To reflect real-field loading conditions, the cavern's fluid and lithostatic pressure values are set such that the fluid pressure does not exceed 24% − 80% of lithostatic pressure 4 . The minimum pressure difference between the cavern's fluid pressure and lithostatic pressure will be at the cavern's roof during the injection; and, the maximum will be at the cavern's bottom during production (see points 1 and 2 in the left Fig. 4). The cavern's fluid pressure, thus, is a function of the overburden rock density, rock salt density and depth of the roof and the bottom of the cavern. Roller boundary condition was imposed at the bottom and the face along the cavern. Traction free Neumann boundary conditions were imposed on the top and the far end face of the geological domain. This allows to observe any subsidence on the top face or any deformation that can cause in the far of the geological domain.
The minimum and maximum cavern pressure and corresponding pressure difference are used to calculate the equivalent surface forces acting on the cavern's wall, as shown in Fig. 4 (right). The density of hydrogen is employed to compute the forces on the cavern walls. These surface forces are then converted into equivalent nodal forces, which are finally used in the numerical model as input parameters.
Triangular mesh with refinement around the caverns was used for the simulations. A mesh convergence study was also conducted. Consistency and 2nd order accuracy of deformation in the implementation for the elastic domains are confirmed.  Table 3. Classical creep methodology is employed in this paper.
The comparison between numerical and experimental results is satisfactory. The numerical results can be further optimised if necessary using optimization algorithms. However, in this work the chosen parameters for the creep constitutive law is obtained from the past literature 10 . Accordingly using these constants, the numerical methodology is validated and further compared with the experiments. All the further numerical experiments are conducted using the parameters from this test case as tabulated in Table 3.  Table 3. A constant fluid pressure of 20% of lithostatic pressure with respect to time for 275 days is imposed on the caverns. The time-dependent solutions are illustrated in Fig. 6, which shows the displacement evolution over time and von Mises stress distribution across the domain. A higher magnitude of displacement is observed around the cavern near larger curvatures at the end of 275 days. This is due to the reason of high stress accumulation near the curvature causing higher creep deformation. Due to the applied loading, the cavern volume shrinks. Figure 6 also demonstrates the evolution of strain ε xx , ε yy and stress σ vM over time, under constant load. Here the evolution of strain is linear because of the chosen formulation for the strain rate, i.e., Eq. (8). Higher rates of strains ( ε xx , ε yy ) are observed for points A and B near the cavern roof and floor due to a high stress distribution caused by low surface area. Points B and F show slightly lower strain rates. Lastly, the points close to the mid-plane of the domain (C, D and E) have very low strain rates compared to the rest of the labeled points. This shows that if there are any cracks or heterogeneity near the curvature, there is a higher probability that it will lead to failure than near the mid plane region. The loss of volume of the cavern given by, The change in the boundary of the cavern is not significant. In this case, the approximate loss of volume is around 3 % after 275 days for a constant load.
Due to the low permeability of salt caverns, they are considered of best use in storage technology. However, due to volumetric strain permeability of salt rocks can increase, causing the stored gases to penetrate. Authors 63,75 propose the permeability of salt caverns which is related to dilatant volumetric strain given by Here, constants α k and β k are model parameters subjected to change due to the material properties and loading conditions. In this work α k = 2.13e − 8 and β k = 3 is employed as suggested by 75 . Khaledi et al. 63 showed the higher permeability near the roof and floor of the caverns. A similar observation is made here, as shown in Fig. 6f. These zones also depict the potential fail zones around the cavern that can lead to gas leakages. However, it is important to see how the contours of permeability are observed in the geological domain with complex shapes as in the real field, which is not shown in the previous studies 63 . This is further elaborated in the future sections.
Test Case 3: Creep under cyclic loading. When excess renewable energy is produced, it is converted to green gas and stored in the subsurface. Depending on the supply and demand of the energy, the gas will be injected and produced. This will result in the cyclic loading of salt caverns. This test case addresses the important aspect of energy storage, i.e., the deformation under cyclic loading. For this reason, the cavern's fluid pressure is assumed to be a function of time. Here, a discrete step function is used as shown in Fig. 7. The maximum and minimum pressure applied during cyclic loading depends on the lithostatic pressure. It varies between P max = 80% and P min = 20% of lithostatic pressure. Figure 8 shows the variation of horizontal and vertical deformation with  www.nature.com/scientificreports/ time for different points located around the cavern. Horizontal displacement is 0 for points A and G due to the imposed roller boundary conditions as shown in Fig. 4. The full cycle of the cavern's loading and unloading in this simulation is evaluated for 6 days. The high peak values are related to the instantaneous elastic response of the rock salt material, after which there is a short period of creep development with a linear slope representing the magnitude of the creep rate. The higher the load, the steeper the slope and the higher the creep rate. It can be seen that the creep deformation is insignificant at P min compared to P max (zero slope). From Fig. 8a it can be seen that the rate of deformation for every cycle is the highest point, F. Followed by the points E, D, C, and B. This trend is observed because of incorporating the lithostatic pressure in the geological domain. So accordingly, highest depth is observed at point F, causing the highest deformation. Similarly, from Fig. 8b it can be seen that the highest magnitude of vertical deformation is observed at point G, and the next highest magnitude is observed at point A. Although these points are symmetrically placed around the cavern along the x axis mid-plane, unequal vertical deformations are observed due to the lithostatic pressure. Accordingly, the magnitude of deformation reduces from highest stress to lowest stress (in the center point D) with the least curvature. In cyclic load, the loss of volume after 275 days is less than 1.5 %.
Test Case 4: Caverns with irregular geometries. In this test case, irregular-shaped caverns are studied for both homogeneous and heterogeneous domains. Figure 9 shows the simulation results of the model with such an irregular cavern shape. From the figures, it can be seen that the maximum displacement is higher than that of the regular cylindrical-shaped cavern. Note that a larger surface is subject to the force exerted by the pore-fluid pressure fluctuations for irregular shapes. Also, it is noted that-for the shown geometry-the displacement is lower in the magnitude at the top and bottom of the cavern when compared with the cylinder-shaped cavern. Point E shows the maximum horizontal deformation in Fig. 9d at the end of 275 days and also has the highest rate of deformation. This is again due to large stresses at point E with high curvature. Figure 9e shows the variation of vertical deformation with time. Since there is no boundary condition imposed on the vertical deformation highest magnitudes of the rate of vertical deformation is seen in Point A and Point G. The potential failure zones in  www.nature.com/scientificreports/ the domain is point around E and the points near the floor and roof of the cavern. Due to the irregular shape of the caverns, the potential local failure zones are more important because of the high-stress zones, which might lead to micro-cracking or lead to damage. The above point is supported by high permeability values near high curvature, as shown in Fig. 9f. Compared to permeability observed in the cylindrical cavern, as shown in the previous test case, higher permeability is observed near point E than near the roof or floor of the cavern. Therefore, high permeability zones are potential failure zones that can reduce the cavern's tightness and lead to stored hydrogen or any other gas leakage. When heterogeneity is incorporated, different local properties or constitutive equations can cause different stress distributions in the geological domain leading to different potential failure zones. This is investigated in the following sections.

Test Case 5: Heterogeneity of elastic material properties.
Heterogeneity is incorporated in this work by introducing a potash interlayer within a homogeneous halite rock salt deposit, as shown in Fig. 10 Fig. 11. The magnitudes of deformation in x and y directions are comparable. Moreover, noticeable is that the introduction of potash interlayer has impacted the stress field and the deformation. For example, the point E within the interlayer (i.e., pink line in Fig. 11) shows higher deformation compared with the homogeneous test case. However, it shows slower creep rates for heterogeneous case, due to the lower local stress at the location of E. The maximum deformation change due to introducing the heterogeneous potash layer for this test case after 275 days is found to be approximately 9.5 % and 38 % for horizontal and vertical deformation, respectively. This is found by comparing the values of Figs. 9 and 11. Figure 11f shows the variation of permeability in the domain, computed from the volumetric strain as stated in Eq. (27). It is observed that the permeability near the cavern wall within the heterogeneous potash is www.nature.com/scientificreports/ slightly smaller, compared with the homogeneous case. Note that here, only elastic parameters are considered to be heterogeneous. These results motivate the next test case to investigate the impact of full heterogeneity in both elastic and plastic properties for the potash interlayer.
Test Case 6: Heterogeneous interlayers. The second, more realistic method is employing different insoluble layers and using their elastic and creep properties. To incorporate heterogeneity in the geological domain surrounding salt caverns, interlayers considering many insoluble impurities (anhydrites, potash salts, shale, gypsum, mudstone, etc.) 35,65,69 are included. The detailed composition of rock salt and interlayers can vary depending on the type of region. Experimental analysis needs to be conducted to understand the structure and the lithological composition of the geological domain 65,79 . In this work, interlayers of different materials are shown in Table 4 are chosen as heterogeneity in the geological domain of halite. The table shows the constitutive parameters of creep formulations for different materials are taken from the literature. Considering i) different laboratory conditions, ii) timescales for conducting experiments, ii) water content in the rock samples, iv) different loading conditions, and v) purity of the samples chosen in the literature, these creep constitutive parameters can change depending on the type of application it is being used for instance, in this case, being energy storage. This study becomes very critical given right parameters of the interlayers are chosen. The experiments were conducted for roughly about 500 hrs, and at 70 • C for both carnallite and bischofite 68 . This work used these constants due to the lack of other suitable literature for energy storage applications. The creep constants are reduced by three orders of magnitude for both carnallite and bischofite because of the difference in temperature the experiments are conducted and the current case study because of their temperature dependence from the Arrhenius formulation as shown in 80 . Figure 12 shows the interlayers in the geological domain. The width of the interlayer is 30m. The interlayers are placed in the two regions near the salt cavern. Figure 12a shows the interlayer around the midplane of the cavern. Locally around that interlayer, the curvature is minimum. Figure 12b shows the interlayer near the roof of the cavern with the highest curvature and more lithostatic pressure compared to the midplane. The entire simulation was conducted for 50 days. Figures 13 and 14 shows various graphs for carnallite and bischofite respectively. Von mises stress distribution, horizontal deformation at initial and final time step, and variation of horizontal and vertical deformation with time for points around the cavern as shown in the schematic Fig. 4 are shown for both the interlayer distributions. The magnitudes of von Mises stress for both carnallite and bischofite materials are the same. However compared to the homogeneous test case of cylindrical cavern Fig. 6a, the minimum value of stress has reduced and the distribution is not smooth anymore. This is due to the discontinuous distribution of Young's modulus and the magnitudes of Young's modulus of interlayers is lower than halite. Due to this, the magnitudes of horizontal deformation at t=0 are similar as shown in Figs. 13b,h, 14b,h. However, the magnitudes of these deformations are higher than the homogeneous test case as shown in Fig. 6b. The deformation magnitude is higher for interlayers in the mid than the interlayer near the floor due to low stress in the central region and low curvature.
When creep is incorporated in these interlayers with separate properties, it's an additional nonlinear physics that depends on the stress distribution. Figures 13c,i and 14c,i show the horizontal deformation distribution at end of simulation after 50 days when interlayer are incorporated comprising carnallite and bischofite respectively. Due to different stress distributions caused by lithostatic pressure and curvature, the accumulated creep deformation also varies. Figures 13c and 14c show a similar distribution of deformation due to the same location of heterogeneity; however, the bischofite shows higher deformation magnitudes compared to carnallite. This can be     Table 4. www.nature.com/scientificreports/ because even though the creep exponent for carnallite is higher than bischofite, the creep constant of bischofite is higher than carnallite by five orders of magnitude. The sensitivity of the creep constant and creep exponent is studied in the later section. Figures 13i and 14i show the horizontal deformation when the interlayer is located near the floor of the cavern. Qualitatively, they look similar; however, the interlayer with bischofite shows higher horizontal deformation than the carnallite test case and a homogeneous test case. Figure 13d,e,j,k show the variation of horizontal and vertical deformation with time at different points as shown in the schematic for carnallite interlayers. Figure 14d,e,j,k show the variation of horizontal and vertical deformation with time at different points as shown in the schematic for bischofite interlayers. Compared to the homogeneous test case as shown in Fig. 6, the horizontal deformation plots for each point are more widely spread when the interlayer is in the midplane. However, when the interlayer is near the floor, Point F shows a much steeper slope than the rest of the points. This shows that given a local heterogeneity in the domain, a large amount of deformation could occur, causing failure of the salt cavern. In the vertical deformation, the plots are magnified proportionately compared to homogeneous test cases where the magnitude is higher but the distribution looks similar when the interlayer is the midplane. When the interlayer is near the cavern floor, the magnitude of vertical deformation at point F with high curvature and heterogeneity has a higher rate of increase in deformation than a homogeneous test case. Permeability of salt caverns which is obtained from volumetric strain, can be seen in the graphs Fig. 13f,l for carnallite with interlayers in the mid strip and near the floor of the cavern. Figure 14f,l show the variation of permeability with interlayer as bischofite in the mid and bottom strip, respectively. When the interlayer is the midplane, we can see that the potential failure zones are the floor, roof (same as without interlayers), and the midplane section with the interlayer. This potential failure zone is created only because of the interlayer heterogeneity. When the interlayer heterogeneity is near the cavern floor, the main potential failure zone is near the floor around point E. Followed by the roof of the cavern location. Here the impact of material and geometrical heterogeneity has combined, resulting in a very high potential failure zone. Also, the magnitude of bottom strip interlayer permeability is higher. Similar permeability distribution is observed for both carnallite and bischofite when the interlayers are in the mid or near floor, respectively.
The maximum increase in percentage for horizontal and vertical deformation in the geological domain when compared to homogeneous test case (Fig. 6) after 50 days of simulation is shown in Table 5. The percentages are higher for bischofite material compared to carnallite. When the interlayers are present near the cavern floor, a higher increase in deformation is observed than the mid-plane interlayer. Also, the percentage increased for both the deformations when the creep properties are incorporated to study heterogeneity is much higher than when only elastic properties are included, as shown in the previous test case. Considering the chosen creep parameters are in the realistic range, the results from this section clearly show that interlayers can fail the cavern locally around the cavern. Test Case 7: Real field test case. Using the field data on cavern shape from echo logs of a salt cavern in Germany 2 , a cross-section of the modeled cavern was generated to be used in the developed simulator. The test case simulation results are shown in Fig. 15. The computed displacement distribution, the simulation output, has a similar distribution to the cylinder-shaped cavern. Due to the increased lithostatic pressure, the maximum horizontal deformation is just below the midplane of the cavern. While the maximum horizontal deformation appears at the cavern roof, lower deformation values are predicted due to the Dirichlet boundary conditions along the bottom boundary. This simulation was run for two years.
Test Case 8: Tertiary creep and material failure. In this test case, the evolution of damage within the material is studied for a given parameter set. Figure 16 show the transition from secondary creep to tertiary state and subsequent material failure. These graphs are plotted for a node inside the domain with a maximum strain, where the material failure would begin. In a homogeneous geological domain with rock salt, the location of that node is observed to be on the cavern's wall. This is because the cavern wall would undergo maximum stress in given time duration. From the results, it can be seen that as the creep forces increase over time, the strain rate becomes steeper, causing the magnitude of the deformation to increase drastically. Depending on the chosen damage parameters, the rate at which the material fails can be calibrated.
The damage constants used in the reference 19 cannot be used in this formulation mainly due to different time scales and different rock salt properties. Figure 17a shows for different B (constant parameter in the formulation) the damage parameter increases significantly. No literature is available to compare the numerical results with damage parameters incorporated with the timeline of 250 days. Hence the damage parameters had to be assumed by considering the same physics involved as explained in the "Tertiary creep".
Damage continuum analysis using Kachanov law involves three material constants B, r, σ . The variation of damage parameter D with time for different B, r, σ is shown in Fig. 17. Figure 17a shows for different B (constant parameter in the formulation) the damage parameter increases significantly. Lower the value of B higher the D. In 19 the parameter B is lower by 100 in magnitude due to the time scale less than ten days. So from this analysis, it can be said that the parameter B depends on the timescale of the simulation significantly. Figure 17b shows the variation of damage parameter for different r. Lower the value of r higher is the damage parameter. Higher the stress higher is the expected damage can be seen in Fig. 17c. Multi-cavern simulations. Multi-cavern simulations are helpful when more than one cavern can be used close to each other in the geological domain for subsurface energy storage. Cavern to cavern distance (CTC) between the caverns is critical to understand the influence of each cavern on the other. In this study, multicavern simulation studies were conducted for regular and irregular-shaped caverns. CTC is the minimum cavern to cavern distance possible in the geological domain. CTC for regular cylindrical-shaped cavern is shown in   Fig. 18. Figure 19 shows the same for irregular-shaped caverns. It can be seen that as the distance between the caverns reduces, causing an increase in Von-mises stress. These simulations can be used when there is heterogeneity in the geological domain which can handle less stress. Accordingly, Drucker Prager or Von Mises failure criteria could be used to determine how critical the distance is between the caverns 81 .
Vertical deformation ( u y ) is compared here since traction-free boundary conditions are applied on the top face. It can be seen that the closer the caverns it higher the vertical deformation is observed. When the caverns have complex geometry, the magnitude of the Von-Mises stress s vm and vertical deformation is higher. This could lead to subsidence 29 in the geological domain or could amplify seismicity. Figure 20 shows the variation of von Mises stress along the horizontal distance in the mid-plane section for regularly shaped caverns and the minimum distance (CTC plane) for complex shaped caverns. CTC mp test case shows the variation of von Mises stress at the mid-plane section of the irregular-shaped cavern. Here from Fig. 20b it can be seen that the von Mises stress distribution is more uniform similar to the Gaussian distribution. As CTC reduces, the minimum stress in the center of the homogeneous geological domain increases.  Figure 15. Test Case 7: The above contours show the horizontal displacement (a), vertical displacement (b) and von Mises stress (c) for real field cavern shape model. The simulation was run for a period of 2 years.  www.nature.com/scientificreports/ Regularly-shaped caverns have minimum surface area due to their smooth surfaces. Hence, there is no higher stress observed near the cavern walls. The stress reduces slowly from peak stress closer to the cavern and drops to a minimum value near the central vertical plane of the geological domain. From Fig. 20c, it can be seen that the stress distribution at the center is more spread out, and closer to the cavern, it peaks out. This is due to the irregular shape. From the curve representing CTC mp it can be seen that the stress increases at first and then reduces again. This is because this section is not the closest distance with the adjacent cavern and not the surface with the most curvature.

Sensitivity analysis.
Creep is a complex phenomenon that involves various parameters that can influence the deformation and stresses over time. To quantify the impact of the critical parameters, sensitivity analysis needs to be conducted. The parameters chosen are cavern depth, creep constant, Young's modulus of halite rock, temperature, and creep exponent. Energy storage technology could involve different operating conditions and heterogeneous properties of rock salt. Due to this, the above parameters are chosen to study their influence on the time-dependent deformation. Figure 21a shows the two points (A and B) in the domain where sensitivity analysis is conducted. The influence of Young's modulus on horizontal deformation is shown in Fig. Fig. 21b. Change in Young's modulus changes the elastic deformation. Higher the Young's modulus, lower is the deformation. The salt cavern's far-field shows a lower slope than the node closer to the cavern because of higher stresses closer to the cavern. The sensitivity of deformation towards shear modulus is not shown here since it can be expected that the deformation is not sensitive towards it. Figure 21c shows the variation of horizontal deformation with different creep multipliers. Higher the creep multiplier, the higher the rate of change in deformation with time, leading to a higher accumulation of plastic deformation. The increase in creep constant is non-linearly increasing the deformation after a certain period.
The influence of temperature on deformation is seen in Fig. 21d. The higher the temperature, the higher is the accumulated creep deformation. Again, a similar trend is observed as the last plot where the temperature change is not linearly changing the change in deformation. The impact of creep exponent n is shown in Fig. 21e. The higher the creep exponent, the higher is the deformation. Compared to creep constant and temperature,  www.nature.com/scientificreports/ creep exponent is very sensitive to the accumulated deformation. For n = 3 and n = 3.25 the variation of deformation with time is almost horizontal, however for creep exponent, n = 3.5 , the rate of change in deformation is very high. Comparing the sensitivity plots for the creep parameters, it can be seen that the influence of creep exponent is biggest, followed by creep constant, and lastly, it is creep temperature. Figure 21f shows the variation of horizontal deformation with time for different depths of the cavern. In the geological domain with a lot of heterogeneity and interlayers, it is important to construct the cavern in the domain with minimum cracks and heterogeneity. Depth of the cavern is a parameter that can allow this argument. The higher the depth, the higher is the lithostatic pressure leading to higher creep deformation.

Conclusion
In this work, the influence of complex shapes and material heterogeneity in the geological domain on salt caverns employed for energy storage technology is studied using a 2D finite element solver. The secondary (i.e., steadystate) creep behavior is introduced in the mathematical model based on power law, and the power-law parameters were taken from the work of Carter et al. 10 . The numerical methodology and the chosen creep constitutive constants are compared with the experimental data. Tertiary creep is introduced by utilizing the damage evolution parameter, allowing to predicting the material failure 8 . The developed simulator allowed for both explicit and implicit time integration schemes. After the consistency check and benchmarking with the experimental data, Various conclusions can be derived from this work.
(1) It is evident from the above results that creep is a slow phenomenon of insignificant magnitude compared to the elastic component for year-long operations. This justifies the present strategy of omitting the low-  www.nature.com/scientificreports/ stress creep mechanism of pressure solution from the current analysis. However, on a scale of several years with a significant amount of heterogeneity, the effect of creep strains on the deformation and stresses can become evident, which would eventually cause the material to reach the tertiary stage leading to failure at zones of critical stress intensity. Less than five % volume changes in the cavern were observed considering cylindrical homogeneous caverns after 275 days. (2) The role of heterogeneity by considering only different elastic parameters is showed for irregular-shaped caverns. Depending on the type of impurity (potash, halite), the location of impurity in the domain, and the type of distribution imposed in the geological model, stress and deformation distributions vary. (3) When heterogeneous interlayers are incorporated in the simulation model with the right constitutive models for the composition, a much higher increase in deformation is observed locally and around the cavern than the homogeneous test case and the heterogeneous test case without considering creep properties (test case 5). In general, bischofite showed higher creep deformations compared to carnallite. When the interlayers are located near the high curvature region, that region has the highest potential to be a failure zone. (4) The influence of curvature of the caverns was studied in this work for both with and without heterogeneity.
Interlayers near the high curvature showed at least five times more deformation than interlayers located near the minimum curvature. Potential failure zones by computing permeability were also identified to be in the regions around high curvature. (5) Tertiary creep was also incorporated in this model with the assumption creep rate depends on the damage state and stress of the system. The underlying physics could be included however the constants involved in the damage law had to be assumed due to unavailable literature for more extended time scales. (6) A detailed sensitivity analysis of various parameters shows their influence on the deformation, which can gauge the deformation in real field test cases. (7) The developed open-source simulator was extensively tested on various test cases, as shown in the paper.
The developed model is publicly available as an open-source simulator in the TU Delft repository of the ADMIRE project located at https://gitlab.tudelft.nl/ADMIRE_Public/Salt_Cavern. The developed code allows the user to use different features simultaneously to understand the underlying field-relevant physics quickly. (8) Multi-cavern simulations were critical to understanding the significance of CTC, irregular shape of the caverns on the stress distributions in the geological domain. The role of additional creep mechanisms such as pressure solution and uncertainties in dislocation creep parameters will be considered in the future, also the effects of gas (hydrogen or green methane) permea- www.nature.com/scientificreports/ tion on damage development in the cavern wall. Future work involves including thermal and (visco)plastic strains to allow for more reliable simulations and sensitivity analyses. Since the field test cases are large-scale of the order (km), efficient formulations such as multi-scale formulation would be developed to reduce the computational costs. Further research would be required to benchmark the tertiary creep model with lab and field data. Other heterogeneous interlayers such as shale rocks need to be studied considering the right creep constitutive laws and underlying physics are chosen.