Long-Term Slope Stability of Abandoned Mine Lake—Numerical Modelling and Risk Assessment "2279

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.


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 largescale 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 ; [3]).

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

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.

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

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.

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

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 < σ c < 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 (q t ), 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][5][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.
with the bearing capacity factor N q = Mater. Proc. 2021, 5, 88 3 of 6 PJIL = sandy clays: 52 kPa < σc < 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][5][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. ( ) 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   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:

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.

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

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. 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 ϕ = 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. 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 φ = 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.

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