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

HAL
open science

# Characteristic features of salt-cavern behavior

Elizaveta Gordeliy, Pierre Berest

► To cite this version:

Elizaveta Gordeliy, Pierre Berest. Characteristic features of salt-cavern behavior. 2022. hal-03882692

HAL Id: hal-03882692

https://hal.science/hal-03882692

Preprint submitted on 2 Dec 2022

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.

Public Domain
# Characteristic features of salt-cavern behavior

*Elizaveta Gordeliy¹*, Pierre Bérest¹

¹Laboratoire de Mécanique des Solides, Ecole Polytechnique, 91128 Palaiseau, France

*elizaveta.gordeliy@polytechnique.edu

## Highlights

* A good fit against both short- and long-term volume evolutions as observed during in-situ tests is provided by the Munson-Dawson constitutive law.

* Following a pressure change in a salt cavern, a long transient evolution is predicted even when the constitutive law predicts no transient evolution during a uniaxial test.

* Transient evolutions of a salt cavern are due both to the transient behavior incorporated in the constitutive law and to the slow redistribution of stresses in the rock mass.

* Reverse creep and onset of tensile effective stresses are observed following an abrupt pressure increase.

* Volume loss rates are much larger when cavern pressure is cycled (rather than constant), especially when the Munson-Dawson constitutive law is selected.

* The Marketos-Spiers constitutive law, which includes pressure solution creep in addition to dislocation creep, predicts very fast cavern-volume loss rates when the cavern is shallow.

* When cavern growth during solution mining is taken into account, volume loss rate and stress distribution are significantly different from what they are when solution mining is described through a linearly decreasing pressure in a cavern with constant volume.

## Key Words

Salt caverns; Closed-form solution; Reverse creep; Pressure-cycled salt cavern; Pressure-solution creep; Moving boundary.
**Abstract**

The behavior of an idealized spherical cavern leached out from a salt formation can be described by a first-order integro-differential equation when the Poisson's ratio is 0.5. Three constitutive laws are considered: Norton-Hoff (N-H), Munson-Dawson (M-D) and Marketos-Spiers (M-S). The parameters of the M-D law can be fitted against both short- and long-term field data. It is shown that, following a pressure change, a "geometrical" transient evolution of the cavern volume is observed even when the constitutive law includes no "rheological" transient behavior (N-H, M-S). Following a large pressure increase, geometrical reverse creep and onset of tensile effective stresses at the cavern wall are observed. In a cavern subject to cyclic pressure, volume loss rate is faster than in a cavern where the average cycle pressure is applied. The loss rate is not extremely sensitive to cycle period; it is much faster when the M-D law (rather than the N-H law) is adopted. The M-S law, which accounts for the effect of pressure solution creep at low deviatoric stresses (in addition to dislocation creep), predicts a higher volume-loss rate (than the N-H law) when caverns are more shallow. A simple solution of the leaching problem (in which cavern radius is a function of time) can be found.

**Introduction**

The objective of this paper is to describe a simple numerical tool used to predict the behavior of hydrocarbon storage caverns leached out from salt formations. These caverns are known to shrink gradually when cavern pressure ($P_c$) is lower than geostatic pressure at cavern depth ($P_\infty$), as salt can be characterized as a Non-Newtonian fluid. Various constitutive laws have been suggested in the literature. For this problem, a couple of closed-form solutions are available — at least in the case of a spherical or cylindrical cavern when **steady state** has been reached [1-4]. A solution of the **transient** problem was obtained by Manivannan and Bérest [5] for an idealized cylindrical cavern in a Norton-Hoff medium. In the following, an idealized spherical shape is assumed. The behavior of the salt formation, both steady state and transient, is described by the Norton-Hoff or Munson-Dawson [6] laws. (The Marketos-Spiers law [7] is also considered.) The Poisson's ratio equals $\bar{\nu} = 0.5$ (i.e., salt is incompressible). The main drawback of such a tool is that the effects of cavern shape cannot be explored; they are known to be significant when the cavern is somewhat flat, for
instance. Its main advantage is that it provides insight in some special features of the
behavior of salt caverns which result from the high non-linearity of the strain-rate vs
stress relation describing salt behavior.

Nomenclature

<figure><img src="image_2.png" /><figcaption>Figure 1. Hollow sphere, main captions</figcaption></figure>

Nomenclature

$a_0$: internal radius of the hollow sphere (Fig. 1)

$A^*$: parameter of the Norton-Hoff power law

$b_0$: external radius of the hollow sphere (Fig. 1)

$c$: parameter of the Munson-Dawson law

$C$: a constant

$D$: salt grain size

$d_{ij}$: strain rate tensor

$E$: Young's modulus

$H$: cavern depth

$J_2$: second invariant of the stress tensor

$K_0$: parameter of the Munson-Dawson law

$P_c(t)$: internal cavern pressure at $r = a$
$P_\infty$: external cavern pressure at $r = b$

$P_c^{-}$: cavern pressure before an abrupt pressure change

$P_c^{+}$: cavern pressure after an abrupt pressure change

$P_c^{av}$: average cavern pressure during pressure cycles

$Q/R$: parameter of the Norton-Hoff power law

$Q_p / R$: parameter of the Marketos-Spiers law

$r, \varphi, \theta$: spherical coordinates (Fig. 1)

$S_{rr}, S_{\varphi\varphi}, S_{\theta\theta}$: radial, hoop and orthoradial deviatoric stresses

$S = \sigma_{rr} - \sigma_{\theta\theta}$: “deviator”

$S^{ss}(r)$: steady-state value of $S(r,t)$

$t$: time

$T_R$: absolute temperature of rock

$v$: radial displacement rate

$\Delta V$: Change in internal volume of the sphere

$V$: volume of the sphere

$\hat{x}(t)$: internal radius of the cavern during leaching

$\alpha_R$: thermal expansion coefficient of the rock mass

$\alpha_w$: parameter of the Munson-Dawson law

$\beta_w$: parameter of the Munson-Dawson law

$\gamma$: ratio ($a/b$) between internal and external radii of the sphere

$\epsilon$: strain

$\epsilon_t$: rheological transient strain

$\epsilon_{el}$: elastic strain

$\epsilon_s$: steady-state strain

$\zeta$: internal parameter of the Munson-Dawson law

$\mu$: shear modulus

$\bar{\nu}$: Poisson's ratio ($\bar{\nu}$ =0.5)

$\sigma_{rr}, \sigma_{\varphi\varphi}, \sigma_{\theta\theta}$: radial, hoop and orthoradial Cauchy stresses (tension-positive)
# 1. Main equations

## 1.1. Constitutive models

The Munson-Dawson (M-D) [6, 8], Norton-Hoff (N-H) and Marketos-Spiers (M-S) [7] models can be defined by the following set of equations:

$$d_{ij} = \frac{1+\bar{v}}{E} \dot{\sigma}_{ij} - \frac{\bar{v}}{E} \dot{\sigma}_{kk} \delta_{ij} + \frac{3}{2} (B(T_R) + A^*(T_R) \sqrt{3J_2})^{n-1} F(\sqrt{3J_2}, \zeta) s_{ij}$$

$$F(\sqrt{3J_2}, \zeta) = \exp\left[\Delta(\sqrt{3J_2})\left(1 - \frac{\zeta}{\varepsilon_i^*(\sqrt{3J_2})}\right)^2\right] \quad \text{when} \quad \zeta(r,t) < \varepsilon_i^*(\sqrt{3J_2})$$

$$F(\sqrt{3J_2}, \zeta) = 1 \qquad \text{when} \quad \zeta(r,t) = \varepsilon_i^*(\sqrt{3J_2}) \qquad (1)$$

$$F(\sqrt{3J_2}, \zeta) = \exp\left[-\delta(\sqrt{3J_2})\left(1 - \frac{\zeta}{\varepsilon_i^*(\sqrt{3J_2})}\right)^2\right] \quad \text{when} \quad \zeta(r,t) > \varepsilon_i^*(\sqrt{3J_2})$$

$$\dot{\zeta}(r,t) = A^*(T_R)(F(\sqrt{3J_2}, \zeta) - 1)\sqrt{3J_2(r,t)}^n$$

$$\varepsilon_i^*(\sqrt{3J_2}) = K_0 \exp(cT_R) \sqrt{3J_2}^m, \text{ and}$$

$$\Delta(r,t) = \alpha_w + \beta_w \log_{10}(\sqrt{3J_2} / \mu),$$

where $E$ and $\bar{v}$ are the Young's modulus and the Poisson's ratio, respectively; $T_R$ is the rock temperature; $\sigma_{ij}$ is the Cauchy stress tensor (in tension-positive sign convention), $d_{ij} = (v_{i,j} + v_{j,i})/2$ is the strain-rate tensor; $v_i$ is the rate of a material point; $\zeta$ is an internal parameter of the Munson-Dawson law; $s_{ij} = \sigma_{ij} - \sigma_{kk} \delta_{ij} / 3$ is the deviatoric stress tensor, and $J_2 = s_{ij} s_{ji} / 2$ is its second invariant.

* The N-H model, which includes no description of transient creep, is obtained when $B = 0$ and $F \equiv 1$ ($\Delta = 0$ and $\delta = 0$).

* The M-D model is obtained when $B = 0$.
* The M-S model, which includes a description of pressure-solution creep and no description of transient creep, is obtained when $B \neq 0$ and $F \equiv 1$.

<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic representation of the Munson-Dawson and Norton-Hoff models</figcaption></figure>

Schematically, in the case of a uniaxial test performed on a cylindrical sample (Fig. 2): a load increase ($\sigma_t$) at $t=0$ triggers both steady-state and transient creep (no transient creep in the case of the N-H model). Transient creep or $\varepsilon_t(t)$ is fast immediately after the load is applied, but exhausts itself after some time, until steady state is reached and $\varepsilon_t(t) = \varepsilon_t^*$. A load decrease (to $\sigma_2$) at $t_1 \to 2$ leads to a long period during which the total (steady state + transient) strain rate is small ("hesitation") before increasing slowly to the steady state value associated to the new stress level. A large load decrease triggers rheological reverse creep. (Sample height increases for a while, though the applied stress is compressive; reverse creep is not included in N-H or M-D models.)

## 1.2. The case of a hollow sphere

Taking into account spherical symmetry, $\underline{\sigma} = \sigma_{rr}(r,t)\underline{e}_r \otimes \underline{e}_r + \sigma_{\theta\theta}(r,t)\underline{e}_\theta \otimes \underline{e}_\theta + \sigma_{\phi\phi}(r,t)\underline{e}_\phi \otimes \underline{e}_\phi$,

$$\sigma_{\theta\theta} = \sigma_{\phi\phi},$$

$$s_{rr} = 2(\sigma_{rr} - \sigma_{\phi\phi})/3,$$

$$s_{\phi\phi} = s_{\theta\theta} = (\sigma_{\phi\phi} - \sigma_{rr})/3,$$

and

$$3J_2 = 3s_{ij}s_{ji}/2 = (\sigma_{rr} - \sigma_{\phi\phi})^2,$$
$\sqrt{3J_2} = |\sigma_{rr} - \sigma_{\phi\phi}|$. It is convenient to set $S(r,t) = \sigma_{rr}(r,t) - \sigma_{\phi\phi}(r,t) = -(r\partial\sigma_{rr}(r,t)/\partial r)/2$ (the "deviator"), or $|S| = \sqrt{3J_2}$ and $3s_{rr}/(2\sqrt{3J_2}) = S/|S|$. In addition, the rate of a material point $M$ can be written $\underline{v}(r,t) = v(r,t) \epsilon_r$ ($v$ is negative in most cases, Fig. 1).

It is convenient to set $I = (B(T_R) + A^*(T_R)|S|^{n-1}) F(|S|, \zeta)S$. The constitutive law (Eq. 2 and Eq. 3), the equilibrium condition (Eq. 4), the boundary conditions and the initial condition (Eq. 5), (Eq. 6), (Eq. 7) can be written as follows:

$$ \frac{\partial v}{\partial r} = \frac{\dot{\sigma}_{rr} - 2\bar{v}\dot{\sigma}_{\theta\theta}}{E} + I \quad (2) $$

$$ \frac{v}{r} = \frac{(1 - \bar{v})\dot{\sigma}_{\theta\theta} - \bar{v}\dot{\sigma}_{rr}}{E} - \frac{I}{2} \quad (3) $$

$$ \frac{r}{2} \frac{\partial \sigma_{rr}}{\partial r} + \sigma_{rr} - \sigma_{\theta\theta} = 0 \quad (4) $$

$$ \sigma_{rr}(a_0, t) = -P_c(t) \quad (5) $$

$$ \sigma_{rr}(b_0, t) = -P_{\infty} \quad (6) $$

$$ \sigma_{rr}(r, -\infty) = \sigma_{\phi\phi}(r, -\infty) = \sigma_{\theta\theta}(r, -\infty) = -P_{\infty} \quad (7) $$

$P_{\infty}$ is the geostatic pressure at cavern depth. Strains and displacements are small. An isothermal problem is considered (except in Section 3.2.2). When $\bar{v} = 0.5$ is assumed (elastic incompressibility), it can be inferred from Eq. (2) and Eq. (3) that $v(r,t) = a_0^2 \dot{a}(t) / r^2$, and Eq. (3) can be rewritten:

$$ \frac{a_0^2 \dot{a}(t)}{r^3} = \frac{r}{4E} \frac{\partial^2 \sigma_{rr}(r,t)}{\partial r \partial t} - \frac{1}{2} I(r,t) \quad (8) $$

Equation (8) can be divided by $r$ and integrated with respect to $r$ from $r=a_0$ to $r=b_0$. Taking into account $\partial\sigma_{rr}(a_0,t)/\partial t = -\dot{P}_c$ and $\partial\sigma_{rr}(b_0,t)/\partial t = 0$,

$$ \frac{1}{3} \frac{\dot{a}(t)}{a_0} (1 - \gamma^3) = \frac{\dot{P}_c(t)}{4E} - \frac{1}{2} \int_{a_0}^{b_0} I(v,t) \frac{dv}{v} \quad (9) $$

Also, when combining Eq. (8) and Eq. (9),
$$ \frac{\partial S(r,t)}{\partial t} = -EI(r,t) - \frac{a_0^3}{r^3} \frac{3}{1-\gamma^3} \left( \frac{1}{2} \dot{P}_c(t) - \int_{a_0}^{b_0} EI(v,t) \frac{dv}{v} \right) \quad (10) $$

and

$$ \dot{\zeta}(r,t) = A^*(F(r,t)-1)|S(r,t)|^n \quad (11) $$

In addition to these equations, initial conditions and boundary conditions must be stipulated. The boundary condition results from the definition of $S(r,t) = -(r\partial\sigma_{rr}(r,t)/\partial r)/2$:

$$ \int_{a_0}^{b_0} \frac{2S(v,t)}{v} dv = P_{\infty} - P_c \quad (12) $$

(The solution is a function of the difference $P_{\infty} - P_c$ only.) However, Eq. (12) implies that $\int_a^b 2\partial S(v,t)/\partial t \, dv/v = -\dot{P}_c$, and, from Eq. (10) it is sufficient that Eq. (12) be satisfied at $t=0$, for instance. In addition, initial conditions must be selected.

An important output of the computations is the volume change rate:

$$ \frac{\dot{V}(t)}{V_0}(1-\gamma^3) = \frac{9\dot{P}_c(t)}{4E} - \frac{9}{2}\int_{a_0}^{b_0} I(v,t) \frac{dv}{v} \quad (13) $$

where $V_0 = 4\pi a_0^3 / 3$ is the volume of the spherical cavern.

A numerical solution of Eq. (10) and Eq. (11) is obtained as follows. The numerical discretization is obtained by choosing a mesh of spatial nodes distributed radially, starting at the internal wall of the cavern. The obtained system of non-linear ODEs for $S(r,t)$ (coupled with non-linear ODEs for $\zeta(r,t)$ in case of the M-D law) is solved using the 4th order Runge-Kutta method. For the examples in this paper, 3000 to 20,000 spatial nodes were used.

### 1.3. Steady state for Norton-Hoff and Munson-Dawson laws

Note that steady state implies $\dot{P}(t)=0$ and $F=1$; the steady state solutions for the M-D and N-H constitutive laws are identical. When steady state is reached, $\partial S / \partial t = 0$; Eq. (10) implies that $S^{ss}(r) = C(a_0 / r)^{3/n}$ where $C$ is a constant which can be computed using Eq. (12):

$$ S^{ss}(r) = \left(\frac{a_0}{r}\right)^3 \frac{3(P_{\infty} - P_c)}{2n(1 - \gamma^{3/n})} \quad \frac{\dot{V}^{ss}}{V} = -\frac{3}{2} A^*(T_R) \left(\frac{3(P_{\infty} - P_c)}{2n(1 - \gamma^{3/n})}\right)^n \quad (14) $$
The steady-state solution for the M-S constitutive law is discussed in Section 5.

## 1.4. Norton-Hoff law, the case of linear viscoelasticity

This case is reached when $B = 0$, $F=1$ (N-H) and $n = 1$. Taking into account Eq. (12), Eq. (10) writes:

$$ \frac{\partial S(r,t)}{\partial t} = -EA^*S(r,t) - \frac{a_0^3}{r^3} \frac{3}{1-\gamma^3} \frac{1}{2} (\dot{P}_c(t) - EA^*(P_\infty - P_c)) \quad (15) $$

the solution for which, when $S(r,0^-) = 0$, is $S(r,t) = 3(a_0/r)^3(P_\infty - P_c(t))/(2(1-\gamma^3))$, which is also an elastic stress distribution. (When $P_c(0^-) = P_\infty$, $P_c(0^+) \neq P_c(0^-)$, and $P_c(t \ge 0^+)$ is fixed, stress distribution is not a function of time, in sharp contrast with the non-linear case; volume loss rate is constant.)

## 2. Fitting Munson-Dawson (M-D) law parameters against tests performed in an actual cavern

<figure><img src="image_4.png" /><figcaption>Figure 3. Daily volumes of liquid expelled from the cavern in Hugout's in-situ test (left), and volume loss obtained by modelling this test using the M-D law (right).</figcaption></figure>

In this section, it is shown that the Munson-Dawson model is able to capture both the transient and long-term features of actual cavern behavior. We use the data from two outflow tests performed by Hugout [9], in the 950-m-deep EZ53 cavern, where geostatic pressure is $P_\infty = 20.9$ MPa, in order to select the parameters for the M-D law. (This cavern currently is considered for hosting a hydrogen pilot plant (HYPSTER project, supported by the European Community).) An outflow test consists of recording the brine expelled from the cavern when the wellhead is opened. This outflow rate is deemed to reflect the effect of cavern creep closure. The effects on brine outflow of
additional dissolution associated with pressure changes and cavern-brine warming are not discussed here (for a complete discussion, see [10]). The first *in situ* test includes a pressure drop from $P_c^-=11.4$ MPa to $P_c^+=8$ MPa on day 93 after the end of leaching. Figure 3 (left) shows the daily volumes of liquid expelled from the cavern during Hugout's test; the transient-expelled flow was as high as 4500 liters/day before decreasing rapidly. The initial pressure, $P_c^+=11.4$ MPa, was restored on day 253 and kept constant for 13 years, after which the second *in situ* test was performed (Fig. 4). The expelled flow then was measured between 7 and 8 liters/day by Brouard, [11].

<figure><img src="image_5.png" /><figcaption>Figure 4. History of EZ53 outflow tests</figcaption></figure>

The M-D law was used to simulate the volume change due to the pressure schedule following Hugout's test and the volume-loss rate observed 13 years later. The leaching phase was simulated by linearly decreasing the cavern pressure. The elastic moduli were set to $E = 17,700$ MPa and $\mu = 5,900$ MPa, in order to fit cavern compressibility, as measured by Brouard [11]. For the parameters of the M-D law, the temperature was set to $T_R = 318$ K, and the power-law exponent was set to $n = 3.1$; the rest of the involved parameters were selected according to Table 1.
<table><caption>Table 1. Parameters values used in the numerical examples</caption><tr><td>Elastic</td><td>E = 17,000 MPa, μ = E/2/(1+ν) = 5667 MPa (in Sections 3 - 6)<br>or<br>E = 17,700 MPa, μ = E/2/(1+ν) = 5900 MPa (in Section 2)</td></tr><tr><td>Norton-Hoff</td><td>A* = A exp[-Q / (RT<sub>R</sub>)] , A = 0.64 /MPa<sup>n</sup>-yr, Q / R = 4100 K,<br>T<sub>R</sub> = 318 K, or T<sub>R</sub> = 283.15 K + (35 K / 1000 m) H<br>n = 3.1, or n = 3 (in Sections 3.3.2, 4 and 5)</td></tr><tr><td>Munson-Dawson</td><td>m = 3, K<sub>0</sub> = 7 × 10<sup>-7</sup> /MPa<sup>3</sup>, c = 0.00902 /K, α<sub>w</sub> = -13.2, β<sub>w</sub> = -7.738, δ = 0.58</td></tr><tr><td>Marketos-Spiers</td><td>B = A<sub>p</sub> exp(-Q<sub>p</sub> / (RT<sub>R</sub>)) , A<sub>p</sub> = 3.76 × 10<sup>-13</sup> / (T<sub>R</sub>D<sup>3</sup>) (MPa<sup>-1</sup>s<sup>-1</sup>) , Q<sub>p</sub> / R =<br>2946.8 K, T<sub>R</sub> = 283.15 K + (35 K / 1000 m) H</td></tr></table>

The parameters for the Norton-Hoff law were taken following Brouard [11]. Parameters of the M-D law were taken following Munson [8]. The M-D model is able to capture the sudden significant increase in the volume loss rate in the few days after the pressure drop (Fig. 3, right). For the second outflow test, the modelling result for the volume loss rate 13 years after the end of leaching is 7.8 liters/day (in comparison to the measured rate between 7 and 8 liters/day). The M-D model provides good agreement with the measurements for both outflow tests, even when the cavern shape is idealized. (Note also that brine warming - which is responsible for a 100 liters/day outflow rate on day 93, a 50 liters/day outflow rate on day 253 and a 0 liter/day during year 13 - is not taken into account in the computation.) The reason for this is that the M-D model is a kind of generalized Kelvin-Maxwell model, which is able to capture both transient and steady-state creep; simpler models (for instance, Lemaitre's model, which describes transient creep only) are less effective in this respect. Note that the M-D model, however, is not able to capture accurately the cavern reverse creep associated with the sudden pressure increase on day 254. (Rheological reverse creep should be incorporated in the model in order to capture it.)
### 3. Effects of an abrupt pressure change

<figure><img src="image_6.png" /><figcaption>Figure 5. Three schematic pressure changes in a cavern</figcaption></figure>

Three examples are discussed (Fig. 5).

#### 3.1. Effect of a cavern pressure drop from geostatic pressure, N-H law

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

<figure><img src="image_8.png" /><figcaption>Figure 6. Evolutions of the transient volume loss rate and the deviator (N-H law) following a pressure drop from geostatic ($P_{\infty}$ = 20.9 MPa) to halmostatic ($P_h$ = 11.4 MPa)</figcaption></figure>

In this example, only the N-H law ($B$ = 0 and $F \equiv 1$) is considered; no transient component is included in this law. Before $t = 0$, cavern pressure was equal to geostatic pressure, $P_c(t < 0) = P_{\infty}$, and $S(r, t < 0) = 0$. At $t = 0$, cavern pressure is decreased abruptly to halmostatic pressure, $P_c(t > 0) = P_h = 11.4$ MPa (Fig. 5, left). The deviator at $t = 0^+$ jumps to (Fig. 6, right):

$$S(r, 0^+) = \frac{a^3}{r^3} \frac{3}{1-\gamma^3} \frac{P_{\infty} - P_c}{2} \quad (16)$$
which is the elastic response to a pressure change. From Eq. (10), the derivative of the deviator is:

$$ \frac{\partial S}{\partial t}(r, 0^+) = EA^* \left( \frac{3}{2} \frac{P_\infty - P_c}{(1-\gamma^3)} \right)^n \left( \frac{1-\gamma^{3n}}{(1-\gamma^3)n} \frac{a_0^3}{r^3} - \frac{a_0^{3n}}{r^{3n}} \right) \quad (17) $$

This derivative is not zero when $n > 1$; after $t = 0$, cavern evolution is transient, even though the constitutive equation includes no transient component. Such a transient cavern evolution is said to be *geometrical*, as opposed to the rheological transient evolution, which results from the constitutive law. The geometrical transient creep is due to the slow redistribution of stresses in the rock mass from its initial elastic distribution to its final steady-state distribution (Fig. 6, right), an effect that does not exist when an idealized triaxial creep test is performed on a rock sample. (Stresses are uniform in the sample.) Note that the deviator decreases at the cavern wall ($\partial S(a_0, 0^+) / \partial t < 0$) and increases at a large distance from the cavern wall ($\partial S(r, 0^+) / \partial t > 0$ when $r / a_0 > [(1-\gamma^3)n / (1-\gamma^{3n})]^{1/(3n-3)}$). It slowly converges to the steady-state distribution, Eq. (14). A dilatant criterion, $\sqrt{3J_2} < \bar{a}[I_1]$, or $|S| < \bar{a}|2S + 3P_c|$ (DeVries et al. [12], Van Sambeek et al. [13]) often is used to assess cavern integrity; $I_1$ is the first invariant of the Cauchy stress tensor; $\bar{a}$ is a constant. As $|S| / |2S + 3P_c|$ is an increasing function of $S$, and $S(a_0, t)$ is a decreasing function of time at the cavern wall, this criterion is more demanding immediately after cavern creation rather than later.

## 3.2. Effect of a pressure decrease

### 3.2.1. Comparison between the N-H and M-D laws

In the second example, geostatic pressure $P_\infty = 20.9$ MPa is applied at $r = b_0$. A constant inner pressure, $P_0 = 10$ MPa, was maintained at the cavern wall during a very long period of time before $t = 0$. At $t = 0^-$, steady state (Eq. 14) was reached.
<figure><img src="image_9.png" /><figcaption>Figure 7. Transient volume rate after abrupt pressure decrease from $P_c^-$ = 10 MPa to $P_c^+$ = 5 MPa.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Behavior of $S(r,t)$ obtained with N-H law or M-D law after an abrupt pressure decrease from $P_c^-$ = 10 MPa to $P_c^+$ = 5 MPa.</figcaption></figure>

At $t = 0^+$, the internal pressure is decreased abruptly from $P_c^-$ = 10 MPa to $P_c^+$ = 5 MPa (Fig. 5, middle), and the deviator experiences an elastic jump by

$$[S]_-^+ = -\frac{a^3}{r^3} \frac{3}{1-\gamma^3} \frac{1}{2}(P_c^+ - P_c^-) \quad (18)$$

After $t = 0^+$, Eqs. (10) and (11) are used. The parameters were taken from the simulations of the *in situ* tests performed on Cavern EZ53 (Section 2, Table 1), except that the elastic moduli were set in this example to $E = 17,000$ MPa and $\mu = 5667$ MPa.
The results were computed using 11,144 spatial nodes. Figure 7 shows the evolution of the relative cavern-volume rate, scaled by the initial cavern volume, $V_0$. The most striking fact is that, as in the first example, a long transient evolution (more than one century — much longer than the transient creep observed during a laboratory test on a salt sample) is observed in both models (M-D and N-H). The steady-state volume-loss rate before the abrupt pressure drop was $\dot{V}^{ss}/V_0 = -4.2 \times 10^{-4} \text{ yr}^{-1}$. The steady-state volume-loss rate at large time after the pressure drop, $\dot{V}^{ss}/V_0 = -1.35 \times 10^{-3} \text{ yr}^{-1}$, is not reached until a century after the pressure drop. Here again, the transient behavior in the Norton-Hoff law results from the slow redistribution of stresses in the rock mass, the so-called geometrical transient behavior. However, the volume rates obtained with the N-H law and with the M-D law become approximately the same about 5 years after the pressure drop — i.e., the largest part of the *rheological* transient creep in the simulation with the M-D law occurs soon after the pressure drop, while the transient behavior can be attributed to the geometrical creep for times much longer than about 5 years after the pressure drop.

Figure 8 shows the deviator, $S(r,t) = \sigma_{rr}(r,t) - \sigma_{\phi\phi}(r,t)$, for both models. At the cavern wall, the deviator stress decreases ($\partial S(r,t)/\partial t < 0$) to reach its steady-state value slowly. At the cavern wall, the deviator equals the opposite of the effective stress, $S(a_0,t) = -P_c - \sigma_{\phi\phi}(a_0,t)$ (see Section 3.3.3). When the cavern pressure is decreased suddenly from $P_c^- = 10 \text{ MPa}$ to $P_c^+ = 5 \text{ MPa}$ at $t=0^+$, the deviator at the cavern wall jumps from 5.3 MPa to 12.8 MPa, as predicted by Eq. (18). For $t > 0$, the effective stress changes much faster in the M-D model than in the N-H model. In the M-D model, the deviator at the cavern wall drops from 12.8 MPa to 9.4 MPa within a day after the sudden pressure decrease, and then to 8.6 MPa within 100 days. In the N-H model, there is almost no change in the deviator within a day after the sudden pressure decrease, and it drops slowly from 12.8 MPa to 9.5 MPa within 100 days. Due to the much faster transient-stress evolution after the pressure drop, the M-D law predicts a much higher volume-loss rate (by two orders of magnitude) in the first few days after the pressure drop, in comparison to the N-H law (Fig. 7).
### 3.2.2. Thermal stresses

Especially in the case of the N-H law, during several days after the initial abrupt pressure drop, rock mass response is almost perfectly elastic. Viscoplastic strains can be disregarded. Moreover, in the case of a gas-filled cavern, gas experiences a severe temperature drop during a rapid (a few days) depressurization, as gas heat capacity is low (much lower than brine's, for instance). Rock temperature change rate during the depressurization period can be written $\dot{\tau}(r,t) < 0$. Onset of thermal stresses must be taken into account and rock mass response during this period is linear thermo-elastic,

$$[S]_{-}^{+} = [S_{EL}^{+}]_{-} + [S_{TH}^{+}]_{-}$$

$[S_{EL}^{+}]_{-}$ is the response to the cavern pressure change when $\dot{\tau}(r,t) = 0$. For simplicity, $b_0 = \infty$ ($\gamma = 0$) is assumed:

$$[S_{EL}^{+}]_{-} = -\frac{3(P_{c}^{+} - P_{c}^{-}) a^{3}}{2 r^{3}} \quad (19)$$

$[S_{TH}^{+}]_{-}$ is the rock mass response to cavern temperature change when $\sigma_{rr}^{TH}(a_0) = 0$, $\sigma_{rr}^{TH}(b_0 = \infty) = 0$ and $r^3 v^{TH}(r) \to 0$ when $r \to \infty$. For this problem, Eq. (2) and Eq. (3) must be rewritten:

$$\frac{\partial v^{TH}}{\partial r} = \frac{\dot{S}^{TH}}{E} + \alpha_{R} \dot{\tau} \quad (20)$$

$$\frac{v^{TH}}{r} = -\frac{\dot{S}^{TH}}{2E} + \alpha_{R} \dot{\tau} \quad (21)$$

Where $\alpha_R$ is the thermal expansion coefficient of the rock mass. Eq. (21) can be subtracted from Eq. (20), the resulting equation can be divided by $r$ and integrated with respect to $r$ between $r=a_0$ and $r=b_0 = \infty$ leading to $[v^{TH}(r)/r + 3\dot{\sigma}_{rr}^{TH}(r)/(4E)]_{a_0}^{\infty} = 0$, from which $v^{TH}(a_0)/a_0 = 0$ (thermal stresses generate no cavern volume change). Integrating Eq. (21) with respect to time from the beginning to the end of the depressurization period leads to:

$$[S_{TH}^{+}(a_0)]_{-}^{+} = 2E\alpha_{R}[\tau(a_0)]_{-}^{+} \quad (22)$$

And
$$[S(a_0)]^+ = [S^{EL}(a_0)]^+ + [S^{TH}(a_0)]^+ = -\frac{3}{2}(P_c^+ - P_c^-) + 2E\alpha_R[\tau(a_0)]^+ \quad (23)$$

The first quantity is positive, the second quantity is negative. The overall deviator at cavern wall is negative (the effective stress is tensile and thermal fractures can be created) when thermal stresses are large enough, i.e., when the thermodynamic evolution is fast enough, as explained in Bérest et al. [14].

### 3.3. Effect of a cavern pressure increase: reverse creep, tensile effective stress

<figure><img src="image_11.png" /><figcaption>Figure 9. Transient volume rate after abrupt pressure increase from $P_c^- = 10$ MPa to $P_c^+ = 15$ MPa.</figcaption></figure>

#### 3.3.1. An example

In this section, the N-H law and the M-D law are used to simulate the case of an abrupt pressure increase (Fig. 5, right). The geostatic pressure is $P_\infty = 20.9$ MPa. The initial conditions and the parameters of the simulations are the same as described in Section 3.2.1, except that at $t = 0^+$ in the present case, the internal pressure is increased abruptly from $P_c^- = 10$ MPa to $P_c^+ = 15$ MPa.

Figure 9 shows the rate of the cavern volume, scaled by the initial cavern volume, $V_0$. The steady-state volume-loss rate before the abrupt pressure change was $\dot{V}^{ss}/V_0 = -4.2 \times 10^{-4}$ yr$^{-1}$. The steady-state volume-loss rate occurring a long time after the pressure increase is $\dot{V}^{ss}/V_0 = -6.3 \times 10^{-5}$ yr$^{-1}$. Here, again, a long transient evolution (more than one century) is observed in both models. After the cavern pressure is
increased from $P_c^-=10$ MPa to $P_c^+=15$ MPa, volume loss rate (which is negative) increases first (its absolute value becomes smaller) before finally decreasing. In the N-H model, volume loss rate starts to decrease after 10 years; this occurs after 30 years in the M-D model. However, even after 100 years, the steady state creep-closure rate, $\dot{V}^{ss}/V_0 = -6.3 \times 10^{-5}$ yr$^{-1}$, is not reached.

### 3.3.2. Geometrical reverse creep

<figure><img src="image_12.png" /><figcaption>Figure 10. Volume rate immediately after an abrupt pressure increase (or decrease) from $P_c^- = 10$ MPa to $P_c^+$. Reverse creep is observed when $P_c^+ > 16$ MPa. The results are obtained with the N-H law, $n = 3.1$.</figcaption></figure>

Note, also, in the example above, volume loss rate remained negative (cavern volume decreased) after the abrupt pressure increase. In fact, it can become positive when the pressure increase is sufficiently large. An example of this in a real cavern was described by Denzau and Rudolph [15]. Indeed, Fig. 10 shows the volume rate immediately after the abrupt pressure jump from $P_c^- = 10$ MPa to the values of $P_c^+$ belonging to the $P_c^+ = 5-20$ MPa domain. The results were computed using the N-H law, and similar results have been obtained using the M-D law. It is seen that for a large enough pressure increase (approximately, $P_c^+ > 16$ MPa in this example), the volume rate becomes positive right after the pressure increase. Immediately after an abrupt (and large) pressure increase, reverse transient volume loss is observed: the cavern volume increases for a while, even though the cavern pressure is lower than the
geostatic pressure. (This is remarkable, as neither M-D nor N-H laws are able to describe rheological reverse creep (see Fig. 2); this behavior can be attributed to geometrical reverse creep.)

In the case of the N-H law, this result can be proved in a more general way, as follows. (It is assumed that steady state was reached at $t = 0$, $P_c = P_c^-$; at $t = 0$, pressure is changed abruptly from $P_c = P_c^-$ to $P_c = P_c^+$, $P_c^- < P_c^+ < P_\infty$. An integer value of $n$ (the exponent of the power law) is selected, $n = 3$. Immediately after the pressure change, from Eq. (14) and Eq. (18):

$$S(r, 0^+) = \frac{a_0}{r} \frac{(P_\infty - P_c^-)}{2(1-\gamma)} - \frac{a_0^3}{r^3} \frac{3(P_c^+ - P_c^-)}{2(1-\gamma^3)} \quad (24)$$

When setting $X = (P_c^- - P_c^+) / (P_\infty - P_c^-)$, integration of Eq. (13) leads to the following result (in which $\gamma = 0$ is set for simplicity),

$$\frac{\dot{V}(0^+)}{V_0} = -\frac{9}{16} A^* (P_\infty - P_c^-)^3 \left[ \frac{1}{3} + \frac{9X}{5} + \frac{27X^2}{7} + \frac{3X^3}{1} \right] \quad (25)$$

where $0 > X > -1$. Note that at $t = 0^-$, $P_c^+ = P_c^-$, $X = 0$, and $\dot{V}(0^-) < 0$. When $X = -1$ (i.e., $P_c^+ = P_\infty$),

$$\frac{\dot{V}(0^+)}{V_0} = \frac{9}{2} A^* (P_\infty - P_c^-)^3 \frac{8}{3 \times 5 \times 7} > 0 \quad (26)$$

i.e., cavern volume increases when $P_c^+$ is close enough to (but smaller than) $P_\infty$.

### 3.3.3. Onset of tensile effective stresses
<figure><img src="image_13.png" /><figcaption>Figure 11. Behavior of S(r,t) obtained with the N-H law or the M-D law after abrupt cavern pressure increase from $P_c^- = 10$ MPa to $P_c^+ = 15$ MPa.</figcaption></figure>

At the cavern wall, the deviator equals the opposite of the effective stress. Immediately after the pressure change:

$$S(a_0, 0^+) = -P_c^+ - \sigma_{\phi\phi}(a_0, 0^+) = \frac{3}{2} \left[ \frac{P_\infty - P_c^-}{n(1 - \gamma^{3/n})} - \frac{P_c^+ - P_c^-}{(1 - \gamma^3)} \right] \quad (27)$$

This becomes negative when $(P_c^+ - P_c^-) > (1 - \gamma^3)(P_\infty - P_c^-) / [n(1 - \gamma^{3/n})]$, a condition that is met easily when $n$ is large enough (typically, $n = 3$ to 6). In the example shown on Fig.11, effective stress is tensile at cavern wall immediately after the pressure jump and remains tensile for a long period. (There is a risk of fracturing, an idea mentioned by Wallner [16] or Wawersik and Stone [17]). This result is similar to that for the cylindrical cavity (Fig. 7 in Bérest and Manivannan [18]).

### 3.3.4. The Pressure Observation Test

A Pressure Observation Test (also called Pressure Vessel Test, or PV Test) consists of increasing cavern pressure from halmostatic to the Maximum Operating Pressure through brine injection in the shut-in cavern; and monitoring further pressure evolution (Bérest et al. [19]). A fast pressure decrease is deemed to be a clear sign of a leak. In fact, evolution of pressure in a shut-in cavern is governed by creep closure, brine thermal expansion and leaks [10]. Brine mass conservation can be written $d(\rho_b V_b) = -\rho_b Q_{leak}$ where brine state equation is $\dot{\rho}_b = \rho_b(\beta_b \dot{P}_c - \alpha_b \dot{T}_c)$, $\beta_b = 2.7 \times 10^{-4}/MPa$ is the brine compressibility, $\alpha_b = 4.4 \times 10^{-4}/°C$ is the thermal expansion coefficient of brine, $Q_{leak}$
are the brine leaks through the cavern wall or the wellbore, and $V_b = V$ is the volume of brine in the cavern. In the following, only the effect of creep closure is considered and $\beta_b \dot{P}_c + \dot{V}/V_0 = 0$. Taking into account Eq. (13):

$$ \left(E\beta_b(1-\gamma^3)+\frac{9}{4}\right)\dot{P}_c(t)-\frac{9E}{2}\int_{a_0}^{b_0}I(v,t)\frac{dv}{v}=0 \quad (28) $$

$$ \frac{\partial S(r,t)}{\partial t} = -EI(r,t) + 3\frac{a_0^3}{r^3}\frac{1}{1-\gamma^3+\frac{9}{4\beta_b E}}\int_{a_0}^{b_0} EI(v,t)\frac{dv}{v} \quad (29) $$

<figure><img src="image_14.png" /><figcaption>Figure 12. Pressure evolution in a 950-m deep shut-in cavern after cavern pressure is increased from 11.4 MPa to 18 MPa (a Pressure Observation Test)</figcaption></figure>

In the M-S case, the independent non-dimensional variables in Eq. (29) are $r/a_0$, $t/(A^*E(P_\infty - P_0)^{n-1})$, $n$, $S/(P_\infty - P_0)$, $B/(A^*(P_\infty - P_0)^{n-1})$, $1-\gamma^3+9/(4\beta_b E)$. Note that in the case of the M-D law ($B=0$) or the N-H law ($B=0, F=1$), a change in $A^*$ is equivalent to a change in the unit of time. As $\beta_b$ is a given physical constant, and $\gamma$ is exceedingly small, the non-dimensional solution of Eq. (29) mainly depends on $n$, $P_\infty - P_0$, $E$, $B/A^*$, and, in the case of the M-D law, on the parameters involved in the definition of $F$.
An example is provided in Fig. 12. The M-D law is considered. Geostatic pressure is $P_{\infty} = 20.9$ MPa. Pressure is increased abruptly at $t = 0$ from 11.4 MPa to 18 MPa before plugging the well. During an 8-year long period, cavern pressure is below the initial pressure even when no actual leak is considered.

# 4. Volume loss during pressure cycles

## 4.1. An example

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Volume change for a cavern under cyclic pressure, $P_{\infty} - P_c^+(t) = 10$ MPa – $\Delta P \sin(2\pi t/T)$. Results are computed with the N-H law or the M-D law for several sets of cycle amplitude $\Delta P$ and period $T$.</figcaption></figure>

In this section, the M-D and N-H laws are used to simulate the case of a cyclic pressure change. The elastic moduli are set to $E = 17,000$ MPa and $\mu = 5667$ MPa. The
parameters of the M-D and N-H laws are taken from the simulations of the in situ tests performed on Cavern EZ53 (Section 2), except that the power law exponent *n* is set to *n* = 3 in these examples. A constant pressure difference, $P_{\infty} - P_{c}^{-} = 10$ MPa, was maintained at the cavern wall over a long period of time before *t* = 0, such that steady state was reached at *t* = 0. At *t* = 0+, the cavern pressure is cycled according to $P_{\infty} - P_{c}^{+}(t) = P_{\infty} - P_{c}^{av} - \Delta P \sin(2\pi t / T)$, where $P_{c}^{av} = P_{c}^{-}$ (average cavern pressure during the cycles), $\Delta P$ (amplitude of the cycles) and *T* (period of the cycles) are three constants. On Fig. 13 (upper left), the M-D law, short cycle periods (1 day and 1 week) and two cycle amplitudes ($\Delta P$ = 3 MPa and 6 MPa) are considered; the average pressure difference is $P_{\infty} - P_{c}^{av} = 10$ MPa. The steady-state volume evolution (when cavern pressure equals average pressure) is drawn for comparison. Volume loss is more than doubled when cycle amplitude is doubled, reflecting the non-linearity of the creep law. Influence of the period is more discreet. Figure 13 (upper right) illustrates the difference between the N-H and M-D (more realistic) creep laws. Volume loss is considerably larger when the M-D law is selected; this is expected, as the M-D law takes into account the rheological transient behavior of salt, which is activated during pressure cycles. Longer periods are considered in the lower pictures (a one-year period is typical of the historical operation mode of natural gas caverns). Lessons are similar to those drawn when the cycle period is shorter: volume loss is much larger when the M-D law is selected; the influence of the cycle period (1 month and 1 year) is relatively minor. (Note, however, that more time is left to restore cavern volume when cavern pressure is higher.)

## 4.2. A simplified model

These facts suggest a possible simplified model. The N-H law ($B = 0, F \equiv 1$) is selected. Deviator evolution can be described roughly as the sum of the elastic response to pressure cycles,

$$S^{el}(r,t) = -\frac{3a_0^3}{2(1-\gamma^3)r^3} \Delta P \sin(2\pi t / T) = -S_{\max}^{el} \sin(2\pi t / T) \quad (30)$$

and an average response $\bar{S}(r,t)$ which varies very slowly during a cycle,
$$S(r,t) = \bar{S}(r,t) - S_{\max}^{el} \sin(2\pi t / T) \quad (31)$$

where $T \partial \bar{S}(r,t) / \partial t \ll \bar{S}$. During one period, $\bar{S}(r,t) = \bar{S}(r)$ is almost constant. The equations become simpler when the exponent of the power law is an integer (for instance, $n=3$). In this case, averaging Eq. (10) for $\partial S / \partial t$ over one cycle leads to the following integro-differential equation for $\bar{S}(r,t)$ (this method was used in Lestringant et al. [20]):

$$\frac{\partial \bar{S}(r,t)}{\partial t} = -EA^* \left[ \bar{S}^3(r,t) + \frac{3}{2} \bar{S}(r,t) (S_{\max}^{el})^2 \right] + \frac{a_0^3}{r^3} \frac{3}{1-\gamma^3} EA^* \int_{a_0}^{b_0} \left( \bar{S}^3(v,t) + \frac{3}{2} \bar{S}(v,t) (S_{\max}^{el})^2 \right) \frac{dv}{v} \quad (32)$$

In addition, averaging Eq. (13) over one cycle leads to

$$\frac{\bar{V}(t)}{V_0}(1-\gamma^3) = -\frac{9}{2}A^* \int_{a_0}^{b_0} \left( \bar{S}^3(v,t) + \frac{3}{2}\bar{S}(v,t)(S_{\max}^{el})^2 \right) \frac{dv}{v} \quad (33)$$

From Eq. (32) and Eq. (33), it is seen that the average response depends on cycle amplitude — not on the period of the pressure cycle — in agreement with Fig. 13 (lower left). In the very long term, $\partial \bar{S}(r,t) / \partial t = 0$ and $\bar{S}^3 + 3(S_{\max}^{el})^2 \bar{S} / 2 = C(a_0 / r)^3$, where $C$ is a constant that can be determined numerically through Eq. (12), $\int_{a_0}^{b_0} 2\bar{S}(v,t)dv / v = P_{\infty} - P_c^{av}$. Finally, from Eq. (33), the volume-loss rate averaged over large time periods can be approximated by

$$\frac{\bar{V}}{V_0} \approx -\frac{3}{2} A^* C \quad (34)$$

Again, constant $C$ depends on cycle amplitude and the average cavern pressure; it does not depend on the period of the pressure cycles.

Figure 14 shows the results obtained using Eq. (34) — namely, the average volume-loss rate as a function of cycle amplitude. As discussed in [20], this large time approximation cannot be mathematically exact, but it can provide an estimate for the volume-loss rate at large periods of time. Indeed, Fig. 15 shows the evolution of the
relative volume (left) and relative average-volume rate (right) for a cavern subject to
the cyclic pressure with $P_{\infty} - P_c^{av} = 10 \text{ MPa}$, $\Delta P = 3$ or $6 \text{ MPa}$, and $T = 1$ year. The transient
results are computed with the N-H law. It is seen that the large time approximation,
obtained from Eq. (34), provides a reasonable estimate for the volume-loss rate after
about 20 years of pressure cycles (in this example).

<figure><img src="image_16.png" /><figcaption>Figure 14: Relative average volume-loss rate for a cavern under cyclic pressure, as a function of cycle amplitude. Results are computed with the N-H law by using the large-time approximation [Eq. (34)].</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15. Relative volume (left) and relative average volume rate (right) for a cavern under cyclic pressure as a function of time. Results are computed with the N-H law.</figcaption></figure>
5. The Marketos-Spiers law predicts much faster volume loss rate than the N-H law when cavern depth is shallower; steady state is reached faster.

It long has been suspected ([21], [22]) that, in the small deviatoric stress domain ($\sigma < 3$ MPa, typically), the governing mechanism for salt creep is pressure solution rather than dislocation creep. A consequence of this should be that the creep rate in this domain is much faster — by several orders of magnitude — than that extrapolated from tests performed in the high stress domain, a view confirmed by creep tests on natural salt performed in the range $\sigma < 4.5$ MPa ([23], [24]). Various attempts were made to take into account both pressure solution and dislocation displacements, which are the governing creep mechanisms in the small (respectively, large) deviatoric stress domains ([25], [7], [3], [4], [5], [18], [26]). The Marketos-Spiers (M-S) law, described above ($\dot{\epsilon} = B\sqrt{3J_2} + A^*\sqrt{3J_2}''$), is a simple such attempt. In this section, the exponent of the non-linear term in the M-S law and in the N-H law is set to $n = 3$. For this case, the steady-state solution for the M-S law is obtained semi-analytically. Constants $A^*$ and $B$ are functions of temperature; in addition, $B$ is proportional to $D^{-3}$, where $D$ is the grain diameter of salt ([21], [7]). Constants $A^*$ and $B$ are selected according to Table 1. The choice for $B$ corresponds to the mean values of the involved parameters ($A_p, Q_p$) given by Spiers et al. [21] and Marketos et al. [7]. Here, the authors use grain diameter $D = 7$ mm. For constant $A^*$, the values of the involved parameters ($A, Q$) are taken from the in-situ tests performed on Cavern EZ53 (Section 2). The values of both $A^*$ and $B$ in this example take into account the temperature dependence on depth, $T_R = 283.15$ (K) + 35 (K)/1000(m) × $H(m)$ (rock temperature is smaller when the cavern is shallower). It can be expected that the addition of the pressure solution term is especially significant when deviatoric stresses in the rock mass are small — i.e., when cavern pressure is high or in the case of a shallow cavern. However, in this last case, this effect may be counterbalanced by the effect of lower rock temperature on creep rate.

We consider first the steady-state solution as a function of the cavern depth $H$. Steady state is defined by $\partial S / \partial t = 0$ in Eq. (10), or

$$0 = -I^{ss} + \frac{a_0^3}{r^3} \frac{3}{1-\gamma^3} \left( \int_{a_0}^{b_0} I^{ss}(v) \frac{dv}{v} \right) \quad (35)$$

From this, since $|S| = S$ when $P_\infty > P_c$:

$$I^{ss}(r) = BS^{ss}(r) + A^*S^{ss}(r)^n = C/r^3 \quad (36)$$

This algebraic equation can be solved numerically, $S^{ss} = S^{ss}(r, C')$, and $C'$ is determined through the condition $\int_{a_0}^{b_0} 2S^{ss}(v)dv/v = P_\infty - P_c$, Eq. (12). (For $n = 3$, the solution $S^{ss} = S^{ss}(r, C')$ is found analytically from the cubic Eq. (36).)
<figure><img src="image_18.png" /><figcaption>Figure 16. Left figures: scaled stress deviator $S/(P_\infty - P_c)$, computed using the N-H law and the M-S law for cavern depths: $H = 200, 500, 1000, 1500, 2000$ m. Right figures: ratio of the steady-state volume rate for the M-S law to the steady-state volume rate for the N-H law, as a function of cavern depth $H$. The results, obtained with the M-S law, correspond to grain size $D = 7$ mm. The results were obtained assuming either fixed temperature 318 K (upper figures), or temperature variation with depth (lower figures).</figcaption></figure>

Geostatic pressure is $P_\infty = \gamma_R H$, $\gamma_R = 0.022$ MPa/m, and cavern pressure is $P_c = \gamma_b H$, $\gamma_b = 0.012$ MPa/m. Figure 16 (left) shows the stress deviator obtained for several cavern depths: $H = 200, 500, 1000, 1500, 2000$ m. To show the effect of temperature variation with depth onto the obtained values, the results are shown for either a constant temperature $T_R = 318$ K (upper figures), or for a temperature varying with depth $T_R = 283.15(\text{K}) + 35(\text{K})/1000(\text{m}) \times H(\text{m})$ (lower figures). The ratio of the volume rate for the M-S law to the rate for the N-H law is larger for shallower depths in both cases, see Fig. 16 (right). It is more than 30 for a very shallow depth ($H = 200$ m, and $T_R = 290.15$ K).
<figure><img src="image_19.png" /><figcaption>Figure 17. Behavior of S(r,t) (left) and transient volume rate (right) after an abrupt pressure drop. Cavern depth is H = 250 m. The black dashed lines (right) show the steady-state volume rate at large time.</figcaption></figure>

We then consider a shallow cavern, $H = 250$ m, in which the geostatic pressure is $P_\infty = \gamma_R H = 5.5$ MPa, and the temperature is $T_R = 291.9$ K. Before $t=0$, steady state was reached with cavern pressure $P_c^- = 3$ MPa. At $t=0$, the cavern pressure is dropped abruptly to $P_c^+ = 2.5$ MPa. Figure 17 (left) shows the behavior of the deviator near the cavern wall after the pressure drop. For the M-S law, the steady-state solution is reached much faster than for the N-H law (when $A^*\sqrt{3J_2}^n$ is much smaller than $B\sqrt{3J_2}$, cavern behavior is close to linear visco-elastic, and stress distribution remains almost constant, see Section 1.4). Indeed, the solution at $t=9$ years for the M-S law is closer to the corresponding steady-state solution than the solution at $t=100$ years for the N-H law. Figure 17 (right) shows that the volume-loss rate, obtained with the M-S law, is much higher than that for the N-H law — roughly by one order of magnitude. The volume-loss rate at 10 years after the pressure drop for the M-S law can be approximated well by the corresponding steady-state rate.
# 6. Solution mining (a moving boundary problem)

## 6.1. Introduction

<figure><img src="image_20.png" /><figcaption>Figure 18. Three descriptions of the leaching phase.</figcaption></figure>

Salt caverns are created through solution mining. A wellbore is tapped to the salt formation. Water is injected in the wellbore through a central string, brine is withdrawn from the annulus between the string and the casing, and the cavity grows until its final shape is reached. Brine pressure is applied to the cavern walls. The cavern wall is a moving boundary. There are different possible methods (see Fig.18) to describe the solution-mining phase: cavern pressure abruptly decreases from geostatic pressure, $P_\infty = \gamma_R H$, to halmostatic pressure, $P_c = \gamma_b H$, (method A); pressure decreases linearly (method B, often adopted when performing numerical computations); and the cavern is created progressively. The stress distribution obtained when using methods A or B cannot be exact, as there is a competition between stress redistribution due to creep closure and cavern growth due to solution mining, an effect that is not accounted for in these two methods.

The equations describing the moving boundary method are discussed below. Only small strains and displacements of material points are considered. All results in this section were obtained using the M-D law. The sphere's external radius, $b_0$, is fixed and submitted to a constant external pressure, $P_\infty$ (geostatic pressure). For the spherical
cavern, methods A and B are based on the assumption that the cavern is created instantaneously with its final radius, $a_0$. For the moving boundary method, evolution of the cavern radius, $\hat{x} = \hat{x}(t)$, results from a physical assumption. For instance, wellbore creation is described as follows: at $t = 0$, $P_c(t)$ experiences a jump from $P_\infty$ to $P_c$; $\hat{x}(0)$ is very small when the initial jump is performed. Later on, solution mining is performed in such a way that cavern growth rate is constant, or $\hat{x}^3 - \hat{x}_0^3 = 3Kt$, $0 < t < t_1$, where K is a constant, and $\hat{x}^2\hat{x} = K$. For $t > t_1$, $\hat{x}(t)$ is constant, $\hat{x}(t > t_1) = \hat{x}_1$.

Equations (2) to (7) still hold, except for Eq. (5). That equation must be replaced by:

$$ \sigma_{rr}(\hat{x}(t)) = -P_c(t) \qquad (37) $$

The Poisson's ratio is assumed to be $\bar{\nu} = 0.5$; from Eqs. (2) and (3), $v(r,t)/r = b_0^2\dot{b}(t)/r^3$ (where $\dot{b}(t) = v(b_0,t)$ is the displacement rate at the external wall of the hollow sphere). Using Eq. (4), equation (3) can be rewritten as

$$ \frac{b^2 \dot{b}(t)}{r^3} = \frac{r}{4E} \frac{\partial^2 \sigma_{rr}(r,t)}{\partial r \partial t} - \frac{I}{2} = -\frac{1}{2E} \frac{\partial S(r,t)}{\partial t} - \frac{I}{2} \qquad (38) $$

Equation (38) can be divided by $r$ and integrated with respect to $r$ from $r=\hat{x}$ to $r=b$. Taking into account $\partial\sigma_{rr}(b,t)/\partial t = 0$,

$$ \frac{d\sigma_{rr}(\hat{x},t)}{dt} = \frac{\partial\sigma_{rr}(\hat{x},t)}{\partial t} + \dot{\hat{x}}\frac{\partial\sigma_{rr}(\hat{x},t)}{\partial r} = \frac{\partial\sigma_{rr}(\hat{x},t)}{\partial t} - \frac{2\dot{\hat{x}}S(\hat{x},t)}{\hat{x}} = -\dot{P}_c \qquad (39) $$

$$ \frac{1}{3} \frac{\dot{b}(t)}{b_0} \left(-1 + \frac{b_0^3}{\hat{x}^3}\right) = \frac{\dot{P}_c(t) - 2S(\hat{x},t)\dot{\hat{x}}/\hat{x}}{4E} - \frac{1}{2} \int_{\hat{x}}^{b_0} I(v,t) \frac{dv}{v} \qquad (40) $$

When combining Eq. (38) and Eq. (40):

$$ \left( \frac{\partial S(r,t)}{\partial t} + EI(r,t) \right) \left( -\frac{1}{b_0^3} + \frac{1}{\hat{x}^3} \right) = -\frac{3}{2} \frac{\dot{P}_c(t) - 2S(\hat{x},t)\dot{\hat{x}}/\hat{x}}{r^3} + \frac{3}{r^3} E \int_{\hat{x}}^{b_0} I(v,t) \frac{dv}{v} \qquad (41) $$

From this, $S(r,t)$, $\partial S(\hat{x},t)/\partial r$ and $d\sigma_{rr}(\hat{x},t)/dt = \partial\sigma_{rr}(\hat{x},t)/\partial t + \dot{\hat{x}}\partial\sigma_{rr}(\hat{x},t)/\partial r$ can be computed.
The rate (creep closure rate) of material point $x$, which coincides with $\hat{x}(t)$ at time $t$, is

$$x^2 \dot{x}(t) = b^2 \dot{b}(t) = -\frac{b^3}{2} I(b,t) \quad (42)$$

## 6.2. Examples

<figure><img src="image_21.png" /><figcaption>Figure 19. Top-left: Volume rate after leaching phase (x-axis shows time from start of leaching). Top-right: Deviator at 0.5 year after start of leaching. Bottom-left: Deviator 1 year after start of leaching. Bottom-right: Deviator at 2 years after start of leaching.</figcaption></figure>

The values of the elastic constants and the parameters of the M-D law are the same as in Section 3.2.1. As the first example, we consider a cavern ($V = 200,000 \text{ m}^3$) leached out in one year. The geostatic pressure is $P_\infty = 20.9$ MPa. At $t = 0$, the initial cavern of radius $\hat{x}_0 = 0.1$ m is created, and the pressure at the cavern wall drops to halmostatic pressure, $P_c = 11.4$ MPa. The final cavern radius after the leaching phase is $\hat{x}_1 = 36.3$ m. Figure 19 shows the evolution of the volume-loss rate after the leaching phase, and the deviator computed at times $t = 0.5$ year (half-way through the leaching phase), $t = 1$ year (at the end of leaching), and $t = 2$ years (one year after the end of leaching). It is seen that during the leaching phase, methods A and B are not able to capture the behavior of the stress deviator obtained using the moving boundary
formulation. At the end of the leaching phase, method B provides an acceptable approximation for the deviator, but the volume rate is still about twice larger than that obtained with the moving boundary. Soon after the leaching phase is completed, method B provides an acceptable approximation for the volume-loss rate as well (in this example, from 1.2 years after the start of leaching). On the contrary, there is a notable difference between the solutions for the stress deviator obtained with method A and with the moving boundary method — even 2 years after the start of leaching. It also takes a longer time for the volume-loss rate obtained with method A to converge to the result of the moving boundary formulation.

In the second example, the results are compared for a 200,000-m³ cavern leached out quickly (in 1 year) or slowly (in 3 years). Figure 20 shows the evolution of the volume-loss rate after the leaching phase. At time $t = 3$ years, the volume rate for the case of slow leaching is more than twice larger than the volume rate for the case of fast leaching. The volume rates converge after about 5 to 6 years from the start of leaching.

<figure><img src="image_22.png" /><figcaption>Figure 20. Volume rate after leaching phase, for a 200,000-m³ cavern leached out in 1 year or 3 years.</figcaption></figure>

## Conclusion

A first-order integro-differential equation allows computing stress distribution in the salt mass and cavern-volume rate when the Poisson's ratio is assumed to be 0.5 and when an idealized spherical shape is considered. This equation allows prediction of several characteristic features of salt cavern behavior: long geometrical transient evolution;
onset of reverse creep and effective tensile stresses; large volume-loss rates when
cavern pressure is cycled or when the effect of pressure solution is taken into account;
and effects of cavern growth during solution mining. It also sheds some light on the
mathematical origin of these features.

Acknowledgement

This project has received funding from the Fuel Cells and Hydrogen 2 Joint
Undertaking (now Clean Hydrogen Partnership) under Grant Agreement No
101006751. This Joint Undertaking receives support from the European Union's
Horizon 2020 Research and Innovation programme, Hydrogen Europe and Hydrogen
Europe Research. Special thanks to Kathleen Sikora and Benoît Brouard.

References

1. Van Sambeek LL, Fossum A, Callahan G, Ratigan J. Salt Mechanics: Empirical and Theoretical Developments. Hidekate Kakihana, Reginald Hardy Jr H, Takeshi Hoshi, Ken Toyokura (eds.), Proc. 7th Symp. on Salt, Vol. I: Elsevier, Amsterdam, 1993, 127-134.

2. Wang L, Bérest P, Brouard B. Mechanical Behavior of Salt Caverns: Closed-Form Solutions vs Numerical Computations. Rock Mechanics and Rock Engineering. 2015; 48(6):2369-2382.

3. Van Sambeek LL, DiRienzo AL. Analytical Solutions for Stress Distributions and creep closure around Open-Holes or Caverns Using Multilinear Segmented Creep Laws. Proc. SMRI Fall Meeting, Salzburg, Austria, 2016, 225-238.

4. Cornet JS, Dabrowski M, Schmid DW. Long-term cavity closure in non-linear rocks. Geophysical Journal International 2017; 210: 1231-1243.
https://doi.org/10.1093/gji/ggx227.

5. Manivannan S, Bérest P. Transient closure of a cylindrical hole in a salt formation considered as a Norton-Hoff medium. Rock Mech. Rock Eng. 2019. https://doi.org/10.1007/s00603-018-1732-6.

6. Munson DE, Dawson PR. Salt constitutive modeling using mechanism maps. Proc. 1st Conf. Mech. Beh. of Salt. Trans Tech Pub, Clausthal-Zellerfeld, 1984, 717-737.

7. Marketos G, Spiers CJ, Govers R. Impact of rock salt creep law choice on subsidence calculations for hydrocarbon reservoirs overlain by evaporite caprocks. J. Geo. Res.: Solid Earth, 2016. 121:4249-4267.

8. Munson DE. Constitutive model of creep in rock salt applied to underground room closure. International Journal of Rock Mechanics and Mining Sciences, Vol. 34 (2), 1997, pp. 233-247, https://doi.org/10.1016/S0148-9062(96)00047-2.
9. Hugout B (1988). Mechanical behavior of salt cavities -in situ tests- model for calculating the cavity volume evolution. Proc. 2nd Conf. Mech. Beh. of Salt, Hannover, Trans Tech Pub., Clausthal-Zellerfeld, Germany, 1984,291-310.

10. Bérest P, Brouard B, Karimi-Jafari M, Van Sambeek L (2007). Transient behaviour of salt caverns. Interpretation of Mechanical Integrity Tests. Int. J. Rock Mech. Min. Sc. 44, 767-786.

11. Brouard B (1998). Sur le comportement des cavités salines, étude théorique et experimentation in situ. PhD Thesis, Ecole Polytechnique, Palaiseau, France, 253 p. [in French]

12. DeVries KL, Mellegard KD and Callahan GD (2003). Laboratory Testing in Support of a Bedded Salt Failure Criterion. Solution Mining Research Institute Fall Meeting, Chester, United Kingdom, October 5-8, 90-113.

13. Van Sambeek LL, Ratigan JL, Hansen FD. Dilatancy of Rock Salt in Laboratory Tests. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, Pergamon Press, 1993, Vol. 30, No. 7, pp. 735-738.

14. Bérest P, Sicsic P, Brouard B (2016). Thermomechanical effects of a depressurization in a CAES. Symp. Am. Rock Mech. Ass, Houston, Tx.

15. Denzau H, Rudolph F (1997). Field test for determining the convergence of a gas storage cavern under load conditions frequently changing between maximum and minimum pressure and its finite element modelling. Proc. SMRI Spring meeting, Cracow, Poland, 71-84.

16. Wallner M. (1988). Frac-Pressure risk for cavities in rock salt. In: Proc. 2nd Conf. Mech. Beh. of Salt:645-658. Clausthal-Zellerfeld: Trans Tech Pub.

17. Wawersik W, Stone CM (1989). A characterization of pressure records in inelastic rock demonstrated by hydraulic fracturing measurements in salt. Int. J. Rock Mech. Min. Sci.:613-627.

18. Bérest P, Manivannan S (2019). Behavior of a salt cavern when creep law is modified to account for low deviatoric stresses. In: Proc. ARMA 2019, New York, Paper 19-477.

19. Bérest P., Brouard B., Hévin G., Réveillère A. Tightness of salt caverns used for hydrogen storage. Proc. ARMA Symposium, Houston, Texas, USA, 20-23 June 2021. ARMA 21-1616.

20. Lestringant C, Bérest P, Brouard B (2010). Thermo-mechanical effects in compressed air storage (CAES). In: Proc. SMRI Fall Technical conference, Leipzig, Germany, 29-44.

21. Spiers CJ, Schutjens PMTM, Brzesowsky RH, Peach CJ, Liezenberg JL, Zwart HJ (1990). Experimental determination of the constitutive parameters governing creep of rocksalt by pressure solution. Geological Society Special Publication 54, Deformation Mechanisms, Rheology and Tectonics, The Geological Society, London, England. 54:215-227, https://doi.org/10.1144/GSL.SP.1990.054.01.21.
22. Urai JL, Spiers CJ (2007). The effect of grain boundary water on deformation mechanisms and rheology of rocksalt during long-term deformation. Proc. sixth Conf. on Mech. Beh. of salt. A. A. Balkema, Netherlands, 149–158.

23. Bérest P, Blum PA, Charpentier JP, Gharbi H., Valès F (2005). Very slow creep tests on rock sample, Int. J. Rock Mech. Min. Sc., 42, 569-576.

24. Bérest P, Brouard B, Brückner D, DeVries K, Gharbi H, Hévin G, Hofer G, Spiers C, Uraï J (2019). Very Slow Creep Tests on Salt Samples. Rock Mech Rock Eng. 2019. https://doi.org/10.1007/s00603-019-01778-9.

25. Bérest P, Brouard B, Karimi-Jafari M (2009). The effect of small deviatoric stresses on cavern creep behavior. In: Proc. 9th Int. Symp. on Salt. 574–589. http://www.brouard-consulting.com/?q=node/47.

26. Sobolik S (2021). Effect of the Addition of a Low Equivalent Stress Creep Mechanism to the Analysis of Geomechanical Behavior of the SPR West Hackberry Site. In: Proc. SMRI Fall Meeting, Galveston, Texas.