Analysis of Spatial and Temporal Characteristics and Spatial Flow Process of Soil Conservation Service in Jinghe Basin of China

Abstract

The supply and demand of ecosystem services and their mutual feedback are important for the formulation of basin ecological environmental policies. Simulation of the spatial flow of ecosystem services can clarify the division of areas and can support policy development. This paper takes the Jinghe Basin in the Loess Plateau of China as the case-study area to simulate the spatial flow of soil conservation service on different scales. The results showed that (1) soil erosion situations in Jinghe Basin improved overall, with a trend of first aggravating and then recovering between 2005 and 2015; (2) the amount of annual soil conservation in the basin accounted for more than 75% of the potential soil erosion and showed a trend of first increasing and then decreasing; and (3) using digital elevation model (DEM) data and ArcGIS software, the experiment divided the basin into sub-basins (58 in total) and hydrological response units (HRUs) (e.g., 2181 HRUs in sub-basin #1), which were used to quantify the spatial flow direction and the corresponding amount of soil conservation service on the “HRU—river-sub-basin” scale. The divided supply and demand helped quantify the spatial flow pattern of soil conservation services from HRU to the sub-basin.

1. Introduction

Ecosystem services [1] are the environmental conditions and utilities formed by ecosystems and ecological processes critical for human survival and development. They are closely related to human life [2,3]. Ecosystem services [4] are generally divided into supply services, regulation services, culture services, and support services. Soil conservation service refers to the erosion control ability of the ecosystem to prevent soil loss and to store sediment and belongs to the scope of regulation services. The assessment of soil conservation services is often based on the revised universal soil loss equation (RUSLE) [5,6,7,8,9]. The RUSLE model is an empirical slope soil erosion prediction model formed by a large number of soil erosion experiments and quantitative observation data, and on this basis, many scholars have applied this model in a lot of research on soil erosion [10,11,12]. Its reasonable structure and strong representative parameters have been widely studied and applied in China [13].

An ecosystem service is used by humans, passed from one area to another so that an ecosystem service flow can be formed. Ecosystem service flow is a vital link connecting the service supply and demand of ecosystem services. Initial research on ecosystem service flow was research on spatial characteristics of the supply services [14]. Later, with an in-depth study of the relationship between the supply and demand of ecosystem services, ecosystem service flow gradually became a research hotspot. For example, based on the Bayesian network model (BBN), Tang et al. constructed an ecosystem service network system in the Jinghe Basin in 2015 and provided spatial optimization suggestions for ecosystem services such as water production, the net primary productivity (NPP), soil erosion, and agricultural productivity in the basin [15]. Li et al. [16] quantified the service flow of freshwater ecosystems in the Beijing–Tianjin–Hebei region of China based on the integrated valuation of ecosystem services and tradeoffs (InVest) model to assess the regional freshwater security pattern. Bagstad et al. [17] used artificial intelligence methods to quantify ecosystem services and simulated the supply areas and demand areas. Xu et al. [18] determined the demand of water supply services (WSSs) from the perspective of ecosystem service flow and simulated the flow path of WSSs. Li et al. [19] simulated spatial matching of supply and demand for carbon sequestration services and flow ratios. Based on the relationship between adjacent grids, they matched carbon sequestration points with carbon emission points to illustrate the ecological process of carbon sequestration services. Most of the articles above used quantitative assessment of the supply and demand relationship of ecosystem services to divide the supply and demand, and to simulate the flow of ecosystem services. There are few articles that consider topographic factors into the spatial flow of ecosystem services [18,20,21]. In addition, there are seldom studies on the specific flow paths of ecosystem services.

This paper holds ecosystem service flow as a vector, and its flow path needs to be represented. The path algorithm for basin network analysis and water catchment area calculation is a type of geographic information system (GIS) digital terrain analysis. The path algorithm describes the high-to-low transmission and flow of surface materials such as water and sediment between topographic units, a critical technical link for automatic basin network extraction. There are many path algorithms currently used for basin network analysis and catchment area calculation. For example, the deterministic eight-direction algorithm (D8) [22], the random eight-direction algorithm (Rho8) [23], the Freeman multiple flow direction (FMFD) algorithm [24], The digital elevation model networks (DEMON) algorithm [25], and the Dinf algorithm [26]. Liu et al. [27] compared the path algorithms above, and the result showed that the single flow direction algorithm is greatly affected by grid structure; that the DEMON algorithm can give an accurate flow direction and catchment area calculation results, but the algorithm design is too complicated and the robustness is low; and that the Dinf algorithm based on these aspects can provide an exact specific catchment area (SCA) value and an ideal SCA distribution. The study in [28] showed that, in a two-dimensional space, surface runoff and the resulting soil loss were determined by the runoff area per unit contour length (referred to as unit catchment area), which is the slope length factor.

This paper focuses on the study of the spatial flow path of soil conservation services in 2005–2015, hoping that the experimental results will provide a certain scientific basis for the formulation of future ecological environmental policies in the basin.

2. Research Methods

2.1. Overview of the Study Area

The Jinghe River is the first tributary of the Weihe River and the second tributary of the Yellow River [15]. The Jinghe Basin is located at the center of the Loess Plateau in China (Figure 1). The basin covers 4.5421 × 10¹⁰ m² area, spanning 31 counties (cities) in three provinces (regions) of Shaanxi, Gansu, and Ningxia in China. The terrain of the Jinghe Basin is higher in the northwest and lower in the southeast, with the highest elevation of 2908 m [29]. The Jinghe Basin is one of the most severely eroded areas on the Loess Plateau [30]; 4.3% of the Jinghe Basin area is mountainous, 41.7% is loess tableland and broken plateau areas, and 48.8% is loess hilly and gully regions [30]. The basin has deep layers (50–80 m) of loess, composed mainly of fine sand, silt, and clay, with silt accounting for up to 50%; has high porosity; and is prone to landslides [30]. The basin is in a transitional area from semi-humid to semi-arid [31]. The annual average temperature of the Jinghe Basin is 8 °C, and the average yearly precipitation is 350–650 mm [15]. The northern part of the basin is a cold-temperate semi-arid zone with an average annual rainfall of less than 400 mm, and the southern part is a cold-temperate semi-humid zone with an average annual rainfall of more than 500 mm [32]. The main vegetation types in the area are forests, shrubs, and grasslands [33]. Among them, forests account for about 28.7% of the total area of the basin, shrubs account for about 10.2% of the total area of the basin, and grassland accounts for about 32.1% of the total area of the basin [34].

2.2. Data Source and Processing

The research datasets include normalized difference vegetation index (NDVI) data, meteorological data, DEM data, soil data, land-use data, sediment data, and socioeconomic data. The data source and main processing processes are summarized as follows: NDVI data was downloaded from the National Aeronautics and Space Administration (NASA) website (https://www.nasa.gov/), using the MODIS Reprojection Tool (MRT) provided by NASA for batch processing of images. NDVI data were in grid format with a resolution of 250 m. In this paper, the rainfall rate data in the China meteorological forcing dataset [35] was selected as the basic rainfall data of the experiment. The data were downloaded from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/zh-hans/). After Python language processing, the daily rainfall data were obtained. The data were in grid form with a resolution of 1.0 × 10⁴ m. The DEM data were downloaded from the geospatial data cloud platform (http://www.gscloud.cn/). DEM data were in grid format with a resolution of 30 m. The soil data were downloaded from the Cold and Arid Regions Science Data Center (http://data.casnw.net/portal/). The soil data were in grid format with a resolution of 1.0 × 10³ m. The land-use data were downloaded from the resource and environmental data sharing platform of the Chinese Academy of Sciences (http://www.resdc.cn/). The land-use data were in grid format with a resolution of 1.0 × 10³ m. The sediment data were the measured annual sediment transport data of Zhangjiashan Hydrological Station in the Jinghe Basin of China, which were collected from the Ministry of Water Resources of the People’s Republic of China (http://www.mwr.gov.cn/sj/tjgb/zghlnsgb/). The socioeconomic data were collected from the relevant records of the National Bureau of Statistics (http://www.stats.gov.cn/). The 1:1 million topographic maps of China were from the Geographic Data Sharing Infrastructure, College of Urban and Environmental Science, Peking University (http://geodata.pku.eud.cn). In order to reduce the deviation caused by spatial data, this paper transformed all spatial data into a unified coordinate system and resampled them to 30 m based on the resampling tool of ArcGIS 10.2.

2.3. Research Methods

2.3.1. Soil Erosion Model

Based on the RUSLE model [14], the amount of potential soil erosion (A_p), the amount of actual soil erosion (A_r), and the amount of soil conservation (A_c) in the Jinghe Basin were calculated from 2005 to 2015. The experiment selected observation data from the Zhangjiashan Hydrological Station at Jinghe Basin in the River Sediment Bulletin of China of the Ministry of Water Resources of the People’s Republic of China to verify the actual soil erosion simulated by the model. $A_p$ refers to soil erosion in the ecosystem without vegetation cover and water and soil conservation measures, that is, C = 1, P = 1. $A_r$ refers to soil erosion in the ecosystem with vegetation cover, and water and soil conservation measures. $A_c$ is the difference between $A_p$ and $A_r$. They can be calculated as follows:

$$A_p=R\times K\times LS$$

(1)

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

(2)

$$A_c=A_p-A_r$$

(3)

where $A_p$ is the amount of potential soil erosion (kg/m²), $A_r$ is the amount of actual soil erosion (kg/m²), $A_c$ is the amount of soil conservation (kg/m²), $R$ is the rainfall erosivity factor (J*mm/m²/s), $K$ is the soil erodibility factor (kg*m²*s/m²/J/mm), $LS$ is the slope length and steepness factor (dimensionless), $C$ is the surface vegetation cover factor (dimensionless), and $P$ is the soil conservation measure factor (dimensionless).

The rainfall erosivity factor R refers to the potential erosivity caused by rainfall. It is the most important factor leading to soil loss. This experiment adopted the empirical formula proposed by Wischmeier et al. [36], which is more suitable for the Loess Plateau region [37], using monthly rainfall. The calculation formula is as follows:

$$R={\displaystyle \sum _{i=1}^{12}(1.73\times {10}^{\left(1.5\times \log \left(\frac{P_i{}^2}{P}\right)-0.8188\right)})}$$

(4)

where $P_i$ is the monthly rainfall (mm) and $P$ is the average annual rainfall (mm).

The soil erodibility factor K refers to the calculation method in the erosion-productivity impact calculator (EPIC) model by Williams et al. [38]. This method was applied in small basins by Cai et al. [39]. The results calculated by this formula are very close to the measured values, and the method is simple. The calculation formula is as follows:

$$K=0.1317\times \left\{0.2+0.3\exp \left[-0.0256\times SAN\times \left(1-\frac{SIL}{100}\right)\right]\right\}\times {\left(\frac{SIL}{SIL+CLA}\right)}^{0.3}\times \left[1-\frac{0.25C}{C+\exp \left(3.72-2.95C\right)}\right]\times \left[1-\frac{0.7SN}{SN+\exp \left(-5.51+22.9SN\right)}\right]$$

(5)

$$SN=1-\frac{SAN}{100}$$

(6)

where $SAN$ is the sand content (%), $SIL$ is silt content (%), $CLA$ is clay content (%), and $C$ is the organic carbon content (%).

The slope length and steepness factor were calculated by the method proposed by Wischmeier et al. [11]. Since the Jinghe Basin is located in the Loess Plateau region and the slope is steep, the slope factor calculation adopted the steep slope calculation method proposed by Mccool et al. [40] and Liu et al. [41]. The formulas are as follows:

$$L={\left(\lambda /22.13\right)}^{\alpha }$$

(7)

$$\alpha =\frac{\beta }{1+\beta }$$

(8)

$$\beta =\left(\sin \theta /0.0896\right)/\left[3.0\times \sin {\theta }^{0.8}+0.56\right]$$

(9)

$$S=\{\begin{array}{ccc}10.8\times \sin \theta +0.03& ,& \theta <9\%\\ 16.8\times \sin \theta -0.50& ,& 9\%\le \theta <14\%\\ 21.9\times \sin \theta -0.96& ,& \theta \ge 14\%\end{array}$$

(10)

where $L$ is the slope length factor, $S$ is the slope steepness factor, λ is the length of the horizontal projection of the slope, α is the slope length exponent, and θ is the slope (%) extracted from DEM.

The vegetation coverage can be estimated more accurately by establishing the relationship between the vegetation index and the vegetation coverage factor [42]. This experiment adopted the calculation formula of the vegetation coverage factor of Cai et al. [39], which can be written as follows:

$$f_c=\frac{NDVI-NDV{I}_{\min }}{NDV{I}_{\max }-NDV{I}_{\min }}$$

(11)

$$C=\{\begin{array}{ccc}1& ,& 0\le fc<0.1\%\\ 0.6508-0.3436\log fc& ,& 0.1\%\le fc<78.3\%\\ 0& ,& fc\ge 78.3\%\end{array}$$

(12)

where $f_c$ is the vegetation coverage (%), C is the vegetation coverage and management factor, $NDVI$ is the normalized vegetation index value, and $NDV{I}_{\max }$ and $NDV{I}_{\min }$ are the maximum and minimum values of $NDVI$ in the study area, respectively.

The p value reflects the soil and water conservation measures, indicating the ratio of soil lost after special measures to the amount lost when planting along the slope. The experiment adopted the method of calculating the P factor proposed by Lufafa et al. [43] in 2003, and the calculation formula is as follows:

$$P=0.2+0.03S$$

(13)

where S is the slope (%).

2.3.2. Hydrological Information Extraction

Based on DEM data and ArcGIS 10.2, the experiment extracted the river network of Jinghe Basin and divided the sub-basins of Jinghe Basin. The specific steps are as follows: ① use the fill tool of ArcGIS to fill the DEM, ② use the flow direction tool of ArcGIS to extract the water flow direction, ③ use the flow accumulation tool of ArcGIS to calculate the basin accumulation, ④ use the raster calculator tool of ArcGIS to obtain the river network of the basin, ⑤ use the stream link tool of ArcGIS to capture the river pour point, and ⑥ use the watershed tool of ArcGIS to divide the basin.

In order to obtain accurate basin river network information, the experiment set 15 different water catchment thresholds between 1 × 10⁴ m²–2.25 × 10⁶ m² in the raster calculator tool of ArcGIS and obtained the river network under the corresponding water catchment threshold. At the same time, the experiment calculated the density of the river network under different water catchment thresholds, which was shown in Figure 2. The fitting relationship between the water catchment threshold and the river network density is y = 0.022 × x^(−0.264) (R² = 0.9954), as shown in Equation (14), where x is the water catchment threshold and y is the river network density. By calculating the second derivative of Equation (14), we found no point in the definition domain where the second derivative of Equation (14) was 0 and where the derivative did not exist. With reference to previous experience [44], the power formula and the linear formula can be used to determine the appropriate catchment threshold for the river network extraction together. Therefore, this experiment calculated the combined solution for Equations (14) and (15) to obtain the catchment threshold suitable for the river basin.

$$y=0.022\times x^{\left(-0.264\right)}$$

(14)

$$y=-x+b$$

(15)

In the following research, this paper refers to Equations (14) and (15) as a formula set. By solving the formula set, we found that, when b = 0.177 and x = 0.14, the equation has a unique solution; when b < 0.177, the formula set had no answer; and when b > 0.177, the formula set had two solutions. Determining the catchment threshold in this way avoided the subjectivity of artificial determination and simplified the complexity of the trial. Compared with the Jinghe Network in the 1:1 million topographic maps of China [45], the river network extracted based on ArcGIS was consistent. Therefore, this paper takes a catchment area of 1.4 × 10⁵ m² (i.e., x = 1.4 × 10⁵ m²) for follow-up research work. In this case, the Jinghe Basin was divided into 58 sub-basins.

2.3.3. Dinf Algorithm

The Dinf algorithm takes the cell to be calculated as the center and diagonally divides the surrounding 8 cells (the diagonal direction of the central cell) in every 45°. The maximum slope and aspect in 8 intervals were calculated [46,47]. Then, the gap with the largest slope was selected and its direction r was regarded as the water flow direction of the current grid cell. The water flow was broken into 2 outflow directions: horizontal and vertical, corresponding to the surrounding 8 cells. An example is presented in Figure 3. The weight of flow distribution was determined by the angle between the slope direction and horizontal or vertical direction and the ratio of angle between the slope direction and diagonal direction. The water flow distribution plan can be seen in Figure 4.

Assuming that the maximum slope of the current grid unit (i, j) is in Δ2, the two downstream units involved are (i + 1, j) (i + 1, j + 1). If the current grid unit (i, j) is A, the flow calculation formula for the two downstream grids can be calculated using Equation (16):

$$\{\begin{array}{c}{A}_{i+1,j}=A_{i,j}\times F_{i+1,j}=A_{i,j}\times \frac{{\alpha }_2}{\alpha }\\ A_{i+1,j+1}=A_{i,j}\times F_{i+1,j+1}=A_{i,j}\times \frac{{\alpha }_1}{\alpha }\end{array}$$

(16)

2.3.4. Determination of Supply and Demand

① Determine the amount of supply and demand.

$$A_c=A_p-A_r$$

(17)

$$D_r=A_r$$

(18)

where $A_c$ (kg/m²) is the amount of supply from soil conservation services, $A_p$ (kg/m²) is the amount of potential soil erosion, $A_r$ (kg/m²) is the amount of actual soil erosion, and $D_r$ (kg/m²) is the amount of demand for soil conservation services.

In addition, this experiment normalized the supply and demand values of soil and water conservation services and divided them into five levels, namely low levels (0–0.2), low level (0.2–0.4), intermediate level (0.4–0.6), and high level (0.6–0.8) and a higher level (0.8–1) [48].

② Determine the supply and demand

$$A_w=A_c-D_r$$

(19)

where $A_w$ (kg/m²) is the difference between the supply and demand of soil conservation services. If the difference is greater than 0, it is a supply area. If it is less than 0, it is a demand area.

3. Results and Discussion

The experimental results showed that the average annual sediment transport modules in the Jinghe Basin from 2005 to 2015 was 1.597 kg/m². With reference to previous research results [49], the sediment transport ratio of the Loess Plateau is approximately 1. This means that the soil erosion modulus of the basin is roughly equal to the sediment transport modules. The RUSLE model calculated the average annual soil erosion modulus of the Jinghe Basin from 2005 to 2015 as 1.441 kg/m² with a ratio of about 1.1 and a correlation coefficient (R²) of 0.7392. The model accuracy verification is presented in Figure 5.

3.1. Temporal and Spatial Distribution Characteristics of Soil Erosion

The spatial distributions of soil erosion in the Jinghe Basin in 2005, 2010, and 2015 are presented in Figure 6. There was a slight degree of erosion in the southwest and southeast of the basin in 2005. The areas with severe erosion in 2010 were concentrated in the north-central and southeastern parts of the basin. The severely eroded areas in the basin in 2015 were located in the western and central-eastern regions of the basin. In 2005, 2010, and 2015, the average annual erosion moduli of the entire basin were 1.826 kg/m², 2.596 kg/m², and 1.298 kg/m² and the amounts of soil erosion were 8.292 × 10¹⁰ kg, 1.179 × 10¹¹ kg, and 5.897 × 10¹⁰ kg, respectively. From 2005 to 2015, soil erosion within the basin showed a trend of first deteriorating and then recovering. Using the scenario simulation method, we analyzed the impacts of rainfall and NDVI on soil erosion in 2010. In quantifying 2010 NDVI data, the 2005 and 2015 rainfall data were used as variables studying the effects of rainfall on the amount of soil erosion. Similarly, in quantifying rainfall data in 2010, the NDVI data in 2005 and 2015 were used as variables to analyze the influence of NDVI on the amount of soil erosion. The results of the scenario simulation method verifying the effect of NDVI data and rainfall data on the amount of soil erosion in 2010 are presented in

Table 1. Scenario simulation method verifies the impact of normalized difference vegetation index (NDVI) data and rainfall data on soil erosion in 2010.

	Rainfall in 2005	Rainfall in 2015	NDVI in 2005	NDVI in 2015
soil erosion amount (kg)	9.522 × 1010	8.882 × 1010	1.448 × 1011	1.400 × 1011

. Overall, rainfall had a more significant impact on soil erosion than NDVI. The average annual rainfall in the Jinghe Basin from 2005 to 2015 was similar. However, the standard deviation of the monthly rainfall in 2010 was large, indicating the variations in rainfall intensity and a large dispersion degree. It resulted in a high amount of soil erosion in the basin in 2010. The R² of soil erosion and rainfall and NDVI from 2005 to 2015 are 0.77 and 0.52, respectively.

According to the standards for classification and gradation of soil erosion SL190–2007 of the Ministry of Water Resources of the People’s Republic of China [50], soil erosion in the Jinghe Basin in 2005, 2010, and 2015 was mild, moderate, and mild, respectively. The proportion of the total erosion area for a strong erosion level in the basin decreased from 6.01% in 2005 to 1.6% in 2015, showing the characteristics of first increasing and then decreasing, with an overall decreasing trend. The moderately eroded area decreased from 21.2% in 2005 to 12.9% in 2015, with a monotonous decreasing trend. The proportion of areas with slight and mild erosion increased from 72.78% in 2005 to 85.47% in 2015, showing a trend of first decreasing and then increasing. It showed that, from 2005 to 2015, about 12.69% of the basin area changed from moderate or above erosion to slight or mild erosion, and the water and soil conservation measures in the basin were successful.

3.2. Temporal and Spatial Distribution Characteristics of Soil Conservation Service

The spatial distribution of soil conservation in the Jinghe Basin in 2005, 2010, and 2015 is presented in Figure 7. The spatial distribution pattern of the amount of soil conservation service was relatively consistent with the spatial distribution pattern of soil erosion. R² between the basin soil conservation rate and NDVI from 2005 to 2015 was 0.60, and that between the soil conservation rate and rainfall for the same time was 0.87. In 2005, the amount of soil conservation in the southwest and southeast of the basin was low, while at other parts, it was relatively high. In 2010, the north-central and southeastern parts of the basin exhibited a high amount of soil conservation. In 2015, these areas with high amounts of soil conservation were located in the western and central-eastern regions of the basin. In 2005, 2010, and 2015, the amounts of soil conservation in the Jinghe Basin were 3.09 × 10¹¹ kg, 4.58 × 10¹¹ kg, and 2.26 × 10¹¹ kg, accounting for 78.86%, 79.54%, and 79.30% of the potential erosion in the corresponding years. During the study period, although the amount of soil conservation fluctuated under the influence of rainfall, the ratio of soil conservation to potential soil erosion increased year by year.

In 2005, 2010, and 2015, the soil conservation amount for different land-use types was in the order of grassland > cultivated land > woodland > construction land > water area > unused land. During the study period, the proportion of soil conservation of grassland to the total soil conservation amount of the basin decreased from 47.1% in 2005 to 46.14% in 2015, showing a trend of first increasing and then reducing overall. The proportion of soil conservation amount of forest land increased from 7.91% in 2005 to 9.44% in 2015, showing a monotonous increase, while the proportion of cultivated land soil conservation amount decreased from 42.82% in 2005 to 41.71% in 2015. The ratio of soil conservation in construction land increased from 1.53% in 2005 to 2.02% in 2015, showing a trend of first decreasing and then increasing. The proportion of soil conservation amount in waters and unused land showed a stable trend. Compared with the ratio of land-use types, it can be found that, during the study period, the area under grassland and cultivated lands decreased by 0.07% and 0.44% while the area of forest and construction land increased by 0.14% and 0.33%, respectively. However, in general, the basin soil conservation measures achieved results during the study period, and the economy of the area was improved.

3.3. Space Flow Path of Soil Conservation Service

This paper defines soil conservation service flow as reducing the amount of sediment carried by surface runoff due to restoration of the basin ecosystem, with improvements in ecological benefits. Soil conservation service flow is a directional ecological service flow [9,51]. The carrier of soil conservation service is runoff or the material or energy carried with runoff as a carrier. For example, the No.1 sub-basin was divided into 2181 hydrological response units (HRUs) based on DEM, land use, and soil. Figure 8 shows the spatial flow directions of soil conservation service at the HRU scale in the No.1 sub-basin of the Jinghe Basin and the corresponding amount of soil conservation service in 2015. Figure 9 shows the spatial flow directions of soil conservation service at the channel scale in the No.1 sub-basin of the Jinghe Basin and the corresponding amount of soil conservation service in 2015.

This section takes the amount of soil conservation in 2015 as an example. The distribution of the supply levels and demand levels of soil conservation service in each sub-basin of the Jinghe Basin in 2015 is presented in Figure 10. Sub-basin No.1 to the north of the basin; No.16 in the middle of the basin; and No.26, No.38, No.45, and No.47 to the southwest of the basin had more quantity supplied than quantity demanded. In this sub-basin, the slope was relatively gentle, the land-use was mainly low-coverage grassland and cultivated land, and the soil type was highly active leached soil. After applying fertilizer, the soil fertility was high and the soil erosion rate was low. The slope in the No.16 sub-basin was relatively large. The land-use types were mainly grassland, woodland, and cultivated land with good vegetation coverage. The sub-basins of No.26, No.38, No.45, and No.47 had gentle slopes with good vegetation coverage. The effect of surface runoff was weak. For the No.2 sub-basin in the northern part of the basin; the No.9 sub-basin to the east; the No.19 sub-basin in the middle; the sub-basins No.24, No.36, and No.55 to the southwest of the basin; and the sub-basins No.27–29, No.31–35, No.37, No.41–44, No.46, No.48–49, No.50–53, and No.56–58 to the southeast of the basin, there was minimal difference between the amount of supply and the amount of demand in these sub-basins. The sub-basins No.2, No.9, No.19, and No.55 belonged to the lower supply (demand) area. The sub-basins No.24, No.27–29, No.31–35, No.37, No.41–44, No.46, No.48–53, and No.56–58 were all low supply (demand) areas. The No.36 sub-basin belonged to the medium supply (demand) area. The sub-basins No.3–8 in the northern part of the basin; No.10–15 sub-basins in the middle part of the basin; the sub-basins No.20, No.22, No.30, and No.40 in the northeast of the basin; and the sub-basins No.17–18, No.21, No.23, No.25, No.39, and No.54 to the southeast of the basin had more quantity demanded than the quantity supplied. The sub-basins No.3–8 and No.10–15 had steep slopes. The land-use types were mainly grassland and cultivated land with low vegetation coverage. The rainfall was high, and the erosion of soil particles by raindrops and runoff was strong, leading to severe soil erosion and nutrient loss. The sub-basins No.20, No.21–23, No.25, No.30, No.39–40, and No.54 received heavy rainfall and had medium vegetation coverage. The main soil types were black with thin layers, which were loose and easy to lose. The sub-basins No.17–18 belonged to the Loess Plateau area, with large undulations and rugged slopes.

According to Equation (19), the distribution of the supply area and demand area of the basin are presented in Figure 11. In 2015, there were five main supply areas for soil conservation service in the basin. The supply amount from these areas was greater than the demand. The remaining regions belonged to the demand area, i.e., the supply amount was less than the demand. There was a surplus of soil conservation service in these supply areas, leading to a flow phenomenon. Following Wang et al. [52], we assumed that the soil conservation service from the supply areas would first flow to the most adjacent area (called connected flow) (represented by the solid arrow in Figure 11). In the demand area, it was a deficit status without the flow phenomenon (called terminating flow), represented by the dotted arrow in Figure 11. In addition, this paper assumed that a sub-basin could receive soil conservation service from multiple sub-basins through the soil conservation service flow. Based on the above assumptions and the consideration of the elevation of each sub-basin and according to the flow law of matter from high to low, we obtained the spatial flow path of the soil conservation service in the Jinghe Basin in 2015 at the sub-basin scale (Figure 11) and the direction and quantity (

Table 2. Spatial flow path of soil conservation service at the sub-basin scale in Jinghe Basin in 2015.

Start Sub-Basin Number	End Sub-Basin Number	Flow Direction (°)	Flow Amount (kg)
1	2	240	1.89 × 109
2	3	130	2.15 × 108
16	19	220	1.29 × 109
26	22	90	5.43 × 108
38	36	87	6.62 × 109
45	38	25	3.42 × 109
47	40	55	7.23 × 108
55	47	315	2.97 × 108
27	31	135	3.41 × 107
31	35	210	1.03 × 108
35	44	155	6.17 × 107
44	46	250	8.32 × 107
48	44	0 (360)	9.39 × 107
50	48	300	2.15 × 108
52	56	120	8.68 × 107
56	57	200	9.88 × 107
57	58	105	1.93 × 107

).

We finally mapped the spatial flow paths of the soil conservation service on different scales in the Jinghe Basin in 2015 (Figure 12). Figure 11 shows the direction of the soil conservation service flow on the HRU-channel scale, and the direction and the amount of soil conservation service trend flow on the scale of the supply area to the demand area. These will provide support to improve the credibility of the soil conservation policy in the basin.

Combined with the socioeconomic data in 2015, we also calculated the value of ecological compensation in the Jinghe Basin. In 2015, the price of N fertilizer was 0.284 USD/kg, of P fertilizer was 0.08 USD/kg, and of K fertilizer was 0.67 USD/kg and the construction cost of the reservoir project was 0.87 USD/m³. The details of the calculation can be found in Reference [53]. According to the 2020 inflation rate (calculated according to the Chinese Consumer Price Index (CPI)), the total ecological compensation from the supply areas for soil conservation services in the Jinghe Basin (a total of 18 sub-basins) was 7.3 × 10⁵ USD. The total ecological compensation from the demand areas (40 sub-basins in total) was 8.1 × 10⁷ USD. The whole ecological compensation from the demand areas paid to the supply areas was 1.67 × 10⁵ USD. This is expected to provide a reference for the formulation of ecological compensation policies in the basin.

4. Discussion

Presently, different experts and scholars have different views on the flow of ecosystem services. Among them, there are two main understandings: (1) with emphasis on the process, the ecosystem service flow is a link between the ecosystem from the supply area to the demand area and the ecosystem service provided by the supply area is transferred to the beneficial area along a direction and path; (2) with emphasis on ultimate utility, the ecosystem service flow is the ecosystem service finally obtained by a human being and the ultimate realization of ecosystem services. Based on the above opinions, this paper conducted a quantitative assessment of soil conservation services in the Jinghe Basin and separated the path and flow of ecosystem service flows to describe the path of ecosystem service flow and to calculate the flow amount of ecosystem service flow. This paper realized the multi-scale analysis of soil conservation service in “HRU-channel-sub-basin” and studied the specific paths and flow amount of soil conservation service flows on different scales.

Soil erosion has always been one of the hot topics studied by scholars from various countries. China’s Loess Plateau is a typical area for soil erosion research. This paper took soil conservation service as the research object and took the Jinghe Basin in the Loess Plateau region of China as the case area to evaluate the soil conservation service of the Jinghe Basin from 2005 to 2015 based on the RUSLE model. This paper took 2005, 2010, and 2015 as representatives to describe the temporal and spatial characteristics of soil erosion and soil conservation services of Jinghe Basin. On this basis, this paper took 2015 as representative to simulate the soil conservation service flow on three scales including “HRU-river-sub-basin”, which is one of the innovations of this paper. Different from previous research, the division of HRUs, river channels, and sub-basins in this research was based on ArcGIS software, which is also one of the highlights of this article. Specifically, this paper used the azimuth to express the flow direction of the soil conservation service flow and gave specific values. This paper carried out a multi-scale expression of “HRU-channels-sub-basin” on the soil conservation services flow in 2015, which showed in more detail the spatial flow process of soil conservation services on different scales. This paper finally gave the general flow trend of soil conservation service throughout the basin, which can provide a scientific basis for formulating local comprehensive management policies and ecological safety policies in the future. Moreover, the rainfall erosivity factor (R) in the experiment was calculated based on Yangkun data. Compared with previous experiments, it avoided the influence of weather station data interpolation on the experiment. The soil conservation measure factor (P) in the experiment was based on the definite calculation formula proposed by LuFaFa [43] et al. Calculating the P factor through a certain formula effectively avoids the subjectivity of artificial assignment. However, the RUSLE model is an empirical model, which cannot fully demonstrate the process of soil erosion. In the future, we will try to simulate the service flow of soil conservation based on a mechanism model.

5. Conclusions

This paper proposed a new approach to the study of ecosystem service flow and realized multi-scale spatial flow simulation of soil conservation services in 2015. From 2005 to 2015, the correlation between the amount of soil erosion simulated by the RUSLE model and the measured value reached 0.74. During the study period, the amount of soil erosion in the Jinghe Basin showed a trend of first deepening and then reducing. Affected by rainfall, the amount of soil conservation showed a trend of first increasing and then decreasing, but the proportion of the amount of soil conservation to the amount of potential soil erosion during the study period showed an increasing trend. This paper took the soil conservation service in the Jinghe Basin in 2015 as an example and separately described the spatial flow path and flow volume of the soil conservation service at the scale of “HRU-channel-sub-basin”. Among them, the division of HRUs, channels, and sub-basins in the Jinghe Basin were all realized through the ArcGIS software. After comparing the supply and demand of soil conservation services, this paper divided the supply area and beneficiary area of soil conservation services in the basin in 2015 and calculated the ecological cost that the soil conservation service benefit area needs to provide to the supply area as 1.67 × 10⁵ USD. The experimental results can provide a certain scientific reference for the formulation of comprehensive management policies for the river basin.