The Pingshuo open-pit coal mining area is situated in the northern Shanxi province of China, and its geographical coordinates are 112°10′ 113°30′ E?39°23′ 39°37′ N (see Figure 1). It is located in the north of the Ningwu Coal Field in Shuozhou city and is linked to the world's largest coal field in the Black Triangle of the Loess Plateau—Dongsheng. It covers a large area [21 kilometers (km) in length from north to south and 22 km in width from east to west] and has a mine life of hundreds of years. The mine has a geological reserve of 12.75 gross tonnages (GT) and a construction scale of 65 metric tons (Mt), including 45 Mt of large national open pit mines and 20 Mt of local coal mines. This national large-scale open-pit coal mining area is divided into three fields: Antaibao open-pit coal mining area, Anjialing open-pit coal mining area, and the east open-pit coal mining area, with an annual average mine output of 20–30 Mt (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.

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%.

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).

The results of the biophysical quantities for the GEP in Pingshuo are shown in Supplementary Table 2. Over the past 30 years, the monetary quantities of GEP increased from RMB¥ 3.34 billion Chinese yuan (RMB¥) to 46.67 billion RMB¥, accompanied by an increase in the production value of coal mining in Pingshuo from 0.65 to 30.07 billion RMB¥ (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 5 shows that the boundaries of the open-pit mining areas were generally within the values of L–L, and the values of L–L gradually increased. The values of H–H were mainly situated in the northeast of the Antaibao and Anjialing open-pit coal mining areas after 1995. There were also H–L and L–H outliers, which were sporadic. The values of H–L were close to the area of L–L, while the values of L–H were close to the area of H–H.

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, low-density and cold spot areas were the main areas in the study area.

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.

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 1995–2000 and 2000–2005, the changes in GDP on GEP were reduced to simple segmented linear functions. However, between 1990–1995 and 2005–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.

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² = 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.

Based on the GEP, this study assessed and monitored the value of ES in Pingshuo, a large open-pit coal mining area in China, from 1990 to 2020. Before 2010, GEP presented a 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). The reasons for these changes in the GEP are relevant to the mining history and construction policies of Pingshuo. Pingshuo is composed of three open-pit mines as a large-scale open-pit coal mining area. The primary open-pit mine of Pingshuo was put into operation in 1987. Before 2009, although the mining policy was less concerned with the ecological environment, the agricultural land, forest, and grassland were the main LULC types in Pingshuo. The wide distribution of vegetated areas ensured a steady increase in GEP (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 2010 for the last three decades. Moreover, the simultaneous growth in the values and growth rates of GEP in 2015 and 2020 resulted from the positive implementation of the “green mine construction” in Pingshuo. In 2014, the China Federation of Mining Industry issued a notice on “National Green Mine Pilot 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” (Yu et al., 2022).

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). 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.

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. 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.