Spatial-Temporal Changes in Ecosystem Service Value and Its Overlap with Coal Mining Intensity in the Yellow River Basin, China, During 2000–2030

Abstract

Understanding the ecosystem services and their interaction with coal resource development is crucial for formulating sustainable development policies. In this study, we focused on the Yellow River Basin, characterized by both rich coal resources and ecological fragility. The key findings are that (1) the ecosystem service value (ESV) in the Yellow River Basin exhibited significant spatial heterogeneity during 2000–2030, decreasing from the southeast to northwest, and decreasing the most notably in the southern part of the upper reaches of the river basin; (2) the high-high clustering area of the ESV shifted from the upper-middle reaches in 2000 to the middle-lower reaches in 2020, while the low-low clustering area remained within Inner Mongolia. By 2030, the high-high clustering area is expected to stabilize in southern Shaanxi Province, and the low-low area will potentially spread eastward; (3) the overall ESV is low, and it experienced a significant decline from 2000 to 2020, with water supply emerging as a major limiting factor, although some policy-supported counties had better ecological service values and trends. (4) From 2000 to 2020, the coal mining intensity (CMI) was concentrated in the upper and middle reaches, particularly at the junctions of Shanxi, Shaanxi, and Inner Mongolia, and the pattern remained stable, but local areas experienced increased mining intensity; (5) the overlap of the CMI and ESV primarily exhibited a low-high clustering pattern in the middle and upper reaches of the Yellow River Basin and eastern Ordos City, and a high-high clustering pattern in the middle reaches of the basin in Shanxi Province, which remained stable and slightly expanded from 2000 to 2030; (6) the trade-off between the ecosystem services in the overlap area intensified, especially between the provisioning and support services, and was significantly impacted by the coal mining activities. The findings indicate that the area that overlaps with the coal mining area in the Yellow River Basin has expanded and has had an increasing negative impact on the ESV. It is also essential to address the trade-offs between the provisioning and support services and to implement ecological restoration measures to mitigate the risk of ESV loss. Future efforts should focus on the regions where the CMI and ESV overlap and have poor coordination and the adverse effects of resource extraction on ecosystem services are becoming more pronounced. The results of this study demonstrate that spatial overlap analysis is effective in identifying the hotspots and provides a foundation for developing sustainable and high-quality policies for ecologically fragile basins.

1. Introduction

Ecosystem services could contribute to human survival, health, livelihood, and well-being. The human welfare generated by the ecosystem services is economically treated as the ecosystem service value (ESV) [1]. The significance of these services has generated interest in safeguarding them. Since the groundbreaking study conducted by Costanza et al. [2], there has been a significant increase in the global research focused on estimating, measuring, and mapping ecosystem services. This growth has been particularly evident since the millennium ecosystem assessment was conducted [3], an international study involving more than 1300 scientists. The millennium ecosystem assessment provided crucial evidence that about 60% of the world’s ecosystems have degraded over the past fifty years [4]. Given the range of the consequences of the human exploitation of ecosystems, it has become vitally important to prevent ecosystem degradation and to sustainably manage ecosystems and landscapes sustainably to ensure future ecosystem supply and meet the growing demands for ecosystem services [5]. This requires an understanding of the spatiotemporal interactions between human activities and ecosystems, enabling decision makers and policymakers to adopt sustainable management strategies and nature-based solutions, thereby enhancing the resilience and adaptability of ecosystems. Such an understanding can be derived from assessing the past land use/land cover dynamics [6], but it also necessitates forecasting the future situation, predicting the potential impacts of current decisions and trends, and revealing the trade-offs [7].

The Yellow River Basin is located in northern China and characterized by a wide range of natural features and socio-economic conditions [8]. It serves as a significant freshwater source for the northern region of China. The basin is commonly referred to as the “energy basin” because of its abundant reserves of coal, oil, natural gas, and non-ferrous metals. This makes it a crucial hub for China’s energy industry, chemical industry, raw materials, and basic industries [9]. Due to the challenging ecological conditions and economic growth pressures, the natural and semi-natural ecosystems in the Yellow River Basin have become more vulnerable [10,11]. The depletion of water resources and contraction of wetlands in the upper Yellow River pose a threat to the river’s resource supply to the upper Yangtze River [12]. The deep loess layer, arid climate, and intricate topography have resulted in significant soil erosion [13]. The ecological fragility of the region has worsened due to the increased resource destruction, environmental pollution, and geological disasters resulting from large-scale mineral resource development. To address these issues, a set of ecological restoration initiatives have been implemented in the Yellow River Basin, aiming to facilitate the rejuvenation of forests and grasslands. However, the effort to strengthen the ecosystem’s ability to withstand disruptions has been rather restricted. The conflict between the exploitation of resources and the preservation of the ecological environment continues to be a significant concern [14]. Hence, examining the correlation between the alterations of the ecosystem and the exploitation of mineral resources is of great importance for enhancing the ecological environment in the Yellow River Basin and achieving superior regional development.

Most previous studies have focused on the coupling and coordination relationship between resource exploitation and ecosystems [15,16,17], often forming a singular relationship at a global scale and neglecting the cluster patterns and spatial dependency of ecosystems and resource exploitation at the local scale. This has limited policymakers in formulating differentiated policies for resource exploitation and ecosystem protection, thereby reducing the efficiency of decision making. Moreover, most research on the impact of mineral resource development on ecosystem services has been based on static analyses of historical periods or single time points [15,18,19]. Assessing the trends of ecosystem services and predicting future ecosystem services are crucial for accurately grasping the future macro trends and formulating sustainable regional development policies.

Research conducted on continuous timescales is of great significant as it not only effectively identifies spatial and temporal changes but also demonstrates resilience in intricate landscapes, thereby closely correlating the results with the spatiotemporal dependencies [20]. An important problem is that time series image data have rarely been utilized to calculate the ecosystem service value despite the continuous land use/land cover changes over time [16,21,22]. This is primarily due to the time-consuming and labor-intensive nature of acquiring and processing time series remote sensing data, which deter most researchers interested in studying the ESV. The Google Earth Engine has been made available to the scientific community, providing an efficient resolution of this issue. The Google Earth Engine has compiled an extensive assortment of commonly utilized spatial datasets that can be analyzed through online tools. In addition, the Google Earth Engine has created a web-based interactive development environment for users. The Google Earth Engine application programming interface (API) and interactive development environment facilitate the swift development of applications and visualization of outcomes [23]. The use of the Google Earth Engine has enabled the analysis of time series images, allowing researchers to obtain a more profound comprehension of the ecosystem services. Hence, the Google Earth Engine provides researchers with a wide range of analytical techniques and resources, facilitating comprehensive investigations of spatial data. This research will foster the sustainable advancement of ecosystems and the preservation and administration of ecosystems.

In summary, the correlation between the alterations of the ecosystem services and the exploitation of mineral resources at the scale of river basins has not been examined. Considering the knowledge gaps and the strategic importance of the Yellow River Basin, in this study, we assessed the evolution of the ESV under the mining activities in the Yellow River Basin and its interactive mechanisms. The specific objectives of this study were (1) to develop a dynamic ecosystem service assessment module based on the Google Earth Engine platform and spatial analysis technique; (2) to analyze the spatiotemporal characteristics of the ecosystem services in the Yellow River Basin from 2000 to 2020; (3) to predict the ESV in the Yellow River Basin in 2030; (4) to analyze the evolution of the trade-offs and synergies between the ecosystem services; and (5) to reveal the overlap between the ESV and coal mining intensity (CMI) at the global and regional scales. The overall goal was to assist policymakers in understanding the spatiotemporal interactions between mineral resource development and the ecosystem services to promote the formulation of refined and sustainable management strategies, thereby providing scientific support for balancing economic development and ecological restoration in large-scale basins.

2. Materials and Methods

2.1. Study Area

The Yellow River Basin (31°31′–43°31′ N, 89°19′–119°39′ E) is the sixth-largest basin in the world, and it has a significant impact on both China’s ecology and economy. It spans nine provinces and regions, covering an area of almost 1.86 million km² at elevations ranging from −64 m to 6819 m [14] (Figure 1). The basin contains 30.3% of China’s population and contributes 26.5% of the national gross domestic product (GDP) [24]. The climate ranges from semi-humid in the southeast to semi-arid and arid in the northwest, with corresponding decreases in temperature and precipitation. The basin functions as a crucial ecological corridor that links the Tibetan Plateau, the Loess Plateau, and the North China Plain. It also serves as an important ecological barrier in Northern China [9,12]. Nevertheless, the ongoing exploitation of resources has resulted in a substantial decrease in water flow, and the basin’s capacity to support resources and the ecosystem is approaching its critical threshold.

As China’s energy basin, the Yellow River Basin has a staggering 1563 coal mines, accounting for almost one-third of the nation’s coal mines. China’s total coal production capacity is roughly 2.97 billion tons per year, and the coal production of the study area accounts for nearly two-thirds of the country’s overall capacity. In 2017, Shanxi, Shaanxi, Inner Mongolia, Ningxia, and Gansu provinces, situated in the middle and upper parts of the basin, collectively generated 2.46 billion tons of coal, accounting for almost 70% of the entire coal production in the country [9]. This region is widely regarded as the most economically valuable and has the greatest potential for coal production and distribution in China. Coal mining in the Yellow River Basin and the rest of the country has contributed to China’s socio-economic growth, but it has also caused severe disturbances to the ecological environment [9,14,16,18,19].

2.2. Research Framework

The primary sequence of activities in this study was as follows. (1) ESV calculation: Utilizing the Google Earth Engine platform, the ESV in the Yellow River Basin in China during 2000–2020 was calculated. Considering the regional heterogeneity, spatial heterogeneity indices such as the net primary productivity and standardized precipitation evapotranspiration index were used to adjust the equivalent weights annually. Furthermore, to illustrate the influence of socio-economic progress on the ecosystem, the ESV was dynamically adjusted using the GDP and Engel’s coefficient. (2) Trend analysis: the Sen+Mann-Kendall test was employed to examine the change trend of the ESV. (3) ESV prediction: the statistical spatiotemporal kriging model was used to predict the spatial distribution of the ESV in the Yellow River Basin in 2030. (4) Autocorrelation analysis: The local indicators of the spatial association (LISA) model were applied to analyze the autocorrelation of the ecosystem services. The bi-variate LISA (BILISA) model was employed to examine the correlation between the ESV and the extent of mining development. (5) Trade-off and synergy evolution analysis: correlation analysis was conducted to explore the evolution of the trade-offs and synergies between the ecosystem services in the overlap area. The research framework is illustrated in Figure 2.

2.3. Data

The land use classification data were obtained from the moderate resolution imaging spectroradiometer (MODIS) Land Cover Type product (MCD12Q1 v061). The net primary productivity (NPP) data were obtained from the MOD17A3HGF v6.1 datasets, which are accessible in Google Earth Engine’s MODIS collection. The standardized precipitation evapotranspiration index (SPEI) data were obtained from the SPEIbase v2.9 datasets, which are hosted on the Google Earth Engine platform. The data for further elements of the ESV, such as the GDP and annual GDP index, were obtained from the statistical yearbooks published by the National Bureau of Statistics. The production data for mineral resources were obtained from the Global Coal Mine Tracker database. Considering that the earliest available MODIS data began in January 2000, and the latest data are updated to early 2021 due to the lag of remote sensing products and statistical yearbooks, we chose a 20-year range (2000–2020) for the research. The data utilized in this study are described in

Table 1. Variables utilized in the computation of the ecosystem service value.

Data Application	Unit	Source
Delineation of Land Systems	Class Index (0–6)	MCD12Q1 v061, 2000–2020 https://doi.org/10.5067/MODIS/MCD12Q1.061 (accessed on 10 September 2024)
Gross Domestic Product (GDP)	Chinese Yuan (CNY)/km2	National Bureau of Statistics, 2000–2020 (https://www.stats.gov.cn/) (accessed on 23 September 2024)
Per Capita Net Income (PCNI)	CNY
Net Primary Productivity (NPP)	kg·C/m2	MOD17A3HGF v6.1, 2000–2020 (https://lpdaac.usgs.gov/resources/data-action/aster-ultimate-2018-winter-olympics-observer/) (accessed on 10 September 2024)
Standardized Precipitation Evapotranspiration Index (SPEI)	mm/month	The Global SPEI database (SPEIbase) v2.9, 2000–2020 (doi:10.20350/digitalCSIC/15470) (accessed on 10 September 2024)
Coal Mining Intensity (CMI)	Continuous Index (0–1)	Global Coal Mine Tracker, Global Energy Monitor, April 2024 release (https://globalenergymonitor.org/projects/global-coal-mine-tracker/tracker-map/) (accessed on 10 September 2024)

.

2.4. Assessment of the ESV

In this study, we used the millennium ecosystem assessment methodology to categorize the ecosystem services into four main types: provisioning services (raw materials, food, and water), regulating services (climate and gas regulation and waste treatment), support services (biodiversity protection and soil formation), and cultural services (recreation and culture). An ecosystem service accounting model based on the Google Earth Engine platform [25] was applied to calculate the ESV. This process involved determining the land use and ecological function parameters followed by static and dynamic calculations incorporating social development variables.

The land use types were classified based on the system developed by the National Committee of the State Ministry of Agriculture (1984), including woodland, grassland, cropland, wetlands, water bodies, unused land, and built-up land. The calculation was based on the equivalence weight factors proposed by Xie et al. [26]. Xie et al. built on the methodology of Costanza et al. [2] by conducting surveys among 200 ecological experts in China, and developing a value equivalent factor table for terrestrial ecosystems. This table, updated in 2015 using the literature data, expert knowledge, and biomass methods, defines one equivalent factor as 1/7 of the market value of grain produced per hectare annually. In this study, we used data from the National Statistical Yearbook, and we determined that the average grain yield was 4405 kg/ha with a price of 3.05 CNY/kg in 2000, resulting in an ESV per hectare of CNY 1919.32.

Given the heterogeneity of the ecosystems in China, directly adopting the national average equivalence coefficient could lead to inaccurate estimates. Therefore, we recalculated these coefficients according to the specific conditions in the study area. Additionally, spatial and temporal adjustments were made by incorporating dynamic factors such as the NPP and SPEI to account for the evolving nature of the ecosystem services over time and across regions. This was achieved using the following formula:

$$F_{nij}=\left\{\begin{array}{c}{P}_{ij}\times F_1\\ R_{ij}\times F_2\end{array}\right\}$$

(1)

where $F_{nij}$ is the equivalent value per unit area of an ecosystem in the nth ecosystem service category in year j in study area i; $P_{ij}$ is the spatiotemporal adjustment factor of the NPP for this type of ecosystem in year j in study area i; and $R_{ij}$ is the spatiotemporal adjustment factor of the SPEI for this type of ecosystem in year j in study area i. F₁ includes gas regulation, raw materials, soil formation and retention, climate regulation, waste treatment, biodiversity protection, and recreation and culture. F₂ includes water supply and food production.

To assure the comparability of all the ESVs at different points in time, it is essential to compute the corresponding economic value, taking into consideration the impacts of the price level, inflation, and other relevant aspects. The studied basin relies on the development of coal mining resources as a key industry, which has a significant influence on the GDP. The equivalent economic value is calculated using the following equation, with 2000 as the constant and the base year for determining the yearly average ESV:

$${ESV}_C={\displaystyle \frac{{ESV}_s}{E_{avg}}}\times E_{an}$$

(2)

$$E_{an}={\displaystyle \frac{E_m}{{\prod }_{i=n}^m{\varnothing }_i}}$$

(3)

where ESV_c is the equivalent ESV. ESVs is the static ESV. E_avg, E_an, and E_m are economic values associated with a single weight factor. E_avg is the average value, and E_an is derived based on the E_m value for the current year m within the research period. The constant year (start year) is denoted by n. ${\varnothing }_i$ is the yearly GDP index in year i, which is a measure of the annual percentage change in the GDP.

As civilization progresses and the quality of life improves, individuals will become more concerned about ecosystem services, resulting in a greater inclination to financially support these services. Hence, it is necessary to modify the fixed ESV by incorporating the socio-economic coefficient. The S-shaped development curve of Pearl can be employed to depict the attributes of the ecological value [27].

$$ESV={ESV}_c\times A_c$$

(4)

$$A_c={\displaystyle \frac{1}{1+{exp}^{-t}}}$$

(5)

$$t={\displaystyle \frac{1}{E_n}}-3$$

(6)

where A_c is the dynamic adjustment coefficient; t is the stage of social development; and E_n is Engel’s coefficient.

To implement the methodology described above, we adapted the Google Earth Engine code developed by Liang et al. [25] to compute the ecosystem services and the ESV in the study area. This was accomplished utilizing the API and the code editor on the Google Earth Engine platform.

To implement the methodology described above, we adapted the Google Earth Engine code developed by Liang et al. [25] to compute the ecosystem services and the ESV in the study area. This was accomplished utilizing the API and the code editor on the Google Earth Engine platform. Google Earth Engine is a cloud computing tool developed by Google, which is characterized by the ability to process remote sensing and geographic information data of large areas online. Compared with standalone processing software such as ENVI and ArcGIS (Environmental Systems Research Institute, USA), it can process data online without downloading it to a local computer. And these data and their processing code are shared in the form of links. The calculation model included the following steps:

(1)

equivalent weight coefficient adjustment

(https://code.earthengine.google.com/5d8a6b6287dddb24181111e7ae06965e)(accessed on 10 September 2024);

(2)

ESV calculation

(https://code.earthengine.google.com/91cbd359658c52bf747cd0298cc0af95) (accessed on 10 September 2024);

(3)

dynamic ESV adjustment

(https://code.earthengine.google.com/91cbd359658c52bf747cd0298cc0af95) (accessed on 10 September 2024).

2.5. Spatial-Temporal Change Analysis of the ESV

2.5.1. Trend of the ESV

The Sen+Mann-Kendall (MK) test can measure the trend of the ESV. The Sen slope estimator can be utilized to compute the trend value of the changes in the ESV. Typically, it has been employed in conjunction with the MK non-parametric test. This involves calculating the Sen trend slope and then using the MK test to assess the significance of the trend [19,28].

Theil–Sen median trend analysis, also known as Sen trend analysis, is a reliable non-parametric statistical technique used to calculate trends. Sen trend analysis has an advantage over linear regression trend analysis as it can mitigate the impact of missing data and the shape of the data distribution in the time series. Additionally, it can reduce the interference caused by outliers in the time series. The computation formula for the Sen trend degree is as follows [29].

$${\beta }_{\mathit{ESV}}=median\left({\displaystyle \frac{{ESV}_j-{ESV}_i}{j-i}}\right),\forall j>i$$

(7)

where ${\mathit{ESV}}_j$ and ${ESV}_i$ are the ESV time series data. ${\beta }_{\mathit{ESV}}>0$ indicates an increasing trend, while ${\beta }_{\mathrm{ESV}}<0$ indicates a decreasing trend.

The Mann-Kendall test is commonly employed in conjunction with Sen trend analysis [29]. This method is a non-parametric statistical test that is robust to missing values and outliers, and it does not require any specific distribution for the sample data. The statistical test approach can be expressed as follows:

$$Z=\left\{\begin{array}{c}{\displaystyle \frac{S-1}{\sqrt{Var\left(S\right)}}} \left(S>0\right)\\ 0 \left(S=0\right)\\ {\displaystyle \frac{S-1}{\sqrt{Var\left(S\right)}}} \left(S<0\right)\end{array}\right.$$

(8)

$$S={\sum }_{j=1}^{n-1}{\sum }_{i=j+1}^{n}sign\left({ESV}_j-{ESV}_i\right)$$

(9)

$$Var(S)={\displaystyle \frac{n(n-1)(2n+5)}{18}}$$

(10)

$$sign\left(\theta \right)=\left\{\begin{array}{c}1 (\theta >0)\\ 0 \left(\theta =0\right)\\ -1 (\theta <0)\end{array}\right.$$

(11)

where ${ESV}_j$ and ${\mathrm{E}\mathrm{S}\mathrm{V}}_{\mathrm{i}}$ are the $ESV$ time series data; the sign function is denoted by the sign symbol; S is the test statistic; and Z is the standardized test statistic, which is employed to assess the presence of a significant trend in the time series data. A higher Z-value indicates a more pronounced trend [30,31]. n is the number of data points. For significance level α, $\left|Z\right|>Z_{1-\alpha /2}$ indicates that there is a significant change trend. In this study, the significance level α was set to 0.05 to determine the significance of the trend shift in the ESV time series.

2.5.2. Prediction of the ESV

Kriging is an optimal spatiotemporal prediction method, and it considers the statistical relationships of the process in both space and time. This technique was developed by Matheron [32], who built upon the research conducted by Krige [33], a mining engineer hailing from South Africa. Kriging was subsequently extended to spatiotemporal applications of the ESV by utilizing the spatiotemporal statistical modeling and interpolation features in the R package gstat.

2.5.3. Spatial Hotspot Analysis of the ESV

Spatial autocorrelation analysis was used to express the spatial correlation and spatial patterns of the ESV across the studied basin. The spatial hotspot analysis was conducted using the local Moran’s I statistic, which has been commonly used for detecting local spatial clustering patterns in data. The local Moran’s I measures the degree of spatial autocorrelation between a given location and the surrounding areas. The formula for the local Moran’s I is as follows:

$$I_i=Z_i{\sum }_{j=1}^n{w}_{ij}{Z}_j$$

(12)

where I_i is the local bi-variate Moran’s I in spatial unit i. $Z_i$ is the standardized value of the ESV at location i. The standardization is achieved by subtracting the mean and dividing by the standard deviation. $Z_j$ is the standardized value of the ESV at neighboring location j. w_ij is the spatial weight matrix, which defines the spatial relationship between locations i and j. In this study, a distance-based spatial weight matrix was employed, in which adjacent areas were assigned a weight of 1 and non-adjacent areas were assigned a weight of 0.

2.5.4. Correlation Analysis

In this study, we utilized Pearson’s correlation analysis to reveal the trade-offs and synergies among various ecosystem services. Pearson’s correlation analysis is generally considered to yield more accurate and consistent results compared to other correlation analysis methods, such as Spearman correlation analysis and Kendall correlation analysis, when examining the relationships between continuous data variables [34]. Therefore, it has been commonly employed to reduce the interactions between ecosystem services [17,35].

2.6. Assessment of the CMI

The mining intensity refers to the degree to which mining activities disturb a specific geographic area. In this study, we calculated the mining intensity of coal mines by considering a single administrative county as the unit of computation. The calculation was based on three indicators: the production of coal mines, the proven coal reserves, and the number of mining areas within a single county. These indicators were weighted to determine the overall mining intensity. Initially, we employed the kernel density method [32] to analyze and describe the geographic arrangement of the coal mine production and proven coal reserves. Furthermore, we applied Equation (13) to normalize the three indicators and subsequently combined them into a comprehensive index, namely the CMI, using the equal-weight summation method, which is expressed by Equation (14).

$$U_{ij}^{\prime }={\displaystyle \frac{U_{ij}-U_{j,min}}{U_{j,max}-U_{j,min}}}$$

(13)

where $U_{ij}^{\prime }$ is the normalized value of mining intensity indicator j in study unit i, $U_{ij}$ is the initial value of mining intensity indicator j, and U_j,_max and U_j,_min are the maximum and minimum values of mining intensity indicator j, respectively.

$${CMI}_i={\sum }_{j=1}^3{U}_{ij}^{\prime }\times w_j$$

(14)

where ${CMI}_i$ is the composite level index of the mining intensity, w_j is the weight of mining intensity indicator j, and its value is one-third in this study.

2.7. Overlap Analysis of the ESV and CMI

Local bivariate spatial autocorrelation analysis, which is expressed by Equation (15), examines the spatial correlation of the attributes inside various geographic units using LISA cluster maps [36]. The exact equation is as follows:

$$I_i^{xy}=Z_i^x{\sum }_{j=1}^n{w}_{ij}{Z}_j^y$$

(15)

where I_i is the local bi-variate Moran’s I in spatial unit i, and $Z_i^x$ and $Z_j^y$ are the normalized values of the variance for the observed ESV and CMI values in spatial units i and j, respectively. Four possible cluster patterns can be generated. (i) High-high aggregation refers to the combination of a high ESV and high CMI values. (ii) Low-low aggregation refers to the combination of a low ESV and low CMI values. (iii) Low-high aggregation refers to the combination of a low ESV and high CMI values. (iv) High-low aggregation refers to the combination of a high ESV and low CMI values.

3. Results

3.1. Spatial Distribution of the ESV from 2000 to 2030

Figure 3 illustrates the spatiotemporal patterns of the ESV from 2000 to 2030. Overall, the ESV exhibits significant heterogeneity, decreasing gradually from the southeastern regions (Shandong Province, Henan Province, the downstream areas of Shanxi Province, and southern Shaanxi Province) to the northwestern regions where the ecological conditions are poorer. Additionally, the ESVs are relatively low in the downstream areas of Henan Province and Shandong Province, as well as the Weihe River Basin in the middle reaches. From the spatiotemporal perspective, the ESV generally exhibits a diffusion-like decline originating from the northwest between 2000 and 2030, and the most pronounced decrease occurred in the southern areas of the upper reaches of the Yellow River where the rate of change is mostly less than −75.71. In these regions, the fragile ecological environment, severe water scarcity, serious soil erosion, strong wind erosion, and expanding desertification and soil salinization are the primary contributing factors. However, the rate of decline has slowed in the southern upper reaches, Weihe River Basin, and Qinling Mountains. As a major beneficiary of the Returning Farmland to Forest project, the increase in forested land in Shaanxi Province has contributed to improvements in both the ecological environment and ecosystem. In contrast, the downstream regions with widespread urban land use have experienced the smallest reduction in the ESV.

3.2. The Hotspots of the ESV

As shown in Figure 4, the LISA cluster and significance analysis revealed the spatial clustering characteristics from 2000 to 2030. The clustering of the ESV in most areas of the Yellow River Basin was not spatially significant, with significant areas primarily characterized by high-high and low-low ESV clusters. In 2000, the southernmost part of the study area was a high-value ESV cluster, mainly including northern Sichuan Province, southern Shaanxi Province, southern Gansu Province, and northern Henan Province. The central-western region of the Inner Mongolia section of the studied basin was primarily a low-value cluster. By 2010, the high-high cluster in northern Sichuan Province had diminished, a high-value cluster began to appear in the central-southern part of Shaanxi Province, and there were no significant changes in the low-value clusters. In 2020, the high-value cluster in Sichuan Province disappeared, and a high-high cluster evolved in the border area between Shanxi and Henan provinces. The low-value cluster shifted toward the northeast. However, according to the predicted ESV’s LISA clustering results for 2030, the high-high clustering pattern will be distributed in the southern part of Shaanxi Province, indicating that only this region’s high ESV will exhibit temporal and spatial stability. The low-low cluster will be distributed in the northwestern part of the study area, suggesting a long-term trend of spreading eastward.

3.3. The Trend of the ESV

According to Figure 5, there are significant differences among the ecosystem services across the studied basin, and each service exhibits a certain degree of spatial autocorrelation. Biodiversity protection, recreation and culture, climate regulation, and soil formation and retention are more prominent in the southeastern and southwestern parts of the study area, and they exhibit relatively strong performances in the central region. However, in the western and northwestern areas, such as the Inner Mongolia Plateau and the Qinghai-Tibet Plateau, these services exhibit low values and an extremely significant decreasing trend. Similar patterns can be observed in the food production and raw material services. Notably, although the values are relatively low in the upper reaches of the Yellow River, specifically in the Alxa League and Bayannur City, they remain stable. The water supply is the most significant limiting factor across the studied basin, exhibiting negative values in most regions. The waste treatment services are primarily concentrated in the upper reaches of the Yellow River and the Qinling Mountains, exhibiting a significant increasing trend, while most of the other regions exhibit a notable decrease. The value of the gas regulation service is generally low and exhibits an extremely significant decreasing trend during the study period, and it only performs well in certain areas near the Taihang Mountains in Shanxi Province. Overall, the spatial heterogeneity of the average ESV and its change trend across the studied basin is not pronounced. The overall value is low and experiences a sharp decline between 2000 and 2020. However, certain counties, such as Qingjian County in Shanxi Province, a national ecotourism county, Wudi County in Shanxi Province, a model county for the Returning Farmland to Forest project, and the Shiguai District in Baotou City, Inner Mongolia, a model autonomous county, exhibit better average ESVs and trends. These exceptional areas clearly benefit from the targeted policy support and incentives.

3.4. CMI During 2000–2030

As shown in Figure 6, from 2000 to 2020, the spatiotemporal distribution of the CMI index was primarily high in the upper and middle reaches of the studied river basin, particularly in Shanxi, Shaanxi, Inner Mongolia, Ningxia, and Gansu provinces. The CMI was notably high in the junction areas between Shanxi, Shaanxi, and Inner Mongolia. During the two decades, the overall spatiotemporal pattern of the CMI remained relatively stable, although the intensity shifted in certain areas. For example, the CMI continued to increase in eastern Ordos City, Inner Mongolia, and Shenmu City, Shaanxi Province. Additionally, the scale and intensity of the coal mining expanded in the southern parts of Shaanxi and Gansu provinces, located in the middle and upper reaches of the Yellow River Basin. In contrast, the CMI gradually weakened in northern Gansu Province, located in the upper reaches.

Coal mining, a key industry in the study area, produced 2.463 billion tons of coal in 2017, accounting for over 69.7% of the national coal production, and this trend is expected to continue. Due to the inability of domestic oil and gas production to meet the demands of national economic development, the focus of coal resource development has increasingly shifted to the upper and middle reaches of the studied river basin. However, these regions are located in arid and semi-arid areas with fragile ecosystems. The inevitable ecological damage caused by coal mining further exacerbates the problems of soil erosion, loss of arable land, and vegetation degradation [37].

3.5. Overlap of the ESV and CMI

As shown in Figure 7, the bivariate LISA clustering results for the ESV and CMI reveal a significant low-high clustering pattern, and the clusters are particularly concentrated in the upper and middle reaches of the studied river basin, especially in the eastern region of Ordos City. This area has hilly gullies and the Mu Us Desert, which experience severe wind erosion. Resource exploitation has further damaged the surface vegetation, intensifying land desertification, dust storms, soil erosion, grassland degradation, and resource depletion, thus making the ecosystem highly fragile [38]. However, the clustering of the ESV difference and CMI in this region exhibits a significant high-high clustering pattern, indicating that policymakers have effectively considered the dual pressures of ecological vulnerability and resource exploitation in their planning and have achieved positive outcomes.

The high-high clustering of the ESV and CMI is mainly distributed in Shanxi Province, in the middle reaches of the study area, and is interspersed with the low-high clustering areas. The region’s ecosystem services are strongly influenced by the topography. The high-high clusters are concentrated in the Taihang and Lüliang Mountains, while the low-high clusters are in the Taiyuan and Xinzhou basins. This is primarily due to the difficulty of developing mountainous areas, which experience less human interference. The different topographic conditions directly impact the land use patterns, creating a locked-in effect on the regional ecosystem service capacity. The mountain areas have a stronger ecosystem service supply capacity than the plains and basins, in which the ecosystem service performance often falls below the resource extraction intensity, leading to nearly imbalanced coupling coordination. However, the high-high clustering pattern of the ESV difference and CMI in this region has significantly contracted, while the low-high clustering pattern has notably strengthened, especially between 2010 and 2020. This change suggests that with increasing resource extraction intensity, particularly in flat and ecologically fragile basin areas, the negative impacts of resource exploitation on the ESV are becoming more pronounced. This trend requires increased attention to prevent the further degradation of the ecosystem and to implement effective measures to alleviate the conflict between resource utilization and ecological protection.

From 2000 to 2030, the spatial patterns of the high-high and high-low clusters of the ESV and CMI are remaining relatively stable, with a slight expansion. The high-high clustering pattern, bordered by the Lüliang Mountains, is expanding from east to west, while the low-low clustering pattern, centered in northeastern Inner Mongolia, is expanding outward. In addition, the high-high clustering pattern of the ESV change and CMI during 2000–2020 is in the Shanxi–Shaanxi–Inner Mongolia adjacent region, where the ESV decreases significantly and the CMI is high. This indicates that the resource exploitation activities in these regions have intensified, potentially placing greater pressure on the already fragile ecosystems. During 2020–2030, the central area of Shannxi Province will exhibit a high-high clustering pattern of the ESV change and CMI. Therefore, future attention should focus on enhancing the ecosystem protection and restoration in the Shanxi–Shaanxi–Inner Mongolia adjacent region and central Shaanxi Province to prevent further ecological degradation.

4. Discussion

4.1. Trade-Offs of the ESV in the Overlap Areas

In 2000, among the 36 pairs of ecosystem services, significant correlations were found between 34 pairs (p ≤ 0.05), and 30 pairs occurred in the overlap areas (p ≤ 0.05), with 30 pairs observed in overlap areas (Figure 8). In the Yellow River Basin, 31 pairs exhibited positive correlations, and 3 pairs exhibited negative correlations. The overlap areas contained 28 pairs with positive correlations and 2 pairs with negative correlations. Of these pairs, seventeen pairs were highly correlated (Pearson coefficient r ≥ 0.7), five pairs were moderately correlated (0.5 ≤ r ≤ 0.7), and eight pairs were weakly correlated (r < 0.7). The synergy between the support services and cultural services was the strongest in the Yellow River Basin. This was mainly because the cultural services are primarily related to economic variables and do not require extensive land resources. For instance, agricultural and forest landscapes can also become leisure spaces that attract tourists. However, the synergy between the support services and provisioning services was weaker, and it even exhibited a trade-off with the food production services. A trade-off existed between the climate regulation and provisioning services. This was likely because the areas with high soil retention had extensive forest vegetation, and there was land use competition between the forested and cultivated lands. Therefore, in this study, the trade-off may have been due to the spatial incompatibility caused by the dependency of the provisioning and regulating services on the specific land use types.

A comparison of the correlation coefficients of the ecosystem services (Figure 8, Figure 9 and Figure 10) revealed that the synergy between the raw materials and services such as water supply, waste treatment, climate regulation, and recreation and culture services decreased significantly in the overlap areas compared to the control counties. Specifically, the difference coefficients for the raw materials and waste treatment were 0.361 (2000) and 0.168 (2010), indicating that coal mining-related ecological damage hindered the promotion of other services by the raw materials [18]. Additionally, the trade-offs between food production and regulating services (e.g., climate regulation) strengthened, suggesting an increase in the agricultural burden due to coal mining [39]. The synergy between the soil formation and retention services and gas regulation services declined by 0.510 (2000) and 0.121 (2010), indicating that the degradation of the soil quality affected the gas regulation capacity [40]. Interestingly, the overlap areas had a stronger synergy between the water supply services and the supporting and regulating services compared to the control counties. This could be because the underground water was influenced by the mining activities. By 2020, the overlap areas shifted from a balanced relationship between food production and regulating services to a synergistic relationship. The waste treatment services also improved, possibly due to the implementation of sustainable management practices [41].

A comparison of the changes in the trade-off and synergy coefficients of the overlap areas throughout the study period revealed that the synergy coefficients between the water supply services and the recreation, waste treatment, and climate regulation services remained consistent. However, the correlation coefficient between the water supply and food production services decreased from −0.184 to −0.226 (by 0.043). The synergy coefficient between the food production and raw materials decreased by 0.215 and that between the food production and soil formation and retention services decreased by 0.269. Furthermore, the trade-offs between the food production services and the waste treatment, climate regulation, and gas regulation services strengthened. In particular, the gas regulation services exhibited a stronger synergy with the other services, especially the soil formation and retention services. This may be because China’s new energy practices have contributed to the reduction in greenhouse gas emissions and thus have improved the soil quality [42]. The synergy between the raw material sand food production continued to decrease, particularly between 2010 and 2020, with a decrease of 0.502. The shift from a synergistic relationship to a trade-off between the raw materials and water supply, as indicated by the change in the correlation coefficient from 0.005 to −0.050 between 2010 and 2020, likely resulted from the insufficient water supply and environmental pressures caused by the mining activities [19]. In conclusion, the supportive role of the water supply in the food production and raw materials in the study area decreased over time. The coal mining activities have affected the equilibrium between the raw materials and supporting, regulating, and cultural services and have impeded the promotion of the regulating services by the soil formation and retention services. The changes in the environmental services within the study area reflect the conflict between economic development and environmental protection [18,43].

4.2. Implications, Limitations, and Future Work

Starting in the twentieth century, the Chinese government has undertaken a range of ecological and environmental conservation initiatives in the study area [44]. As shown in Figure 3, these have somewhat improved the regional ecological environment. However, as shown in Figure 6, the synergies between the ecosystem services in coal mining-affected counties, specifically between raw materials and provisioning services and between gas regulation and supporting services, have continuously decreased. Moreover, the gap with the control countries has widened. Simple re-vegetation appears not to have formed a well-functioning ecosystem with multiple services and functions working in synergy. This may be due to coal mining disrupting the surface and causing subsidence and fissures, which directly affect the soil structure and hydrological conditions, causing soil nutrient loss and subsequently impacting the vegetation growth and soil nutrient cycling [41,45]. Additionally, the re-planted vegetation struggles to replicate the structure of the native communities, and its recovery may require a long period of time [46]. This suggests that the ecological restoration should focus on the synergistic effects of soil formation and conservation, diversity protection, and gas regulation services to enhance the ecosystem stability, thereby improving the recovery efficiency and capacity of the ecosystem.

Ecological restoration should be tailored to the local conditions. Figure 7 shows that the issue of uneven regional development is particularly pronounced in the coal mining provinces. Taking Shanxi Province as an example, according to the announcement by the Shanxi Provincial Energy Bureau regarding the production capacity of coal mines in Shanxi Province, there were a total of 668 coal mines in Shanxi at the end of December 2020 [47]. Long-term intensive coal mining has left the ecological environment fragile and the ecological carrying capacity low. The analysis of the spatiotemporal evolution of the ESVs in counties in Shanxi Province revealed that the ecosystem service capacity of the study area is mainly influenced by the topography, and the areas with better ecosystem service capacities are concentrated in the Taihang and Lüliang mountain regions. However, the ecosystem service capacity in the Lüliang Mountain area is significantly lower than that in the Taihang Mountain area. This is mainly because the Lüliang Mountains are located in a typical loess hilly and gully region, where severe soil erosion has occurred. The water conservation service per unit area is only 320.9 m³/ha, 20% lower than the provincial average, and the mountainous region is adjacent to the area where mineral resources have been found. Human activities are regularly conducted in this area. In 2021, the Lüliang Mountain area had a total of 91 coal mines, accounting for 13.62% of the province’s total number of coal mines. The extensive operations and insufficient management measures have amplified the strain on the ecosystem services and diminished the quality of the regional ecosystem services in this area [48]. The northwestern part of the study area in Inner Mongolia has been identified as a location with limited value in terms of ecosystem services. This is mostly due to the unfavorable natural climate conditions and frequent occurrence of sandstorm catastrophes in this area. This area has a mid-temperate continental semi-arid monsoon climate and experiences frequent strong winds and little rainfall, resulting in a low upper limit of its ecosystem service capacity [49,50,51]. Therefore, it is strongly recommended that policymakers and decision makers consider these factors during the implementation of ecological restoration, develop more targeted restoration policies, and consult experts during the restoration process to enhance the effectiveness of the restoration measures. New technology of coal mining such as water-preserved mining, and backfill mining, should be widely applied to protect the ecosystem services such as water supply.

In this study, we used the administration cells, namely the counties within the Yellow River Basin, to assess the changes in the ESV and its relationship with the CMI. It is beneficial to assess the official statistical data and make management decisions for each administration cell. Nevertheless, the spatial division of the towns and cities is relatively imprecise [52], which hinders the comprehension of the cause-and-effect relationships between the ecosystem services [53]. Through a survey, we found that the mean size of the towns and cities was 27,979 km², ranging from 1656 to 253,015 km². To mitigate the impact of the changes in the geographic unit area, the indicators of ecosystem services and socio-ecological factors pertaining to the towns and cities were measured by normalizing them with respect to the area. Nevertheless, the spatial unit areas analyzed in this study were significantly larger compared to those used in previous studies that employed the same ecosystem service valuation approach [54,55,56]. This technique was unsuccessful in capturing the detailed ecological and human-related factors within these huge spatial units, which are especially crucial for understanding and predicting the outcomes regarding the ecosystem services [56]. Additionally, the dynamic ESV adjustment variables used in this study, such as the GDP and Engel’s coefficient, have a national-scale precision, which may lead to the overestimation of the ESV in some of the western provinces due to the uneven economic development in these areas [57]. Therefore, future research should assess and analyze the ESV across the entire study area at multiple spatial scales such as the geographic grid, individual coal mine, and watershed scales.

5. Conclusions

In this study, we utilized abundant data from the Google Earth Engine platform, comprehensively assessed the ecosystem services in the Yellow River Basin from 2000 to 2020, analyzed the evolution patterns and trends of the ecosystem services, predicted the future ESV, and determined the interactive stress mechanisms between the ESV and coal mining. The findings of this study are summarized below:

• From 2000 to 2020, the overall spatial distribution of the ESV increased from east to west; however, from 2000 and 2030, the area with a decreasing ESV trend is expanding from the initial position in the northwest. The most significant decreases occurred in the southern upstream region of the Yellow River Basin where most of the rates of change were less than −75.71.
• From 2000 to 2030, the high-high clustering areas of the ESV in the Yellow River Basin are shifting from Sichuan and southern Shaanxi in 2000 to the central-southern part of Shaanxi and the border area between Shanxi and Henan, while the low-low clustering areas are remaining stable in the central-western part of Inner Mongolia, with a trend of spreading eastward. By 2030, the high-high areas are expected to be stable and to be in the southern part of Shaanxi, while the low-low areas may continue to spread eastward.
• The ecosystem services in the study area have exhibited significant spatial heterogeneity. The gas regulation, biodiversity protection, and recreation and culture services are stronger in the southeastern and southwestern regions and weaker in the northwestern regions such as Inner Mongolia and the Tibetan Plateau. The provision of the ecosystem services has decreased over the past 20 years, but certain services along the Yellow River, such as water supply and waste treatment services, have exhibited increasing trends.
• From 2000 to 2030, the CMI in the study area is mainly concentrated in the upstream and midstream areas, especially at the junction of Shanxi, Shaanxi, and Inner Mongolia, where the overall pattern is stable but the intensity of the mining continues to increase in localized areas such as Ordos and Shenmu cities. From 2000 to 2020, the ESVs of these areas decreased significantly and the CMI intensified.
• The development of the coal resources in the Yellow River Basin must take into account the value of the ecosystem services. It is necessary to avoid setting up coal mines in highly sensitive areas during the planning phase. During the mining phase, methods such as water-preserved mining should be adopted to reduce the impact on the hydrological system, while ecological restoration measures should be taken to restore the ecosystem functions, thereby reducing the trade-off among ecosystem services and enabling more ecosystem services to be preserved during mining.