ABSTRACT
This study employs the InVEST model to conduct a quantitative assessment of water-related ecosystem services in five geomorphological units of the Chishui River Basin, over the period 2000–2020. Results indicate that mean annual water yield increased from 537.6 to 667.8 mm, with high-yield zones migrating markedly downstream, while low-yield areas persist in the western mid-elevation highly undulating mountains (MEHUM). Soil conservation improved from 269.3 to 400.3 t a−1, featuring a ∼15% surge during 2005–2010 attributable to the Grain-for-Green program. Nitrogen export rose from 0.44 to 0.48 t a−1, exhibiting a ‘higher near-water, lower farther’ spatial pattern. Geodetector analysis reveals land use as the primary driver for water yield (WY) and nitrogen export (NE) (q > 0.60), whereas slope and topographic roughness explain 20%–40% of soil conservation (SC) variance. In the low-elevation hills (LEH) unit, precipitation accounts for 36.9% of WY variance, with LD ∩ PRE interaction boosting it to 82.7%; population density ∩ elevation interaction explains 69.7% of SC variance. Under dual-factor constraints in LEH, peak WY, SC, and NE reach 1,331.4 mm, 4,566.4 t a−1, and 2.04 t a−1, respectively. These findings underscore the amplified coupling of anthropogenic and climatic controls in low-altitude terrains, highlighting the necessity of terrain-specific, adaptive management within surveying and geospatial frameworks to sustain WES.
HIGHLIGHTS
The resolution of relevant parameters was improved, effectively estimating the WES of small watersheds with surface heterogeneity.
The Chishui River Basin was divided into five different geomorphological regions to analyze the driving patterns and constraint effect differences of WES under varying geomorphological conditions, providing more targeted basis for the potential enhancement and optimization of WES.
INTRODUCTION
Ecosystem services (ES) refer to the various functions and services provided by ecosystems that directly affect the sustainable development of society and human well-being (Shen et al. 2023). However, with the acceleration of urbanization, drastic changes in land use, and the continuous development of natural resources, the structure and processes of ecosystems have undergone significant alterations, leading to the degradation of ES and triggering a series of environmental problems (Kaczorowska et al. 2016). These changes have not only intensified ecological crises such as soil erosion and biodiversity loss but have also weakened the stability and adaptive capacity of regional ecosystems, posing a severe threat to ecological security (Pedrinho et al. 2024). In the context of current territorial ecological governance, achieving comprehensive protection of ecosystems requires a systematic analysis of the complex interactions between society and ecosystems, underpinned by the coordination of multiple ecological factors (Gao, Wang et al. 2024). Fully understanding the driving mechanisms of ES and exploring collaborative optimization pathways for enhancing these services are crucial for achieving sustainable regional development (Li et al. 2024). Therefore, scientifically quantifying ES, revealing their spatial patterns and driving factors, and establishing a reasonable ecological management system are of great practical significance for promoting the stability and sustainable utilization of ecosystem functions (Peng et al. 2023).
In research on water resource management and ecological protection, the hydrological processes and water-related ecosystem services (WES) within a watershed play a pivotal role. The intensification of climate change has caused significant alterations in precipitation patterns across spatial and temporal scales, profoundly affecting the reorganization of watershed hydrological networks and ecological functions. Previous studies have utilized optimal channel network (OCN) models to explore the evolution of river network structures under different precipitation gradients, finding that uneven precipitation leads to river migration, expansion or contraction of watershed areas, and the formation or disappearance of waterways and catchment areas (Elmdoust et al. 2016). This structural reorganization not only changes the hydromorphological characteristics but also reflects the deep impacts of precipitation changes on watershed hydrological responses. Moreover, identifying key nodes in river networks has become a hot topic in watershed hydrological studies, with researchers using algorithms to identify the most influential nodes in fractured networks, revealing their crucial role in sustaining ecosystems (Scheffer et al. 2001; Sarker et al. 2019).
WES primarily encompasses services such as water supply, water regulation, and water quality (WQ) purification (Brauman et al. 2007), arising from the interactions between aquatic systems and surrounding terrestrial ecosystems (Ren 2021), including water yield (WY), soil conservation (SC), and WQ (Valente et al. 2021). WES can have both positive effects, such as enhancing water supply, and negative effects, such as exacerbating soil erosion (Cao et al. 2021; Lan et al. 2021). Therefore, changes in WES are closely linked to the sustainable development of a watershed. Although some research has been conducted, there is still a lack of a unified definition and classification standard, with most studies focusing on the role of water resources in ES and how water resources can be optimized for ecological management (Schmalz et al. 2016).
In complex karst regions, the spatial heterogeneity and regional differentiation of WES are particularly evident. The strong heterogeneity in karst areas leads to complex driving mechanisms for WES, making it difficult to relate local-scale research results to regional-scale trends (Xiong & Chi 2015). Especially in karst areas with severe rocky desertification, ecosystem degradation is a serious issue that limits regional socio-economic development (Wang et al. 2004; Wang & Li 2007). Previous studies have primarily employed remote sensing and empirical model methods, such as the Revised Universal Soil Loss Equation (RUSLE) and the InVEST ecosystem service assessment model, to analyze the spatial patterns of ecosystems and their relationship with environmental factors (Zeng et al. 2017; Lang et al. 2018). The InVEST model, with its spatially explicit modeling capabilities, ability to evaluate multiple services comprehensively, and flexibility across scales, has distinct advantages in WES studies in karst regions (Zhou et al. 2025). Based on remote sensing and land use data (LD), the model can simulate the spatial distribution of ES and explore the interactions between natural and human factors, providing scientific support for understanding the complex ecological processes in karst regions. Unlike process models that rely on field observations (e.g., Soil and Water Assessment Tool (SWAT)), InVEST can conduct reliable ecological assessments even in data-sparse situations (Tang et al. 2020). Furthermore, InVEST can simulate the dynamic changes of WES under different management scenarios, providing quantitative support for watershed ecological restoration and sustainable development. However, research on the spatial heterogeneity and regional differences of WES in karst watersheds remains limited.
The driving mechanisms of WES vary significantly across scales, with driving factors changing across spatial and temporal dimensions (Scheffer et al. 2001; Bennett et al. 2009; Duraiappah et al. 2014). In karst mountainous areas, different geological backgrounds and topographical features lead to pronounced spatial differences, making it difficult to apply regional research conclusions to local areas (Cai 2009). For example, Ma & Zhang (2018) pointed out that there are significant differences in soil erosion patterns across different slopes, watersheds, and regional scales. Therefore, understanding the formation mechanisms of ES across different scales is a major challenge in current research. In particular, for example, Hou et al. (2018) found that the runoff and normalized difference vegetation index (NDVI) were positively correlated in the upper reaches of the Sancha River but negatively correlated in the lower reaches, indicating the need to consider topographical differences when analyzing the spatial variation of WES drivers (Schirpke et al. 2019; Sun et al. 2020).
Therefore, when studying WES, it is essential to consider the diversity of geomorphological features and their influence on ES. This is particularly important in complex terrains such as karst watersheds, where the spatial heterogeneity of ES is especially pronounced. To address this, the present study conducts geomorphological zoning of the Chishui River Basin, dividing it into five major geomorphological types: low-elevation hills (LEH), low-elevation slightly undulating mountains (LESUM), low-elevation moderately undulating mountains (LEMUM), mid-elevation moderately undulating mountains (MEMUM), and mid-elevation highly undulating mountains (MEHUM). Based on these geomorphological differences, the InVEST model is employed to quantitatively assess three core WES variables: WY, SC, and nitrogen export (NE). By analyzing the spatial distribution of WES across different geomorphological regions, this study reveals the varying effects of different geomorphological types on WES. In addition, the geographical detector (GD) method is used to analyze the key driving factors and their interactions in different geomorphological zones, aiming to explore the dominant factors influencing changes in WES. The main objective of this study is to uncover the spatial heterogeneity and driving mechanisms of WES in the Chishui River Basin, especially across its various geomorphological zones, by identifying key influencing factors and their interactions. The findings provide a scientific foundation for regional ecological protection and socio-economic coordination, offering theoretical support for watershed ecological management and policy-making, particularly in regions with complex terrain, and serve as a valuable reference for ecological management in similar areas.
DATA AND METHODS
Study area
The Chishui River is the only first-order tributary of the Yangtze River in its upper reaches that has not been dammed along its main course and maintains a natural flow regime (Figure 1). It originates in Zhenxiong County, Zhaotong City, Yunnan Province, and flows through four cities and 13 counties (cities and districts) in the provinces of Yunnan, Guizhou, and Sichuan. The coordinates of the river are 104°44′–106°59′E and 27°14′–28°50′N. The total length of the main stream is 436.5 km, and the total watershed area is 19,600 km2, with elevations ranging from 192 to 2,487 m. The entire basin falls under a subtropical monsoon climate, characterized by cold and dry winters and hot and humid summers. The average annual temperature ranges from 13.1 to 17.6 °C, with annual precipitation between 749 and 1,286 mm, and the average annual runoff is approximately 9.74 × 109 m3. The Chishui River gets its name from the high sediment content and reddish-yellow water color, as it flows through karst and Danxia landform areas. The river is also known for its good ecological status, beautiful environment, and suitability for brewing. It is home to important red cultural heritage sites such as the ‘Four Crossings of Chishui’ and the ‘Zunyi Conference.’ The Chishui River Basin is located within the Yangtze River Upper Reaches National Nature Reserve for Rare Fish Species, making regional WQ protection a top priority. Therefore, studying the WES of the Chishui River Basin is not only crucial for scientifically assessing the ecological function and value of the region but also provides essential support for achieving green development and ecological protection goals.
Research data
The input data for the WES assessment using the InVEST model mainly includes the digital elevation model (DEM), precipitation (PRE), potential evapotranspiration (PET), LD, soil database (HWSD), depth to bedrock (DTB), plant available water content (PAWC), biophysical parameter tables, and watershed data. Among these, DTB is derived from bedrock depth data, and PAWC is obtained from soil data. For the driving mechanism analysis, the selected factors are PRE, PET, slope, NDVI, land surface temperature (LST), DEM, temperature (TEM), population density (POP), rock outcrop index (RFI), and landscape pattern indices such as number of patches (NP), patch density (PD), largest patch index (LPI), Shannon's diversity index (SHDI), edge density (ED), and contagion index (CONTAG). Precipitation and temperature data are interpolated using the ANUSPLIN method based on meteorological station data from the Chishui River region. Finally, the resolution of the raster data for each factor is resampled to 30 m, and the projection coordinate system is standardized to WGS_1984_UTM_48N. The specific data and their sources are provided in Table 1.
Required data and sources
Required data . | Data source . | Resolution . |
---|---|---|
LD | https://zenodo.org/record/8176941 | 30 m |
PRE | Chishui River Meteorological Station Data | |
PET | https://www.earthdata.nasa.gov/ | 1 km |
TEM | Chishui River Meteorological Station Data | |
NDVI | Institute of Resources and Environmental Sciences, Chinese Academy of Sciences (IRES, CAS) Data Center | 30 m |
LST | 250 m | |
HWSD | Harmonized World Soil Database version 1.2 (HWSD) | 1 km |
DTB | http://globalchange.bnu.edu.cn/research/cdtb.jsp. | 30 m |
DEM | Geospatial Data Cloud Platform (https://www.gscloud.cn) | 30 m |
POP | https://www.worldpop.org/ | 1 km |
Required data . | Data source . | Resolution . |
---|---|---|
LD | https://zenodo.org/record/8176941 | 30 m |
PRE | Chishui River Meteorological Station Data | |
PET | https://www.earthdata.nasa.gov/ | 1 km |
TEM | Chishui River Meteorological Station Data | |
NDVI | Institute of Resources and Environmental Sciences, Chinese Academy of Sciences (IRES, CAS) Data Center | 30 m |
LST | 250 m | |
HWSD | Harmonized World Soil Database version 1.2 (HWSD) | 1 km |
DTB | http://globalchange.bnu.edu.cn/research/cdtb.jsp. | 30 m |
DEM | Geospatial Data Cloud Platform (https://www.gscloud.cn) | 30 m |
POP | https://www.worldpop.org/ | 1 km |
Research methods
Overview of the study area. Note: LEH represents low-elevation hills, LESUM represents low-elevation slightly undulating mountains, LEMUM represents low-elevation moderately undulating mountains, MEMUM represents mid-elevation moderately undulating mountains, and MEHUM represents mid-elevation highly undulating mountains.
Overview of the study area. Note: LEH represents low-elevation hills, LESUM represents low-elevation slightly undulating mountains, LEMUM represents low-elevation moderately undulating mountains, MEMUM represents mid-elevation moderately undulating mountains, and MEHUM represents mid-elevation highly undulating mountains.
Water yield





Soil conservation
Nitrogen export


ANUSPLIN data interpolation









Geographical detector



RESULTS
Spatiotemporal variation trends of WES in the Chishui River Basin
Spatial pattern of WES in the Chishui River Basin. Note: (a), (b), and (c) represent the spatial distribution maps of WES in 2000; (d), (e), and (f) represent the spatial distribution maps of WES in 2005; (g), (h), and (i) represent the spatial distribution maps of WES in 2010; (j), (k), and (l) represent the spatial distribution maps of WES in 2015; and (m), (n), and (o) represent the spatial distribution maps of WES in 2020.
Spatial pattern of WES in the Chishui River Basin. Note: (a), (b), and (c) represent the spatial distribution maps of WES in 2000; (d), (e), and (f) represent the spatial distribution maps of WES in 2005; (g), (h), and (i) represent the spatial distribution maps of WES in 2010; (j), (k), and (l) represent the spatial distribution maps of WES in 2015; and (m), (n), and (o) represent the spatial distribution maps of WES in 2020.
In terms of SC, high-value areas were predominantly found in LESUM, LEMUM, MEHUM, and downstream MEMUM regions, while low-value areas were mainly concentrated in LEH regions. SC increased significantly from 269.28 t/a in 2000 to 400.25 t/a in 2020. The spatial pattern of SC remained relatively stable during 2000–2005 and 2010–2020, but a notable increase was observed between 2005 and 2010, which is closely linked to the implementation of the ‘Grain-for-Green’ project in the basin. The predominance of SC in the upstream region can be attributed to steep slopes in mid-altitude moderately undulating areas, which facilitate deep-rooted vegetation growth, thereby enhancing soil retention. Additionally, the region's high forest-cover plays a crucial role, as forests and grasslands effectively intercept rainwater, slow down runoff, and stabilize the soil through root systems, significantly reducing soil erosion.
The spatial pattern of NE in the Chishui River Basin followed a ‘high near-water, low far from water’ distribution. High NE values were concentrated in downstream LEMUM, MEMUM, and MEHUM. These areas are key agricultural hotspots, where the intensive use of fertilizers contributes to excessive nitrogen input. Fertilizers, rich in nitrogen, are washed into rivers during rainfall and irrigation, significantly increasing nitrogen export. Moreover, the flat topography of downstream regions, combined with slow water flow and prolonged water retention times, hampers the dilution and removal of nitrogen in surface runoff, leading to greater accumulation and export of nitrogen in the river system. Over the past two decades, nitrogen export has increased from 0.44 t/a in 2000 to 0.48 t/a in 2020.
Driving mechanisms of WES in the Chishui River Basin
Univariate driving mechanism analysis
Explanatory power of major influencing factors on WY
. | Factor . | 2000 . | 2005 . | 2010 . | 2015 . | 2020 . |
---|---|---|---|---|---|---|
Chishui River Basin | PRE | 0.1365 | 0.2736 | 0.1263 | 0.2264 | 0.1313 |
LD | 0.8588 | 0.7090 | 0.8188 | 0.7091 | 0.8178 | |
NDVI | 0.0809 | 0.1299 | 0.1935 | 0.1559 | 0.2873 | |
LEH | PRE | 0.1073 | 0.3687 | 0.2277 | 0.4896 | 0.3081 |
LD | 0.9724 | 0.9390 | 0.9741 | 0.9236 | 0.9802 | |
NDVI | 0.0645 | 0.1554 | 0.3075 | 0.3436 | 0.3455 | |
LESUM | PRE | 0.2111 | 0.1430 | 0.3450 | 0.1711 | 0.1341 |
LD | 0.9802 | 0.9681 | 0.9779 | 0.9545 | 0.9753 | |
NDVI | 0.1568 | 0.0727 | 0.3513 | 0.1911 | 0.3521 |
. | Factor . | 2000 . | 2005 . | 2010 . | 2015 . | 2020 . |
---|---|---|---|---|---|---|
Chishui River Basin | PRE | 0.1365 | 0.2736 | 0.1263 | 0.2264 | 0.1313 |
LD | 0.8588 | 0.7090 | 0.8188 | 0.7091 | 0.8178 | |
NDVI | 0.0809 | 0.1299 | 0.1935 | 0.1559 | 0.2873 | |
LEH | PRE | 0.1073 | 0.3687 | 0.2277 | 0.4896 | 0.3081 |
LD | 0.9724 | 0.9390 | 0.9741 | 0.9236 | 0.9802 | |
NDVI | 0.0645 | 0.1554 | 0.3075 | 0.3436 | 0.3455 | |
LESUM | PRE | 0.2111 | 0.1430 | 0.3450 | 0.1711 | 0.1341 |
LD | 0.9802 | 0.9681 | 0.9779 | 0.9545 | 0.9753 | |
NDVI | 0.1568 | 0.0727 | 0.3513 | 0.1911 | 0.3521 |
Explanatory power of major influencing factors on SC
. | Factor . | 2000 . | 2005 . | 2010 . | 2015 . | 2020 . |
---|---|---|---|---|---|---|
Chishui River Basin | DEM | 0.0239 | 0.0362 | 0.0401 | 0.0286 | 0.0307 |
Slope | 0.1761 | 0.1539 | 0.1563 | 0.1906 | 0.1741 | |
TRI | 0.1670 | 0.1389 | 0.1546 | 0.1880 | 0.1777 | |
LEH | DEM | 0.1513 | 0.2146 | 0.2528 | 0.1668 | 0.2463 |
Slope | 0.3364 | 0.5308 | 0.4179 | 0.4395 | 0.4353 | |
TRI | 0.3164 | 0.4803 | 0.4353 | 0.4132 | 0.4646 | |
LESUM | DEM | 0.3015 | 0.2714 | 0.2590 | 0.3349 | 0.2044 |
Slope | 0.3281 | 0.2902 | 0.3868 | 0.2656 | 0.4075 | |
TRI | 0.2610 | 0.3272 | 0.3589 | 0.2514 | 0.3911 |
. | Factor . | 2000 . | 2005 . | 2010 . | 2015 . | 2020 . |
---|---|---|---|---|---|---|
Chishui River Basin | DEM | 0.0239 | 0.0362 | 0.0401 | 0.0286 | 0.0307 |
Slope | 0.1761 | 0.1539 | 0.1563 | 0.1906 | 0.1741 | |
TRI | 0.1670 | 0.1389 | 0.1546 | 0.1880 | 0.1777 | |
LEH | DEM | 0.1513 | 0.2146 | 0.2528 | 0.1668 | 0.2463 |
Slope | 0.3364 | 0.5308 | 0.4179 | 0.4395 | 0.4353 | |
TRI | 0.3164 | 0.4803 | 0.4353 | 0.4132 | 0.4646 | |
LESUM | DEM | 0.3015 | 0.2714 | 0.2590 | 0.3349 | 0.2044 |
Slope | 0.3281 | 0.2902 | 0.3868 | 0.2656 | 0.4075 | |
TRI | 0.2610 | 0.3272 | 0.3589 | 0.2514 | 0.3911 |
Explanatory power of major influencing factors on NE
. | Factor . | 2000 . | 2005 . | 2010 . | 2015 . | 2020 . |
---|---|---|---|---|---|---|
Chishui River Basin | LD | 0.6620 | 0.6801 | 0.6304 | 0.6752 | 0.6681 |
NDVI | 0.1293 | 0.2260 | 0.2264 | 0.2414 | 0.2932 | |
LEH | LD | 0.5776 | 0.6302 | 0.5763 | 0.6650 | 0.6429 |
NDVI | 0.0724 | 0.0982 | 0.2700 | 0.3382 | 0.3358 | |
LESUM | LD | 0.5422 | 0.7076 | 0.7087 | 0.5995 | 0.6447 |
NDVI | 0.0446 | 0.1172 | 0.3553 | 0.2550 | 0.2817 |
. | Factor . | 2000 . | 2005 . | 2010 . | 2015 . | 2020 . |
---|---|---|---|---|---|---|
Chishui River Basin | LD | 0.6620 | 0.6801 | 0.6304 | 0.6752 | 0.6681 |
NDVI | 0.1293 | 0.2260 | 0.2264 | 0.2414 | 0.2932 | |
LEH | LD | 0.5776 | 0.6302 | 0.5763 | 0.6650 | 0.6429 |
NDVI | 0.0724 | 0.0982 | 0.2700 | 0.3382 | 0.3358 | |
LESUM | LD | 0.5422 | 0.7076 | 0.7087 | 0.5995 | 0.6447 |
NDVI | 0.0446 | 0.1172 | 0.3553 | 0.2550 | 0.2817 |
Driving forces of individual factors on WES in the Chishui River Basin for different topographical types. Note: Specifically, (a, b, c) represent the Chishui River Basin; (d, e, f) represent the LEH; (g, h, i) represent the LESUM; (j, k, l) represent the LEMUM; (m, n, o) represent the MEMUM; (p, q, r) represent the MEHUM. X1: DEM, X2: Slope, X3: TRI, X4: LT, X5: PRE, X6: LD, X7: PET, X8: TEM, X9: LST, X10: NDVI, X11: POP, X12: NP, X13: PD, X14: LPI, X15: SHDI, X16: ED, X17: CONTAG, X18: RFR (rock fraction rate).
Driving forces of individual factors on WES in the Chishui River Basin for different topographical types. Note: Specifically, (a, b, c) represent the Chishui River Basin; (d, e, f) represent the LEH; (g, h, i) represent the LESUM; (j, k, l) represent the LEMUM; (m, n, o) represent the MEMUM; (p, q, r) represent the MEHUM. X1: DEM, X2: Slope, X3: TRI, X4: LT, X5: PRE, X6: LD, X7: PET, X8: TEM, X9: LST, X10: NDVI, X11: POP, X12: NP, X13: PD, X14: LPI, X15: SHDI, X16: ED, X17: CONTAG, X18: RFR (rock fraction rate).
Table 2 presents the temporal variations in the explanatory power of PRE, LD, and NDVI for the Chishui River Basin and its subregions (LEH and LESUM) from 2000 to 2020. In the Chishui River Basin, LD consistently exhibits high explanatory power, ranging from 0.7090 in 2005 to 0.8588 in 2000, indicating the strong influence of landscape structure on ES. PRE showed some fluctuations, with a noticeable peak in 2005 (0.2736) and a significant drop in 2010 (0.1263), suggesting the variable impacts of precipitation on ecosystem functioning across the years. NDVI steadily increased over time, from 0.0809 in 2000 to 0.2873 in 2020, reflecting the enhancement of vegetation cover and its increasing role in ES.
In the LEH subregion, LD maintained exceptionally high explanatory power, peaking at 0.9802 in 2020, reflecting the region's landscape complexity and stability. PRE also exhibited notable fluctuations, with a significant rise in 2015 (0.4896), while NDVI increased progressively, highlighting improved vegetation cover.
The LESUM subregion showed lower values for both PRE and NDVI, but LD remained high, indicating that landscape diversity plays a dominant role in ecosystem functioning in this region as well.
Overall, LD consistently demonstrates a strong influence, while NDVI shows increasing importance over time, particularly in more recent years.
Table 3 illustrates the explanatory power of topographic factors – DEM, slope, and topographic roughness index (TRI) – on WES across the Chishui River Basin and two representative geomorphic units (LEH and LESUM) from 2000 to 2020. At the basin scale, topographic factors generally showed low to moderate explanatory power, with DEM consistently weak (0.0239–0.0401), indicating that elevation alone exerts minimal constraint on ES at the regional level. Slope and TRI showed relatively stable influence (0.15–0.19), reflecting their modest but consistent control over hydrological and soil processes.
In contrast, the LEH region demonstrated substantially higher values, particularly for slope and TRI. Slope peaked at 0.5308 in 2005, and TRI consistently explained over 0.41 after 2010, suggesting that terrain heterogeneity significantly influences ecosystem processes in hilly areas, likely due to impacts on runoff, erosion, and vegetation distribution.
The LESUM region also exhibited strong topographic influence, especially from DEM and slope. DEM reached 0.3349 in 2015, and both slope and TRI maintained high explanatory power throughout the period. These results underscore the increasing importance of terrain complexity in smaller spatial units, highlighting the need to incorporate fine-scale topographic variation when managing ES in mountainous and karst environments.
Table 4 presents the explanatory power of LD and NDVI on WES in the Chishui River Basin, as well as in the LEH and LESUM subregions from 2000 to 2020. At the basin scale, LD consistently shows high explanatory power (ranging from 0.6304 to 0.6801), indicating its dominant role in shaping ES. This suggests that landscape composition and configuration – such as land use heterogeneity and connectivity – have a strong influence on ecological processes like water retention and nutrient cycling. Meanwhile, NDVI showed a moderate increasing trend from 0.1293 in 2000 to 0.2932 in 2020, reflecting vegetation's growing role in ecosystem service dynamics, likely due to reforestation and ecological restoration efforts.
In the LEH region, LD also maintained relatively high explanatory power, especially in 2015 (0.6650), affirming the structural complexity of landscapes as a key driver in hilly terrains. NDVI in this region increased significantly after 2010, peaking at 0.3382 in 2015, indicating improved vegetation cover contributed notably to ecosystem service enhancement.
The LESUM subregion exhibited similar trends. While LD remained the dominant factor, NDVI's explanatory power saw a sharp increase in 2010 (0.3553), underscoring the role of vegetation in small watershed-scale ecosystem functioning. These results highlight the synergistic influence of landscape pattern and vegetation condition on sustaining WES.
From Tables 2–4, it is evident that in the LEH and LESUM regions of the Chishui River Basin, the influence of various driving factors on WES is greater than that observed at the watershed scale. This can be attributed to the relatively homogeneous environmental conditions in these areas (e.g., soil depth, moisture distribution, and solar radiation), which amplify the impact of natural factors. For instance, PRE exhibits a more uniform spatial distribution in regions with low topographic relief, thereby exerting a stronger influence on WY. Additionally, flat terrains, which are more susceptible to development and land use changes, tend to concentrate agricultural expansion, urbanization, and infrastructure construction. Consequently, land use patterns and human activities exert a more direct and pronounced impact on WES in these regions.
Dual-factor driver study
The driving forces of single factors on WES in different geomorphological types of the Chishui River Basin. Note: (a) represents the entire watershed; (b) represents the LEH; and (c) represents the LESUM.
The driving forces of single factors on WES in different geomorphological types of the Chishui River Basin. Note: (a) represents the entire watershed; (b) represents the LEH; and (c) represents the LESUM.
In LEH terrain, the interaction effects between PRE and other factors are significantly stronger than at the basin scale, with notable interactions observed for PRE ∩ PET (0.473), PRE ∩ TEM (0.516), PRE ∩ LST (0.539), and PRE ∩ PD (0.614). Additionally, the interactions PET ∩ LST (0.495) and PET ∩ PD (0.529) exhibit substantially higher values than at the basin scale. The most pronounced interaction, excluding those involving LD, is LST ∩ PD (0.671). In LESUM, the interaction patterns largely align with those in LEH terrain, except for PET ∩ PD (0.829), which emerges as the most enhanced interaction.
For SC, the dominant interacting factors are slope ∩ DEM (0.234) and slope ∩ PD = 0.360. In both LEH and LESUM, interactions involving slope, PD, and other factors exhibit a marked increase. In LEH, the most notable interactions include slope ∩ DEM (0.573), slope ∩ TRI (0.573), slope ∩ LD (0.573), slope ∩ PRE (0.573), slope ∩ NDVI (0.573), slope ∩ PD (0.573), PD ∩ TRI (0.697), PD ∩ LT (0.497), and PD ∩ PRE (0.584). In LESUM, PD interactions become even more significant than in LEH, with key interactions including PD ∩ DEM (0.692), PD ∩ slope (0.607), PD ∩ LT (0.780), PD ∩ PRE (0.514), and PRE ∩ POP (0.535).
Similar to WY, as LD is the primary driver of NE, interactions between LD and other factors demonstrate varying degrees of enhancement, with the most significant being LD ∩ PD (0.652). Additionally, NDVI ∩ PD (0.569) shows a notable increase. In low-altitude hilly terrain, strong interactions are observed for DEM ∩ LD (0.629), DEM ∩ PD (0.659), and LST ∩ PD (0.671). In low-altitude undulating terrain, the interaction effects of PD with other factors surpass those in the LEH, particularly for DEM ∩ PD (0.760) and LD ∩ PD (0.740). Among these, DEM ∩ NDVI (0.701) exhibits an improvement exceeding 20% compared with the basin scale.
Overall, dual-factor interactions significantly enhance the explanatory power of WY, SC, and NE. Across different landforms, the interaction trends remain consistent with the basin scale, where WY and NE are predominantly influenced by interactions involving LD, although the role of PD should not be overlooked. The interaction effects of PD with other factors are considerably stronger than those of individual factors alone. For SC, while the interactions between slope and other factors show moderate improvement, the most pronounced enhancements remain those between PD and other factors. These findings suggest that interactions between dominant drivers and secondary factors further amplify their explanatory power. To alter future patterns of WES, modifications to land use patterns are essential. Moreover, the strong influence of PD in combination with other factors underscores the critical role of human activities in shaping WES, highlighting the joint effects of anthropogenic and natural processes.
Analysis of the dual-factor constraint effect
Dual-factor constraint effects on WY in different geomorphological types of the Chishui River Basin. Note: (a) represents the Chishui River Basin, (b) represents the LEH, and (c) represents the MEMUM.
Dual-factor constraint effects on WY in different geomorphological types of the Chishui River Basin. Note: (a) represents the Chishui River Basin, (b) represents the LEH, and (c) represents the MEMUM.
Double-factor constraint effects on SC in different geomorphic types of the Chishui River Basin. Note: (a) represents the entire Chishui River Basin, (b) and (c) represent the LESUM, and (d) and (e) represent the LEMUM.
Double-factor constraint effects on SC in different geomorphic types of the Chishui River Basin. Note: (a) represents the entire Chishui River Basin, (b) and (c) represent the LESUM, and (d) and (e) represent the LEMUM.
Double-factor constraint effects on nitrogen output (NE) in different geomorphic types of the Chishui River Basin. Note: (a, b, c) represent the Chishui River Basin, and (d, e) represent the LEH.
Double-factor constraint effects on nitrogen output (NE) in different geomorphic types of the Chishui River Basin. Note: (a, b, c) represent the Chishui River Basin, and (d, e) represent the LEH.
Based on the constraint results (Figure 7), at the full basin scale, the interaction between DEM and slope and the interaction between slope and TRI are identified as the primary factors influencing SC. Under the combined influence of these factor interactions, SC shows significant potential for improvement. In the LEH, the interaction between slope and TRI (slope = 32.76, TRI = 476.00, SC = 4,566.37) is the most effective in enhancing SC, with its maximum value exceeding that of the entire basin scale. In the LESUM the interaction between slope and TRI exerts a much stronger constraint effect on SC than at the basin scale. Additionally, the interaction between slope and PRE also contributes to higher SC values. To alter SC in the Chishui River Basin, it is essential to comprehensively consider the combined effects of DEM, slope, and TRI, while also giving due attention to the influence of PRE.
DEM, LD, NDVI, and PD are critical factors influencing nitrogen export, each playing a significant role in different geomorphic regions. At the basin scale, the interactions between the following pairs of factors (Figure 8(a)–8(c)) – DEM ∩ LD (DEM = 220.00, LD = 1, NE = 2.06), LD ∩ NDVI (LD = 8, NDVI = 0.86, NE = 2.12), and NDVI ∩ PD (NDVI = 0.66, PD = 14.71, NE = 1.74) – can lead to peak NE.
In the LEH (Figure 8(d)–8(e)), the interaction between LD ∩ NDVI (LD = 1, NDVI = 0.83, NE = 2.04) and NDVI ∩ PD (NDVI = 0.83, PD = 14.76, NE = 1.69) results in the optimal nitrogen export. Compared with the entire basin, peak nitrogen output in the low-altitude hilly region is smaller and more easily reached under these interaction conditions.
The primary sources of nitrogen output include agricultural activities (such as fertilizer application and livestock waste), industrial wastewater discharge, domestic sewage, and nitrogen deposition from atmospheric settling. Nitrogen from these sources enters water bodies and soil through runoff, discharge, and deposition, thereby impacting the nitrogen cycle and ecosystem quality. Therefore, policies such as converting farmland to forest, reducing agricultural land, and enforcing stricter standards for industrial wastewater discharge are essential to reduce nitrogen output and achieve environmental protection.
DISCUSSION
Differences in driving factors between the Chishui River Basin and different geomorphic regions
In past ecosystem service research, many scholars have focused on exploring the impacts of factors such as climate, topography, land use, and human activities on ES. Existing studies are often based on large-scale watersheds or ecosystems and seldom address small-scale watersheds with strong surface heterogeneity. Research on larger-scale watersheds usually assumes that the driving factors of ES are spatially homogeneous, neglecting the effects of topographic differences. Especially in small-scale watersheds, the interactions between topography, climate change, and human activities lead to complex spatial heterogeneity in ecological processes, which has not been sufficiently addressed. For example, Huang et al. (2023) studied the driving mechanisms of ES in large-scale watersheds but did not delve into the effects of surface heterogeneity on the spatial distribution and driving mechanisms of ES in small-scale watersheds. Xue et al. (2024) focused on ES at a larger scale but similarly ignored the influence of topographic differences on the spatial distribution of ES. Dee et al. (2025) discussed the impact of climate change on ES but did not consider the interaction between topography and human activities in small-scale watersheds. Additionally, Zhang et al. (2023) pointed out that while climate factors have a significant impact at larger scales, the interaction between topography and human activity factors may play a more important role in small-scale watersheds. Gao, Song et al. (2024) also noted that the constraint effects of human activities on ES in small-scale watersheds deserve further investigation, but their studies mainly focus on single-factor effects, neglecting the interactions of multiple factors. Therefore, there is a gap in the existing research regarding the spatial distribution and driving mechanisms of ES in small-scale watersheds. This study aims to address this gap by focusing on the Chishui River Basin, a typical small-scale watershed, exploring the role of topographic heterogeneity and its interaction with climate and human activities in driving WES, thus filling the research gap in small-scale watershed studies.
This study takes the Chishui River Basin as a representative small-scale watershed and, through geomorphic regionalization, systematically analyzes how climate, topography, and human activities drive WES. Results show that PRE and LD are the dominant drivers of WY, SC, and NE, especially in LEH and LESUM, underscoring the dominant role of climate in low-relief areas (Forootan & Sadeghi 2023). Furthermore, interactions between slope and TRI enhanced SC capacity, while in the LESUM, the interaction between PD and natural factors became increasingly significant, highlighting the role of human activities (Gopalakrishnan et al. 2016; Watson et al. 2019). These findings indicate that the driving mechanisms of ES vary across geomorphic zones, and management strategies should be tailored to regional characteristics.
Constraint effect analysis further revealed the amplifying impact of multi-factor interactions on ES. In the LEH region, the combined effect of precipitation and land use resulted in a maximum WY of 1,331.37 mm, far exceeding the watershed-scale maximum of 1,041.65 mm, illustrating the synergistic effect of land use change and climatic factors (Bai et al. 2019). In addition, the combination of slope and terrain relief significantly boosted SC capacity (4,566.37 t/a) in this region. Meanwhile, the interaction between land use and vegetation index effectively reduced nitrogen export in both LEH and LESUM zones (Metzger et al. 2017). These results suggest that ecosystem management in small-scale watersheds should emphasize the synergistic effects of multiple factors, promoting spatially targeted strategies for water regulation, SC, and WQ protection. This provides both theoretical insight and practical guidance for future ecosystem service research.
In conclusion, this study, through a detailed analysis of different topographic regions of the Chishui River Basin, demonstrates the key role of topographic heterogeneity and multi-factor interactions in small-scale watershed ES, providing specific recommendations for regional ecological and environmental management. Unlike traditional large-scale watershed studies, this paper not only reveals the individual effects of climate and human activities but also deeply explores the complex interactions between various factors through multi-factor interaction analysis. It suggests that watershed management should comprehensively consider the roles of topography, climate, and human activities, providing a new theoretical framework for ecosystem service research in small-scale watersheds.
Uncertainty analysis and future prospects
This study is of significant importance in the fields of environmental protection and climate change. By quantifying the WES of the Chishui River Basin, including WY, SC, and NE, the research reveals the impact of various environmental factors on the basin's hydrological processes. It provides a scientific basis for water resource management, SC, and pollution control, and contributes to optimizing land use planning to reduce the impact of extreme climate events on the basin's ecosystem. Additionally, the study uses the InVEST model to quantitatively assess WES and employs the geographic detector to analyze their driving mechanisms and constraint effects. This not only enhances the accuracy of ecosystem service assessments at the watershed scale but also offers insights for the sustainable management of other karst regions. However, this study has certain limitations. For instance, the input data for the InVEST model (such as PRE, soil, and LD) has inherent uncertainties, which may affect the accuracy of the WES estimates. Furthermore, the study primarily focuses on natural factors and does not delve into the influence of socio-economic factors (such as policy interventions and human activities) on WES. Additionally, the research relies on static data and does not adequately consider the dynamic response of WES to future climate change scenarios, which limits its ability to predict long-term trends. These limitations may affect the generalizability of the findings, reducing their applicability in different climate and social contexts. Therefore, future research could improve the accuracy of WES assessments by integrating high-resolution remote-sensing data and long-term monitoring data, while also incorporating socio-economic variables to comprehensively explore the effects of human activities. Moreover, by combining climate models to predict future changes in precipitation and temperature, WES dynamics under various climate scenarios can be assessed, and the reliability of simulation results can be enhanced through multi-model comparison approaches (such as SWAT and LPJ-GUESS). These improvements will help provide a more comprehensive understanding of the evolution of ES and offer scientific support for regional ecological protection and global climate change mitigation.
CONCLUSION
This study focuses on the driving mechanisms and constraint effects of WES in different geomorphological regions of the Chishui River Basin. Using the Geodetector model, single-factor, two-factor explanatory power, and constraint effect analyses, it quantitatively reveals the spatial heterogeneity of the interactions between natural factors and human activities. The results show:
LD is the primary single-factor driver of WES, with an average explanatory power of 0.7627 for WY between 2000 and 2020, especially in LESUM, where the q-value reaches as high as 0.9802. Furthermore, the explanatory power of NDVI shows an increasing trend, rising from 0.0645 in 2000 to 0.3455 in 2020 in the LEH region, indicating that vegetation recovery has had an increasingly positive impact on WES supply. Terrain factors such as slope and TRI show significantly higher explanatory power in the LEH region compared with the overall basin average, with slope having a q-value of 0.5308 and TRI 0.4803, reflecting the important regulatory role of terrain in services like SC.
Two-factor interaction analysis shows that the interaction between LD and PRE in the LEH region has the highest explanatory power (0.9791 in 2015), significantly higher than the single factors of LD (0.9236) and PRE (0.4896), indicating that the synergy between climate and human activities more effectively explains the spatial distribution of ES. Constraint effect analysis shows that the maximum WY in LEH reaches 1,331.37 mm, much higher than the basin-wide maximum of 1,041.65 mm, and is mainly driven by the interaction between LD and PRE. Additionally, the interaction between slope and TRI in LEH results in an SC value of 4,566.37 t/a and the lowest nitrogen output of 2.04 t/a, emphasizing the importance of rational land use and vegetation recovery in controlling non-point source pollution and enhancing ES.
In conclusion, the driving mechanisms of WES in the Chishui River Basin exhibit significant regional differences and complex interactions. The quantitative analysis results indicate that zoning management strategies should be developed based on geomorphological characteristics to improve land use efficiency and promote ecological restoration, thereby enhancing the spatial resilience and supply capacity of ES.
ACKNOWLEDGEMENTS
This research was supported by the Geological Research Project of Bureau of Geology and Mineral Exploration and Development Guizhou Province (Qian Di Kuang Ke He (2020) No.27), the Guizhou Provincial Science and Technology Project (no. Qian Ke He Zhi Cheng [2022] General 199), and the Guizhou Provincial Science and Technology Project (no. Qian Ke He Zhi Cheng [2023] General 169). Acknowledgment for the data support from Resource and Environment Science and Data Center (https://www.resdc.cn/), National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/), and National Earth System Science Data Center (http://www.geodata.cn/).
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.
CONFLICT OF INTEREST
The authors declare there is no conflict.