*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). The calculation method of WESV can be expressed formally as follows:

$$V\_w = \sum\_{n=1}^{9} V\_n \tag{1}$$

where *Vw* is the value of WESV; *Vn* is the *n*-th sub-category of WESV; *n* = 1, 2, . . . , 9.


**Table 1.** Classification of WESV.

#### 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 = \sum\_{j=1}^{3} Q\_j \times P\_{j\prime} \tag{2}$$

where *V*<sup>1</sup> is the value of water supply services; *Qj* and *Pj* 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–48] as the value *V*<sup>2</sup> 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\_{pcr} = \frac{1}{7} \times \frac{(\chi\_1 \times P\_1 + \chi\_2 \times P\_2)}{A\_1 + A\_2},\tag{3}$$

where *Vper* is the value of ESs per unit area; corn and potato are mainly planted in the study area, *Y*<sup>1</sup> and *Y*<sup>2</sup> are the yield of corn and potato, respectively, which comes from "Shaanxi Statistical Yearbook" [45–48]; *P*<sup>1</sup> and *P*<sup>2</sup> are the corresponding grain prices, which comes from "Compilation of National Agricultural Product Cost and Benefit Data" [50–52]; *A*<sup>1</sup> and *A*<sup>2</sup> 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/hm2, 1318.81 CNY/hm2 and 1199.63 CNY/hm2, 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). 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\_{\rm m} = V\_{\rm per} \times A\_{\rm w} \times E\_{\rm m} \tag{4}$$

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


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

<sup>1</sup> The equivalent coefficient was quoted from Ref. [43].

#### 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(\mu\_{i\prime}\nu\_i) + \sum\_k \beta\_k(\mu\_{i\prime}\nu\_i)\mathbf{x}\_{i,k} + \varepsilon\_{i\prime} \tag{5}$$

where *yi* is the explained variable; (*μi*, *νi*) are the coordinates of the target area *i*; *β*<sup>0</sup> (*μi*, *νi*) is the intercept; x*i,k* is the value of the explanatory variable *xk* on the target area *i*; the value of the function *β<sup>k</sup>* (*μi*, *νi*) at geographic location *i* is *β<sup>k</sup>* (*μi*, *νi*); *k* is the number of explanatory variables; *ε<sup>i</sup>* 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:

$$\beta(u\_{i\prime}, v\_i) = \left(X^T W(u\_{i\prime}, v\_i)X\right)^{-1} X^T W(u\_{i\prime}, v\_i) Y\_{\prime} \tag{6}$$

where *β*ˆ(*ui*, *vi*) 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*(*ui*, *vi*) 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.
