Ecological Impact Prediction of Groundwater Change in Phreatic Aquifer under Multi-Mining Conditions

Abstract

In aeolian sandy grass shoal catchment areas that rely heavily on groundwater, mining-induced geological deformation and aquifer drainage are likely to cause irreversible damage to natural groundwater systems and affect the original circulation of groundwater, thus threatening the ecological environment. This study aimed to predict the impact of groundwater level decline on vegetation growth in the Hailiutu River Basin (HRB), which is a coal-field area. Based on remote-sensing data, the land use/cover change was interpreted and analyzed, and the central areas of greensward land in the basin were determined. Subsequently, the correlation between groundwater depth and grassland distribution was analyzed. Then, the groundwater system under natural conditions was modeled using MODFLOW, and the groundwater flow field in 2029 was predicted by loading the generalized treatment of coal mine drainage water to the model. The change in groundwater depth caused by coal mining and its influence on the grassland were obtained. The results show that coal mining will decrease the groundwater depth, which would induce degradation risks in 4 of the original 34 aggregation centers of greensward land that originally depended on groundwater for growth in HRB because they exceeded the groundwater threshold. The prediction results show that the maximum settlement of groundwater level can reach 5 m in the northern (Yinpanhao), 6 m in the eastern (Dahaize), and 10 m in the southern (Balasu) region of HRB. Attention should be paid to vegetation degradation in areas where groundwater depth exceeds the minimum threshold for plant growth.

1. Introduction

As the transfer base of coal mining in China, the Inner Mongolia Autonomous Region and Shaanxi province are gaining importance [1,2]. Controlled by the burial characteristics, more than 90% of coal resources in China are extracted by underground mining [3]. However, high-intensity underground coal mining has brought about a series of water resource and ecological environmental problems [4,5]. Continuous development of underground coal mining would lead to the deformation and collapse of the overlying strata, resulting in fissures and subsidence, which will destroy the groundwater aquifer system [6], reduce the regional groundwater level, and affect the surface ecological environment. Unlike other regions, coal deposits in Inner Mongolia and Shaanxi are located in arid and semi-arid regions [7,8,9], which are areas of serious soil erosion in China, and their ecological environment is very fragile. Consequently, there is a major contradiction between groundwater loss through drainage by coal mining and the lack of water resources for vegetation growth. Therefore, it is meaningful to study the ecological changes caused by the decline of groundwater level due to coal mining.

Coal burial in China is characterized by strong geological reconstruction after coal formation due to which water inrush accidents occur frequently during coal mining [10]. Such accidents pose serious threats to the safety of miners and mine production equipment, hinder the development of coal industry and mining operations [11], and even bring a series of ecological and environmental problems. Therefore, it is necessary to study the groundwater system of coal mining areas. For example, Yu et al. [12] carried out field trials on coal mine geological engineering in the Bucun coal mine, China, based on theoretical calculation and numerical simulation, and proposed a combined controlling measure considering underground coal mining and the water environment. Andres and Jose [13] presented a conceptual and numerical model of two linked mines in Spain for the assessment of possible environmental risks following closure. Wu and Wang [14] presented a conceptual framework with multilayered groundwater flow systems for the characterization of water bursting and discharge in underground mines, based on the characteristics of regional geological conditions, and they analyzed mine water inrush and flooding events in the north China coal basin.

With the gradual improvement of governmental efforts and increase in people’s awareness on the protection of the ecological environment [15], the impact of coal mining on the ecological environment is receiving increasingly more attention [4]. For example, in order to evaluate related effects in mining areas and non-mining areas from an ecological viewpoint, Xiao et al. [16] constructed an ecosystem service measurement and evaluation index system in a typical high-intensity coal mining area using remote sensing and geographic information system (GIS) and obtained good application results. Aiming to investigate the contamination level, distribution, and possible sources of polycyclic aromatic hydrocarbons (PAHs) in six coal mines in Xuzhou, China, and assess their ecological risks, Chen et al. [17] determined the concentrations of 16 PAHs in 26 underground samples, including coal, mine water, and underground sludge. Quaranta et al. [18] studied the influence of coal slurry injection in two underground coal mines on the potential water quality of surface water and groundwater in the central Appalachians and predicted the standard excess of surface water and groundwater using a dilution analysis model.

Recent studies have defined three stages of impact (damage phase, post-damage phase, and restoration phase) based on the effect of mining-induced groundwater changes on vegetation conditions [19]. Xie et al. [20] established a hydrogeological conceptual model and a groundwater flow numerical model to predict mining-induced changes in the groundwater flow field over five years (30 April 2015 to 30 April 2020). Based on the vegetation index determined using remote-sensing data, Zhang and Wang [21] reported that groundwater depth has a prominent effect on the spatial variation of vegetation distribution in the Ordos Plateau of China.

The border between Inner Mongolia and Shaanxi feature mega-coal fields, which are slated for mining in the next decade [1,22]. However, these coal fields occur below the Mu Us Sandy Land, which has the most vulnerable environment in China [23,24]. Therefore, changes in the groundwater and natural environment induced by coal mining are bound to affect the sustainable development of the region. These effects are mainly manifested in the reduction in groundwater, deterioration of vegetation that relies heavily on groundwater for growth, and regional vegetation shifting toward more drought-tolerant species [25], which are not suitable for grazing. More seriously, when groundwater shortage exceeds the basic threshold for maintaining vegetation growth, desertification will occur, making the area unsuitable for human habitation. Therefore, when groundwater discharge by coal mining exceeds natural recharge, groundwater resources will inevitably decrease sharply, ultimately threatening the healthy development of regional agriculture, animal husbandry, and industry. In particular, the Hailiutu River Basin (HRB) is a typical representative of this area. HRB has multiple planned coal fields, and its ecological vegetation and agricultural development mainly depend on groundwater recharge [26]. The relationship between the groundwater and ecology of HRB and its surrounding areas has been widely studied. In areas with groundwater depth less than 10 m within the HRB, the mean and standard deviation of normalized difference vegetation index (NDVI) values have been found to decrease with increasing groundwater depth [27]. Mata-González et al. [28] reported that vegetation growth and plant diversity are closely related to groundwater depth, especially in arid and semi-arid regions, and areas with shallow groundwater have higher vegetation coverage and plant diversity. Accordingly, the decline of groundwater level caused by coal mining will inevitably affect the status of vegetation growth, especially in the aeolian sand grass shoal catchment area, which is heavily dependent on groundwater [29]. Previous studies have investigated and predicted the influence of coal mining on groundwater mainly by establishing groundwater models [30,31]. In these models, the drainage by coal mining is used as input [32]. Subsequently, future changes in the groundwater flow field are predicted after model calibration and verification [33]. However, the prediction of the impact of groundwater change on the ecological environment has not been comprehensively analyzed thus far. In particular, detailed information on the specific impact of model predictions on the corresponding spatial distribution of a certain vegetation type after water level changes is lacking.

To address this gap, this study predicted the grassland degradation risk associated with groundwater decline induced by coal mining. Selecting HRB as a representative of the coal fields of Inner Mongolia and Shaanxi, mining-induced changes in groundwater circulation and their effects on the ecological environment were investigated. The groundwater depth of HRB was determined by measuring submersible wells, and the distribution of greensward land in the basin was interpreted using remote-sensing satellite images. Finally, changes in the groundwater flow field in the next 10 years (2020–2029), after the commencement of mining in first three mining areas in the basin, were simulated using a numerical model.

2. Materials and Methods

2.1. Study Area

HRB is located at the junction of Inner Mongolia and Shaanxi (Figure 1), which is an important coal mining base in northwest China. This basin belongs to the middle course of the Yellow River Basin. HRB is a small multiple bottomland basin in the transition area from the Ordos Plateau to the northern Shaanxi Loess Plateau [34]. The entire basin lies above the Mu Us Sandy Land, and the surface is mainly covered by Quaternary aeolian sand and lacustrine sand [35] (Figure 1). The basin covers an area of approximately 2600 km² and has a surface elevation ranging from 1020 m in the southeast to 1480 m above mean sea level in the northwest [36]. The perennial Hailiutu River flows from the northwest to the southeast of the basin and enters the Wuding River, the main tributary of the middle Yellow River. The Hanjiamao hydrological station is located at the outlet of watershed catchments, which had an average annual runoff of 0.8 × 10⁶ m³ from 1957 to 2014. HRB has a typical steppe environment with arid, semi-arid sandy land and beach land. It has a temperate continental monsoon climate, and atmospheric precipitation is the main source of water supply. The annual precipitation is 370 mm, and the annual average potential evaporation is 2000 mm [37]. The precipitation is mainly concentrated between July and September, and the wind direction is mainly northwest.

2.2. Geology and Hydrogeology

The study area has flat terrain, gentle strata, simple geological structure, and no distinct faults. The tectonic activity is not strong [38], and the watershed of HRB is composed of denudation loess ridges and high dunes [39]. Along at the river channel, the topography is steep, mainly due to river scouring and cutting of Cretaceous strata through flow erosion at the outlet of the basin. In addition to psammophytes, grass, and farmland, most of the land surface is covered by Quaternary sand dunes, with scattered spots of small ponds formed by groundwater. The Quaternary layer has various diagenetic features, but the thickness is generally small. Mesozoic and Cenozoic strata are widely exposed and generally water bearing, and the Lower Cretaceous Huanhe Formation of the Baoan Group is exposed only in some areas [40]. Within the basin, the Quaternary and Cretaceous aquifers exhibit distinct vertical variation, and their thickness varies widely.

The average thickness of the upper Quaternary loose rock soil is approximately 50 m. The thickness of the Cretaceous sandstone aquifers is approximately 350 m. Between the main 2^# coal seam and the Cretaceous floor is a special aquifer group formed by the alternate deposition of sandstone aquifers and impervious layers, which is the main source of mine water inrush [41]. The thickness of this special aquifer group is approximately 180 m. The Jurassic 2^# coal seam in the study area is the first mining coal seam with an average thickness of 5 m (Figure 2). Restricted by the paleotopography during deposition, the thickness of each stratum in the basin varies greatly [42]. Generally, the strata are the thickest in the paleochannel and low-lying center, gradually thinning to both sides, and nip-out at the watershed [43]. However, the thickness of Quaternary strata gradually decreases from downstream to upstream, and the thickness of some beach land in the middle of the basin increases. The thickness of Cretaceous strata gradually increases from downstream to upstream, forming a large aquifer. Among them, there is no continuous and effective aquifuge between Quaternary and Cretaceous aquifers, but there are multiple aquifuge and aquifer alternations in the Jurassic strata [44]. Therefore, the aquifer above Cretaceous strata is regarded as a water-rich layer, and the aquifer below is regarded as a weak aquifer. Based on the above, the Quaternary and Cretaceous strata in the study area can be regarded as a unified phreatic aquifer group. The Quaternary and Cretaceous aquifer groups have different thicknesses, but they generally exhibit a continuous distribution. Stratigraphic integration is relatively gentle among the Jurassic strata, and the 2^# coal seam is also evenly distributed in the basin because of the stress action of the upper strata [45].

2.3. Land Use/Cover Classification

Land use/cover change (LUCC) can provide crucial information for global environmental change and human environmental science [46]. LUCC can directly reflect the way humans utilize the land and the influence of human activities on the ecological environment [47]. In this study, we used Landsat images (TM/ETM+/OLI) of moderate resolution (30 m), remote-sensing data, and optical images from October 2019 acquired from the United States Geological Survey (USGS). We selected remote-sensing images with no or less cloud cover to improve the accuracy of land cover classification. For HRB, only one-scene images are required. Regarding the images, the cloud cover was 0.49%, and the Landsat path/row numbers were 128/33. After pretreatment, which involved radiometric calibration and atmospheric correction, land cover information was extracted using a supervised classification method [48]. In this study, land cover was classified into seven classes: greensward, psammophyte, water bodies, roads, buildings, farmland, and bare sand. The classification was performed according to the quantity and distribution characteristics in land use/cover types of the study area determined through field surveying and the research objectives of this study. We selected training samples from field investigations and Google Earth Engine for land use classification. A total of 287 samples were selected from evenly distributed sampling points in the watershed, including 58 greensward, 23 psmmophyte, 6 water bodies, 20 roads, 45 buildings, 68 farmlands, and 67 bare sandy lands. For pixel sample training, 70% of the samples were employed, and the remaining samples were used for testing (

Table 1. Description of land use/cover classification.

LUCC	Training	Testing	Total	Description
Greensward	144,330	61,855	206,185	Emergence of green
Psammophyte	857,465	367,485	1,224,950	No obvious geometric features, mostly dark green or brown in color
Water bodies	3331	1427	4758	Lakes and reservoirs, and some wide rivers
Roads	7146	3063	10,209	Usually near buildings, in a regular straight line
Buildings	45,126	19,340	64,466	Regular geometry exists, concentrated distribution, with red or blue roofs
Farmland	348,413	149,320	497,733	Emergence of green, with clear geometric boundaries
Bare sand	622,805	266,916	889,721	Yellow, with the feeling of sand flowing

).

The distribution of greensward area in HRB was processed in the ENVI software, which was used to validate the classification results. Although greensward area accounts for only 7.1% of the basin, its sensitivity to changes in groundwater depth is the highest among all land use types.

2.4. Interpolation of Groundwater Depth

Spatial interpolation is an effective tool for studying groundwater distribution and has been widely used in the field of groundwater resources [49]. As only limited point-scale groundwater depth information can be obtained through monitoring stations in HRB, we used spatial interpolation to process the measured 100 wells and obtain spatially continuous groundwater depth data. To ensure that the wells were spatially representative, we selected 21 wells for continuous monitoring from July to early September (Figure 3). The monitoring showed that the average variation of groundwater level was only 0.1 m. To gain a more realistic view of groundwater depth in the watershed, we used the measured groundwater depth as the main variable and some small lakes and river node locations in the watershed as auxiliary variables, and we used the Co-kriging method of the ArcGIS software for interpolation [50]. Co-kriging is an extension of ordinary kriging, in which the best estimation method of regionalized variables is developed from a single attribute to two or more synergistic regionalized attributes, and one or more auxiliary variables are applied. These auxiliary variables, which are interpolated and estimated, are related to the main variables, and the correlations between the variables are assumed to be used to improve the accuracy of the main prediction target [51]. The interpolation of groundwater depth is mainly based on the following equation:

$$Z^{\ast }\left(x_0\right)={\displaystyle \sum }_{i=1}^n{\lambda }_{1i}{Z}_1\left(x_i\right)+{\displaystyle \sum }_{j=1}^m{\lambda }_{2j}{Z}_2\left(x_j\right)$$

(1)

where x is variable location, $Z^{\ast }\left(x_0\right)$ is the predicted value for the x position, is the measured value of the primary variable, $Z_2\left(x_j\right)$ is the measured value of the covariate, λ is Co-kriging weights for Z, n is the number of the primary variable, and m is the number of the covariate.

2.5. Model Description and Setup

Groundwater was simulated to predict changes in groundwater level in the basin under the combined exploitation of multiple underground mines and to study the change in groundwater resources using the simulated water level. The circulation of the groundwater system was simulated using the MODFLOW groundwater model, based on Darcy’s law and the foundation of Dupuit, which are widely used in groundwater research [52,53]. MODFLOW was implemented in the GMS software platform to simulate the watershed. Through MODFLOW finite difference numerical simulation, HRB was subdivided into square meshes of 200 rows and 200 columns, on average. Finally, the watershed was divided into 20,986 grids. To establish the model, hydrological data, hydrogeological information, meteorological data, and groundwater level data were acquired from various sources. Geological drilling information and topographic elevation data are needed in the initial stage of the model. Hydrogeological parameters, water level, evaporation, and rainfall were applied as the driving forces of the hydrological cycle. For the initial groundwater flow field, important data were acquired through a groundwater survey from April to May 2019. Simulated and measured values were compared according to various evaluation metrics, and the model was optimized by fitting the actual values. The optimized model was finally used to predict the future conditions. In the groundwater simulation, the input factor values required adjustments within a certain range for adaptation to the ground conditions.

In order to establish an effective groundwater model, the aquifer system was generalized into three layers consisting of Quaternary aquifers, Cretaceous aquifers, and Jurassic strata above the 2^# coal seam considering the strata lithology, groundwater occurrence, and hydrodynamic characteristics of the study area. In addition, the boundary conditions and internal structure of the study area were generalized according to the actual relationship between water supply and discharge (Figure 4). The phreatic aquifer in the basin is directly related to the atmosphere and vertical water exchange through surface processes such as precipitation infiltration, evaporation and drainage, and farmland irrigation. Horizontal groundwater movement is determined by the topography and geological characteristics of the study area. Based on the accurate description of the conceptual hydrogeological model of groundwater system, a mathematical model of the basin was established. This three-dimensional and finite difference groundwater flow model is based on the groundwater flow equation, combined with the three-dimensional mathematical model of non-steady flow. The model was visualized in GMS. The final mathematical model of the basin is as follows:

$$\frac{\partial }{\partial x}\left(K_h\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_h\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_z\frac{\partial h}{\partial z}\right)+W={\mu }_s\frac{\partial h}{\partial t}$$

(2)

$${h|}_{t=0} = h_0(x,y,z) (x,y,z) \in \Omega$$

(3)

$$K\frac{\partial h}{\partial n}{|}_{\Gamma } = q(x,y,z,t) (x,y,z) \in \Omega$$

(4)

$$h\left(x,y,z,t\right){|}_{И} = h (x,y,z) \in \Omega$$

(5)

$$\frac{K_{r}A}{M_r}\left(H_r-h\right)=Q_r$$

(6)

where ${\mu }_s$ is the water storage coefficient, h is the groundwater head or elevation (m), $K_h$ is the hydraulic conductivity in the horizontal direction (m/d), $K_z$ is the hydraulic conductivity in the perpendicular direction (m/d), W is the strength of source and sink (m³/d), q(x,y,z,t) is the inflow or outflow function from a unit area in a unit time of the second boundary condition, Γ is the flow boundary, И is the head boundary, Ω is the scope of the model, $K_r$ is the river permeability coefficient (m/d), $M_r$ is the sedimentation thickness of the river channel (m), $H_r$ is the water level in the river (m), and $Q_r$ is the stream flow (m³/d).

According to the characteristics of aquifers in the watershed and the distribution of the groundwater flow field, shallow water aquifers in the watershed were divided into flow boundary and closed boundary aquifers. Along the streamflow direction of groundwater, the north of the watershed is the recharge boundary, and the south is the discharge boundary. The Hailiutu River receives water from continuous groundwater discharge, and it was thus simulated as a drain [21,36]. In order to characterize the discharge process of the river in detail, the river course was divided into 7 sections, and the actual flow was measured in each section. The measured average flow was input into the model. The phreatic aquifer in HRB is directly recharged by atmospheric precipitation. As the sandy land is characterized by a high infiltration rate, scarce vegetation on the surface, and low evaporation capacity, the free surface of the phreatic aquifer is regarded as the boundary surface of vertical recharge. In the model setting, precipitation directly recharges groundwater, and evaporation is affected by groundwater depth. Moreover, farmland irrigation wells are evenly distributed in the basin, which can be regarded as the surface discharge of groundwater resources. For the simulation of recharge and evaporation processes, the Recharge Package and Evapotranspiration Package of MODFLOW were implemented [54], with data obtained from the measured precipitation and evaporation at the meteorological stations of Yulin City near the basin. The groundwater flow conditions of all aquifers in the study area are essentially similar, but the hydrogeological parameters are spatially different. The hydrogeological parameters were mainly determined through pumping tests, but the initial hydrogeological parameters simulated by the model may not closely match the actual situation. This is because the parameters obtained by the pumping test are only point data, while the model requires area data. Considering the characteristics of the groundwater flow field and hydrogeology data collected during the geological survey, based on hydrogeological parameters, the phreatic aquifer was divided into 4 zones (Quaternary) and 3 zones (Cretaceous), and the confined aquifer (Jurassic) was divided into 4 zones (Figure 5). Groundwater in the study area is mainly recharged by atmospheric precipitation, and the infiltration coefficient of precipitation is determined by the surface geotechnical properties. The infiltration coefficients of precipitation are 0.34 (aeolian sand), 0.26 (lacustrine sand), 0.29 (valley sediment sand), 0.07 (eolian loess), and 0.09 (Cretaceous sandstone) for the five zones of the top layer (Figure 1c). The depth limit of phreatic evaporation was defined as 4 m according to preliminary knowledge of the investigator. Water quotas for farmland irrigation were recorded from actual observations.

2.6. Model Calibration and Validation

Through the accurate calibration of the measured and calculated water levels, the groundwater flow in HRB was successfully simulated. The actual observed groundwater level and flow field were mainly fitted by trial-and-error calibration, such that the simulated groundwater levels are reasonably coordinated with the observed values. The observed groundwater level data from July 2019 to December 2019 (184 days in total) were used for model calibration. To ensure that the dynamic change in the simulated groundwater level was as close as possible to the observed data, the hydrogeological parameters were constantly adjusted [55]. At the end of the simulation (184 days), no significant difference was observed between the simulated non-pressurized aquifer flow field and the measured groundwater flow field.

As shown in Figure 6 and Figure 7, the simulated groundwater level and the observed groundwater level were in line with the groundwater hydraulic properties of the basin. The groundwater levels of our long-term observation wells in Huhetala and Xinmiaotan are in good agreement with the simulated values. According to the centralized survey of the 21 selected wells at the beginning of each month from July to September, the water levels in the observation wells agree with the simulated values (Figure 7). Moreover, the maximum groundwater level error between the observed data and the simulated data of the two observation wells was 0.5 m, which is negligible compared to the natural variation of groundwater level. In the fitting process, hydrogeological parameters, as important data, were repeatedly adjusted manually to meet the actual situation, and the spatial distribution of hydrogeological parameters was finally determined. The hydrogeological parameters determined after model verification are shown in

Table 2. Hydrogeological parameters of aquifers.

Stratum Area		Hydraulic Conductivity (m/d)	Specific Yield	Storage Coefficient (m−1)
Q	Ⅰ	7	0.26	-
	Ⅱ	3.5	0.2	-
	Ⅲ	1.5	0.06	-
	Ⅵ	0.8	0.01	-
K1	Ⅰ	0.25	0.05	-
	Ⅱ	0.3	0.05	-
	Ⅲ	0.3	0.06	-
J2	Ⅰ	0.003	-	1.5 × 10−7
	Ⅱ	0.004	-	1.9 × 10−8
	Ⅲ	0.008	-	3.3 × 10−7
	Ⅵ	0.007	-	4.0 × 10−8

. In addition, as shown in

Table 3. Results of groundwater budget in the study area.

Balance Project		Annual Amount (×106 m3/a)
Recharge	Precipitation	361.83
	Lateral recharge	3.23
Discharge	Evaporation	197.58
	Lateral discharge	4.17
	River discharge	81.99
	Artificial exploitation	80.43
	-	0.89

, the discharge term in the basin was basically equal to the recharge term, indicating that the water cycle of HRB is in a dynamic equilibrium state, further indicating that the generalization of the water resources cycle in the study area is valid. Groundwater recharge mainly occurs through meteoric precipitation and lateral recharge, and drainage occurs mainly through surface evaporation, overflow to surface water, lateral discharge of groundwater, and artificial exploitation. The final verification results show that the groundwater model of HRB meets the accuracy requirements, which can reflect the hydraulic characteristics of the groundwater system in the basin and can be used to predict the change of groundwater level in the basin after coal mining.

2.7. Mine Inflow

In order to ensure the safety of underground mining, groundwater in the upper aquifer of the coal seam and water inrush will be drained, which will significantly affect the groundwater flow field of the overlying aquifer and the natural circulation of groundwater resources in the region. Using the established numerical model of groundwater flow, the relationship between groundwater level and drainage in the mining area can be quantitatively described and predicted, and the possible impact can be analyzed. For the prediction, coal mine drainage should be preliminarily input to the model, while other hydrological data remain unchanged. After coal seam mining, the deformation and failure of the overlying strata would generate caved, fractured, and continuous deformation zones [6,56], which is the main reason for the change in the groundwater system. Considering that the study area has thick coal seams for mining, the combined thickness of caved and fractured zones corresponds to the major drained aquifers. Through the field investigation of the study area and surrounding coal mines, the water conductive fracture zone of the coal mine was found to develop only in the area above the working face, and the strata outside the mining face remained largely unaffected. Therefore, the effective thickness of the drainage layer is defined as the distance between the top boundary of the fractured zone with water flow and the bottom boundary of the coal seam. Considering the current situation of deep buried coal seams in HRB [45], with coal mining, the largest water-conducting layer would include only the Jurassic formation, and it would be separated from the Cretaceous bottom by a certain distance. The distance would be 100.87 m in the Yinpanhao (Yph) coal mine, 94.22 m in the Dahaize (Dhz) coal mine, and 110.9 m in the Balasu (Bls) coal mine. Pumping wells for simulating mine inflow were expanded to the drainage layer previously defined on the basis of the first mining area. According to the prediction report of the collected mine water inflows, the drainage water of the three coal mines was 13,549.8 m³/d (Yph), 35,520 m³/d (Dhz), and 37,000.8 m³/d (Bls), and the planned mining volume of the three coal mines was 15 Mt/a (Yph), 15 Mt/a (Dhz), and 10 Mt/a (Bls). The drainage water of Bls was equal to that of Dhz, and the drainage water of Yph was the smallest. The drainage water increased gradually from upstream to downstream, according to the location of the coal mines in the basin.

3. Results and Discussion

3.1. Relationship between Vegetation and Groundwater Depth

Compared with the use of only groundwater depth as the interpolation factor, Co-kriging interpolation provides results with smoother spatial distribution without any sudden point aggregation characteristics, which can better reflect the spatial distribution of groundwater [57]. Therefore, during the field survey of 100 wells (Figure 3) in October 2019, the groundwater depth of the basin (Figure 8a) was acquired by Co-kriging interpolation. The results show that groundwater depth exhibits wide spatial variability across the study area. Specifically, groundwater depth is relatively large in the south and northwest regions, mainly distributed in the loess ridge and high dunes. On the basin scale, areas with suitable groundwater depth for vegetation are mainly distributed in the upper and middle reaches. According to the results of the land use type distribution (Figure 8b), vegetation coverage is relatively high in the areas with suitable groundwater depth, concentrated in the middle and upper reaches of the basin. Comparing Figure 8a,b, it can be ascertained that the location of shallow groundwater depth coincides with the dominant greensward area, which exhibits central aggregation points. Accordingly, the frequency graph of each pixel was plotted cumulatively, and the value of central aggregation points was determined based on the following conditions: greensward land use type with a resolution of 1500 m and dominant greensward pixels over other land use classifications. The dependency of greensward land on groundwater depth was investigated by organizing pair data of greensward land and groundwater depth on the same base map (Figure 8c). The diagram of the final results reveals 34 central aggregation points of greensward land, of which 24% are distributed in areas with groundwater depth less than 3 m and 11.17% in areas with groundwater depth more than 10 m.

The dominant spatial pattern of the distribution of greensward land in HRB is presented in Figure 8b. Greensward land clearly occupies a relatively independent niche in space, which is related to groundwater depth, and indirectly reflects its different dependences on different groundwater depth conditions. The pixel number of greensward land was counted from LUCC data, and 82% of the greensward area was found to be distributed in areas with groundwater depth of 0–4 m (Figure 9). From the relationship between the distribution of groundwater depth and the number of greensward pixels, the groundwater depth of 4 m can be ascertained to be an important threshold. In aeolian sandy soil, grass can grow normally only when the deep roots of grass have access to groundwater. If the groundwater depth exceeds the root depth, the survivability of grass significantly declines. For example, greensward land is unlikely to develop in sand dunes because the groundwater level exceeds 8 m. The results suggest that grass growth is seriously limited in areas where the groundwater depth exceeds 4 m. It is noteworthy that groundwater irrigation also determines the coverage of green space. Figure 8c also confirms this, but this aspect is not the focus of our study.

The relationship between groundwater and ecological vegetation has always been an important issue for research [58,59]. In particular, the relationship between vegetation and groundwater is very complex in arid and semi-arid areas, where plant roots have a strong response to changes in groundwater [60]. It is not sufficient to discuss the relationship between vegetation and groundwater depth by taking only a specific area of HRB as an example. To comprehensively present the characteristics of vegetation coverage in arid and semi-arid regions suitable for groundwater burial depth, we compared studies on vegetation and groundwater in similar regions to confirm the credibility of our conclusions [61,62,63,64,65,66]. As shown in

Table 4. Relationship between vegetation growth and groundwater level in different areas.

Research Objectives	Number of Studies	The Link between Vegetation and Groundwater Levels
The Owens Valley of California, USA (1991–2004)	30 permanent monitoring sites	About 3.5 m is the groundwater threshold for normal vegetation coverage in most areas
The riparian area in southeastern Arizona, USA (2003)	Three sites	The average annual groundwater depth of the grassland is 2.6 m
The lower section of the Tarim River, China (2005)	Nine 50 m wide transects	When the groundwater depth is −3.14 m, the density and diversity of above-ground plants are significantly improved
The hinterland of the Badain Jaran Desert, China (2014)	10 long-term observation sites	In the area where the groundwater depth is 0–2 m, the growth of vegetation is more vigorous
Al Qunfudah City and the surrounding coastal plain in southwest Saudi Arabia (2018)	Six locations	Vegetation coverage is higher in areas where the groundwater depth is less than 5 m
The Manasi River riparian zone, China (2018–2019)	13 sites	The water-table depth appropriate for herbs is 1–1.5 m

, some arid and semi-arid regions around the world have a limited range of groundwater levels for the healthy growth of vegetation. If the minimum threshold of groundwater required for vegetation growth is exceeded, the growth of vegetation will be affected. Attributable to the difference in hydrogeology, different study areas exhibit distinct differences in the threshold. Nevertheless, comparing the conclusions of previous studies, the threshold of groundwater burial depth for the normal growth of herbs was found to be approximately 2–5 m, which is essentially the same as the results obtained in this study. This can be explained by the limited root systems of herbs.

3.2. Groundwater Level Prediction

According to the groundwater model established after verification, aquifer drainage caused by coal seam mining was separated, and mine gushing water was set as the output of the model. Finally, the model was used to predict future changes in groundwater induced by coal mining. In the model, wells were used to simulate mine water inflow. A simulation was performed for each grid of the first mining area, and the inflow was then interpolated within the effective mining drainage layer of the Jurassic strata. The value of total mine water inflow of the grids should be close to the value of the previously mentioned prediction report. Therefore, to ensure accurate simulation results, all pumping wells were set in the mesh grids. These wells, representing mine water inflow, were negative and were evenly distributed in the first mining area of Yph, Dhz, and Bls. In the model, the first mining area of the three coal mines was assumed to satisfy the designed production volume in the future. By running the model, the groundwater flow field of the aquifers in the next 10 years (2020–2029) could be predicted. Figure 10a shows the simulated flow of the phreatic aquifer in December 2029.

The change of the flow field curve represents the influence of coal mine drainage. Therefore, the flow field curve of the phreatic aquifer above the first mining area of the coal mine is different from that of other areas of the basin. It can be observed that the maximum drawdown may be as high as 5 m (Yph), 6 m (Dhz), and 10 m (Bls), indicating three cones of depression centered on the three mines. The cone of depression exerts a distinct influence in Yph, but the change of flow field in Dhz and Bls coal mines is more complicated because of the influence of groundwater watershed on the hydraulic gradient. Overall, mining activities affected groundwater in phreatic aquifers within a certain range [2]. The size of groundwater depression is closely related not only to the amount of drainage by the coal mine but also to the geographical location of coal mines in the basin. An obvious drop of groundwater levels could be observed for floor leakages in the first mining area of the three coal mines (Figure 10b), and the groundwater depth shows no significant change in the upstream of HRB because it is located far away from coal mining sites. It is essentially certain that, as coal mining continues, groundwater depression will continue to expand, thereby affecting the hydrological cycle of the basin.

3.3. Vegetation Change Prediction

Given the single vegetation structure and sensitive environment of the study area, the ecological function of the vegetation community is susceptible to interference from climate change and human activities, which may even completely destroy the original vegetation cover [21]. Because HRB is the recharge source area of the Wuding River Basin, and precipitation is the main source of groundwater recharge in the basin, groundwater resources are important ecological factors. The greensward land in the basin is mainly supported by groundwater. Therefore, groundwater depth is the decisive factor controlling the distribution of greensward land. Grass coverage and occurrence frequency are related to groundwater depth. With the increase in groundwater depth, the coverage, density, and population richness of grass decrease. The growth and development of grass is higher in areas with suitable groundwater depth, with higher frequency and coverage. On the contrary, in water deficit areas, the growth and development of grass are relatively poor, with low occurrence frequency and corresponding coverage. The most suitable groundwater level for grass growth in HRB is 0–4 m, and this groundwater depth range is also the key guarantee to maintaining the sustainable development of the ecosystem in arid and semi-arid greensward areas [67,68]. If the groundwater level continues to decline below 4 m, the growth and development of greensward vegetation in pastoral areas will be limited, and the number of biological populations will gradually decrease. As the groundwater depth continues to drop below 7 m, arid greensward areas will undoubtedly undergo irreversible desertification. Therefore, the upper limit for the stable survival of the existing greensward vegetation can be determined to be 0 m and the lower limit to be 4 m.

Using MODFLOW to predict the groundwater flow field in the basin after 10 years, we obtained the groundwater depth of the future basin, and central aggregation points of greensward land over the spatial distribution of groundwater depth were determined (Figure 10c). Compared with Figure 8b, the groundwater level in the first mining area of Yph decreased significantly, and the groundwater depth increased. The change of the flow field above the first mining area causes the increase in groundwater depth near the river, and this area is the central aggregation point of greensward land in the basin. The increase in groundwater depth will seriously affect the normal growth of grass vegetation. Therefore, these regions require serious attention. As mentioned earlier, the groundwater depth of 4 m is an important threshold. However, in order to take necessary measures before reaching the warning line, attention should be paid to affected areas when the groundwater depth drops below 3 m. Accordingly, the central aggregation point of greensward land in this region is marked as red to alert the government and local population of the degradation risk. By red tagging, four greensward degradation risk areas were found to be concentrated in the area and surroundings of Yph in the midstream of the basin. As clearly shown in Figure 11, the predicted groundwater depth in this area after coal mining exceeds 3 m, increasing the proportion central aggregation points with groundwater depth greater than 3 m to 35%. Although the greensward coverage near Dhz and Bls is small, and there is no major central aggregation point of the greensward land, attention should also be paid to the impact of groundwater level decline on ecological vegetation diversity. Finally, the prediction of groundwater and grassland coverage under multi-mining conditions shows that the overall ecological prospect of the basin is negative without the protection of the ecological environment. From the final research results, the change of groundwater in the next 10 years will affect only the growth of vegetation within a certain range. However, the influence of groundwater depth change caused by coal mining on natural vegetation should be carefully examined at the watershed scale.

Although vegetation can adapt to changes in the environment, some vegetation types [69] require a certain amount of groundwater to maintain normal growth. With changes in groundwater, groundwater-sensitive plants may completely disappear from the corresponding area. Although the environmental crisis caused by coal mining is globally widespread [70], the coal industry is particularly important for regional economic development and national energy security in developing regions. Unless a breakthrough is achieved in the technology of renewable energy, and renewable energy is made affordable, coal mines will continue to operate, which will lead to the deterioration of the ecological environment after the drainage of groundwater [71]. Indeed, the best solution is to prevent mining and stop the destruction, but this solution is not feasible for developing countries and poor countries. We believe that coal mining activities in developed countries should be gradually discontinued, and clean energy measures should be increasingly implemented, such as in Germany [72]. For developing countries, the mining of coal resources should be limited by environmental protection. Under the condition of maintaining environmental sustainability, underground resources should be gradually exploited, and part of the profits from coal mining should be used for the management of hazards caused by coal mining, such as in China [73]. The contradiction between coal mining and environmental protection requires the preliminary intervention of government departments, for which the prediction of potential environmental changes caused by coal mining and subsequent active formulation of environmental protection policies are necessary.

4. Conclusions

In this study, variations in the flow field of phreatic groundwater caused by the joint action of multiple coal mines in a basin were explored, and the groundwater flow field in the basin in the next few years was predicted. The results show that the MODFLOW model exhibits good performance in simulating the variation law of the groundwater flow field. The law of groundwater movement in the natural state of the basin was first determined by simulating groundwater movement in the basin under non-mining conditions, and the model as well as the accuracy of the completed system were verified. The drainage attributable to underground coal mining was input to the established MODFLOW model, and potential changes in groundwater depth throughout the basin were simulated. Through the correlation between groundwater depth and distribution of greensward land obtained by remote-sensing interpretation, the influence of mining on the future distribution of greensward land was predicted. The following conclusions can be drawn:

• The topographical features of HRB comprise arid and semi-arid dunes and denuded loess ridges. The distribution of greensward land is closely related to aquifer thickness and groundwater depth. Greensward land is mainly distributed in the upstream and midstream of the basin, especially in areas with shallow groundwater depth and thicker aquifers. The sporadic distribution of greensward land in some areas is related to artificial pumping irrigation. There are 34 central aggregation points of greensward land in HRB. Among them, 24% is distributed in areas with groundwater depth less than 3 m and only 8.82% in areas with groundwater depth more than 10 m. Furthermore, grassland accounts for only 11.11% of the total area in the downstream of the basin. The relationship between groundwater depth and central aggregation points of greensward land in the basin contributes to deeper understanding of the distribution of vulnerable vegetation in the Mu Us Sandy Land.
• The combined action of multiple underground mines has a strong impact on the phreatic aquifer system in the basin, and coal mine drainage changes the flow field and increases the depth of groundwater. The MODFLOW simulation of the groundwater aquifer flow field in 2020–2029 suggests that coal mining will have distinct effects on the groundwater aquifer above the first mining area. However, variations in the groundwater hydraulic gradient induced by coal mining would not be prominent in the next 10 years compared with the hydraulic gradient caused by differences in terrain. The groundwater model reveals three cones of groundwater depression centered on three mines, and the maximum drawdown may be as high as 5 m (Yph), 6 m (Dhz), and 10 m (Bls).
• According to the prediction of changes in groundwater depth induced by the combined mining of multiple underground mines in HRB, the groundwater level will prominently decrease around Yph, posing degradation risks to four central aggregation points of greensward land in the basin. It is necessary to implement timely measures to counter the changes. Although the coverage of greensward land is small in the first mining area of Dhz and Bls, the depth of groundwater decline is greater than that of Yph, and its impact on plant growth is more serious. This simulation approach provides a reference for the prediction of the distribution of vegetation with changes in groundwater depth.