# Long-term slope stability of abandoned mine lake - numerical modelling and risk assessment

Vincent Renaud, Marwan Al Heib, Roxane Muller, Jan Burda

► To cite this version:

Vincent Renaud, Marwan Al Heib, Roxane Muller, Jan Burda. Long-term slope stability of abandoned mine lake - numerical modelling and risk assessment. International Conference on Raw Materials and Circular Economy (RawMat 2021), Sep 2021, Athens, Greece. pp.art. 88, 10.3390/mater proc2021005088. ineris-03539344

HAL Id: ineris-03539344

https://ineris.hal.science/ineris-03539344

Submitted on 21 Jan 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 Attribution 4.0 International License
Proceeding Paper

# Long-Term Slope Stability of Abandoned Mine Lake—Numerical Modelling and Risk Assessment †

Vincent Renaud <fnref n="1" />, Marwan Al Heib <fnref n="1" />,*, Roxane Muller <fnref n="1" /> and Jan Burda <fnref n="2" />

**Abstract:** Almost all post exploitation open pit mines in the world are shaped as a final reservoir intended to be filled with water. In Europe, the creation of water lakes is the most common way of reclaiming post open pit mines. The safety and the security of mine lakes is one of the priorities of mine regions. One of the main hazards identified is the slope stability of lake banks. To develop a reliability methodology for assessing the long-term stability of flooded open pit mines, a large-scale numerical model of the lake was carried out and was applied on Lake Most, which is one of the largest mining lakes in Europe (Czech Republic). The large-scale numerical model was built, based on the site observations, large scale LiDAR data and geotechnical data. The results highlighted the reliability of the methodology to combine the geometric model with the geological model to create a large-scale numerical model, and to identify local and potentially instable zones.

**Keywords:** mining lake; slope stability; risk assessment; LiDAR; numerical modelling

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

check for
updates

**Citation:** Renaud, V.; Al Heib, M.; Muller, R.; Burda, J. Long-Term Slope Stability of Abandoned Mine Lake—Numerical Modelling and Risk Assessment. *Mater. Proc.* **2021**, *5*, 88. https://doi.org/10.3390/materproc2021005088

Academic Editor: Evangelos Tzamos

Published: 29 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction and Objective

Almost all post exploitation open pit mines in Europe and in the world are shaped as a final reservoir intended to be filled with water. These artificial lakes are currently (and in the future) dedicated to recreational purposes including energy installation purposes (SUMAD RFCS project). To ensure safe utilization of these localities by the public, it is necessary to assess their potential risk of instability. Slope instability is the major long-term hazard. The RFCS RAFF project aims to research issues related to pit lakes. The stability of pit lake slopes after flooding remains an area of uncertainty. Examples of geotechnical failures in slopes and banks of pit lakes are quite well documented, for example those at pit lake Pątnów ([1]), Zülpich Mitte and Concordia Lake near Nachterstedt ([2]). Back analysis is generally carried out to understand the cause of these phenomena and to suggest mitigation strategies.

The Lake Most site (Czech Republic) was selected, as a case study, to perform a large-scale stability analyses based on in-situ observation, a 3D geometric model and a large-scale numerical model using strength reduction method. Lake Most (located below the hill Hněvín) was formed on the site of the former royal town of Most, which had to give way to brown coal mining in the second half of the 20th century. Coal mining was terminated in August 1999. The lake was established in the open pit mine and the dumps. The flooding up began in October 2008 and was finished in September 2014 (surface = 309 ha, volume = 70 hm³; [3]).

## 2. Long-Term Slope Stability

To develop a reliability methodology for assessing the long-term stability of flooded open pit mines, a slope stability study of the lake area was carried out with a large-scale 3D

<footnote n="1">¹ Institut National de L'Environnement Industriel et des Risques (Ineris), 54042 Nancy, France; vincent.renaud@ineris.fr (V.R.); marwan.alheib@ineris.fr (M.A.H.); roxane.muller@ineris.fr (R.M.)</footnote>

<footnote n="2">² Výzkumný Ústav pro Hnědé Uhlí (VUHU), tř. Budovatelů 2830/3, 434 01 Most, Czech Republic</footnote>

<footnote>* Correspondence: burda@vuhu.cz</footnote>

<footnote>† Presented at International Conference on Raw Materials and Circular Economy, Athens, Greece, 5–9 September 2021.</footnote>
numerical model using Flac3D, utilizing an explicit finite difference formulation that can model complex soil and rock behaviors. The 3D model is based on the site observations, large LiDAR data, geological, hydraulic and geotechnical data.

## 2.1. Geometric and Geologic Data

A 3D numerical model is time consuming, however the study area must, at least, integrate the areas of ground movement already identified and must consider the heights of geometric anomalies (valleys or hills). The maximum depth of the lake is 75 m, the highest hill in the immediate vicinity of the lake has a height of 70 m and the boundaries of the 3D model were positioned at more than 6 times 70 + 75 m from the shores of Lake Most. Additionally, a large LiDAR campaign was carried out in August 2019. The final data point cloud was used to create the digital terrain model to build the 3D volumetric mesh.

To avoid too many meshes and therefore prohibitive calculation a finer mesh corresponding to 5 m was adopted based on 2D calculation for the area of interest whose horizontal perimeter includes all records and observed ground movements in the Most lake slopes. The elements at the boundaries are cubes of 50 m edge. The 3D geometric model thus designed has 11,826,069 elements. A total 14 non tabular geologic formations are identified, and faults are present. The geology of the site is based on thousands of boreholes (carried out between 1867 and 2018) that have been used to build DMT of the interfaces of the different formations. Since the spatial distribution of these boreholes is not homogeneous and does not always extend to the limits of the model, it was necessary to interpolate and extrapolate the position of these interfaces. The analysis of the numerous data of this campaign has helped to determine both the position of these interfaces as well as the values of the geomechanical parameters. Figure 1 represents the different geologic units used for building the 3D numerical model.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** 3D geological model of lake Most (11,826,069 elements); rectangular limit: area where the FoS is minimal.</figcaption></figure>

## 2.2. Geotechnical Data

The geotechnical data are based on 2 groups: the data relating to geological formations and the data characterizing the dump units ([3]). The first group concerns all the units described in Figure 2 except NS-TV1 and TV2 anthropogenic units. Among these geological units, the 3 least resistant units (quaternary gravel, JIL1 = plastic soft clays and PJIL = sandy clays: 52 kPa < σ<sub>c</sub> < 143 kPa) will play a role in the calculation of slope stability, especially in areas where these units are near the topographic surface. The second group of data is very rich in geotechnical information due to a recent measurement campaign (CPT). There
are 9538 CPT measures available, including pore pressure (u2), sleeve friction (fs), depth, cone resistance (qc), corrected cone resistance (qt), soil weight (γ), density (ρ), pore water pressure (u), total stress (σv) and Young's modulus (E). The last 6 parameters are classically calculated from the data recorded every 10 cm and parameters characterizing the used cone (net area ratio: a = 0.8 and penetrometer diameter B = 0.0357 m). As materials are assumed to have a Mohr-Coulomb elastoplastic constitutive model, the cohesion (C) and the friction angle (φ) were calculated from these data. The CPT data interpretation theory manual gives 3 empirical relationships ([4–6]) for assessing φ. The disadvantage is that these relationships do not give access to cohesion and that their validity is limited to sands or to fine-grained soils. Therefore, we preferred to use the system of 2 equations proposed by Motaghedi & Eslami ([7]) which gives access to C and φ after numerical resolution.

$$u_2 + \gamma B \tan \phi + (\sigma_v - u) N_q + \gamma B N_q \tan^2 \phi + C \frac{N_q - 1}{\tan \phi} = q_t \quad (1)$$

$$C = f_s - 0.000789(1 - \sin \phi)(\sigma_v - u) \tan \frac{2\phi}{3} \left( \frac{q_t - \frac{\sigma_v}{3}(1 + 2(1 - \sin \phi))}{\frac{\sigma_v - u}{3}(1 + 2(1 - \sin \phi))} \right)^{1.44} \quad (2)$$

$$\text{with the bearing capacity factor } N_q = \frac{(\sin \phi + 1)^2}{\cos \phi} \exp \left[ \left( 0.0061 \frac{q_t}{P_a} + 2.4 \right) \tan \phi \right] \quad (3)$$

<figure><img src="image_5.png" /><figcaption>Figure 2. (a) Cohesion and friction angle computed for profile P1_2. (b) Linear variation of cohesion and friction angle for all CPT profiles for NS-TV1 unit; minimum and maximum limits in red and mean value in cyan.</figcaption></figure>

The parameters of the CPT campaign were recorded for 23 exploitable penetrometric profiles whose depth of investigation varies between 10 m and 101 m. This depth was calibrated to the estimated thickness of the mining deposits. These 23 profiles have been statistically analyzed with an automated method based on the evolution of the coefficient of determination (between depth and cohesion). The position of the interface between different dump units is obtained by analyzing the cohesion profiles (Figure 2). This completes the determination of the geological model shown in Figure 1 with a homogeneous structured mesh. To gain accuracy in the areas of interest, this mesh will be replaced by a heterogeneous unstructured mesh. In addition, an assumption was made concerning a
contact layer. This layer was not detected from the CPT campaign measurements. This layer may play a main role on the slope stability and in situ observations.

Two calculation scenarios have been built from the average values and the minimum bounds of C, φ, ρ and E. A third scenario is possible by modeling the properties with the most representative distribution of the measurements. For the TV2 unit (properties not very sensitive to depth), we considered all the data attributed to this unit (2843 measurements) to find the best distribution. To account for the variation of the statistical parameters with the depth in the NS-TV1 unit, we gathered the measurements in 2 m depth. The measurements are numerous enough for this unit (population = 6629) to have enough individuals to analyze. The criterion for selecting the best distribution is a minimization of the square of difference of the theoretical and measured frequencies. As one of the aims of the RAFF project is to develop an operational methodology, we have chosen, for modeling, only the easily implementable distributions: normal distribution, lognormal and Birnbaum-Saunders. The results of this statistical analysis are reported in Table 1. Note that the lognormal distribution is the one that best represents variations in cohesions and Young's modulus. The same is true for the friction angle except for the TV2 unit which correlates better with a Birnbaum-Saunders distribution. The normal distribution is the best fit for spatial distribution of the density masses in the 23 CPT profiles. The probability density function of Birnbaum-Saunders distribution is:

$$pdf = \frac{\sqrt{\frac{x}{\beta}} + \sqrt{\frac{\beta}{x}}}{2\sqrt{2\pi}\gamma x} e^{-\frac{\frac{x}{\beta} + \frac{\beta}{x} - 2}{2\gamma^2}} \quad \text{with} \quad \beta = \frac{\mu}{3}(4 - \sqrt{3c^2 + 1}), \gamma = \sqrt{\frac{2(c^2 - 1 + \sqrt{3c^2 + 1})}{5 - c^2}} \quad \text{and} \quad c = \text{COV} = \frac{\sigma}{\mu} < \sqrt{5} \quad (4)$$

<table><caption>**Table 1.** Statistical properties of dump units. N, LN & BS dist: Normal, Log-normal and Birnbaum-Saunders distributions; *d*: depth.</caption><thead><tr><th>Dump Unit/Property</th><th>NS + TV1</th><th>TV2</th><th>Contact Layer</th></tr></thead><tbody><tr><td>cohesion: C (kPa)</td><td>mean: 4.59 d + 46.15<br>min: 3.5 d + 10<br>LN dist ∈ [0.6, 650]: μ<sub>LN</sub> = 0.0324 d + 4.034<br>σ<sub>LN</sub> = -0.0063 d + 0.598</td><td>mean: 247.4<br>min: 4 d - 60<br>LN dist ∈ [6, 1098]: μ<sub>LN</sub> = 5.31 & σ<sub>LN</sub> = 0.626</td><td>6.0</td></tr><tr><td>friction: φ (°)</td><td>mean: -0.323 d + 30.69<br>min: -0.13 d + 20<br>LN dist ∈ [7, 44]: μ<sub>LN</sub> = -0.241 d + 28.73<br>σ<sub>LN</sub> = -0.071 d + 6.074</td><td>mean: 22.7<br>min: 16.8<br>BS dist ∈ [8.2, 38.6]: β = 22.019 & γ = 0.244</td><td>6.0</td></tr><tr><td>Young's modulus: E (MPa)</td><td>mean: 3.67 d + 18.51<br>min: 2.96 d + 5<br>LN dist ∈ [1, 354]: μ<sub>LN</sub> = 0.037 d + 3.577<br>σ<sub>LN</sub> = -0.004 d + 0.377</td><td>mean: 193.9<br>min: 3.6 d - 52<br>LN dist ∈ [29.6, 636.5]: μ<sub>LN</sub> = 3.855 d + 32.046<br>σ<sub>LN</sub> = 0.6195 d + 28.295</td><td>70</td></tr><tr><td>saturated density: ρ (t/m³)</td><td>mean: 0.0124 d + 1.739 (d ≤ 12.5 m)<br>0.0035 d + 1.869 (d > 12.5 m)<br>min: 0.0035 d + 2.0<br>N dist ∈ [1.3, 2.2]: μ = 0.0053 d + 1.813<br>σ = -0.0008 d + 0.0919</td><td>mean: 2.023<br>min: 0.0017 d + 2.079<br>N dist ∈ [1.55, 2.285]: μ = 2.023 & σ = 0.098</td><td>2.0</td></tr></tbody></table>

## 2.3. Hydraulic Data

Typically, the position of the water table is poorly known and it is often determined by a few points where the pore pressure is known, by the hydromechanical properties of the studied layers (permeability, porosity, Biot coefficient, dynamic viscosity) and by boundary conditions (null fluid flux, prescribed pore pressure field, inflow, outflow, leakage). Fortunately, the Lake Most site is very well instrumented: the piezometric heights of the water table have been recorded for 93 wells since 2014. The water table DMT has been built from 93 piezometric water levels and lake surface. For slope stability calculations, the position of the water table has been used to calculate the pore pressure and effective stress fields.
## 3. Results

The use of the strength reduction method produces one global minimum stability state per simulation. This method is applied with the Mohr-Coulomb failure criterion ([8]) by progressively reducing the shear strength of the material to bring the slope to a state of limiting equilibrium. However, along a complex slope profile (as the profile of the NNW-SSE cross section: Figure 3), it is interesting to be able to compute multiple minimum states. Instead of excluding different regions of the slope when performing the strength reduction calculation, we have used the ability of Flac3D ([9]) to compute multiple local stability surfaces in a single simulation.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** FoS contours for the 4629 m length NNW-SSE cross section with distributions of C, E, $\phi$ and $\rho$ in the dump units (a) without contact layer. (b) with contact layer.</figcaption></figure>

Six calculations were made to estimate the safety factor of the Most site in its current (short-term) situation. The 6 calculations were established from the scenarios of geomechanical properties of dump units (mean values, minimum bounds or statistical distribution) and the presence or not of the contact layer at the bottom of the dump bodies. The contact layer is characterized by very low geomechanical parameters ($C = 6$ kPa and $\phi = 6$). Some results of these scenarios are shown in Figure 3. First, the 3 scenarios without contact layer produce the same isovalues of FoS because these 6 scenarios differ only in the properties of the dump units and because dump units are not the weakest geologic units. On the other hand, the 3 scenarios including the contact layer change the stability of all dump units (NS-TV1, TV2 and contact layer). This is explained by the very low properties of the contact layer. Areas with FoS of 2.75 (NNW) and 2.9 (SSE) without contact layer have a safety factor of between 1.14 and 1.53 at NNW and between 1.5 and 1.6 at SSE when contact layer is present. The 3D calculations give results compatible with this 2D cross section: FoS = 2.2 and 1.38 without or with contact layer respectively. The global 3D FoS is located on the north bank of Lake Most (Figure 2) at the location where the majority of slope failure stabilization operations took place in the past. In those zones earth and stabilization work were carried out to insure long-term stability.

## 4. Conclusions

In order to establish a reliability methodology for assessing the long-term stability of flooded open pit mines, a large-scale numerical model of Lake Most was carried out with Flac3D. The model integrated the complex geology of the mine and the dumps as well as the surface of the water table interpolated from 93 piezometric levels. The results of the 2D and 3D numerical modelling were analyzed as large scale by calculating global and local safety factors. The results highlighted the reliability of the methodology to combine the geometric, geological and hydraulic models to create a large-scale numerical model, and to identify local and potentially instable zones. The hypothesis of the presence of a very weak contact layer (at the bottom of dump bodies) is therefore a strong hypothesis, capable of questioning the stability of the slopes of the site (Most Lake). It should be noted that the contact layer was not detected from the CPT campaign measurements. But neither the absence nor the presence of the contact layer can be confirmed because only 2 profiles are
deeper than the dump units on an area to be investigated of more than 8 km². The most realistic hypothesis is probably to consider a partial presence of the contact layer. These new calculations could show whether the reality lies between the minorant (with contact layer) and majorant (without contact layer) scenarios. In this case, the Lake Most site could then be concluded to be stable in the short to medium term.

**Author Contributions:** Methodology, V.R. and M.A.H.; investigation, V.R. and J.B.; data acquisition, R.M. and V.R.; writing—V.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has received funding from the European Union's Research Fund for Coal and Steel under the project "Risk Assessment of final Pits during Flooding (RAFF)" No. 847299. Financial assistance by the European Commission is much appreciated.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Szafranski, C.; Stachowski, P.; Kozaczyk, P. Influence of meteorological conditions on water balance of soils created from post-mining grounds. In Proceedings of the 5th National Polish Scientific Conference on Complex and Detailed Problems of Environmental Engineering, Warsaw, Poland, 24–27 May 2001.

2. Vinzelberg, G.; Dahmen, D. The Inden Residual Lake and the Federal Autobahn A44n—Geotechnical Requirements to Be Met by Large-Scale Structures in the Rhenish Lignite Mining Area. In *Proceedings of the 12th International Symposium Continuous Surface Mining—Aachen 2014*; Springer: Cham, Switzerland, 2014; pp. 91–100. [CrossRef]

3. Žižka, L.; Burda, J. Post Exploitation Lakes. In *Zpravodaj Hnědé Uhlí*; Brown Coal Research Institute: Budovatelů, Czech Republic, 2021; pp. 11–46.

4. Robertson, P.K.; Campanella, R.G. Interpretation of cone penetration tests—Part I (sand). *Can. Geotech. J.* **1983**, *20*, 718–733. [CrossRef]

5. Kulhawy, F.H.; Mayne, P.W. *Manual on Estimating Soil Properties for Foundation Design*; Report EL-6800; Electric Power Research Institute, EPRI: Washington, DC, USA, 1990.

6. Mayne, P.W. In situ test calibrations for evaluating soil parameters. In Proceedings of the Characterization and Engineering Properties of Natural Soils II, Singapore, 29 November–1 December 2006.

7. Motaghedi, H.; Eslami, A. Analytical Approach for determination of soil shear strength parameters from CPT and CPTu data. *Arab. J. Sci. Eng.* **2014**, *39*, 4363–4376. [CrossRef]

8. Zienkiewicz, O.C.; Humpheson, C.; Lewis, R.W. Associated and non-associated visco-plasticity and plasticity in soil mechanics. *Géotechnique* **1975**, *25*, 671–689. [CrossRef]

9. Itasca. *FLAC3D (Fast Lagrangian Analysis of Continua in Three Dimensions) Version 7.0, User's Guide and Command Reference*; Itasca Consulting Group Inc.: Minneapolis, MI, USA, 2019.