Spatial–Temporal Variations of Water Ecosystem Services Value and Its Influencing Factors: A Case in Typical Regions of the Central Loess Plateau

Abstract

Water resources provide indispensable ecosystem services, which are related to human well-being and sustainable social development. Accurately measuring the water ecosystem services value (WESV), and then grasping its changing characteristics, is particularly important for solving water problems. In this study, the typical area of the central Loess Plateau location is taken as the research area. Based on remote sensing images and statistical data, the direct market method combined with the equivalent factor method was used to calculate the WESV including groundwater and surface water, which is of greatest originality. The temporal and spatial variation characteristics in 2010, 2015 and 2020 were analyzed. Then, four WESV driving factors including per capita GDP, population density, proportion of water areas, and water consumption were selected, and the geographically weighted regression (GWR) model was used to analyze the spatial distribution pattern and temporal variation of WESV’s response to the influencing factors. The results showed that WESV experienced a process of first decreasing and then increasing, which was mainly caused by Yulin. For the composition of WESV, the proportion of provisioning services value has increased, which caused the proportion of regulating services value to decrease. The correlations between four factors and WESV were different. The distribution pattern of the influences was spatially heterogeneous, which showed regular variations over time. These results indicate the necessity of WESV’s independent research and provide a realistic basis for ecological compensation in the Yellow River Basin.

1. Introduction

Ecosystems provide basic and necessary services for human survival and social functioning, namely ecosystem services (ESs). ESs are the continuous provision of ecosystem goods and services by ecosystems and their ecological processes [1]. However, factors such as population growth, industrialization and urbanization have led to a rapid increase in the demand for ESs, and ecosystems are facing unprecedented pressure [2]. A comprehensive and reasonable quantitative assessment of ecosystem services value (ESV) is necessary to alleviate the contradiction between supply and demand of ESs, manage effectively and formulate relevant policies. Water resources are an essential part of ecosystems. Whereas, due to the existence of water pollution, waste of water resources and climate change, the water ecosystem is facing greater pressure than other types of ecosystems. Therefore, it is crucial to study the water ecosystem services value (WESV) and the various factors that affect WESV.

Monetization of ecosystem service value is the most recognized and practical form of ESV, and the calculation methods can be divided into two categories. One is to adopt the relevant methods of traditional ecological economics or environmental economics. Most of these methods obtain the output of physical quantity based on statistical data or ecological model, and then they calculate the value of ecosystem services by combining the market value method, willingness survey method (CV) and revealed preference method. This kind of method has high requirements on data, parameters, model accuracy and method applicability, etc. Specific to the value of water resources ecosystem services, this kind of method is suitable for calculating the value of water resources with commodity attributes. The other is the equivalence factor method, which constructs the economic value equivalent per unit area of different ecosystems based on the division of ecosystem service functions and quantifies ESV in combination with the distribution of ecosystems. Compared with traditional methods, the equivalence factor method requires relatively fewer data, and the evaluation of ESVs is more comprehensive. Costanza et al. divided the global ESs into 17 species and calculated that the ESV of 16 biomes in the world was US$33 trillion per year [3], and rivers or lakes were one of the 16 biomes. Meanwhile, another contribution of this study is to provide the equivalent table of services value per unit area of 17 ecosystems in each biome at the global scale, which provides a reference for future research. In 2003, the Millennium Ecosystem Assessment (MA) conducted by the United Nations classified global ESs into four primary categories, namely Provisioning Services, Regulating Services, Supporting Services, and Cultural Services [4]. At the regional scale, Xie et al. improved Costanza’s research and obtained the equivalent factor of ecosystem service value suitable for China. Then, the ESV can be calculated combined with the area of each ecosystem. This study was accepted by a large number of Chinese scholars [5,6]. Whether it is a global scale or a regional scale, in the study of using the equivalent factor method to evaluate ecosystems, most of the rivers/lakes or waters are divided into one type of ecosystem for discussion. However, this method does not clarify the service value of groundwater, which is an important component of water resources. The most direct service value of groundwater is that it provides a part of water for production and domestic use, accounting for only 70% of the total global water use for agriculture, and 40% of agricultural water comes from groundwater [7]. Therefore, it is more comprehensive and accurate to combine the two methods to assess the WESV including surface water and groundwater.

It is not the ultimate goal of scholars to study the value of WESV; the more important goal is to study the relationship between WESV and the interaction of various factors. In order to deal with more severe water and environmental problems, exploring the influence of various factors on WESV has gradually become one of the research hotspots. Due to the social nature of water resources, social and economic factors have become one of the components that affect WESV. The characteristics of cities, populations, communities, and cultures [8] all have an impact on water resources and water ecology [9,10]. The urbanization level is a concentrated expression of the social and economic development degree, which profoundly affects the spatial distribution and potential functions of ecosystem services [11]. Many researchers have conducted related research in North China [12,13], Yangtze River Delta [14], Pearl River Delta [15], Southwest Mountainous [16,17] and Northwest arid regions of China [18,19]. WESV also has a significant response to changes in natural factors. Climate conditions [20,21,22], ecosystem types [23], and environmental quality [24] are the main influencing factors. Meanwhile, the coupling of many factors, such as nature, social economy and human activities, has a more realistic impact on WESV [25,26]. Many of the above studies have fully considered various factors and provided important references for explaining the changes caused by EVS or WESV. Unfortunately, the collinearity among some influencing factors and the nonstationarity in space have not attracted enough attention.

However, EVS exhibit spatial heterogeneity and spatial dependence with changes in geographic space due to differences in the socioeconomic development degree, natural resources, and geographic environment. Thus, incorporating geospatial aspects into the research scope is the key to addressing spatial heterogeneity. Geographically Weighted Regression (GWR) model is an effective tool for dealing with spatial heterogeneity, which is improved on the basis of ordinary least square [27]. The model incorporates the spatial location information as a coefficient into the regression equation and explores to eliminate the nonstationarity caused by spatial changes based on the fitted values of geographic element parameters [28]. The GWR model has been widely used in the fields of natural resources and ecological environment [29]. The water footprint has been extensively researched, including concepts, methods and applications, for better management of water resources and water ecology [30,31]. With the deterioration of ecological environment problems, it is necessary to analyze the evolution of ecological footprint and the spatial differences of influencing factors from the perspective of spatial heterogeneity [29]. In the related research on land use variation and ESV, the GWR model is used to solve the problem of spatial heterogeneity and compare with the OLS model [32]. In coastal counties of Mississippi and Alabama (U.S.), GWR was used in the estimation of the monetary value of distance to different waterfront types, in the extension to a traditional hedonic pricing method, and in analyzing the value of ecosystem services associated with waterfronts differed geospatially [33]. However, few researches utilize the GWR model to study WESV, which becomes the main content of this study.

In arid and semi-arid regions, water resources are very precious, which means that the ecological services provided by water resources play a vital role. Therefore, on the basis of accurately measuring the value of water resources ecosystem services, analyzing the impact of various factors on WESV is of great significance to effectively manage water resources and alleviate the contradiction between supply and demand of water resources ecosystem services.

Many studies have been conducted on the value of ecological services and their influencing factors. However, this study is more relevant. Specifically, the ecological service value of water resources is the object of this study. In addition, the scope of water resources is broader to include groundwater and surface water. This study serves the ecological compensation policy in China.

In this study, the typical area in the central Loess Plateau of China was taken as the research object, and the evaluation of water resources ecosystem service value and the analysis of the temporal and spatial changes of the influencing factors were carried out. The main works are as follows:

(1)

Combining the environmental economics method and the equivalent factor method, the WESVs of 25 counties (districts) in the study area including groundwater and surface water in 2010, 2015 and 2020 were calculated, and the distribution characteristics were analyzed;

(2)

Selecting the representative factors of nature, economy and society, and the applicability of OLS model and GWR model in studying the impact of each indicator on WESV was analyzed;

(3)

Using the more applicable GWR model, the spatial heterogeneity and spatial and temporal distribution of the effects of various factors on WESV were shown.

2. Materials and Methods

2.1. Study Area

The Loess Plateau (33°43′–41°160′ N, 100°54′–114°33′ E) covers an area of 640,000 km² in the upper and middle reaches of China’s Yellow River (Figure 1a) [34,35]. Most areas of the Loess Plateau belong to arid and semi-arid areas, with fragile environment, scarcity of water resources and serious soil erosion [36,37]. The study area of this paper is Yulin and Yan’an Cities (Figure 1b), with an area of about 79,957 km², accounting for 12.49% of the area of the Loess Plateau. In terms of location, the study area is located in the middle of the Loess Plateau, which is a representative area in the middle of the Loess Plateau. From the perspective of administrative division, the study area belongs to the north of Shaanxi Province. The study area includes 25 county-level administrative units (Figure 1c). The study area with large topographic relief and hilly gully, is the core area for controlling soil and water loss in the Yellow River Basin. The annual average temperature of Yulin City in the north of the study area is about 10.5 °C, and the average annual precipitation is about 400.00 mm. From north to south, the landform gradually transits from sandy land to gullies and hills [38,39]. It is an important energy and chemical base in China. Yan’an City in the south of the study area belongs to the hilly area of the Loess Plateau, which is high in the northwest and low in the southeast. The annual average temperature is about 7.70~10.60 °C, and the annual average precipitation is about 500.00 mm [40,41,42].

2.2. Data Sources

The data involved in this study include economic indicators, social indicators, and natural indicators of Yulin and Yan’an Cities in 2010, 2015 and 2020. It should be noted that due to the impact of COVID-19, some economic and social indicators in 2020 are replaced by data from 2019, including population, population density, urbanization rate, per capital GDP, and gross product of primary industry.

The water consumption comes from the water resources bulletin of Yulin and Yan’an Cities (2010, 2015 and 2020). Through reclassification, water consumption is divided into three categories: agricultural water consumption, residential water consumption and nonresidential water consumption. The price of water is obtained from the research of the water supply department.

GDP, GDP per capita, population, population density and urbanization rate are from the statistical yearbooks of Yulin and Yan’an (2010, 2015 and 2019).

The water area and the proportion of water area are calculated through the statistical calculation of the spatial distribution data of the national land use type remote sensing monitoring provided by the Resource and Environmental Science and Data Center, Chinese Academy of Sciences (http://www.resdc.cn/, accessed on 2 July 2018). The resolution of the data is 30 m × 30 m, which is generated by manual visual interpretation using the remote sensing images of Landsat TM of various phases of the US Landsat as the main data source.

2.3. Methods

2.3.1. Classification and Calculation Method of WESV

According to the classification of ESs by the Millennium Ecosystem Assessment (MA 2005) carried out by the United Nations [4], combined with the research practice in China [43] and the practicability of ecological compensation in the study area, this study divided WESV into three categories. The reason why the cultural services value was not considered is that at this stage, ecological compensation in China rarely involves such services value. Considering the ways of WES and the availability of data, the three categories were subdivided into 9 specific services (

Table 1. Classification of WESV.

Categories of WESV	Work Mode	Content	Symbol
Provisioning Services Value	Direct supply	Water supply	V1
	Indirect supply	Aquatic product	V2
		Material	V3
Regulating Services Value	Direct regulation	Hydrological regulation	V4
		Water purification	V5
	Indirect regulation	Climate regulation	V6
		Air regulation	V7
Supporting Services Value	Indirect support	Soil conservation	V8
		Biodiversity	V9

). The calculation method of WESV can be expressed formally as follows:

$$V_w={\displaystyle \sum _{n=1}^9{V}_n}$$

(1)

where V_w is the value of WESV; V_n is the n-th sub-category of WESV; n = 1, 2, …, 9.

2.3.2. Market Value Method

Whether groundwater or surface water, the most direct and important service is water supply for living and production. The value of water services is reflected in the water price [44]. Therefore, this study adopts the market value method to calculate the water supply services value of water resources. According to the type of water price, water is divided into agricultural water, residential water and nonresidential water, among which nonresidential water mainly refers to industrial, business service water, administrative institution water, municipal water, etc. The formula for calculation is as follows:

$$V_1={\displaystyle \sum _{j=1}^3{Q}_j\times }{P}_j,$$

(2)

where V₁ is the value of water supply services; Q_j and P_j are the water consumption and water price of the j-th type of water; j = 1, 2, 3 represent agricultural water, residential water and nonresidential water, respectively.

The market value method is also used to calculate the indirect supply value of aquatic products provided by water resources. Due to the variety of aquatic products and the different prices, this study uses the fishery output value in the statistical yearbook [45,46,47,48] as the value V₂ of aquatic product supply services, and the calculation method also adopts the market value method.

2.3.3. Equivalent Factor Method

Water resources are the most basic elements to maintain the normal operation of ecosystems and social systems, and most of its ecosystem services cannot be measured by market value. The ESV equivalent factor method is obtained from the study of Xie G. [5] by calculating the ecological service value per unit area of each ecosystem in China, which has been widely recognized and applied [49]. The equivalence factor is defined as a relative quantity, which represents a relative value of ecosystem services relative to the economic value of grain output in that year. In the case of the study area, due to the difference in planting structure, the economic value of grain output will change accordingly. The ecological service value of the study area can be calculated by the follows:

$$V_{per}=\frac{1}{7}\times \frac{\left(Y_1\times P_1+Y_2\times P_2\right)}{A_1+A_2},$$

(3)

where V_per is the value of ESs per unit area; corn and potato are mainly planted in the study area, Y₁ and Y₂ are the yield of corn and potato, respectively, which comes from “Shaanxi Statistical Yearbook” [45,46,47,48]; P₁ and P₂ are the corresponding grain prices, which comes from “Compilation of National Agricultural Product Cost and Benefit Data” [50,51,52]; A₁ and A₂ are the planting area of corresponding grain; 1/7 is the ratio of the economic value provided by natural ecosystems without artificial inputs to the economic value provided by existing farmland. The calculated values of ESs per unit area in 2010, 2015 and 2020 were 1457.16 CNY/hm², 1318.81 CNY/hm² and 1199.63 CNY/hm², respectively.

The equivalence factor adopts the corresponding part of the water ecosystem in the equivalence coefficient table of ecosystem service value per unit area calculated by Xie G. [43] (

Table 2. Equivalent coefficients of WESV, and the WESV per unit area in different years.

	Provisioning Services	Regulating Services				Supporting SERVICES
	Material	Hydrological Regulation	Water Purification	Climate Regulation	Air Regulation	Soil Conservation	Biodiversity
Equivalent coefficient (Em)	0.23 1	102.24 1	5.55 1	2.29 1	0.77 1	0.93 1	2.55 1
2010 WESV per unit area	335.15	148,979.93	8087.23	3336.89	1122.01	1355.16	3715.76
2015 WESV per unit area	1055.04	303.33	134,834.65	7319.37	3020.06	1015.48	1226.49
2020 WESV per unit area	275.92	122,650.62	6657.97	2747.16	923.72	1115.66	3059.07

¹ The equivalent coefficient was quoted from Ref. [43].

). At the same time, the value of WESV per unit area in 2010, 2015 and 2020 was calculated (Table 2).

Then, the value of services other than Water supply and Aquatic product could be calculated as follows:

$$V_m=V_{per}\times A_w\times E_m,$$

(4)

where m = 3, 4, 5, …, 9; V_m is the value of other services besides Water supply and Aquatic product; A_w represents the watershed area for the study year; E_m is the equivalent coefficient of different service in Table 2.

2.3.4. Geographically Weighted Regression

The first law of spatial geography shows that the correlation between ground objects gradually increases as the distance decreases [53,54]. Inevitably, spatial correlation and spatial heterogeneity coexist. GWR achieves better results when using local smoothing to deal with the problem of spatial heterogeneity [55,56]. GWR was based on kernel-weighted regression. Instead of estimating global values for regression parameters, GWR allows these parameters to be derived for each location separately [57]. The model can be expressed as

$$y_i={\beta }_0\left({\mu }_i,{\nu }_i\right)+{\displaystyle \sum _k{\beta }_k\left({\mu }_i,{\nu }_i\right)}{x}_{i,k}+{\epsilon }_i,$$

(5)

where y_i is the explained variable; (μ_i, ν_i) are the coordinates of the target area i; β₀ (μ_i, ν_i) is the intercept; x_i,k is the value of the explanatory variable x_k on the target area i; the value of the function β_k (μ_i, ν_i) at geographic location i is β_k (μ_i, ν_i); k is the number of explanatory variables; ε_i is the random disturbance term, i.e., the error. The coefficient of each sample point is a parameter estimate obtained by weighting the adjacent observations [58], and the expression is as follows:

$$\widehat{\beta }(u_i,v_i)={(X^{T}W(u_i,v_i)X)}^{-1}{X}^{T}W(u_i,v_i)Y,$$

(6)

where $\widehat{\beta }(u_i,v_i)$ is the parameter estimate of the local coefficient of the i-th sample with coordinates (μ_i, ν_i); X and Y are the vectors of the explanatory and the dependent variables; $W(u_i,v_i)$ is the weight matrix, which is usually calculated by the Gaussian function of the distance decay function (kernel function).

ArcGIS software has developed the related functions of OLS and GWR into convenient tools. Therefore, this study uses ArcGIS 10.4 to perform related calculations.

3. Results

3.1. WESVs and Characteristics in 2010–2020

3.1.1. Temporal and Spatial Variation Analysis of WESV

Taking the county-level administrative region as the minimum calculation unit, and adopting appropriate methods according to different service types, the WESVs of typical areas in the middle of the Loess Plateau in 2010, 2015 and 2020 were calculated (

Table 3. WESV of the typical area in the central Loess Plateau. unit: million CNY.

City	County	2010	2015	2020
Yulin	Yuyang	866.28	1071.75	1386.48
	Shenmu	3336.85	2765.08	3237.86
	Fugu	1294.82	900.46	1194.77
	Hengshan	589.99	686.95	811.72
	Jingbian	1086.09	1070.41	1133.69
	Dingbian	918.33	905.37	1036.14
	Suide	264.05	249.26	235.08
	Mizhi	114.52	119.55	125.89
	Jiaxian	404.83	401.90	363.29
	Wubu	161.24	147.49	143.38
	Qingjian	340.77	309.77	317.19
	Zizhou	54.41	64.33	77.94
	Subtotal	9432.18	8692.32	10,063.43
Yan’an	Baota	352.40	366.65	427.71
	Yanchang	215.76	333.62	302.61
	Yanchuan	324.42	228.90	233.01
	Zichang	92.63	101.52	117.77
	Ansai	126.53	161.89	152.78
	Zhidan	234.91	251.12	270.24
	Wuqi	281.47	319.77	345.85
	Ganquan	213.12	202.37	183.24
	Fuxian	507.99	468.73	542.46
	Luochuan	221.64	221.65	274.95
	Yichuan	96.23	96.43	129.86
	Huanglong	15.50	21.05	24.65
	Huangling	255.62	280.84	310.91
	Subtotal	2938.22	3054.54	3316.04
Total		12,370.42	11,746.86	13,379.48

). The changes in counties among the three years have also been shown (Figure 2). Overall, the WESV in the study area showed a trend of first decreasing and then increasing. The WESV decreased from 12,370.42 million CNY in 2010 to 11,746.86 million CNY in 2015, which showed a drop range of 5.04%. Then, it increased to 13,379.48 million CNY in 2020. The increasing range was 13.90% from 2015 to 2020, and 8.16% during the entire study period of 2010–2020.

It can be seen from

Table 4. Structural of WESV in 2010, 2015, 2020.

WESVs	Content	2010		2015		2020
		WESVs	Proportion	WESVs	Proportion	WESVs	Proportion
Provisioning Services	Water supply	1260.94	10.19%	1963.99	16.72%	3729.11	27.87%
	Aquatic production	45.08	0.36%	122.14	1.04%	188.06	1.41%
	Materials	22.21	0.18%	19.40	0.17%	19.00	0.14%
Regulating Services	Hydrological regulation	9874.52	79.83%	8621.80	73.40%	8444.72	63.12%
	Water Purification	536.03	4.33%	468.03	3.98%	458.41	3.43%
	Climate Regulation	221.17	1.79%	193.11	1.64%	189.15	1.41%
	Air Regulation	74.37	0.60%	64.93	0.55%	63.60	0.48%
Supporting Services	Soil Conservation	89.82	0.73%	78.43	0.67%	76.82	0.57%
	Biodiversity	246.28	1.99%	215.04	1.83%	210.62	1.57%

that, from 2010 to 2020, some districts and counties in Yan’an City in the southern part of the study area experienced a decrease in WESV, but the trend was rising, with a rate of 3.96%. The leading area that led to the decreasing trend in 2015 was Yulin City. From 2010 to 2015, the WESV in Yulin City decreased by 7.84%, which was significantly greater than the entire study area decreasing rate. Only four counties of Yuyang, Hengshan, Mizhi and Zizhou showed an increase in WESV, and the remaining two-thirds of the counties presented a decreasing trend; the decreasing degree in Fugu reached the highest 30.46%. These phenomena indicated that the ecological service function provided by water resources in Yulin City was in the stage of degradation from 2010 to 2015, which also indirectly reflected the trend of ecological environment deterioration.

From 2015 to 2020, WESV in Yan’an City was still in a growth trend as a whole, and the growth rate increased to 8.56%. Only Yanchang, Ansai and Ganquan experienced negative growth. During the same period, WESV in Yulin City increased rapidly, with an overall increasing rate of 15.77%, and 75% of the counties in the jurisdiction were on the rise, mainly due to the low value of WESV in 2015, which showed that the water ecosystem and its service functions in Yulin City were in a significant recovery situation in 2015–2020.

From the whole study period, the increasing rate of WESV in the southern part of the study area was significantly greater than that in the northern part. In 2020, WESV in Yan’an increased by 12.86% compared with 2010, while Yulin increased by only 6.69% in this decade.

3.1.2. Variation Analysis of WESV Structural Characteristics

Table 4 and Figure 3 exhibited the changes in the value and proportion of nine subservices that comprise WESV from 2010 to 2020. Among the three first-level classifications of WESV, Regulating Services accounted for the largest proportion, but the proportion gradually decreased from 86.55% in 2010 to 68.44% in 2020. The Provisioning Services proportion has increased year-by-year, from 10.73% to 29.42% in the past decade.

The Supporting Services proportion has decreased steadily, from 2.72% in 2010 to 2.14% in 2020.

In the subcategories, only the proportions of Water Supply and Food Production have increased year-by-year, and the proportions of other ESV have decreased to varying degrees, which showed that the development of economy and society has caused a sharp increase in the demand for water resources and aquatic products, and at the same time, the regulating function of water resources has gradually weakened.

3.2. Analysis of Spatial Heterogeneity Characteristics of Each Impact Factor

3.2.1. Identification of Influencing Factors

This paper selected 12 influencing factors from two aspects of society-economy, and natural environment. The values of NDVI, Forest and grass area, Forest and grass coverage ratio, and Proportion of water area were directly or indirectly obtained from the Resource and Environmental Science and Data Center, Chinese Academy of Sciences (http://www.resdc.cn/, accessed on 2 July 2018). All other influence factors were obtained from the statistical yearbook [45,46,47,48]. For the three different years of 2010, 2015 and 2020, the collinearity investigation of influencing factors and the significance test between them and the explained variables were carried out respectively, which is the first step in building a model. The specific work was to construct models for different parameter combinations by the OLS model in Arcgis 10.4. Then, parameters, that did not meet the criteria for variance inflation factor (VIF) and significance tests were excluded one by one gradually. The results (

Table 5. Statistical test results and the selection of influencing factors.

Classification	Influence Factor	Unit	Statistical Test Results
Social and economy	Population	104	VIF > 7.5
	Population density	person/km2	VIF < 7.5, significant
	Urbanization rate	%	VIF < 7.5, no significant
	Per capital GDP	yuan	VIF < 7.5, significant
	Gross product of primary industry	yuan	VIF > 7.5
	Water consumption	104 m3	VIF < 7.5, significant
Natural environment	Rainfall	mm	VIF > 7.5
	NDVI	/	VIF > 7.5
	Forest and grass area	hm2	VIF > 7.5
	Forest and grass coverage ratio	%	VIF < 7.5, no significant
	Proportion of water area	‰	VIF < 7.5, significant
	Area of soil erosion control	hm2	VIF < 7.5, no significant

) showed that in the three years, only four influencing factors satisfy the two conditions of being VIF < 7.5 and significant at the same time. Therefore, we selected per capita GDP, population density, the proportion of water areas, and water consumption as the research objects.

3.2.2. Effect of the GWR Model

In order to illustrate the applicability of GWR, this study first used the OLS model to simulate. As shown in (

Table 6. Statistical test of OLS and GWR in 2010, 2015, and 2020.

	2010		2015		2020
	OLS	GWR	OLS	GWR	OLS	GWR
R2	0.758	0.856	0.792	0.953	0.850	0.906
Adjusted R2	0.709	0.789	0.751	0.915	0.820	0.865
AICc	377.181	373.591	364.539	351.095	365.338	362.201

), the R² and Adjusted R² of the GWR model in the three years are larger than those of the OLS model, indicating that the simulation effect is more accurate and representative. Meanwhile, the AICc of the GWR model is smaller than that of the OLS model, and the difference is greater than 3.0, which demonstrated that the GWR model is more applicable. Therefore, the GWR model can accurately explain the relationship between WESV and each explanatory variable and the spatial heterogeneity.

3.2.3. Per Capita GDP

Figure 4 displayed the spatiotemporal distribution of the per capita GDP impact on WESV. The darker the color, the greater the positive effect. In 2010 (Figure 4a) and 2020 (Figure 4c), WESV and Per Capita GDP were positively correlated. The magnitude of the coefficient indicated that, in general, the 2020 Per Capita GDP had a stronger impact on WESV. In 2010, the correlation relationship gradually increased from west to east, with Huanglong, Huangling and Yanchang in the southeast, Wubao in the east and Fugu in the northeast being the most influential. In 2020, the impact of Per Capita GDP on WESV gradually increased from southwest to northeast, and WESV in Fugu, Shenmu and Jia County was more sensitive to Per Capita GDP. Because Yulin is the base of energy and chemical industry in China, the economy had developed rapidly after 2008, and the demand for water had increased. The combined effect of the two led to changes in the spatial distribution from 2010 to 2020. Nevertheless, in 2015 (Figure 4b), the WESV of 11 districts and counties in the northern and central parts of the study area showed a negative correlation with Per Capita GDP, and the negative correlation effect was strongest in the northeast. This illustrated that during the period from 2010 to 2015, the protection of water resources and water environment in the northern region was neglected due to the great economic development, which was consistent with the reality of the reduction of the water area in this region. It also showed that from 2015 to 2020, the water resources condition in the northern part of the study area and the ecological services provided had been greatly restored.

3.2.4. Population Density

In a certain area, population density determines the demand for WES. Figure 5 displayed the spatiotemporal distribution of WESV response to Population Density. Throughout the study period, WESV exhibited a negative correlation with population density, and the negative correlation gradually weakened with time. From 2010 to 2015, the more northerly the geographical location is, the more sensitive WESV is to changes in population density. Compared with Figure 5a,b, although the overall negative correlation was slightly enhanced, the area of region with the lowest level of negative correlation was increased. As shown in Figure 5c, the spatial distribution pattern of the impact of population density on WESV had fundamentally changed in 2020, and the coefficient increased from the southeast to the outside. The larger the coefficient, the weaker the negative correlation.

3.2.5. Proportion of Water Areas

WESV is largely determined by the water areas. The spatial distribution pattern of the proportion of water areas affecting WESV and the temporal change of this distribution pattern were analyzed (Figure 6). Looking at the entire study area, although the pattern distribution has changed, all coefficients showed that the proportion of water area has a positive impact on WESV. As shown in Figure 6c, the impact pattern in 2010 is that the degree of impact decreased from north to south. In 2015, the impact pattern evolved into that population density had the lowest impact on the southeastern region, and gradually increased in the west and north directions. In 2020, the gravity center of the impact continued to shift eastward, showing the phenomenon that the west was large and the east was small. The two counties with the least impact appeared in Shenmu and Fugu in the northeast, and these two districts and counties belonged to the most affected areas in 2010 and 2015.

3.2.6. Water Consumption

Water consumption is a key factor affecting the service value of Water Supply, and it is also the most direct manifestation of the ecological service value of water resources. Figure 7 exhibited the change in the spatiotemporal distribution of water consumption affecting the degree of WESV and showed that there is a positive correlation between water consumption and WESV. In 2010, the correlation coefficient between WESV and water consumption in the study area showed an increasing state from west to east. In 2015, the degree of influence showed a circular increase from the western to the eastern region (Figure 7b). The increase in the value range of the coefficient indicated that the differences between districts and counties were expanding. WESV in the southern region was gradually affected by changes in water consumption. With the passage of time, in 2020 (Figure 7c), the spatial distribution of WESV affected by water consumption had a strong regularity, decreasing in a stepwise manner from south to north. However, the coefficient ranges were less discrete.

4. Discussion

4.1. Necessity to Assessing WESV

The significance of evaluating the value of ecosystem services is to better manage the ecological environment and natural resources to achieve sustainable development. The water issue is prominent today, so it is more practical to discuss the value of water resources ecosystem services. In China, the value of ecosystem services is always in the form of the upper limit of the ecological protection compensation standard. The water ecosystem services value is often calculated as part of a comprehensive ecosystem service value and is rarely discussed in isolation. Even in the basin ecological compensation, only the fluctuation of the direct use value caused by the change of water quantity is calculated, and how the ecosystem service value of water resources including groundwater changes is not fully explored. In 2019, the ecological protection and high-quality development of the Yellow River Basin was established as one of China’s major national strategies. Ecological compensation is one of the key tasks of the strategy. The Loess Plateau, located in the middle reaches of the Yellow River, is the main source of sediment in the Yellow River. The region has scarce water resources, a large population, and an urgent need for development. Hence, from the perspective of the integrity of water resources, this study selects typical regions to assess WESV and analyzes the temporal and spatial variation characteristics of WESV, which can provide a basis for the ecological compensation development.

4.2. The Spatiotemporal Distribution of WESV Response to Different Influencing Factors

Previous studies have shown that the value of water resources is affected by multiple factors, such as water quantity, water quality, use of water resources, economic development, and educational level of residents [59,60].

With the deepening of relevant research, the water resources value has been extended to the value of ecosystem services that water resources can provide for human well-being. Correspondingly, WSVE is also affected by economic, social, natural and cultural factors. Many researches have explored the mechanism by which ESV is affected by various driving factors, including land use change, socioeconomic development indicators, and human acceptance willingness. Furthermore, the distribution characteristics of the sensitivity of ESV to various influencing factors at different temporal and spatial scales were obtained. WESV is an important part of ESV and is also affected by various factors. As mentioned above, the extent to which WESV is affected by external factors should be studied separately.

With the rapid economic development and the sharp expansion of cities, the differences between different regions in natural resources, economy, society and culture are gradually increasing. In this study, the central part of the Loess Plateau was selected as the research area, and four main influencing factors were selected. The GWR model was used to analyze the influence of per capita GDP, population density, the proportion of water areas, and water consumption on WESV in different regions.

Per Capita GDP is an important parameter to measure the development degree of a region. The rise in prices and the massive consumption of water resources brought about by economic development will lead to an increase in WESV. Meanwhile, WESV will greatly decreased because of the ecological degradation caused by development, that is, the reduction in water area. Therefore, the relationship between water conservation and economic development should be balanced. The effect of per capita GDP on ESV was confirmed by Song F. [61] in a study on the value of wetland ecosystem services. Dai X. concluded that per capita GDP was negatively correlated with ESV in Chengdu [62]. As can be seen from Figure 4, per capita GDP showed time instability. The study area has undergone a process of environmental damage and recovery during a 10-year development period, which is consistent with the process in China. During 2015–2020, the construction of ecological civilization has been raised to a new level. At last, a pattern in which the degree of economic development roughly matches WESV is formed. Many studies have found that per capita GDP has an unstable impact on ESV, by influencing the level of awareness, willingness to protect the environment, and investment in environmental protection [63,64,65].

It is generally believed that the increase in population density is accompanied by the expansion of human activity areas, which will encroach on other land-use types. The water area is also experiencing depletion of rivers and lakes due to the massive depletion of surface water while the water area was occupied. Under the dual effect of the two, WESV is bound to decrease, which also explains the negative correlation between population density and WESV. It was confirmed in Chen Y.’s study of the relationship between population density and agro-ecosystem services [66]. Jiang Z. also obtained a consistent conclusion in his study of the South Four Lakes that population density was negatively correlated with Esv [67]. In recent years, the population growth rate in the study area has slowed down as in the whole country. WESV is less sensitive to population density.

In the same region, WESV will increase with the increase in water area, which can also be explained by the calculation method adopted in this study. The equivalent factor method is to calculate the corresponding ESV according to the type of land use. Figure 6 showed that the proportion of water area has the greatest influence on WESV. During the initial stage of the study, the drier regions in the north were more sensitive. With the comprehensive management of the Mu Us Desert, the environment in the northern region has been gradually improved, and the scarcity of water resources and the vulnerability of the water environment have been alleviated. The center of influence then transferred to the west, where the pace of governance was relatively slow.

In 2010 and 2015, the impact of water consumption on WESV had a certain regularity in space, but the spatial instability that occurs cannot be ignored. Figure 7a,b can well demonstrate this phenomenon. From 2015 to 2020, the state strictly controlled the regional water consumption, and the strictest water resources management system was successfully implemented locally [68], which well explained the phenomenon that the spatial distribution pattern in 2020 was more regular.

There is spatial heterogeneity in the sensitivity of WESV to explanatory variables. Using the GWR model can more clearly show the spatial pattern of the influence degree and its evolution trend over time.

4.3. Limitations

Based on the water area obtained from remote sensing image statistics and the water consumption obtained from statistical data, this study proposed a method to quickly calculate WESV and analyze the temporal and spatial distribution of driving factors. However, the impact of water quality on WESV does not depend on water consumption or water area. Deterioration of water quality not only reduces the quality of services provided by the water resource but also incurs additional remediation costs [69]. This study did not include water quality as the basic data in the WESV assessment for three reasons. First, different types of WESVs have different sensitivity to water quality, and its mechanism of action and degree of influence are not clear, which will cause uncertainty in the calculation of WESV. Second, the equivalent factor method used in this paper is calculated from the ecological service value and the economic value of grain, and the output of grain is affected by water quality. That is to say, this study considered the impact of water quality on WESV to a certain extent. Third, the spatial and temporal scales of this study determined that water quality will not have an essential impact on WESV. In space, the county-level administrative region is the smallest research unit; in terms of time, a year is the minimum span of time. From this point of view, the water quality in the region is relatively stable, and this is also verified by the water quality data released by the environmental department. Therefore, the WESV calculated in this study is still representative and reliable without considering water quality.

5. Conclusions

In this study, taking a typical area in the central Loess Plateau as the study area, WESV was explored and assessed, and the spatiotemporal relationship between WESV and driving factors was analyzed. The conclusions obtained are as follows:

(1)

Considering surface water and groundwater as a whole, a WESV calculation method based on multi-method fusion of multiple data sources such as remote sensing images is proposed.

(2)

In total, the WESV in 2020 is 8.16% higher than that in 2010. However, during the study period, the WESV in the typical area of the central Loess Plateau experienced a process of first decreasing and then increasing. The main factor leading to this phenomenon was that the WESV of Yulin City in the north of the study area decreased by 7.84% in 2015 compared with 2010. From the perspective of the WESV composition structure, the proportion of regulating services has decreased by 18.11% in the past 10 years, but its proportion is still the largest; the proportion of provisioning services has increased year by year to 29.42%; the proportion of supporting services has decreased steadily, from 2.72% in 2010 to 2.14% in 2020.

(3)

Considering the spatial heterogeneity, the GWR model has better applicability. Among the four influencing factors, the proportion of water area and water consumption showed a positive correlation with WESV throughout the study period. Population density is negatively correlated with WESV. Per capita GDP was positively correlated with WESV in both 2010 and 2020. In 2015, there was a negative correlation between per capita GDP and WESV in the northern and western districts and counties of the study area.

(4)

With the passage of time, the spatial distribution pattern of the influence of the four factors on WESV has changed, and the evolution directions are different. The response degree of WESV to per capita GDP has evolved gradually from low in the west and high in the east to low in the southwest and high in the northeast, and the laddering nature is stronger. The water area plays a leading role in the size of WESV, and the center of gravity of the influence of the water area proportion shifts from the north to the west, with a strong decreasing law. The distribution of the influence of water consumption on WESV has experienced three evolutionary stages: small in the west and large in the east, increases in a circular shape from the west to the outside, small in the north and large in the south.