# Effects of Stress Extrapolation and Dilatancy on Salt Cavern Design: An Application to Underground Hydrogen Storage

Jubier A Jiménez-Camargo, Laura Blanco-Martín, Ahmed Rouabhi, Faouzi Hadj-Hassen

► To cite this version:

Jubier A Jiménez-Camargo, Laura Blanco-Martín, Ahmed Rouabhi, Faouzi Hadj-Hassen. Effects of Stress Extrapolation and Dilatancy on Salt Cavern Design: An Application to Underground Hydrogen Storage. Rock Mechanics and Rock Engineering, 2025, 10.1007/s00603-025-04897-8. hal-05252182

HAL Id: hal-05252182

https://hal.science/hal-05252182v1

Submitted on 12 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effects of Stress Extrapolation and Dilatancy on Salt Cavern Design:
An Application to Underground Hydrogen Storage

Jubier A. Jiménez-Camargo<fnref n="1" /> · Laura Blanco-Martín<fnref n="1" /> · Ahmed Rouabhi<fnref n="1" /> · Faouzi Hadj-Hassen<fnref n="1" />

Received: 20 December 2024 / Accepted: 25 August 2025

© The Author(s) 2025

## Abstract

The long-term behavior of rock salt is commonly predicted using creep laws fitted on laboratory data at medium-high deviatoric stresses (≥ 5 MPa), which do not fully capture the relevant stress range for underground applications. Indeed, it is a common practice in numerical modeling to extrapolate high-stress experimental data to estimate material behavior at lower stress levels. Additionally, the analysis of the simulation results is often combined with specific design criteria, with particular focus on dilatancy. While some dilatancy criteria are embedded at the constitutive level, others are used only during post-processing, or are of empirical nature. This study simulates hydrogen storage in salt caverns using two in-house finite element codes and the RTL constitutive model. This model accounts for deviatoric and volumetric viscoplastic strains, incorporates dilatancy at the constitutive level and has been calibrated with recent multi-stage confined creep tests over a differential stress range of 0.5–20 MPa and temperatures of 20–50 °C. The results of the thermomechanical simulations show that extrapolating high-stress data to low-stress regimes can lead to significant underestimations of cavern volume loss and to misrepresentation of stress and strain magnitudes and profiles. Likewise, incorporating dilatancy at the constitutive level offers more accurate predictions of the onset and extent of dilatancy zones. Such zones stabilize over time, remaining smaller than the region influenced by very low deviatoric stresses. As emphasized in this study, accurate assessments of stability require that laboratory experiments and the calibration of constitutive models be conducted as much as possible within the ranges of stress, temperature, strain and strain rates relevant for the underground application.

## Highlights

* Extrapolation from high to low differential stresses can lead to underestimation of cavern volume loss and misrepresentation of stresses and strains.

* Creep laws should rely highly on experiments covering ranges of stress, temperature, strain and strain rates relevant to underground applications.

* Incorporating dilatancy criteria at the constitutive level and across a wide stress range improves predictions of dilatancy onset and extent.

* From a design perspective, considering a wider range of deviatoric stresses enables more efficient engineering decisions on cavern spacing.

**Keywords** Rock salt creep · Dilatancy · Hydrogen storage · Stress extrapolation · Thermomechanical modeling

Laura Blanco-Martín, Ahmed Rouabhi, Faouzi Hadj-Hassen have
contributed equally to this work.

<figure><img src="image_2.png" /></figure>

Jubier A. Jiménez-Camargo
jubier-alonso.jimenez camargo@minesparis.psl.eu

Laura Blanco-Martín
laura.blanco_martin@minesparis.psl.eu

Ahmed Rouabhi
ahmed.rouabhi@minesparis.psl.eu

Faouzi Hadj-Hassen
faouzi.hadj_hassen@minesparis.psl.eu

<footnote n="1">¹ Centre de Géosciences, Mines Paris, Université PSL, 35 rue Saint-Honoré, Fontainebleau 77300, France</footnote>
## List of symbols

<table><tr><td>σ</td><td>Cauchy stress tensor</td></tr><tr><td>σ'</td><td>Deviatoric component of σ</td></tr><tr><td>ε̇<sub>vp</sub></td><td>Viscoplastic strain rate tensor</td></tr><tr><td>I, J, N</td><td>Unit tensors forming an orthonormal basis</td></tr><tr><td>q = √(3/2) ||σ'||</td><td>von Mises equivalent stress</td></tr><tr><td>p = -1/3 tr(σ̄)</td><td>Mean pressure</td></tr><tr><td>θ = 1/3 arccos (√6tr(J̄<sup>3</sup>))</td><td>Lode angle</td></tr><tr><td>γ<sub>vp</sub></td><td>Equivalent deviatoric visco-plastic strain</td></tr><tr><td>ζ<sub>vp</sub></td><td>Equivalent volumetric visco-plastic strain</td></tr><tr><td>ρ</td><td>Mass density</td></tr><tr><td>W</td><td>Current volume</td></tr><tr><td>Ψ</td><td>Thermal energy</td></tr><tr><td>ḡ</td><td>Gravity acceleration vector</td></tr><tr><td>ψ̄</td><td>Heat flux vector</td></tr><tr><td>T</td><td>Temperature</td></tr><tr><td>C<sub>σ</sub></td><td>Specific heat capacity at constant stress</td></tr><tr><td>λ</td><td>Thermal conductivity</td></tr><tr><td>α<sub>th</sub></td><td>Linear thermal expansion coefficient</td></tr><tr><td>E</td><td>Young's modulus</td></tr><tr><td>ν</td><td>Poisson's ratio</td></tr><tr><td>α, k<sub>r</sub>, β, A, T<sub>r</sub>, χ<sub>1</sub>, χ<sub>2</sub></td><td>Viscoplastic rheological parameters (deviatoric part)</td></tr><tr><td>V, v, N, n, M, m</td><td>Viscoplastic rheological parameters (volumetric part)</td></tr></table>

# 1 Introduction

The time-dependent deformation of rock salt has been an important research topic in subsurface engineering, particularly in the fields of energy storage, mining, and nuclear waste disposal. Rock salt exhibits very low permeability and porosity in the undisturbed state, along with healing/sealing potential, ductility and hardening up to significant irreversible strain, and relatively high thermal conductivity as compared to other shallow-crustal rock types. These assets make rock salt an ideal geological barrier for containment, ensuring both environmental and operational safety (Bérest and Brouard 2003; Wolters et al. 2012; Hunsche and Hampel 1999; Bechthold et al. 1999; Stormont 1997; Hou 2003).

To predict the long-term behavior of rock salt surrounding underground facilities, it is crucial to use constitutive models calibrated on reliable data sets (phenomenological approach). These models must capture the complex

time-dependent behavior of rock salt when subjected to relevant thermal and mechanical loadings. However, it is not possible to establish a unique, comprehensive phenomenological constitutive model that accurately describes the behavior of a material such as rock salt (Noll 1974). In this regard, approximate laws are developed to closely reflect the actual behavior of the rock under specific loading conditions. This explains the large number of constitutive models for rock salt found in the literature.

Many existing models are based on experimental data obtained at deviatoric stresses exceeding 5 MPa and different temperatures (e.g., Bérest et al. (2019); Herchen et al. (2018); Lüdeling et al. (2022)). While these stress levels provide useful insights, they do not represent the full range of stresses encountered in underground openings, particularly those at lower stress levels (temperature is in general less extrapolated, except if very high frequency cycles are considered in salt cavern storage). Unfortunately, experimental challenges associated with low-stress tests limit the number of studies on rock salt creep in the low stress regime. The primary technical challenge in studying low differential stresses in the laboratory is that viscoplastic strains and strain rates are very small under such conditions (Bérest et al. 2008), requiring highly stable temperature and relative humidity (for unjacketed samples) to ensure accurate measurements. Even minor fluctuations in the applied stress can lead to significant changes in strain and strain rate, complicating the interpretation of results (see Blanco-Martín et al. (2023) and references therein). As a consequence, high-stress behavior is often extrapolated towards low-stress scenarios, which can lead to underestimations of strain rates in critical applications. In the case of salt caverns, accurate creep rate predictions are essential for understanding processes such as cavern convergence, pressurization, and the thermomechanical response far away from the salt cavern or during abandonment (van Gent et al. 2024).

Over the years, several models have been proposed to describe the creep behavior of rock salt, starting with early models like Norton-Hoff (Bérest et al. 2008), Munson-Dawson (Munson and Dawson 1981) or Lemaitre (Tijani et al. 1983), and evolving to incorporate more complex phenomena such as recovery strain (e.g., Munson (1997)), volumetric deformation, dilatancy, damage and healing (e.g., Cristescu (1993); Hunsche and Hampel (1999); Wolters et al. (2012); Labaune et al. (2018); Günther et al. (2015)). Despite these advancements, many models still rely heavily on empirical approximations and lack robust calibration at low-stress data (Lüdeling et al. 2024). Addressing this gap, Bérest et al. (2005, 2019, 2023) conducted uniaxial creep tests on natural salt samples under controlled temperature and humidity conditions in remote drifts in salt mines. In agreement with microstructural investigations by Spiers et al. (1990) and Urai et al.
(1986), their results suggest that there is a transition from high-n power-law to linear creep at low stress levels and that extrapolation of creep rates at moderate-high stress levels towards low levels leads to an underestimation of the creep rates by several orders of magnitude. However, in these studies low deviatoric stress levels were investigated with uniaxial tests, while high levels were commonly examined with confined tests. The tests were also performed at different temperatures in the different differential stress domains, and different scales of displacement measurements were employed. Additionally, even for tests using salt from the same origin, sample preparation and preconditioning were not always equivalent in the two domains. Despite these limitations, the experimental evidence underscores the importance of reducing stress extrapolation. More recently, Blanco-Martín et al. (2023, 2024) conducted multi-stage confined creep experiments on salt samples from the same origin and using similar preparation and preconditioning procedures, allowing for the study of a wider range of deviatoric stresses (between 0.5 and 20 MPa). Their results revealed marked differences in creep rates between high and low deviatoric domains. In parallel, recent modeling works have started to incorporate low-stress (linear or low-n power-law) components, improving predictions of the closure of underground openings and their long-term stability (e.g., Cornet et al. (2018); Hunfeld et al. (2022); Sobolik and Ross (2022); Reedlunn et al. (2022); Speirs et al. (2022)).

In addition to creep behavior, another critical factor is dilatancy, which is associated with micro-cracking and the development of flow paths (DeVries et al. 2002). Traditional stress-based dilatancy criteria (e.g., Van Sambeek et al. (1993); DeVries et al. (2005)), often based on (fast-rate) short-term tests, are fixed in the stress space and are not always consistent with the parameters derived from (slow-rate) long-term creep experiments (such criteria represent the stress state at the onset of dilatancy, expressed in terms of the equivalent deviatoric stress as a function of the mean pressure and possibly the Lode angle). Recent advancements, such as those proposed by Wolters et al. (2012), Lux et al. (2018), Labaune et al. (2018) and Rouabhi et al. (2019), have addressed this issue by integrating dilatancy directly into constitutive models. These enhanced models account for both short- and long-term conditions, improving the accuracy of stability predictions and providing a more comprehensive framework for understanding rock salt mechanical behavior under varying loading rates. Although progress has been made, most dilatancy criteria continue to disregard the microphysical mechanisms driving deformation, a key shortcoming that future research should tackle to extend our understanding of rock salt behavior beyond laboratory constraints.

This study examines the creep behavior of rock salt across a wide range of deviatoric stresses (relevant for underground applications), analyzing the impact of stress extrapolation for cavern stability and design. Additionally, we analyze the effects of incorporating the volumetric component of the viscoplastic strain tensor and compare dilatancy predictions with those of fixed stress-based criteria. Our findings reveal that, due to the highly nonlinear and stress-dependent nature of creep in rock salt, the viscoplastic response occurring at low stress levels is not accurately captured when extrapolating moderate-high stress behavior to low-stress conditions, leading to underestimations of the cavern volume loss and a misrepresentation of the magnitude of deviatoric stresses and strains. The impacts are visible during the cycling phase of a cavern, well before the abandonment phase. These insights are essential for improving long-term stability predictions of salt caverns under cyclic loading, contributing valuable data for underground energy storage and subsurface engineering in rock salt formations.

The structure of this paper is as follows. Section 2 details the numerical tools and governing equations considered in the simulations, along with the approach to solve the thermodynamic behavior within the cavern and the thermo-mechanical problem of the surrounding rock. This section also introduces the constitutive model, calibrated over a large range of deviatoric stresses and incorporating dilatancy, followed by a description of the generic scenario of cyclic hydrogen storage in a salt cavern and the rock salt properties used in the simulations. Section 3.1 outlines the criteria for evaluating cavern stability. Section 3.2 discusses the effects of stress extrapolation and their implications for stability, as well as the role of dilatancy. Section 3.3 provides insights into cavern spacing and cavern field design considerations. Finally, the key findings are summarized, and directions for future research are proposed in Sect. 4.

# 2 Materials and Methods

## 2.1 Numerical Framework

The numerical simulations were conducted using two in-house finite element method (FEM) codes: Demether (Rouabhi et al. 2017) and Phimef (Rouabhi 2022), developed at the Centre for Geosciences and Geo-engineering of Mines Paris. Demether is designed for the resolution of thermodynamic problems related to fluid storage in underground caverns, from the leaching phase to the abandonment phase. Fluid circulation along the well is also considered. In addition, the software is intended for the resolution of the thermo-mechanical problem of the rock mass around a spherical or a cylindrical cavern with finite
transformations and with strong coupling between the thermal and the mechanical components (i.e., a two-way coupled problem in which the output of one component serves as the input for the other until convergence is achieved). In Demether, the host rock is assumed to be a homogeneous, infinite medium with uniform initial temperature and uniform initial stress; such values are respectively the averages over the cavern surface of the geothermal temperature and the lithostatic stress. The mass and energy balance equations that govern the thermodynamical state of the stored fluid are solved assuming uniform pressure and temperature conditions (averaged over the cavern surface), and their temporal evolution can be expressed as:

$$ \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & -TA_{11} \end{pmatrix} \begin{pmatrix} \dot{T} \\ \dot{P} \end{pmatrix} = \begin{pmatrix} \dot{W} \\ \dot{\Psi} \end{pmatrix} \qquad (1) $$

where the values for the coefficients $A_{11}$, $A_{12}$ and $A_{21}$ depend on the current volume of the cavern and, respectively, on the isobaric thermal expansion coefficient, the isothermal compressibility factor and the isobaric specific heat capacity of each separate phase present in the cavern (e.g., stored fluid, residual brine). The term $\dot{W}$ is the current cavern volume variation, which is corrected by the variations in the partial volumes of the different species considered, and $\dot{\Psi} = \dot{\Psi}_{\sigma} + \dot{\Psi}_{l}$ is the heat exchanged with the rock salt formation and the heat resulting from the fluid mass entering the cavern (for further details, see Rouabhi et al. (2017)). The thermodynamic behaviour of the different species is taken into account using accurate equations of state (Rouabhi and Tijani 2014). For hydrogen, the GERG-2008 EOS is used (Kunz and Wagner 2012).

The thermo-mechanical problem for the rock salt around the cavern is solved with balance equations for momentum and energy, which read respectively

$$ \begin{gathered} \vec{\nabla} \cdot \underline{\underline{\sigma}} + \rho \vec{g} = \vec{0} \\ \rho C_{\sigma} \dot{T} + \vec{\nabla} \cdot \vec{\psi} = \vec{0} \end{gathered} \qquad (2) $$

where $\underline{\underline{\sigma}}$ is the Cauchy stress tensor, $\rho$ is the mass density, $\vec{g}$ the gravity acceleration vector, $C_{\sigma}$ is the specific heat capacity of the salt rock at constant stress, $\dot{T}$ is the material derivative of the temperature, $\vec{\psi} = -\lambda \vec{\nabla} T$ (Fourier's law) is the heat flux vector with thermal conductivity $\lambda$, and $\vec{\nabla} \cdot \vec{\psi}$ is the divergence operator. The mechanical constitutive model used will be described in Sects. 2.2 and 2.3. Boundary conditions at the cavern wall are provided by the thermodynamic state of the cavern, so that Eqs. (1) and (2) are solved iteratively within a time step until convergence is reached.

The second numerical tool used in this study, Phimef, is capable of solving coupled Thermal-Hydro-Mechanical (THM) processes in one, two, or three dimensions. Here, we consider only TM coupling, as its impact will be the

most important for the purposes of this work (note that we do not investigate the abandonment phase in the simulations, and also that the fluid flow mechanisms in rock salt should be better understood to properly implement THM coupling [e.g., entry pressure, Biot coefficient]). Phimef solves the momentum and energy balance equations (see Eq. (2)) using a finite strain framework. As compared to Demether, a more accurate solution is obtained for complex geometries and/or geological layers with different properties. Here, the simulations are 2D axisymmetric.

The workflow of our simulations is as follows. First, we use Demether to model the thermodynamic evolution of hydrogen inside the cavern and the thermo-mechanical response of the host rock. This preliminary model serves two purposes: (1) ensuring that the cavern dimensions meet well-established stability criteria (more details will be given in Sect. 3.1) using low computing-intensive simulations, and (2) generating the cavern pressure and temperature profiles that will be used as boundary conditions at the cavern wall during the thermo-mechanical simulations with Phimef (Phimef does not solve Eq. (1)). The leaching (assumed isothermal), debrining, and filling phases are modeled prior to the cycling phase (30 years) to account for the loading history and the stress state changes within the rock salt formation (Blanco-Martín et al. 2021). Next, we use Phimef to analyze the same scenario (considering a more realistic cavern geometry) and using the Demether-computed pressure and temperature histories at the cavern wall during cycling (boundary conditions applied at the cavern wall).

## 2.2 Constitutive Model

To solve the thermomechanical problem of the rock salt surrounding the cavern, we selected a constitutive law introduced by Labaune et al. (2018) that captures the time-dependent behavior of rock salt and reproduces both short- and long-term responses, using a unified set of parameters. The RTL model, which has been employed in previous studies on the mechanical behavior of rock salt (e.g., Azabou et al. (2021); Blanco-Martín et al. (2024)), serves as the foundation for the present analysis. It is based on the Lemaître model (Tijani et al. 1983) (which is purely deviatoric) and incorporates inelastic volumetric strain as a key feature (Labaune et al. 2018). It also provides a dilatancy criterion embedded at the constitutive level, which offers a more realistic alternative to the traditional criteria applied during post-processing (e.g., DeVries et al. (2005)). Let us consider the decomposition of the Cauchy stress tensor $\underline{\underline{\sigma}}$ into its deviatoric and spherical parts:

$$ \underline{\underline{\sigma}} = \sqrt{2/3} q J - \sqrt{3} p I \qquad (3) $$
where $q = \sqrt{3/2} ||\underline{\underline{\sigma}}'||$ is the von Mises equivalent stress, $J = \underline{\underline{\sigma}}'/||\underline{\underline{\sigma}}'||$ is a unit tensor in the direction of $\underline{\underline{\sigma}}'$, $I = \underline{\underline{1}}/\sqrt{3}$ ($\underline{\underline{1}}$ is the identity tensor), and $p = -\frac{1}{3}\text{tr}(\underline{\underline{\sigma}})$ is the mean pressure. Let us assume as well that the total strain rate tensor $\underline{\underline{\varepsilon}}$ can be decomposed into elastic ($\underline{\underline{\varepsilon}}_{el}$), thermal ($\underline{\underline{\varepsilon}}_{th}$) and viscoplastic ($\underline{\underline{\varepsilon}}_{vp}$) parts, as follows:

$$ \underline{\underline{\varepsilon}} = \underline{\underline{\varepsilon}}_{el} + \underline{\underline{\varepsilon}}_{th} + \underline{\underline{\varepsilon}}_{vp} \quad (4) $$

where $\underline{\underline{\varepsilon}}_{el} = \frac{(1+v)}{E}\underline{\underline{\sigma}} - \frac{v}{E}\text{tr}(\underline{\underline{\sigma}})\underline{\underline{1}}$ ($E$ is the Young's modulus and $v$ is the Poisson's ratio), $\underline{\underline{\varepsilon}}_{th} = \alpha_{th}\dot{T}\underline{\underline{1}}$ ($\alpha_{th}$ is the linear thermal expansion coefficient) and $\underline{\underline{\varepsilon}}_{vp}$ is given by a viscoplastic rheological model. Note that similarly to the stress tensor, the viscoplastic strain rate tensor can be decomposed into its distortional and volumetric components:

$$ \underline{\underline{\varepsilon}}_{vp} = \sqrt{3/2} \dot{\gamma}_{vp} \underline{\underline{N}} - \sqrt{1/3} \dot{\zeta}_{vp} \underline{\underline{I}} \quad (5) $$

where $\dot{\gamma}_{vp} = \sqrt{2/3} ||\underline{\underline{\varepsilon}}'||$ is the rate of the viscoplastic distortion (the deviatoric component of a tensor $\underline{\underline{a}}$ is denoted by $\underline{\underline{a}}' = \underline{\underline{a}} - \frac{1}{3}\text{tr}(\underline{\underline{a}})\underline{\underline{1}}$) and $\dot{\zeta}_{vp} = -\text{tr}(\underline{\underline{\varepsilon}}_{vp})$ is the rate of viscoplastic volumetric strain. Note that the unit tensor $\underline{\underline{N}}$ defines the direction of the deviatoric flow (for further details, see Labaune et al. (2018)).

The RTL model assumes that the past history of the material is represented by $\gamma_{vp}$, and also that $d\zeta_{vp}/d\gamma_{vp}$ depends only on the thermodynamic state, which is a function of $\underline{\underline{\sigma}}$, $T$ and $\gamma_{vp}$. Therefore, $\dot{\gamma}_{vp} = \dot{\gamma}_{vp}(\underline{\underline{\sigma}}, T, \gamma_{vp})$ and $\dot{\zeta}_{vp} = \varphi(\underline{\underline{\sigma}}, T, \gamma_{vp})\dot{\gamma}_{vp}$. The evolution of the viscoplastic distortion is governed by an extended Lemaitre model, adding the influence of stress triaxiality through the Lode angle, $\theta \in [0, \pi/3]$, defined as

$$ \theta = \frac{1}{3} \cos^{-1}(\ell), \quad \ell = \sqrt{6} \text{tr}(\underline{\underline{J}^3}) \quad (6) $$

Note that $\ell \in [-1, 1]$. The evolution law of the hardening variable is defined by

$$ \frac{d\gamma_{vp}^{1/\alpha}}{dt} = \left( \frac{q/\varphi(\theta)}{K(T)} \right)^{\beta/\alpha} \quad (7) $$

where $K(T) = k_r \exp(A \frac{\alpha}{\beta} \left( \frac{1}{T} - \frac{1}{T_r} \right))$ ($T_r$ is a reference temperature, at which $k_r$ is determined) and function $\varphi$ defines the shape of the deviatoric section of the deviatoric potential (Labaune et al. 2018):

$$ \varphi(\theta) = \varphi(\ell(\theta)) = \frac{\cos(\frac{\pi}{6}\chi_2 - \frac{1}{3}\arccos(-\chi_1))}{\cos(\frac{\pi}{6}\chi_2 - \frac{1}{3}\arccos(\chi_1\ell))} \quad (8) $$

Material parameters are, therefore, $\alpha, \beta, k_r, A, \chi_1$ and $\chi_2$ (all positive). Note that $0 < \alpha < 1$, $\chi_1 \in [0, 1]$ and $\chi_2 \in [0, 2]$.

Finally, the function $\varphi$ is based on experimental data (mostly derived from short-term triaxial tests conducted under constant strain rate or constant mean pressure conditions) and determines the evolution of the viscoplastic volumetric strain:

$$ \varphi(\sigma, T, \gamma_{vp}) = \min \left\{ v \frac{(\langle p \rangle / N)^n - \gamma_{vp}}{(\langle p \rangle / M)^m + \gamma_{vp}}, V \right\} \quad (9) $$

with positive parameters $V, v, N, n, M$ and $m$. Note that $\langle \cdot \rangle$ are the Macaulay brackets, such that $\langle x \rangle = (x + |x|)/2$. Eq. (9) defines a dilatancy threshold through parameters $N$ and $n$: when $\gamma_{vp} = (\langle p \rangle / N)^n$, the viscoplastic volumetric strain rate changes sign, transitioning between contracting and dilating behavior. That is, $\dot{\zeta}_{vp} = 0$ represents the onset of dilatancy, and $\dot{\zeta}_{vp} < 0$ corresponds to dilating conditions.

## 2.3 Model Fitting to Consider a Wide Range of Deviatoric Stresses

To demonstrate the consequences of extrapolating the behavior observed at high deviatoric stresses ($q \ge 5$ MPa) towards lower stress levels ($q < 5$ MPa), the deviatoric component of the RTL model was fitted using recent multi-stage confined creep tests, covering a differential stress range between 0.5 and 20 MPa (see for instance Blanco-Martín et al. (2024)). In that previous work, the authors proposed a modification to the Lemaitre law to better fit the experimental data (double power law of the deviatoric stress). Notwithstanding, in the generic scenario investigated here, we use the classic formulation of the Lemaitre law with a single power law of the deviatoric stress for the sake of simplicity (Eq. (7)). Even if a better fit of the experimental data could be captured with a double power law (showing the transition from power law to linear creep) at the expense of having more parameters, a single power law fitted over a large stress range can capture significant differences in the creep response at low deviatoric levels, see Fig. 1. As can be seen, the relative differences in the hardening parameter between the two fittings (data at $T = 323.15$ K) are more significant for low differential stresses than for high differential stresses. Note that the effect of temperature has been primarily investigated at high differential stresses ($T = 323.15$ K and $T = 293.15$ K), and should be evaluated in a larger, relevant deviatoric range
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Data analysis (Eq. (7)) showing the comparison between the two fittings for the deviatoric component of the RTL model in a differential stress range between 0.5 and 20 MPa (dashdotted line), and between 5 and 20 MPa (dashed lines). The effect of temperature is also displayed in the range 5-20 MPa. Parameters are: $\alpha = 0.326$, $\beta = 1.569$, $k_r = 0.065$ for the range 0.5-20 MPa, and $\alpha = 0.326$, $\beta = 2.496$, $k_r = 0.396$ for the range 5-20 MPa. $T_r = 323.15$ K, $A = 5,600$ K, $\chi_1 = \chi_2 = 0$. The units for the parameter $k_r$ are MPa·d$^{\alpha/\beta}$($\mu\text{m/m}$)$^{-1/\beta}$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Cavern geometry (2D axisymmetric model) and mesh refinement close to the cavern wall. Boundary conditions: no cross-boundary movement or thermal flux at the bottom and lateral boundaries; free surface at the top; fixed temperature at top and bottom boundaries (see text)</figcaption></figure>

in the future to confirm the formulation in Eq. (7) (for now, we assume the same temperature effect for low and high differential stresses).

Based on the fitting lines in Fig. 1, two simulations of underground hydrogen storage in salt caverns are presented in Sect. 3.2. The first simulation uses parameters calibrated from experimental data across the deviatoric range 0.5-20 MPa. This case represents the ideal case, in which experiments and models consider a large, relevant range of deviatoric stresses. The second simulation uses parameters fitted using data limited to the upper deviatoric range, 5-20 MPa. This case represents current standard practice, with experiments being conducted at high deviatoric stresses, constitutive models being fitted on those data, and their predictions being extrapolated to lower deviatoric values. Note that although the first case reduces stress extrapolation, uncertainties persist beyond the range of the experimental data because the constitutive model used is phenomenological.

## 2.4 Modeling Scenario

For illustrative purposes, we have considered a generic operational scenario of hydrogen storage in salt caverns. In this context, integrity issues arise over time due to frequent pressure and temperature fluctuations during injection and withdrawal cycles, which can affect cavern stability (Bérest et al. 2012; Anderson et al. 2022). Therefore, it is crucial to

investigate the thermomechanical response of rock salt to the rapid variations in the cavern occurring under fast cyclic loading conditions, as seen in short-term daily or weekly demand scenarios.

The modeling scenario presented here contemplates 7 days of withdrawal followed by 7 days of injection over a 30-year cycling period. The cylindrical-shaped cavern is located at a depth of 1,000 m (measured at the cavern mid-height) and has radius of 32 m, a height of 100 m, and a volume of 153,600 m³. The Demether computations assume a perfect cylindrical shape (see Sect. 2.1). To avoid stress concentrations at the corners and to comply with standard solution mining practice, in Phimef the cavern shape is modified at the top and bottom as shown in Fig. 2. The salt **body** extends vertically to the ground surface and has a lateral extension of 100 times the cavern radius.

The design parameters are based on engineering calculations from Hadj-Hassen et al. (2020) and optimized to meet the energy demands of a mid-sized city with a renewable electricity consumption of 200 MW. The withdrawal/injection rate is 232,258 kg/day, with an injection pressure of 13.75 MPa and an injection gas temperature of 30 °C. The ground surface temperature is 15 °C, and a geothermal gradient of 0.025 °C/m is considered. The initial stress state is assumed isotropic and equal to the lithostatic gradient (generic scenario with no particularities regarding salt tectonics, so that the focus is kept on the role of stress extrapolation and dilatancy). Accordingly, the initial in situ stress and temperature are -21.58 MPa and 40 °C, respectively, at cavern-mid depth (compression is negative). The leaching
<figure><img src="image_5.png" /><figcaption>Fig. 3 Imposed pressure at the cavern wall</figcaption></figure>

<table><caption>**Table 1** Properties of the rock salt used for the simulations</caption><tbody><tr><th colspan="4">Thermo-physical parameters</th></tr><tr><td>ρ (kg/m<sup>3</sup>)</td><td>C<sub>σ</sub> (kJ/kg/K)</td><td>λ (MJ/d/m/K)</td><td></td></tr><tr><td>2,200</td><td>0.8</td><td>0.499</td><td></td></tr><tr><th colspan="4">Thermo-elastic parameters</th></tr><tr><td>α<sub>th</sub> (K<sup>-1</sup>)</td><td>E (MPa)</td><td>ν (-)</td><td></td></tr><tr><td>25</td><td>22,000</td><td>0.3</td><td></td></tr><tr><th colspan="4">RTL, deviatoric part fitted on 5 ≤ q ≤ 20MPa</th></tr><tr><td>α</td><td>k<sub>r</sub></td><td>β</td><td></td></tr><tr><td>0.326</td><td>0.396</td><td>2.496</td><td></td></tr><tr><th colspan="4">RTL, deviatoric part fitted on 0.5 ≤ q ≤ 20MPa</th></tr><tr><td>α</td><td>k<sub>r</sub></td><td>β</td><td></td></tr><tr><td>0.326</td><td>0.065</td><td>1.569</td><td></td></tr><tr><th colspan="4">RTL, deviatoric part (additional parameters)</th></tr><tr><td>A</td><td>T<sub>r</sub></td><td>χ<sub>1</sub></td><td>χ<sub>2</sub></td></tr><tr><td>5,600</td><td>293.15</td><td>0</td><td>0</td></tr><tr><th colspan="4">RTL, volumetric part</th></tr><tr><td>V</td><td>ν</td><td>n</td><td>N m M</td></tr><tr><td>-1</td><td>0.11</td><td>2.34</td><td>0.24 4.99 2.80</td></tr></tbody></table>

For RTL, the unit system is such that strain is in µm/m, stress in MPa, temperature in K and time in days

stage, which lasts 1,174 days, is modeled by linearly reducing the cavern pressure ($P_{cav}$) from the lithostatic pressure ($P_{lith}$) to the halmostatic pressure ($P_{brine}$) at each cavern depth, followed by 130 days of debrining and first filling, during which $P_{cav}$ increases to the maximum cavern pressure ($P_{max}$). From that point, the (uniform) cavern pressure oscillates between $P_{max}$ and $P_{min}$. At the beginning of the cycling phase, it is assumed that the temperature of the fluid ($T_{fluid}$) inside the cavern has reached thermal equilibrium with the surrounding rock mass ($T_{∞}$). Fig. 3 illustrates the pressure history applied at the cavern wall and Table 1 lists the parameters used in the simulations (note that the temperature effect requires further experimental investigation as explained above). Constitutive parameters of the deviatoric part of the RTL model are obtained as described in Sect. 2.3

and those of the volumetric part are obtained from strain-controlled tests at $2·10^{-6}$ s$^{-1}$ and confining pressures in the range 5-20 MPa (Blanco-Martín et al. 2022). Boundary conditions are as follows. The bottom and lateral boundaries of the domain have no-flow and no-normal displacement. Displacement is allowed in the upper boundary, and the temperature is imposed on the top and bottom boundaries (15 and 42.5 °C, respectively).

# 3 Results and Discussion

## 3.1 Stability Criteria

To establish a robust basis for evaluating the numerical results it is crucial to define the necessary stability criteria for the underground facilities. For salt cavern storage, there are well-established criteria that have been used for the last decades (e.g., Heusermann et al. (2003); Bérest (2013); Bérest et al. (2020); Nguyen Minh et al. (1993); Minkley et al. (2015)). A first important criterion stipulates that the internal pressure of the cavern should be maintained between 20% and 80% of the lithostatic stress in absolute value ($σ_{∞}$) to respectively limit excessive deviatoric stresses and prevent both micro-fracture formation and fluid leakage through the cavern wall. The results show that $P_{min}$ ranges between 5.2 and 5.9 MPa, while $P_{max}$ fluctuates between 13.9 and 15.9 MPa. These values are within the predefined stability limits of 0.2 $σ_{∞}$ and 0.8 $σ_{∞}$, as depicted in Fig. 4 (bottom curves). Note that the overall temperature trend (top curves) decreases over time due to the injection temperature (30 °C), which is lower than $T_{∞}$ (40 °C). Conversely, the overall increase in pressure over time can be attributed primarily to cavern closure due to the creep of the host rock. As can be seen, the temperature range falls within the experimental range used to fit the temperature effect in Eq. (7) (20-50 °C, see Fig. 1, although more data points are required at low differential stresses to confirm the trends).

The second stability criterion requires that the stress state of the rock mass surrounding the cavern remains compressive to avoid the formation of tensile stresses that could surpass the low tensile strength of the rock salt host rock. The simulation results show that the least compressive principal stress ($σ_1$) remains negative (indicating compression) and becomes increasingly compressive over time as illustrated in Fig. 5. The evolution of $σ_1$ is similar for the two parameter sets in Table 1, and the temperature variations (see Fig. 4) do not induce any tensile stresses. In addition, following Blanco Martín et al. (2015) and references therein, we evaluated whether the applied pressure at the cavern wall ($P_{cav}$) locally exceeds the least compressive principal stress $σ_1$, which would indicate conditions favorable to fluid infiltration (Fokker 1995; Kenter et al. 1990; Wolters et al. 2012).
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Evolution of the least compressive principal stress at the cavern wall during the cycling phase. For interpretation purposes, the left plot shows only the maximum and minimum values of each cycle. The right plot is a close-up to the cycles during the first year</figcaption></figure>

<figure><img src="image_7.png" /></figure>

Indeed, in these past works hydraulic fracturing is assumed to occur if $\sigma_1 + P_{cav}$ exceeds a critical threshold $\sigma_T$, typically ranging from 2 to 4 MPa. During the cycling period, in our simulations $\sigma_1 + P_{cav}$ remains close to 0 MPa, indicating that hydraulic fracturing conditions are not reached within the simulations time frame.

Another crucial stability criterion is the change in cavern volume over time. Although all solution-mined caverns and mined drifts tend to decrease in volume gradually due to

creep behavior, a closure rate of up to 1% per year is considered acceptable for maintaining stability. This gradual volume reduction is important not only because it affects the storage capacity of the cavern, but also due to its potential to cause surface subsidence. In the next section, we analyze the volume loss of the cavern while considering the effects of stress extrapolation towards low deviatoric stress regimes.

Finally, the development of dilatancy zones around the openings is an important criterion for the operation of

<figure><img src="image_8.png" /></figure>
underground facilities. As explained before, dilatancy (non-
elastic volumetric strain), is related to the development of
micro-fracturing and possibly flow paths, and its extent
around the excavations should be limited, if any. In Sect. 3.2
we will compare the predictions of dilatancy from the crite-
ria implemented at the constitutive level (they evolve as the
stresses and strains evolve), with predictions from criteria

fixed in the stress space and used during post-processing
(they define a boundary in the stress space between the con-
tracting and dilating behavior of rock salt).

**3.2 Impact of Stress Extrapolation and Dilatancy**

Figure 6 illustrates the evolution of the cavern contour at different stages: the start of the simulation (0 years), the end of leaching (1174 days), the beginning of cyclic loading (1303 days), and after 1, 10 and 30 years of cycling. Starting from the leaching stage, the model fitted on the wider range of deviatoric stresses predicts greater cavern convergence than the model fitted only at high deviatoric stresses. After 30 years of cycling, the maximum horizontal displacements (negative, indicating closure) at the cavern wall predicted by the two simulations are 2.4 m and 1.2 m, respectively. Despite the differences in the magnitude of creep deformation, both cases show intensified convergence with depth (compare the deformation at the top of the cavern with that at the bottom), highlighting the critical role of deviatoric stresses in driving the creep behavior of rock salt.

As previously discussed, a closure rate of up to 1% per
year is considered acceptable for maintaining long-term
cavern stability. Fig. 7 shows that the modeled annual clo-
sure rates remain well below this threshold in both cases,
demonstrating compliance with the stability criterion. By
the end of the cycling phase, the case fitted only at high
deviatoric stresses predicts a volume loss of approximately
-8%, whereas the wider-range case forecasts a volume loss
close to -15%. This greater volume loss can be attributed

<figure><img src="image_9.png" /><figcaption>**Fig. 6** Evolution of the cavern contour over 30 years for the model calibrated on experimental data within the range $5 \le q \le 20$ MPa (left) and the model fitted in the range $0.5 \le q \le 20$ MPa (right). An amplification coefficient of 10 has been applied to highlight the geometry effects on the cavern convergence</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7** Volume loss of the cavern over time. For the model calibrated on experimental data within the range $0.5 \le q \le 20$ MPa, the volume loss is calculated with and without activating the volumetric component of $\dot{\epsilon}_{\parallel} = \nu_p$ (parameter $\nu$ in Eq. (9))</figcaption></figure>
to the higher levels of distortion occurring at low deviatoric stresses (refer to Fig. 1). Similar to Fig. 6, this figure shows that stress extrapolation is important since an early stage of a cavern lifecycle.

An additional simulation using the wider-range param-
eter set was conducted to evaluate whether the activation
of the volumetric component of the viscoplastic strain rate
tensor influences the cavern volume loss. In this scenario,
the parameter v = 0 was set in Eq. (9) to exclude the volu-
metric component. When the volumetric component is dis-
regarded, the predicted volume loss at the end of cycling is –
15.13%, compared to –15% when the volumetric component
is included. This negligible difference aligns with the results
obtained by Labaune et al. (2018), although our analysis
is thermomechanical and accounts for greater viscoplastic
strains at low deviatoric stresses. As can be seen, adding
a volumetric component to the viscoplastic strain rate ten-
sor has minimal impact on the predicted volume loss. This
observation also supports the notion that cavern closure
is not primarily governed by the dilatancy of the rock salt
(mostly confined to the near field as will be seen below), but
is largely influenced by deviatoric stresses promoting the
creep of the surrounding rock mass.

Underground excavations disrupt the initial stress equi-
librium of the rock formation, leading to a redistribution
of deviatoric stresses around the openings. To assess how
the stress distribution changes when the creep response at
low deviatoric stresses is extrapolated, Fig. 8 presents the
profiles of *q* at the end of the last withdrawal stage. The dis-
tribution of deviatoric stresses differs between the two cases.
For the wider-range parameter set, the region influenced by
very low deviatoric stresses is significantly larger, due to the
faster stress evolution at higher creep rates. This is consist-
ent with our previous observations of greater volume loss
associated with higher viscoplastic strains at low stresses.
Additionally, higher deviatoric stresses are observed at the

cavern wall and near the bottom of the cavern. These results
are consistent with those reported by Lüdeling et al. (2024),
further supporting that the accelerated creep rate of the rock
mass surrounding the cavern increases the stress magnitude
at the cavern wall as a result of the faster closure rate. Bear-
ing Fig. 1 in mind, Fig. 8 reflects that the higher cavern
convergence in the wider-range case is mostly due to the
far field: for both parameter sets, *q* is large at the cavern
wall and small in the far field, but more so for the wider-
range case. At higher values of *q*, γ̇<sub>vP</sub><sup>1/α</sup> is slightly higher for
the model fitted on the upper range of deviatoric stresses;
however, for lower values of *q*, γ̇<sub>vP</sub><sup>1/α</sup> is much higher for
the wider-range parameter set. Note that the wider-range
parameter set predicts minimum deviatoric stresses around
0.004 MPa, which are below the experimental range shown
in Fig. 1 (on the other hand, *q*<sub>max</sub> < 20 MPa). Therefore,
some stress extrapolation is still present, and would require
additional tests at very low deviatoric stresses under highly
controlled experimental conditions. At the same time, more
research is needed on the creep of rock salt in this stress
range (experiments and microstructural observations).
Indeed, recent studies by Hunfeld et al. (2022) and van Oost-
erhout et al. (2022) suggest the existence of a differential
stress threshold (approximately 0.07–0.7 MPa) below which
the creep response may differ as a result of a change in the
dominant deformation mechanism. Note that this hypothesis
still requires microstructural validation of grain boundary
behavior, and confined creep tests at very low deviatoric
stress to confirm if the stress threshold should be accounted
for in future numerical models.

As previously discussed, dilatancy should be avoided
or minimized to prevent the formation of flow paths and
the associated potential increase in porosity and permea-
bility. To analyze the occurrence and distribution of dila-
tancy zones, we refer to the profiles in Fig. 9. These iso-
lines represent the loci of ζ̇<sub>vP</sub> = 0, indicating the onset of

<figure><img src="image_11.png" /><figcaption>**Fig. 8** Distribution of the von Mises equivalent stress for the last day of cycling (a withdrawal stage). Model fitted on experimental data in the range 5 ≤ *q* ≤ 20 MPa (left) and model fitted in the range 0.5 ≤ *q* ≤ 20 MPa (right). The dashed-dotted line indicates the initial cavern contour at the start of the cycling phase</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 9** Limits of the dilatancy regions at different moments of the simulation. The isolines represent the loci of $\zeta_{vp} = 0$ (see Eq. (9))</figcaption></figure>

dilatancy. At a given time, dilatancy occurs within the region bounded by the cavern wall and the corresponding isoline. Note that the distance between isolines equally spaced in time (e.g., from +10 years to +30 years) decreases, demonstrating that the extent of the dilatancy-affected zone stabilizes over time and does not continue to expand indefinitely (this is also consistent with Labaune et al. (2018)). For the upper-range parameter set (fitted within $5 \le q \le 20$ MPa), dilatancy begins one month after the cycling phase starts, primarily affecting the region near the cavern midsection and the bottom of the cavern, with the lateral width of the dilatancy zone reaching approximately 20 m (from the cavern midsection) after 30 years of cycling. In comparison, the wider-range parameter set (fitted within $0.5 \le q \le 20$ MPa) shows dilatancy occurring after 14 days of cycling, with the width of the dilatancy zone reaching 26 m by the end of the cycling period. These results are obtained with the parameters of the volumetric part of RTL displayed in Table 1 (which were calibrated on triaxial tests performed at constant strain rates as explained above), and might be considered conservative with a view to the width of the predicted dilatancy zones (notwithstanding, field evidence is hard to obtain from convergence or fluid mass balance estimates). Here, our purposes are only illustrative rather than site-specific. We highlight that the wider-range set predicts an earlier onset of dilatancy and a larger damage zone compared to the model calibrated solely at high deviatoric stresses (recall that $\zeta_{vp} = \varphi(\sigma, T, \gamma_{vp})\dot{\gamma}_{vp}$). Interestingly, the differences in the spatiotemporal extent of dilatancy between the two models are more limited than those displayed in Figs. 7-8.

<figure><img src="image_13.png" /><figcaption>**Fig. 10** Evolution of the deviatoric stress vs. the mean pressure at the bottom of the cavern wall (where the deviatoric stresses are higher). Note that no dilatancy occurs according to the stress-based criterion $q = 1.23p + 8.18$ in any of the two scenarios</figcaption></figure>

Next, for the purposes of illustration, we consider a linear stress-based dilatancy criterion following Labaune et al. (2018). For the rock salt considered in Table 1, dilatancy occurs at points in stress space that satisfy $q = 1.23p + 8.18$ (expressed in MPa). As shown in Fig. 10, no dilatancy is predicted for either of the two cases analyzed using this stress-based criterion, which is in clear contrast with the results obtained when dilatancy is considered at the constitutive level (Fig. 9). These results highlight the limitations of stress-based criteria, derived from the onset of dilatancy during short-term triaxial tests, and analyzed independently of factors such as the loading rate. In contrast, the dilatancy criteria embedded at the constitutive level evolve concurrently with the stresses and strains. The latter formulation, such as that described in Sect. 2.2, is believed to be more realistic than a fixed post-processing approach based on discrete experimental data points that lack direct strain-rate correspondence.

## 3.3 Spatial Extent of Dilatancy vs. Creep at Low-Stress Levels and Implications for Cavern Spacing

In Fig. 9, we showed that the damage zone extends less than one cavern radius after 30 years of cycling when dilatancy is incorporated at the constitutive level, which is well within the typical engineering requirement for pillar width, usually two to three times the cavern diameter. It is useful to compare the lateral extent of the dilatancy-affected area with the lateral extent of the region subjected to very low deviatoric stresses in both cases (model fitted on both low-high stress deviators and model fitted only on high stress
<figure><img src="image_14.png" /><figcaption>**Fig. 11** Profile at the end of the last withdrawal stage showing the equivalent deviatoric stress as a function of the radial distance to the cavern wall. The inset zooms in on the region in the far field affected by very low deviatoric stresses, emphasizing that these stresses do not reach zero</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12** Profile at the end of the last withdrawal stage showing the magnitude of the viscoplastic distortion as a function of the radial distance to the cavern wall. The zoomed-in subplot emphasizes the differences in the magnitude of the viscoplastic distortion predicted by the two parameter sets in the far field</figcaption></figure>

deviators), to provide more insights for the design of a multi-cavern configuration (thus, cavern spacing).

Figure 11 shows the evolution of the deviatoric stress at cavern mid-height as a function of the lateral distance from the cavern wall. The stress at the cavern wall is larger when the parameter set fitted on $0.5 \le q \le 20$ MPa is used, and decreases with the distance to the cavern wall in both cases (note the higher decay rate when the wider-range parameter set is considered). Observe that, despite the similar trends, the wider-range case predicts lower values of $q$ from about 40 m due to enhanced creep. Regardless of the parameter set used, the lateral extent of the region subjected to very low deviatoric stresses is notably larger than the width of the dilatancy zones in Fig. 9 (20–26 m). Note that $q > 0$ far away from the cavern wall. In the far field, and for a given deviatoric stress level such as $q = 1$ MPa, the extent of rock volume affected by low deviatoric stresses varies significantly between the two parameter sets. For the wider-range set, this region extends 138 m from the cavern wall, whereas it extends beyond 234 m for the upper-range set. This difference suggests that reducing stress extrapolation can effectively identify rapid spatial decreases of deviatoric stresses in low-stress zones, thereby improving the accuracy of stress field predictions. These factors are critical to consider in the design process, as they imply faster deviatoric stress reduction accompanied by accelerated creep rates.

Now, let us examine the evolution of viscoplastic distortion as a function of the radial distance to the cavern wall. As shown in Fig. 12, stress extrapolation towards low $q$ leads to underestimation. At the cavern wall, where high deviatoric stresses dominate (see Figs. 8 and 11), viscoplastic distortion after 30 years of cycling reaches about 9.4% with the wider-range parameter set, compared to 5% with the

upper-range set. Although the magnitude of the viscoplastic distortion decreases with distance in both models, the differences remain significant. At 100 m from the cavern wall, the viscoplastic distortion with the wider-range parameter set is about twice that with the upper-range set, at 0.39% and 0.20%, respectively. Similar to the area affected by very low deviatoric stresses, the region influenced by viscoplastic distortion extends farther than the dilatancy-affected damage zone, which remains confined at closer distances from the cavern wall.

# 4 Conclusions and Perspectives

The thermomechanical simulations presented in this study build on previous work by delving deeper into the effects of stress extrapolation and dilatancy on predicting the creep behavior of rock salt. Our findings underline some problems associated with neglecting the effects of low differential stress regimes, offering insights that could lead to optimized storage designs. These numerical results reveal that relying on stress extrapolation can lead to erroneous estimations of cavern volume loss, stresses and strains. In this sense, they enhance our understanding of how rock salt properties influence spatiotemporal stress and strain distributions. By reducing stress extrapolation, we obtain more accurate predictions of the effect of the accelerated closure rates at low deviatoric stresses, with direct implications on stability, performance and storage capacity. From a design perspective, the consideration of a wider range of deviatoric stresses could allow for more efficient engineering decisions on cavern spacing (i.e.,

<figure><img src="image_16.png" /></figure>
pillar widths). Indeed, we provided additional evidence that supports the dominance of creep under low deviatoric stresses in the far field, leading to enhanced cavern closure effects. An important conclusion of this work is therefore that the laboratory experiments on which creep laws are fitted should be performed as much as possible on the thermomechanical ranges (e.g., stresses, temperature, strain and strain rates) that are relevant for a particular underground application (obtained from modeling and/or from field data), so that extrapolation is reduced. This normally requires combining different types of experiments (e.g., triaxial, creep, tensile) to fit the constitutive laws. Here, we focused on extrapolation towards low stress levels, but the same applies if the deviatoric stresses expected in the field are higher than the ones investigated in the laboratory, and also to other variables such as temperature and strain rates (very small in the far field, but orders of magnitude larger close to the cavern wall during cycling).

Dilatancy criteria incorporated at the constitutive level and applied across a broad range of deviatoric stresses enable more precise predictions of the onset of dilatancy and the extent of the damaged zone. Notably, this affected zone is limited compared to the larger region influenced by very low deviatoric stresses, emphasizing the need for enhanced creep laws that can reliably capture the creep response in all relevant ranges for the field. These findings suggest that current design practices could benefit from a reassessment that incorporates a broader spectrum of differential stresses, as well as addresses the limitations of stress-based dilatancy criteria used during post-processing of simulation results.

Building on these findings, several avenues for future research stem from the unresolved questions about extrapolation and dilatancy. Future work will incorporate a more comprehensive Thermal-Hydraulic-Mechanical framework, including the effects of fluid flow through rock salt (fluid permeation in the dilated zone accounting for a possible distinction between total and effective stresses) and predictions for the cavern abandonment phase. The constitutive models will be calibrated on data from laboratory experiments conducted, whenever feasible, within stress, temperature, strain and strain rate ranges relevant to the specific underground application, and complemented with microphysical considerations for creep and dilatancy to constrain extrapolations beyond laboratory-accessible conditions. Future research could also consider heterogeneous rock salt formations, a broader variety of cavern geometries and operational conditions. Additionally, calibration of the models with long-term field data would enhance the reliability of the numerical predictions, addressing challenges associated with time extrapolation and further clarifying rock salt behavior for safer, more effective underground storage.

**Acknowledgements** Funding for this study has been provided by the ANR (French National Research Agency) under contract number ANR-21-CE05-0027-05 (RockStorHy project). The authors express their gratitude to the project partners (Laboratoire de Mécanique des Solides - École Polytechnique, Laboratoire Navier-ENPC, Laboratoire Angevin de mécanique, procédés et innovation - ENSAM and Brouard Consulting) for the constructive discussions which have helped improve the quality of the work presented here.

**Author Contributions** JAJC, LBM, AR and FHH designed the study. JAJC and LBM performed the numerical simulations, analyzed the data and prepared the manuscript with contributions from all co-authors.

**Funding** Open access funding provided by Mines Paris - PSL. Funding for this study has been provided by the ANR (French National Research Agency) under contract number ANR-21-CE05-0027-05 (RockStorHy project).

**Data Availability** Simulation files can be made available upon request.

## Declarations

**Conflict of interest** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Anderson S, Barnes L, Doyle R, Frey F, Garcia E, Laszlo J, Marshall M, Ozah R (2022) Solutions for storage of green hydrogen in salt caverns. In: SMRI Spring Technical Conference, Rapid City, SD, USA

Azabou M, Rouabhi A, Blanco-Martín L, Hadj-Hassen F, Karimi-Jafari M, Hévin G (2021) Rock salt behavior: from laboratory experiments to pertinent long-term predictions. Int J Rock Mech Min Sci 142:104588. https://doi.org/10.1016/j.ijrmms.2020.104588

Bérest P, Brouard B (2003) Safety of salt caverns used for underground storage blow-out; mechanical instability; seepage; cavern abandonment. Oil and Gas Science and Technology-revue De L'Institut Français Du Petrole: Oil Gas Sci Technol 58, 361-384 https://doi.org/10.2516/ogst:2003023

Bérest P, Blum P, Charpentier J-P, Gharbi H, Valès F (2005) Very slow creep tests on rock samples. Int J Rock Mech Min Sci 42(4):569-576. https://doi.org/10.1016/j.ijrmms.2005.02.003

Bérest P, Brouard B, Karimi M, Réveillère A (2020) Maximum admissible pressure in salt caverns used for brine production and hydrocarbon storage. Oil Gas Sci Technol Revue de l'IFP. https://doi.org/10.2516/ogst/2020068
Bérest P, Djizanne H, Brouard B, Hévin G (2012) Rapid depressurizations: can they lead to irreversible damage? In: Proceedings SMRI Spring Technical Conference, Regina, SK, Canada

Bérest P (2013) The mechanical behavior of salt and salt caverns. In: ISRM international symposium - EUROCK 2013, Wrocław, Poland, pp 17–29. https://doi.org/10.1201/b15683-5

Bérest P, Gharbi H, Brouard B, Brückner D, DeVries K, Hévin G, Hofer G, Spiers CJ, Urai JL (2019) Very slow creep tests on salt samples. Rock Mech Rock Eng 52(9):2917–2934. https://doi.org/10.1007/s00603-019-01778-9

Bérest P, Gharbi H, Blanco-Martín L, Brouard B, Brückner D, DeVries K, Hévin G, Hofer G, Spiers CJ, Urai JL (2023) Salt creep: transition between the low and high stress domains. Rock Mech Rock Eng. https://doi.org/10.1007/s00603-023-03485-y

Bérest P, Karimi-Jafari M, Brouard B (2008) Comportement mécanique à très long terme des mines et cavernes dans le sel gemme: loi de Norton-Hoff ou loi de Lemaître? Revue Française de Géotechnique 124, 45–59. https://doi.org/10.1051/geotech/2008124045

Blanco-Martín L, Azabou M, Rouabhi A, Hadj-Hassen F, Hévin G (2022) Insights from an experimental campaign on rock salt and implications at the cavern scale. In: de Bresser, Drury, Fokker, Gazzani, Hanx, Niemeijer, Spiers (eds.) 10^th^ Conference Mechanical Behavior of Salt (Saltmech X), Utrecht, The Netherlands, pp 119–129

Blanco-Martín L, Azabou M, Rouabhi A, Hadj-Hassen F (2023) Laboratory experiments on rock salt and phenomenological observations. Int J Rock Mech Min Sci 170:105452. https://doi.org/10.1016/j.ijrmms.2023.105452

Blanco-Martín L, Jiménez-Camargo JA, Gharbi H, Dimanov A, Bornert M, Brouard B (2023) Rock salt creep: a cross-check test to cover the full relevant range of deviatoric stresses. Rev Fr Géotech. https://doi.org/10.1051/geotech/2024017

Blanco-Martín L, Rouabhi A, Hadj-Hassen F (2021) Use of salt caverns in the energy transition: application to Power-to-Gas-Oxyfuel. J Energy Storage 44:103333. https://doi.org/10.1016/j.est.2021.103333

Blanco-Martín L, Rouabhi A, Hadj Hassen F, Jaworowicz J, Azabou M, Hévin G, Labaune P (2024) Creep of rock salt under a large range of deviatoric stresses. Rock Mech Rock Eng. https://doi.org/10.1007/s00603-024-03841-6

Blanco-Martín L, Rutqvist J, Birkholzer JT (2015) Long-term modeling of the thermal-hydraulic-mechanical response of a generic salt repository for heat-generating nuclear waste. Eng Geol 193:198–211. https://doi.org/10.1016/j.enggeo.2015.04.014

Bechthold W, Rothfuchs T, Poley A, Ghoreychi M, Heusermann S, Gens A, Olivella S (1999) Backfilling and sealing of underground repositories for radioactive waste in salt, Phase I (BAMBUS project). Technical report EUR 19124 EN, Nuclear Science and Technology Commission of the European Communities

Cornet JS, Dabrowski M, Schmid DW (2018) Long term creep closure of salt cavities. Int J Rock Mech Min Sci 103:96–106. https://doi.org/10.1016/j.ijrmms.2018.01.025

Cristescu ND (1993) A general constitutive equation for transient and stationary creep of rock salt. Int J Rock Mech Min Sci Geomech Abstr 30(2):125–140. https://doi.org/10.1016/0148-9062(93)09070-5

DeVries K, Mellegard K, Callahan G (2002) Salt damage criterion proof-of-concept research. Technical Report DE-FC2600NT41026, RESPEC, Rapid City, SD

DeVries K, Mellegard K, Callahan G, Goodman W (2005) Cavern roof stability for natural gas storage in bedded salt. Technical report DE-FG26-02NT41651, RESPEC, Rapid City, SD

Fokker PA (1995) The behavior of salt and salt caverns. Ph.d. thesis, Delft University, Delft, The Netherlands

Günther R-M, Salzer K, Popp T, Lüdeling C (2015) Steady-state creep of rock salt: improved approaches for lab determination and

modelling. Rock Mech Rock Eng 48(6):104452. https://doi.org/10.1007/s00603-015-0839-2

Hunfeld L, Breunese J, Wassing B (2022) The influence of a threshold stress for pressure solution creep on cavern convergence and subsidence behavior - A FEM study. In: de Bresser, Drury, Fokker, Gazzani, Hanx, Niemeijer, Spiers (eds.) 10^th^ Conference Mechanical Behavior of Salt (Saltmech X), Utrecht, The Netherlands, pp 577–589

Hunsche U, Hampel A (1999) Rock salt - the mechanical properties of the host rock material for a radioactive waste repository. Eng Geol 52(3):271–291. https://doi.org/10.1016/S0013-7952(99)00011-3

Hadj-Hassen F, Abuaisha M, Rouabhi A (2020) Simulation de scénarios types de stockage en cavités salines (in French). Mines Paris - Centre de Géosciences

Hou Z (2003) Mechanical and hydraulic behavior of rock salt in the excavation disturbed zone around underground facilities. Int J Rock Mech Min Sci 40(5):725–738. https://doi.org/10.1016/S1365-1609(03)00064-9

Herchen K, Popp T, Düsterloh U, Lux K-H, Salzer K, Lüdeling C, Günther R-G, Rölke C, Minkley W, Hampel A, Yildirim S, Staudtmeister K, Gährken A, Stahlmann J, Reedlunn B, Hansen FD (2018) WEIMOS: Laboratory investigations of damage reduction and creep at small deviatoric stresses in rock salt. In: Fahland H, Hansen H, Lux M (eds) 9^th^ conference mechanical behavior of salt (Saltmech IX). Hanover, Germany, pp 175–192

Heusermann S, Rolfs O, Schmidt U (2003) Nonlinear finite-element analysis of solution mined storage caverns in rock salt using the LUBBY2 constitutive model. Comput Struct 81(8):629–638. https://doi.org/10.1016/S0045-7949(02)00415-7

Kenter C, Doig S, Rogaar H, Fokker P, Davies D (1990) Diffusion of brine through rock salt roof of caverns. In: Proceedings of SMRI fall meeting, Paris

Kunz O, Wagner W (2012) The GERG-2008 wide-range equation of state for natural gases and other mixtures: an expansion of GERG-2004. J Chem Eng Data 57(11):3032–3091. https://doi.org/10.1021/je300655b

Lüdeling C, Günther RM, Hampel A, Sun-Kurczinski J, Wolters R, Düsterloh U, Lux K-H, Yildirim S, Zapf D, Wacker S, Epkenhans I, Stahlmann J, Reedlunn B (2022) WEIMOS: creep of rock salt at low deviatoric stresses. In: de Bresser, Drury, Fokker, Gazzani, Hanx, Niemeijer, Spiers (eds.) 10^th^ conference mechanical behavior of salt (Saltmech X), Utrecht, The Netherlands, pp 130–140

Lux K-H, Lerche S, Dyogtyev O (2018) Intense damage processes in salt rocks—a new approach for laboratory investigations, physical modelling and numerical simulation. In: 9^th^ Conference Mechanical Behavior of Salt (Saltmech IX), Hanover, Germany, pp 47–67

Labaune P, Rouabhi A, Tijani M, Blanco-Martín L, You T (2018) Dilatancy criteria for salt cavern design: a comparison between stress- and strain-based approaches. Rock Mech Rock Eng 51(2):599–611. https://doi.org/10.1007/s00603-017-1338-4

Lüdeling C, Spiers C, Urai J, Zapf D (2024) Creep of Rock Salt at Low Deviatoric Stress. Technical Report RR2023-3, SMRI

Munson DE, Dawson PR (1981) Salt constitutive modeling using mechanism maps. In: Reginald Hardy Jr., Langer, M. (eds.) 1^st^ Conference Mechanical Behavior of Salt (Saltmech I), University Park, PA, United States, pp 717–737

Minkley W, Knauth M, Fabig T, Farag N (2015) Stability and integrity of salt caverns under consideration of hydro-mechanical loading. In: Roberts, Mellegard, Hansen (eds.) Proceedings Saltmech8, Rapid City, SD, pp 217–227

Munson DE (1997) Constitutive model of creep in rock salt applied to underground room closure. Int J Rock Mech Min Sci 34(2):233–247. https://doi.org/10.1016/S0148-9062(96)00047-2
Nguyen Minh D, Brahams S, Durup JG (1993) Surface subsidence over deep solution mined storage cavern field. In: 3rd international conference on case histories in geotechnical engineering. University of Missouri, St. Louis, Missouri

Noll W (1974) The foundations of mechanics and thermodynamics - selected papers. Springer, Heidelberg

Reedlunn B, Argüello JG, Hansen FD (2022) A reinvestigation into Munson's model for room closure in bedded rock salt. Int J Rock Mech Min Sci 151:105007. https://doi.org/10.1016/j.ijrmms.2021.105007

Rouabhi A, Hévin G, Soubeyran A, Labaune P, Louvet F (2017) A multiphase multicomponent modeling approach of underground salt cavern storage. Geomech Energy Environ 12:21–35. https://doi.org/10.1016/j.gete.2017.08.002

Rouabhi A, Labaune P, Tijani M, Gatelier N, Hévin G (2019) Phenomenological behavior of rock salt: on the influence of laboratory conditions on the dilatancy onset. J Rock Mech Geotech Eng 11(4):723–738. https://doi.org/10.1016/j.jrmge.2018.12.011

Rouabhi A (2022) Phimef users' guide (in French). Mines Paris - Centre de Géosciences

Rouabhi A, Tijani M (2014) Problèmes Thermodynamiques liés aux stockage de fluides dans des cavités souterraines. Notice de conception de Demether (in French). Mines Paris - Centre de Géosciences

Speirs DCD, Bere A, Roberts D (2022) Geomechanical modelling of salt caverns under operational loading from hydrogen storage. In: ARMA/DGS/SEG international geomechanics symposium, pp 2022–036. https://doi.org/10.56952/IGS-2022-036

Sobolik SR, Ross TS (2022) Geomechanical analysis of oil storage caverns in salt domes with a low stress creep mechanism added to the M-D model. In: de Bresser, Drury, Fokker, Gazzani, Hanx, Niemeijer, Spiers (eds.) 10th Conference mechanical behavior of salt (Saltmech X), Utrecht, The Netherlands, pp 662–673

Spiers CJ, Schutjens PMTM, Brzesowsky RH, Peach CJ, Liezenberg JL, Zwart HJ (1990) Experimental determination of constitutive

parameters governing creep of rocksalt by pressure solution. Geol Soc Lon Spec Publ 54(1):215–227. https://doi.org/10.1144/GSL.SP.1990.054.01.21

Stormont JC (1997) In situ gas permeability measurements to delineate damage in rock salt. Int J Rock Mech Min Sci 34(7):1055–1064. https://doi.org/10.1016/S1365-1609(97)90199-4

Tijani M, Vouille G, Hugout B (1983) Le sel gemme en tant que liquide visqueux. In: Congrès International de Mécanique des Roches, Melbourne, Australia, pp 241–246

Urai JL, Spiers CJ, Zwart H, Lister GS (1986) Weakening of rock salt by water during long-term creep. Nature 324:554–557. https://doi.org/10.1038/324554a0

van Gent H, Hoedeman G, Bakker E (2024) The abandonment of large caverns and underground hydrogen storage in salt, views and lessons from a mining authority. In: SMRI Spring Technical Conference, Krakow, Poland

van Oosterhout B, Hangx S, Spiers CJ (2022) A threshold stress for pressure solution creep in rock salt: Model predictions vs. observations. In: de Bresser JH, Drury MR, Fokker PA, Gazzani M, Hangx SJ (eds.) 10th Conference Mechanical Behavior of Salt (Saltmech X), Utrecht, The Netherlands, pp. 57–67

Van Sambeek LL, Ratigan JL, Hansen FD (1993) Dilatancy of rock salt in laboratory tests. Int J Rock Mech Min Sci Geomech Abstr 30(7):735–738. https://doi.org/10.1016/0148-9062(93)90015-6

Wolters R, Lux K-H, Düsterloh U (2012) Evaluation of rock salt barriers with respect to tightness: influence of thermomechanical damage, fluid infiltration and sealing/healing. In: 7th conference mechanical behavior of salt (Saltmech VII), Paris, France, pp 417–426

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.