# Numerical modelling of leached cavern fields using mixed BEM-FEM Method

José Eduardo Tavares Quintanilha de Menezes, Duc Nguyen-Minh

► To cite this version:

José Eduardo Tavares Quintanilha de Menezes, Duc Nguyen-Minh. Numerical modelling of leached cavern fields using mixed BEM-FEM Method. 4th Conference on the Mechanical Behaviour of Salt, 1996, Montréal, Canada. pp.417-426. hal-00116431

HAL Id: hal-00116431

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

Submitted on 28 Oct 2022

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
NUMERICAL MODELLING OF LEACHED
CAVERN FIELDS USING MIXED
BEM-FEM METHOD

J. Eduardo Q. de Menezes and Duc Nguyen-Minh

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

ABSTRACT

More than ten years of surface subsidence observation on Tersanne and Etrez storage cavity fields in France, is subject to numerical analysis, using a mixed boundary and finite element approach. The consideration of large domains to discretize and non linear behaviours to take into account is simplified as the bedded rock mass is modelled in two horizontal strata : an elastic clay cover and a viscous salt bed, the first one discretized in boundary elements and the second in finite elements. Although the problem is typically three-dimensional, an axisymmetric model for a single cavity is used throughout the numerical analysis, and a superposition method, previously validated, is applied to a 14-cavity field. There is a reasonable agreement between in situ data and numerical results.

INTRODUCTION

The search for a good characterization of long term behavior of leached rock salt cavities is justified by the need to minimize eventual damages resulting from subsidence on the surface.

In two locations in France, Tersanne and Etrez, subsidence control has been
performed during more than ten years. Simulation by numerical methods allows to
perform some predictions on future evolution, and also to test specific numerical
procedures. The issue is to obtain a good accuracy in the far field region, e.g. the
surface subsidence. The data from Tersanne site have been already subject to
numerical analyses (Nguyen-Minh et al, 1993, 1994), where a two step finite element
procedure was used, by taking advantage of the two-layered medium and the elastic
response of the cover-layer. Herein, such procedure in not needed as the far field
region can be efficiently discretized with boundary elements and specific infinite
elements (Bettess, 1992). A special procedure is used to speed up numerical
calculations (Menezes & Nguyen-Minh, 1993).
TERSANNE FIELD DATA AND CAVITY ARRAYS

At Tersanne and Etrez, in France, subsidence has been measured periodically
for more than ten years (Durup, 1990, 1991). The distribution of cavities is shown in
Figure 1. Etrez Cavities are more regularly spaced than in Tersanne and distance
between cavities is larger.

A quasi axisymmetric subsidence pattern can be observed at the surface at Tersanne (Figure 2). The volume loss at the surface to the total caverns volume loss shows a 60% ratio. As rock salt is believed to behave as an incompressible material, one can suspect a possible dilatancy inside the rock mass, which could affect its tightness.

<figure><img src="image_2.png" /><figcaption>Figure 1 - Pattern Distribution of Storage Cavities in Tersanne and Etrez</figcaption></figure>

The observation zone (Figure 3) covers an area of 15 km², where cavities are leached 1500 m depth, in a thick 650 m bedded salt layer and located 150 m beneath the clay cover. Beyond 2000 m depth there is a rigid sandstone base. There are fourteen cavities spaced 600 m in a hexagonal arrangement, with an equivalent volume of 268 000 m³. We will assume that all the cavities are leached at the same time and that they have the regular distribution as seen in Figure 3. Comparing with the real pattern distribution of cavities (Figure 1), this hypothesis is not far from reality in Etrez case.
<figure><img src="image_3.png" /></figure>

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

Figure 2 - Subsidence data from Durup (1991): (a) Subsidence Sinkhole; (b) Comparison between Subsidence and Cavity Volume Losses

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

Figure 3 - The Cavity Array: Subsidence Observation Zone and Idealized Cavity Field
# NUMERICAL MODEL DESCRIPTION

## Two-layered mixed rock mass model

The numerical analysis of the subsidence phenomena of a layered unbounded field problem is easier to accomplish if one uses a mixed method. The infinite domain is divided in two or more layers, with linear and elastoviscoplastic mechanical behavior. Boundary elements provide the conditions to represent a large surface line where one needs to compute displacements over quite a large distance from the axis of symmetry. Finite elements are easier to deal with non linear materials (salt). The discretization in finite elements makes easier the computation of stresses inside the salt layer needed in every time step.

The boundary element region is treated as a larger finite element. For the FE region the finite element matrices are $KU=F$, where $K$ is the stiffness matrix, $F$ the equivalent nodal forces and $U$ represents the displacements. The boundary elements matrices are different: $HU=GP$, where $U$ also represents the displacements and $P$ the surface traction's (Figure 4.a). To combine the two equations one must reduce the BE expression to a finite element form: $G^{-1}HU=P$. Then the traction values must be converted by means of the distribution matrix $M$, that is $F=MP$. Finally the BE equation will be $M(G^{-1}H)U=MP=F'$ or $K'U=F'$. The two matrices of BE and FE regions can now be assembled to form the global stiffness matrix (Brebbia, 1984; Pande, 1990).

A layered rock mass model was chosen with an high behaviour contrast: the cover layer, consisting of non viscous material, and the viscous rock salt layer. There were also different primary states of stress in each layer: isotropic initial state in the salt strata and anisotropic above in the clay strata. Figure 4.b shows the discretization of a single isolated cavity.

To simulate the infinite conditions, characteristic to this kind of problems, the border elements are replaced by mapped infinite elements (Figure 4.b). These elements have particular shape functions in the infinite direction. They provide good results in the far field region if they remain inside a non plastic yielding zone (Zienkiewicz, 1981; Bettess, 1992).

## Constitutive models for rock salt

Two alternative models are used herein: the Norton-Hoff (N-H), and Lemaitre-Menzel-Schreiner (L-M-S) laws. N-H law considers the existence of a secondary or steady state creep (Langer et al, 1984), while L-M-S law supposes a continuously strain rate decreasing creep (Vouille et al, 19843). These creep laws can be described by the same formula:

$$ \begin{gathered} \dot{\varepsilon} = \dot{\varepsilon}_{el} + \dot{\varepsilon}_{vp} \qquad \dot{\varepsilon}_{el} = \text{elastic strain rate} \qquad \varepsilon_{vp} = A \left(\frac{\sigma}{\sigma_o}\right)^n t \\ \dot{\varepsilon}_{vp} = \text{viscoplastic strain rate} \end{gathered} $$
Tersanne (N-H)

$$A = 5.68 \times 10^{-8} \text{ day}^{-\alpha}$$

$$n = 3$$

$$\alpha = 1$$

Tersanne

$$A = 1.795 \times 10^{-6} \text{ day}^{-\alpha}$$

$$n = 3.6$$

$$\alpha = 0.5$$

Etrez

$$A = 1.70 \times 10^{-3} \text{ day}^{-\alpha}$$

$$n = 2.98$$

$$\alpha = 0.35$$

Salt layers: $E = 25000 \text{ MPa}$

$$\nu = 0.25$$

$$\sigma_v = 1 \text{ MPa}$$

$$\gamma = 22 \text{ kNm}^{-3}$$

(reference stress)

Elastic layers: $E = 6000 \text{ MPa}$

$$\nu = 0.4$$

$$K_o = \frac{\nu}{1-\nu} = \frac{\sigma_h}{\sigma_v}$$

$$\gamma = 23 \text{ kNm}^{-3}$$

These parameters were obtained by both laboratory and *in situ* observations (Fauveau, Le Bitoux, 1986). The equivalent N-H model for Tersanne salt was determined on the basis of a same volume loss for a single cavity at 1500 m depth, with an inner pressure of 16.5 MPa, during 10 years (Nguyen-Minh et al, 1993). For the moment, such procedure was not followed with Etrez case.

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

<figure><img src="image_7.png" /><figcaption>Figure 4 - a) Layered Unbounded Problem Simulation; b) Infinite Elements</figcaption></figure>

## SINGLE CAVITY PROBLEM AND SUPERPOSITION PROCEDURE

Figure 5 shows evolution of a single cavity volume loss and the subsequent surface subsidence volume. These were computed for the two constitutive laws with Tersanne parameters. After a period of 50 years, numerical results with N-H law are approximately double from those of L-M-S law. However, such discrepancy should be expected.

Initially, the same 60% ratio surface/cavity volume loss is obtained, but subsidence volume variation becomes progressively larger than cavern volume loss probably due to numerical errors as no geometry change is considered here.
<figure><img src="image_8.png" /><figcaption>a) N-H law</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 5 - Cavity and Subsidence Volume Variations</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6 - Evolution of Subsidence Surface above a Single Cavity</figcaption></figure>
For the 14-cavities storage array (Figure 3), a plain superposition of the single cavity solution (Figure 6) is performed. Theoretically, this should not be valid as non-linear deformations are involved. But for an infinite array, and for the spacing under concern (600 m), the superposition method has been checked to be accurate when compared to a direct calculation of a unit axisymmetric cell. Thus, the method can be applied with confidence to any array of finite number of cavities (Nguyen-Minh et al., 1993).

The resulting subsidence for the 14 cavity array was determined along its two axes of symmetry. The difference between these two subsidence values is very small, thus showing that on the surface the subsidence is approximately axisymmetric with a radius of 3000 m. The subsidence values along the vertical axis of symmetry of the 14-group cavity are represented on Figure 7. After a period of 50 years there would be a subsidence of more than 22 cm in the centre of the array. However, by N-H law, such provision should be rather pessimistic.

There is a good agreement between the subsidence values corresponding to 5 and 10 years interval time. The vertical lines in Figure 8 are proportional to the subsidence in situ for periods of 5 and 10 years. The figures correspond to a superposition of surface deformations from the real arrays (Figure 1) and the model array (Figure 3). It must be emphasised that the numerical model used herein admits that all 14 cavities have the same age.

Subsidence at Etrez is extremely small. As Figure 1.b shows, cavities are more largely spaced than in Tersanne which leads to a bigger zone of influence, but the main reason is due to the different viscous nature of Etrez salt. Numerical calculations for the Etrez case are more intricate than in Tersanne. There are four layers, which results in two BE domains, and besides, cavities are closer to the interface between salt and elastic (clay) layer. This particular situation has not been analysed in detail until now.

## CONCLUSIONS

Subsidence and volume losses obtained depend on the choice of constitutive models. With Norton-Hoff law larger values (50% plus) are obtained after 50 years, as this model describes a stationary creep, while the Lemaitre-Menzel-Schreiner is a continuous decaying creep rate model.

The mixed method used to analyse the subsidence problem exposed herein appears to be quite efficient. Compared to a plain Finite Element approach, the discretization process is much simpler, as large elastic domains are represented by boundary elements, and the results are obtained more rapidly, with a similar or better accuracy.

For large periods of time, a regular change in structure geometry must be done in order to take account of finite displacements, leading to a final closure of the cavity.
<figure><img src="image_11.png" /><figcaption>Figure 7 - 14 Cavity Field Surface Subsidence (N-H law)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8 - Comparison between Subsidence Measured In Situ and Computed Values in Tersanne (a) and Etrez (b)</figcaption></figure>
ACKNOWLEDGEMENTS

This work would not have been possible without the valuable contribution from
Mr. J.G. Durup and Mr. C. Rolin from GDF - France. The first author would like to
express his acknowledgement to JNICT - Portugal, for financial support.

REFERENCES

Bettess, P., (1992), "Infinite Elements", Penshaw Press.

Brebbia, C.A., Telles, J.C.F., and Wrobel, L.C., (1984), "Boundary Element Techniques", Springer-Verlag.

Durup, J.G., (1990), "Surface subsidence measurements on Tersanne cavern field", Solution Mining Research Inst., Paris, October 14-19.

Durup, J.G., (1991), "Relationship between subsidence and cavern convergence at Tersanne (France)", Solution Mining Research Inst., Atlanta, April 28-30.

Fauveau, M. and P. Le Bitoux, (1986), "Progrès récents de la connaissance du comportement mécanique des cavités salines", Assoc.Tecnh. Ind. du Gaz en France, Congrès 1986.

Langer, M., (1981), "The rheological behavior of rock salt", Mechanical Behavior of Salt I, Trans.Tech.Pub., 1984, Pag. 201-240.

Nguyen-Minh, D.; Braham, S., and Durup, J.G., (1993), "Modelling Subsidence of the Tersanne underground gas storage field", 3rd Conf. Mechanical Behavior of Salt, Palaiseau, France.

Nguyen-Minh, D.; Maitournam, H.; Braham, S., and Durup, J.G., (1993), "A tentative interpretation of long-term surface subsidence measurements over a solution mined cavern field", 7th Symposium on Salt, Vol.1, pag. 441-449, Elsevier.

Menezes, J.E.Q., and D. Nguyen-Minh, (1993), "Incompressible numerical modelling for long term convergence evaluation of underground works in salt", 3rd Conf. Mechanical Behavior of Salt, Palaiseau, France.

Pande, G.N., Beer, G., and Williams, J.R., (1990), "Numerical Methods in Rock Mechanics", J.Wiley & Sons.

Vouille, G., S.M. Tijani and B. Hugout, (1983), "Le sel gemme en tant que liquide visqueux", Proceedings 5th Int. Cong. Rock Mechanics, pp. D241-D246, Melbourne, Australia.

Zienkiewicz, O.C., P. Bettess, T.C. Chiam and C. Emson, 1981, "Numerical methods for unbounded field problems and a new infinite element formulation", ASME, AMD, 46, Pag.115-148, New York.