Subsidence Management and Prediction System: A Case Study in Potash Mining

: Subsidence is an important environmental and safety issue in the mining sector, yet there remain voids in knowledge in terms of management and prediction. This study aims to improve knowledge on the impact of mining operations on the surface, reducing their effect on the environment, increasing the safety of mining operations, monitoring stress behavior and predicting rock mass. Therefore, an analysis was carried out to process and analyze the measured subsidence data and, subsequently, create a numerical model to predict the surface subsidence of a case study mine. The model was developed based on a ﬁnite element method (FEM). It was achieved by considering the geological characteristics of the area, the design features of the mine, the surface subsidence measured over twelve years and the time-dependent behavior of the geological layers. The correlation obtained between the measured subsidence and the modelling results was very satisfactory, with a 90% conﬁdence level, over the years analyzed. Hence, the efﬁciency of the system was conﬁrmed, enabling the evaluation and the prediction of potential surface effects, and therefore improving the safety and environmental levels of the mining area.


Introduction
Mineral resources are crucial to providing welfare to society, so it is necessary to find and extract them on a sustainable basis. Nowadays, the mining industry represents a key factor in the development of societies and communities around ore deposits [1,2]. However, mining activity must be managed properly, since the effects of the extraction process can produce serious environmental problems in the vicinities of the mining area, such as acid drainage, heavy metal contamination of soils or subsidence [3,4]. In recent decades, various initiatives have emerged to construct the legal frameworks and legislate on the sustainability of extractive activities [5][6][7].
In the case of subsidence, it is well known that when mineral extraction is not properly planned, it has a high probability of having a direct impact on the surface, affecting not only mining infrastructures but also urban infrastructures or the communities that live close to the mining area [8][9][10]. In addition, subsidence has become of special importance as it is closely linked to the overexploitation of groundwaters [11][12][13]. Groundwaters represent one of the most important natural resources for ecosystems, as stated in the Sustainable Development Goals (SDGs) declared by the United Nations [14]. Various mining and geological factors, including the quality and characteristics of the underground rock and the surface conditions, can influence the subsidence magnitude, shape, mode, and extent [15]. Therefore, it is necessary to consider the usage of less aggressive mining methods as part of corporate social responsibility (CSR) [16].
The need for mineral resources pushes the industry towards extracting ores at greater depths, with buildings and infrastructures on the surface or with challenging bedrocks This study was developed based on data from a mine located in northeastern Spain which extracts potash using a room and pillar method at depths between 500 and 700 m, varying the extraction rate depending on the arrangement of the stratum. The ore deposit is made of two different mineable layers, Bed A and Bed B, with an additional intermediate layer of varying thickness made of salt and clay minerals. Two different drifts are used to extract the potash ore when the salt layer is more than 5 m thick, while it is mined jointly when the thickness is less than 5 m. Bed A and Bed B are interstratified layers of sylvinite rock salt and clay minerals. The average production during the period 2007-2019 was 2,870,570 tons.

Data Management
The case study has approximately 700 control points in an area of 50 km 2 , with measures taken from 2008 to 2021. These control points are distributed throughout the geographical area potentially affected by subsidence, obtaining the X, Y, and Z coordinates of each point using GPS. The post-process method is employed to obtain the coordinates, using four static differential GPS with double frequency receptors. Two receptors are placed at two points of known and well-established coordinates, while the other two receptors are used to perform measurements in the control points, with a minimum measurement time per point of 12 min. This methodology guarantees an accuracy of 1-1.5 cm in planimetric and altimetric coordinates, which implies that the maximum possible error between measures from the same point is 2-3 cm. Consequently, only higher differences are relevant to indicate if there is subsidence. Regular quality control measurements are taken to verify the reliability of the measurements.
InSAR images are also used, since 2016, as an additional control system. Nonetheless, the surface characteristics of the area, mainly forested and agricultural, still require a considerable amount of GPS control points to obtain reliable subsidence displacement. In this case, it was determined that a minimum of 200 points are required. The accuracy of this type of technology was widely proved [50][51][52]. The combination of both systems allows for obtaining extremely reliable information regarding X, Y, and Z surface displacements [39].

Data Processing
The case study area is defined by the fully exploited zones from 2008 to 2018 based on the subsidence data collected. The initial year, 2008, is determined as it is the first year in which surface subsidence data are obtained. The final year, 2018, is determined according to Sanmiquel et al., 2018 [39], where it is stated that 90% of the surface subsidence for an active mine occurs within the initial 5 years.
Five profiles are defined within the area of interest. Each profile crosses the area completely, every 10 m in the direction of the mine advance, east-west, and Metric Points (MPs) are listed. These MPs show the displacement in the Z coordinate for each subsidence calculation period. According to Sanmiquel et al. [39], a 400-meter protection area is determined in order to separate the targeted area from the area exploited before 2008, as no data is available ( Figure 1).  After analyzing all the raw data collected, the MPs are reduced and re-listed every 100 m, since it is considered an adequate distance to appreciate the evolution of subsidence. Figure 1 shows the targeted area in the enclosed orange perimeter and the mining drifts excavated during the years 2008-2018 are marked in light green. The light blue perimeter delimits the area exploited before 2008. Pink lines depict the 5 profiles used to determine the average subsidence of the targeted area, while dark blue marks the 12 mining zones selected to test the model, as they have been sorted from north to south and east to west, respectively. On the other hand, Figure 2 shows the targeted area amplified, so that the five profiles and the 12 selected mining zones can be seen in more detail as well as the 400-meter protection area.  The analysis is initially performed individually, defining how subsidence evolves over annual and biannual measurements. Subsequently, subsidence values are gathered together to obtain the profile of maximum subsidence for each MP and perform a statistical analysis to determine the average subsidence profile of the targeted area, as well as the level of confidence of the subsidence measured. Figure 3 represents the methodology to The analysis is initially performed individually, defining how subsidence evolves over annual and biannual measurements. Subsequently, subsidence values are gathered together to obtain the profile of maximum subsidence for each MP and perform a statistical analysis to determine the average subsidence profile of the targeted area, as well as the level of confidence of the subsidence measured. Figure 3 represents the methodology to follow and obtain the average surface subsidence profile. The characteristics of the case study are in brackets.

Numerical Model Definition
The numerical model was generated using indirect time-dependent behavior modelling, using the RS2 v11.003 software (Toronto, ON, Canada). After analyzing the geometries of the drifts and the extraction ratio, an initial model was run. The total size of the model is 225 m by 666 m, with a minimum of 37.5 m of walls on each side; therefore, 150 m of exploitation is modelled. The model size is based on the stratigraphic column size proposed by Orellana et al., 1996 [53]. The models developed by Xu et al. [31] and Parmar

Numerical Model Definition
The numerical model was generated using indirect time-dependent behavior modelling, using the RS2 v11.003 software (Toronto, ON, Canada). After analyzing the ge-ometries of the drifts and the extraction ratio, an initial model was run. The total size of the model is 225 m by 666 m, with a minimum of 37.5 m of walls on each side; therefore, 150 m of exploitation is modelled. The model size is based on the stratigraphic column size proposed by Orellana et al., 1996 [53]. The models developed by Xu et al. [31] and Parmar et al. [30] have been used to define the wall width. The drifts are designed with an average geometry section, 6 m in height and 8.157 m in width, based on real mine planning. The mesh consists of 6 nodal triangles, while it is graded and densified around 150 m of the mined area, as well as in the areas of lithological variations to adequately capture displacements adjacent to the drift wall. No support material was placed in the drifts. The Mohr-Coulomb constitutive model was adopted for the analysis.
The model was designed based on the stratigraphic column taken from Campos de Orellana [53] and subsequently modified by Cendon et al. [54], adding two horizons of Anhydrite and Halite on top of those, known as roof salt. RS2 model is presented in Figure 4, where all the materials can be seen. The first lithology consists of sedimentary material such as sandstones, limestones and lutites, the second group of lithologies consists of two layers of anhydrite. These two groups of stratums correspond to the overburden part of the model. The following layers are defined by the potassic seams, formed by carnallite, Bed B, intermediate salt, Bed A, and footwall salt, and finally comes the last lithology formed by the lower salt. Figure 4 is not to scale, the top and bottom have been cut off, only showing the middle part of the model, which is considered the representative part of it. The geomechanical values of the lithologies were determined by laboratory and field testing, while some data were taken from Campos de Orellana [53]. Table 1 shows the values of geomechanical properties for each lithology.   The materials with time-dependent behavior are the following lithologies: carnallite, Bed B, intermediate salt, Bed A and footwall salt. Seven stages have been defined, each one representing a two-year period. As can be seen in Figure 5, in the first stage, the original properties of the geological materials have been used, while Young's modulus of the materials with time-dependent behavior was reduced from the second stage onwards.  The materials with time-dependent behavior are the following lithologies: carnallite, Bed B, intermediate salt, Bed A and footwall salt. Seven stages have been defined, each one representing a two-year period. As can be seen in Figure 5, in the first stage, the original properties of the geological materials have been used, while Young's modulus of the materials with time-dependent behavior was reduced from the second stage onwards. Figure 6 illustrates the methodology to define Young's modulus decrease. The first modulus reduction was defined according to Paraskevopoulou et al. [21], while the following reductions are determined by an iterative approach (Figure 6), comparing the actual subsidence data and RS2 values. The factor used in each stage to reduce Young's modulus is included in Figu weakening the material properties of the model as a time-dependent behavior. The d mination of the geomechanical characteristics of each stage is verified using the su subsidence values. Once all the steps of the initial model were verified, it was tested selected mining zones, all of them included in the target area. The selection of these z was performed randomly in order to achieve better representativity of the procedure posed and its usage.

Results and Discussion
The technique to indirectly assess the time-dependent behavior allows for obta the following outcome. Results related to the definition of the average subsidence p are presented, followed by the results obtained with RS2. Figures 7 and 8 show how sidence evolves, at point MP 03 (which corresponds to the first 300 m), in the five pr tested between 2008 and 2021. The relative values are represented in Figure 7, sho the magnitude of surface subsidence year by year. On the other hand, Figure 8 dis The factor used in each stage to reduce Young's modulus is included in Figure 5, weakening the material properties of the model as a time-dependent behavior. The determination of the geomechanical characteristics of each stage is verified using the surface subsidence values. Once all the steps of the initial model were verified, it was tested in 12 selected mining zones, all of them included in the target area. The selection of these zones was performed randomly in order to achieve better representativity of the procedure proposed and its usage.

Results and Discussion
The technique to indirectly assess the time-dependent behavior allows for obtaining the following outcome. Results related to the definition of the average subsidence profile are presented, followed by the results obtained with RS2. Figures 7 and 8   The factor used in each stage to reduce Young's modulus is included in Figure 5, weakening the material properties of the model as a time-dependent behavior. The determination of the geomechanical characteristics of each stage is verified using the surface subsidence values. Once all the steps of the initial model were verified, it was tested in 12 selected mining zones, all of them included in the target area. The selection of these zones was performed randomly in order to achieve better representativity of the procedure proposed and its usage.

Results and Discussion
The technique to indirectly assess the time-dependent behavior allows for obtaining the following outcome. Results related to the definition of the average subsidence profile are presented, followed by the results obtained with RS2. Figures 7 and 8 show how subsidence evolves, at point MP 03 (which corresponds to the first 300 m), in the five profiles tested between 2008 and 2021. The relative values are represented in Figure 7, showing the magnitude of surface subsidence year by year. On the other hand, Figure 8 displays the cumulative evolution of subsidence at point MP 03 during the period studied, obtaining a representation of the maximum values of subsidence and how surface subsidence evolves from year to year. In both figures, it can be observed how the maximum vertical displacement occurs during the years 2012 to 2016, with the period 2012-2014 the one with the steepest slope, hence the one with the maximum subsidence. This behavior is in accordance with previous studies performed in the area [39]. From 2016 onwards, it can be observed how the vertical subsidence is progressively reduced until it became stable during the period 2018-2021. the steepest slope, hence the one with the maximum subsidence. This behavior is in accordance with previous studies performed in the area [39]. From 2016 onwards, it can be observed how the vertical subsidence is progressively reduced until it became stable during the period 2018-2021. Profile 4 shows the clearest subsidence evolution since it can be observed how the subsidence velocity evolves progressively in the first stages, increasing until reaching its maximum value of it in the period 2014-2016, followed by a reduction in the period 2016-2020, until reaching the stabilization of the subsidence velocity in the period 2020-2022. This subsidence trend is possible due to the fact that it is the profile located in the south part of the mine (Figure 1), being the least influenced by the whole group of mining drifts and excavations. Its analysis shows how the subsidence increases until reaching the maximum vertical displacement in 2016. From this date forward, vertical displacement is gradually reduced, reaching a stabilization stage during the period 2018-2021. On the other hand, the other profiles are located in areas where mining was more sustained in the period 2019-2020 and forward, making it more difficult to see the stabilization progress. Profile 4 shows the clearest subsidence evolution since it can be observed how the subsidence velocity evolves progressively in the first stages, increasing until reaching its maximum value of it in the period 2014-2016, followed by a reduction in the period 2016-2020, until reaching the stabilization of the subsidence velocity in the period 2020-2022. This subsidence trend is possible due to the fact that it is the profile located in the south part of the mine (Figure 1), being the least influenced by the whole group of mining drifts and excavations. Its analysis shows how the subsidence increases until reaching the maximum vertical displacement in 2016. From this date forward, vertical displacement is gradually reduced, reaching a stabilization stage during the period 2018-2021. On the other hand, the other profiles are located in areas where mining was more sustained in the period 2019-2020 and forward, making it more difficult to see the stabilization progress.
The average subsidence profile of the area of interest can be determined by means of the average subsidence profile from each MP. Figure 9 shows the average subsidence profile with a 90% confidence level. It can be seen that the maximum variability of the analyzed values is detected around the years 2016 and 2020 since these are the years in which the raw subsidence data have the greatest variability. The data detailed in Figure 9, and processed, have been used to define the decreasing Young's modulus ( Figure 5). As can be seen in Figure 6, the reduction in Young's modulus in each stage varies from that observed by Paraskevopoulou et al. [21], since it was adapted to the subsidence conditions of the case study. Nevertheless, the fundamental procedure applied remains the same. The average subsidence profile of the area of interest can be determined by means of the average subsidence profile from each MP. Figure 9 shows the average subsidence profile with a 90% confidence level. It can be seen that the maximum variability of the analyzed values is detected around the years 2016 and 2020 since these are the years in which the raw subsidence data have the greatest variability. The data detailed in Figure 9, and processed, have been used to define the decreasing Young's modulus ( Figure 5). As can be seen in Figure 6, the reduction in Young's modulus in each stage varies from that observed by Paraskevopoulou et al. [21], since it was adapted to the subsidence conditions of the case study. Nevertheless, the fundamental procedure applied remains the same. Surface subsidence values have been used to calibrate the FEM model by comparing them with vertical displacement values located on the surface of the model. Figure 10 shows one example of the stages obtained after the FEM post-process calculation, where it can be seen how the vertical displacement evolves throughout the stratigraphy that forms the ore deposit. These values, in all the evaluated mining areas, are homogeneously Surface subsidence values have been used to calibrate the FEM model by comparing them with vertical displacement values located on the surface of the model. Figure 10 shows one example of the stages obtained after the FEM post-process calculation, where it can be seen how the vertical displacement evolves throughout the stratigraphy that forms the ore deposit. These values, in all the evaluated mining areas, are homogeneously shown on the surface of the model, as can be seen in Figure 10. Yet, the actual subsidence values show variations in all the evaluated MPs. Thus, the average subsidence profile must be used to verify and calibrate the model (Figure 9). The model is proposed to enable the analysis of subsidence scenarios depending on the extraction ratio, depth and size of the drifts in a very user-friendly approach. The idea is that both the parameters to be introduced and those returned by the model are as simple as possible to be applied in the day-to-day life of a mining company, which represents the main advantage of the model. Therefore, the predicted values of the vertical surface displacements are considered average values. The model was tested in 12 different cases, and as a result, all tests were verified with real subsidence data.  Table 2 shows the values obtained for the twelve mining zones tested, the average subsidence profile with a 90% level of confidence. It can be seen that in most of the stages, the values obtained with the RS2 model are within the level of confidence (Figure 11). In the first stage, the vertical displacement obtained by the FEM model ranges between −14.2 and −15.2 cm, being slightly out of the interval of −8.7 and −9.7 cm defined by the level of confidence of the actual values. This may be due to the values used as original material  Table 2 shows the values obtained for the twelve mining zones tested, the average subsidence profile with a 90% level of confidence. It can be seen that in most of the stages, the values obtained with the RS2 model are within the level of confidence ( Figure 11). In the first stage, the vertical displacement obtained by the FEM model ranges between −14.2 and −15.2 cm, being slightly out of the interval of −8.7 and −9.7 cm defined by the level of confidence of the actual values. This may be due to the values used as original material properties from laboratory tests could have some variation from their real in situ characteristics. However, the variation is still quite small and corrected over time with the subsidence evolution model. The main limitation is detected in the final stage, where subsidence stabilization should be modelled. However, achieving this trend requires recovering the initial properties in the last stage, which leads to an inconsistency with Young's modulus reduction technique. Another disadvantage of the software used is the need to simplify the model design. Therefore, if much detail is required, this fact can represent a limitation. In the case of the twelve models tested, only the main geological structures were detailed. limitation. In the case of the twelve models tested, only the main geological structures were detailed.  Another factor to consider is the subsidence measured on the surface, using GPS. As mentioned in the methodology, it has a potential error between 2 and 3 cm. This means Figure 11. Average subsidence profile with the results of the sections analyzed.
Another factor to consider is the subsidence measured on the surface, using GPS. As mentioned in the methodology, it has a potential error between 2 and 3 cm. This means that the differences between real and theoretical values could be accepted when a range of ±5 cm is not exceeded and can, therefore, be considered accurate.
Based on the results from the numerical modeling, the 2D FEM method allows handling the nonlinear and inhomogeneous properties of rock masses that produce the timedependent behavior. The model is a useful tool for the prediction of subsidence in mining infrastructures [30,31], as well as for analyzing and testing the time-dependent behavior [21]. It was detected that the most affected areas by subsidence are those with the widest mining galleries. This type of gallery is usually built when it is necessary to connect two parts of the mine infrastructure. In these cases, it is recommended to use underground supports and to have sufficiently wide pillars to limit their effect on the surface. In order to gain knowledge and improve the approach proposed, it would be interesting to test other areas, either potassic deposits or other rock types, with time-dependent behaviors.

Conclusions
The approach proposed gives the possibility to determine the evolution of subsidence over time. The system allows the generation of average profiles at different points and an average subsidence profile of an area. The system also provides crucial information on mining activity behavior and gives the possibility to analyze the subsidence velocity of a representative period in a very simple and straightforward way, which represents one of the major strengths of the assessed technique. The tested numerical method has sufficient capacity to evaluate the time-dependent behavior of the potash and salt layers, using an indirect method-a 2D FEM model.
Although the model has some limitations in achieving subsidence stabilization in the final stage, it was proven as a particularly useful tool, especially during the years of mining activity. The model is able to analyze potential extraction scenarios and measures taken to reduce the subsidence whenever it is necessary and, therefore, reduce the environmental impact of the mining extraction and increase the operational safety levels.