Issue 
Sci. Tech. Energ. Transition
Volume 77, 2022



Article Number  6  
Number of page(s)  12  
DOI  https://doi.org/10.2516/stet/2022001  
Published online  04 May 2022 
Regular Article
Modelling the operation of gas storage in salt caverns: numerical approaches and applications
Storengy, Immeuble Djinn, 12 rue Raoul Nordling, 92274 BoisColombes Cedex, France
^{*} Corresponding author: laurent.jeannin@storengy.com
Received:
14
September
2021
Accepted:
24
January
2022
This paper examines numerical approaches to model operation of gas storage in salt caverns. The emphasis is on taking into account the thermal exchanges between the well, the cavern, and the rock mass, as well as the modelling of the volume losses of the cavern, while keeping the model fast, simple, and operationally usable. The use of this model in the context of monitoring storage operations is illustrated.
Key words: Underground energy storage / Salt caverns / Numerical modelling
© The Author(s), published by EDP Sciences, 2022
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Underground gas storage in salt caverns, especially natural gas storage, is a mature technique of more than 50 years. Caverns are leached out in salt formations at a depth ranging between a few hundreds to around 1000 m and filled with gas. Due to the evolution of the gas market, caverns are now often operated with cycles shorter than an annual cycle. It is then very important to accurately predict the thermodynamic and mechanical behavior of caverns in order to optimize their exploitation. Salt caverns could be used, as part of the energy transition, to store other gases such as H_{2}, which may be produced by water electrolysis from renewable electricity [1].
The caverns represent volumes of a few tens to a few hundred thousand cubic meters. Injection and withdrawal of gas cause the gas to compress and expand respectively in the caverns. Gas is injected at a prescribed rate and temperature at the wellhead and withdrawn at a given rate. However, the gas pressure in the cavern is set between a minimum value to limit salt creep and a maximum value to remain well below the lithostatic pressure and avoid mechanical damage to the cavern wall. Pressure and temperature measurements show that the gas does not behave adiabatically in caverns; moreover, the volume loss and convergence of the cavern can be significant.
Numerical modelling [2–6] is an approach widely used and constitutes a support for the operation of caverns (to inventor the gas volume, to determine the risk of hydrates, etc.). The existing approaches are based on a 1D modeling of the well and 0D of the cavern and propose a simplified modeling of the heat exchanges with the rock mass and of volume losses of the cavern. This work differs from the previous approaches because it explicitly describes the diffusion of heat in the rock mass and the heat exchanges with wellbore (solving the problem with a decomposition domain algorithm). Moreover, geomechanical modeling to account for cavern volume losses is based on solving the mechanical equilibrium equation in the rock mass using a 1D finite element formulation and sequential coupling.
The physical and numerical choices made in order to have a usable and predictive tool are detailed in this paper. Its remainder is organized as follows. The physical model describing gas in the cavern and in the well connecting the cavern and the surface is first presented. The flow and heat exchange modelling and associated algorithms for the well are described in Section 3. A simplified description of heat losses is adopted where the heat transfer problem is approximated by an axisymmetric problem around the well. A specific Robin–Robin decomposition domain algorithm provides a robust and efficient convergence to describe the well flow for the range of physical parameters of interest including low flow rates and very small time steps. In the next section, the mechanical problem related to salt creep is analyzed. The fluid pressure in the cavern is always smaller than the lithostatic stress, which causes the cavern to converge. The volume losses are estimated from a spherical mechanical model in large transformation. An updated Lagrangian formulation of the mechanical problem solved in finite elements is used to estimate the cavern volume at each time step. Finally, the value of this approach is illustrated with various examples of applications such as the modelling of operation on a CH_{4} cavern and the study of H_{2} storage.
2 Modelling of thermodynamical gas behaviour of salt caverns
Inside the cavern and well, the gas does not behave adiabatically: it exchanges heat with the surrounding rock. It is not possible to accurately represent the cavern shape and heat exchanges in detail, partly because of the imperfect knowledge of the geological environment and partly because a detailed threedimensional numerical model would be far too timeconsuming. Therefore, an intermediate approach is adopted with a simplified description of the geometry and heat exchange of the well and cavern.
In the following, the gas verifies:(1)where T is the gas temperature, P the gas pressure, Z(P, T) the compressibility factor and ρ(P, T) denotes the gas density in the well.
2.1 Cavern description
In this study, the cavern is represented by a spherical volume V_{c} and exchange surface S_{c} with the rock mass. Gas pressure P_{c} and temperature T_{c} inside the cavern are assumed uniform. Equations describing the thermodynamical evolution of the gas are integrated on the cavern volume.
The mass conservation is given by:(2) Q denotes the mass flow rate of gas injected (resp. positive) or withdrawn (resp. negative) in the cavern.
Energy balance equation in the cavern can be written:(3)with,where γ denotes a shape surface factor, as the real cavern is not spherical. T_{wall,c} represents the temperature at the cavern wall, T_{shoe} the gas temperature at the well bottom. c_{ p } is the gas specific heat capacity at constant pressure and h′ a heat transfer surface coefficient.
The last term of the r.h.s of equation (3) represents the heat input during injection, while the first term describes heat exchange with the rock mass.
In the rock mass, the temperature field T_{rock} verifies the heat equation:(4)where ρ_{ r } the rock density, the specific heat capacity of the rock per unit mass and λ_{rock} the rock thermal conductivity.
The temperature distribution is assumed to be spherically symmetric around the cavern. The rock mass around the cavern is then discretized in spherical shells (Fig. 1). The cavern problem with the rock mass around is solved as a one dimensional coupled nonlinear problem.
Fig. 1 Simplified model considered. 
2.2 Well flow
The well is taken vertical with radius r_{ w } and length L and identified by a vertical coordinate z oriented downward in this section.
The flow is one dimensional along the well and mass flow rate is supposed constant along the well (acoustic waves along the well are not considered). The gas velocity u verifies:(5)
The momentum equation along the well takes the following form:(6)where the last term of the r.h.s represents pressure drop. λ is a friction coefficient estimated by the Colebrook equation and .
Denoting e the gas internal energy per mass unit, energy conservation, and momentum conservation are combined. It yields:(7)where ϕ_{ v } is the dissipation rate related to viscous effects and ϕ the heat flux.
This equation (7) becomes:(8)with,(9)and,(10) is gas heat capacity of the gas at constant volume per unit of mass. T_{wall}(z, t) denotes the temperature on the wellbore, at the interface between the well and the rock mass. T_{wall} as well as the rock mass temperature are assumed axisymmetric.
Note that the typical order of magnitude of gas velocity is in the range 1–15 m s^{−1}. Gas flow is thus turbulent.
h is the heat convective transfer coefficient given by:(11)where λ_{ g } is the gas conductivity and N_{ u } the Nusselt Number. Nu is estimated from correlation with Reynolds and Prandtl numbers.
Heat exchanges between the surrounding axisymmetric rock mass domain (around the vertical well) may not be neglected. The temperature field T_{rock} verifies the axisymmetric heat equation.
At the interface between the rock mass and the well (Fig. 2), the temperatures and thermal fluxes must be continuous:(12)
Fig. 2 The wellrock mass coupled problem assuming axisymmetry. 
The numerical method used for solving this rock masswell coupled problem is detailed in Section 3.
2.3 Solving the well and cavern
Assuming pressure equilibrium between well bottom pressure and cavern pressure, it becomes:(13) R_{c} is the cavern radius. The well + cavern problem is discretized by adopting finite volume discretization along the well and in the rock mass (see [7]) and formulated as a fixed point problem solved iteratively, when the gas is injected. During withdrawal, cavern problem is solved first; P_{shoe} = P_{c} + ρ(P_{c}, T_{c})gR_{c} and T_{shoe} = T_{c} are then imposed as boundary conditions at the well bottom.
Figure 3 summarizes the strategy chosen for modelling gascavern storage operations. The domain decomposition algorithm between the well and the rock mass is described in Section 3, while Section 4 gives an overview of the geomechanical model.
Fig. 3 Different coupling and modules used in the software. 
3 Modelling heat exchange between well and rock mass
We present in this section domain decomposition algorithm to solve the coupled wellrock mass model [7]. Domain decomposition methods solve iteratively at a given time the well and rock mass subproblems with transmission conditions chosen in order to obtain wellposed subproblems and a quick convergence to the coupled solution. Compared with a fully coupled algorithm, domain decomposition approaches have the advantage to allow the use of adapted solvers for each subproblem. Morevoer, it is modular. For example, for long withdrawals, the thermal influence of the rock mass can be neglected in the well (Sect. 5) and a a reasonable approximation consists of decoupling the well from the rock mass.
In order to explain the implemented algorithm, we consider a simplified version of the well problem. Assuming the coupled problem is axisymmetric, the rock T_{rock}, well T and wall T_{wall} temperatures verify the following linear system of equations, the gas velocity u being assumed constant:with continuity conditions at the wellbore:(14)
The two last equations are the continuity of temperature and heat flow at the wellbore (denoted wall).
In practice, the well and rock mass subproblems given by the two first equations are solved iteratively, until convergence to the fully coupled solution (for which temperature and heat flux of each domain at the wellbore are equal). The rate of convergence of domain decomposition methods is highly dependent on the boundary conditions applied at the interface between the two subdomains. A Robin–Robin algorithm achieves an efficient convergence over the range of physical parameters of interest for well modeling (including low flow rates and small time steps) [7]. It exchanges linear combinations of pressure and heat flux at the well interface with wellchosen coefficients.
Considering an Euler implicit integration scheme, we define the current time step Δt and:
Having defined T_{rock}, T_{wall}, T the current time step temperatures, , will denote rock and well temperatures at the previous time step.
The Robin–Robin algorithm with Robin coefficients β_{rock} ≥ 0 and β_{ff} ≥ 0, updates the temperatures (T ^{ k−1}, , ) at iteration k – 1 ≥ 0 of the decomposition domain algorithm by the temperatures (T ^{ k }, , ) at iteration k ≥ 0 of the decomposition domain algorithm verifying:(15)followed by the solution (T ^{ k }, ) of the following well subproblem:(16)
The convergence analysis [7] shows that a simple and efficient choice of Robin’s coefficients is:and,where K_{0} is the second kind modified Bessel function and K_{1} = −.
This algorithm studied in [7] allows a robust and modular approach to account for heat exchange between the well and the rock mass.
4 Finite transformation and mechanical modelling
The change in volume of the cavern during operation is assessed from a mechanical model in large transformation. In general, few volume measurements are available (made by echometry) and even a simple mechanical model allows the mechanical effects to be taken into account and ensures satisfactory modelling of the long term behaviour of caverns.
During the operation of the storage, the salt creeps, because the gas pressure is lower than the lithostatic stress, the deviatoric stresses are nonzero in the rock in the vicinity of the cavern.
To account for this volume variation, a spherically symmetric problem is considered and the associated large transformation problem is solved. This model allows the variation in volume of the cavern to be evaluated at each time step. In the following paragraph, we recall the numerical approach used. It is based on a large deformation approach detailed in [8] in the case of a poroplastic medium and implemented here for viscoplastic (and poroviscoplastic in Sect. 5) behaviour. One of the interest of the large strain formulation used is that it manipulates not the Piola–Kirchhoff tensor, but the Cauchy tensor which has a direct physical meaning [8].
4.1 Finite transformation
The displacement field of a material point of the solid rock mass is defined as the difference of vector position in the reference configuration frame between a final configuration of the solid and the reference configuration. The transformation gradient is defined as the gradient of the vector position in the reference configuration.
Following [8], it is assumed that the transformation gradient (which is a second order tensor) may be written as:(17)where is the elastic component and the viscoplastic component. Assuming elasticity is isotropic, the unstress configuration may be defined [8]. The viscoplastic deformation is then associated with a residual deformation, when the elementary volume considered is unloaded. Elastic deformation is assumed infinitesimal related to the unstress configuration. In other words, large deformations are considered viscoplastic in nature in our modelling. The deformation gradient is thus linearized and is given by:(18) is the second order identity tensor, the linear elastic deformation tensor with respect to the unstressed configuration. The deformation rate is defined by:(19)
Considering a constant elastic stiffness tensor , the macroscopic stress evolution is given:(20)which is usually written in term of the Jaumann derivative:(21)with,(22) , which denotes the Jaumann derivative, may be interpreted as the derivative in the reference frame rotating at velocity relative to the observer’s reference frame. It merges with the usual derivative in the one dimensional spherical case considered in this study. is related to Cauchy stress using a viscoplastic behaviour law.
4.2 Updated Lagrangian formulation
The Updated Lagrangian method commonly used consists in using the last calculated configuration Ω_{ t } (and not the initial one Ω_{0}), as the reference configuration for resolution of the mechanical problem. The mechanical equilibrium, as well as behavior laws, at each time are written on a final unknown configuration (Fig. 4).
Let us now define u as the displacement vector on Ω_{ t } between the configuration at time t and t′ (Fig. 4). The gradient of the transformation between the configuration at time t and t′, estimated on the configuration at time t is given by:(23)
The mechanical equilibrium at time t′ verifies on the configuration :(24) ρ_{ r } is the rock density (for a porous medium see [8]).
Following [8], we define T the stress vector acting on S ^{ T }, the part of ∂Ω_{ t } where the stress vector is imposed. The weak form of the equilibrium equation takes the final form:(25)where,(26)and represents the part of the boundary surface where surface forces are imposed, v is a virtual displacement field defined on the configuration (see [8] for more details).
is estimated from the behaviour law, which will take the following form in our spherical case:(27)where represents the stress value on the configuration at time t transported on the configuration at time t′.
, the rate deformation tensor is then approximated by:(28)
This approach is implemented using first order finite elements on the 1D spherical mesh around the cavern defined in Section 2 (Fig. 1). It takes a particularly simple form in the 1D spherical case. The geometric nonlinearities are solved by a fixed point method at each time step [8].
The displacement field, solution of the mechanical problem, is therefore used sequentially at each time iteration of the well+cavern thermodynamic model. Volume variation of the cavern is estimated from the cavern wall displacement. At iteration in time n, the cavern model calculates the pressure at the wall. This wall pressure is used as a boundary condition in the mechanical model between times n and n + 1 to estimate a volume variation, which will be used in the cavern model between times n and n + 1.The time step coupling is imposed by the well + cavern thermodynamic model. This approach makes it possible to account for volume variations related to large withdrawals over short periods (a few days) in the cavern, the timestep being of the order of a few hours.
4.3 Viscoplastic behaviour law
Many creep laws have been developed for rock salt to relate as a function of the stress history. The standard models describing transient and/ or stationary creep strain are Norton–Hoff and Lemaître behaviour laws, usually used for cavern design.
The Norton law is defined by:(29) denotes the deviatoric component of the Cauchy stress. n, K are the parameters of the behaviour law. J_{2} is the second invariant of the deviatoric stress tensor.
Lemaître law is an hardening law, which takes the following form:(30)where ζ is an hardening variable defined by,(31)
The Norton–Hoff law is indeed a specific case of the Lemaître law with α = 1 and β = n.
5 Applications to operations modelling
5.1 Modeling operations on a methane cavern
From the operator’s point of view, questions arise throughout the exploitation of caverns concerning the thermodynamic evolution of the gas, volume losses and safety issues. In the short term, the thermodynamic conditions in the well and cavern during operation must be predicted [9]. Questions asked by the operator are for example:

what is the working gas volume (gas volume which can be withdrawn from a storage cavern) and what flow rates can be injected without exceeding the maximum allowable pressure in the cavern?

what is the maximum flow rate that can be withdrawn during withdrawal without the risk of hydrate formation on the surface?

how much and how quickly does the working gas volume increase after withdrawal due to heating from heat exchange with the rock mass?
Longterm questions concern for example studies of exploitation scenarios and evolution of the free gas volume; time at low pressure should be minimized to avoid volume losses due to saltcreep.
To answer these questions numerical models are created with geological, geometrical, and flowrates data as input. Those models are historymatched with continuous or oneoff measurements as wellhead or cavern pressure and temperature and free volume of the cavern.
From a practical point of view, each cavern has a historymatched model, which allows shortterm predictions to be made. The calculations are launched and interpreted using a web interface usable by the operator.
5.1.1 History matching and prediction
Figure 5 shows an example of historymatch during gas storage operation. The cavern shape factor γ times h′ the heat transfer surface coefficient and the cavern volume can be used as calibration parameters. Wellhead pressure and temperature are monitored and modelled; they are plotted on Figure 5 for a methane cavern in operation. Note that the temperature estimated by numerical modeling is compared to measurements only when gas is withdrawn. Indeed, when the cavern is at rest, wellhead temperature depends upon the weather, while during injection, it is an imposed value.
Fig. 5 Wellhead pressure and temperature – available data (grey marks) – estimated by model (black marks) – flow rate (black line) is imposed in the model. 
The historymatched model may be used to verify gas inventory and to accurately determine cavern performance. For example estimates of temperature and pressure at the wellhead during withdrawal are used to prevent hydrate formation near the wellhead (Fig. 6), by adapting flowrate.
Fig. 6 Hydrate risk prevention: hydrate formation stability curve (black line) – Estimate temperature and pressure at the wellhead during a withdrawal (blue points). 
5.1.2 Numerical examples: Mechanical coupling and thermal exchange
The next two examples are numerical models which illustrate respectively the interest of accounting for heat exchanges with the rock mass in the well for short time tests and of having a mechanical model to update the cavern volume.

Accounting thermal exchange in well
For relatively shorttime tests, it is important to take into account the heat exchange between the well and the rock mass. In this example, a cavern with a volume of 275 000 m^{3}, an initial temperature of 62 °C and an initial pressure of 240 bar in the cavern is considered. The rock thermal properties in the cavern model are the λ = 5.5 W·K^{−1}·m^{−1}, c_{ p } = 920 J·kg^{−1}·K^{−1}, ρ_{rock} = 2100 kg·m^{−3}, while along the well λ = 5 W·K^{−1}·m^{−1}, c_{ p } = 850 J·kg^{−1}·K^{−1}, ρ_{rock} = 2100 kg·m^{−3}.
An operating sequence consisting of a withdrawal (during three days), a cavern rest (two days) and an injection (three days) is considered (Fig. 7); creep is neglected. The temperature at the wellhead, when heat exchange with the rock mass is considered in the well, differs appreciably from an adiabatic solution (no heat exchange with rock mass) and a solution with an infinite rock heat capacity during the phases of withdrawal and rest. The injected gas temperature is 20 °C at the wellhead; the temperatures of these three solutions at the casing shoe (well bottom) are also different.

Mechanical coupling
Fig. 7 Wellhead temperatures during an operating sequence of injection and withdrawal – flow rate: grey line – temperature for coupled case: dotted – temperature for adiabatic case: dashed line – temperature for infinite rock mass conductivity: long dashed line. 
Some caverns can lose more than half their volume in 30 years of operation depending on the mechanical properties of the salt rock mass. Volume losses are particularly important, when the cavern is subjected to strong withdrawals and remains at low pressure.
We compare in this numerical example the effect of updating the cavern volume annually or at each iteration of the thermodynamic evaluation (see Sect. 3). Initial free gas volume of the cavern is 270 000 m^{3}, depth casing shoe at 1200 m. Initial lithostatic stress is 367.8 bar.
Young’s modulus of salt rock is 2 × 10^{5} bar, Poisson’s ratio is 0.2. The parameters of Lemaître law used for salt rock are α = 0.6, β = 0.53, K = 8.6 bar·day^{0.6}, while thermal behaviour is characterized by λ = 7 W·K^{−1}·m^{−1}, c_{ p } = 800 J·kg^{−1}·K^{−1}, ρ_{rock} = 2200 kg·m^{−3}.
One year leaching process is assumed. The following scenario of flowrate is used for both update approaches of the cavern volume: first 100 days, the cavern is at rest then every year, two cycles of withdraw (20 days), rest (14 days) and injection (60 days) are simulated. Cavern losses are of the order of 1% of volume per year.
Figure 8 shows caverns pressure calculated for both approaches. The first year the pressure difference is bigger than during the rest of the scenario because creep is higher at the beginning of the exploitation of the salt cavern. Afterwards, the pressure difference between the scenarios varies from 2 to 4 bar: an annual update does not allow for an accurate description of the effect of large withdrawals. In practice, the cavern volume is updated at each iteration of the wellcavern thermodynamic calculation.
Fig. 8 Geomechanical coupling: annual update of the cavern volume (dashed line), update of the volume at each time step (dotted line). 
5.2 Hydrogen storage
In this section we focus on the performance of hydrogen storage in a salt cavern. We compare hydrogene storage performance with that of two different gas storages, respectively methane and a mixture of 80 vol% methane and 20 vol% hydrogen, over the same withdrawal and injection cycle.
The cavern volume considered in this section is 190 000 m3, with a casing shoe at a depth of 1300 m. The initial rock mass temperature at cavern depth is 53 °C. The rock thermal properties in the cavern model are the λ = 5.5 W·K^{−1}·m^{−1}, c_{ p } = 920 J·kg^{−1}·K^{−1}, ρ_{rock} = 2100 kg·m^{−3}, while along the well λ = 5 W·K^{−1}·m^{−1}, c_{ p } = 850 J·kg^{−1}·K^{−1}, ρ_{rock} = 2100 kg·m^{−3}.
The minimum and maximum allowable pressures at the casing shoe (which is the top of the cavern) are 60 and 240 bar respectively. These are determined by geomechanical considerations.
5.2.1 Initial performance at the start of the storage
Two cycles were considered for a total simulated period of about two years (720 days). Each cycle consists in a first withdrawal phase of 180 days, from the maximum admissible pressure to the minimum pressure, followed by an injection phase of 180 days to reach the maximum pressure again. Volume losses are neglected. The injected gas was considered at a temperature of 30 °C.
For each of the gases studied, the constant flow rate of withdrawal (and injection), the working and cushion gas volumes (gas volume used to ensure the minimum storage pressure), the pressure evolutions at the well head and the temperature variations in the cavern are estimated using modelling (Tab. 1).
Working and cushion gas volumes.
Figure 9 represents well head pressure and cavern temperature. It is worth noting that the difference between wellhead pressure and cavern pressure is around 5–20 bar for methane or mix (Fig. 9) and around 0.5–2 bar for hydrogen. On one hand, hydrogen is less dense than methane. On the other hand, pressure losses by friction are lower for hydrogen. Finally, it is important to highlight that less energy is stored in a cavern filled with hydrogen than in a cavern filled with methane or a mixture. H_{2} has a higher calorific value per mass compared to natural gas. But the density of hydrogen is about 8 times lower (under normal conditions, 0.090 kg/m^{3} for hydrogen versus 0.789 kg/m^{3} for natural gas). In energy per volume at a given pressure, natural gas is about 3–4 times more energetic than hydrogen [1]. However, it should be noted that an hydrogen cavern will not be operated in the same mode, i.e. with the same frequencies and quantities withdrawn and injected as for a natural gas cavern [10].
Fig. 9 Top: Variations of wellhead pressure for methane (solid line), hydrogen (dotted line) and mix (dashdotted line). Down: Variations of cavern temperature for methane (solid line), hydrogen (dotted line) and mix (dashdotted line). 
5.3 Other applications
5.3.1 Cavern conversion
The conversion from Lgas to Hgas is a large infrastructure project in the gas industry in Europe and concerns some storages for example in North Europe. The Lgas (with low calorific value) will be replaced by Hgas (with high calorific value). To model this conversion, the thermodynamic behaviour of the gas mixture is estimated from Dalton’s law and the respective thermodynamic properties of Lgas and Hgas. It is assumed that the mixing during the injection in the cavern is instantaneous, which is a good approximation for caverns where the gas is in a convective state [11]. The standstill and withdrawal phases are characterized by unchanged gas composition. The cavern volume considered in this section is 554 000 m^{3}, with a casing shoe at a depth of 1251 m. The rock mass temperature at cavern depth is equal to 51 °C. The rock thermal properties in the cavern model are the λ = 6 W·K^{−1}·m^{−1}, c_{ p } = 920 J·kg^{−1}·K^{−1}, ρ_{rock} = 2100 kg·m^{−3}, while along the well λ = 5 W·K^{−1}·m^{−1}, c_{ p } = 850 J·kg^{−1}·K^{−1}, ρ_{rock} = 2100 kg·m^{−3}. Creep is assumed negligible.
The minimum and maximum wellhead pressures are 60 and 207 bar respectively during routine operation of the cavern. The pressure has been lowered to 45 bar exceptionally for the last Lgas withdrawal before the Hgas injection starts. Wellhead pressure and temperature and gas composition are predicted by the numerical model. Figure 10 shows the evolution of the gas composition during the L/H gas conversion. At the end of the first injection cycle, the cavity contains a mass fraction of 70% of Hgas (Fig. 10, point A); the first Hgas injection cycle makes the most significant change in composition. The mass fraction exceeds 90% by the next two cycles (Fig. 10, point B). Note that the inventory of Hgas is higher than that of Lgas at the same pressure (due to their different thermodynamic properties).
Fig. 10 Hgas mass fraction as function of the wellhead pressure during conversion. 
5.3.2 Abandonment and fluid transfer modeling in the rock mass
For storage operations, rock salt can be considered as an impermeable rock. However, when very longterm behavior is considered, the picture is modified: insitu abandonment tests of salt caverns demonstrate that the very low permeation of salt may be beneficial for the abandonment phase [12, 13]. In sum, salt permeability can prevent large pressure buildup in an abandoned cavern.
The generalization of this mechanical approach presented in Section 4 to the poromechanical case (in the case of cavern abandonment) is also available in the software; it is based on the formulation given in [8] for poromechanical medium. To account for small transfer of fluid through the cavern walls, the problem of transfer in the rock salt, considered a weakly permeable and porous medium, in finite transformation is solved with the poromechanical problem. A one way coupling method is used [14].
6 Conclusion
In this paper we have presented a tool for thermodynamic modelling of gas behaviour in caverns. Particular emphasis was placed on certain numerical aspects, namely the modelling of heat transfer between the fluid flowing in a well and a cavern and the description of volume losses.
This approach has been implemented in order to support operations on natural gas storage caverns and study the feasibilty of hydrogen storage in salt caverns.
This software is now a tool used on a daytoday basis to monitor, predict the behaviour of caverns, prevent operational issues and also to define commercial strategies in terms of storage capacity.
References
 Louvet F., Charnavel Y., Hevin G. (2017) Thermodynamic studies of hydrogen storage in salt caverns, in: SMRI Spring 2017, Albuquerque, April 2017. [Google Scholar]
 Tijani M. (2008) Contribution à l’étude thermomécanique des cavités réalisées par lessivage dans des formations géologiques salines, Habilitation à diriger des recherches, Université Pierre et Marie Curie, Paris VI. [Google Scholar]
 Rouabhi A., Hevin G., Soubeyran A., Labaune P., Louvet F. (2017) A multiphase multicomponent modeling approach of underground salt cavern storage, Geomech. Energy Environ. 12, 21–35. [CrossRef] [Google Scholar]
 KarimiJafari M. (2016) Thermodynamic simulation of gas caverns with GUSTS V.2, in: SMRI Fall 2016, Zalzburg, Austria, Sept. 2016. [Google Scholar]
 Nieland J.D. (2004) Salt cavern thermal simulator version 2.0 user’s manual, RESPEC Topical Report RSI1760, Ed. by RESPEC, Rapid City, South Dacota. [Google Scholar]
 Nieland J.D. (2008) Salt cavern thermodynamics – comparison between hydrogen, natural gas, and air storage, in: SMRI Fall 2008, Galveston, TX, USA, October 2008. [Google Scholar]
 Masson R., Jeannin L., Louvet F., Vuddamalay A. (2020) Domain decomposition methods to model heat between a well and a rock mass, Comput. Geosci. 24, 1377–1392. [CrossRef] [MathSciNet] [Google Scholar]
 Bernaud D., Deudé V., Dormieux L., Maghous S., Schmidt D.P. (2002) Evolution of elastic properties in finite poroplasticity and finite element analysis, Int. J. Numer. Anal. Methods Geomech. 26, 845–871. [CrossRef] [Google Scholar]
 Rodriguez P. (2021) Effect of heat transfer on gas storage cavern operations, in: Virtual Technical Conference, SMRI Spring 2021, April 2021. [Google Scholar]
 Evans J., Shaw T. (2021) Storage of hydrogen in solution mined salt caverns for long duration energy storage, in: Virtual Technical Conference, SMRI Spring 2021, April 2021. [Google Scholar]
 Berest P., Louvet F. (2020) Aspects of the thermodynamic behavior of salt caverns used for gas storage, OGST 75, 57. [CrossRef] [Google Scholar]
 Bérest P., Brouard B. (2003) Safety of salt caverns used for underground storage blow out; mechanical instability; seepage; cavern abandonment, Oil Gas Sci. Technol. 58, 3, 361–384. [CrossRef] [Google Scholar]
 KarimiJafari M. (2007) Comportement transitoire des cavités salines profondes, Thèse de doctorat, Ecole Polytechnique, Palaiseau. [Google Scholar]
 Jeannin L., Mainguy M., Masson R., VidalGilbert S. (2007) Accelerating the convergence of coupled geomechanicalreservoir simulations, Int. J. Anal. Numer. Methods Geomech. 31, 1167–1183. [Google Scholar]
All Tables
All Figures
Fig. 1 Simplified model considered. 

In the text 
Fig. 2 The wellrock mass coupled problem assuming axisymmetry. 

In the text 
Fig. 3 Different coupling and modules used in the software. 

In the text 
Fig. 4 Updated Lagrangian Formulation ([8]). 

In the text 
Fig. 5 Wellhead pressure and temperature – available data (grey marks) – estimated by model (black marks) – flow rate (black line) is imposed in the model. 

In the text 
Fig. 6 Hydrate risk prevention: hydrate formation stability curve (black line) – Estimate temperature and pressure at the wellhead during a withdrawal (blue points). 

In the text 
Fig. 7 Wellhead temperatures during an operating sequence of injection and withdrawal – flow rate: grey line – temperature for coupled case: dotted – temperature for adiabatic case: dashed line – temperature for infinite rock mass conductivity: long dashed line. 

In the text 
Fig. 8 Geomechanical coupling: annual update of the cavern volume (dashed line), update of the volume at each time step (dotted line). 

In the text 
Fig. 9 Top: Variations of wellhead pressure for methane (solid line), hydrogen (dotted line) and mix (dashdotted line). Down: Variations of cavern temperature for methane (solid line), hydrogen (dotted line) and mix (dashdotted line). 

In the text 
Fig. 10 Hgas mass fraction as function of the wellhead pressure during conversion. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.