Quantitative Analysis about the Spatial Heterogeneity of Water Conservation Services Function Using a Space–Time Cube Constructed Based on Ecosystem and Soil Types

Abstract

Precisely delineating the spatiotemporal heterogeneity of water conservation services function (WCF) holds paramount importance for watershed management. However, the existing assessment techniques exhibit common limitations, such as utilizing only multi-year average values for spatial changes and relying solely on the spatial average values for temporal changes. Moreover, traditional research does not encompass all WCF values at each time step and spatial grid, hindering quantitative analysis of spatial heterogeneity in WCF. This study addresses these limitations by utilizing an improved water balance model based on ecosystem type and soil type (ESM-WBM) and employing the EFAST and Sobol’ method for parameter sensitivity analysis. Furthermore, a space–time cube of WCF, constructed using remote-sensing data, is further explored by Emerging Hot Spot Analysis for the expression of WCF spatial heterogeneity. Additionally, this study investigates the impact of two core parameters: neighborhood distance and spatial relationship conceptualization type. The results reveal that (1) the ESM-WBM model demonstrates high sensitivity toward ecosystem types and soil data, facilitating the accurate assessment of the impacts of ecosystem and soil pattern alterations on WCF; (2) the EHSA categorizes WCF into 17 patterns, which in turn allows for adjustments to ecological compensation policies in related areas based on each pattern; and (3) neighborhood distance and the type of spatial relationships conceptualization significantly impacts the results of EHSA. In conclusion, this study offers references for analyzing the spatial heterogeneity of WCF, providing a theoretical foundation for regional water resource management and ecological restoration policies with tailored strategies.

1. Introduction

Water conservation service function (WCF) falls within the research scope of eco-hydrology and fundamentally signifies that the biotic and abiotic elements that constitute an ecosystem in a particular environment can offer tangible services to humanity as natural resources [1]. Among the diverse array of ecosystem services, water conservation stands out as a pivotal service function of regional ecosystems [2], exerting a direct influence on climate, hydrology, vegetation, productivity, and soil nutrient cycling within the region and its surrounding areas [3]. Accurately assessing and analyzing this service is particularly important in areas with significant disparities in natural resources and ecosystem service distribution.

With the advancement of comprehensive assessment models and ongoing research in WCF, scholars are increasingly inclined to employ more intricate hydrological models or integrate multiple models for related studies. Prominent hydrological models include the SCS model [4,5,6], the TVDI model [7,8], the SWAT model [9,10,11], and the SWIM model [12]. Likewise, common methods for integrating multiple models include the InVEST model [13], the WaSSI-C model [14], the BEPS-TerrainLab model [15,16], and cellular automata [17]. These approaches enable more comprehensive assessments and evaluations, contingent on the availability of sufficient regional data. However, in practical applications, conducting relevant research encounters challenges often due to inadequate foundational data and extensive research areas. Consequently, water balance models are extensively utilized to guide ecological protection policies at the national and regional levels, owing to their applicability across spatiotemporal scales and their capacity to comprehensively analyze the input (precipitation) and output (runoff and evapotranspiration) of the study area [18,19]. Moreover, this model has become a widely adopted general assessment method, implemented through standards and specifications in China’s National Ecological Status Assessment, China’s Ecological Red-line Delineation, and other regular remote sensing surveys and land spatial planning [20,21].

In the application of water balance models, a significant challenge lies in obtaining accurate runoff data that align with the data format, spatial scale, and resolution. Current studies [20] primarily estimate the spatial distribution of runoff depth by employing runoff coefficients specific to different ecosystem types. Traditional water balance models typically represent the runoff depth of a region as the product of a runoff coefficient and precipitation. However, this approach gives rise to three issues. Firstly, the absence of actual hydrological observation data hinders the representation of the ecosystem’s true water regulation capacity in the results. Secondly, the runoff coefficient used in the traditional method fails to account for watershed geomorphology and physical–geographical characteristics, thereby obscuring variations in the runoff production capacity among watersheds of the same ecosystem type. Lastly, the traditional approach assumes a fixed constant for the runoff coefficient, disregarding changes in vegetation growth and ecosystem type. To address these uncertainties, Zhai et al. proposed a runoff depth calculation method in 2023 based on measured runoff data from hydrological stations within specific watersheds. This method has been proven to provide a more accurate description of runoff depth at the watershed scale.

Additionally, there are limitations in the specifications of relevant regular survey assessments and general scientific research on the assessment of WCF. These limitations include using only the multi-year average value to describe spatial changes and relying solely on the spatial average value to study temporal changes. In addition, this research method does not encompass all WCF values at each time step and spatial grid. To overcome these limitations, Emerging Hot Spot Analysis (EHSA) is considered a viable approach. Currently, EHSA has been utilized in studying temporal and spatial changes in diverse fields, such as surface evapotranspiration rate [22], prediction of hydrological drought risk [23], and surface deformation [24]. Furthermore, Liu [25] investigated the spatial heterogeneity of WCF in the Yangtze River Basin using EHSA. However, the aforementioned studies solely employed this method for basic interpretation and analysis, and the impact of associated parameter factors on the output results of EHSA remains unclear. The principle of EHSA typically encompasses neighborhood distance and the conceptualization of spatial relationships. The configuration of neighborhood distance and the choice of spatial relationship conceptualization should precisely depict the distribution characteristics of the parameters. Importantly, unlike vector or drought index data studied in other disciplines, WCF is a composite evaluation index that integrates various factors, including meteorology, hydrology, and geography. Its distribution pattern is influenced by multiple conditions such as ecosystem type, soil type, precipitation, and evaporation. Its distribution pattern and heterogeneity differ significantly from those of other data types. Thus, further exploration of the influence of neighborhood distance and spatial relationship conceptualization on the representation of WCF spatial heterogeneity using EHSA is necessary.

This study utilizes observed data from hydrological stations to analyze the region’s WCF. To achieve this, a water balance model improved based on ecosystem type and soil type (ESM-WBM) is introduced. The model assumes a closed system, treating the forest as a black box and focusing solely on the input and output of water. By evaluating the differences between water input and consumption within the study area, the overall water-saving capacity of the region is determined. At the same time, EFAST and Sobol’ method for area are employed for parameter sensitivity analysis. Furthermore, a space–time cube of WCF is constructed using remote sensing data of ecosystem and soil types. The space–time cube is further explored using EHSA (Emerging Hot Spot Analysis). Additionally, this study investigates the impact of two core parameters, neighborhood distance and spatial relationship conceptualization type, on the expression of spatial heterogeneity in WCF.

2. Materials and Methods

2.1. Research Area

The study area is located within the Yellow River Basin (also named Huang He) and encompasses a total of 58 tributaries (Figure 1). Geographically, it lies at the convergence of Shaanxi, Shanxi, and Inner Mongolia provinces, situated in the central part of the Loess Plateau. The area spans approximately 121,241.92 square kilometers with an average elevation of 1227.60 m. In terms of vegetation coverage, the region primarily consists of grassland (45.69%) and cultivated land (28.10%). Additionally, forest land constitutes 13.60%, while urban and rural residential areas, as well as construction land, account for 2.34%. The study zone’s entry point in the Yellow River Basin is the Toudaoguai station, while the exit point is the Longmen station, both serving as national-level hydrological observation sites in China. Situated on the Loess Plateau, the region boasts an average elevation of 1227.60 m and is predominantly characterized by meadows and forests. Serving as a critical transitional area between the Yellow River’s upstream river source and downstream water demand zones, the region plays a pivotal role in the water supply–demand dynamics between these regions. The northern segment is comprised of residential and urban areas, agricultural land for irrigation, and limited industrial zones, while the southern part is distinguished by extensive grasslands and forests. The ecological divergence between the northern and southern areas renders the region highly conducive for research purposes. After combining existing boundary data, satellite remote-sensing imagery with a resolution of 2 m or higher was employed and underwent manual verification and calibration processes to ensure precision.

2.2. Data Sources

Based on various data sources, this study classifies the required data into five main types: hydro-meteorological station measurements, meteorological data, hydrological data, remote sensing images, and hydrological soil group data (

Table 1. Data types and sources used in this study.

Category	Spatial Data	Description	Data Sources	Data Source Websites
Hydro-meteorological Station Measurements	Measured runoff	Value	China Gazette of River Sedimentation	http://www.mwr.gov.cn/sj/tjgb/zghlnsgb/ (accessed on 10 December 2023)
	Measured precipitation	Value	China Meteorological Data Service Centre	http://data.cma.cn/data/detail/dataCode/A.0053.0002.S007.html (accessed on 11 December 2023)
Meteorological	MODIS Evapotranspiration product MOD16A2	Raster 500 m (2012–2022)	Level 1 and Atmosphere Archive and Distribution System DAAC (LAADS DAAC)	https://www.earthdata.nasa.gov/eosdis/daacs/laads (accessed on 14 December 2023)
	Monthly precipitation	Raster 1 km (2012–2022)	National Tibetan Plateau Data Center Third Pole Environment Data Center TPDC	https://data.tpdc.ac.cn/zh-hans/data/faae7605-a0f2-4d18-b28f-5cee413766a2 (accessed on 12 December 2023)
Hydrological	Runoff efficient	Value	Existing Research [26]	https://doi.org/10.3390/w15081475 (accessed on 11 December 2023)
Remote Sensing	DEM, digital elevation model	Raster 250 m	Spatial distribution data of multi-period ecosystem types in China	https://www.resdc.cn/ (accessed on 10 December 2023)
	Ecosystem Type	Raster 16 m	Landsat series satellite images	https://www.resdc.cn/ (accessed on 13 December 2023)
Soil	Hydrological Soil Group and Soil Data	Raster 1 km	Harmonized World Soil Database ver 1.2 by Food and Agriculture of the United Nations	https://www.fao.org/soils-portal/data-hub/soil-maps-and-databases/harmonized-world-soil-database-v12/en/ (accessed on 13 December 2023)

). The hydro-meteorological station measurements encompass measured evapotranspiration, runoff, and precipitation data. Meteorological data consist of MODIS evapotranspiration products (MOD16A2) and monthly precipitation datasets with a 1 km resolution for China. Runoff data includes runoff coefficient data for the study area. Land remote-sensing images mainly pertain to the ecosystem-type data derived from satellite remote-sensing images for the study area. Hydrological soil data refers to soil cover types and corresponding runoff curve number (CN) values extracted from the HWSD Ver1.2 global soil database for the study area.

2.3. Methods

The primary framework of this study comprises four main components (Figure 2): data collection, calculation of WCF, variation trend and spatiotemporal heterogeneity analysis, and driver factor analysis. The detailed methods are outlined below. The construction of space–time cubes and EHSA analysis are conducted using Python code. The EFAST and Sobol’ sensitivity analyses are performed using Matlab R2016a and the SAlib. All other analyses were conducted using the ArcGIS Pro 3.02 platform.

2.3.1. Runoff Depth Calculation Method Based on Ecosystem Type and Soil Type (ESM)

This study adopts a runoff depth calculation method based on ecosystem and soil type (ESM). In contrast to the lookup table method employed in traditional water balance approaches, this calculation method takes into account the ecosystem type and soil type of the study area. As a result, the comprehensive runoff index of ecosystem type and soil (R_ist) determined using this method can more accurately portray the actual eco-hydrological conditions of the study area, thereby offering a more precise estimation of the water conservation capacity of the study area. The detailed calculation method is outlined below:

$$R={\displaystyle \frac{\mathsf{\Delta }R}{S}}\times R_{ist}(j,K)$$

(1)

In Formula (1), R is the runoff depth of the YBR study region (mm), ΔR is the study basin’s variation of total runoff (mm), R_ist(j, K) is the comprehensive runoff index of ecosystem type j of HSG classification K, and S is the sum of all pixel values of the R_ist(j, K) gride data.

$$\mathsf{\Delta }R=R(outlet)-R(inlet)$$

(2)

Regarding the processing of hydrological station data, the runoff data used in this study is from the China Gazette of River Sedimentation [27]. This report encompasses the annual observed values of the effluent hydrological stations within the basin for the period spanning 2012 to 2022. The alteration in total runoff (ΔR) is derived by computing the disparity between the annual runoff recorded at the outlet hydrological station and the inlet station, as indicated in Formula (2).

$$\begin{array}{cc}& \mathrm{HSG}\mathrm{Classification}\\ \mathrm{EcosystemType}& \overline{\begin{array}{cc}& \begin{array}{ccc}{R}_{is}(HSG1)& \cdots & R_{is}(HSGK)\end{array}\\ \begin{array}{c}{R}_{it}(Type1)\\ ⋮\\ R_{it}(Typej)\end{array}& Rist(j,K)\end{array}}\end{array}$$

(3)

$$Rist(j,K)=0.5\times R_{it}(j)+0.5\times R_{is}(j,K)$$

(4)

$$R_{it}(j)={\displaystyle \frac{R_c(j)}{R_c(\min )}}$$

(5)

$$R_{is}(j,K)={\displaystyle \frac{CN(j,K)}{CN(\min )}}$$

(6)

Formulas (3) and (4) shows the calculation principle of the comprehensive runoff index of ecosystem type and soil, and R_ist (j, K), which adopts an equal weight manner to add both of the runoff indexes. A schematic diagram of the logical principle is presented in Figure 3. In Formula (5), R_it (j) is the runoff index of ecosystem type j, R_c(j) is the runoff coefficient of type j, and R_c(min) is the minimum value of the runoff coefficient of all ecosystem types in the study region. In the same way, in Formula (6), R_is(j, K) is the runoff index of the value of ecosystem type j of HSG classification K. CN(j, K) is the CN value of ecosystem type j of HSG classification K. CN(min) is the minimum value of them (

Table 2. Runoff coefficients, index, and CN of different ecosystem types.

Ecosystem Types	Runoff Parameters		CN of Different HSG and Runoff Index of Soil (Ris)
	Runoff Coefficient (Rc)	Runoff Index (Rit)	CNA	CNB	CNC	CND	Ris(A)	Ris(B)	Ris(C)	Ris(D)
Broad-leaved forest	3.33	1.55	36	62	75	81	1	1.72	2.08	2.25
Coniferous forest	2.15	1.00	37	62	75	81	1	1.68	2.03	2.19
Mixed forest	2.4	1.12	38	62	75	81	1	1.63	1.97	2.13
Sparse forest	16.02	7.45	72	82	83	87	1	1.14	1.15	1.21
Broadleaf shrub	3.58	1.67	45	65	75	80	1	1.44	1.67	1.78
Coniferous shrub	3.41	1.59	49	69	79	84	1	1.41	1.61	1.71
Open shrubland	16.02	7.45	72	82	83	87	1	1.14	1.15	1.21
Marshy grassland	9.11	4.24	49	69	79	84	1	1.41	1.61	1.71
Steppe	12.34	5.74	49	69	79	84	1	1.41	1.61	1.71
Tussock	10.18	4.73	49	69	79	84	1	1.41	1.61	1.71
Sparse grassland	16.02	7.45	72	82	83	87	1	1.14	1.15	1.21
Marsh	0	0	0	0	0	0	0	0	0	0
Lakes	0	0	0	0	0	0	0	0	0	0
River	0	0	0	0	0	0	0	0	0	0
Farmland	49.69	23.11	67	78	85	89	1	1.16	1.27	1.33
Plantation	4.62	2.15	52	69	79	84	1	1.33	1.52	1.62
Settlement	90	41.86	80	85	90	95	1	1.06	1.13	1.19
Urban green space	7.91	3.68	52	70	79	84	1	1.35	1.52	1.62
Industrial, mining and transportation	73.33	34.11	80	85	90	95	1	1.06	1.13	1.19
Desert	30	13.95	72	82	83	87	1	1.14	1.15	1.21
Glacier/Permanent snow cover	0	0	0	0	0	0	0	0	0	0
Bare soil	25	11.63	72	82	83	87	1	1.14	1.15	1.21

Note: This study assumes the runoff coefficients of 0 for wetlands, lakes, rivers, and glaciers/permanent snow.

).

The runoff coefficient (R_c) data utilized in this study is derived from the aggregated data on runoff coefficients corresponding to various ecosystem types within the Yellow River Basin, as documented by Zhai et al. [26]. The study employed meta-analysis to consolidate the R_c data applicable to the Yellow River Basin, drawing from comparable regions or natural geographical environments. Furthermore, this paper categorizes the 12 common soil textures in the Yellow River Basin into four groups denoted as A, B, C, and D (

Table 3. Hydrological soil group (HSG) classification according to soil properties.

HSG	USDA Soil Texture Class	Category	Proportion of YRB Surface (%)
A	1, 2, 3	Sand, loamy sand, or sandy loam	30.90
B	4, 5, 6	Loam, silt loam, or silt	62.82
C	7	Sandy clay loam	2.46
D	8, 9, 10, 11, 12	Clay loam, silty clay loam, sandy clay, silty clay, or clay	3.50

), based on the research conducted by Zeng et al. [28].

2.3.2. Evaluation of Water Conservation Service Function (WCF)

This study conducts a quantitative analysis of the WCF using the water balance principle. The water balance principle evaluates the overall WCF by comparing the input and consumption of water in the study area. This widely used calculation method serves to assess and quantify the overall state of water conservation in the study area, providing insight into the regional water conservation condition [29]. The specific calculation formula is as follows:

$$WCF=P-E-R$$

(7)

In this formula, WCF refers to the amount of WCF (mm/a), P is the annual precipitation (mm/a), E is the annual evapotranspiration (mm/a), and R is the runoff depth (mm/a). The detailed calculation method of R is detailed in Section 2.3.1.

The precipitation data was acquired from the 1-km monthly precipitation dataset [30] for China, which is available from the TPDC (National Tibetan Plateau Data Center). This dataset has been downscaled using the Delta spatial downscaling method for the China region, based on the global 0.5° climate dataset by CRU and the high-resolution global climate dataset by World Clim. Furthermore, validation of the dataset was conducted using 496 independent meteorological observation points, and the results of the validation are deemed reliable.

2.3.3. Construction of the Space–Time Cube

The space–time cube is a form of multidimensional spatiotemporal data that integrates time and space [31], enabling visualization of spatiotemporal data, time-series forecasting, and analysis of spatiotemporal patterns [22]. In this study, the space–time cube takes the form of a netCDF file with three dimensions: x, y, and z. The x and y dimensions consist of grid data for the WCF, while the z dimension, representing the time dimension, is aligned with the research time step. Figure 4 illustrates the construction principle of the space–time cube based on ecosystem and soil type data.

Initially, remote-sensing data for ecosystem and soil types, alongside additional meteorological data for the study area, were inputted to acquire yearly WCF grid data for the period from 2012 to 2022 (Section 2.3.1/Section 2.3.2). The two-dimensional WCF grid data were then stacked in chronological order and augmented with temporal dimension properties to form a multidimensional grid layer. Subsequently, the neighborhood distance and time interval for the space–time cube were established in accordance with the experimental requirements to align with the multidimensional grid layer. Finally, these two components were combined to construct the necessary space–time cube, with the z-axis direction representing the progression of the time scale from bottom to top.

2.3.4. Emerging Hot Spot Analysis (EHSA)

EHSA is a spatial analysis tool that combines spatial and temporal information using a space–time cube. It incorporates the temporal dimension to identify the spatiotemporal heterogeneity of WCF more accurately and comprehensively. In this study, EHSA is applied to examine the statistically significant clustering trend of WCF in the study area from 2012 to 2022.

The EHSA process involves several steps. Initially, a space–time cube is created based on the preliminary steps (Section 2.3.3). This space–time cube, representing the WCF data, is then used as input for EHSA operations. The analysis requires setting parameters such as neighborhood distance, time step, and spatial relation conceptualization type. These parameters define the spatiotemporal domain (local range) considered in the analysis. The spatial relation conceptualization type determines the spatial scale of the bins used in the computation and offers four specific options: Fix Distance, K Nearest Neighbors, Contiguity Edges Only, and Contiguity Edges Corners (Figure 5). After the above, the Getis-Ord Gi* statistic values are calculated for each bin in the space–time cube. These values indicate the intensity of clustering or dispersion of the WCF data. Finally, the Mann–Kendall trend test is applied to analyze the trend of hot or cold spots. Based on the EHSA analysis, the spatiotemporal trend of WCF is classified into 17 patterns, including new, persistent, sporadic, intensifying, consecutive, diminishing, oscillating, historical, or no pattern detected.

2.3.5. Perturbation Analysis Method

To investigate the key parameters in the ESM-WBM model, this study conducted a quantitative evaluation of sensitivity indexes (SI) for five critical input parameters using two global sensitivity analysis methods: Sobol’ and EFAST (Extended Fourier Amplitude Sensitivity Test). Unlike traditional local sensitivity analysis methods, global sensitivity analysis methods can simultaneously evaluate the influence of multiple parameters on the model output, including the interaction effects between parameters. Prior research has indicated that, given the diverse characteristics and intricate usage scenarios of ecological and hydrological models, the Sobol’ and EFAST methods offer the most reliable and stable global sensitivity analysis [32]. The Sobol’ method, a classic algorithm, facilitates the quantitative determination of variance-based global sensitivity, enabling accurate quantification of both overall sensitivity and primary sensitivity of model parameters. On the other hand, the EFAST method combines the strengths of the Sobol’ algorithm in assessing parameter interactions with the efficient sampling of the FAST algorithm. To ensure the robustness of the results, this study employed both methods to calculate SI, thereby mitigating potential systemic errors associated with a single method.

3. Results

3.1. The Characteristics of Changes in Water Conservation Service Function (WCF)

3.1.1. Spatial Patterns of Multi-Year Average WCF

Figure 6 presents the spatial distribution of the WCF in the study area. In 2012, the disparity between the maximum and minimum values of regional WCF stood at 580.01 mm, whereas in 2022, the difference decreased to 494.55 mm. The spatial distribution patterns of average annual precipitation (2012–2022) and actual evapotranspiration (AET) demonstrate a striking resemblance, both exhibiting a gradual increase from the northwest to the southeast, characterized by a distinct stepped pattern (Figure 6a,b). The southeastern part of the RBW basin and the HLS basin encompass the majority of areas with high rainfall, where agricultural activities are relatively abundant, resulting in increased runoff compared to regions with lower precipitation levels. Observing Figure 6c,d, it becomes apparent that high WCF value areas (>300 mm) have shifted toward the southeast since 2012, accompanied by a notable decline in the RAW and HLN basins. Likewise, the northern and central areas of the RBW basin also experienced a significant decrease in WCF.

This study utilized the Natural Breaks Classification Method (JENKS) to categorize the WCF in the YRB study area into five levels of importance. The regional distribution and area proportions of WCF importance levels were visually represented in a spatial visualization format (Figure 7). The analysis revealed that, across the entire study area, the highest proportions were observed at levels V and III, accounting for 28.24% and 25.95%, respectively, reflecting a distinct stepped distribution.

For a more detailed assessment of the WCF distribution in each sub-basin, Figure 7 also illustrates the detailed classification of WCF importance levels within each sub-watershed. The upstream area of the study region primarily comprises levels Ⅲ and Ⅳ WCF areas, with level Ⅲ accounting for 45.17% and 33.25% in the RAW and HLN basins, respectively. What’s more, in the HLN basin, the proportions of level I and II areas are notably higher at 15.32% and 27.78%, respectively, which is attributed to the distribution of settlements and industrial land in the northern part of HLN. As the area progresses downstream, the proportion of level V areas steadily increases, reaching 31.78% in the RBW basin and peaking at 69.01% in the HLS basin, influenced by the widespread distribution of broad-leaved forests and shrubs in the southeastern part.

In summary, the substantial proportion of the study area in the downstream basins is classified as level V, which necessitates prioritized protection. Conversely, the low-level regions in the northern part of the HLN basin, the central and northwestern parts of the RAW basin, and the northwestern part of the RBW basin experience more strained relationships between water supply and demand, thus requiring specific attention to water resource utilization and allocation issues.

3.1.2. Inter-Annual Variation of WCF

Figure 8 presents the inter-annual variation of average WCF trends and its four sub-basins from 2012 to 2022. The annual average WCF in the study area exhibits a declining trend, with an average annual growth rate of −0.9377 mm/a. It reached a peak of 403.25 mm in 2020 and a trough of 272.94 mm in 2015. In contrast to the inter-annual variation trends in runoff depth, as presented in Section 3.1.1 (2), the WCF in the HLS and RBW basins displays an increasing trend, with growth rates of 1.691 and 1.096 mm/a, respectively. However, the WCF in the HLN and RAW basins demonstrates a declining trend, with rates of −4.8499 and −4.8788 mm/a, respectively, based on the trend lines from linear regression. The absolute values of their slopes are notably higher than the absolute values of the increasing slopes in the other two basins (HLS and RBW). Notably, the inter-annual variation in WCF in each sub-basin is not substantial, with frequent occurrences of elevated values observed only in 2020.

The use of annual change slope analysis enabled the examination of the temporal trend of WCF from 2012 to 2022 (Figure 9). In the downstream of RBW and HLS basins, the annual change slopes surpassed 3.10 mm/year, signifying a substantial increase in WCF in these regions over the past decade. Furthermore, specific sections in the southern areas of these two basins exhibited annual change slopes of WCF exceeding 6 mm/a, which can be attributed to ecological compensation measures such as land reforestation and eco-relocation initiatives implemented in the Yellow River Basin. Despite the ascending trend in these areas, the majority of the region experienced a declining trend. Particularly, the RAW basin and the majority of the HLN basin demonstrated a decreasing trend, with change slopes plummeting to −4.85 mm/a in the northern part of the RAW basin and approximately −6 mm/a in most regions of the HLN basin.

3.2. Parameter Sensitivity Analysis

The EFAST and Sobol’ results of the global sensitivity analysis for the five key parameters of the ESM-WBM model are showcased in Figure 10. Figure 10a,b illustrate the Major Sensitivity Index (MSI) and Total Sensitivity Index (TSI), respectively. The MSI measures the impact of each input variable on the variance of the dependent variable and is utilized to gauge the sensitivity of individual variables within the model. Conversely, the TSI reflects the variance of each input variable and its interaction with other input variables, indicating the degree of interaction between parameters.

The similarities between the sensitivity indices obtained from the two methods suggest the reliability of the analysis results. The most sensitive parameters are R_ist, ΔR, and P. Overall, excluding the TSI corresponding to R_ist, the SI calculated by the Sobol’ method slightly exceeds that of the EFAST method. Specifically, both methods yield high values for the MSI and TSI of R_ist, at 0.471/0.472 and 0.541/0.540, respectively. This indicates the model’s emphasis on considering the ecosystem and soil types of the study area in quantitatively analyzing WCF. Precipitation P and runoff change ΔR are also highly sensitive parameters. The TSI of the runoff change ΔR (TSI = 0.274) is notably higher than its first-order sensitivity (MSI = 0.205), suggesting a relatively strong interaction with other parameters.

3.3. Spatiotemporal Heterogeneity of Water Conservation Service Function (WCF)

3.3.1. Evaluation Results of Emerging Hot Spot Analysis (EHSA)

To delineate the spatiotemporal heterogeneity of the WCF, we employed the emerging spatiotemporal hot spot analysis. In this analysis, clusters with high WCF values within the designated spatiotemporal neighborhood were identified as spatiotemporal hot spots, while areas with clustered low WCF values were classified as spatiotemporal cold spots. Figure 11 delineates (a) the spatial distribution and (b) statistical analysis results for each pattern, respectively. As depicted in Figure 11, a total of 16 EHSA patterns were identified in this investigation, except the No Pattern Detected area. In addition, to enhance the presentation of the EHSA results,

Table 4. Definition of 17 patterns in EHSA.

Legend	Pattern Name	Definition
	Undetected Patterns	Does not fall into any of the hot or cold spot patterns defined below.
	New Hot Spot	A location that is a statistically significant hot spot for the final time step and has never been a statistically significant hot spot before.
	Consecutive Hot Spot	A location with a single uninterrupted run of at least two statistically significant hot spot bins in the final time-step intervals. The location has never been a statistically significant hot spot prior to the final hot spot run, and less than 90 percent of all bins are statistically significant hot spots.
	Intensifying Hot Spot	A location that has been a statistically significant hot spot for 90 percent of the time-step intervals, including the final time step. In addition, the intensity of clustering of high counts in each time step is increasing overall, and that increase is statistically significant.
	Persistent Hot Spot	A location that has been a statistically significant hot spot for 90 percent of the time-step intervals with no discernible trend in the intensity of clustering over time.
	Diminishing Hot Spot	A location that has been a statistically significant hot spot for 90 percent of the time-step intervals, including the final time step. In addition, the intensity of clustering in each time step is decreasing overall, and that decrease is statistically significant.
	Sporadic Hot Spot	A statistically significant hot spot at the last time step interval has also been a recurring hot spot throughout history. Up to 90% of time-step intervals are already statistically significant hot spots, and none of the time-step intervals are statistically significant cold spots.
	Oscillating Hot Spot	A statistically significant hot spot at the last time step interval that has a history of a statistically significant cold spot in the previous time step. Up to 90% time-step intervals are already hot spots of statistical significance.
	Historical Hot Spot	The most recent period is not a hot spot, but at least 90% of the time step intervals are already statistically significant hot spots.
	New Cold Spot	This location is a statistically significant cold spot for the last time step and has never been a statistically significant cold spot before.
	Consecutive Cold Spot	This position has a single uninterrupted run of a cold spot bar with two statistical significances in the last time step interval. It is never a statistically significant cold spot until the final cold spot runs, and up to 90% of all bars are statistically significant cold spots.
	Intensifying Cold Spot	This position is already a statistically significant cold spot for 90% of time-step intervals, including the last time step. In addition, the smaller number of clusters in each time step increased in strength overall, and the increase was statistically significant.
	Persistent Cold Spot	This location is already a statistically significant cold spot with a 90% time-step interval, and there is no clear trend indicating that the clustering strength has changed over time.
	Diminishing Cold Spot	This position is already a statistically significant cold spot for 90% of time-step intervals, including the last time step. In addition, the smaller number of clusters in each time step was reduced overall, and the reduction was statistically significant.
	Sporadic Cold Spot	A statistically significant cold spot for the final time-step interval with a history of also being an on-again and off-again cold spot. Less than 90 percent of the time-step intervals have been statistically significant cold spots, and none of the time-step intervals have been statistically significant hot spots.
	Oscillating Cold Spot	A statistically significant cold spot for the final time-step interval that has a history of also being a statistically significant hot spot during a prior time step. Less than 90 percent of the time-step intervals have been statistically significant cold spots.
	Historical Cold Spot	The most recent time period is not cold, but at least 90 percent of the time-step intervals have been statistically significant cold spots.

Tips: The main content of this table is derived from the How Emerging Hot Spot Analysis Works section of the ArcGIS Pro official website.

provides definitions for the 17 patterns mentioned above. The distribution of hot spots and cold spots appeared relatively independent. Specifically, hot spots were predominantly concentrated in the southeast, while cold spots were mainly distributed in the northwest. Furthermore, a small-scale hot spot cluster was observed in the southern part of the HLN basin.

Figure 11b provides a detailed representation of the proportion and spatial distribution of different types, facilitating an in-depth analysis of the spatial heterogeneity of WCF. In terms of hot spot types, the study area primarily exhibits three patterns: Intensifying Hot Spot (4.1178%), Persistent Hot Spot (6.5228%), and Sporadic Hot Spot (11.2880%). The Intensifying Hot Spot is predominantly observed in the southern parts of the RBW and HLS basins, indicating statistically significant increasing trends in WCF at 90%-time step intervals, suggesting a significant annual growth trend in WCF in these areas. The central parts of the HLS and the southern parts of the RBW basin exhibit Persistent Hot Spots, signifying sustained high WCF levels from 2012 to 2022 without significant temporal changes. The Sporadic Hot Spot is widely dispersed in the southern central region of the study area, encompassing the downstream of RBW and the entire HLS basin, which indicates intermittent dispersion across various areas of the region. In addition, in these areas, at most 90% of the time-step intervals are statistically significant hot spots, highlighting strong spatial heterogeneity in the distribution of WCF without significant spatial clustering and temporal fluctuations.

The areas with low WCF in the study area are primarily Persistent Cold Spot (11.5%), Sporadic Cold Spot (6.9%), and Oscillating Cold Spot (19.9%). Persistent Cold Spot is mainly distributed in the northern regions of HLN and RBW, with a few scattered in the areas of RAW, indicating persistent low WCF, possibly related to urban and industrial land use in the region. In contrast to the concentrated distribution of Persistent Cold Spot, Sporadic Cold Spot is dispersed in the northern regions of HLN and RAW, indicating that at most 90% of the time-step intervals are statistically significant cold spots. Oscillating Cold Spot covers the northern part of the RAW basin, the northeastern part of the RBW basin, and the southern part of the HLN basin, manifesting a concentrated and clustered distribution in this study area. This indicates that these areas are statistically significant cold spots in the last time step interval and have a history of statistically significant hot spots in previous time steps, implying unstable fluctuations in WCF in this area, with previously high values decreasing to low values recently. It is worth noting that a New Cold Spot appears in the central part of RAW, representing a statistically significant cold spot in the last time step, indicating a recent decrease in WCF in this area, necessitating targeted mitigation.

3.3.2. Impact of Neighborhood Distance/the Type of Spatial Relationship Conceptualization on Expression of WCF Spatial Heterogeneity

(1)

Influence of neighborhood distance on the expression of WCF spatial heterogeneity.

In this study, we defined neighborhood distance (ND) at intervals of 100 m, 150 m, 200 m, 250 m, 500 m, 750 m, and 1000 m. Subsequently, we calculated the Getis-Ord Gi* statistics to analyze the EHSA spatial distribution (Figure 12) and investigate the impact of neighborhood distance on WCF spatial heterogeneity. With the increase in ND from 100 m to 150 m, a pronounced distribution of hot spot types was observed in the southern part of the study area, accompanied by notable changes in the cold spot types in the northern region. When ND was increased to 200 m, the pattern of hot spot distribution in the southern part displayed even greater heterogeneity. However, with a further increase in ND, there was a significant expansion in the area of “Sporadic Cold Spot” and “Oscillating Cold Spot” in the northern part. The same situation occurred in the southern part, where “Intensifying Hot Spot” and “Persistent Hot Spot” gradually occupied the entire downstream area of the HLS and RBW basins. Nevertheless, once ND exceeded 750 m, its impact on the heterogeneity expression became almost negligible. This unfavorably affects the portrayal of WCF spatial heterogeneity expression, suggesting that it weakens as ND increases and exhibits a declining correlation trend.

To analyze the transition of patterns at various neighborhood distance gradients, transition matrices were introduced to assess the relevant characteristic points at intervals of 100 m, 150 m, 200 m and 250 m (Figure 12). With the increase in ND from 100 m to 150 m, the areas of “No Pattern Detected”, “Sporadic Cold Spot”, and “Sporadic Hot Spot” decreased by 15.270%, 6.778%, and 2.862%, respectively, resulting in total area changes of −18512.22 km², −8217.64 km², and −3470.12 km². Concurrently, “Oscillating Cold Spot” and “Persistent Cold Spot” increased by 11.217% and 4.253%, respectively, equating to total area changes of 13598.62 km² and 5156.43 km². Notably, “No Pattern Detected”, “Sporadic Cold Spot”, and “Sporadic Hot Spot” underwent the most significant changes. “No Pattern Detected” mainly transformed into “Sporadic Hot Spot” (4.57%), “Sporadic Cold Spot” (2.99%), “Oscillating Hot Spot” (1.28%), and “Oscillating Cold Spot” (5.47%). The main transformation direction for “Sporadic Cold Spot” was “Oscillating Cold Spot” (4.85%), “Persistent Cold Spot” (3.78%), and “Sporadic Cold Spot” (3.90%). The primary transformation directions for “Sporadic Hot Spot” were “Intensifying Hot Spot” (1.96%) and “Persistent Hot Spot” (5.37%). In addition, the transition patterns from 150 m to 200 m and from 200 m to 250 m gradients exhibited consistency.

(2)

Influence of conceptualizing neighborhood relationships on WCF spatial heterogeneity expression.

Four different conceptual types of spatial relationships were utilized to examine the spatial heterogeneity of WCF in the study area (Figure 12a). Results from three of the methods are largely consistent, with the exception of the “Contiguity Edges Only” type, which demonstrates notably greater spatial disparities. Furthermore, the results obtained from the “Contiguity Edges Only” method closely approximate those from the “Fixed Distance” method (ND = 100 m). It is worth noting that the “Fixed Distance” method applied a neighborhood distance of 150 m in Figure 12a.

Figure 13b illustrates the proportions (≥1%) of the main patterns obtained from the four methods. This aligns with the findings in Figure 13, wherein the results from the “Fixed Distance” (ND = 150 m), “K Nearest Neighbors”, and “Contiguity Edge Corners” methods display high similarity, while the “Contiguity Edge Only” method demonstrates larger errors. Specifically, the areas of the “Diminishing Hot Spot”, “Consecutive Hot Spot”, “Intensifying Hot Spot”, “Sporadic Hot Spot”, “Persistent Hot Spot”, and “Persistent Cold Spot” obtained from the “Contiguity Edge Only” method are underestimated by 0.651%, 2.550%, 11.217%, 1.542%, 4.253%, and 5.429%, respectively. In contrast, “Historical Hot Spot”, “Oscillating Hot Spot”, “Oscillating Cold Spot”, and “No Pattern Detected” are overestimated by 0.426%, 15.270%, 6.778%, and 2.862%, respectively.

4. Discussion

4.1. The Characteristics of Changes in WCF

The study concentrates on a sub-basin within the Yellow River Basin to quantitatively analyze the WCF. This analysis is conducted utilizing the ESM-WBM model, employing actual observed data obtained from hydrological stations. The results show that, in the study area, the WCF displays a spatial pattern with lower values in the northwest and higher values in the southeast, with a slight overall decrease from 2012 to 2022. It is worth noting that the study area is located in the Loess Plateau region, which has influenced the overall spatial pattern of WCF distribution. In the HLN basin, located on the eastern side of the study area, the higher altitude and fragile ecosystem structure, combined with the cold climate, result in a weaker WCF due to the high loss of surface water. On the other hand, the downstream region of the study area, with lower altitude, more abundant precipitation, and favorable meteorological conditions, supports diverse vegetation types due to favorable thermal conditions and gentle terrain, which aids in water resource conservation [33]. Areas with higher WCF values, such as the southern parts of the RAW and RBW basins, as well as almost the entire HLS basin, are characterized by abundant rainfall and relatively high vegetation coverage. These areas consist of vast forests, grasslands, broad-leaved shrubs, and marshy grasslands, which provide strong capabilities to intercept and retain precipitation [34]. Conversely, regions with lower WCF values, such as the RAW and HLN basins, have inferior precipitation and vegetation coverage, along with shallow soil depth, making them prone to surface runoff rather than retention, consistent with the study’s findings on runoff depth. Furthermore, these areas have a high proportion of human settlements and industrial land, contributing to the weaker WCF in the region.

4.2. Spatiotemporal Heterogeneity of WCF

EHSA analyzes the WCF space–time cube and categorizes it into 17 patterns, which in turn allows for adjustments to ecological compensation policies in related areas based on each pattern. In the study area, hot spots are mainly concentrated in the downstream region, while cold spots are primarily situated in the northern and northwestern parts of the upper reaches. Specifically, the “Intensifying Hot Spot”, located at the southern end of the HLS and RBW basin in the downstream area, signifies a significant increase from 2012 to 2022. It indicates a substantial improvement in the ecological environment and highlights the effectiveness of existing restoration or compensation measures. The “Persistent Hot Spot” in the middle of HLS and RBW indicates that this area has maintained a desirable WCF since 2012 and has not significantly changed over time, thereby requiring the continuation of existing protection policies without additional compensation or restoration. The distribution of cold spots in the study area primarily includes “Persistent Cold Spots”, “Oscillating Cold Spots”, and “Sporadic Cold Spots”. The distribution of “Persistent Cold Spots” is mainly in the northern part of HLN and the northwestern part of RBW, signifying that significant ecological compensation or governance in this area may present challenges, including high input costs and slow effects. “Oscillating Cold Spots” are widespread in the northern part of RAW and the middle and southern parts of HLN, indicating fluctuations in the WCF in these areas and emphasizing the importance of ecological protection and restoration in these areas. Additionally, a “New Cold Spot” has emerged in the middle of RAW, representing a statistically significant cold spot at the last time step, necessitating specific mitigation measures due to recent damage to the ecosystem service function in this area.

4.3. EHSA Influencing Factors Analysis

The influence of neighborhood distance and the conceptualization type of spatial relationships on the manifestation of spatial heterogeneity in WCF analysis is explored in this study. The outcomes reveal noteworthy disparities among the various spatial relationship conceptualization types, with Contiguity Edge Only exhibiting significant aberration. This finding suggests that this particular spatial relationship model inadequately represents the spatial dispersion of WCF within the study area. Conversely, Fixed Distance, K Nearest Neighbors, and Contiguity Edge Corners demonstrate effectiveness in capturing the spatial heterogeneity of WCF distribution. The study establishes seven gradients (100 m/150 m/200 m/250 m/500 m/750 m/1000 m) based on the space–time cube characteristics to investigate neighborhood distances. Analysis indicates that significant variations in WCF spatial heterogeneity expression occur at neighborhood distances of 100/150/200/250 m, with a discernible escalating trend as the distance increases. Conversely, neighborhood distances exceeding 500 m categorize numerous areas as exhibiting a similar pattern, hindering the accurate representation of WCF spatial heterogeneity in the region. The impact of neighborhood distance on WCF spatial heterogeneity manifestation is thus deemed equivocal.

Moreover, considering the aforementioned findings, it is hypothesized that a threshold may exist for the influence of neighborhood distance on the spatial heterogeneity expression of WCF. Adjusting the neighborhood distance beyond this threshold is anticipated to substantially influence the spatial heterogeneity expression outcomes. To pinpoint this threshold precisely, a recommendation is made to employ a Monte Carlo model by conducting an extensive sample range of neighborhood distances following a normal distribution pattern. Subsequently, statistical analysis can be performed on the corresponding Getis-Ord Gi* values, and the EHSA grid results can be summarized and subjected to linear regression fitting. Unfortunately, due to current limitations in computational resources and algorithm efficiency, the practical implementation of this methodology is currently unfeasible in the context of this study.

4.4. Limitations and Uncertainties

The limitations and uncertainties of this study primarily stem from three main aspects:

(1)

First, one main source of uncertainty is the meteorological data, which is a crucial input. While ESM uses high-resolution ecosystem type data to accurately describe the spatial patterns of WCF, the precipitation and evapotranspiration data used still have uncertainties. On the one hand, the precipitation data is downscaled in China based on the global 0.5° climate dataset published by CRU and the global high-resolution climate dataset published by World Clim using the Delta spatial downscaling scheme. This data was validated using actual data from 496 meteorological stations across China, but it still inevitably has a certain systematic error at the sub-basin scale. On the other hand, the evapotranspiration data used in this study comes from MOD16A2, which ignores the evapotranspiration of water bodies such as rivers and lakes, hence also introducing uncertainties.

(2)

Limitations arise from the research design of this method. ESM neglects the influence of terrain slope, which can lead to errors in the runoff capacity of ecosystems at smaller watershed or landscape scales. Similarly, the presence of cropland and residential areas in the northern part of the study area may underestimate water consumption and thus lead to a higher WCF value. The next step in this research is to consider terrain slope more comprehensively and scientifically consider the supply-and-demand relationship of certain land use types.

(3)

This study provides a preliminary discussion of the factors affecting the EHSA analysis results based on the WCF space–time cube, but specific mechanism explanations still require further research. For example, in exploring the neighborhood distance’s expression of WCF spatial heterogeneity, a Monte Carlo model can be introduced, using large-scale sampling of neighborhood distances, such as through a normal distribution and regression analysis, finally converging to obtain the best-fit parameter values. However, this study is limited by factors such as software underlying code and algorithm efficiency and only carried out related analyses using seven representative gradients. Hence, how to improve the efficiency of software or algorithms, such as using parallel computing, to achieve large-scale simulation is a possible direction for future researchers.

5. Conclusions

This study aims to quantitatively analyze the WCF of ecosystems and explore their spatial heterogeneity. To achieve this, the study employed an evaluation method for WCF based on ecosystem and soil types and conducted a sensitivity analysis of this evaluation method’s parameters. Based on this, the study constructed a space–time cube using remote sensing data of ecosystem and soil types and used EHSA to conduct data mining, exploring the impact of neighborhood distance and conceptual spatial relationship types. The application of EHSA can detect the WCF status and its spatiotemporal patterns at specific locations, providing scientific support for differentiated management measures for water resource allocation and ecological protection in related areas. The main conclusions are as follows.

(1)

The annual variation of WCF in the study area is significant. During the period from 2012 to 2022, a decreasing trend was observed in the WCF of the study area basin, with rates of −0.9377 mm/year. HLN and RAW showed declines at rates of −4.8499 mm/year and −4.8788 mm/year, respectively, necessitating corresponding mitigation measures. On the other hand, HLS and RBW exhibited increases in WCF at rates of 1.691 mm/year and 1.096 mm/year, respectively, indicating significant ecological improvement in these areas.

(2)

The ESM-WBM comprehensively depicts changes in ecosystem patterns, providing a more accurate representation of the spatial distribution of WCF. Global sensitivity analysis using EFAST and Sobol’ methods shows that the output of ESM is highly sensitive to R_ist, indicating its ability to capture ecosystem and soil pattern changes effectively. Additionally, the results are sensitive to runoff and precipitation variations, demonstrating the method’s capacity to integrate observed runoff data from hydrological stations. Moreover, the ESM method, relying on high-precision remote sensing data of ecosystem types, inherits the accuracy of this data, enabling a detailed description of the spatial pattern of WCF.

(3)

The Emerging Hot Spot Analysis (EHSA) based on the space–time cube allows for accurate and comprehensive identification of the spatial heterogeneity of WCF. By visually representing the spatiotemporal aggregation patterns of WCF, EHSA can more effectively assist in the development of differentiated ecological protection management policies. EHSA analyzes the WCF space–time cube and categorizes it into 17 patterns, which in turn allows for adjustments to ecological compensation policies in related areas based on each pattern. Based on the results of the analysis, the “Intensifying Hot Spot” in the research area shows significant improvement in the ecological environment due to the effectiveness of current restoration and compensation measures. Continuing current conservation efforts without additional restoration investment are recommended for the “Persistent Hot Spot”. However, addressing challenges related to high costs and slow effectiveness for ecological management is crucial in the “Persistent Cold Spot” area. The “Oscillating Cold Spot” should be a priority for ecological conservation and restoration efforts, while the “New Cold Spot” calls for immediate mitigation measures to counter recent damage to ecosystem services.

(4)

The conceptualization of neighborhood distance and spatial relationships significantly impacts the results of EHSA. In this study, seven gradients of neighborhood distance (ND) were defined, and EHSA was conducted separately for each. The results revealed a dual characteristic of ND’s impact on the results, with the existence of a threshold. Specifically, when ND is less than 500 m, increasing ND benefits the expression of the spatial heterogeneity of WCF, while surpassing 500 m has an adverse effect, and once ND exceeds 750 m, its impact diminishes. Except for “Contiguity Edge Only”, the results of “K Nearest Neighbors” and “Contiguity Edge Corners” were similar in their conceptualized types of spatial relationships, with significant deviations. Additionally, these results were highly consistent with “Fixed Distance (ND = 150 m)”, indicating similar performances in describing the inherent spatial relationships of the space–time cube of WCF. This provides a reference for exploring key parameters best suited for describing the spatial relationships of the WCF spatiotemporal cube in future research. This study offers references for analyzing the spatial heterogeneity of WCF, providing a theoretical foundation for regional water resource management and ecological restoration policies with tailored strategies.