A new pattern to quantitatively evaluate the value of ecosystem services in the large-scale open-pit coal mining area

Introduction Open-pit coal mining could disrupt the ecosystem and lead to the loss of service values for the ecosystem through direct occupation or indirect impacts on adjacent ecosystems. Methods In this research, we combined a new accounting system, gross ecosystem product (GEP), with spatial–temporal analyses to quantify the ecological variation and explore its driving factors in Pingshuo, a large-scale open-pit coal mining area in China. GEP is an aggregate accounting system that can summarize the value of provisioning, regulating, and cultural ecosystem services (ES) in a single monetary metric. The spatial–temporal approaches used in our study were known as exploratory spatial data analyses and interpretable models in machine learning. Both spatial and non-spatial data, including remote sensing images, meteorological data, and official statistics, were applied in the research. Results The results indicated the following: (i) From 1990 to 2020, the annual average growth rates of GEP decreased from 30.78 to 9.1%. Furthermore, the classified results of GEP revealed that the regions with rich ES quality rapidly reduced from 51.90 to 32.18%. (ii) Spatial correlation of GEP was significant, and the degree of spatial clustering was relatively high in the mining areas. Moreover, the mining areas also continually presented concentrated high-density and hot spot areas of GEP changes. (iii) The spatial–temporal effects were notable in the relationship between GEP and three socioeconomic factors, i.e., the mining effects, human activity intensity, and gross domestic product (GDP). (iv) The win–win development for both the economy and ecological environment in Pingshuo could be realized by restricting the annual growth rate of mining areas to between 4.56 and 5.03%. Discussion The accounting results and spatial–temporal analyses of GEP will contribute to the future regional sustainable development and ecosystem management in Pingshuo.


. Introduction
Coal is the most important source of energy and sustains social and economic progress for human wellbeing (Qian et al., 2018;Zhang et al., 2021). While mineral resources can provide the impetus for the rapid development of the economy, mining and processing of coal could result in ecological and environmental problems (Worlanyo and Jiangfeng, 2021;Xiang et al., 2021). It is significant to note that, among the various types of mining, open-pit coal mining is one of the most destructive activities in terms of its impact on the surface landscape and the surrounding environment (Kibria, 2014;Balasubramanian, 2016;Lechner et al., 2017). The reason for this is that the construction, exploitation, and rehabilitation of an open-pit mining site essentially involve large-scale earth movements in space and time. Furthermore, these movements often lead to massive and irreversible changes to landscape, such as removing vegetation cover, stripping topsoil, and reshaping topography . Because of the spatial and temporal features of open-pit mining, the negative impacts on the ecosystem are also spatial and temporal. Specifically, the adverse effects occur not only in the mining area but also in the surroundings (Sonter et al., 2017). Moreover, mining projects generate unique and diverse environmental impacts at every stage of the process (Lechner et al., 2017;Wang Z. Y. et al., 2018). The precise measurement of ecosystem services (ES) in open-pit coal mining areas, with spatial and temporal dimensions, is the first step in seeking out methods to minimize the negative impacts on the ecosystem, thus achieving sustainable development of mining areas, and it will help provide great convenience to further realizing the necessity of green mines (Hu et al., 2020).
Ecosystem services (ES) are the benefits people directly or indirectly derive from functioning ecosystems (MA, 2005). Since the 1990's, researchers have been studying ecosystem functions and attempting to quantify the ES they provide for human survival and development (Daily, 1997;Costanza et al., 1998). Ecosystem service value (ESV) and gross ecosystem product (GEP) were then proposed, in turn, to provide the fundamental theories and methodologies for ES accounting (Daily et al., 2009;Ouyang et al., 2013). Due to the relative ease of ESV accounting, most researchers have assessed the ES of coal mining areas through ESV Qian et al., 2018;Xiao et al., 2020). Initially, studies of ES in mining areas were carried out following the environmental impact assessment requirements (Rigina, 2002;Gomez and Herbert, 2015;Preece et al., 2016). Later, researchers increasingly considered the spatial and temporal characteristics of the coal mining areas in their assessment of the ESV (Qian et al., 2018;Liu and Chen, 2021). Wang et al. (2020) assessed ES by mapping and quantifying changes at multiple spatial and temporal scales, confirming that ES in mining areas has distinct spatial patterns in different years. Yang et al. (2022) further focused on the spatial and temporal evolution of past mining and exposed the cumulative ecological losses due to open-pit mining. In addition, the research on the driving factors of ES changes has gradually increased. Furthermore, the relationship between changes in land use/land cover (LULC) and ESV was the most extensively researched topic (Bian and Lu, 2013;Kumar et al., 2022;Zhang L. J. et al., 2022). However, studies still need to be conducted on quantitative analysis of ES changes, especially on the driver of economic development and physical conditions within mining areas. More research on how these factors affect ES and how they change over time needs to be conducted.
While the unit value-based approach is more convenient for ES accounting in ESV (Gashaw et al., 2018;, the method cannot differentiate the value of the same type of ES in different regions (Zou et al., 2020). The accounting methods for ES measurement models in the GEP can address this problem by considering the variation in ecosystem types and quality conditions in different locations. It was in 1985 that Hannon first came up with this theory (Hannon, 1985), and the basic concept of the GEP has been further defined by Chinese scholars (Ouyang et al., 2013;Ma et al., 2015). GEP is a statistical and accounting system that translates the ES to the economy into monetary terms (Ouyang et al., 2013). Currently, the National Development and Reform Commission of China and other related departments have launched pilot projects for GEP in 15 provinces and 23 cities (https://seea. un.org/home/Natural-Capital-Accounting-Project). Moreover, the International Union for Conservation of Nature (IUCN) and its partners have conducted national pilot studies of GEP accounting projects in Ordos city, Tonghua city, Xing'an county, and Xishui county in China (IUCN, 2017). At the same time, an increasing number of researchers have investigated GEP for global and national terrestrial ecosystems (Jiang H. Q. et al., 2021;, provinces (Ouyang et al., 2013, and megacities (Zou et al., 2020). Moreover, there has been a general recognition that the way we measure development needs to change to include nature's contribution to human wellbeing, i.e., ES, in decision-making (Ma et al., 2015;Nations, 2015;Hao et al., 2022). However, there were a few examples of ES accounting and measuring the development of mining areas through GEP. The spatial and temporal explicit characterization of GEP concerning such heterogeneity has yet to be well-studied and therefore necessarily limits our ability to manage ecosystems in mining areas sustainably.
In this research, we focused on Pingshuo and evaluated the value of ES in Pingshuo by GEP. As the largest open-pit coal mining area in China , Pingshuo is a typical comprehensive land use area for mining, with the influence of coal mining and high-speed urbanization (Zhang et al., 2019a). The specific goals are as follows: (i) calculating the GEP with the help of spatial and non-spatial data to verify its feasibility in the large-scale open-pit mining area; (ii) applying the exploratory spatial data analysis methods to analyze the spatial-temporal dynamics of GEP; (iii) identifying the relationship between the four types of driving forces, including socioeconomic factors, degree of human interference, neighborhood effects of LULC, and physical conditions, and the variations of GEP applied with GeoDetector and the interpretable models in machine learning; (iv) proposing a practical measure that is beneficial to the sustainable development of Pingshuo.
. Materials and methods

. . Study area
The Pingshuo open-pit coal mining area is situated in the northern Shanxi province of China, and its geographical coordinates are 112 Figure 1).  (Bai et al., 2000). Pingshuo has a typical temperate semi-arid continental monsoon climate, with few drops of rain, cold, dry, and windy winters, and cool and less windy summers. It is located in the eastern part of the Loess Plateau (see Figure 1), a highly eroded and ecologically fragile system. The ecological environment of Pingshuo is weak. It is sensitive to environmental changes and has less plasticity to maintain stability. However, the natural secondary forest has been destroyed, and the extensive grassland communities, which are rare, are scattered with low vegetation cover due to a long development history and a high reclamation index (Bai and Yun, 2018). Therefore, it is vital for the sustainable development of Pingshuo to maintain regional land and environmental security.

. . Data sources
Spatial and publicly available official non-spatial data are the primary data used in this research. In particular, spatial data were obtained from data product service websites, namely, the National Tibetan Plateau Scientific Data Center, the Resource and Environment Science and Data Center of the Chinese Academy of Sciences, and NASA EarthData. The official non-spatial data, including statistics and monitoring data, have been maintained by the product service websites in China. Moreover, Table 1 summarizes the sources, types, resolutions, and functions of data in the GEP accounting process.
The data for LULC were obtained from Landsat 1990-2010 Thematic Mapper (TM), Landsat 2015, and 2020 Operational Land Imager (OLI) with the strip numbers 126/33 and 126/32 in the Geospatial Data Cloud, and Landsat images were chosen because of their extended temporal coverage. For all these images, atmospheric correction and radiation normalization are essential to remove the effects of smoke, dust, and haze, as well as changes in the sun's angle (Chavez, 1996). Radiometric and atmospheric corrections were performed for these 14 scenes over seven periods using The Environment for Visualizing Images (ENVI) 5.3. The geometric correction was then performed on the above seven periods of images using the previously corrected Landsat TM images, and the geometric accuracy after the correction was 0.5 pixels (Luo et al., 2016). Moreover, the images with different periods also need to be aligned by polynomial and nearest neighbor methods, and the root mean square error after alignment can be guaranteed within 0.5 pixels (Zhang et al., 2017). Finally, the simultaneous and adjacent images were stitched and cropped using the vector boundary of the study area.
A combination of supervised classification and visual interpretation was utilized to classify the LULC types into 10 categories: shrubs, forests, grasslands, farmlands, water areas, industrial and mining lands, urban lands, rural lands, bare lands, and roads. The classification process is carried out by referencing the remote sensing monitoring dataset of China's LULC released by the Resource and Environment Science Data Center of the Chinese Academy of Sciences (Xu, 2023) and applying Google Earth to verify the LULC types. The confusion matrix method was used to assess the accuracy of the interpretation of these remote sensing images, and the overall accuracy was >85%.

. . Methodology
The whole framework of our studies consists of three parts: the first part is the data preparation and pre-processing; the second part is the GEP accounting with temporal and spatial features; and the last part is the spatial-temporal exploration of the GEP in Pingshuo (see Figure 2).

. . . Evaluation methods for GEP
The ecosystem types in Pingshuo are mainly comprised of farmland, forest, and grassland ecosystems. Moreover, after working in the field, we observed that no ecological area in Pingshuo currently provides cultural services. Based on these actual situations, the availability of data, the connotation of GEP, the values of provisioning service, and eight secondary indicators of regulating service accounted for GEP in Pingshuo.
. . . . Ecosystem product value (EPV) The ES that humans directly obtain from the ecosystem to meet the material needs of human life and to finish the production of derived products that can be traded on the market are called provisioning services. The value of the provisioning services in GEP is known as the ecosystem product value (EPV).
The annual output of the agriculture, forestry, livestock, fisheries, and other ecological energy sources such as fuelwood, straw, and biogas was taken as the biophysical quantities for EPV. Moreover, the annual production values were the monetary quantities. The formulas for the biophysical and monetary quantities of EPV are shown in Supplementary Table 1. The output values of agricultural and livestock products contain the value of labor and other human inputs . The EPV has been modified to include only the share of natural inputs; the share of labor and other human inputs have been deducted. Based on the relationship between soil fertility contribution and crop yield in China, we set a transfer rate of 32.91% for agricultural products (Tang and Huang, 2009). Based on the factor contribution constituted of livestock, we set a transfer rate of 36.5% for livestock products (Liu and Dong, 2014).

. . . . Ecosystem regulating value (ERV)
Ecosystems can provide not only direct provisioning services to humans but also functions to improve the human living and survival environment (Sannigrahi et al., 2018;Zou et al., 2020), such as pollutant degradation and carbon dioxide fixation. Moreover, the value of the regulating services in GEP is known as the . /fevo. .  ecosystem regulating value (ERV). Considering the real situation in Pingshuo, eight indicators were identified for ERV: water retention, soil retention, sandstorm retention, carbon sequestration, oxygen supply, air purification, water purification, and climate regulation.
The biophysical quantities of each indicator were valued using different models and methods, and the monetary quantities, where market values were unavailable, were usually valued using the replacement cost method. The formulas for the biophysical and . /fevo. .  It is worth noting that we applied the geographic information system (GIS) concept in the ERV accounting process. Specifically, the inverse distance weighted interpolation method was used to create a spatial distribution model of the monitoring data. The relationship between the indicators of GEP and LULC was mapped to the corresponding characteristics using raster calculations, thus completing the spatialization of the statistics. The resampling method mainly addressed the problem of inconsistent scales in spatial data. Consequently, all data resolutions were consistent with LULC at 30 m.
. . . Exploration methods for spatial-temporal heterogeneity with GEP Spatial heterogeneity can be revealed by the Local Moran's I (LISA), one of the exploratory spatial data analysis (ESDA) methods. It is a decomposed form of the Global Moran's I (Anselin, 2010), and its formula is given in Equation 1. The range of LISA is from −1 to 1, and the closer the value of LISA is to 1, the stronger the positive spatial correlation is; the closer the value of LISA is to −1, the stronger the negative spatial correlation is; the closer the value of LISA is to 0, the more it is irrelevant to the spatial dimension. Accordingly, the LISA clustering maps can be divided into five categories, i.e., the positive correlation is represented by the high-high (H-H) type and the low-low (L-L) type in clustering maps; the negative correlation is represented by the high-low (H-L) type and the low-high (L-H) type; and the insignificant type represents the spatial decorrelation. In addition, the spatial characteristics of the variations for GEP within every 5 years were studied and observed through the density and hot spot analysis in ArcGIS 10.2 with the Getis-Ord Gi * index, whose equations are shown in Equations 2, 3.
where n is the number of grids; y i and y j are the attribute values of the spatial object y at the i and j points; y is the average value of y; W ij is the spatial weight matrix element, appraising the spatial connection between the point i and the other points, which can be . /fevo. .
where E(G * i ) and Var(G * i ) are represented as the mathematical expectation and variance, respectively, and W ij is the spatial weight.
. . . Exploring the driving forces of spatial-temporal heterogeneity with GEP . . . . Driving factors The potential driving factors influencing the GEP changes and their spatial heterogeneity in this study were selected based on the experimental conditions of Pingshuo, literature reviews Cheng et al., 2019;Han et al., 2020), and the principle of data availability. Moreover, the driving factors fall into four classes: socioeconomic factors, degree of human interference, neighborhood effects of LULC, and physical conditions. Specifically, (i) socioeconomic factors included population density (POP), human activity intensity, gross domestic product (GDP), and the extent of mining area effects; and in this research, human activity intensity was represented by night lights; (ii) the formula for the integrated indicator of human disturbances (Hemeroby) is shown in Equation 4 (Han et al., 2020); (iv) the formula for the neighborhood effects of the six LULC types  is given in Equation 5; and (v) physical conditions including annual average temperature, precipitation, wind speed, and NDVI.
where D is the Hemeroby; m is the number of land types; HI i is the coefficient of land type i (Table 2); S i is the area of land type i; and S is the total study area.
where F i,j,d is the enrichment of neighborhood d of location j with land use type i; n i,d,j is the number of cells with the i th land use type in the neighborhood d of grid j; n d,j is the number of cells in the neighborhood d; N i is the number of cells with the i th land use type in the whole raster, and N represents all cells of the study . /fevo. . area. We set d = 2 with a 5 * 5 rectangular window to figure up the neighborhood interactions. The neighborhood interactions of croplands, grasslands, forests, water areas, built-up lands, and bare lands were obtained by the "Focal Statistics" tool in ArcGIS 10.2.
. . . . Driving factor analyses The initial driving factors used for driving force analyses may be highly correlated, which can lead to a faulty system analysis (Dormann et al., 2013). Multicollinearity analyses and Pearson's correlation coefficient were applied to eliminate this problem. Specifically, multicollinearity, which uses variance inflation factors (VIFs) and tolerance VIF, can be used to determine the standard errors of these factors (Allison, 2001;Liao and Valliant, 2012). The higher the standard error, the greater the multicollinearity. If the VIF was more remarkable than 10 or the tolerance was < 0.1, it was confirmed that multicollinearity existed in the data set (Hair et al., 2006). The correlation coefficient of the two factors was evaluated by Pearson's correlation coefficient (Booth et al., 1994), and the formula is defined as follows. When the value between the two factors is >0.7, it means that there is a correlation between them.
where cov(X, Y) is expressed as the covariance of the conditioning factors X and Y; σ X and σ Y are the standard deviations of X and Y, respectively.
. . . . Driving forces exploration methods GeoDetector is a statistical approach for detecting spatial heterogeneity and uncovering its potential driving forces. By analyzing the distinctions between the regional and total variances, the factor detection in GeoDetector can be used to determine the extent to which the potential driving factor explains the spatial divergence (Wang et al., 2010;Wang and Hu, 2012). The value of the q statistic can reflect how the dependent variable was affected by the values of one or more independent variables, which are between 0 and 1. When q is 0, the driving factor is independent of the dependent variable. When q is 1, the independent variable controls the spatial distribution of the dependent variable (Wang and Xu, 2017). The value of q is calculated by the formula as follows: where SSW and SST are the sums of variance within the classified region and the whole region, respectively; h is the classification of the variable or independent variable; N h and N are the number of rasters in the classified area h and the whole area of the study area, and σ 2 are the variances of the classified region h and the whole region of the study area. The permutation feature importance and partial dependency plot (PDP) methods in machine learning interpretable models were used to investigate the drivers of GEP changes. In the permutation feature importance method, a factor is "important" if shuffling its values increases the model error because, in this case, the model relied on the factor to make the prediction (Wieder et al., 2014;Fisher et al., 2019). Moreover, PDP can reveal the relationship between the GEP and the factors by marginalizing the model output over the feature distribution (Greenwell et al., 2018).
. /fevo. .   (Table 3). Moreover, the monetary quantities of EPV increased from 13.46 to 128.67 million RMB , and the monetary quantities of ERV increased from 3.33 to 46.54 billion RMB . The proportion of EPV in Pingshuo remained below 0.5% (see Figure 3). This indicates that the contributions of ES to Pingshuo were mainly derived from the regulating service.
As the units of the biophysical quantities of GEP were inconsistent (Supplementary Table 2), the following research focused on the monetary quantities of GEP, abbreviated as GEP. According to its value range, the visualization of the GEP was carried out using the natural break method in ArcGIS 10.2, and the classified results are shown in Figure 4 in descending order: Level 1 (F), Level 2 (E), Level 3 (D), Level 4 (C), Level 5 (B), and Level 6 (A). In different years, the spatial distribution of GEP showed different characteristics (see Figure 4). Moreover, the proportion of the lowest grade (Level 1) increased from 4.67 to 15.33%, while the higher grades (Level 4 and above) decreased from 51.90 to 32.18% (Table 4). Therefore, it was essential to study the spatial and temporal heterogeneity of GEP.
Four significant correlations were found in Pingshuo, of which H-H and L-L were the main types (see Figure 5).  Figure 6 is the LISA significance map, which can reflect the significance level of the GEP, and it could be observed that, within the L-L and H-H areas, the significance level of the GEP was 0.001, with some values at the significance level of 0.01. In other words, the GEP of these regions was strongly positively related to the GEP of neighboring regions. In the mining areas, the significance levels were substantially >0.05, which meant that the spatial distribution of GEP was firmly convergence and showed a strong correlation. Moreover, the areas of significance (p = 0.01) and the higher significance (p = 0.001) manifested an increasing trend before 2015.
Before analyzing the density and hot spot, the spatialized GEP was resampled through the 300 * 300 m fishnet and applied using the spatial linkage method in ArcGIS 10.2. The detailed spatial distributions of the GEP variation density and hot spot are presented in Figures 7, 8, supplemented by changes in the boundaries of adjacent non-contemporaneous mining areas. Areas of high density and hot spots were concentrated in and around mining areas, and these areas gradually increased before 2015. Compared with the shrunken mining areas, the extended areas were likely to lead to the variation of GEP (see Figures 7, 8). While the low-density and cold spot regions were mainly distributed within other areas, the areas of low violent artificial disturbance and, excluding 2010-2015, lowdensity and cold spot areas were the main areas in the study area.

. . Results of driving factor analyses
The results revealed that temperature (tolerance = 0.018, VIF = 56.941) and POP (tolerance = 0.025, VIF = 40.412) had minor tolerances and the most significant VIF values (Table 5). Neither the tolerances of the two driving factors nor the VIF values meet the threshold values (VIF <10 or tolerance> 0.1). The maximum correlation value was between POP and night lights (0.875), and it was followed by the correlation value between GDP and Hemeroby (0.753; see Figure 9), and both of the two values were >0.7, indicating that the POP and night lights, GDP, and Hemeroby are strongly collinear. Therefore, we removed the driving factors of temperature, POP, and Hemeroby before the processes of driving force analyses.
. /fevo. .    . . Driving forces of the spatial-temporal heterogeneity for GEP Gross domestic product (GDP; q = 0.348), NDVI (q = 0.297), night lights (q = 0.295), precipitation (q = 0.249), mining effects (q = 0.164), and wind speed (q = 0.147) had stable and dominant solid effects on the spatial heterogeneity of GEP (Table 6). However, the values of q for the neighborhood effects of LULC were all below 10%, indicating that they played a small role in the spatial distribution heterogeneity of GEP, and the impacts of neighborhood effects of LULC on GEP tended to remain constant compared with other factors (see Figure 10).
Gross domestic product (GDP; mean weight = 0.400), mining effects (mean weight = 0.325), and night lights (mean weight = 0.270) were the most prominent factors contributing to the variation in GEP (see Figure 11, Table 7). However, the various factors influencing how GEP varied (see Figure 12). For example, GDP and night lights demonstrated negative trends in the variations of GEP. At the same time, the mining effects were divided into two different trends by the boundary line of the mine. Indeed, the effects of different factors on GEP varied greatly, and even the effects of the same factors could change over time. For example, the GDP and GEP were negatively correlated, with the degree of impact varying yearly. In 1995In -2005, the changes in GDP on GEP were reduced to simple segmented linear functions. However, between 1990However, between -1995However, between and 2005However, between -2010, the effects became more complex and required the intervention of non-linear functions to explain them.
As for the bivariate effects, adding new variables usually resulted in a change in the trend of the single factor, although this effect could be greater or lesser (see Figure 12). For instance, night lights were more likely to cause a change in the trend of GDP than the mining effects. Generally, the bivariate results would be more varied and complex, and specific analyses would be needed for each year and situation.

. . Exploratory research on economic development and ecological balance in Pingshuo
While the growth of GEP from 1990 to 2020 increased from 3.34 to 46.67 billion RMB , the annual average growth rates of GEP decreased from 30.78 to 9.1% (see Figure 13), and the areas of GEP with high levels were decreasing year by year (see Table 4). Furthermore, the exploitation of coal resources played an essential role in the spatial distribution and the spatial-temporal heterogeneity of GEP. Therefore, we sought to investigate how to maximize the mining economy while ensuring that the impact on the ecosystem in the open-pit coal mining area was minimized by controlling the exploration seed of mining areas.
The regression analysis was carried out for the mining areas and the production values of coal mining from 1990 to 2020, and the results showed that the mining areas had a 98.2% explanatory rate for the production values (R 2 = 0.982), which is the reason why we seek the balance between the economy and the ecosystem in Pingshuo with the help of mining areas. Until 2010, the increasing annual growth rate of the mining area was associated with the decreasing growth rate of the GEP. After 2010, when the annual growth rate of the mining area fell to 3.89%, the GEP experienced a simultaneous increase in both growth rates and values (Table 8). The growth rates of the mining areas in Table 8 were calculated as the ratio of the growth of the mining area to the initial area.
However, when the annual growth rates of the mining areas were <3.89%, the economic development of Pingshuo did not achieve the most optimal solution. To obtain the optimal solution, we accounted for the GEP and the production values from 1990 to 2020 after setting the annual growth rates of the mining areas within a range of 3.89-7.19%. It was found that, when the annual growth rate of the mining areas was restricted to between 4.56 and 5.03%, the variations in the annual growth rates of the GEP could be maintained between 0.02 and 0.08%. At that time, the GEP could reach 48.46-51.55 billion RMB , which is 1.04-1.10 times that of the current one, while the output values of mining could reach 23.91-26.79 billion RMB , which is 0.80-0.89 times that of the current ones. /fevo. .

FIGURE
Pearson correlation coe cient between the pairs of driving factors. NDVI represents the normalized di erence vegetation index. NE represents neighborhood e ects. Hemeroby represents an integrated indicator of human disturbances. GDP represents gross domestic product. POP represents population density.
gradually increasing tendency, rising from 3.34 to 31.91 billion RMB (Table 3). However, the GEP in 2010 was 13.97% lower than that in 2005 (Table 3). Following that, GEP and the growth rate of GEP showed a steady growth trend from 2015 to 2020 (Table 8) (Zhang et al., 2019a). After 2009, the overall production of the three open-pit mines entered a period of rapid development (Xu et al., 2021). The ecosystem had the lowest stability due to the intensification of mining activities and inadequate monitoring during 2009-2013. During this period, mining and dumping areas became the most crucial land types in Pingshuo, along with cultivated land and forest (Zhang et al., 2019a), which is the reason why GEP was lower in 2010 than in 2005 and 2015. After further research, we also found that mining effects, GDP, and night lights all had a decisive and negative impact on the GEP (see Figure 13), which reached its maximum in . /fevo. .  Unit Acceptance Measures (for Trial Implementation)" (China Federation of Mining Industry, 2014). As the third batch of pilots, Pingshuo has paid particular attention to mining area governance and green reclamation, actively implementing "green mine construction" .

FIGURE
Results of the driver factors' importance assessment in Pingshuo for GEP changing. NDVI represents the normalized di erence vegetation index. NE represents neighborhood e ects. Hemeroby represents an integrated indicator of human disturbances. GDP represents gross domestic product. POP represents population density. Regarding components, ERV was an essential component of GEP in Pingshuo (see Figure 3). However, the percentages of ERV were more significant compared with other types of regions. Pingshuo experienced poor production conditions, fewer crop varieties, and lower yields due to natural constraints, such as the environment and soil. It also limited the output and production of livestock, all of which resulted in a relatively low share of EPV. At the same time, the share of the different ES in ERV can be divided into two types. The years 1990 and 1995 have similar profiles and differ from the latter two decades (Supplementary Table 3).

FIGURE
Partial dependency diagrams of the first three important driving factors. NDVI represents the normalized di erence vegetation index. NE represents neighborhood e ects. Hemeroby represents an integrated indicator of human disturbances. GDP represents gross domestic product. POP represents population density.
Specifically, in 1990 and 1995, water retention and oxygen supply services were considered essential ERV sources. In contrast, carbon sequestration replaced water retention as a vital part of ERV after 2000. Finally, it should be noted that the ERV was not included in our research due to the uniqueness of the regional functions and the different ecological aptitudes. Some of the results of this study were consistent with previous research (Ma et al., 2017;Miao et al., 2019), while others were unique to this mining area. Both of these indicate that it was practicable to monitor the value of ES in the large open-pit coal mining area by producing an estimate of GEP using existing data and methods.

. . E ects of mining on GEP in Pingshuo
We used the spatial-temporal analysis to quantitatively identify the driving factors for the variation and spatial distribution of GEP. The results indicated that mining effects had essential influences on the spatial-temporal variation (Table 7) and spatial distribution (Table 6) of GEP in Pingshuo.
From the point of view of variation, although the monetary quantities of the GEP were basically in the ascending states, the growth rates had a general tendency to fall (Supplementary Tables 3, 4). In addition, the number of regions with a high quality of ES decreased rapidly from 1990 to 2020. Specifically, the areas of higher levels (C-A) of GEP decreased from 51.90 to 32.18%, and the lowest level (F) increased from 4.67 to 15.33% (Table 4). Moreover, the expanded mining areas in Pingshuo were the main areas where higher values were converted to lower values. In comparison, the shrunken areas were the main areas where lower values were converted to higher values (see Supplementary Figure 1). We analyzed that the exclusive functional areas of open-pit coal mining, such as stripping areas and waste dumps, could damage forests and reduce vegetation cover and biomass (Zhang et al., 2019b). This resulted in the decline of ES with water retention, soil retention, sandstorm retention, carbon sequestration, oxygen supply, and climate regulation in Pingshuo (Supplementary Table 3). Moreover, the land degradation caused by the exploitation of coal resources cannot be ignored (Bian et al., 2020), which leads to the deduction of the ES with sandstorm retention.
The spatial correlation was significant in terms of spatial distribution, and the degree of spatial clustering was relatively high in Pingshuo. The H-H values were located in the northwest of Pingshuo, where no mining exploration had taken place.
. /fevo. .  Meanwhile, L-L was situated in the coal mining areas (see Figure 5), with significant levels >0.05 (see Figure 6). As these areas have often been disturbed by mining exploration and other relevant human activities, the quality of the surrounding ecological environment has deteriorated. Moreover, the mining areas are always accompanied by drastic spatial variations in GEP (see Figures 7, 8). This is due to the inherent characteristics of duration, spatial extent, and frequent ecological disturbance that affect ES in open-pit mining areas (Gvozdkova et al., 2019;Xu et al., 2019). It is clear that, when mining development is chosen in ecologically fragile areas, economic gains are associated with corresponding ecological pressure and loss of ES value. As the exploitation of mineral resources has contributed significantly to the spatial heterogeneity and variability of the GEP (Tables 7, 8), we investigated the balance between economic and ecosystem development in Pingshuo by controlling the rate of exploration in mining areas. Moreover, we found that the value of ES could be improved when the annual growth rate of the mining areas in Pingshuo was restricted to between 4.56 and 5.03%. Previous experiments have also confirmed that the rational usage of the land, the manageable exploitation of mineral resources, and the reducing land occupation were all beneficial for the local ecosystem (Xu et al., 2021). This may have practical implications for future resource development and utilization in Pingshuo, which could ensure the balance between mineral supply and ecosystem sustainability.

. Conclusion
In this study, GEP was implemented in a large-scale open-pit coal mining area to quantitatively assess and monitor the value of ES from 1990 to 2020. The study presents the following conclusion: 1) The GEP changes in Pingshuo have been significant over the past 30 years, and the mining areas also continuously showed concentrated high-density and hot spot areas of GEP changes. While the GEP increased from 3.34 to 46.67 billion RMB , the annual average growth rate of GEP decreased from 30.78 to 9.1%, and the number of regions with high ES quality fell rapidly from 51.90 to 32.18%. 2) The spatial correlation of GEP was pronounced, the degree of spatial clustering was relatively high, and the significance levels of GEP in the mining areas were significantly higher than 0.05. The areas of significance (p = 0.01) and the . /fevo. .
higher significance (p = 0.001) showed an increasing trend before 2015. 3) Three socioeconomic factors such as GDP (q = 0.348, mean weight = 0.400), night lights (q = 0.295, mean weight = 0.270), and mining effects (q = 0.164, mean weight = 0.325) had significant effects on the spatial heterogeneity and spatial-temporal variation of GEP, and the impact on GEP varies with factors and time. 4) The GEP of the mining areas can be guaranteed to grow in value and growth rate by limiting the annual growth rate of the mining areas to between 4.56 and 5.03%.
Since the accounting methods for GEP belong to the interdisciplinary subject, the processing is relatively diverse in indicators, model algorithms, and parameter selection. Moreover, GEP accounting is still at an early stage of development. Apart from the services we mentioned and accounted for in this research, many services still need more clarity in the formation mechanisms of ES (Pandeya et al., 2016). A typical example is preventing and reducing the possibility of mass movements by ecosystems, such as debris flows, landslides, mudflows, rock falls, etc. (Kirschbaum and Stanley, 2018). Although extremely important for human survival and economic development, the quantitative relationships between ecosystems and mass movement regulation have yet to be thoroughly investigated (Bruzzese et al., 2021). Moreover, even for the types of services generally recognized by humans, it is difficult to quantify the values using immature technical methods and unreliable data. In addition, ERV indicators not available in market values are typically measured by replacement cost in this study. Even though this direct market-based method has a relatively low level of uncertainty in assessing the value of the ERV (TEEB, 2010), there are still inherent uncertainties in the valuation methodologies. For example, some common pricing mechanisms, such as carbon and industrial oxygen and the cost of reservoir construction, have different price ranges and can vary widely from different sources.
This study provides a novel way of investigating ES in largescale open-pit coal mining areas, which can clarify the value differences of the same type of ecosystem in different regions. Despite these limitations, GEP can reflect the contribution of an open-pit coal mining area's ecosystems to human wellbeing, and the spatial and temporal trend of GEP can reflect changes in the quality of ecosystems in the area and the reasons for changes. Our research findings may be instructive not only for the future development of Pingshuo but also for policies aimed at promoting sustainability in other mining areas.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.