Spatiotemporal Distribution Characteristics and Their Driving Forces of Ecological Service Value in Transitional Geospace: A Case Study in the Upper Reaches of the Minjiang River, China

Abstract

Ecosystem service value (ESV) is a key indicator for evaluating ecosystem services. Thus, a unique quantitative assessment instrument that comprehensively and objectively evaluates ESV is of great significance for protecting regional ecosystems and achieving sustainable development. Based on data for meteorology, hydrology, soil use, and land use, this paper comprehensively employs the InVEST model, spatial autocorrelation, and geographic detectors to study the spatiotemporal characteristics and driving forces of spatial variations in ESV in the upper reaches of the Minjiang River. The results indicate the following: (1) The ecosystem service capacity of the study area has continuously improved, with the ecosystem service value (ESV) increasing by USD 4.078 billion over 20 years. Soil conservation has made the most significant contribution to the growth of ESV, accounting for over 85%. (2) The distribution of ESV exhibits a “lower in the northwest, higher in the southeast” trend. The Moran’s I value for each year exceeds 0.7, indicating characteristics of High–High and Low–Low aggregation. (3) Slope plays a dominant role in causing the spatial differentiation of ESV, contributing 30.9%. Slope is followed by HAI at 19.7% and the urbanization rate at 16.8%. Rainfall has the least impact at 4%. (4) The results from the multi-factorial interactions reveal that all factors experience synergistic enhancement effects when interacting. The spatiotemporal differentiation of ESV is the result of multiple factors acting in conjunction, underscoring the importance of coordinated efforts in ecological restoration and comprehensive environmental management in the upper reaches of the Minjiang River. The methodology of this research could be applied to assess the impact of natural changes and human activities on ESV. The findings offer theoretical support for regional resource and environmental management, as well as ecological compensation decision making.

1. Introduction

Ecosystem assessments over the years have indicated that ecosystem services play an irreplaceable role in the process of sustainable development for humankind, the ecological environment, and the economy [1]. As complex ecosystems, basins contribute to a variety of ecosystem services. Basins also feature a high level of human activity [2]. While human beings benefit from ecosystem services, they often ignore the destruction of watershed ecosystems, resulting in water scarcity [3], air pollution [4], soil erosion [5], and vegetation degradation [6]. The resulting climate change will directly impact the capacity of ecosystem services [7]. The Global Environment Outlook 6 (GEO-6) highlights that the Earth has suffered extreme damage. Without protection of the environment, the Earth’s ecosystems and the sustainable development endeavors of humanity will face increasingly severe threats. To avoid excessive exploitation of the ecosystem and negative results of ecological service functions, it is critical to study the spatiotemporal characteristics and driving forces of ESV in such a basin.

Quantitative assessment of ecosystem service value (ESV) is of the utmost importance. Previously, the most widely used method for this task was the equivalent factor method [8]. The equivalent factor method used is based on the average ESV per unit area and the corresponding area [9]. However, current equivalent factors are based on global or national averages, necessitating some adjustments based on the actual situation [10]. In recent years, using traditional empirical equations for ESV calculations can no longer fully meet research requirements. Relying on ecological benefit calculation models to depict and simulate the functionalities of ecosystem services has become a notable assessment approach [11]. For instance, models such as VEST, ARIES, Eco AIM, EPM, and NAIS are being increasingly utilized in ESV assessments [11,12].

In studying the driving factors in the evolution of ESV, scholars have explored the impact of multiple factors such as urbanization, population, economic development, industrial structure, and climate change on changes in ecosystem service value [13]. However, most scholars have focused on exploring how land use change and climate change affect ecosystem services [13,14,15], which is insufficient to comprehensively describe the spatiotemporal evolution pattern of ecosystem service values. With the advancement of remote sensing (RS), geographic information systems (GIS), and other technologies, scholars have begun to conduct interdisciplinary and more diversified research. Increasingly, more models and computational methods, such as correlation analysis, multiple regression analysis, principal component analysis, and structural equation models [16,17], have been introduced into the study of driving factors of ESV changes. However, these methods do not reflect the interplay between driving forces, nor do they reflect the heterogeneities of different driving forces in different regions. In recent years, methods such as Geographically Weighted Regression and geographic detectors have been introduced into the study of driving forces of spatial variations in ESV. These methods can interpret the interplay between driving forces and dominant regions of ESV distribution and avoid the strong interplay between forces, which can affect the accuracy of the research results. As a result, these methods have gradually become important topics of discussion [18,19,20].

Transitional geospace represents a complex area of interplay between nature and mankind, with obvious ecological, industrial, and social transitions, presenting the uniqueness and complexity of the territorial system of human–environment interactions. Transitional geospace is, therefore, a key region for the high-quality development and governance of territorial space [21]. Scholars have conducted research at different scales, including global [22], national [23,24], regional [16,25], and watershed [17,26,27] scales. Authors have also conducted research in different ecosystems, such as forests, farmland, and wetlands [28,29,30]. Alternatively, the research can focus on specific areas such as natural reserves [31,32], urban groups [33], and urban economic zones [34,35]. However, there is still a lack of research on the ESV of transitional geospace, as well as its evolution and driving forces.

The upper reaches of the Minjiang River are located in the transitional area from the Qinghai–Tibet Plateau to the Sichuan Basin. This area is a typical region of ecological vulnerability in western China, with high mountains, steep valleys, frequent geological hazards, and obvious climate change. In this study, we focus on research gaps related to ESV in transitional geospaces, such as empirical equations struggling to meet the needs of decision-makers and the difficulty of quantitatively detecting the driving forces as well as their interplays of ESV spatial variations through correlation analysis and structural equation models. The objectives of this study are as follows. (1) Using the upper reaches of the Min River as a case study, this paper attempts to assess the region’s ESV based on multi-source data and the InVEST model. (2) Global Moran’s I and Getis-Ord Gi* are employed to analyze the spatiotemporal evolution characteristics of ESV. (3) A geographic detector is used to investigate the key driving factors and their intensities behind the spatiotemporal differentiation of ESV. This study breaks through the limitations of administrative scales and more intuitively reflects the changes in the ESV of each grid unit. This model identifies regional ecosystem issues and their causes, not only offering insights for the assessment of ESV in transitional geographical spaces but also providing a scientific basis for maintaining the regional ecological balance and promoting sustainable development.

2. Materials and Methods

2.1. Study Area

The upper reaches of the Minjiang River (north latitude 30°45′~33°10′, east longitude 102°35′~103°57′) are located in southwest China (Figure 1), with a total length of 330 km and covers a total area of 2.48 × 10⁴ km. It is a transitional area between the Qinghai–Tibet Plateau and the Sichuan Basin, as well as an important water source and ecological barrier for the Chengdu Plain [36]. This region contains high mountains and steep valleys, with an elevation ranging from 780 m to 6200 m. This area features rough terrain, fold and fault structure development, and active new tectonic movements, making it one of the most geologically disaster-prone areas in China [37]. Consequently, this region is ecologically fragile and vulnerable. This region has a plateau mountain monsoon climate, with an annual rainfall of 500 mm~850 mm [38]. Here, soil and vegetation have an obvious vertical distribution rule. The main soil types include subalpine meadow soil, frigid desert soil, and yellow earth. The main vegetation types are broad-leaved forests, mixed coniferous and broad-leaved forests, coniferous forests, shrubs, and meadows. These different types create biodiverse ecosystems [39,40].

2.2. The Data Source

Reliable data, including DEM and administrative boundary data, were obtained from the Data Center for Resources and Environmental Science, CAS. Data for temperature, rainfall, and decreases in evapotranspiration were sourced from the National Earth System Science Data Center, National Science and Technology Infrastructure of China. Soil data were sourced from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn). Data on the available soil water in plants, coefficient of velocity, slope, length, soil saturated hydraulic conductivity, and soil root depth were obtained through basic data processing and calculations. Land use data were sourced from the Resource and Environment Science Data Center of the Chinese Academy of Sciences, with a spatial resolution of 30 m. With reference to the Chinese land use classification standard, land use types were reclassified into six categories: arable land, forest land, grassland, water area, construction land, and unutilized land. The land use data for the study area were obtained through processes such as cropping and projection. Economic data such as population density, level of urbanization, regional GDP, per capita GDP, and the proportion of primary industry were generated from various levels of statistical yearbooks and bulletins. For the main data sources, see

Table 1. Main data sources and types.

The Data Source	Type	Resolution	Year	Data Source
DEM	Raster	30 m	-	http://www.geodata.cn, (accessed on 12 December 2022)
Administrative boundaries	Vector	-	2020	http://www.geodata.cn, (accessed on 25 October 2022)
Land use	Raster	30 m	2000, 2005, 2010, 2015, 2020	http://www.resdc.cn, (accessed on 25 October 2022)
Potential evapotranspiration decreases	Raster	1 km	2000, 2005, 2010, 2015, 2020	http://www.geodata.cn/, (accessed on 5 January 2023)
soil dataset	Raster	1 km	2014	https://data.tpdc.ac.cn/, (accessed on 5 January 2023)
Temperature	Raster	1 km	2000, 2005, 2010, 2015, 2020	http://www.geodata.cn/, (accessed on 10 March 2023)
Average annual rainfall	Raster	1 km	2000, 2005, 2010, 2015, 2020	http://www.geodata.cn/, (accessed on 10 March 2023)
Annual average temperature	Raster	1 km	2000, 2005, 2010, 2015, 2020	http://www.geodata.cn/, (accessed on 10 March 2023)
NDVI	Raster	30 m	2000, 2005, 2010, 2015, 2020	http://www.resdc.cn/, (accessed on 11 March 2023)
Per kilometer of GDP spatial distribution (GDP)	Raster	1 km	2000, 2005, 2010, 2015, 2020	https://www.resdc.cn/, (accessed on 11 March 2023)
Per kilometer of population spatial distribution (PD)	Raster	1 km	2000, 2005, 2010, 2015, 2020	http://www.ornl.gov/sci/landscan/, (accessed on 11 March 2023)
Other socio-economic data	Text	-	2000, 2005, 2010, 2015, 2020	Sichuan Statistical Yearbook, A BA Prefecture Yearbook, and National Economic and Statistical Communiqué of the People’s Republic of China on the National Economic and Social Development

.

2.3. Methods

The paper employs the InVEST model combined with valuation assessment methods. Using ArcGIS, we compute the ecological service value of the research area. Following this process, the spatial clustering characteristics of ecological system service values at different times are examined using global autocorrelation analysis and hotspot analysis. Finally, a geographic detector is applied to analyze the main driving factors and forces behind the spatial differentiation of ecological system service value. The technical order of operations of this research is shown in Figure 2.

2.3.1. Calculation of ESV Based on InVEST

The InVEST model quantifies the impacts of land/water changes on ecosystem services through the production function approach combined with GIS tools to estimate the material or economic value of ecosystem services [41]. The InVEST model facilitates the visual representation of assessment results and is sensitive to data variations, which is why this model has been widely adopted for the evaluation of ecosystem service functions. In this paper, the InVEST model is applied to assess the carbon storage, water conservation, and soil conservation values in the upper reaches of the Min River. The calculation principles are explained as follows:

(1)

Calculation of Carbon Storage Value

The carbon storage model of the InVEST model divides the carbon storage of the ecosystem into four basic carbon pools: above-ground carbon, underground carbon, soil carbon, and dead organic matter carbon [19].

The calculation formula for total carbon storage is as follows:

$$C_{total}=C_{above}+C_{below}+C_{soil}+C_{dead}$$

(1)

where C_total is the total carbon storage, C_above is the above-ground carbon storage, C_below is the underground carbon storage, C_soil is the soil carbon storage, and C_dead is the dead organic matter carbon storage.

In the results of the calculation, the carbon storage value is calculated by applying the market value method:

$$V_C=P_C\times C_{total}$$

(2)

where V_C represents the value of carbon sequestration in the ecosystem in USD, C_total refers to the total carbon storage (in tons), and P_C represents the value of carbon stocks per unit in USD·t⁻¹.

Here, the value of carbon storage is primarily assessed based on transaction data from the EU carbon trading market, carbon tax law, and afforestation cost method [42,43,44]. Given the continuous improvement of China’s carbon trading market over the past few decades and the release of the “First Compliance Period Report on China’s Carbon Emission Trading Market” in December 2022, the carbon storage value in this paper is based on the average transaction price data published in this report, which is USD 5.96 per ton. Through the conversion coefficient of 0.2727 between C and CO₂, the unit carbon storage value is calculated as follows [12]:

$$P_C=\frac{S}{0.2727}$$

(3)

where S is the average trading price in China’s carbon market in USD·t⁻¹.

(2)

Calculation of Water Conservation Value

The water source conservation component of the InVEST model is based on the principle of hydrologic cycle, and the water yield of the watershed is obtained through parameters such as rainfall, plant transpiration, and surface evaporation. The water yield data are then corrected using the topographic index, soil saturated hydraulic conductivity, and flow velocity coefficient to determine the water source conservation [45]:

$$Y_{jx}=\left(1-\frac{AE{T}_{xj}}{P_x}\right)\times P_x=\left(1-\frac{1+{\omega }_x{R}_{xj}}{1+{\omega }_x{R}_{xj}+\frac{1}{R_{xj}}}\right)\times P_x$$

(4)

where $Y_{jx}$ represents the annual water yield (m³), $P_x$ refers to the annual average rainfall (mm) of grid cell x, and $AE{T}_{xj}$ serves as the annual average of evapotranspiration decreases in grid cell x on land use type j. In addition, $R_{xj}$ is the aridity index of grid cell x on land use type j, which represents the ratio of potential decreases in evaporation to rainfall, and ${\omega }_x$ is the corrected ratio of vegetation’s annual available water to rainfall.

Based on the calculated amount of water conservation, a shadow project is selected to further calculate the water conservation value. The calculation formula is as follows [46]:

$$V_W=Y_{jx}\times A$$

(5)

where V_W is the value of water conservation in USD, $Y_{jx}$ is the annual water yield (m³), and A is the unit volume of the reservoir cost. This study adopts a unit storage value cost of 0.88 USD·m⁻³ with reference to the Zipingpu Hydropower Station in the research area.

(3)

Calculation of Soil Conservation Value

Next, we calculate the amount of soil erosion (RKLS) under bare ground conditions using the SDR module in the InVEST model, which is the potential amount of soil erosion, as well as the actual amount of soil erosion (USLE) under the influence of vegetation cover and soil conservation. Then, soil conservation (SD) is taken as the difference between the two aforementioned conditions:

$$USLE=R\times K\times LS\times C\times P$$

(6)

$$RKLE=R\times K\times LS$$

(7)

$$SD=RKLE-USLE$$

(8)

where R is the factor of rainfall erosivity (MJ·mm (ha·hr)⁻¹), K refers to the factor of soil erodibility (t·ha·hr (MJ·ha·mm)⁻¹), LS serves as the factor of slope and length, C represents the factor of vegetation cover and crop management, and P stands for the factor of soil and water conservation measures.

Soil conservation seeks to achieve three main goals: producing solid soil, fertilizer maintenance, and avoiding siltation. The calculation formulas are as follows [47]:

$$V_S=V_i+V_j+V_k$$

(9)

where V_S means the value of soil conservation, V_i stands for the value of solid soil, V_j embodies the value of fertilizer maintenance, and V_k refers to the value of avoiding siltation. The calculation process for V_i, V_j, and V_k was detailed by Xu Wenxiu et al. [47].

The soil fixation value is assessed using the opportunity cost method. The abandoned land area caused by soil erosion can be determined based on the amount of soil conservation and the average thickness of the soil surface layer. The economic value of the loss is then estimated using the annual average benefit of local natural vegetation [47]. The calculation formulas are as follows:

$$V_i=\frac{Q}{\rho h}M$$

(10)

where Q represents the soil conservation quantity, measured in tons (t), ρ stands for soil bulk density, in units of t·m⁻³, h denotes the average soil thickness, and M signifies the annual average benefit of natural vegetation in USD·ha⁻¹. Within the study area, the average soil bulk density is 1.195 t·m⁻³. According to the 2020 Sichuan Province Statistical Yearbook, the annual average benefit of natural vegetation is 22,021.91 USD·ha⁻¹.

Soil erosion results in the loss of nutrients, including nitrogen (N), phosphorus (P), and potassium (K). The conservation value of soil fertility can be estimated using the price of chemical fertilizers; the soil conservation quantity; the content of N, P, and K in the soil; and the nutrient content in a given fertilizer [12]. The calculation formula is as follows:

$$V_j=Q\times \left(a_N\times b_N\times M_N+a_P\times b_P\times M_P+a_K\times b_K\times M_K\right)$$

(11)

where a represents the nutrient content in the soil, b is the nutrient content in the fertilizer, and p is the price of the fertilizer. Based on field research, the fertilizers chosen for this study are diammonium phosphate and potassium chloride. The price of diammonium phosphate is USD 445.1 per ton, with a nitrogen content of 18% and phosphorus content of 46%, while potassium chloride has a market price of USD 528.56 per ton and contains 60% potassium. The relevant literature on soil research in the upper reaches of the Min River [48] indicates that the soil in this area contains 0.29% N, 0.02% P, and 0.29% K.

The value of reducing sediment deposition is generally assessed using the shadow project method. The typical sediment transport pattern in China’s major rivers suggests that over one-third of soil erosion sediment remains at the source, roughly one-third reaches the sea, and approximately 24% is deposited in river channels, lakes, and reservoirs [47]:

$$V_k=\frac{Q}{\rho }\times F\times 0.24$$

(12)

where F stands for the cost of sediment removal in USD·m⁻³. Referring to pertinent studies [47,49], this removal is valued at 0.91 USD·m⁻³.

(4)

Ecosystem Service Value

By adding up the values of carbon storage, water conservation, and soil conservation, the ecosystem service value of the study area can be obtained. The calculation formula is as follows:

$$ESV=V_C+V_W+V_S={\displaystyle \sum }_{i=1}^n(V_C{}_i+V_W{}_i+V_S{}_i)$$

(13)

where V_i, V_j, and V_k represent the value of carbon storage, water conservation, and soil conservation in region i.

2.3.2. Analysis of Global Moran’s I

This study uses exploratory spatial analysis to determine the spatial agglomeration types and autocorrelation of ESV. The global Moran’s I index is used for global autocorrelation measurements and examinations, and its calculation formula is as follows [17,50]:

$$\begin{array}{cc}\mathrm{Moran}’\mathrm{s}& I\end{array}=\frac{n{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{\omega }_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{{{\displaystyle \sum }}_{j=1}^n\left(x_i-\overline{x}\right){{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{\omega }_{ij}}$$

(14)

where n is the number of grids; x_i and x_j are the attribute values of grid i and j, respectively; x is the average value of attributes; and w_ij is the spatial weight matrix. A Moran’s I range of [–1, 1] indicates spatially significant irrelevance when I = 0. This result corresponds to a positive correlation with space and represents a spatial feature aggregation when I > 0 and a negative correlation with space and represents a discrete spatial feature when I < 0.

Using a z-test to conduct a statistical test on Moran’s I, the spatial correlation of features is determined by comparing the values of p and z obtained via global autocorrelation. The criteria for judgment can be found in

Table 2. Level of confidence for global spatial autocorrelation.

Z Score	The Value of p	Level of Confidence
<−1.65 or >+1.65	<0.10	90%
<−1.96 or >+1.96	<0.05	95%
<−2.58 or >+2.58	<0.01	99%

[17].

2.3.3. Analysis of Getis-Ord Gi*

Getis-Ord Gi* can identify the high-value spatial agglomeration and low-value spatial agglomeration of regional spatial distribution, thereby reflecting the spatial distribution of Getis-Ord G* areas at the regional response level. The formula is as follows [17]:

$$\begin{array}{c}{G}_i^{*}\left(d\right)=\frac{{{\displaystyle \sum }}_{i=1}^n{w}_{ij}\left(d\right)X_i}{{{\displaystyle \sum }}_{i=1}^n{X}_i}\end{array}$$

(15)

where $X_i$ is the observed value of region unit i, and $w_{ij}$ is the spatial weight matrix. When the value of $G_i^{*}\left(d\right)$ is significantly positive, it indicates that the value near region i is relatively high, which is the region with the highest ESV. When the value $G_i^{*}\left(d\right)$ is significantly negative, it indicates that the value near the region i is relatively low, which is the area with the lowest ESV.

2.3.4. Geo-Detector

Geo-detectors can determine the impact of individual driving factors on ESV and can also explore the interplay between factors, including factor detection, risk detection, interplay detection, and ecological detection. Factor detectors can detect the influence of various potential influencing factors on the value of the ecosystem, with an influence value q ranging from 0 to 1. Here, a larger q value indicates a stronger explanatory power for the dependent variable [51]:

$$q=1-\frac{{{\displaystyle \sum }}_{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}$$

(16)

where h = 1, 2…; L is the classification or partition of variable (Y) or factor (X); N_h and N are layer h and the regional number units, respectively; and ${\sigma }_h^2$ and ${\sigma }^2$ are the variance of layer h and the regional value Y, respectively.

The variance of the regional value Y is calculated as follows:

$${\sigma }^2=\frac{{{\displaystyle \sum }}_{i=1}^n{(Y_i-\overline{Y})}^2}{N-1}$$

(17)

where Y_j and $\overline{Y}$ are the mean value of sample j and the region Y, respectively:

$${\sigma }^2=\frac{{{\displaystyle \sum }}_{i=1}^{n_h}{(Y_{h,i}-\overline{Y_h})}^2}{N_h-1}$$

(18)

where Y and $\overline{Y}$ are the value and mean value of sample i in layer h, respectively.

Interaction detection: it is used to identify the interactions between different impact factors Xs, that is, to evaluate whether the combined action of X₁ and X₂ will increase or weaken the explanatory power of vegetation coverage Y or whether the influence of these factors on Y occurs independent of each impact factor. The specific relationship between the two factors can be divided into the following five categories, as shown in

Table 3. Result types of the interplay between two factors.

The Basis for Judging	Result Types of the Interplay between Two Factors
q(x1∩x2) < min(q(x1), q(x2))	Nonlinear Weakening Trend
min(q(x1), q(x2)) < q(x1∩x2) < max(q(x1), q(x2))	Nonlinear Weakening Trend with One Factor
q(x1∩x2) > max(q(x1), q(x2))	Increasing Trend with Two Factors
q(x1∩x2) = q(x1) + q(x2)	No Mutual Interference
q(x1∩x2) > q(x1) + q(x2)	Nonlinear Increasing Trend

:

Next, we combine relevant research results [52,53,54], and progress is made to combine them with the actual situation of the study area. For the natural environment, we consider five factors: NDVI, elevation, slope, average annual temperature, and average annual rainfall. Socially and economically, the five factors of PD, urbanization rate (RA), GDP, Human Influence Index (HAI), and the share of primary industry output in GDP (SPIO-GDP) are considered. Data on slope and elevation were obtained from DEM. The sources of annual average rainfall, annual average temperature, NDVI, UR, GDP, SPIO-GDP, and PD are shown in Table 1. The HAI calculation method is based on the research in [55], and the driving factor data are reclassified by natural breaks (Figure 3).

3. Analysis and Results

3.1. Ecological Service Value in the Upper Reaches of the Minjiang River, China

Figure 4 shows values of carbon stock, water conservation, and soil conservation in the upper reaches of the Minjiang River from 2000 to 2020.

The value of carbon stock is lower in high-altitude areas in the north and west than the value of carbon stock in the Zagunao River Basin and Yuzi Creek Basin. The value of water conservation, however, is generally higher in the west and south than that on both sides of the river valley. Additionally, the distribution characteristics for the value of soil conservation are higher in the northwest corner and lower in the southeast corner. The key reason for this value distribution is that the northwest region is located at a high altitude with relatively flat terrains and no obvious undulation. However, the southeast region is located at a high altitude with steeper terrain, leading to frequent geological disasters. Here, the capability of soil conservation is very weak, decreasing the soil retention value, which is relatively low. From the perspective of spatial distribution, the spatial variations in the values of carbon stock, water, and soil conservation are quite apparent. This factor is related to the vertical distribution of land use and land cover with increasing altitude.

The spatial distribution of ESV (ecosystem service value) in the upper Minjiang River region exhibits a characteristic of being ‘lower in the northwest and higher in the southeast. The ecosystem service value of the Zhenjiangguan upstream and Xiaoxinggou basin is relatively low; however, the ecosystem service value was found to be relatively high in the Heishui River, Zagunao River Basin, and Yuzi Creek Basin, primarily due to the large forest cover in these areas, which has better effects on carbon stock, soil conservation, and the capacity of the soil to retain water. Simultaneously, there are significant differences in climatic conditions, altitude, and other factors between the northwest and southeast parts of the study area, which are crucial contributors to the uneven spatial distribution of the ESV.

From the perspective of changes in value (

Table 4. Changes in ecosystem service value in the study area.

Year	2000	2005	2010	2015	2020	Changes in 2000–2020
Carbon Stock/Billion USD	6.131	6.130	6.132	6.131	6.128	−0.0026
Water Conservation/Billion USD	1.631	1.837	1.910	1.697	2.120	0.489
Soil Conservation/Billion USD	5.275	5.297	7.308	6.703	8.867	3.592
ESV/Billion USD	13.036	13.264	15.350	14.530	17.114	4.078

), the total value of carbon stock in the upper reaches of the Minjiang River from 2000 to 2020 remained overall stable, decreasing by only USD 2.6 million. The primary reasons for this result are the increase in built-up land, bare land, and water areas in the region. The capacity of water conservation also continued to increase, with an annual increase in the water conservation value of USD 24.45 million. In addition, the value of soil conservation increased significantly, with an average annual increase of USD 179.6 million. However, the value of soil conservation was linked to high volatility. The ecosystem in the research area is fragile, making it easily susceptible to mountainous disasters and disturbances from human activities. From the perspective of ESV in the upper reaches of the Minjiang River, the ESV over the five years was USD 13.04 billion, USD 13.26 billion, USD 15.35 billion, USD 14.53 billion, and USD 17.11 billion, in that order. The ESV increased by 31.28% in 2020 compared to 2000, with an average annual increase of USD 203.9 million. The increase in the region’s ESV primarily stems from improvements to soil conservation functions, with soil conservation contributing up to 88.08%. Water conservation follows closely, accounting for nearly 12% of the contributions. Research indicates that the long-term implementation of policies, including ecological protection and the Project of Natural Forest Conservation, have greatly improved ecosystem services in the upper reaches of the Minjiang River, especially the capacity for soil and water conservation.

3.2. Analysis of the Spatiotemporal Distribution Characteristics of ESV

To reflect the variations in spatial changes in ESV in the study area on a smaller scale, the grid method was used to divide the study area. In total, 25,503 grids were created with grid cells of 1 km × 1 km to calculate the ecosystem service value in each grid. The Global Moran’s I was employed to analyze the spatial distribution characteristics of ESV in the study area. The Moran’s I index of ESV in 2000, 2005, 2010, 2015, and 2020 was 0.7168, 0.7178, 0.7243, 0.7141, and 0.7114, respectively. The Z scores were 55.9255, 54.8715, 56.5106, 55.7092, and 55.5003, respectively. When p < 0.001, the ESV in the upper Minjiang River region possesses a certain degree of spatial correlation, primarily characterized by “High-High” and “Low-Low” clustering patterns. Concurrently, as time progresses under the influences of socio-economic policy adjustments and various natural factors, the spatial aggregation of the ESV presents a weakening trend.

Getis-Ord Gi* is used to analyze ESV local autocorrelation in the study area (Figure 5). The spatial distribution characteristics of ESV from 2000 to 2020 feature no significant changes. The ESV is very high in the Heishui River basin, the Zagunao River basin, and the Mao County of the upper reaches of the Minjiang River–Dujiangyan City. Low-value areas are mainly clustered in the north at higher elevations, where the land type is dominated by grassland and bare ground. Judging from the position of ESV distribution, ESV high-value areas are mostly located in Wenchuan County, Li County, and Mao County, whose land use types are mainly woodland, with favorable hydrothermal conditions and abundant natural resources. Across the five time periods, the areal sizes for different spatial clustering patterns presently displayed the following sequence: non-significant > hotspot area > coldspot area. Both hotspot and coldspot areas have steadily increased. The proportion of the coldspot areas increased from 16.70% in 2000 to 17.53% in 2020, while the hotspot areas increased from 18.39% to 19.38%. This result indicates that the spatial aggregation characteristics of high and low ESV values in the upper Minjiang River region are continuously intensifying.

3.3. Driving Forces Analysis of ESV Spatial Variations

3.3.1. Single Factor Detection

By exploring the factors’ impacts on the spatial variations in ESV in the upper reaches of the Minjiang River, the values of q from different driving forces are obtained (

Table 5. Results of driving factor’s impacts on the spatial variations of ESV.

Factors		qVC Statistic	qVw Statistic	qVS Statistic	qESV Statistic
Natural factors	Slope (X1)	0.11	0.037	0.333	0.309
	Rainfall (X2)	0.099	0.229	0.045	0.04
	Temperature (X3)	0.189	0.402	0.29	0.11
	NDVI (X4)	0.18	0.214	0.1	0.058
	Elevation (X5)	0.203	0.385	0.31	0.13
Socio-economic factors	Urbanization rate (X6)	0.027	0.039	0.224	0.168
	GDP (X7)	—	0.032	0.171	0.093
	HAI (X8)	0.542	0.164	0.173	0.197
	SPIO-GDP (X9)	0.023	—	0.107	0.104
	PD (X10)	—	0.051	0.113	0.052

Note: “—” indicates that it fails the correlation test when p > 0.01.

).

Based on the values of q, it can be seen that the dominant factor affecting the spatial variations in the value of carbon stock is HAI, the dominant factor affecting the spatial variations in the value of water conservation is temperature, and the value of soil conservation is the slope. These are the effects of individual factors. The contribution of slope to ESV in the study area is 30.9%, indicating that slope is the dominant factor causing spatial variations in ESV. The contributions of HAI and urbanization to ESV are 19.7% and 16.8%, respectively, indicating that these factors have an important influence on regional ESV spatial differentiation. Therefore, in addition to natural factors, the intensity of human activities also caused greater disturbances to impacted ecosystems. In addition, elevation, temperature, and SPIO-GDP all feature contributions to the ESV of over 10%, making them secondary influencing factors behind spatial variations in ESV. The contribution of other factors is less than 10%, and their impacts on the spatial variations in ESV in the study area are relatively small. These results indicate that the spatial variations in ESV in the upper reaches of the Minjiang River are the result of a combination of natural and social factors.

3.3.2. Exploration of the Interplay between Factors Underlying Spatial Variations in ESV

As shown in Figure 6a–c, based on the results of spatial variation in carbon stock value, the interplay between HAI and other factors is significantly greater than the interplay between other factors, with the strongest driving forces being X₈∩X₅, with 0.582. X₃∩X₆ has the strongest driving influence on the value of water conservation, with 0.478. In addition, the interplay force between temperature and other factors is larger than that between other factors. In the results related to the value of soil conservation, X₁∩X₅ has the strongest interplay force, with 0.493. The interplay between slope and other factors is larger than that between other factors. Therefore, natural factors and social factors have the strongest impact on ESV, and the interplay between factors leads to different driving forces on ESV.

As shown in Figure 6d, the driving forces of all factors on ESV after interplay are greater than the driving force of a single factor, indicating a close relationship between the spatial variations in ESV in the upper reaches of the Minjiang River and all factors. X₁∩X₈, X₁∩X_5, and X₁∩X₃ have a strong influence on the driving force of ESV, which are 0.391, 0.387, and 0.371, respectively. In other words, this suggests that slope plays a leading role in the spatial variations in ESV in the upper reaches of the Minjiang River. Meanwhile, rainfall and PD have a weak driving force on the spatial variations in ESV, but their influence is significantly increased after being combined with slope, which further demonstrates the influence of slope on the spatial variations in ESV. However, the driving force of PD and GDP on ESV did not increase significantly after they combined these factors with other factors, showing that PD and GDP are not the main influencing factors behind the spatial variations in ESV in the upper reaches of the Minjiang River.

4. Discussion

This paper used the InVEST model to quantitatively assess the physical quantity and value of ecosystem carbon storage, water conservation, and soil retention in the region. We also used Global Moran’s I and Getis-Ord Gi* to analyze the characteristics of the spatial and temporal evolution of the ESV, as well as geographic probes to study the key drivers and magnitude of the driving forces underlying the spatial and temporal divergence of the ESV. Those efforts not only reduced biases in the results caused by subjective judgments but also provided support for research on the influencing mechanisms of various geographical factors on the spatial patterns of rural settlements. These results can provide a clear scientific basis for resource and environmental management and ecological compensation in the upper reaches of the Minjiang River and similar alpine valleys all over the world.

Judging from the spatial distribution characteristics of ESV, features in the upper reaches of the Minjiang River are “lower ESV in the north and higher ESV in the south”. Additionally, the ESV is low on both sides of the valley and in the high-altitude areas. Conversely, areas with higher ESVs are concentrated in the Heishui River and Mao County of the upper reaches of the Minjiang River–Dujiangyan City, and areas with lower ESV are mainly distributed in the northwest area of Songpan County. These results are generally consistent with those of previous studies [17,56,57] because the western and southern parts of the study area are mainly covered by forests and grasslands [56,58], which have stronger water conservation, carbon storage, and soil erosion prevention capabilities [55,59]. The steep slopes on both sides of the valley, with minimal vegetation cover and relatively severe soil erosion, contribute to the lower value of water conservation in this area. The total ESV has continued to increase in the past 20 years, especially in recent years. This increase indicates that soil erosion in the study area has been effectively contained and that the overall trend of ESV is enjoying sound development [60]. The value of soil conservation contributed most to the increase in ESV. However, the fluctuations in the value of soil conservation during the past 20 years were the largest. This result shows that a series of ecological restoration projects in the upper reaches of the Minjiang River basin have played a positive role in preventing soil erosion but that stability over the long term does not exist in the ecosystem. Therefore, relevant measures need to be continuously implemented [61].

In this paper, we selected 10 factors related to the natural environment and socio-economic aspects to study the drivers of ecosystem service value. There are differences in the leading factors for the spatial variations in carbon storage value, water conservation, and soil conservation. Slope remains the key factor leading to the spatial variations in ESV, followed by HAI and the urbanization rate. These results indicate that the spatial distribution of ecological service value in the upper reaches of the Minjiang River is affected by a combination of factors, including natural and social factors. Natural factors affect the spatial variation of ESV in the study area, and socio-economic factors will further change ecosystem services [62,63,64]. Therefore, enhancing ecological service value in this area should be guided by the law of geographical differentiation based on the ecological function positioning and adhering to relevant zoning and piecemeal promotion. Additionally, the results of observing the interactions of the various factors revealed that the driving force of ESV spatial variations between the two factors could be enhanced, with 66.67% of the factors multiplying the corresponding forces and 33.33% presenting nonlinearity. All drivers interacted with greater explanatory power than the explanatory power of a single factor, showing that a combination of natural and social factors significantly increases the value of ecological services. Therefore, in the process of ecological and environmental protection and risk control in the Minjiang River Basin, we must fully recognize the action mechanisms of various influencing factors, especially their mutual promotion and strengthening effects.

The changes in ESV in the upper reaches of the Minjiang River are closely related to forest land and grassland. The valley lowland is the main area occupied by cultivated land and construction land, and the main area has a lower ESV. Therefore, it is suggested to strengthen the control of ecological protection red lines in combination with regional territorial spatial planning and carrying out ecosystem restoration. In this context, efforts should be made to repair and protect the ecosystems in areas with severe ESV loss, thereby improving the service capabilities of the entire regional ecosystem [56]. In response to issues such as vegetation destruction, land degradation, and decreased wetlands in some regions, especially frequent natural disasters [37], it is necessary to adopt differentiated regulatory policies by exploiting ecological advantages. This measure can enable one to select the appropriate development mode of land use that matches the regional development level. By preventing unreasonable human interference, continuous progress can be attained to promote the construction of livelihood projects and ecological engineering and promote high-quality regional development [63,65].

By simulating changes in the quality and value of ecosystem services under different land cover types, InVEST can be used to spatially and visually express the quantitative assessment of ecosystem service function value, which has certain advantages over the equivalent factor method of Costanza [66] and Xie Gaodi [8]. However, InVEST is very sensitive to data changes, and the lack of verification for the simulation results in this paper will have a certain impact on the scientific nature of the conclusions. In the long run, it will be necessary to further verify the theoretical results, predict the development trends of ESV evolution, and further improve the accuracy of the research results.

5. Conclusions

In this study, the upper reaches of the Minjiang River were taken as the research area. Multi-source data and the InVEST model were applied to quantify the value of water conservation, carbon storage, and soil conservation. These values were obtained to comprehensively assess the ESV of this basin and analyze the characteristics of its spatiotemporal variations. After employing a geographic detector to explore the main driving forces leading to spatial variations in ESV, the results are as follows.

(1)

From 2010 to 2020, the ESV in the study area increased from USD 13 billion to USD 17 billion, showing a general upward trend. This increase in ESV is mainly due to improvements in the value of soil and water conservation, while the value of carbon stock contributes little to this phenomenon. ESV generally presents a spatial distribution trend of ‘low in the north and high in the south’. The ESV in the upstream areas of Zhenjiangguan and Xiaoxinggou watershed is low, while the ecosystem service value in the Heishui River, Zagunao River, and Yuzi Creek watershed is relatively high;

(2)

ESV features intense spatial clustering. The ESV is obviously high in three places: the Heishui River basin, the Zagunao River basin, and the Mao County in the upper reaches of the Minjiang River–Dujiangyan City area. There is also one area in the north of Songpan County and the northwest of Heishui County where the ESV is visibly low;

(3)

The spatial variations in ESV in the upper reaches of the Minjiang River are the result of both natural and social factors. Here, slope is the main factor causing spatial variations in ESV. The driving force of all factors on ESV after interplay is greater than that of a single factor, and the interplay between slope and HAI has the greatest impact on ESV. The key factors causing the spatial variations in ESV vary from place to place. Natural factors affect the spatial variation of ESV in the study area, and socio-economic factors will further change ecosystem services. Therefore, continuous progress should be made in the ecological protection and restoration of the upper reaches of the Minjiang River by enhancing land use management and the stability of the ecosystem.