© The Author(s), published by EDP Sciences, 2022

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₂, 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₄ cavern and the study of H₂ 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 three-dimensional numerical model would be far too time-consuming. 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⁻¹. 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 well-rock mass coupled problem assuming axisymmetry.

The numerical method used for solving this rock mass-well 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 gas-cavern 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 well-rock 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 well-posed 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 well-chosen 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₀ is the second kind modified Bessel function and K₁ = −.

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 non-zero 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 Ω₀), 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).

Fig. 4 Updated Lagrangian Formulation ([8]).

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 time-step 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₂ 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?

Long-term 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 salt-creep.

To answer these questions numerical models are created with geological, geometrical, and flowrates data as input. Those models are history-matched with continuous or one-off measurements as wellhead or cavern pressure and temperature and free volume of the cavern.

From a practical point of view, each cavern has a history-matched model, which allows short-term 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 history-match 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 Well-head pressure and temperature – available data (grey marks) – estimated by model (black marks) – flow rate (black line) is imposed in the model.

The history-matched 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 short-time 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³, 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⁻¹·m⁻¹, c _p  = 920 J·kg⁻¹·K⁻¹, ρ_rock = 2100 kg·m⁻³, while along the well λ = 5 W·K⁻¹·m⁻¹, c _p  = 850 J·kg⁻¹·K⁻¹, ρ_rock = 2100 kg·m⁻³.

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³, depth casing shoe at 1200 m. Initial lithostatic stress is 367.8 bar.

Young’s modulus of salt rock is 2 × 10⁵ 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⁻¹·m⁻¹, c _p  = 800 J·kg⁻¹·K⁻¹, ρ_rock = 2200 kg·m⁻³.

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 well-cavern 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⁻¹·m⁻¹, c _p  = 920 J·kg⁻¹·K⁻¹, ρ_rock = 2100 kg·m⁻³, while along the well λ = 5 W·K⁻¹·m⁻¹, c _p  = 850 J·kg⁻¹·K⁻¹, ρ_rock = 2100 kg·m⁻³.

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).

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₂ 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³ for hydrogen versus 0.789 kg/m³ 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 well-head pressure for methane (solid line), hydrogen (dotted line) and mix (dash-dotted line). Down: Variations of cavern temperature for methane (solid line), hydrogen (dotted line) and mix (dash-dotted line).

5.3 Other applications

5.3.1 Cavern conversion

The conversion from L-gas to H-gas is a large infrastructure project in the gas industry in Europe and concerns some storages for example in North Europe. The L-gas (with low calorific value) will be replaced by H-gas (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 L-gas and H-gas. 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³, 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⁻¹·m⁻¹, c _p  = 920 J·kg⁻¹·K⁻¹, ρ_rock = 2100 kg·m⁻³, while along the well λ = 5 W·K⁻¹·m⁻¹, c _p  = 850 J·kg⁻¹·K⁻¹, ρ_rock = 2100 kg·m⁻³. 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 L-gas withdrawal before the H-gas 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 H-gas (Fig. 10, point A); the first H-gas 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 H-gas is higher than that of L-gas at the same pressure (due to their different thermodynamic properties).

Fig. 10 H-gas mass fraction as function of the well-head 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 long-term behavior is considered, the picture is modified: in-situ 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 build-up 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 day-to-day basis to monitor, predict the behaviour of caverns, prevent operational issues and also to define commercial strategies in terms of storage capacity.

References

All Tables

All Figures

	Fig. 1 Simplified model considered.
In the text

	Fig. 2 The well-rock 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 Well-head 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 well-head pressure for methane (solid line), hydrogen (dotted line) and mix (dash-dotted line). Down: Variations of cavern temperature for methane (solid line), hydrogen (dotted line) and mix (dash-dotted line).
In the text

	Fig. 10 H-gas mass fraction as function of the well-head pressure during conversion.
In the text