Assessment of the Impact of the Spatial Extent of Land Subsidence and Aquifer System Drainage Induced by Underground Mining

: The environmental impact assessment of underground mining usually includes the direct e ﬀ ects of exploitation. These are damage to rock mass and land subsidence. Continuous dewatering of the aquifer system is, however, necessary to carry out underground mining operations. Consequently, the drainage of the aquifer system is observed at a regional scale. The spatial extent of the phenomenon is typically much wider than the direct impact of the exploitation. The research presented was, therefore, aimed at evaluating both the direct and the indirect e ﬀ ects of underground mining. Firstly, the spatial extent of land subsidence was determined based on the Knothe theory. Secondly, underground mining-induced drainage of the aquifers was modeled. The 3D ﬁnite-di ﬀ erence hydrogeological model was constructed based on the conventional groundwater ﬂow theory. The values of model hydrogeological parameters were determined based on literature and empirical data. These data were also used for model calibration. Finally, the results of the calculations were compared successfully with the ﬁeld data. The research results presented indicate that underground mining’s indirect e ﬀ ects cover a much larger area than direct e ﬀ ects. Thus, underground mining requires a broader environmental assessment. Our results can, therefore, pave the way for more e ﬃ cient management of groundwater considering underground mining. exploitation 3B and It consists of 56 piezometers located in the Lower Cretaceous-Upper Jurassic and Middle Jurassic aquifers. continuous monitoring of the groundwater table in the Quaternary-Upper Cretaceous aquifer is carried out in 34 wells distributed irregularly


Introduction
The environmental impacts of mining can occur through direct and indirect mining practices at local, regional, and global scales [1]. Impacts can result in erosion [2], sinkholes [3], loss of biodiversity, or soil [4], groundwater [5], and surface water contamination from the chemicals emitted from mining processes [6]. Nevertheless, one of the most significant direct impacts of underground mining is land subsidence, which can cause surface and underground infrastructure damage, and even threat to surface user safety [1,7,8]. The spatial extent of land subsidence and its value depends primarily on the depth of exploitation, the seam thickness, rock mass geomechanical conditions, and post-mining void reclamation method [7]. During recent years, a lot of research has been presented documenting the phenomenon of land subsidence resulting from mining [9][10][11][12][13][14][15]. It can therefore be stated that the spatial extent of land subsidence can cover an area of approximately several km 2 [9] to tens of km 2 [16], whereas the phenomenon values may vary from a few mm/year [17] to several cm/day [14].
Given the increasing importance of issues related to groundwater resources [66] and the environmental impact of underground mining [1], the main objectives of the paper comprise the following: • To determine the spatio-temporal extent of land subsidence induced by underground mining; • To model aquifer system drainage resulting from the studied phenomenon; • To analyze and discuss the conjugated regional-scale impact of underground mining on the environment.
The present research is conducted at a hard coal extraction mine. The mining operation in this area is carried out at a considerable depth of about 1 km, in complex geological and mining conditions, changing over time and space (varying hydrogeological parameters of the rock mass). In the case of aquifer deformation resulting from rock mass drainage induced by underground mining, the advanced numerical groundwater flow models are currently not commonly used. The lack of a uniform algorithm for modeling such deformations and their negative impacts demonstrate the need for continuous optimization and the search for new solutions that will allow more effective management of the indirect effects of mining operations on the aquifer system.

Study Area
The research works presented in this paper were carried out at Bogdanka Hard Coal Mine in the eastern part of Poland, approximately 30 km NE from the city of Lublin ( Figure 1A). The topography of the study area is relatively flat and gently sloping. The surface morphology was shaped mainly by the fluvioglacial processes associated with the Quaternary glaciations. The altitude ranges from 170 to 191 masl (meters above sea level). The research area is mainly covered by arable fields, forests, and lakes ( Figure 1B). The buildings comprise small villages and towns of between 500 and 18,600 inhabitants [67]. The surface and underground infrastructures are relatively poorly developed.
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 28 hydrogeological models of the underground mines are partially based on literature data and the results of previous research carried out at mining sites [60]. Moreover, the hydrogeological parameters of the aquifer systems can be determined, e.g., based on InSAR (Interferometric Synthetic-Aperture Radar; [5,35,[61][62][63][64][65]). Given the increasing importance of issues related to groundwater resources [66] and the environmental impact of underground mining [1], the main objectives of the paper comprise the following: • To determine the spatio-temporal extent of land subsidence induced by underground mining; • To model aquifer system drainage resulting from the studied phenomenon; • To analyze and discuss the conjugated regional-scale impact of underground mining on the environment.
The present research is conducted at a hard coal extraction mine. The mining operation in this area is carried out at a considerable depth of about 1 km, in complex geological and mining conditions, changing over time and space (varying hydrogeological parameters of the rock mass). In the case of aquifer deformation resulting from rock mass drainage induced by underground mining, the advanced numerical groundwater flow models are currently not commonly used. The lack of a uniform algorithm for modeling such deformations and their negative impacts demonstrate the need for continuous optimization and the search for new solutions that will allow more effective management of the indirect effects of mining operations on the aquifer system.

Study Area
The research works presented in this paper were carried out at Bogdanka Hard Coal Mine in the eastern part of Poland, approximately 30 km NE from the city of Lublin ( Figure 1A). The topography of the study area is relatively flat and gently sloping. The surface morphology was shaped mainly by the fluvioglacial processes associated with the Quaternary glaciations. The altitude ranges from 170 to 191 masl (meters above sea level). The research area is mainly covered by arable fields, forests, and lakes ( Figure 1B). The buildings comprise small villages and towns of between 500 and 18,600 inhabitants [67]. The surface and underground infrastructures are relatively poorly developed.

Regional Geology and Hydrogeology
Rock strata grouped in the Quaternary, Cretaceous, Albian, and Jurassic formations dominate in the analyzed area (Table 1; [68,69]). They constitute a caprock for the extracted Carboniferous coal beds. Its thickness totals 650 m in the eastern part of the research area to monoclonal increase to about 745 m in the west. The hard coal deposit is located in a synclinorium characterized by a flat bottom and asymmetrical system of wings-small dip angles in the NE wing and a steep course in the SW direction. The stratigraphic throw rates of recognized faults are relatively small, i.e., less than 2.6 m ( Figures 1C and 2). Table 1. Lithology and thickness (m) of the geological strata in the study area in selected boreholes B1-B4 (see Figure 1 for boreholes location); source of data: the Polish Geological Institute-National Research. formations are very thick and are characterized by a very low hydraulic conductivity, they comprise an almost impermeable hydrogeological barrier that hinders the infiltration of groundwater from the nearby surface strata. The Quaternary-Upper Cretaceous aquifer is, therefore, isolated from the confined aquifers. As a result, no groundwater infiltration into the confined aquifer is observed from the Quaternary-Upper Cretaceous aquifer.

Figure 2.
Schematic cross-section of geological and hydrogeological formations in the study area (see Figure 1C for the location of the profile); modified after [69]; source of data: the Polish Geological Institute-National Research.

Damage to Rock Mass and Groundwater Extraction
The extraction of hard coal from the Bogdanka mine began in 1982. Hard coal was extracted from panels 2500 m in length and 350 m in width. By 2020, the extraction was concentrated mainly in the central and northwestern portions of the analyzed area. At that time, the mining operation was conducted in the depths of approximately 820 m to 1005 m from a deposit of 1.6 m to 3.5 m thick. In 2020, the area of the exploitation amounted to approximately 51 km 2 ( Figure 3A). According to Witkowski et al., the total value of the land subsidence resulting directly from hard coal mining was 5.5 m in the years 1982-2013 [70]. The future scenario of mining operations assumes that hard coal mining will also take place in the southern part of the analyzed area by the end of 2030. Mining operations will be conducted in the depths of approximately 900 m to 1000 m from a 1.3 m to 2.9 m thick deposit. The total area of exploitation will, therefore, be approximately 76 km 2 in 2030 ( Figure 3A).
The groundwater pumping started in June 1976 and is continuously carried out in four drainage wells located in the central part of the mining area ( Figure 3A,B). The mine drainage scheme implemented assumes that the boreholes are drilled from the roof of the mining operations and from the shafts to the Middle Jurassic aquifer. These boreholes aim to extract groundwater from the aquifer. It can therefore be inferred that the pumping wells are located directly in the Middle Jurassic aquifer and not at the mine excavation level. As expected, the groundwater from the Jurassic complex  Figure 1C for the location of the profile); modified after [69]; source of data: the Polish Geological Institute-National Research.
Four aquifers can be distinguished in the study area ( [69]; Figure 2). There are unconfined Quaternary-Upper Cretaceous and confined Lower Cretaceous-Upper Jurassic, Middle Jurassic, and Carboniferous aquifers. The thickness of these aquifers varies between 2-17 m, 0-3 m, 45-102 m, and 2-9 m, respectively (Table 2). Therefore, the Middle Jurassic aquifer is the dominant one. It is made of loose sands and strongly karstic limestones. The aquifer is completely separated from the roof by a series of practically impermeable carbonate layers of the Upper Jurassic as well as very thick Middle and Lower Cretaceous beds. The thickness of the series in the entire study area exceeds 350 m, to a maximum of 680 m.
There are many backyard wells drilled in the Quaternary-Upper Cretaceous aquifer that provide drinking water for residents in the study area. The other aquifers are not exploited for domestic purposes. However, the groundwater pumped out of the Middle Jurassic aquifer is partly used for technological processes carried out during mining operations. Nevertheless, much of the groundwater extracted from this aquifer is discharged directly into the surface water system. Table 2. Lithology and thickness (m) of the aquifers in the study area in selected boreholes B1-B4 (see Figure 1 for boreholes location); source of data: the Polish Geological Institute-National Research. Pressure in the rock mass caused a fracture zone resulting from a mining operation to disappear completely at a depth of about 300 m [69]. With an increase in depth, a constant decrease in the average hydraulic conductivity of the aquifer system can be observed. As the Jurassic and Cretaceous formations are very thick and are characterized by a very low hydraulic conductivity, they comprise an almost impermeable hydrogeological barrier that hinders the infiltration of groundwater from the nearby surface strata. The Quaternary-Upper Cretaceous aquifer is, therefore, isolated from the confined aquifers. As a result, no groundwater infiltration into the confined aquifer is observed from the Quaternary-Upper Cretaceous aquifer.

Damage to Rock Mass and Groundwater Extraction
The extraction of hard coal from the Bogdanka mine began in 1982. Hard coal was extracted from panels 2500 m in length and 350 m in width. By 2020, the extraction was concentrated mainly in the central and northwestern portions of the analyzed area. At that time, the mining operation was conducted in the depths of approximately 820 m to 1005 m from a deposit of 1.6 m to 3.5 m thick. In 2020, the area of the exploitation amounted to approximately 51 km 2 ( Figure 3A). According to Witkowski et al., the total value of the land subsidence resulting directly from hard coal mining was 5.5 m in the years 1982-2013 [70]. The future scenario of mining operations assumes that hard coal mining will also take place in the southern part of the analyzed area by the end of 2030. Mining operations will be conducted in the depths of approximately 900 m to 1000 m from a 1.3 m to 2.9 m thick deposit. The total area of exploitation will, therefore, be approximately 76 km 2 in 2030 ( Figure 3A).
The groundwater pumping started in June 1976 and is continuously carried out in four drainage wells located in the central part of the mining area ( Figure 3A,B). The mine drainage scheme implemented assumes that the boreholes are drilled from the roof of the mining operations and from the shafts to the Middle Jurassic aquifer. These boreholes aim to extract groundwater from the aquifer. It can therefore be inferred that the pumping wells are located directly in the Middle Jurassic aquifer and not at the mine excavation level. As expected, the groundwater from the Jurassic complex flows through the fracture zone to the mine and is drained out of the workings through the shafts. However, the volume of such groundwater is relatively small in comparison to the direct drainage of the Middle Jurassic complex.  Initially, only rock strata were drained near the drilled mine shafts. As the opening-up work progressed and hard coal extraction developed, the aquifer drainage was extended to the mineinduced area and the cracked zones overlying mining panels. The total amount of groundwater pumped out in the years 1976-2016 exceeded 180 million m 3 ( Figure 4). However, the average inflow of groundwater into the underground mining voids varied over time. In the years 1976-1978, when the shafts were drilled, there was a total of about 2800 m 3 /day. In the period 1978-1990, when the deposit was intensely cut with the opening-out, the groundwater inflow rapidly increased to about 16,000 m 3 /day in 1990. In successive years this value has slightly decreased and is now maintained at a constant level of approximately 14,500-16,000 m 3 /day. The observation network was established mostly before groundwater pumping to monitor changes in groundwater resulting from aquifer system drainage induced by underground mining Initially, only rock strata were drained near the drilled mine shafts. As the opening-up work progressed and hard coal extraction developed, the aquifer drainage was extended to the mine-induced area and the cracked zones overlying mining panels. The total amount of groundwater pumped out in the years 1976-2016 exceeded 180 million m 3 ( Figure 4). However, the average inflow of groundwater into the underground mining voids varied over time. In the years 1976-1978, when the shafts were drilled, there was a total of about 2800 m 3 /day. In the period 1978-1990, when the deposit was intensely cut with the opening-out, the groundwater inflow rapidly increased to about 16,000 m 3 /day in 1990. In successive years this value has slightly decreased and is now maintained at a constant level of approximately 14,500-16,000 m 3 /day.  Initially, only rock strata were drained near the drilled mine shafts. As the opening-up work progressed and hard coal extraction developed, the aquifer drainage was extended to the mineinduced area and the cracked zones overlying mining panels. The total amount of groundwater pumped out in the years 1976-2016 exceeded 180 million m 3 ( Figure 4). However, the average inflow of groundwater into the underground mining voids varied over time. In the years 1976-1978, when the shafts were drilled, there was a total of about 2800 m 3 /day. In the period 1978-1990, when the deposit was intensely cut with the opening-out, the groundwater inflow rapidly increased to about 16,000 m 3 /day in 1990. In successive years this value has slightly decreased and is now maintained at a constant level of approximately 14,500-16,000 m 3 /day. The observation network was established mostly before groundwater pumping to monitor changes in groundwater resulting from aquifer system drainage induced by underground mining  The observation network was established mostly before groundwater pumping to monitor changes in groundwater resulting from aquifer system drainage induced by underground mining exploitation ( Figures 3B and 5). It consists of 56 piezometers located in the Lower Cretaceous-Upper Jurassic and Middle Jurassic aquifers. Also, continuous monitoring of the groundwater table in the Quaternary-Upper Cretaceous aquifer is carried out in 34 wells distributed irregularly in the study area. The results of the observations revealed that the Lower Cretaceous-Jurassic and the Middle Jurassic aquifers are continually dewatered ( Figures 3B and 5). The Middle Jurassic aquifer, however, is exposed to the most intense drainage. The aquifer is permanently dewatered by mine shafts and boreholes drilled from the mining panels. Before groundwater pumping began, the initial groundwater head in the Lower Cretaceous-Upper Jurassic and Middle Jurassic aquifers varied between 140-170 masl and 130-150 masl, respectively [69]. The observed decrease in the Middle Jurassic aquifer reached a maximum of about 550 m, leading to a significant drop in the hydrostatic pressure. The depression cone has an irregular shape due to the geological configuration of the rock mass (see Section 2.1). It assumed the form of a large dewatering cone, extending 7.5 km southwest to approximately 20 km to the northeast of the drainage center ( Figure 3B). In the Lower Cretaceous-Upper Jurassic aquifer, a relati ely small decrease (up to 30 m) in the groundwater head is observed only in the mine shaft drainage center. This is the result of a rapid flow of quicksand into underground mine voids at the initial stage of the mining operation. The data on the change in the groundwater head in the Carboniferous aquifer are fragmentary, therefore the spatial extent and the value of the depression cone in this aquifer are only assumed.  Figures 3B and 5). The Middle Jurassic aquifer, however, is exposed to the most intense drainage. The aquifer is permanently dewatered by mine shafts and boreholes drilled from the mining panels. Before groundwater pumping began, the initial groundwater head in the Lower Cretaceous-Upper Jurassic and Middle Jurassic aquifers varied between 140-170 masl and 130-150 masl, respectively [69]. The observed decrease in the Middle Jurassic aquifer reached a maximum of about 550 m, leading to a significant drop in the hydrostatic pressure. The depression cone has an irregular shape due to the geological configuration of the rock mass (see Section 2.1). It assumed the form of a large dewatering cone, extending 7.5 km southwest to approximately 20 km to the northeast of the drainage center ( Figure 3B). In the Lower Cretaceous-Upper Jurassic aquifer, a relatively small decrease (up to 30 m) in the groundwater head is observed only in the mine shaft drainage center. This is the result of a rapid flow of quicksand into underground mine voids at the initial stage of the mining operation. The data on the change in the groundwater head in the Carboniferous aquifer are fragmentary, therefore the spatial extent and the value of the depression cone in this aquifer are only assumed.

Assessment of the Environmental Impact of Underground Mining
The methodology of the research presented consisted of three stages. The spatial extent of land subsidence resulting from underground mining was determined in the first part of the research work (Section 3.1). The maximum supposed spatial extent of the depression cone in the Middle Jurassic aquifer was calculated in the second part of the research work (Section 3.2). Subsequently, the computed spatial extent of the depression cone was used to determine the boundaries of the hydrogeological model. Finally, groundwater head change modeling in the Middle Jurassic aquifer was carried out in the third, main part of the research (Section 4).
Calculations of the spatial extent of the land subsidence were made for the period from the beginning of the mining operation, i.e., 1982 to 2030, as the exploitation is already planned for the next 10 years to be carried out. The maximum supposed spatial extent of depression cone in the Middle Jurassic aquifer was therefore also determined for the year 2030. Given the above, the numerical modeling of the changes in the groundwater head for the period 1975-2020 was carried out. Based on the calibrated model, the prediction of the groundwater head in the aquifer in the years 2020-2030 was also made.

Determination of the Spatial Extent of Land Subsidence
The determination of the spatial extent of land subsidence induced by underground hard coal exploitation was carried out with the use of the stochastic model based on the Knothe influence function [21,71,72]. Models based on influence function use a limited number of parameters. Therefore, these methods can serve as a comprehensive tool, especially when describing phenomena that require a significant degree of parameterization [73]. Many studies have shown that the application of influence function-based methods might indicate similar accuracy in predicting the spatial extent and values of land subsidence as the results of numerical modeling based on the most commonly used theoretical models [74].
The model based on the Knothe influence function is based on the division of the geological structure into elementary cuboids, i.e., reservoir elements [22]. In such models, the exploitation of the elementary cuboids and the horizontal plane generates the elementary subsidence at point A, expressed in the following Equation (1): where dS A is the elementary subsidence at point A; f is the influence function; x, y are the horizontal coordinates of point A and dP is the elementary area of the reservoir element. Implementing the principle of superposition, it has been assumed that the subsidence at point A is the sum of the elementary subsidence derived from all the elementary volumes in the area of Exploitation (2): where S A is the subsidence at point A.
Assuming the limited range of influence of the exploitation and introducing of the parameter R, the generalized influence function in Knothe model has the form of Equation (3): where S max is the maximum subsidence and R is the radius of the influence range, which take the form of Equation (4): where H is the depth of exploitation and tgβ is the parameter corresponding to geomechanical properties of the rock mass. The spatial extent of land subsidence in the hard coal mine of Bogdanka was calculated based on Equation (4). To assess the past impact of the mining operations, as well as predict the future exploitation scenario, the calculations were carried out between 1982 and 2030. The depth of the exploitation values was obtained from the geodetic measurements carried out and made available by the Mining and Surveying Department of the mine. Due to the inability to conduct field research and to obtain information on geomechanical rock mass parameters in the study area, the value of the tgβ parameter was obtained based on research conducted by [38,70]. It was set to 1.6.

Estimation of the Supposed Spatial Extent of the Depression Cone
The radius of influence is defined as the maximum distance at which groundwater drawdown can be detected with the usual field measuring devices. The size of the depression cone is influenced by the hydrogeological conditions of the aquifer system, the volume of the pumped groundwater, and the pumping time. Although the spatial extent of the depression cone is currently determined primarily based on numerical modeling, to quickly estimate its spatial extent, simple empirical formulas can be used [75].
In the research presented, the supposed spatial extent of the depression cone in the Middle Jurassic aquifer was determined with the use of empirical Weber Equation (5): where h 0 is the thickness of the confined aquifer; k is the hydraulic conductivity; t is time; and n is the porosity. The values of the parameters used in the calculations were determined based on field measurements carried out by the Mine and Surveying Department of the Bogdanka hard coal mine as well as literature data [69].

Mathematical Model
Groundwater flow in a porous medium is governed by two laws: the conservation of mass and the conservation of momentum [75,76]. These laws are expressed in two basic physical quantities: the hydraulic head (6): where p = p(x, y, z) is the hydrostatic pressure at a point with the coordinates x, y, z; ρ = ρ(x, y, z) is the density of the water and g is the gravity acceleration; and the specific discharge, defined as the volume of groundwater flowing in 1 s through 1 m 2 of the surface perpendicular to the direction of flow. It takes the form of Equation (7): where q x , q y , q z is the specific discharge in the x, y, z coordinates direction. Specific discharge is directly related to the superficial flow velocity (8): Sustainability 2020, 12, 7871 10 of 28 The law of conversation of groundwater mass in a porous medium expresses the continuity equation. It can take the form of Equation (9): The left side of Equation (9) represents the rate of accumulation of groundwater per unit volume of the porous medium, whereas the right-hand side of this equation expresses the excess volume of groundwater flowing into the porous medium per unit of time concerning the unit volume of the porous medium. Noticeably, the accumulation of groundwater in a saturated 3D porous medium is always associated with an increase in hydrostatic pressure. Empirically, the change in hydrostatic pressure induced by groundwater accumulation was found to be proportional to the rate of groundwater accumulation per unit volume of the porous medium [76,77]. It is expressed in Equation (10): where S s is the specific storage.
If the water density is constant, the 3D groundwater flow continuity equation most often assumes the following form (11): The law of conservation of momentum of groundwater in the porous medium takes the form of Darcy law. According to this law, there is a linear relationship at each point in the x = (x, y, z) space between the specific discharge and the force moving the groundwater to flow per unit volume of the medium. The Darcy law has the form (12): where κ is the permeability of a medium and µ is the dynamic viscosity of water. If the water density is constant, the right-hand side of the Equation (12) can be expressed by the hydraulic head. Thus, the Darcy law takes a simple and most commonly used form of Equation (13): where k is the hydraulic conductivity expressed in the Equation (14): The equation describing the transient 3-D flow in a saturated porous medium can be derived by substituting the Darcy law for the continuity Equation (15): The groundwater flow equation in the 3D saturated porous medium is, therefore, expressed in the form of (16): The determination of the solution Φ(x, t) = Φ(x, y, z, t) of Equation (16) is the starting point for all quantitative and qualitative analyses of hydrogeological changes, in particular those considered to be man-made induced impacts on the aquifer system. The dependence of the hydraulic conductivity on the spatial coordinates x, y, z, k = k(x, y, z) in Equation (16) enables this parameter to be differentiated in aquifers. The groundwater flow equation can, therefore, be applied successfully to heterogeneous aquifers.
In general, it is not possible to find an unequivocal solution Φ(x, t) to Equation (16). For this reason, approximate methods are used to solve the groundwater flow equation. The necessary and sufficient condition of the Φ(x, t). solution is that it simultaneously meets the following conditions: • Initial conditions Φ(x, 0) = Φ 0 (x) characterizing the spatial variation of the hydraulic head within the entire groundwater flow area Ω at a time t = 0 using known, assumed, or estimated function; • Boundary conditions describing the variation in the time of the hydraulic head or its derivative on all segments of the edge of the groundwater flow area ∂Ω.
The groundwater flow in the specific discharge is defined by differential equations using continuous variables. To build a numerical model, these equations must be transformed into a discreet form [78,79]. The discreet form of equations means that their approximate solutions are determined only at discrete points in space and at discrete moments in time.
In the research presented, the solution of the groundwater flow equation was carried out using the finite difference method. It is based on the approximation of the derivatives in the groundwater flow equation employing the differential quotients of the values of the desired function in the nodes of the discretization grid. For the first and second derivatives of function f concerning x, the simplest discrete approximation has the form of Equations (17) and (18), respectively: where f i+1 , f i−1 , f i are the desired values of the function in the discreet nodes of the discretization grid. Replacing derivatives with difference quotients (17) and (18) leads to a system of algebraic equations for N unknowns ( f i , . . . , f N ), where N is the number of discretization nodes. Usually, the values of the sought function are defined in the mesh nodes and the parameters in the respective equations are defined in the spaces between the nodes. An approximate solution is obtained by solving N algebraic equations with N unknown for discreet moments in time. An appropriate discreet scheme for variable t enables the transition from solutions obtained for a certain time to solutions for a subsequent time. As a result, the spatial-temporal characteristics of the approximate solution are achieved. If approximate solutions meet the desired accuracy criteria, the results obtained can be the basis for hydrogeological analysis.
A broader mathematical description of the groundwater flow theory and numerical methods used in hydrogeology is provided by [75,[78][79][80].

Numerical Modeling
IMOD software, the OpenSource GUI version of MODFLOW [79], was used to simulate the change in the groundwater head in the study area [81]. IMOD is a simulation system that models groundwater flow with rapid, flexible, and consistent sub-domain modeling techniques that implement the modular 3D finite-difference groundwater models MODFLOW 2005 and MODFLOW 6 of the U.S. Geological Survey. Application of MODFLOW to the description and prediction of the hydrogeological conditions of the aquifer systems has increased significantly and covers various types of research problems that have occurred worldwide [5,41,42,[82][83][84][85][86][87][88].

Numerical Representation of the Aquifer System
•

Model Conceptualization
The construction of a numerical hydrogeological model that would allow the assessment of the past and current impact, as well as predict future changes in the aquifer system in the conditions of underground hard coal mining, required many simplifications, reflecting the existing hydrogeological conditions. These simplifications are the result of the lack of sufficient recognition of the hydrogeological conditions and the complexity of these conditions, as well as the lack of need for their accurate mapping due to the modeling objectives assumed.
The hydrogeological conceptualization assumed the following: the Quaternary-Upper Cretaceous aquifer is completely isolated from the underlying aquifers by a thick series of almost impermeable Cretaceous strata; mining-induced drainage includes Lower Cretaceous-Upper Jurassic, Middle Jurassic, and Carboniferous aquifers; the main aquifer undergoing mining-induced drainage is the Middle Jurassic aquifer. Due to the low thickness of the Lower Cretaceous-Upper Jurassic and Carboniferous aquifers (see Table 1 and Figure 3), the inflow of groundwater into the mine voids from these aquifers was considered to be insignificant and therefore negligible; the crack zone occurred in the area of mining operations. Due to the presence and spatial propagation of the zone, the infiltration of groundwater from the Middle Jurassic and the Lower Cretaceous-Upper Jurassic aquifers increased; there is a direct hydraulic conductivity between the Lower Cretaceous-Upper Jurassic and the Middle Jurassic aquifers at the bottom of the Lower Cretaceous-Upper Jurassic aquifer, the Namurian complex outcrops and the crack zone over the hard coal mining area.
Given the foregoing, the following hydrogeological scheme for the research area was implemented: it was decided to map the Lower Cretaceous-Upper Jurassic and the Middle Jurassic aquifers as one groundwater complex; due to the hydrogeological complexity of the Quaternary-Upper Cretaceous aquifer, its isolation from the underlying aquifers and, therefore, no impact of mine-induced drainage on that aquifer, the Quaternary-Upper Cretaceous aquifer was not modeled in the research presented; regional groundwater inflow into the Jurassic aquifer was assumed to be not fully recognized due to a lack of empirical data; the surface water infiltration was mapped by the recharge of the top of the Cretaceous aquitard. •

Model Area and Finite Element Mesh
The model consisted of three layers: Cretaceous aquitard at the top (aquitard 1), Jurassic aquifer at the middle (aquifer 1), and Carboniferous aquitard at the bottom (aquifer 3; Figure 6). The model boundaries were determined based on the supposed maximum spatial extent of the depression cone in the Middle Jurassic aquifer in 2030 (see Section 3, Figure 7). Consequently, the impact of mining-induced drainage did not influence the model boundaries. The model covered the entire mining area, whereas all three layers covered the entire modeling area. The model was of regular shape with a dimension of 60 km by 60 km. It covered an area of 3600 km 2 and was divided into a regular mesh of 240 rows and 240 columns. The size of the cell was 250 m by 250 m, whereas the total number of cells was 57,600.

•
Geological and Initial Groundwater Head Input Data The elevation data of the top and the bottom of model layers were determined based on the results of laboratory tests of the geological cores carried out by the Polish Geological Institute-National Research Institute. They were obtained from 532 geological profiles which were distributed almost-regularly across the entire study area (Figure 1). Elevation maps of the model layers were generated by kriging interpolation. Due to the low value of Carboniferous stratigraphic throws (Figure 1), they were not included in the interpolation. Simultaneously, the initial groundwater head in aquifer 1 was obtained by interpolating the observational data from 41 piezometers located in the Middle Jurassic aquifer (see Figure 3 for the location of piezometers) and literature data [69] with the use of kriging. Sustainability 2020, 12, x FOR PEER REVIEW 13 of 28 According to the assumed conceptualization scheme, the infiltration of surface water at the top of aquitard 1 was assumed. The rate of the infiltration was assumed based on the literature data provided by Bocheńska et al. [69] and estimated as a single value of approximately 1.4 mm/year. The bottom of aquitard 2 is the impermeable bedrock. Hence, it was regarded as the no-flow boundary. Therefore, the boundaries of the layer were treated as Dirichlet boundary. •

Hydraulic Parameters
The hydraulic properties of the model layers were determined by horizontal hydraulic conductivity , vertical horizontal conductivity , and storage coefficient .
According to the literature data, several zones with different horizontal hydraulic conductivity values were distinguished for aquifer 1 and aquitard 2 ([69]; Figure 7). Due to the determination of these parameters only in parts of the model (hard coal deposits), the data were extrapolated to all blocks of layers 2 and 3 using the kriging method. For aquitard 1, due to data scarcity, a single horizontal hydraulic conductivity value was assigned. The values of the vertical hydraulic conductivity of all model layers were estimated as 1-2 orders of magnitude smaller than the values of the horizontal hydraulic conductivity [75].
The storage coefficient of aquifer 1 was estimated by comparing the volume of groundwater pumped from the Middle Jurassic aquifer with the observed values of the depression cone (19): where is the volume of groundwater released from storage; ∆ℎ is the unit drop in the groundwater head; and is the unit area of the hydrostratigraphic unit. The storage coefficient of aquitards 1 and 2 was estimated based on literature data using Equation (20): where is the thickness of the hydrostratigraphic unit; is the vertical skeletal compressibility; and is the compressibility of water. •

Boundary Conditions
The Dirichlet boundary conditions for the aquifer 1 were assumed. Thus, groundwater was simulated as moving in or out of the aquifer at a rate sufficient to maintain the specified head. It was assumed that the values of the groundwater head in the aquifer beyond the depression cone are constant and correspond to the initial values.
Mining-induced drainage of the aquifer system was modeled by setting four drainage wells in aquifer 1. The volume of pumped groundwater was determined based on data provided by the Mine and Surveying Department of the Bogdanka hard coal mine (see Section 2.2 and Figure 4). Due to data scarcity describing the volume of groundwater pumped after 2016, a constant drainage value of 15,000 m 3 /day was assumed in 2017-2030. The extraction values of the individual wells were calculated by comparing the total volume of water pumped out of the rock mass with the values of the depression cone in the Middle Jurassic aquifer at a given time.
According to the assumed conceptualization scheme, the infiltration of surface water at the top of aquitard 1 was assumed. The rate of the infiltration was assumed based on the literature data provided by Bocheńska et al. [69] and estimated as a single value of approximately 1.4 mm/year. The bottom of aquitard 2 is the impermeable bedrock. Hence, it was regarded as the no-flow boundary. Therefore, the boundaries of the layer were treated as Dirichlet boundary. •

Hydraulic Parameters
The hydraulic properties of the model layers were determined by horizontal hydraulic conductivity k h , vertical horizontal conductivity k v , and storage coefficient S.
According to the literature data, several zones with different horizontal hydraulic conductivity values were distinguished for aquifer 1 and aquitard 2 ([69]; Figure 7). Due to the determination of these parameters only in parts of the model (hard coal deposits), the data were extrapolated to all blocks of layers 2 and 3 using the kriging method. For aquitard 1, due to data scarcity, a single horizontal hydraulic conductivity value was assigned. The values of the vertical hydraulic conductivity of all model layers were estimated as 1-2 orders of magnitude smaller than the values of the horizontal hydraulic conductivity [75].
The storage coefficient of aquifer 1 was estimated by comparing the volume of groundwater pumped from the Middle Jurassic aquifer with the observed values of the depression cone (19): where V w is the volume of groundwater released from storage; ∆h is the unit drop in the groundwater head; and S a is the unit area of the hydrostratigraphic unit. The storage coefficient of aquitards 1 and 2 was estimated based on literature data using Equation (20): where b is the thickness of the hydrostratigraphic unit; α is the vertical skeletal compressibility; and β w is the compressibility of water.

Model Calibration, Verification, and Future Scenario
Model calibration and validation were achieved using a manual test-and-error procedure by matching the observed groundwater head with the simulated values to determine the set of best-fit parameters. The hydraulic head was fitted by varying horizontal and vertical hydraulic conductivity and storage coefficient. The calibration of each model should first set a target with an acceptable error.

Model Calibration, Verification, and Future Scenario
Model calibration and validation were achieved using a manual test-and-error procedure by matching the observed groundwater head with the simulated values to determine the set of best-fit parameters. The hydraulic head was fitted by varying horizontal and vertical hydraulic conductivity and storage coefficient. The calibration of each model should first set a target with an acceptable error. The range of the error depends mainly on the purpose of the model. During calibration, the user should focus on parameters that have been determined with less accuracy or are assumed and only slightly modify those parameters that are more certain. Most other parameters are less sensitive and need only be changed within a specific realistic range [89]. Therefore, during the calibration of the model presented, we tried to distinguish mostly the parameters assigned to aquitard 1 and 3. These layers were relatively poorly recognized. Firstly, the model was calibrated under unsteady conditions. Such calibration typically involves groundwater heads recorded in piezometers during long-term aquifer dewatering. Storage properties of the aquifer are critical for transient calibration. For this purpose, the groundwater head data recorded in 30 piezometers located in the Middle Jurassic aquifer were used from January 1976 to January 1995. This interval was divided into 20 time steps with a length of 1 year. The groundwater head estimated for the period before groundwater pumping and hard coal mining started was used as the initial condition. Numerous executions of the model were carried out to obtain optimal hydrogeological parameters values (approximately 170 executions). There was a good agreement between the observed and the calculated groundwater heads. The final values of the hydromechanical parameters selected for the model layer simulation are summarized in Table 3. Secondly, the model was verified for the period from January 1995 to January 2020. Therefore, for model validation, 26 time steps were used and each time step covered one year. The groundwater head, estimated for January 1995 during the calibration process, was used as the initial condition. Similarly, the groundwater heads observed in 34 piezometers in the Middle Jurassic aquifer were used for model verification. The calculated and observed groundwater heads in piezometers showed similar trends in general, and the calculated heads were mostly consistent with the observed data. However, a quantification of the model error for the different statistical parameters was carried out. These include the root mean square error (RMSE) and the determination coefficient (R 2 ; [90,91]). The results of the study were RMSE = 0.59 and R 2 = 0.64. Due to data scarcity and complex hydrogeological conditions, the accuracy of the model was not very high. These values were, however, assumed to be acceptable and reliable [91].
To illustrate the results of the model calibration and validation, two piezometers were selected for which quasi-continuous groundwater head data were available. These piezometers were located near the drainage center and in the outlying area of the depression cone (see Figure 3 for piezometer location). The time-series presented show acceptable agreement between the modeling results and the empirical data (

Spatio-Temporal Extent of Land Subsidence Induced by Underground Mining
The calculated values of the radius of the influence range, based on the Knothe theory, varied between 521 and 625 m, which is related to the variable depth at which hard coal is exploited. The total area of direct effects of underground mining in the study area was therefore approximately 16, 33, 53, and 82 km 2 , respectively, in 1990, 2000, 2010, and 2020 ( Figure 9A-D). The expected spatial extent of land subsidence will reach approximately 86 km 2 in 2030 ( Figure 9E). The spatial extent of the land subsidence resulting from mining in the successive years is consistent with the current boundaries of the mining area. This indicates that the direct effects of mining activities are restricted only to the immediate vicinity of mining shafts and mining panels. This is consistent with theories explaining land surface movements as a result of underground mining [7].
Due to the lack of access to more detailed data, computation of land subsidence resulting directly from mining activities was not introduced in the article presented. However, based on the precision leveling data provided by the Mine and Surveying Department and the study of such measurement data carried out by Witkowski et al. [92], it can be stated that the total value of the land subsidence resulting from hard coal mining was 5.5 m in the years 1982-2013. Based on this work, it can be concluded that land subsidence has occurred within the limits of the current mining area. It should be noted that this area is of agricultural use mostly ( Figure 9F). It is covered by a dense network of farmland, meadows, pastures, and forests. This area covers a few small villages with a population of not more than a few hundred inhabitants. The only town-Puchaczów-is located on the border of the mining area. The mining area includes over a dozen industrial buildings, as well as an industrial railroad used for the transport of hard coal.
Although the land subsidence values determined by Witkowski et al. are significant, owing to the relatively large area and the long period in which they were observed, no discontinuous deformations and substantial compressive or tensile strains occurred in the study area. For this

Spatio-Temporal Extent of Land Subsidence Induced by Underground Mining
The calculated values of the radius of the influence range, based on the Knothe theory, varied between 521 and 625 m, which is related to the variable depth at which hard coal is exploited. The total area of direct effects of underground mining in the study area was therefore approximately 16, 33, 53, and 82 km 2 , respectively, in 1990, 2000, 2010, and 2020 ( Figure 9A-D). The expected spatial extent of land subsidence will reach approximately 86 km 2 in 2030 ( Figure 9E). The spatial extent of the land subsidence resulting from mining in the successive years is consistent with the current boundaries of the mining area. This indicates that the direct effects of mining activities are restricted only to the immediate vicinity of mining shafts and mining panels. This is consistent with theories explaining land surface movements as a result of underground mining [7]. stability of the construction site. Therefore, the protective pillars are established as a rock mass where no mining is carried out. Due to the lack of access to information on the established mining pillars, we have not been able to include this issue in our research. The results obtained are, therefore, partially incorrect. This should be taken into account when the calculations are carried out in the case of an occurrence in the analyzed area of a building which is particularly sensitive to mining damage and therefore intended to be protected.

Aquifer System Drainage Resulting from Underground Mining
The results of the numerical modeling carried out in the study presented indicate that the Middle Jurassic Aquifer was heavily drained between 1976 and 2020 (Figures 8 and 10A). The maximum decrease in the groundwater head was observed at the drainage center and was 604 m. The depression cone has been developing most intensively in the north-east and south-east direction. As a result, a large drainage center related to drainage shafts was identified in 2020. The groundwater head in this area was approximately 400 mbsl. The spatial extent of the depression cone is Due to the lack of access to more detailed data, computation of land subsidence resulting directly from mining activities was not introduced in the article presented. However, based on the precision leveling data provided by the Mine and Surveying Department and the study of such measurement data carried out by Witkowski et al. [92], it can be stated that the total value of the land subsidence resulting from hard coal mining was 5.5 m in the years 1982-2013. Based on this work, it can be concluded that land subsidence has occurred within the limits of the current mining area. It should be noted that this area is of agricultural use mostly ( Figure 9F). It is covered by a dense network of farmland, meadows, pastures, and forests. This area covers a few small villages with a population of not more than a few hundred inhabitants. The only town-Puchaczów-is located on the border of the mining area. The mining area includes over a dozen industrial buildings, as well as an industrial railroad used for the transport of hard coal.
Although the land subsidence values determined by Witkowski et al. are significant, owing to the relatively large area and the long period in which they were observed, no discontinuous deformations and substantial compressive or tensile strains occurred in the study area. For this reason, due to the direct effects of hard coal mining, no significant threat to the surface infrastructure is expected. It should be noted, however, that the observed land subsidence caused a slight change in the watercourse profiles. However, this is related also to the extremely low land nature of the analyzed area. As a result, the local marshes occurred. These changes have been continually remedied by the mining entrepreneur. Nevertheless, owing to the further progress of the planned mining operation in the south-east by 2030, similar issues are expected to occur in the area under analysis.
It should be noted that the spatial extent of the land subsidence determined in subsequent years includes the mining shaft area ( Figure 9A-E). However, mining shafts are highly sensitive to the impacts of mining operations and are therefore subject to special protection. This protection shall be implemented by the establishment of the pillars of protection. The protection pillar is a rock mass area where mining is not expected to have any impact that could adversely affect the geomechanical stability of the construction site. Therefore, the protective pillars are established as a rock mass where no mining is carried out. Due to the lack of access to information on the established mining pillars, we have not been able to include this issue in our research. The results obtained are, therefore, partially incorrect. This should be taken into account when the calculations are carried out in the case of an occurrence in the analyzed area of a building which is particularly sensitive to mining damage and therefore intended to be protected.

Aquifer System Drainage Resulting from Underground Mining
The results of the numerical modeling carried out in the study presented indicate that the Middle Jurassic Aquifer was heavily drained between 1976 and 2020 (Figures 8 and 10A). The maximum decrease in the groundwater head was observed at the drainage center and was 604 m. The depression cone has been developing most intensively in the north-east and south-east direction. As a result, a large drainage center related to drainage shafts was identified in 2020. The groundwater head in this area was approximately 400 mbsl. The spatial extent of the depression cone is approximately 8 km to the south-west, 16 km to the north-east, 11 km to the north-west, and 14 km south-east of the drainage center. The depression cone covers an area of around 548 km 2 in 2020. To present the dynamics of the development of the depression cone in the Middle Jurassic aquifer, two piezometers were selected. In point J/5443, it can be noted that the greatest decrease in the groundwater head was observed only during the initial calculation period, in the years 1975-1995, and amounted to 230 m ( Figure 8A). From 1995 to 2020, the time series of the groundwater head showed a slighter decrease in the piezometric level. Such dynamics indicate practically complete depletion of the Middle Jurassic aquifer in the drainage center in 2020, i.e., approximately 45 years after the groundwater extraction began. The analysis of changes in the groundwater head in the J/5451 piezometer showed that the depression cone in this area had been developing intensively in the first years of the mining exploitation ( Figure 8B). The piezometric level decreased by 72 m between 1975 and 1995. In the period 1995-2020, the dynamics of the decrease in the groundwater head were much lower. Based on the modeling results, it can be concluded that the development of the depression cone has a strong time-based relationship with the exploitation of hard coal and, to a limited extent, with the geological structure of the deposit. In the area of the damage of rock layers caused by mining operations characterized by the highest hydraulic conductivity, the decrease in the groundwater head was the greatest (compare Figures 7 and 10). Furthermore, the cone of depression is characterized by a certain asymmetry. It is most likely related to the synclinical geological structure of the hard coal deposit and higher hydraulic conductivity in the Middle Jurassic aquifer to the north and north-east from the drainage center. However, the asymmetry of the depression cone can also arise from regional groundwater flow, which is not completely recognised due to lack of data. The depression To present the dynamics of the development of the depression cone in the Middle Jurassic aquifer, two piezometers were selected. In point J/5443, it can be noted that the greatest decrease in the groundwater head was observed only during the initial calculation period, in the years 1975-1995, and amounted to 230 m ( Figure 8A). From 1995 to 2020, the time series of the groundwater head showed a slighter decrease in the piezometric level. Such dynamics indicate practically complete depletion of the Middle Jurassic aquifer in the drainage center in 2020, i.e., approximately 45 years after the groundwater extraction began. The analysis of changes in the groundwater head in the J/5451 piezometer showed that the depression cone in this area had been developing intensively in the first years of the mining exploitation ( Figure 8B). The piezometric level decreased by 72 m between 1975 and 1995. In the period 1995-2020, the dynamics of the decrease in the groundwater head were much lower.
Based on the modeling results, it can be concluded that the development of the depression cone has a strong time-based relationship with the exploitation of hard coal and, to a limited extent, with the geological structure of the deposit. In the area of the damage of rock layers caused by mining operations characterized by the highest hydraulic conductivity, the decrease in the groundwater head was the greatest (compare Figures 7 and 10). Furthermore, the cone of depression is characterized by a certain asymmetry. It is most likely related to the synclinical geological structure of the hard coal deposit and higher hydraulic conductivity in the Middle Jurassic aquifer to the north and north-east from the drainage center. However, the asymmetry of the depression cone can also arise from regional groundwater flow, which is not completely recognised due to lack of data. The depression cone observed in 2020 and the geological structure of the study area suggest that the regional groundwater flow could be from the south-east to the north-west. However, this direction of water flow may be greatly influenced by the drainage of the main aquifer. Nonetheless, there are no comprehensive empirical data that could support the claim. Therefore, further research is needed to determine the impact of underground mining operations on the groundwater flow direction.
The results of the model prediction indicate that the depression cone will continue to develop along with the progress of hard coal mining in the years to come (Figures 8 and 10B). The maximum decrease in the groundwater head will be 644 m in 2030 in the area of drainage center. The depression cone will extend approximately 15 km north, 20 km north-east, and 12 km south of the drainage center. The total area of the depression cone will cover around 756 km 2 in 2030. It will, therefore, be about 208 km 2 greater than in 2020. The depression cone predicted for 2030 would be similar in shape to the currently observed decrease in the groundwater head in the Middle Jurassic aquifer.

Analysis and Discussion of the Conjugated Regional-Scale Impact of Underground Mining on the Environment
The spatial extent of the direct effects of underground coal mining is much less than the depression cone (see Figures 9 and 10C). In 2020, the area of the depression cone in the Middle Jurassic aquifer is about 6.7 times larger than the spatial extent of the land subsidence directly related to the mining operation. This factor will amount to approximately 8.8 in 2030. Significant land subsidence has occurred within the boundaries of the current mining area, resulting directly from the exploitation of hard coal and the formation of a post-mining void in the rock mass. However, the depression cone in the Middle Jurassic aquifer caused a substantial drop in hydrostatic pressure and is therefore probably related to the progressive drainage-induced compaction of the aquifer system. Although the mine-induced drainage does not affect the surface water and the Quaternary-Upper Cretaceous aquifer, the underlying aquifers have been intensively dewatered. The decrease in the groundwater head modeled for the years 2020 and 2030 is the greatest near the drainage center (Figures 11 and 12).
Based on the research presented, it can be concluded that the area of approximately 2.3 km 2 near the mine shafts was significantly drained by 2020 ( Figures 11A and 12). In this region, the groundwater in the Middle Jurassic aquifer was up to 24 m below the top surface of the aquifer. The thickness of the aquifer in the area is approximately 123 m. It is therefore expected that the depression cone center will be further drained in the years to come as mining works progress. In 2030, the modeled groundwater head will be 482 mbsl, i.e., 66 m below the upper surface of the Middle Jurassic aquifer ( Figures 11B and 12). The total spatial extent of the groundwater head, which is less than the upper surface of the Middle Jurassic aquifer, will be approximately 5.6 km 2 .
depression cone in the Middle Jurassic aquifer caused a substantial drop in hydrostatic pressure and is therefore probably related to the progressive drainage-induced compaction of the aquifer system. Although the mine-induced drainage does not affect the surface water and the Quaternary-Upper Cretaceous aquifer, the underlying aquifers have been intensively dewatered. The decrease in the groundwater head modeled for the years 2020 and 2030 is the greatest near the drainage center (Figures 11 and 12).  As a result of such a significant drop in the groundwater level of the Middle Jurassic aquifer, part of the drainage center area has been completely drained. Subsequently, it should be assumed that the rock mass could be disintegrated. This phenomenon is associated with irreversible aquifer system compaction. It should be noted, however, that the overall size of this area is relatively small (Figures 11 and 12).
A number of studies have shown that most of the developed groundwater basins are subject to inelastic deformation that occurs over the years (elastic deformation is usually seasonal) [93][94][95][96]. This is due to over-extraction of groundwater. Reducing hydrostatic pressure in the pores and cracks of the aquifer system is inevitably accompanied by some deformation of the aquifer system. Almost all permanent subsidence is due to irreversible consolidation during the generally slow process of aquifer drainage [97]. The relationship between changes in ground-water heads and compaction of the aquifer system is based on the principle of effective stress proposed by Terzaghi [98]. Under this theory, when the support provided by the fluid pressure is reduced, such as when groundwater is lowered, the support previously provided by the hydrostatic pressure is transferred to the skeleton of the aquifer system which compresses to a certain degree. The maximum level of past stress of the skeletal of the aquifer is referred to as pre-consolidation stress. When the load on the aquitard skeleton exceeds the pre-consolidation stress, the aquitard skeleton may undergo major, permanent rearrangement, resulting in irreversible compaction. As drainage-induced compaction of the susceptible aquifer was not a key point of the article presented, for a more detailed discussion of the issue please refer to [16,35,36].
upper surface of the Middle Jurassic aquifer, will be approximately 5.6 km 2 .
As a result of such a significant drop in the groundwater level of the Middle Jurassic aquifer, part of the drainage center area has been completely drained. Subsequently, it should be assumed that the rock mass could be disintegrated. This phenomenon is associated with irreversible aquifer system compaction. It should be noted, however, that the overall size of this area is relatively small (Figures 11 and 12).  Figure 11 for the location of profiles); a comparison of the altitude of the groundwater head in the Middle Jurassic aquifer, the terrain surface, the upper and lower surfaces, as well as the terrain surface and initial groundwater head in the aquifer under analysis.
A number of studies have shown that most of the developed groundwater basins are subject to inelastic deformation that occurs over the years (elastic deformation is usually seasonal) [93][94][95][96]. This is due to over-extraction of groundwater. Reducing hydrostatic pressure in the pores and cracks of the aquifer system is inevitably accompanied by some deformation of the aquifer system. Almost all permanent subsidence is due to irreversible consolidation during the generally slow process of aquifer drainage [97]. The relationship between changes in ground-water heads and compaction of the aquifer system is based on the principle of effective stress proposed by Terzaghi [98]. Under this  Figure 11 for the location of profiles); a comparison of the altitude of the groundwater head in the Middle Jurassic aquifer, the terrain surface, the upper and lower surfaces, as well as the terrain surface and initial groundwater head in the aquifer under analysis.
There has also been a decline in the groundwater head beyond the existing mining area. Therefore, land subsidence caused by aquifer drainage can also be assumed to occur in the entire analyzed region. This has been shown by the research carried out by Witkowski et al. [38], who predicted surface deformations due to groundwater extraction in the Bogdanka hard coal mine using a stochastic model. Based on the findings of this study, the maximum value of vertical displacement in the study area reached 0.56 m in the years 1976-2013. Nevertheless, the estimated value of the land subsidence resulting from the main aquifer drainage in the study area is not known in the following years.
It should be noted, however, that although the expected land subsidence values resulting from the drainage of the Middle Jurassic aquifer are much smaller than the land subsidence directly linked to hard coal mining, the spatial extent of the aquifer drainage is much wider. Therefore, it can also have a much wider effect on the environment. It should be noted that, under Polish law, the mining contractor is responsible for the negative effects of mining activities only in the mining area. Given that the decrease in the groundwater head extends significantly outside the mining area, it would be advisable to carry out more comprehensive research that would confirm and accurately determine the values of the drainage-induced land subsidence, especially outside the mining area established based on the direct impacts of mining workings.
Owing to the lack of hydraulic connections between the Lower Quaternary-Upper Createcous aquifer and the underlying aquifers, the drainage of the Middle Jurassic aquifer would not affect the depletion of the surface watercourses. However, due to the minor variations in land elevation (see Section 2), the drainage-induced land subsidence, which may occur outside the mining area, could in the years to come interrupt the surface runoff of the watercourses and alter the lake boundary. Noticeably, the depression cone modeled for 2020 and 2030 covers many protected areas ( Figure 12). These include 'Poleski' National Park and three landscape parks: 'Poleski', 'Nadwieprzański', and 'Pojezierze Łęczyńskie'. The ecosystems which dominate the landscape of the parks comprise swamps and peat-bogs. Such ecosystems are considered very delicate and can easily be influenced by outside factors. Nonetheless, there is no threat of draining water from protected areas due to the complete lack of hydraulic connections between the Middle Jurassic aquifer and surface water. However, the drainage-induced compaction of the aquifer system can lead to an increase of floodplains in protected areas, thereby enriching the local ecosystem [99]. Nevertheless, it is advisable to carry out more comprehensive research in the field of assessing the overall effect of rock mass drainage and the aquifer system compaction on the natural environment, including the land surface.

Concluding Remarks
Underground mining results in major and long-term changes in the conditions of the water environment. These changes could also be irreversible. It is therefore important to properly predict the impact of mining operations on the aquifer system through model tests to mitigate damage to that ecosystem, improve rock mass drainage processes, and reduce groundwater-related hazards. Reliable representation of processes taking place in aquifers involves the selection of appropriate tools to allow accurate mapping of several complex processes, rapid integration with other numerical models, and adaptation of the software to individual needs. As a result, a comprehensive assessment of the environmental impact of underground mining is a complex research issue. It requires the identification of the direct and indirect effects of mining operations. These primarily include land subsidence and aquifer system drainage.
Many methods may be used to determine the spatial extent of land subsidence. The function of influence-based method is most frequently used. This method allows for obtaining satisfactory results of such calculations. Besides this, due to low parametrization, this method is relatively easy to implement. However, the modeling of groundwater head changes in underground mines remains a challenge. The main problems associated with the construction of hydrogeological models are linked to the determination at the considerable depth of the hydrogeological parameters of the aquifer system. The presence of structural discontinuities, rock strata dips, or local thickness fluctuations has a significant impact on groundwater flow. For this reason, the authors searched for methods that would optimize the modeling of the phenomenon being studied.
Our process-based approach, using the estimation of the spatial extent of land subsidence resulting from underground mining as well as the spatial extent and values of the depression cone in the Middle Jurassic aquifer, has enabled us to compare the direct and indirect effects of underground mining on the regional environment. This is an important step towards a thorough recognition of the combined effects of underground mining. The approach also facilitates the analysis of the dynamics of the groundwater head decrease in the Middle Jurassic aquifer and the corresponding development of the mining operations. The local hydrogeological model was built based on the integration of qualitative and quantitative data from different sources.
As mining continues to develop in the study area, groundwater exploitation is likely to continue in the coming decades. The results of the model calibration and forecasting of the groundwater head in the Middle Jurassic aquifer could be considered to be acceptable and at least consistent with the empirical data. Predicted values of the depression cone formed as a result of the post-mining drainage assumed the shape of a large dewatering trough. This was consistent with the expected decrease in the groundwater head in the area during the analysis period. Our modeling results indicate that the hydrogeological system has been transformed from undisturbed to human-impacted aquifer depletion due to an increase in groundwater extraction over the past almost 45 years. Depletion of the aquifer system has exceeded the limits of the mining area. Although the values of the drainage-induced land subsidence are generally much lower than the direct impact of mining operations, the spatial extent of the aquifer depletion may influence the aquifer system far away from the drainage center and the mining area.
However, there are also significant limitations to the approach presented. The obtained depression cone values may show very good congruence with the empirical data, as the study area was well recognized geologically. Detailed altitude maps of geological strata were generated based on a wide range of data. It should also be noted that the degree of reliability of the model calibration and the accuracy of the forecast analysis should be spatially diversified. The results of the research have the highest level of certainty close to the observation points. Notably, the hydrogeological model was parameterized based on unknown groundwater flow direction values and partly based on literature data not directly describing local hydrogeological conditions. However, the reliability of the results of the calculation can be fully determined based on the analysis of changes in the groundwater head, accurate analyses, and uncertainty assessments. Where sufficient hydrological and geological data are available, this modeling approach can be applied to other underground mining areas worldwide facing local and regional impacts of groundwater-related environmental impacts.
The results of the research presented contribute to a more reliable environmental impact assessment of the underground hard coal mining. They point out that the environmental impact is much wider than the direct land subsidence. Research results clearly show that the depression cone extends considerably beyond the mining area. The combined spatial scale of land subsidence and aquifer dewatering should therefore be included in a comprehensive analysis of the impact of underground mining on the environment. The 3D numerical model presented in this study has the potential to provide highly relevant predictions of the mining-induced drainage of the aquifer system to be computed. This type of analysis can be used to understand the regional-scale processes that occur in the aquifer system due to mine-induced drainage. It should be noted that the recognition of hydrogeological phenomena and their dependence on the local characteristics of the aquifer is a prerogative of proper assessment of the effects of hydrogeological hazards, rational planning of mineral extraction, and protection in underground mining areas.
With future research in mind, it is worth paying attention to the difficulties associated with distinguishing land subsidence caused directly by mining from drainage-induced land subsidence. This problem concerns, in particular, the areas above exploitation, near post-mining voids, where these influences overlap. In the case of drainage-induced land subsidence, InSAR and GRACE (Gravity Recovery and Climate Experiment) seem to be useful tools for detecting this type of phenomenon. They have been successfully used to determine the full spatial extent of the direct and indirect environmental effects groundwater pumping in many parts of the world. However, further research is still needed to better implement such measurement techniques in terms of deep, underground mining.