The Challenges of Open-Pit Mining in the Vicinity of the Salt Dome (Bełchatów Lignite Deposit, Poland)

Abstract

The extraction of the Bełchatów lignite deposit located in the vicinity of the Dębina salt dome requires careful planning that considers the influence of mining projects on the slope and salt dome stability conditions. The instability problem is directly related to horizontal and vertical displacement, as well as the complex geological and mining conditions. These conditions are very unique with regard to the co-occurrence of the salt dome and lignite deposits in the same area, as well as the large scale of the pit wall slope. Thus, predicting rock mass behavior and ensuring the safety of mining operations are important issues. The presented analysis focused on the influence of long-term lignite extraction on the western pit wall slope of the Bełchatów field and the salt dome’s stability conditions. This study offers a comprehensive approach to a complex geotechnical problem defined by large-scale, complex geometry, and geological conditions. The rock mass behavior and stress conditions are simulated in numerical modelling. The results of the presented analysis will be useful not only for present mining activities but also for future developments related to post-mining and recultivation plans.

1. Introduction

As the global economy is transforming, pushing it towards renewable and green energy, the world’s demand for energy is also increasing at its fastest level since 2010 [1]. While renewable energy sources and gas will satisfy a portion of this demand, coal, as the largest fuel-source for global industrial energy, is not predicted to be phased out for many decades, partially due to the dependency of many processes on coal as fuel, partially due to its rising popularity in developing countries, which use it to fuel their growth, and partially due to coal representing too large of a portion of the energy market to be substituted by another source completely [2]. Therefore, issues such as slope stability still need to be investigated and solved to improve the safety of mines and the industry as a whole.

Long-term stability plays an important role in the design of slope engineering and rock mass extraction in open-pit mines. The slope stability in open-pit mines is an extremely complex issue that requires a special approach. Due to the variety of factors that can affect the safety of open-pit mining operations, proper maintenance with proper monitoring and prediction tools is necessary. Providing effective lignite extraction and ensuring the safety of mining operations are key factors in the mining of the Bełchatów deposit. Another important issue is the prediction of slope stability and the behavior of the rock masses affected by mining activity. In the Bełchatów lignite deposit (BLD), the geological and mining conditions are difficult due to the complex geological structure of the BLD and the presence of the Dębina salt dome (DSD). This complex geological structure distinguishes the Bełchatów deposit from other lignite open-pit mines [3,4,5]. Moreover, the scale of mining operations in the Bełchatów field is very large, with lignite extraction accomplished over 300 m high pit wall slopes. These factors are the reasons for the application of geotechnical tools and numerical modeling in the monitoring and prediction of rock mass behavior and ensuring the safety of mining operations. The geological conditions and their influence on slope stability have been described in a few papers. Arikan et al. [3] investigated the geological and structural features of the rock units that govern the stability of slopes. Kayabasia and Gokceoglu [4] presented the slope instability problems in the Can (Canakkale—northwestern Turkey) coal mine characterized by difficult geological conditions. The application of numerical modelling indicated that the weak contacts between the lithological units are important factor underlying this instability. The paper by Zevgolisa et al. [5] provides a review of the slope failure incidents in lignite mines in Greece. Here, the authors focused on the role of geological factors in the slope failure mechanism.

The issue of slope stability in open-pit lignite mining, with the aid of numerical modeling techniques, was investigated by several authors. Numerical modeling methods such as the finite difference method (FDM), the finite element method (FEM), and the limited equilibrium method (LEM) provide the opportunity to conduct a comprehensive study of slope stability, its failure mechanisms, and landslide prediction. Moreover, numerical analyses can be applied to both two-dimensional (2D) and three-dimensional (3D) approaches to study the scale and complexity of simulated problems. Ozbay and Cabalar [6] described two large landslides that took place in the Çöllolar the lignite mine (Turkey). A series of 2D numerical analyses in geological cross-sections were carried out using FEM and LEM methods to evaluate the factors that cause landslides, indicating that the impact of the water table and the weak layers trigger landslides. Deliveris et al. [7] studied the slope stability of an open pit lignite mine in Greece using 2D numerical methods (FEM and FDM) by means of three popular software packages (PLAXIS 2D, FLAC 2D, and Phase2), whose performance was compared. All three programs demonstrated good agreement in the determination of safety factors and the kinetics of the collapse mechanisms. However, the authors stated that FLAC obtained more realistic results and was the most consistent in terms of the agreement between the localization of the maximum shear strains and the plasticity indicators. The paper by Lem et al. [8] utilized FEM and LEM analysis methods to investigate the stability of slopes with soft (weak material) bands, a postquake slope, and rock slopes in 2D. The authors found that results from the FEM and LEM analyses differed significantly in their failure mechanisms. The results obtained using the FEM analysis were more realistic than those from LEM because they did not require the failure surface to be pre-determined. The same methods were applied by Vanneschi et al. [9] to two-dimensional stability analyses of a lignite open-pit mine affected by landslides in North Bohemia (Czech Republic). This work notes that the complementary use of both approaches is recommended for routine checks of model responses and the interpretation of the associated results. Tutluoglu et al. [10] described the results of a 2D and 3D slope stability analysis in the FDM code of the Elbistan–Collolar lignite mine (Turkey). The aim of this study was to compare 2D and 3D modeling techniques. The authors observed that 3D methods can properly simulate the confining effect of 3D slope geometry and field stresses. The 3D method simulated the geometry of the slide mass more accurately. Moreover, the 3D FDM analysis provided the correct solution for the back analysis of slope failures and a better solution for the shear strength parameters of the different slope materials. Wang et al. [11] analyzed the stability of the inner dump slopes in the Shengli # Surface Coal mine located in China. They developed a 2D numerical model based on field observations and a borehole core quality analysis. They performed calculations with FDM and LEM and re-calibrated 2D models to a 3D case. The rock slope deformation stability for the Antaibao open pit coal mine (China) was analyzed by He et al. [12] using the FDM and FLAC 3D software. An optimal excavated scheme with steeper slope angle was successfully simulated and implemented in this analysis, taking into account the 3D effects of the nonlinear large deformation of the slope. Pinnaduwa et al. [13] incorporated field deformation monitoring data with 3D numerical modelling (with 3DEC) in order to predict the rock mass deformation in an open pit mine. The study conducted by Havaej et al. [14] offers an excellent practical example of calibrating a 3D numerical model (using 3DEC code) of a slate quarry. The authors calibrated their numerical model with the data obtained from the laser scanning of a portion of the mine.

Most of aforementioned investigations applied 2D numerical modeling to slope stability and analyzed this problem in chosen cross-sections. The 3D methods offer a complex approach for studying geological and mining conditions, as well as stability issues. Another factor that increases the complexity of the case presented in this paper is the presence of the salt dome and its influence on rock mass behavior. The numerical modeling of such a unique case and the challenges associated with the related mining operations have not been not discussed in literature.

The analysis presented in this paper focuses on the western pit wall slope of the Bełchatów field (WSBF) located in the vicinity of the DSD. The co-occurrence of rock salt, lignite, limestones, and other rocks affected by tectonic, sedimentation, and halokinetic processes make the geological structure of the BLD unique. However, these complex conditions produce new challenges for mining and geotechnics. This paper presents a comprehensive approach to solving the compound geotechnical and mining problems defined by their large scale, as well as their difficult geological and mining conditions, their complex geometry, and their rock mass’s mechanical behavior.

The present analysis aims to predict the impact of long-term lignite extraction on pit wall slope and salt dome stability. Numerical modeling and the FLAC3D v. 6.0 software were used to simulate rock mass behavior under the influence of mining projects. The geometry of the mining operation’s front-line, its changeability over time, and its position in relation to the salt dome, as well as the volume of its extracted rocks, were included in the numerical simulations. The period between January 2010 and October 2018 was analyzed. Due to a long timeframe of the analysis, the results are presented as two-years periods.

In this paper, the rock mass behavior in the open-pit mine adjacent to the salt dome is projected by its vertical and horizontal displacements and stress conditions. Moreover, to reflect the rock salt mass’s behavior in detail, the stress measurements from the AWID system installed in the borehole inside the DSD are included in the analysis. The results of the numerical modeling presented in this paper comply with the measurements of the subsurface deformations from the inclinometer system and the surface deformations from the geodetic network, which are the parts of the open-pit mine monitoring system. This problem has been broadly described in other works [15,16]; thus, it is not discussed in this paper.

2. Outline of the Studied Area’s Geology

The BLD is located in central Poland and is situated in the deepest part of the Kleszczów graben, which is an E–W trending segmented and subsided tectonic structure limited by vertical faults. The orientation of these faults complies with the direction of the regional faults. Thus, in the eastern part of the Kleszczów graben, the faults trend WNW–ESE, and in the western part, they trend WSW–ENE (Figure 1A). The main bounding faults are USB no. 1 (southern bounding fault) (in the south), UNB no. 2 (northern bounding fault), and UNB no. 3 (in the North) (Figure 1B). There are also strike-slip faults that cross the Kleszczów graben: the Faustynów faults, the Dębina fault, and the Chabielice fault [17,18].

The BLD is divided by the DSD into two mining fields (Figure 2): the Bełchatów field in the east and the Szczerców field in the west [17,19,20,21]. The DSD is one of eleven diapirs piercing the Cenozoic surface known as the Polish part of the Zechstein basin (Middle to Late Permian). The DSD is situated in the central part of the Kleszczów graben (Figure 1A), which is the easternmost part of the Neogene graben system formed in the Late Alpine phases (Early Miocene).

In addition to the complex tectonics related to regional trends, the processes connected to the formation of the DSD have resulted in the uplifting and brecciation of the surrounding sediments [22,23]. The lignite seams are inclined at an angle of 25–85° and folded in the area adjacent to the DSD [21]. Therefore, the study area is characterized by a complex geological structure related to the polyphasic tectonic activity and dynamics of salt structures.

The lithostratigraphic profile of the Bełchatów field in the vicinity of the DSD contains (from the bottom to the top) Jurassic and Cretaceous marls, sandstones, limestones, mudstone, a thick Miocene sequence, and Quaternary sands, loams, and clays (Figure 2). The Miocene sediments are subdivided into three units: an undercoal formation, a coal formation, and a clay–coal complex. The undercoal formation is mainly loams, sands, and clays. The coal formation in the northern part of the Bełchatów field is represented by several seams 5–30 m thick; however, in the southern part, there are up to two seams that are 80–100 m thick. The maximum thickness of the coal complex reaches 250 m in the eastern part of the Bełchatów field [24]. The clay–coal formation includes clays, sands, and one seam up to 5 m thick.

The structure of the DSD was initially recognized via drilling and geophysical methods [17,19,21]. The salt mirror was found at a depth of 170–200 m under the surface (Figure 2). The horizontal shape of the dome is an irregular ellipse elongated in the NW–SE direction. The maximum dimensions of the DSD axes reach 960 m (SW–NE) and 1040 m (NW–SE), with an area of about 0.5 km². The top of the salt dome is elevated 60–100 m above the top of Mesozoic sediments in the northern and southern flanks of the graben. However, at deeper levels, it extends parallel to the graben elongation and dips gently to the East [18,25].

The salt series consist mainly of Zechstein sediments (the first cycle/cyclothem PZ1) represented by medium, coarse, and mixed grained white rock and grey rock salt, with intercalations of fine-grained potassium-magnesium salt and rock salt with clay impurities. The rock salt beds are covered by a gypsum-anhydrite-clay cap. The thickness of the cap ranges from 30 to 118 m and is a residuum of the dissolution of the rock salt layers [17,21,25,28].

3. Methods and Conditions of the Stability Analysis

The stability analysis of the DSD was performed for the period 2011–2018 and based on numerical modeling.

3.1. Methods

Stability analysis is an effective method to predict the behavior of the rock mass under the given geological and mining conditions. The use of the numerical modeling in stability analysis enables one to describe the behavior of a rock mass through mathematical models, while the physical processes are represented by differential equations. The relevant mechanisms and constitutive laws with the associated variables and parameters are applied during the numerical modeling. The stability conditions of the DSD were evaluated with the use of the FLAC3D v. 6.00 software (Itasca Consulting Group) and the finite difference method (FDM), which is a valuable tool for solving geomechanical and geotechnical problems. In order to perform a stability analysis, a numerical model of the study area was built. The geological structure in the numerical model was projected based on 28 cross-sections covering the studied area (Figure 3). The cross-sections included lithostratigraphy and the main tectonic features, such as regional discontinuities, faults, folds, and brecciations (Figure 4). All these elements were projected into the 3D numerical model and reflected by the mechanical parameters. The numerical model is a cuboid with a length of 4000, a width of 2000, and a height of 2600 m (Figure 5). The model extends vertically from the ground surface (approx. +200 m.a.s.l.) to 2400 m below sea level (b.s.l.), which is the base of the salt dome. The boundaries of the numerical model were established based on previous experience and calculations performed in the 2D models. The bottom of the numerical model was fully fixed, and the left and back sides of the model were constrained by vertical rollers.

The mesh contained about 5 million hexahedral elements. To ensure that the complex geological structure of the area was accurately represented, the element size ranged from 5 to 80 m depending on the depth. Moreover, the geometry of the mining operations’ front line and its changeability based on time and its position relative to the salt dome, as well as the volume of the extracted rocks, were included in the numerical model.

3.2. Lignite Extraction Progress in 2011–2018

The progress of lignite extraction in the WSBF was analyzed to examine the impact of mining activity on the DSD. The starting point for the analysis presented in this paper is January 2010. At that time, a pit slope was formed in five benches +198/+176/+148/+120/+109 m.a.s.l. The elevation of the pit floor was +85 m.a.s.l. and reached +63 m.a.sl. at a distance of 1200 m from the pit toe (Figure 6A). The future pit slope design involved the formation of five more benches to the end of 2018 (+88/+74/+60/+41/+22 m.a.sl.). Another five were planned to be formed after 2018 (+2/−12/−30/−50/−65 m.a.s.l.). The designed benches are marked in blue and orange in Figure 3, while the benches that had been formed before 2010 are marked black with red captions.

Due to the long timeframe of the analysis, the geometrical changes are presented as two-year periods from October 2010 to October 2018. Between January and October 2010, lignite extraction was carried out mostly in the northern part (corner) of the Bełchatów pit slope located on the north-eastern side of the DSD. Benches +176 and +148 m.a.s.l. were formed on the north of the DSD in a specific elliptical shape outlined in the horizontal projection of the extraction area. Moreover, the northern part of bench +109 m.a.s.l. (specific triangular shape in the horizontal projection) was deepened to an elevation of +97 m.a.s.l. The first designed bench of +88 m.a.s.l. was formed, and an elevation of +74 m.a.s.l. was reached at the pit toe. The elevation of the pit floor was still +85 m.a.s.l., but reached +63 m.a.s.l. at a distance of 1100 m from the pit toe and was 100 m closer to the pit wall slope than in January (Figure 6B).

In the following two years (between October 2010 and October 2012), the pit floor was deepened. An elevation of +63 m.a.s.l. was reached at a distance of 800 m from the pit toe and was 300 m closer to the pit wall slope than in October 2010. The northern part of the pit floor was deepened from an elevation of +85 to +72 m.a.s.l. This section is located at a distance of 200 m from the pit toe. Furthermore, the pit floor reached an elevation of +41 m.a.s.l. at a distance of 1200 m from the pit toe (Figure 6C).

Over the next two years, the intensity of mining works was very high, and the front-line of the mining operations was positioned parallel to the slope edge. Lignite extraction moved from the west to the east. The second designed bench +74 m.a.s.l. was finally formed, and an elevation of +60 m.a.s.l. was reached at the pit toe. The north corner was deepened from +97 to +85 m.a.s.l. The pit floor was excavated into four elevations: +67/+41/+22/−7 m.a.s.l., which were located at a distance of 100, 400, 800, and 1100 m, respectively, from the pit wall slope (Figure 6D).

In October 2016 the third designed bench at +60 m.a.s.l. was fully formed. An elevation of +41 m.a.s.l. was achieved at the pit toe, and elevations of +22, −7, and −12 m.a.s.l. were reached at a distance of 300, 700, and 900 m, respectively, from the pit wall slope (Figure 6E). The front-line of mining operations was still parallel to the slope edge, and the lignite extraction moved from west to east.

Over the last two years of the study period (2017 and 2018), the position of the mining operation’s front-line changed after the formation of the fourth and fifth designed benches (+41 and +22 m.a.s.l.). These two benches were formed through the western pit wall slope to the northern pit wall slope of the Bełchatów field. The pit floor was formed into three elevations (−7, −30, and −40 m.a.s.l. (Figure 6F)), which were planned to be deepened in the following years (after 2018).

3.3. AWID Measurement System

The changes in stress inside the DSD are described based on data from the AWID measurement system. This system was placed in the borehole located in the central part of the DSD. This type of AWID system is equipped with six pressure pads and designed to determine the stress tensor (Figure 7). The inner pressure in the pads is produced hydraulically or pneumatically. Initially, a small defined inner pressure was injected according to the hydrostatic pressure in the rock mass. Measurements determined the difference between the inner and outer pressure [31,32,33]. Each pad was arranged independently and represented a standard component of the stress vector effecting the pad’s surface. Pads 2, 3, and 4 measured the horizontal pressure in the directions rotated 120°. Pad 1H0° measured vertical pressure; however, pads 5 and 6 were inclined +45° and −45° from vertical. The following directions of pressure were monitored by the AWID system (Figure 7):

✓

SSE–NNW by the 2V0° pad,

✓

W–E by the 3V120° pad,

✓

NNE–SSW by the 4V240° pad,

✓

NNW–SSE by the 5V0°+45° pad,

✓

SSE–NNW by the 6V0°−45° pad.

The AWID system was installed in the borehole in 2007. However, it was impossible to determine the stress tensor because the six pads were never active together. The changes in pressure registered by the pads resulted from a convergence of rock salt. The measurement data were collected four times a year.

3.4. Materials, Models, and Parameters

The mechanical behavior of the rock mass surrounding the DSD was simulated based on the constitutive models and the mechanical parameters of the rocks. Generally, the constitutive models of the materials describe the stress–strain behavior in response to the applied loads. In this paper, the Mohr–Coulomb model was chosen to reflect the complexity of mining and the geological conditions and simulate the elastic–plastic response of rocks. Rheology and dilatancy was not considered in numerical modeling due to the low value of displacements in the salt dome.

The physical and mechanical parameters of the rock salt and non-salt sediments (

Table 1. Physical and mechanical parameters used in numerical modelling.

Rock Types	Unit Weight γ [kN/m3]	Cohesion c [MPa]	Friction Angle φ [0]	Uniaxial Compressive Strength σc [MPa]	Tensile Strength σt [MPa]	Young’s Modulus E [GPa]	Poisson’s Ratio ν [-]
Zechstein rock salt	2220	5.080	46.08	25.820	2.056	4.482	0.27
Gypsum–anhydrite–clay cap	2202	0.161	24.00	12.800	0.176	8.050	0.40
Brecciated Mesozoic rocks over the DSD	2243	0.129	30.20	3.500	0.104	1.113	0.40
Brecciated Mesozoic rocks surrounding the DSD	2243	0.129	30.20	3.500	0.104	1.113	0.40
Triassic formation	1998	0.021	25.50	0.300	0.19	0.097	0.27
Jurassic formation	2700	2.000	43.50	44.200	1.006	12.000	0.40
Cretaceous formation	2345	0.616	43.50	44.20	0.320	12.000	0.40
Miocene undercoal formation	2059	0.500	22.00	0.500	0.250	0.110	0.40
Miocene coal and clay–coal formation	1590	0.149	16.50	3.400	0.219	0.162	0.32
Tertiary clay–sand complex	1998	0.070	25.50	0.300	0.019	0.097	0.27
Quaternary formation	1998	0.070	25.50	0.300	0.019	0.097	0.27

) were determined based on laboratory tests, geotechnical mapping of the cores and outcrops, and a geotechnical evaluation based on the Hoeck–Brown classification [15,16,34,35]. Based on this information, the 2D and 3D numerical models were built and calibrated by the authors for over 20 years. These assumed mechanical properties are the result of the continuous process of determining the most reliable simulation under real conditions. Relative to the complex geology and tectonics of the study area, 11 geotechnical layers were distinguished for the needs of numerical modeling. The physical and mechanical properties of each geotechnical layer were calculated as the weighted average based on the lithological diversity of each layer. Calculations were performed considering the percentage contribution of the layer thicknesses and the values of adequate parameters determined in the laboratory tests. Moreover, tectonic features such as brecciation and fracturing were reflected in the definitions of the mechanical parameters for each geotechnical layer.

3.5. Factors Analysed in Numerical Modelling

Analyzing the impact of open-pit mining and developing a stability assessment of the DSD included determining the vertical and horizontal displacements, as well as the vertical and horizontal stress. The vertical and horizontal displacements were calibrated based on the geodetic measurements performed over the salt dome and at the WSBF, as well as the measurements of the subsurface deformations from the inclinometer system [15,16]. The geodetic network is well-developed and contains 69 points. These points are distributed in four observation lines (Figure 8): axial lines No. 1, 2, and 3 in directions NW–SE, W–E, and NE–SW, with Line No. 4 encircling the DSD. Measurements were performed twice a year. The state for January 2010 was set as the foundation for the displacement calculations. The displacements before January 2010 were reset. Thus, the states and values of the displacements presented in this paper refer to the period between January 2010 and October 2018. In the vertical displacement analysis, the variations in the level of the water table had to be considered. The dewatering system in the Bełchatów mine was adapted to the geological structure of the deposit. The DSD is protected by a barrier consisting of 41 wells that are centrically distributed around the salt dome and the monitoring system (150 piezometers). During the analyzed period, the average level of the water table registered by the piezometers was stable, and the changes between successive years only reached up to 1 m. Consequently, it was determined that the changes in the water table had an insignificant impact on the vertical displacements. The vertical stress and horizontal stress were calculated to evaluate the impact of lignite extraction on the salt dome. The volume of extracted rock mass was introduced gradually (according to the mining progress) into the numerical model.

4. The Results of the DSD Stability Analysis

The results of the stability analysis are presented in the form of maps. These maps show, with the help of a color scale, the values of the analyzed factors, including the vertical and horizontal displacements and the vertical and horizontal stress in the relevant time period. Different colors represent different values. For instance, the red represents a value of +5 cm (uplifting) on the map of vertical displacements and a value of 5 cm (movement in W direction) on the map of horizontal displacements. Dark blue corresponds to a value of −4 cm (subsidence) on the map of vertical displacements and a value of −5 cm (movement in E direction) in the map of horizontal displacements. The maps present the distribution of the analyzed factors in cross-sections through the model. The locations of the W–E and S–N cutting planes within the numerical model are presented in Figure 9. Maps of the horizontal and vertical displacements were prepared for the E–W cross-section and oriented towards NE to detail all the changes in the eastern part of the model (Bełchatów field). However, the maps of vertical and horizontal stress were prepared for the W–E and S–N cross-sections to show changes in the stress field. In addition, during the analysis period (January 2010–October 2018), the mining operations in the vicinity of the DSD were performed only in the Bełchatów field. The mining activity in the Szczerców field was concentrated at a distance over 2 km from the DSD. Thus, the maps of the horizontal displacements were prepared only for the W–E cross-section. The displacements maps were prepared for the two-year periods considering the scale of the model. However, the maps of the vertical and horizontal stress are presented only for the end-period of the analysis (October 2018).

The horizontal displacements were analyzed in two directions: towards W (positive values) and towards E (negative values). Consequently, in the maps that present the W–E cross-section through the model, the x-axis corresponds to the W–E direction, and y-axis corresponds to the S–N direction. The vertical displacements correspond to the z-axis in the maps. Generally, the horizontal displacements in the W–E direction prevail over the S–N direction in the whole analysis period. For this reason, horizontal displacements in S–N direction are not presented.

In the first period (from January 2010 to October 2010), uplifting dominated in the Bełchatów field, with a maximum value of 4 mm. The vertical displacements were concentrated in two areas at distances of about 500 m and 1500 m (Figure 10A) from the DSD (each distance was calculated from the central axis of the DSD). The horizontal displacements ranged from 1 to 3 cm and comprised the same areas as the vertical displacements (Figure 11A). Both the vertical and horizontal displacements had an effect at a depth of about 300 m below the surface.

The vertical displacements largely expanded over the next two years to form a zone in the eastern part of the model. This zone comprises the major area in the eastern part of the model in October 2012. The zone of the largest displacements is about 700 m long, 200 m wide, and extends to a depth of about 400 m below the surface. The largest uplift (max. +5 cm) occurred in two areas separated by the subsidence region (Figure 10B). The largest area reaches the eastern border of the model, but a small area is located at the distance of 500 m from the DSD. The vertical displacements of 1–2 cm affected the brecciated Mesozoic sediments surrounding the salt dome. The largest subsidence of −5 cm characterized the region between the uplifting areas. The lowest subsidence (−2 cm) was observed in the surface over the salt dome. Similarly, the horizontal displacements expanded between October 2010 and October 2012. The horizontal displacements spread in the eastern direction and formed two elongated displacements in the E–W direction zones (Figure 11B). These zones expanded to a depth of about 400 m and reached the brecciated sediments surrounding the salt dome.

Over the next two years, further progress was made in expanding the vertical and horizontal displacements in the eastern part of the model and the DSD. In the map for October 2014, uplifting occurred in the entire eastern part of the model with different magnitudes distributed concentrically (Figure 10C). The largest area is observable in the zone of the largest vertical displacements (+5 cm), which extended at a distance of about 1300 m in the W–E direction and to a depth of 1400 m below the surface. The salt dome was affected by uplifting of 1–2 cm. Subsidence occurred in a small area between the uplifting zone and the brecciated sediments surrounding the DSD. In October 2014, the area of the largest horizontal displacements was an extensive zone that could be a potential slide surface (Figure 11C). As shown, this displacement is inclined according to the geometry of the Bełchatów western slope. The zone extends at a distance of 1750 m and reaches the area directly over the DSD. The displacements in this zone are directed towards the E. There is also another zone in the shape of a circle that influences the salt dome at a depth of 600–2000 m. The horizontal displacements directed towards the W occur at the surface and form a small belt adjacent to the potential slide surface.

The progress in the described trends continued over the following years. In October 2016, the vertical displacements expanded to the west. The zone of the largest uplift (5 cm) spread in a western direction and downwards (Figure 10D). As a result, the zone reached the brecciated sediments adjacent to the DSD and a depth of about 1900 m. The lower uplifting (1–3 cm) influenced the salt dome and brecciated sediments. The subsidence area also expanded to the west and increased in depth. The zone of the largest horizontal displacements (the potential slide surface) underwent the same westward tendency. However, the circular zone expanded in all directions and changed in an irregular ellipse. The value of the horizontal displacements within this ellipse ranged from 1 to 3 cm and affected the salt dome from a depth of about 400 to 1900 m (Figure 11D). Moreover, the belt of horizontal displacements directed towards the W expanded to the south.

In the last period, the same tendencies were observed in the development of the vertical and horizontal displacements. The largest uplifting zone continued to expand, and, in October 2018, comprised the entire eastern part of the model (Figure 10E). Consequently, the zones characterized by lower uplifting values extended in a western direction. The horizontal displacements underwent similar changes. The elliptical area spread further to the east, and the maximum value of displacements increased to 4 cm. The horizontal displacements directed toward the E expanded to a depth of 2200 m (Figure 10E).

The stress state was stable during the analysis period (January 2010–October 2010), and there were no significant changes in the stress field around the DSD. Thus, the presented maps depict the stress state in October 2018. The vertical stress depends on the overburden pressure and increases with depth from 0 on the surface to 59 MPa in the lowest parts of the model (Figure 12). The stress field is disturbed in the upper part of the model. The low weight density of rock salt and the halokinetic movements caused the distribution of the vertical stress in the DSD to be different from that in the surrounding rocks. In the upper part of the DSD, the disturbances of the stress field are larger than those in the surrounding rocks. Horizontal stress is related to vertical stress, thus, both values increase with depth. However, the irregular distribution of horizontal stress in the upper part of the model is larger than that for the vertical stress (Figure 12). These disturbances are visible in the N–S cross-sections (Figure 12) and occur mostly in the salt dome and adjacent area.

The presented maps of vertical and horizontal stress show the stress state in the DSD and the surrounding areas on a large-scale. Consequently, the stress changes in the range of several MPa cannot be detected at this scale. Thus, an analysis of data from the AWID measurement system is presented below.

The data registered by the AWID system give information about the value of stress, as well as the directions of and changes in stress (Figure 7). From 2010 to 2016, four pads measured the pressure in the DSD; however, since the end of 2016, only three pads were active (Figure 13 and Figure 14). During 2010–2018, only the 3V120° pad was active for the whole period (Figure 13 and Figure 14). The periods of inactivity and their durations are different for each pad. The least amount of data was registered by the 2V0° pad, which was active only temporarily from September to December 2013 (Figure 13C). On the contrary, the 5V0° + 45° pad displayed only one period of inactivity from January to June 2010 (Figure 13A). The horizontal pad 1H0° was inactive between April 2010 and July 2011 (Figure 13A,B). Inactivity periods were registered for the 6V0°–45° pad from January to April 2013 and from June 2013 to March 2018, as well as for the 4V240° pad from December 2016 to March 2018 (Figure 13A–E).

The values of the pressures registered by the pads equal 1–3 MPa and underwent constant fluctuations (Figure 13). Pressure above 2 MPa was registered by the 5V0° + 45° pad since July 2011 and by the 4V240° pad from May 2012 to March 2016 (Figure 13B–E). The other active pads measured lower pressure values. In addition, two major changes in pressure were registered for all active pads. Firstly, an increase (+1.0 MPa) was recorded between April 2010 and June 2010 (Figure 13A and Figure 14A). Secondly, a sudden increase (+1.0 MPa) followed by a drop (−1.5 MPa) was registered from September 2013 to June 2014 (Figure 13C and Figure 14C). The recorded data characterize the high variability of values and the large differences between the two consecutive measurements for a single pad. The same trends are related to the measurements from all pads. These changes are related to the mining activity in the WSBF.

5. Stability Analysis in Comparison to the Mining Progress

The intensity and distribution of the mining operations in the WSBF from January 2010 to October 2018 were examined to assess their impact on the behavior of the rock mass, including the rock salt. This assessment was divided into two-year periods in accordance with the analysis of displacements. As non-renewable resources continue to be key in the world’s demand for energy [1], it is crucial to continue to ensure the safety and stability of mining regions.

5.1. Mining Progress Influence on the Analysed Factors

In the first period (January–October 2010), lignite extraction was performed in close (200 m) vicinity to the DSD at a distance of 800–1800 m (Figure 6A,B). Mining operations were concentrated in the northern and northeastern surroundings of the DSD (Figure 6A). This mining activity resulted in the development of horizontal and vertical (uplifting) displacements in the Bełchatów field. The area of the largest horizontal and vertical displacements was located at a distance of about 500 and 1500 m (Figure 10A and Figure 11A) from the DSD (this distance was calculated from the central axis of the DSD). The mining activity and the occurrence of displacements had an influence on the salt rock mass of the DSD and was registered by the AWID system. In the first part of 2010, a pressure increase was recorded in all four active AWID pads. However, in the next three months, the pressure decreased (Figure 13A and Figure 14A).

In the following years (2011–2018), the mining operations were directed to the deepening of the pit floor in the analyzed area. Mining operations moved from the east to the west in order to reach the design slope’s geometry. In the period from October 2010 to October 2012, lignite extraction was carried out first in the north–eastern vicinity of the DSD at a distance of about 400 m and 1000–1800 m (Figure 6C and Figure 14B). Then, the mining operations were concentrated in the eastern surroundings of the salt dome at a distance of about 800 m. These works resulted in the formation of two uplifting zones separated by the subsidence region (Figure 10B), as well as the occurrence of horizontal displacements (Figure 11B). In comparison to the previous period, the zones of the vertical and horizontal displacements expanded mainly in the W–E direction (and in depth), reaching the brecciated Mesozoic sediments in the surrounding salt dome. Two decreases in pressure were recorded by the AWID pads: once at the end of 2010 and once in the second part of 2012 (Figure 13B and Figure 14B). The largest decrease was registered by the 4V240° pad that measured pressure in the NNE–SSW direction and the 1H0° pad that measured the vertical pressure (Figure 14B). At the beginning of 2012, the increasing trend was noticed for all active pads.

The two following years (from October 2012 to October 2014) showed the development of the tendencies described above. Mining activity was undertaken in the east at a distance of 800 m from the DSD and in the north–east at a distance of 600–1200 m from the DSD (Figure 6D). The intensity of mining works was very high. As a result of this intensive extraction, the zones of the vertical and horizontal displacements largely expanded (Figure 10C and Figure 11C). At the end of this period, the uplifting zone comprised most of the analyzed area. However, the largest vertical displacements covered an area about 1300 m long and 1000 m wide. Moreover, the zone reached a depth of about 1400 m below the surface. The subsidence zone enlarged and reached the brecciated Mesozoic sediments over the salt dome (Figure 10C). The area of the largest horizontal displacements formed an extensive zone, making it a potential slide surface (Figure 10C). The horizontal and vertical displacements affected the salt dome and brecciated the Mesozoic sediments surrounding the DSD (Figure 10C and Figure 11C). These processes were reflected in the AWID pad activity. A major decrease in pressure was noted at the end of 2013 and registered by the 1H0° pad that measures the vertical pressure and the 3V120° pad whose measurements are directed W–E. Moreover, the 2V0° pad that measures pressure in the SSE–NNW direction was activated for 6 months, and the fluctuations in the value of pressure were recorded by the other active pads (Figure 14C).

Further progress in the described trends occurred over the next two years (between October 2014 and October 2016). Lignite extraction was concentrated in the northern and northeastern vicinity of the salt dome at a distance of about 600–2000 m (Figure 6E). Mining activity contributed to the enlargement and expansion of the vertical and horizontal displacement zones and the areas with the highest displacement values (Figure 10D and Figure 11D). The uplifting zone covered the whole analyzed area of the Bełchatów field and the salt dome (Figure 10D). Moreover, the subsidence zone located in the south–eastern surroundings of the DSD extended over the salt dome. The active AWID pads recorded two decreases in pressure: one at the start of 2014 and 2015 and another during the first months of 2016. The largest decreases were registered by the 1H0° pad that measures the vertical pressure, the 3V120° pad that measures pressure in the W–E direction, and the 4V240° pad that measures pressure in the NNE–SSW direction (Figure 14D).

In the last analyzed period (from October 2016 to October 2018), the mining operations were located at a distance of about 1200–1800 m in the east and north–east directions from the DSD (Figure 6F). A further development was observed in the vertical and horizontal displacements according to the described trends (Figure 10F and Figure 11F). These processes were followed by the AWID pad activity. At the end of 2016 and during the first six months of 2017, a decrease in pressure was recorded by the pads: 3V120° (W–E), 1H0°, 4V240° (NNE–SSW), and 5V0° + 45° (NNW–SSE). However, after 2017, only three pads were active. Fluctuations in the values of pressure occurred over the following months of 2017, as well as in 2018 (Figure 14E).

For the pressure changes registered by the AWID measurement system, the general tendency in the stress field of the DSD and its surroundings was observable in the maps (Figure 12). These maps clearly outline the differences in the distribution of vertical and horizontal stress between the W–E and S–N cross-sections (Figure 12). In the W–E cross–sections, vertical and horizontal stress were distributed almost evenly inside the DSD and its surroundings. Some disturbances in the stress field are observable in the upper and eastern parts of the maps (Figure 12). These disturbances are associated with the Mesozoic brecciated rocks surrounding the DSD and result from mining operations. However, in the S–N cross sections, the distribution of the vertical and horizontal stress inside the DSD is asymmetric (Figure 12). The disturbances in stress field occurred inside the DSD, and their intensity decreased with the depth. This asymmetry is related to the lignite extraction performed in the eastern and north–eastern vicinity of the DSD.

5.2. General Trends and Observations

The results of the 3D numerical modeling indicate that the lignite extraction and mining activity in the WSBF had an impact on the Debina salt dome and the surrounding rocks. The stability of the WSBF and the salt dome was not affected for the assumed conditions and time. The numerical analysis revealed the development of horizontal and vertical displacements, as well as the disturbances in the stress field that resulted from mining operations. Lignite extraction required unloading of the rock mass in the mining field. This rock mass relaxation triggered the uplifting (vertical displacements) in most of the analyzed area. The vertical displacements reached a large depth (about 2 km below the surface) resulting from the large volume of rock mass that was extracted. Moreover, this process was accelerated by the present halokinetic activity of the salt dome [36,37,38]. The deepening of the pit floor contributed to the enlargement and expansion of the vertical and horizontal displacements, mainly in the eastern direction. In addition, the development of these displacements was accelerated by the geological structure of the area, especially the occurrence of the brecciated Mesozoic sediments that surround the salt dome, the regional discontinuities, and the local faults, as well as the lithological variability of the rocks. The trends for the development of displacements comply with the direction of the graben border faults, whose quaternary activity was indicated by the measurements and observations [39]. Moreover, the unloading of the upper layers caused by lignite extraction may contribute to the reactivation of these regional discontinuities. The lithological and petrological diversity of the rocks in the analyzed area is associated with the mechanical parameters of the rocks. The co-occurrence of rock salt, lignite, limestone, marls, clays, sand, and mudstone, as well as the petrological variations within all these types of rocks, contributed to the high diversity in the mechanical parameters. However, the largest value of the displacements under numerical modeling was low. Thus, the indicated displacements are not associated with a risk of stability loss or landslide occurrence.

The described processes were followed by changes in the stress measured by the AWID system that was installed in the borehole inside the salt dome. The rock salt reaction on the mining operations and the rock mass relaxation that resulted from the lignite extraction were recorded by the AWID pads as a pressure increase followed by a decrease. It should be noted that the reaction of the rock salt was delayed over time and registered by the AWID pads after several months. The value of the pressure in the AWID pads was low for the values of the largest displacements. Moreover, the active pads provided information about the approximate area of mining activity. For example, the 4V240° pad that measured pressure in the NNE–SSW direction was active when mining works were concentrated in the northeastern vicinity of the salt dome. The activity of the AWID pads was associated with the intensity of mining works. All six pads were active between October 2012 and October 2014 when the lignite extraction was the most intensive. The salt domes are generally considered to be homogenous. However, taking into consideration their genesis, the results indicate that the salt’s movement upward over time took place as a series of salt spines. Bordering shear zones occurred between these differentially moving salt spines. These features of the salt dome’s internal structure impacted its behavior [40,41]. Consequently, the internal structure of the DSD may affect its impact on mining works concentrated in its vicinity. In addition, some irregularities and asymmetry in the stress distribution in the salt dome and surrounding area are indicated by the data. These disturbances in the stress field are a result of the lignite extraction process in the western slope of the Bełchatów field.

6. Conclusions

The results of the numerical simulations performed for the large-scale pit wall slope are characterized by very unique geological conditions indicating complex rock mass behavior. The numerical simulations presented in this paper were confirmed by measuring the subsurface deformations using an inclinometer system, as well as from data from the geodetic network. The numerical model of the western slope of the Bełchatów field and the Dębina salt dome, as well as the surrounding area, has a broad scope of applications covering not only the simulation and prediction of rock mass behavior, but also the planning of mining operations and different technological solutions. These applications are not limited to mining operations but also apply to the future utilization of the analyzed area, such as the technological development and trends in underground energy storage, including salt domes, as well as the various methods of post-mining area recultivation. The world is transitioning towards renewable energy sources, but as this transition takes place, raw materials will continue to be key in achieving the Sustainable Development Goals [42]. Thus, the most accessible future-oriented approach for existing mining projects is advisable. Taking into account the progress of mining works in the Szczerców field and the deepening of the pit floor in the western slope of the Bełchatów field, lignite’s extraction impact on the salt dome is expected to increase. Consequently, the rheological processes in the salt dome should be considered in future numerical simulations.