# A multi-step calibration strategy for reliable parameter determination of
salt rock mechanics constitutive models

Hermínio T<EMAIL_ADDRESS>Maartje Houben  Kevin
Bisdom  Arjan van der Linden  Karin de Borst  Lambertus J. Sluys  Hadi
Hajibeygi  Faculty of Civil Engineering and Geosciences, Delft University of
Technology, Stevinweg 1, 2628CN, Delft, The Netherlands  Shell Global
Solutions International B.V., Grasweg 31, 1031 HW Amsterdam, The Netherlands

###### Abstract

The storage of renewable hydrogen in salt caverns requires fast injection and
production rates to cope with the imbalance between energy production and
consumption. Such operational conditions raise concerns about the mechanical
stability of salt caverns. Choosing an appropriate constitutive model for salt
mechanics is an important step in investigating this issue, and many
constitutive models with several parameters have been presented in the
literature. However, a robust calibration strategy to reliably determine which
model and which parameter set represent the given rock, based on stress-strain
data, remains an unsolved challenge. For the first time in the community, we
present a multi-step strategy to determine a single parameter set based on
many deformation datasets for salt rocks. Towards this end, we first develop a
comprehensive constitutive model able to capture all relevant nonlinear
deformation physics of transient, reverse, and steady-state creep. The
determination of the single set of representative material parameters is then
achieved by framing the calibration process as an optimization problem, for
which the global Particle Swarm Optimization algorithm is employed. Dynamic
data integration is achieved by a multi-step calibration strategy for a
situation where experiments are included one at a time, as they become
available. Additionally, our calibration strategy is made flexible to account
for mild heterogeneity between rock samples, resulting in a single set of
parameters that is representative of the deformation datasets. As a rigorous
mathematical analysis for the presented method and the lack of relevant
experimental datasets, we consider a wide range of synthetic experimental
data, inspired by the existing sparse relevant data in the literature. The
results of our performance analyses show that the proposed calibration
strategy is robust. Moreover, the model accuracy becomes increasingly better
as more data is included for calibration.

###### keywords:

Salt rock, Transient creep, Particle swarm, Viscoplasticity, Model
calibration, Underground hydrogen Storage

††journal: International Journal of Rock Mechanics and Mining Sciences

##  1 Introduction

Solution-mined salt caverns have been used for storage of hydrocarbons [1, 2,
3, 4], compressed air [5], and even feed-stock hydrogen mainly for chemical
industry [6, 7]. More recently, salt caverns are being considered as viable
options for large-scale (e.g., in the orders of TWh) storage of renewable
hydrogen to support energy transition [8]. However, due to the intermittent
nature of renewable sources and variable consumption rates, injections and
productions of renewable hydrogen are expected to occur at much more
unpredictable patterns than the current feed-stock-based storage systems. This
can result in fast pressure fluctuations inside the caverns, which naturally
increases the associated uncertainties related to the mechanical stability of
the caverns. To perform stability analyses, selecting an appropriate
constitutive model is the first crucial step. The constitutive model, used in
the numerical simulation framework, is expected to allow for accurate
prediction of the mechanical behavior of salt rock [9]. Constitutive models
entail many parameters, which are expected to be tuned for a specific
scenario. Therefore, an equally important step is the calibration of these
material parameters such that the combination of the appropriate model and
parameters leads to reliable safety assessments. Model calibration becomes
especially challenging when many experimental data sets are available, for
which the determination of a single representative set of parameters is
required. This particular point is the main focus of this study.

Salt rocks entail a complex mechanical behavior, with different yet co-
existing time-dependent deformation mechanisms [10]. Capturing all these
complexities in a single constitutive model is not a trivial task. In this
category, the Hou/Lux [11] and CDM [12] models are developed for both
transient and steady-state creep. In addition, they are also developed for
tertiary creep in the dilatancy zone and damage healing in the compressibility
zone. Conversely, other models disregard the compressibility/dilatancy regions
and focus on specific deformation mechanisms, such as transient, steady-state,
and tertiary creep [13, 14, 15, 16, 17, 18, 19]. However, as pointed out in
the literature [20], full constitutive models are not always necessary
depending on the specific application and the stress regimes. For instance,
transient creep might be disregarded for disposal operations, since internal
pressure is expected to be fairly stable. However, this is not the case for
hydrogen storage operations, where fast cyclic injections and productions are
expected to take place. Although some researchers consider transient creep as
a viscoelastic process [14, 13, 17], others develop constitutive models that
assume time-dependent inelastic deformations [21, 22, 11, 23]. When rapid
pressure swings are expected to occur, some studies in the literature also
include reverse creep in their modeling framework [24, 9].

The complexity of salt rock mechanics is often reflected in the constitutive
models, which tend to depend on many material parameters. As an example, the
viscoplastic part of the Hou/Lux-ODS and Hou/Lux-MDS [11] models depend on 11
and 18 material parameters, respectively. This makes the calibration strategy
a very challenging task, which also requires many experimental data sets to
characterize the material behavior. The viscoplastic model of Desai, for
instance, requires a set of at least six experiments for determining all
material parameters [21]. An additional challenge appears by recognizing the
potential inherent differences between samples and difficulties in
experimental controls [25]. In other words, different material parameters
might be found for the same batch of samples, thus requiring many sets of
experiments for an ideal calibration procedure. In this case, a set of
material parameters that is representative of all experiments should be
pursued through the calibration process. Moreover, it should be noted that
performing creep experiments in salt rocks is a time-consuming task [26, 27,
28], which means that the size of the dataset increases at a low pace.
Therefore, it would be convenient to perform partial calibrations with the
experimental data currently available and include new experiments as they
become available. During this process, the quality of the model results should
increase as more experiments are used for calibration, and it should stabilize
when a sufficient number of experiments is reached.

There are different approaches for determining the material parameters of a
constitutive model, some of which are complementary. As pointed out in [29],
some material parameters, such as Young’s modulus and Poisson’s ratio, can be
determined by well-defined procedures, while others are obtained either by
manual (trial and error) or automated procedures, with the latter being
preferred in face of a large number of parameters. A common approach for
automated calibration consists of defining an appropriate objective function
to be minimized in an optimization process. This is called direct calibration
[30], and although some researchers choose gradient-based algorithms [29],
meta-heuristic optimization approaches [25, 31] are more common choices as
they are more likely to avoid local minima. The direct calibration approach
treats the constitutive model as a black-box, so the optimization process is
agnostic to the field of application. Nevertheless, such an approach is not
found in the literature for constitutive models related to salt rocks.
Furthermore, a procedure for incrementally including new experiments in the
calibration process is also not reported in the literature.

In this context, the present work proposes a calibration strategy for a salt
rock constitutive model that includes deformation mechanisms relevant to
cyclic operations in salt caverns, such as transient, steady-state, and
reverse creep. The calibration of the material parameters related to each
deformation mechanism is presented. A multi-step direct calibration procedure
is developed to include one experiment at a time by properly adjusting the
objective (loss) function and solving it as a multi-objective function
optimization problem. Moreover, a regularization term is added to the loss
function to roughly favor the same fitting quality for all experiments. The
Particle Swarm Optimization (PSO) algorithm is employed for solving the
optimization problems. Moreover, we investigate a situation where each
experiment is performed on salt samples with slightly different material
properties. As discussed before, the goal is to obtain a single set of
representative material parameters. For this purpose, synthetic experiments
are employed, such that quality assessment of the proposed calibration
strategy is possible. Finally, a sensibility analysis is carried out to
investigate the influence of some material parameters and to clarify the
possible difficulties when performing the optimization process.

The remainder of this paper is organized as follows. Section 2 presents
experimental results obtained from a triaxial test on a salt rock. The
experimental results serve as a guide for the choice of our constitutive
model, which is presented in Section 3. Once the model is defined, the
calibration strategy adopted in this work is discussed in Section 4. In
Section 5, the model validation against an experimental result is presented,
and the calibration strategy is tested and investigated in different
situations. Finally, Section 7 closes our presentation.

##  2 Experimental methodology and results

To give additional insights into the mechanical behavior of salt rock and to
provide a benchmark for the constitutive model employed in this work, we
briefly discuss in this section a triaxial test performed on a salt rock. The
experiment was performed on the Z3 Leine rock salt provided by the Institut
fur Gebrismechanik GMBH (IFG), Germany. The Z3 Leine rock salt originated from
the Bernburg mine in Germany, it is very homogeneous ($>98\%$ NaCl),
colorless, and no obvious layering is visible. The grain size is less than 10
mm and the core sample was prepared as a 1-inch diameter and 2-inch long
cylindrical shape. Additionally, the sample was stored inside an airtight bag
and kept in a watertight container.

We show here the methods and results of a cyclic loading test on the Z3 Leine
rock salt. The experimental data is generated as testing data to assess the
ability of our constitutive computer model to describe salt behavior under
cyclic loading conditions. The experiment was carried out in a standard
triaxial cell operating at room temperature (21°C). The choice for room
temperature testing was made to keep the experiment as simple as possible. The
triaxial cell was built in-house (Energy Transition Campus Amsterdam, ETCA)
with an axial and radial stress limit of 95 MPa. The axial and radial
displacements were measured during the test. We first loaded the cylindrical
sample isostatically to 12.9 MPa and let the samples consolidate until the
measured axial and radial strain equilibrated. After equilibrium, the radial
stress was kept at 12.9 MPa for the duration of the experiment, whereas the
axial stress was increased to 23.5 MPa (at 2.0 MPa$/$hr) and kept at this
stress for 4 days and thereafter, decreased back down to 12.9 MPa (at 2.0
MPa$/$hr) followed by a consolidation for another 2 days. The 4 days of high
axial stress and 2 days of low axial stress were repeated for other four axial
stress steps, namely, 26.0/12.9 MPa, 27.5/12.9 MPa, 31.0/12.9 MPa, 33.5/12.9
MPa (see Fig. 1). Figure 1 also shows axial and radial strains, both measured
during the experiment, and volumetric strain (calculated) versus time. In this
case, all five differential stress steps show axial shortening and volumetric
shrinkage, whereas the radial strain shows dilation in all differential stress
steps.

Figure 1: Experimental results obtained from a triaxial test performed on the
Z3 Leine salt rock sample.

A closer look at this experiment can provide valuable insights into the
modeling of salt mechanics. For this purpose, the left graph of Fig. 2 shows a
zoomed view of the last loading stage (i.e. from 12.9 to 33.5 MPa). Notice
that the black dashed line indicates the previous maximum stress of 31.0 MPa.
When the axial stress goes from point A to point B, the corresponding strain
goes from point A’ to B’, as indicated by the dashed red curve. As can be
seen, the strain increment is relatively small along the stress path A-B.
However, when the axial stress exceeds the maximum stress ever applied to the
sample during the experiment, which at this point is 31.0 MPa, an abrupt
increase in strain is observed from point B’ to C’. This suggests that an
additional inelastic strain has been triggered during the stress path B-C.
This can be also concluded from the right graph of Fig. 2, which shows the
corresponding stress-strain relation during the loading path A-B-C. To
facilitate the visualization, notice that the light and dark gray areas in the
two graphs of Fig. 2 correspond. The right graph shows that from point A” to
B”, the salt sample behaves as a (visco)elastic material. However, from point
B” to C” there is a large increase in strain for a relatively small stress
increment. This is a characteristic of plastic (or viscoplastic) deformation,
and it can be observed in every unloading/reloading cycle during the
experiment.

Figure 2: Zoomed-in view of a reloading step of the experiment that allows for
identification of both viscoplastic and viscoelastic deformations.

Finally, the right graph of Fig. 2 shows a hysteretic effect during the
unload/reload path. The unload path goes from point B” to point A” and the
load path goes from point A” back to B”. The fact that these two paths are not
the same is an indication of the hysteretic effect. It is well-known that
viscoelastic materials present this type of behavior during loading/unloading
cycles. This indicates that, in addition to viscoplasticity, the salt rock
also presents a viscoelastic behavior. This hysteretic effect is actually a
manifestation of the well-known reverse creep.

##  3 Model formulation

Following the infinitesimal deformation assumption, the total strain can be
decomposed into different independent contributions by different mechanisms.
The first contribution accounts for the instantaneous elastic response.
Reverse creep, which is observed after an unloading step, is interpreted as a
time-dependent recoverable strain (i.e. viscoelastic) and is described by a
Kelvin-Voigt element [24]. Transient creep, on the other hand, is regarded as
a time-dependent inelastic deformation (i.e., viscoplastic) [23]. Finally,
steady-state dislocation creep is captured by a dashpot element following a
power-law function. In this manner, the total strain can be written as

$\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}_{e}+\boldsymbol{\varepsilon}_{ve}+\boldsymbol{\varepsilon}_{vp}+\boldsymbol{\varepsilon}_{cr},$
(1)

where $\boldsymbol{\varepsilon}_{e}$, $\boldsymbol{\varepsilon}_{ve}$,
$\boldsymbol{\varepsilon}_{vp}$ and $\boldsymbol{\varepsilon}_{cr}$ denote the
elastic, viscoelastic, viscoplastic and steady-state creep strain tensors,
respectively. Next, the specific models for each of these contributions are
described.

###  3.1 Elastic strain

For convenience, we define a function of Poisson’s ratio, $\nu$, as
$\mathbb{C}=\mathbb{C}(\nu)$, such that
$\mathbb{C}:\mathbb{R}\rightarrow\mathbb{R}^{3\times 3\times 3\times 3}$ is a
4th-order tensor with elements given by

$\mathbb{C}_{ijkl}(\nu)=\frac{\nu}{(1+\nu)(1-2\nu)}\delta_{ij}\delta_{kl}+\frac{1}{2(1+\nu)}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}).$
(2)

The time-independent elastic strain tensor at any time $t_{i}$ can be
calculated with the linear elasticity, assuming isotropic properties, i.e.,

$\boldsymbol{\varepsilon}_{e}(t_{i})=\mathbb{C}_{1}^{-1}:\boldsymbol{\sigma}(t_{i}),$
(3)

where $\mathbb{C}_{1}$ is the constitutive 4th-order tensor, defined as

$\mathbb{C}_{1}=E_{1}\mathbb{C}(\nu_{1}).$ (4)

Here, $E_{1}$ and $\nu_{1}$ denote the Young’s modulus and Poisson’s ratio
associated to the elastic deformation, respectively.

###  3.2 Viscoelastic strain

The solution for the time-dependent elastic (i.e. viscoelastic) strain,
represented by a Kelvin-Voigt element, for continuous stress functions has
been studied in the literature [32]. However, experimental data sets are
discrete, not continuous. Therefore, in this work, we present a solution
strategy for discrete stress conditions.

Figure 3: Discrete stress path in time.

Towards this end, consider a discrete stress function as depicted in Fig. 3.
The viscoelastic strain (i.e., for the Kelvin-Voigt element) subjected to this
discrete stress condition reads

$\boldsymbol{\varepsilon}_{ve}(t_{i})=\mathbb{C}_{2}^{-1}(t_{i}):\boldsymbol{\sigma}_{0}+\sum_{j=0}^{i-1}\mathbb{C}_{2}^{-1}(t_{i}-t_{j}):\Delta\boldsymbol{\sigma}_{j+1},$
(5)

where

$\mathbb{C}_{2}^{-1}(t)=\frac{1}{E_{2}}\left(1-e^{-\frac{E_{2}}{\eta_{2}}t}\right)\mathbb{C}^{-1}(\nu_{2})$
(6)

holds. Here, $\boldsymbol{\sigma}_{0}$ represents the stress tensor at
$t=t_{0}$ and $\mathbb{C}(\nu_{2})$ is defined by Eq. (2). Moreover, $E_{2}$
and $\nu_{2}$ denote the Young’s modulus and Poisson’s ratio of the spring in
the Kelvin-Voigt element, and $\eta_{2}$ is the dashpot viscosity of the same
element.

###  3.3 Viscoplastic strain

Following Perzyna’s theory [21, 23], the viscoplastic strain rate is expressed
as

$\dot{\boldsymbol{\varepsilon}}_{vp}=\mu_{1}\left\langle\dfrac{F_{vp}}{F_{0}}\right\rangle^{N_{1}}\dfrac{\partial
Q_{vp}}{\partial\boldsymbol{\sigma}},$ (7)

where $\mu_{1}$ and $N_{1}$ are material parameters, $F_{0}$ is a reference
value, and the operator $\left\langle\cdot\right\rangle$ denotes the ramp
function. The yield function [21] is given by

$F_{vp}(\alpha,\boldsymbol{\sigma})=J_{2}-(-\alpha I_{1}^{n_{1}}+\gamma
I_{1}^{2})\left[\exp{(\beta_{1}I_{1})}+\beta\dfrac{\sqrt{27}J_{3}}{2\sqrt{J_{2}^{3}}}\right]^{m}$
(8)

in which $n_{1}$, $\gamma$, $\beta_{1}$, $\beta$ and $m$ are material
parameters. In addition, an associative viscoplastic deformation mechanism is
considered, implying that the potential function $Q_{vp}$ in Eq. (7) is equal
to the yield function $F_{vp}$ of Eq. (8).

Moreover, denoting the deviatoric stress tensor as $\mathbf{s}$, the stress
invariants in Eq. (8) can be written as

$I_{1}=\text{tr}(\boldsymbol{\sigma})+\sigma_{t},\quad
J_{2}=\dfrac{1}{2}\mathbf{s}:\mathbf{s},\quad\text{and}\quad
J_{3}=\text{det}(\boldsymbol{\sigma}),$ (9)

where $\sigma_{t}$ is the tensile strength of the rock.

The hardening parameter in Eq. (8) characterizes the behavior of the yield
function. The hardening rule adopted in this work is

$\alpha=\dfrac{a_{1}}{\left[\left(\dfrac{a_{1}}{\alpha_{0}}\right)^{\frac{1}{\eta}}+\xi\right]^{\eta}},$
(10)

where $\eta$, $a_{1}$, $\alpha_{0}$ are material parameters. In particular,
$\alpha_{0}$ denotes the initial configuration of the yield surface.
Additionally, $\xi$ is the accumulated viscoplastic strain, which is given by

$\xi=\int_{t_{0}}^{t}\sqrt{\dot{\boldsymbol{\varepsilon}}_{vp}:\dot{\boldsymbol{\varepsilon}}_{vp}}\text{d}t.$
(11)

Note that, before any viscoplastic deformation takes place, Eq. (11) results
in $\xi=0$ and Eq. (10) leads to $\alpha=\alpha_{0}$, as expected.
Furthermore, $\alpha\rightarrow 0$ when $\xi\rightarrow\infty$, which causes
the yield surface to assume the position of the short-term failure boundary.
Figure 4 shows the yield surface (i.e. $F_{vp}(\alpha,\boldsymbol{\sigma})=0$)
for different values of hardening parameter $\alpha$.

Figure 4: Yield function for the viscoplastic model for different hardening
parameters.

Once the viscoplastic strain rate is obtained, the viscoplastic strain at time
$t_{i}$ is obtained by performing a time integration and using backward Euler
scheme, i.e.,

$\boldsymbol{\varepsilon}_{vp}(t_{i})=\boldsymbol{\varepsilon}_{vp}(t_{i-1})+\Delta
t\dot{\boldsymbol{\varepsilon}}_{vp}(t_{i}),$ (12)

in which $\Delta t=t_{i}-t_{i-1}$. Note that the viscoplastic strain rate is
also computed at $t_{i}$ due to the implicit formulation.

###  3.4 Dislocation creep strain

In this study, we neglect pressure solution creep and consider dislocation
creep only, as most of the available data in the literature are conducted
under deviatoric stresses exceeding 5 MPa [33]. Dislocation creep is often
modeled as using a power law [34], which reads

$\dot{\boldsymbol{\varepsilon}}_{cr}=A\exp\left(-\dfrac{Q}{RT}\right)q^{n-1}\mathbf{s},$
(13)

where $A$ and $n$ are material parameters, $Q$ is the activation energy, $T$
is temperature and $R$ is the universal gas constant ($R=8.32$ JK-1mol-1).

Similar to the viscoplastic strain, the creep strain at time $t_{i}$ is
obtained by

$\boldsymbol{\varepsilon}_{cr}(t_{i})=\boldsymbol{\varepsilon}_{cr}(t_{i-1})+\Delta
t\dot{\boldsymbol{\varepsilon}}_{cr}(t_{i}),$ (14)

where implicit time integration is also employed.

##  4 The model calibration strategy based on ensemble of data sets

In this section, we present the developed calibration strategy adopted to
determine the material parameters of the multi-physics constitutive model. As
summarized in Table 1, the multi-physics model entails 19 material
parameters222The variable $\alpha_{0}$ defines the initial position of the
yield surface. Therefore, it is not strictly a material parameter but rather
an initial condition. Nevertheless, for the purposes of this work,
$\alpha_{0}$ is also regarded as a material parameter. in total. Any manual
strategy to fit them would be certainly inapplicable. Moreover, performing
optimization in a 19th-dimensional space is cumbersome and prone to suffer
from the local optimum traps or too many possibilities to match the given data
set. As a remedy to this significant challenge, in this study, we propose a
strategy to adjust the material parameters separately in three different
groups as follows. The first group, i.e., group 1, comprises the parameters
associated with the dislocation creep model. The second group, i.e., group 2,
is the union of the material parameters representing both elastic and
viscoelastic models. Finally, the parameters related to the viscoplastic model
are considered in the third group, i.e., group 3. The sequence of calibration
goes from group 1 to 3. A brief discussion of the calibration steps of each
group is provided below.

Table 1: Material parameters to be calibrated. Group | Model | Parameters
---|---|---
1 | Dislocation creep |  $A$, $n$
2 | Elastic |  $E_{1}$, $\nu_{1}$
Viscoelastic |  $E_{2}$, $\nu_{2}$, $\eta_{2}$
3 | Viscoplastic |  $\mu_{1}$, $\eta$, $N_{1}$, $n_{1}$, $\beta$, $m$
$a_{1}$, $\beta_{1}$, $\gamma$, $\alpha_{0}$, $k$, $\sigma_{t}$

###  4.1 Group 1: Dislocation creep model

As expressed in Eq. (13), the creep model depends on 2 material parameters:
$A$ and $n$. Typical values of $n$ for dislocation creep range between 3 to 7
[35]. Based on this information, we adjust parameters $A$ and $n$ such that
the slopes of the creep strain curve match those from the experimental data.
For example, Figure (5) shows the axial strain measured from the experimental
results and the target slopes (stationary creep strain rates) represented by
the magenta lines. For this step, only the solution provided by Eq. (14) is
necessary, since only the slopes of the curves are required.

Figure 5: Steady state creep rates (slopes) represented by the magenta lines.

The material parameters found for the dislocation creep model are summarized
in Table 2.

Table 2: Material parameters for Group 1. Parameter | Value | Unit
---|---|---
$n$ | 3.0 | -
$A$ | $1.1\times 10^{-21}$ | Pa-ns-1

###  4.2 Group 2: Elastic and viscoelastic models

For the elastic model, the material parameters are the Young modulus and
Poisson’s ratio. The latter can be determined by the ratio between radial and
axial strains. Additionally, is it assumed that the Poisson’s ratio for
elastic and viscoelastic elements are the same, i.e., $\nu_{1}=\nu_{2}$. The
Young’s modulus, $E_{1}$, can be obtained by the slope $\Psi$ of the
unloading/loading path, as illustrated in Fig. 6. Finally, the values of
$E_{2}$ and $\eta_{2}$ are adjusted all together to match the distance $D$ and
the angle $\phi$ of the average unloading/loading slope.

Figure 6: Features ($\phi$ and $D$) observed when adjusting the elastic and
viscoelastic parameters.

For the elastic and viscoelastic contributions of the constitutive model, the
material parameters obtained based on the above described procedure are shown
in Table 3.

Table 3: Material parameters for Group 2. Parameter | Value | Unit
---|---|---
$E_{1}$ | 102 | GPa
$\nu_{1}$ | 0.32 | -
$E_{2}$ | 42 | GPa
$\nu_{2}$ | 0.32 | -
$\eta_{2}$ | $2.5\times 10^{5}$ | GPa.s

###  4.3 Group 3: Viscoplastic model

As summarized in Table 1, the viscoplastic model requires the definition of 12
parameters. Although different types of salt rocks and the presence of
impurities may produce different values for these material properties, typical
values for salt rocks can be found in the literature [23, 21]. The parameter
$\gamma$ is associated with the short-term failure boundary, whereas $\beta$,
$\beta_{1}$ and $m$ account for the influence of Lode’s angle. The position of
the dilatancy boundary can be adjusted mainly through parameters $\beta_{1}$
and $n_{1}$. Also, $F_{0}$ is a reference value and $k$ is regarded as zero,
meaning that an associative flow rule is adopted in this work. These
parameters are summarized in Table 4.

The remaining parameters $\mu_{1}$ and $N_{1}$ relate to the rate-dependent
behavior of the salt rock and should be calibrated for specific samples.
Similarly, the hardening rule described by Eq. (10) depends on $a_{1}$ and
$\eta$, which should also be calibrated. Additionally, the initial position of
the yield surface is defined by the parameter $\alpha_{0}$, which depends on
the stress history of the rock sample. The definition of $\alpha_{0}$ is very
important for the accuracy of the predictions, in the presence of
viscoplasticity.

Based on the above discussion, from the 12 material parameters associated with
the viscoplastic model, only $\mu_{1}$, $N_{1}$, $a_{1}$, $\eta$ and
$\alpha_{0}$ are to be determined. The remaining parameters are taken from
typical values found in the literature. Next, the strategy developed in this
work for obtaining these sets of parameters is described.

Table 4: Fixed material parameters for the viscoplastic element. Parameter | Value | Unit
---|---|---
$n_{1}$ | $3$ | $-$
$\beta_{1}$ | $4.8\times 10^{-3}$ | MPa-1
$\beta$ | $0.995$ | $-$
$m$ | $-0.5$ | $-$
$\gamma$ | $0.095$ | $-$
$F_{0}$ | $1.0$ | MPa2
$k$ | $0.0$ | $-$

###  4.4 Optimization strategy

As shown in the previous subsections, the constitutive model depends on 19
material parameters. Let us denote the entire set of material parameters as
$\textbf{m}\in\mathbb{R}^{d\times 1}$, where $d=19$. Furthermore, for the
purpose of the optimization procedure, let us split m in two subsets:
$\textbf{w}\in\mathbb{R}^{w\times 1}$, denoting the material parameters to be
optimized; and $\textbf{k}\in\mathbb{R}^{k\times 1}$ comprising the material
parameters that are fixed. Mathematically,
$\textbf{m}=\textbf{w}\cup\textbf{k}$, which implies that $w+k=d$. As
discussed before, the material parameters to be optimized are

$\textbf{w}=\begin{bmatrix}\mu_{1}&N_{1}&a_{1}&\eta&\alpha_{0}\end{bmatrix}^{T},$
(15)

whereas k comprises all the remaining material parameters, which are fixed
(Table 4). Note that the set of material parameters to be optimized
($\mathbf{w}$) belong exclusively to the viscoplastic contribution of the
constitutive model.

The strain tensor for a certain experiment $i$ is then obtained by

$\boldsymbol{\varepsilon}_{i}(\textbf{w},t)=f(\textbf{w},\textbf{k},\boldsymbol{\sigma}_{i}(t))\quad\text{for}\hskip
4.2679pti\in[1,\cdots,N],$ (16)

where $N$ is the number of experiments and $\boldsymbol{\sigma}_{i}(t)$ is the
specific stress condition applied to experiment number $i$.

For each experiment $i$, we define a loss function $F_{i}(\textbf{w})$ by
taking the Mean Absolute Percentage Error (MAPE) of the axial
($\varepsilon_{a}$) and radial ($\varepsilon_{r}$) strains, i.e.,

$F_{i}(\textbf{w})=\dfrac{1}{N_{t}}\sum_{j=1}^{N_{t}}\left(\left|\left|\frac{\varepsilon_{a,i}(\textbf{w},t_{j})-\varepsilon_{a}^{*}(t_{j})}{\varepsilon_{a}^{*}(t_{j})}\right|\right|+\left|\left|\frac{\varepsilon_{r,i}(\textbf{w},t_{j})-\varepsilon_{r}^{*}(t_{j})}{\varepsilon_{r}^{*}(t_{j})}\right|\right|\right),$
(17)

where $||\cdot||$ represents the absolute value, $N_{t}$ is the number of time
steps in experiment $i$, and $\boldsymbol{\varepsilon}^{*}(t_{j})$ is the true
deformation measured at time $t_{j}$. The MAPE metric represents an average of
the percentage error at each time $t_{i}$. By considering the percentage
error, the same importance is given to the errors in both the beginning and
the end of the experiment.

When more than one experiment is considered, say $i$ experiments, the loss
functions associated with each experiment are combined by taking the Mean
Squared Error (MSE) plus a regularization term, i.e.,

$L_{i}(\textbf{w})=\dfrac{1}{i}\sum_{j=1}^{i}F_{j}^{2}(\textbf{w})+\phi\sigma^{2},$
(18)

in which $\phi$ is the regularization factor and $\sigma^{2}$ is the variance
of the errors in all $i$ experiments. The incorporation of a regularization
term is designed to mitigate the risk of the optimization process becoming
overly specialized to a single experiment, thereby compromising its
generalizability to other experiments. This is further discussed in subsection
4.4.2.

Every time a new experiment is added, an optimization problem is solved to
obtain the best set of material properties for that particular experiment
only. For a given experiment $\\#i$, for example, this is achieved by
minimizing the corresponding loss function $F_{i}$. In other words, the set of
material parameters resulting from this process is defined as

$\textbf{w}_{i}^{F}\approx\text{argmin}\hskip
1.42262ptF_{i}(\textbf{w})\quad\forall\hskip 4.2679pt\textbf{w}\in\mathbb{P},$
(19)

where $\mathbb{P}\subset\mathbb{R}^{w\times n}$ is the initial population
distribution, with $n$ denoting the initial population (swarm) size. The
approximate sign “$\approx$” is used in Eq. 19 to admit the possibility of
obtaining sub-optimal or, at least, approximate solutions.

Once the single experiment optimization has been performed for experiment
$\\#i$, the set of parameters $\textbf{w}_{i}^{F}$ is added to the population
of the next full optimization, which considers all experiments from 1 to $i$.
That is,

$\textbf{w}_{i}^{L}\approx\text{argmin}\hskip
1.42262ptL_{i}(\textbf{w})\quad\forall\hskip
4.2679pt\textbf{w}\in\mathbb{P}\cup\mathbb{F}_{i}\cup\mathbb{L}_{i}.$ (20)

Here, $\mathbb{F}_{i}$ and $\mathbb{L}_{i}$ represent the union of the
solutions obtained from optimizing Equations (19) and (20), respectively. In
other words,

$\mathbb{F}_{i}=\bigcup_{j=1}^{i}\textbf{w}_{j}^{F}\quad\text{and}\quad\mathbb{L}_{i}=\bigcup_{j=1}^{i-1}\textbf{w}_{j}^{L}.$
(21)

Notice that the full optimization (Eq. 20) is performed with the initial
population plus the individuals contained in sets $\mathbb{F}_{i}$ and
$\mathbb{L}_{i}$. This means that the set of parameters obtained from the
previous full optimizations is also considered for the current full
optimization. In this manner, the calibrations of subsequent experiments do
not start from scratch, as they consider the results from previous
calibrations.

The complete procedure is described in Algorithm 1. Notice that $\mathbb{F}$
and $\mathbb{L}$ are initially empty sets. Moreover, when only one experiment
is available ($i=1$), there is no full optimization, so $\mathbb{L}$ remains
empty. After the second experiment has been added, full optimizations are
performed and $\bar{\textbf{w}}_{i}^{L}$ are added to $\mathbb{L}$. The
algorithm loops over $N$ experiments, including them one by one to the
calibration process.

Algorithm 1 Strategy for simultaneous calibration of a set of experiments.

$\mathbb{F}=\emptyset$

$\mathbb{L}=\emptyset$

for $i\leftarrow 1\cdots N$ do

$\textbf{w}_{i}^{F}\approx\text{argmin}\hskip
1.42262ptF_{i}(\textbf{w})\quad\forall\hskip 4.2679pt\textbf{w}\in\mathbb{P}$

$\mathbb{F}\leftarrow\mathbb{F}\cup\textbf{w}_{i}^{F}$

if $i>1$ then

$\bar{\textbf{w}}_{i}^{L}\approx\text{argmin}\hskip
1.42262ptL_{i}(\textbf{w})\quad\forall\hskip
4.2679pt\textbf{w}\in\mathbb{P}\cup\mathbb{F}\cup\mathbb{L}$

$\mathbb{L}\leftarrow\mathbb{L}\cup\textbf{w}_{i}^{L}$

end if

end for

####  4.4.1 Optimization algorithm

The solutions of the optimization problems (Equations 19 and 20) are obtained
by the Particle Swarm Optimization (PSO) algorithm [36, 37]. The PSO is a
meta-heuristic optimization algorithm based on population intelligence. An
initial population is created within the search space, and each member of the
population represents a possible solution. Each individual is assigned an
initial velocity so that they can move. For each time step (iteration), the
velocity of a certain particle $i$ is computed as

$\mathbf{v}^{t+1}_{i}=\omega\mathbf{v}^{t}_{i}+c_{1}\mathbf{r}_{1}\left(\mathbf{p}_{i}^{t}-\mathbf{x}_{i}^{t}\right)+c_{2}\mathbf{r}_{2}\left(\mathbf{g}_{i}^{t}-\mathbf{x}_{i}^{t}\right),$
(22)

where $\omega$, $c_{1}$ and $c_{2}$ are weights given to the previous velocity
(inertia), particle’s best position ($\mathbf{p}_{i}^{t}$) and swarm’s best
position ($\mathbf{g}_{i}^{t}$), respectively. Moreover, $\mathbf{r}_{1}$ and
$\mathbf{r}_{2}$ are random vectors between 0 and 1. With this new velocity,
the position of particle $i$ is updated by

$\mathbf{x}^{t+1}_{i}=\mathbf{x}^{t}_{i}+\mathbf{v}_{i}^{t+1}.$ (23)

Once all particles have their positions updated, the algorithm proceeds to the
next iteration. This process is repeated until a stop criterion is reached.

####  4.4.2 Multi-objective function optimization

Solving an optimization problem for multiple experiments simultaneously (Eq.
20) consists of multi-objective function optimization, since the goal is to
minimize the loss functions $F_{i}$ for all experiments. This type of problem
admits an infinite number of solutions that are not dominated by other
feasible solutions, which means that they are all acceptable solutions to the
minimization problem. This set of solutions is called the Pareto optimal front
and is illustrated in Fig. 7 for a situation where only two experiments are
considered.

Figure 7: Multi-objective optimization and Pareto front. Experiments #1 and #2
in this figure are only meant for illustration purposes.

Although all solutions along the Pareto front are admissible, we notice that
the solution at point A provides a good fit for experiment #1 and a sub-
optimal solution for experiment #2. Similarly, point C provides a better fit
for experiment #2 than for #1. On the other hand, point B provides a solution
that is equally good for both experiments, since $F_{1}\approx F_{2}$, and
therefore should be preferred.

To “encourage” the solution to stay along the dashed line of Fig. 7 we
penalize the loss function $L_{i}$ when $F_{1}$ $\cdots$ $F_{N}$ are too
different from each other, that is, when the variance is high (see Fig. 8).
Mathematically, this can be written as

$L_{i}(\mathbf{w})=\dfrac{1}{i}\sum_{j=1}^{i}F_{j}^{2}(\textbf{w})+\phi\underbrace{\dfrac{1}{i}\sum_{j=1}^{i}\left[F_{j}(\mathbf{w})-\tilde{F}(\mathbf{w})\right]}_{\sigma^{2}}=\dfrac{1}{i}\sum_{j=1}^{i}F_{j}^{2}(\textbf{w})+\phi\sigma^{2},$
(24)

where $\tilde{F}(\mathbf{w})$ is the average error of all experiments, as
indicated in Fig. 8, and $\phi$ is an arbitrary regularization factor. When
variance is high (Fig. 8-a), the regularization term is large, thus penalizing
the loss function $L_{i}$. When variance is close to zero (Fig. 8-b), which
means that the fit quality for all experiments is about the same, the
regularization term tends to vanish from Eq. (24).

Figure 8: Variance of $F_{i}$. (a) High variance, (b) low variance.

##  5 Results

In the first part of this section, we investigate the capabilities of the
proposed constitutive model in describing the mechanical behavior of salt
rocks observed in the laboratory. In the sequence, the optimization procedure
proposed for the model calibration is explored. We first start by presenting
the synthetic data employed throughout our analysis. Sensitivity analysis is
performed to identify the impact of each material parameter on the model
behavior. The PSO algorithm is first employed to calibrate each synthetic
experiment individually. Then, the calibration strategy is tested for fitting
the entire set of experiments as they become available, as described in
Algorithm 1.

###  5.1 Model validation

The primary purpose of this subsection is to show that the constitutive model
presented in Section 3 is able to accurately capture the time-dependent
behavior of salt rocks operating in the compressibility zone (i.e., no
tertiary creep) and under cyclic loading conditions. Additionally, we intend
to highlight the importance of each term considered in the proposed
constitutive model. Of particular interest is the description of transient
creep, as it is permanently present during cyclic operations. For this
purpose, we consider three models with different characteristics regarding the
transient creep description, as shown in Figure 9. In Model 1, the transient
creep deformations are considered to be fully recoverable (elastic). They thus
can be described by the Kelvin-Voigt element presented in subsection 3.2 (see
Fig. 9). On the other extreme, the hypothesis that all the deformation
observed during the primary stage of creep is fully inelastic is assumed by
Model 3, in which the viscoplastic model presented in subsection 3.3 is
employed. Finally, a combination of both hypotheses is assumed by Model 2,
where the viscoelastic element (Voigt element) is combined with the
viscoplastic element. It should be stressed that the terminologies “fully-
elastic”, “hybrid” and “fully-inelastic”, employed in Fig. 9 to designate the
models, refer exclusively to the underlying assumption of each model regarding
the transient creep stage. For instance, Model 1 is not fully elastic as a
whole, but the transient creep stage is regarded as such.

Figure 9: Models with different assumptions regarding the transient creep.

The models illustrated in Fig. 9 are now used to fit laboratory experimental
data. The experimental setup and salt rock details are presented in Section 2.
It can be verified from Fig. 10 that Model 1 can be calibrated to properly
describe the transient creep during the first step load. However, the
assumption that transient creep strains are fully recoverable implies that
most of the deformation should be reversed when the load is removed. This is
precisely what we observe in the zoomed-in graph on the left of Fig. 10, where
the unloading takes place at around 100 hours. As expected, the total strain
obtained with Model 1 is almost fully recovered after unloading, but this is
not observed in the experimental results. In fact, only a small strain
decrease is observed when the load is removed. Furthermore, the stress-strain
graph in Fig. 11 shows that Model 1 provides a completely wrong behavior. This
shows that transient creep is not fully elastic, and therefore cannot be
described by a simple Kelvin-Voigt element.

Figure 10: Total strain curves obtained with three different models.

It can be verified in Fig. 11 that Models 2 and 3 provide almost identical
results, and both of them are able to describe the experimental data. In spite
of that, the zoomed-in graph in Fig. 11 clearly shows that Model 3 is unable
to capture the reverse creep (hysteretic effect) observed during the
unloading/reloading step. By contrast, the Kelvin-Voigt element of Model 2 can
be tuned to properly describe reverse creep.

Figure 11: Stress-strain curves showing the unloading/loading step (reverse
creep).

##  6 Calibration procedure

In this section, we investigate the calibration procedure described in Section
4. A set of synthetic experiments is designed to test the proposed strategy.
The idea is to simulate a situation in which the experiments are provided one
at a time, and the number of experiments available is not sufficient for a
proper model calibration. In this context of high uncertainty, the goal is to
perform a model calibration as good as possible with the experiments available
and also to quantify the associated uncertainties.

###  6.1 Synthetic experiments

To overcome the lack of available laboratory experiments, in this work we
consider 6 synthetic experiments in which the correct model parameters are
known, thus providing the possibility of quantifying the calibration accuracy.
The idea is to simulate a situation, where the experiments are made available
one at a time, and the calibration takes place with the available data. This
means that the calibration is first performed on Exp. #1. Then, Exp. #2 is
made available and calibration is performed considering experiments #1 and #2.
This is repeated until all six synthetic experiments are available.

The synthetic experiments are shown in Fig. 12. The stress conditions were
chosen as shown in the left column of the figure. The graphs in the middle
column show the stress path of each experiment with the colors representing
time. In these same graphs, the initial yield surface of the viscoplastic
model is also plotted for reference, according to Eq. (8). Finally, the
corresponding strain responses are shown in the right column of Fig. 12. In
all stress conditions shown in Fig. 12, the axial load is larger than the
confining pressure, which means that the same Lode angle ($\theta=60^{o}$) was
considered in all experiments.

Figure 12: Synthetic experiments. Graphs on the left show the axial and radial
stresses for each experiment. The middle column shows the corresponding stress
path in the $I_{1}-\sqrt{J_{2}}$ plane, and the initial yield surface is also
indicated according to Eq. (8) for $\alpha=\alpha_{0}$. The right column
presents the resulting deformations of the synthetic sample.

All experiments shown in Fig. 12 use the same set of material properties. The
dislocation creep (Group 1) and viscoelastic (Group 2) parameters are shown in
Tables 2 and 3, respectively. For the viscoplastic element, the material
parameters that are fixed (i.e., not obtained by calibration) are shown in
Table 4. Finally, Table 5 presents the material parameters meant to be
obtained through the calibration process. Even though these values are known,
they are regarded as unknown so that we can test the optimization procedure.

Table 5: Material properties for the synthetic experiments. Parameter | Value | Unit
---|---|---
$\mu_{1}$ | $6.898\times 10^{-12}$ | $s^{-1}$
$N_{1}$ | 3 | $-$
$a_{1}$ | $1.8\times 10^{-5}$ | MPa2-n
$\eta$ | $0.82$ | $-$
$\alpha_{0}$ | $0.002$ | $-$

###  6.2 Global sensitivity analysis (GSA)

Before proceeding with the model calibration, we investigate how the material
parameters presented in Table 5 affect the model results. For this purpose, we
perform global sensitivity analysis (GSA), where all parameters are varied at
the same time so that we can identify not just how each parameter affects the
results individually, but also the interactions between them [38].

The strategy consists of generating a large population ($\sim$40 thousand) of
material parameter sets. We build random uniform distributions for each
material parameter within specified ranges, as shown in Table 6. For each set
of parameters, we compute a loss function by comparing the resulting strains
with a base solution. The loss function is computed as in Eq. (17), with the
set of parameters summarized in Table 5 regarded as the base solution. In this
manner, we can analyze how the changes in the material parameters affect the
values of the loss function.

Table 6: Ranges of variation of each material property for the uniform distributions. Parameter | $\log{\mu_{1}}$ | $N_{1}$ | $\log{a_{1}}$ | $\eta$ | $\log{\alpha_{0}}$
---|---|---|---|---|---
Max. | -10 | 4 | -4 | 1.5 | -1.9
Min. | -12 | 2 | -6 | 0.6 | -2.8

There are many ways GSA can be performed. One of them is computing statistical
correlations between the material parameters and the loss function333We
actually use the logarithm of the loss function as it spans many different
orders of magnitude.. The Pearson and Spearman correlations are commonly
applied to identify linear and non-linear relationships, respectively, between
variables [39]. Additionally, univariate F-test is also able to capture linear
relations, whereas mutual information can capture any type of dependency. It
is important to stress that both Pearson’s correlation and F-test assume
normal distributions, which is not the case here. Nevertheless, the four
methodologies are employed to identify relationships between the material
parameters and the loss function. According to these methodologies, Fig. 13
reveals a strong relationship between the loss function and parameters
$a_{1}$, $\eta$ and $\alpha_{0}$, the first two related to the hardening rule,
and the last one defining the initial yield surface. All four methodologies
point to a weak dependency on the rate-dependent parameters $\mu_{1}$ and
$N_{1}$. These results are also confirmed by analyzing feature importance
through a machine learning model and column permutation, as presented in A.

Figure 13: Table of correlation of each material parameter with the loss
function. Darker and lighter colors refer to weaker and stronger correlations,
respectively. Figure 14: (Left) Loss function distribution. (Right) Spearman
correlations for loss function values lesser than 0.1.

The dataset is created by choosing random values within the ranges specified
in Table 6. As a consequence, the correlation between the material parameters
is strictly zero. When the simulations are performed on the dataset, the
resulting probability distribution of loss function values is shown in the
left graph of Fig. 14. We can see that there is a high probability of loss
function values between 0.1 and 1.0, which suggests the existence of local
minima. To further investigate this issue, we filter the loss function values
larger than 0.1 out of our dataset and compute Spearman’s correlations between
all material parameters. In this manner, we only consider parameter sets that
produce reasonably good results (loss function values smaller than 0.1). These
correlations are shown in the right graph of Fig. 14, where a strong
relationship between $a_{1}$ and $\eta$ is identified. Weaker dependencies are
also revealed between $\mu_{1}$, $N_{1}$ and $\alpha_{0}$. This shows that
certain combinations between these variables can produce relatively good
results.

A better way to visualize this effect is by computing Spearman’s coefficient
as we continuously filter out the loss function values from $10^{5}$ down to
$10^{-2}$. This is shown in Fig. 15, which clearly shows a significant
increment of dependency between $a_{1}$ and $\eta$ for loss function values
smaller than 1.0. To a lesser degree, similar behavior can be verified for the
pairs $\mu_{1}-N_{1}$ and $\mu_{1}-\alpha_{0}$, as already revealed by the
right graph of Fig. 14. It is important to emphasize that all combinations of
material parameters were obtained by random sampling, thus ensuring zero
correlation between them. What the results of Figures 14 and 15 reveal is that
all sets of parameters that produce a relatively good solution present a
certain (still unknown) correlation between $a_{1}$ and $\eta$. The sets of
material parameters that do not present this correlation do not produce good
results compared to the base case.

Figure 15: Spearman’s correlation with filter out loss values.

To understand the origin of the correlation between $a_{1}$ and $\eta$ we use
the hardening rule expressed by Eq. (10). In this equation, we consider
$\alpha_{0}=0.02$, $\eta=0.82$ and $a_{1}=1.8\times 10^{-5}$. Then we vary the
accumulated viscoplastic strain ($\xi$) from 0 to 0.0075, which are typical
values. The obtained values for the hardening parameter are shown by the black
curve on the left graph of Fig. 16. This curve is regarded as a base case.
Then we can choose different values of $a_{1}$ and $\eta$ and check how it
affects the hardening parameter $\alpha$ with respect to $\xi$. For each pair
$a_{1}-\eta$ we compute the MAPE and plot the color surface shown in the right
graph of Fig. 16. The true $a_{1}-\eta$ pair is indicated by $P^{*}$ in both
graphs of Fig. 16. We can see a region in which the values of $a_{1}$ and
$\eta$ ($P_{1}$ to $P_{6}$) produce hardening curves very similar to the
original $P^{*}$. In other words, whichever pair of values we take in this
region, a good result may be produced by the constitutive model. Moreover, the
existence of this region explains the high correlation between $a_{1}$ and
$\eta$ in the right graph of Fig. 14 and Fig. 15.

Figure 16: The left graph presents the behavior of $\alpha$ as a function of
$\xi$ for different values of $\eta$ and $a_{1}$. The right graph shows the
MAPE for each pair $\eta-a_{1}$ compared to the true value $P^{*}$.

The results obtained from the global sensitivity analysis, particularly the
ones in Fig. 16, give an idea about the behavior of the loss functions for the
optimization procedure. These results show the presence of many local minima
in which the optimization algorithm might be trapped, which allows the
possibility of obtaining sub-optimal solutions.

###  6.3 Single calibrations

The goal of this subsection is to test the Particle Swarm Optimization (PSO)
algorithm for finding the true material parameters of the synthetic
experiments presented in Subsection 6.1 (Fig. 12). The idea is to run
individual optimizations for each synthetic experiment and check whether the
calibration process can indeed find the global solution, which is shown in
Table 5. The hyperparameters used in the PSO algorithm (Eq. 22) are
$\omega=0.8$, $c_{1}=0.8$ and $c_{2}=0.4$, the algorithm is set to perform 100
iterations, and the population size is 1280. To maximize coverage of the
search domain, Latin Hypercube Sampling (LHS) is adopted and the ranges of
each material property are those presented in Table 6.

Figure 17: Synthetic experiments. (Left) Cobweb representation of the material
parameters. (Right) Comparison between the true synthetic experiments and the
calibrated model.

The calibration outcomes obtained for each synthetic experiment are shown in
Fig. 17. The left column of this figure shows a cobweb representation of the
material properties. The gray lines represent the initial population, the red
line denotes the set of properties obtained from the optimization process, and
the true values are plotted in blue for comparison. The rightmost vertical bar
in the cobwebs represents the loss function $F_{i}$ (MAPE), according to Eq.
17, and its corresponding values for the best solution are also indicated. The
right column shows the total strains obtained with the calibrated properties
for each experiment, and they all agree very well with the (synthetic)
experimental results. For experiments #3 and #5 the values of $a_{1}$ and
$N_{1}$ obtained from optimization differ from the true values. However, as
revealed by the GSA, these two parameters have a secondary importance compared
to the others, so these disparities do not have much impact on the behavior of
the constitutive model. The values obtained from the calibrations of each
experiment are summarized in Table 7, along with the corresponding errors.

Table 7: Material parameters obtained from the calibration process. | $\log{\mu_{1}}$ | $N_{1}$ | $\log{a_{1}}$ | $\eta$ | $\log{\alpha_{0}}$ | Loss ($F_{i}$)
---|---|---|---|---|---|---
Reference | -11.1613 | 3.0000 | -4.7447 | 0.8200 | -2.6990 | 0.0
Exp. #1 | -11.1414 | 2.9815 | -4.7450 | 0.8201 | -2.6991 | $6.02\times 10^{-8}$
Exp. #2 | -11.3425 | 3.1648 | -4.8521 | 0.8772 | -2.6993 | $1.47\times 10^{-5}$
Exp. #3 | -11.9118 | 3.9593 | -4.7436 | 0.8202 | -2.6951 | $2.39\times 10^{-4}$
Exp. #4 | -11.1942 | 2.9800 | -4.7543 | 0.8256 | -2.6987 | $7.13\times 10^{-6}$
Exp. #5 | -11.6245 | 3.4175 | -4.8789 | 0.8925 | -2.6972 | $7.54\times 10^{-5}$
Exp. #6 | -11.2350 | 3.0760 | -4.6827 | 0.7875 | -2.6978 | $1.33\times 10^{-6}$

###  6.4 Full calibrations

In this test case, we aim to test the proposed calibration strategy in a more
realistic scenario. In laboratory conditions, even if the rock samples are
taken from the same region with a similar appearance, it is not possible to
guarantee they will present exactly the same mechanical properties.
Additionally, in general, the experiments are performed sequentially, meaning
that they are not all available simultaneously. Considering that each creep
experiment takes at least 20 to 30 days, approximately, it would be convenient
to have partially calibrated models with the experiments at hand and include
new experiments to the calibration as they become available.

To reproduce such a scenario, we take the reference values presented in Table
5 and randomly disturb the material properties of each rock sample. In this
manner, all experiments are performed with rock samples with slightly
different material properties. The levels of perturbations are $\pm 5\%$ for
$\log{\mu_{1}}$ and $N_{1}$, $\pm 2\%$ for $\log{a_{1}}$ and $\eta$, and $\pm
1\%$ for $\log{\alpha_{0}}$. The material properties for the rock samples of
each experiment are presented in Fig. 18. The graphs in this same figure show
the total strains obtained with the reference properties (i.e. those of Table
5) and the modified properties. It can be verified by visual comparison that
the strain curves do not differ much from the reference case, but those small
property variations do cause noticeable differences. Since the material
properties of each experiment are different from each other, it is impossible
to perfectly fit all experiments. Nevertheless, the goal of the calibration is
to find a single set of material properties that fits the experiments as good
as possible. Additionally, the calibration is performed in steps, including
one new experiment at a time, according to Algorithm 1.

Figure 18: The black curves (Ref.) are obtained with the reference properties
of Table 5, whereas the blue curves (Mod.) are created by perturbing
(modifying) the reference properties. The modified material parameters are
indicated below each graph.

The result of the calibration procedure is presented in Fig. 19. The first
line of this figure represents step 1, in which only experiment #1 is
available for calibration. Once the model is calibrated against experiment #1,
the model is used to predict the behaviors of the other experiments444Even
though this would never be possible in real life, we use the remaining
experiments for prediction in order to check predictability.. Above each graph
is indicated the $F_{i}$ error, as in Eq. (17), for each experiment. The loss
function values $L_{i}$, given by Eq. (18), are presented on the left for the
calibrated (blue) and predicted (red) experiments. In the first step, only
experiment #1 is available, so a single optimization (Eq. 19) is performed
using the initial population $\mathbb{P}$ to obtain $\mathbf{w}^{F}_{1}$. When
experiment #2 is available in the second step, another single optimization,
with $\mathbb{P}$ as the initial population, is performed on it, which
provides the set of parameters $\mathbf{w}^{F}_{2}$. Still on the second step,
a combined optimization (Eq. 20) is performed upon experiments #1 and #2
including $\mathbf{w}^{F}_{1}$ and $\mathbf{w}^{F}_{2}$ among the initial
population $\mathbb{P}$. The combined optimization between experiments #1 and
#2 provides $\mathbf{w}^{L}_{2}$, which should be included in the combined
optimization of the subsequent steps as well. This whole process is summarized
in Table 8 for all calibration steps.

Figure 19: Calibration process including one experiment at a time. Table 8: Population sets $\mathbb{P}$, $\mathbb{F}$ and $\mathbb{L}$ during the calibration steps. Step | Initial population | $\mathbb{F}$ | $\mathbb{L}$ | Output
---|---|---|---|---
#1 | $\mathbb{P}$ | $\emptyset$ | $\emptyset$ | $\mathbf{w}^{F}_{1}$
#2 | $\mathbb{P}$ | $\mathbf{w}^{F}_{1}$ | $\emptyset$ |  $\mathbf{w}^{F}_{2}$, $\mathbf{w}^{L}_{2}$
#3 | $\mathbb{P}$ |  $\mathbf{w}^{F}_{1}$, $\mathbf{w}^{F}_{2}$ | $\mathbf{w}^{L}_{2}$ |  $\mathbf{w}^{F}_{3}$, $\mathbf{w}^{L}_{3}$
#4 | $\mathbb{P}$ |  $\mathbf{w}^{F}_{1}$, $\mathbf{w}^{F}_{2}$, $\mathbf{w}^{F}_{3}$ |  $\mathbf{w}^{L}_{2}$, $\mathbf{w}^{L}_{3}$ |  $\mathbf{w}^{F}_{4}$, $\mathbf{w}^{L}_{4}$
#5 | $\mathbb{P}$ |  $\mathbf{w}^{F}_{1}$, $\mathbf{w}^{F}_{2}$, $\mathbf{w}^{F}_{3}$, $\mathbf{w}^{F}_{4}$ |  $\mathbf{w}^{L}_{2}$, $\mathbf{w}^{L}_{3}$, $\mathbf{w}^{L}_{4}$ |  $\mathbf{w}^{F}_{5}$, $\mathbf{w}^{L}_{5}$
#6 | $\mathbb{P}$ |  $\mathbf{w}^{F}_{1}$, $\mathbf{w}^{F}_{2}$, $\mathbf{w}^{F}_{3}$, $\mathbf{w}^{F}_{4}$, $\mathbf{w}^{F}_{5}$ |  $\mathbf{w}^{L}_{2}$, $\mathbf{w}^{L}_{3}$, $\mathbf{w}^{L}_{4}$, $\mathbf{w}^{L}_{5}$ |  $\mathbf{w}^{F}_{6}$, $\mathbf{w}^{L}_{6}$

It can be verified from Fig. 19 that the combined loss function for the
calibrated experiments, $L_{cal}$, substantially increases as more experiments
are included in the calibration procedure. At the same time, the $L_{pred}$ is
observed to decrease, which suggests that the predictability of the calibrated
constitutive model improves by adding new experiments. Figure 20 plots these
two quantities by the number of experiments available for calibration. On the
left graph of this figure, we see that the loss function indeed increases in
the beginning but it remains more or less constant as more experiments are
added to the calibration. The reason for this is, as mentioned before, it is
impossible to simultaneously fit all experiments since they have different
material properties, so after a certain point the error is expected to
stabilize at a constant level. On the left graph, however, we see that the
more experiments we consider for calibration, the better the predictions of
the constitutive model.

Figure 20: Loss function values for calibration and predictions.

Equation (18) includes the parameter $\phi$ as a regularization in the loss
function for the combined optimization. The results presented in Fig. 19 were
obtained with a regularization factor of $\phi=10$. To investigate how the
regularization term affects the calibration process, the same test case is run
considering $\phi=0$, $\phi=10$ and $\phi=100$. The loss function values
$L_{i}$ are shown in Fig. 20 as a function of the number of experiments
considered for calibration. The left graph shows the results of the
calibration experiments, where the loss function increases until the third
experiment is added and stabilizes at approximately constant values. As
expected, the larger $\phi$ the larger $L_{i}$ since the regularization term
becomes larger. The graph on the right shows the loss function values for the
experiments left for prediction. It can be seen that the loss function values
decrease as more experiments are considered for calibration. In other words,
it shows that the constitutive model can make better predictions when more
experiments are considered during calibration. Additionally, it can be
verified that for 5 experiments considered during calibration, the prediction
of the 6th experiment is better when $\phi=10$. Since it might be difficult to
compare the different curves in Fig. 20 because of the regularization term, in
Fig. 21 only the MSE is shown. In this case, the comparisons are unambiguously
and we can see from the right graph that indeed better predictions are
obtained with $\phi=10$.

Figure 21: Mean squared errors (MSE) for calibration and predictions.

As explained before, the sets of properties chosen for each experiment were
obtained by randomly disturbing the reference properties of Table 5.
Therefore, it would be reasonable to expect that the reference properties
themselves would also provide good performance in fitting all six experiments
as well. In Fig. 22 we compare the total strains for all experiments
considering the calibrated (red curves) and the reference (blue curves) set of
material properties. The curves in black represent the true values of the
synthetic experiments. It can be verified that experiments #2, #3 and #4 are
better fitted by the calibrated parameters, whereas the reference properties
perform slightly better for experiments #1, #5 and #6. However, the mean
squared errors (MSE) shown at the bottom of this figure reveal that the
calibrated material parameters present an overall better fit.

Figure 22: Constitutive model results obtained with the material properties
obtained from calibration and the reference values from Table 5. The overbar
refers to the quantities related to the reference parameters.

The calibrated results presented in Fig. 22 are obtained with a regularization
parameter of $\phi=100$. Different regularization parameters, however, may
lead to different calibrated material properties. The obtained values for
$\phi$ equals to 0, 10 and 100 are summarized in Table 9, along with the
reference properties for comparison. As shown in this table, the MSE values of
the calibrated properties are always smaller than the reference properties.
Interestingly, the material properties obtained through the calibration
process are significantly different from each other and the reference values.
Apart from $\alpha_{0}$, the other material properties even lie outside the
limits of the perturbed properties, as it can be verified in Fig. 18. This can
be attributed to the non-linearity of the viscoplastic model.

Table 9: Calibrated material properties obtained with different values of $\phi$, and the reference properties from Table 5 for comparison. The last line shows the mean squared error (MSE) for each set of properties. | $\phi=0$ | $\phi=10$ | $\phi=100$ | Reference
---|---|---|---|---
$\mu_{1}$ | $4.51\times 10^{-12}$ | $1.83\times 10^{-12}$ | $9.77\times 10^{-12}$ | $6.89\times 10^{-12}$
$N_{1}$ | 2.98 | 3.99 | 3.47 | 3.00
$a_{1}$ | $4.33\times 10^{-5}$ | $3.12\times 10^{-6}$ | $2.19\times 10^{-5}$ | $1.80\times 10^{-5}$
$\eta$ | 0.60 | 1.19 | 0.74 | 0.82
$\alpha_{0}$ | 0.00207 | 0.00205 | 0.00202 | 0.0020
MSE | 0.0615 | 0.0635 | 0.0609 | 0.5167

##  7 Conclusion

Geomechanics of salt caverns for safe cyclic hydrogen storage requires a
detailed investigation of the associated creep mechanisms. For this purpose,
we assemble a constitutive model able to capture transient creep, reverse
creep and steady-state dislocation creep. This is achieved by including
elastic, viscoelastic, viscoplastic, and power-law contributions to the
constitutive model. The viscoelastic contribution can capture the hysteretic
effect that appears during unloading/reloading stages, also known as reverse
creep. Moreover, it is shown that the viscoelastic element alone is not able
to appropriately describe transient creep, for which a viscoplastic model is
indispensable. With this constitutive model, the experimental results
presented in this paper can be successfully described.

An important step in devising and deploying a constitutive model relates to
model calibration. Due to the time scales involved in creep deformations,
laboratory experiments in salt rocks are usually very time-consuming
(approximately one month long or even more). Additionally, it takes a lot of
experiments to properly calibrate a constitutive model. There is also the
possibility that mechanical properties might vary from sample to sample,
making the calibration process more challenging. In this context, the main
contribution of this paper is to devise a multi-step calibration strategy that
takes place on the fly, in the sense that the constitutive model is calibrated
against experiments as they are made available. Determination of the material
parameters corresponding to the elastic, viscoelastic, and dislocation creep
can be done independently because the presented procedure is general and
flexible. Another significant contribution of the present study is on the
calibration of the viscoplastic material parameters, which have to be all
tuned simultaneously. Moreover, the choice of the material parameters to be
calibrated is discussed, and a systematic analysis for the impact of these
parameters on the constitutive model behavior is also presented. The
calibration strategy consists of sequentially solving multi-objective
optimization problems as new experiments are included.

It is shown that the proposed strategy provides increasingly better results as
more experiments are incorporated into the calibration process. In this
manner, calibrated constitutive models are provided much more efficiently and
reliably than the existing procedures in the literature. Additionally, it is
shown that including a regularization term into the loss function is
beneficial to promote equal quality fit of the constitutive model against all
experimental data sets. The proposed methodology can be applied to any
constitutive model, and it provides representative material parameters for
salt rocks. As such, this work contributes to improving the accuracy of
constructing a suitable constitutive model with reliable material parameter
sets to study the mechanics of salt caverns under cyclic operations. Future
studies will incorporate this characterization procedure in field-scale 3D
cavern simulations.

## CRediT authorship contribution statement

H.T.H.: Conceptualization, Methodology, Software, Validation, Formal analysis,
Visualization, Writing – Original Draft. M.H.: Conceptualization,
Investigation, Writing – Review & Editing. K.B.: Conceptualization, Writing –
Review & Editing. A.L.: Investigation. K.B.: Conceptualization, Writing –
Review & Editing. L.S: Conceptualization, Methodology, Writing – Review &
Editing H.H.: Conceptualization, Methodology, Writing – Review & Editing.

## Declaration of competing 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.

## Data availability

Digital datasets of the results and input data are available upon request.

## Acknowledgments

This research was partly supported by Shell Global Solutions International B.V
within the project ‘SafeInCave’. The authors also acknowledge members of the
ADMIRE and DARSim research groups at TU Delft for the fruitful discussions
during the development of this work. Additionally, the authors acknowledge the
fruitful discussions with the technical staff of Shell Global Solutions
International B.V and for providing the experimental data set used in this
work.

##  Appendix A Feature importance for sensitivity analysis

Another interesting strategy to measure the influence of each material
parameter on the model behavior is through feature importance, as usually
performed in machine learning (ML) applications. Some ML models, such as
linear regression and decision trees, naturally provide feature importance
rankings. A more general approach can be achieved by column permutation. This
strategy consists of training a ML model to predict the target variable (in
our case, the loss function). After the training step, the model performance
is evaluated according to a certain metric (e.g. MSE). This establishes a base
metric for the ML model. Then, we shuffle the values of a certain feature
(i.e., material parameter), run the model again without retraining, and check
the new metric. Based on how much the new metric is reduced compared to the
base metric we can infer the importance of that feature. A great reduction
means the ML model heavily relies on that feature. Conversely, if the metric
barely changes, it means the ML model does not find that feature very useful.
This procedure is performed using the Random Forest algorithm and the results
are shown in Fig. 23. We can see that parameters $\mu_{1}$ and $N_{1}$ are of
less importance compared to $a_{1}$, $\eta$ and $\alpha_{0}$. This is in
complete agreement with the correlations found in Figure 13.

Figure 23: Feature importance according to the column permutation method.

## References

  * [1] JW Joachim Koenig.  Preparing motor gasolines for salt cavern storage of up to 10 years.  SMRI Fall Mtg, 425, 1994.
  * [2] Radosław Tarkowski, Barbara Uliasz-Misiak, and Piotr Tarkowski.  Storage of hydrogen, natural gas, and carbon dioxide–geological and legal conditions.  International Journal of Hydrogen Energy, 46(38):20010–20022, 2021\.
  * [3] Kermit Allen.  Eminence dome-natural-gas storage in salt comes of age.  Journal of Petroleum Technology, 24(11):1299–1301, 1972.
  * [4] HL Gentry.  Storage of high pressure natural gas in underground salt or rock caverns.  In Symposium on Salt, pages 604–608, 1963.
  * [5] Fritz Crotogino, Klaus-Uwe Mohmeyer, and Roland Scharf.  Huntorf caes: More than 20 years of successful operation.  Orlando, Florida, USA, 2001.
  * [6] Dilara Gulcin Caglayan, Nikolaus Weber, Heidi U Heinrichs, Jochen Linßen, Martin Robinius, Peter A Kukla, and Detlef Stolten.  Technical potential of salt caverns for hydrogen storage in europe.  International Journal of Hydrogen Energy, 45(11):6793–6805, 2020\.
  * [7] Radoslaw Tarkowski.  Underground hydrogen storage: Characteristics and prospects.  Renewable and Sustainable Energy Reviews, 105:86–94, 2019.
  * [8] Leila Hashemi, Martin Blunt, and Hadi Hajibeygi.  Pore-scale modelling and sensitivity analyses of hydrogen-brine multiphase flow in geological porous media.  Scientific reports, 11(1):8348, 2021.
  * [9] Michel Aubertin, Michel R Julien, Stéphane Servant, and Denis E Gill.  A rate-dependent model for the ductile behavior of salt rocks.  Canadian Geotechnical Journal, 36(4):660–674, 1999.
  * [10] Nicolae Cristescu and Udo Hunsche.  Time effects in rock mechanics, volume 350.  Wiley New York, 1998.
  * [11] Zhengmeng Hou.  Mechanical and hydraulic behavior of rock salt in the excavation disturbed zone around underground facilities.  International Journal of Rock Mechanics and Mining Sciences, 40(5):725–738, 2003.
  * [12] A Hampel.  The cdm constitutive model for the mechanical behavior of rock salt: Recent developments and extensions.  In Proceedings of the 7th Conference on Mechanical Behaviour of Salt, Paris, pages 16–19, 2012.
  * [13] Stefan Heusermann, Olaf Rolfs, and Uwe Schmidt.  Nonlinear finite-element analysis of solution mined storage caverns in rock salt using the lubby2 constitutive model.  Computers & structures, 81(8-11):629–638, 2003.
  * [14] HW Zhou, CP Wang, BB Han, and ZQ Duan.  A creep constitutive model for salt rock based on fractional derivatives.  International Journal of Rock Mechanics and Mining Sciences, 48(1):116–121, 2011.
  * [15] Qinghui Jiang, Yajing Qi, Zhijian Wang, and Chuangbing Zhou.  An extended nishihara model for the description of three stages of sandstone creep.  Geophysical Journal International, 193(2):841–854, 2013.
  * [16] Lin-jian Ma, Xin-yu Liu, Qin Fang, Hong-fa Xu, Hui-min Xia, Er-bing Li, Shi-gang Yang, and Wen-pei Li.  A new elasto-viscoplastic damage model combined with the generalized hoek–brown failure criterion for bedded rock salt and its application.  Rock mechanics and rock engineering, 46:53–66, 2013.
  * [17] Jianqiang Deng, Yaoru Liu, Qiang Yang, Wei Cui, Yinbang Zhu, Yi Liu, and Bingqi Li.  A viscoelastic, viscoplastic, and viscodamage constitutive model of salt rock for underground energy storage cavern.  Computers and Geotechnics, 119:103288, 2020.
  * [18] Pedro ALP Firme, Nuno B Brandao, Deane Roehl, and Celso Romanel.  Enhanced double-mechanism creep laws for salt rocks.  Acta Geotechnica, 13:1329–1340, 2018.
  * [19] Pedro ALP Firme, Deane Roehl, and Celso Romanel.  Salt caverns history and geomechanics towards future natural gas strategic storage in brazil.  Journal of Natural Gas Science and Engineering, 72:103006, 2019\.
  * [20] MB Dusseault, S Bachu, and L Rothenburg.  Sequestration of co2 in salt caverns.  Journal of Canadian Petroleum Technology, 43(11), 2004.
  * [21] CS Desai and A Varadarajan.  A constitutive model for quasi-static behavior of rock salt.  Journal of Geophysical Research: Solid Earth, 92(B11):11445–11456, 1987.
  * [22] Jishan Jin and ND Cristescu.  An elastic/viscoplastic model for transient creep of rock salt.  International Journal of Plasticity, 14(1-3):85–107, 1998.
  * [23] Kavan Khaledi, Elham Mahmoudi, Maria Datcheva, and Tom Schanz.  Stability and serviceability of underground energy storage caverns in rock salt subjected to mechanical cyclic loading.  International journal of rock mechanics and mining sciences, 86:115–131, 2016.
  * [24] Darrell E Munson, Kerry L DeVries, Arlo F Fossum, and Gary D Callahan.  Extension of the md model for treating stress drops in salt.  Technical report, Sandia National Labs., 1993.
  * [25] Joseph D Robson, Daniel Armstrong, Joseph Cordell, Daniel Pope, and Thomas F Flint.  Calibration of constitutive models using genetic algorithms.  Mechanics of Materials, 189:104881, 2024.
  * [26] Chunhe Yang, JJK Daemen, and Jian-Hua Yin.  Experimental investigation of creep behavior of salt rock.  International Journal of Rock Mechanics and Mining Sciences, 36(2):233–242, 1999.
  * [27] Hakan Özşen, İhsan Özkan, and Cem Şensöğüt.  Measurement and mathematical modelling of the creep behaviour of tuzköy rock salt.  International journal of rock mechanics and mining sciences, 66:128–135, 2014.
  * [28] P Bérest, JF Béraud, H Gharbi, B Brouard, and K DeVries.  A very slow creep test on an avery island salt sample.  Rock Mechanics and Rock Engineering, 48:2591–2602, 2015.
  * [29] A Anandarajah and D Agarwal.  Computer-aided calibration of a soil plasticity model.  International Journal for Numerical and Analytical Methods in Geomechanics, 15(12):835–856, 1991.
  * [30] Andreas Tolk, John Fowler, Guodong Shao, and Enver Yücesan.  Advances in modeling and simulation.  Springer International Publishing. Epub ahead of print, 10:978–3, 2017.
  * [31] Robin Schulte, Cavid Karca, Richard Ostwald, and Andreas Menzel.  Machine learning-assisted parameter identification for constitutive models based on concatenated loading path sequences.  European Journal of Mechanics-A/Solids, 98:104854, 2023.
  * [32] J Argyris, I St Doltsinis, et al.  Constitutive modelling and computation of non-linear viscoelastic solids. part i: Rheological models and numerical integration techniques.  Computer methods in applied mechanics and engineering, 88(2):135–163, 1991.
  * [33] Kishan Ramesh Kumar and Hadi Hajibeygi.  Influence of pressure solution and evaporate heterogeneity on the geo-mechanical behavior of salt caverns.  The Mechanical Behavior of Salt X, pages 407–420, 2022.
  * [34] Kishan Ramesh Kumar, Artur Makhmutov, Christopher J Spiers, and Hadi Hajibeygi.  Geomechanical simulation of energy storage in salt formations.  Scientific Reports, 11(1):19640, 2021.
  * [35] Luuk Hunfeld, Jaap Breunese, and Brecht Wassing.  The influence of a threshold stress for pressure solution creep on cavern convergence and subsidence behavior–an fem study.  In The Mechanical Behavior of Salt X, pages 577–589. CRC Press, 2022.
  * [36] Russell Eberhart and James Kennedy.  A new optimizer using particle swarm theory.  In MHS’95. Proceedings of the sixth international symposium on micro machine and human science, pages 39–43. Ieee, 1995.
  * [37] Ahmed G Gad.  Particle swarm optimization algorithm and its applications: a systematic review.  Archives of computational methods in engineering, 29(5):2531–2561, 2022.
  * [38] Francesca Pianosi, Keith Beven, Jim Freer, Jim W Hall, Jonathan Rougier, David B Stephenson, and Thorsten Wagener.  Sensitivity analysis of environmental models: A systematic review with practical workflow.  Environmental Modelling & Software, 79:214–232, 2016.
  * [39] Alessandro Rovetta.  Raiders of the lost correlation: a guide on using pearson and spearman coefficients to detect hidden correlations in medical sciences.  Cureus, 12(11), 2020.