*2.1. Study Area*

The study area is the upper reach of the Beidao Hydrological Station in the Wei River (Figure 1). The study area is located at 103◦970−106◦42<sup>0</sup> E and 34◦170−36◦19<sup>0</sup> N, with a catchment area of about 25,000 km<sup>2</sup> , accounting for about 19% of the total area of the Wei River Basin. The study area is a Loess hilly and gully area with an altitude of 1079–3934 m. The climate is a continental monsoon climate, with a cold and dry winter and a hot and rainy summer. The average annual precipitation from 1965 to 2017 was 491 mm, and the precipitation was unevenly distributed throughout the year, mostly concentrated in July-September. The annual average temperature is 7.8–13.5 ◦C, the extreme maximum temperature is 42.8 ◦C, the extreme minimum temperature is −28.1 ◦C and the annual potential evapotranspiration is 700–1200 mm.

## *2.2. Data*

In this study, the ground-measured meteorological data include the daily data of 6 meteorological stations and 16 rainfall stations in the study area. The data of the meteorological stations include daily rainfall, temperature, sunshine duration, relative humidity, air pressure and average wind speed. The measured daily runoff data from 2001 to 2017 come from the Beidao Hydrological Station. The location of the stations is shown in Figure 1. The measured flow data and meteorological data come from the *Yellow River Basin Hydrological Yearbook* and the China Meteorological Data Network, respectively. The vegetation data from 2001 to 2017 were derived from the NDVI dataset provided by the MODIS/Terra website. The dataset was the MOD13Q1 product with a spatial resolution of 250 m and a temporal resolution of 16 days.

The DEM data used to build the hydrological model came from the ASTER GDEM data provided by Geospatial Data Cloud (http://www.gscloud.cn, accessed on: 1 October 2022), with a spatial resolution of 30 m. The land use data in 2010 came from the Resource and Environment Data Center of the Chinese Academy of Sciences (http://www.resdc.cn, accessed on: 1 October 2022), and the soil attribute data came from the Harmonized World Soil Database (HWSD), with a spatial resolution of 1000 m. The detailed model construction steps refer to Wu et al. [25].

**Figure 1.** Location of the study area and the stations. **Figure 1.** Location of the study area and the stations.

#### *2.2. Data*  **3. Methods**

In this study, the ground-measured meteorological data include the daily data of 6 meteorological stations and 16 rainfall stations in the study area. The data of the meteorological stations include daily rainfall, temperature, sunshine duration, relative humidity, air pressure and average wind speed. The measured daily runoff data from 2001 to 2017 come from the Beidao Hydrological Station. The location of the stations is shown in Figure 1. The measured flow data and meteorological data come from the *Yellow River Basin Hydrological Yearbook* and the China Meteorological Data Network, respectively. The vegetation data from 2001 to 2017 were derived from the NDVI dataset provided by the MODIS/Terra website. The dataset was the MOD13Q1 product with a spatial resolution of 250 m and a temporal resolution of 16 days. The DEM data used to build the hydrological model came from the ASTER GDEM data provided by Geospatial Data Cloud (http://www.gscloud.cn, accessed on: 1 October 2022), with a spatial resolution of 30 m. The land use data in 2010 came from the Resource and Environment Data Center of the Chinese Academy of Sciences (http://www.resdc.cn, accessed on: 1 October 2022), and the soil attribute data came from the Harmonized Since there is a lack of recorded solar radiation data at all meteorological stations in the study area, the solar radiation in the watershed was approximated by the number of sunshine hours. The measured meteorological data collected in this study are relatively comprehensive, and the Penman–Monteith method was chosen for the calculation of potential evapotranspiration. For the whole basin average precipitation, potential evapotranspiration, temperature, runoff, vegetation index and other variables in the study area, the linear regression method and the cumulative anomaly method were used to analyze the trend and abrupt change, respectively. The linear regression method is adopted to construct a univariate linear regression equation between the hydrometeorological variable sequence and the time series. The regression coefficient *b* can intuitively reflect the changing rate of the hydrometeorological sequence, and its positive and negative values indicate the increasing or decreasing trend of the hydrometeorological sequence, respectively. The *t*-test method (*a* = 0.05) is used to test the significance of the regression coefficient of the equation; the significance of the change trend of the hydrometeorological sequence is evaluated [26]. The formula of linear regression is as follows:

$$\mathbf{x}(t) = a + bt\tag{1}$$

where, *x* is the variable, *t* is time series, *a* and *b* are parameters.

**3. Methods** Since there is a lack of recorded solar radiation data at all meteorological stations in the study area, the solar radiation in the watershed was approximated by the number of sunshine hours. The measured meteorological data collected in this study are relatively comprehensive, and the Penman–Monteith method was chosen for the calculation of The cumulative anomaly method is a method for evaluating the trend of change according to the increase and decrease in the cumulative anomaly curve. The long-term evolution trend of the hydrometeorological sequence and the approximate time of the sudden change can be evaluated [27]. The formula to calculate the cumulative anomaly value *cd* at time *t* is:

$$cd = \sum\_{t=1}^{n} (\mathbf{x}(t) - \mathbf{x}\_{\text{mean}}) \tag{2}$$

where, *x*(*t*) is the variable at *t* and *xmean* is the mean of *x*.

The SWAT model is a physically based distributed hydrological model and is widely used to simulate the hydrological process of the watershed all over the world. It has great advantages in simulating long-term continuous hydrological processes in large-scale complex watersheds. Therefore, since the model was launched, it has been used in Europe, North America, Canada, Asia and other regions with good results [28,29]. The basin SWAT model was constructed, and the DEM data of the basin were used for catchment analysis and river network extraction. The minimum catchment area threshold was set to 400 km<sup>2</sup> , and the study area was divided into 39 sub-basins. Sensitivity analysis was performed using the SUFI-2 algorithm, in which six sensitive parameters were closely related to runoff, such as the number of runoff curves (CN2), the effective hydraulic conductivity of the main channel bed (CH\_K2), the soil evaporation compensation coefficient (ESCO), the baseflow *α* factor of riparian regulation and storage (ALPHA\_BNK), the Manning coefficient of the main channel (CH\_N2) and the saturated hydraulic conductivity of the soil (SOL\_K). Since CH\_K2 and SOL\_K are parameters that characterize the soil characteristics of the basin, they are not suitable to establish a relationship with the meteorological characteristics, and only CN2, ESCO, ALPHA\_BNK and CH\_N2 are analyzed.

The Nash–Sutcliffe efficiency coefficient (NSE), the coefficient of determination (*R* 2 ) and the Kling–Gupta coefficient (KGE) [30–33] are selected to evaluate the simulated daily flow. The formulas are as follows:

$$\text{NSE} = 1 - \frac{\sum\_{i=1}^{n} (Q\_{o,i} - Q\_{s,i})^2}{\sum\_{i=1}^{n} \left(Q\_{o,i} - \overline{Q}\_o\right)^2} \tag{3}$$

$$\mathcal{R}^2 = \frac{\sum\_{i=1}^{n} \left[ \left( Q\_{s,i} - \overline{Q}\_s \right) \left( Q\_{o,i} - \overline{Q}\_o \right) \right]^2}{\sum\_{i=1}^{n} \left( Q\_{o,i} - \overline{Q}\_o \right)^2 \sum\_{i=1}^{n} \left( Q\_{s,i} - \overline{Q}\_s \right)^2} \tag{4}$$

$$\text{KGE} = 1 - \sqrt{(r-1)^2 + (\beta - 1)^2 + (\gamma - 1)^2} \tag{5}$$

$$\beta = \frac{\mu\_{\rm s}}{\mu\_{\rm o}}, \gamma = \frac{\sigma\_{\rm s}/\mu\_{\rm s}}{\sigma\_{\rm o}/\mu\_{\rm o}} \tag{6}$$

where, *n* represents the length of the runoff sequence; *Q<sup>s</sup>* and *Q<sup>o</sup>* represent the simulated flow and measured flow; *µ* and *σ* are the mean and variance, respectively; *r* is the linear correlation coefficient between the simulated and measured values. The closer NSE and *R* 2 are to 1, the better the simulation effect is. The KGE coefficient is a comprehensive index including Cv, correlation coefficient and mean value. The closer the value is to 1, the higher the simulation accuracy.

## **4. Results**

#### *4.1. Spatial Distribution of Meteorological Factors*

According to the precipitation, air temperature and potential evapotranspiration data calculated by the Penman–Monteith formula from 1965 to 2017 at the six meteorological stations in the upper reach of the Wei River, the spatial distribution of the annual average precipitation, temperature and potential evapotranspiration was plotted by inverse distance weighting interpolation (IDW), as shown in Figure 2. It shows that the annual average precipitation in the basin has obvious differences in spatial distribution, and generally shows a decreasing trend from the southwest of the basin to the northeast. The spatial distribution of the average temperature in the basin shows an increasing trend from the northwest to the southeast. The temperature in about two-thirds of the area is between 6 and 11 ◦C. From the perspective of the river flow, the high temperature area is located at the outlet of the downstream watershed, and the temperature in the middle and upper reaches of the watershed is lower. It is closely related to the changes of elevation in the basin. The spatial distribution of the annual average potential evapotranspiration in the basin is extremely uneven, and the overall distribution pattern is that the potential evapotranspiration in the

northern part of the basin is slightly larger than that in the southern part of the basin. The spatial distributions of potential evapotranspiration and precipitation are basically opposite; the area with the largest potential evapotranspiration has less precipitation and, on the contrary, the area with the smallest potential evapotranspiration has more precipitation. that in the southern part of the basin. The spatial distributions of potential evapotranspiration and precipitation are basically opposite; the area with the largest potential evapotranspiration has less precipitation and, on the contrary, the area with the smallest potential evapotranspiration has more precipitation.

*Water* **2022**, *14*, x FOR PEER REVIEW 6 of 17

the higher the simulation accuracy.

*4.1. Spatial Distribution of Meteorological Factors*

**4. Results**

*R<sup>2</sup>* are to 1, the better the simulation effect is. The KGE coefficient is a comprehensive index including Cv, correlation coefficient and mean value. The closer the value is to 1,

According to the precipitation, air temperature and potential evapotranspiration data calculated by the Penman–Monteith formula from 1965 to 2017 at the six meteorological stations in the upper reach of the Wei River, the spatial distribution of the annual average precipitation, temperature and potential evapotranspiration was plotted by inverse distance weighting interpolation (IDW), as shown in Figure 2. It shows that the annual average precipitation in the basin has obvious differences in spatial distribution, and generally shows a decreasing trend from the southwest of the basin to the northeast. The spatial distribution of the average temperature in the basin shows an increasing trend from the northwest to the southeast. The temperature in about two-thirds of the area is between 6 and 11 °C. From the perspective of the river flow, the high temperature area is located at the outlet of the downstream watershed, and the temperature in the middle and upper reaches of the watershed is lower. It is closely related to the changes of elevation in the basin. The spatial distribution of the annual average potential evapotranspiration in the basin is extremely uneven, and the overall distribution pattern is that the potential evapotranspiration in the northern part of the basin is slightly larger than

**Figure 2.** Spatial distribution map of annual average precipitation, temperature and potential evapotranspiration. (**a**) precipitation, (**b**) temperature, (**c**) potential evapotranspiration. **Figure 2.** Spatial distribution map of annual average precipitation, temperature and potential evapotranspiration. (**a**) precipitation, (**b**) temperature, (**c**) potential evapotranspiration.

The linear regression method was used to analyze the trend of the annual precipita-

further tested. The results are shown in Figure 3. It shows that the annual precipitation in the basin shows an insignificant decreasing trend, with a decreasing rate of 0.59 mm/a. The mean annual precipitation in the basin was 491.37 mm, the annual precipitation reached the maximum in 1967, about 704.04 mm, and the minimum in 1982. The annual average temperature in the basin shows a significant growth trend, with an increasing rate of 0.03 °C/a. The mean annual temperature in the basin is 6.62 °C. The annual average temperature reached the maximum in 2016, about 7.98 °C, and was the minimum in 1967. The annual potential evapotranspiration in the basin showed a significant increasing trend, with an increasing rate of 1.87 mm/a. The mean annual potential evapotranspiration in the basin was 1150.72 mm, and the annual potential evapotranspiration reached the maximum in 2016, about 1280.31 mm, and was the minimum in 1967. The annual runoff depth of the Beidao Hydrological Station showed an insignificant trend with an increasing rate of 0.42 mm/a. The mean annual runoff depth of the Beidao Hy-

*4.2. Interannual Changes of the Environmental Factors*

#### *4.2. Interannual Changes of the Environmental Factors*

The linear regression method was used to analyze the trend of the annual precipitation, temperature, potential evapotranspiration and runoff depth series at the whole basin scale in the study area from 1965 to 2017, and the significance of the trend results was further tested. The results are shown in Figure 3. It shows that the annual precipitation in the basin shows an insignificant decreasing trend, with a decreasing rate of 0.59 mm/a. The mean annual precipitation in the basin was 491.37 mm, the annual precipitation reached the maximum in 1967, about 704.04 mm, and the minimum in 1982. The annual average temperature in the basin shows a significant growth trend, with an increasing rate of 0.03 ◦C/a. The mean annual temperature in the basin is 6.62 ◦C. The annual average temperature reached the maximum in 2016, about 7.98 ◦C, and was the minimum in 1967. The annual potential evapotranspiration in the basin showed a significant increasing trend, with an increasing rate of 1.87 mm/a. The mean annual potential evapotranspiration in the basin was 1150.72 mm, and the annual potential evapotranspiration reached the maximum in 2016, about 1280.31 mm, and was the minimum in 1967. The annual runoff depth of the Beidao Hydrological Station showed an insignificant trend with an increasing rate of 0.42 mm/a. The mean annual runoff depth of the Beidao Hydrological Station was 26.20 mm, and the annual runoff depth reached the maximum in 2013, about 49.01 mm, and was the minimum in 2002. *Water* **2022**, *14*, x FOR PEER REVIEW 8 of 17 drological Station was 26.20 mm, and the annual runoff depth reached the maximum in 2013, about 49.01 mm, and was the minimum in 2002.

*y* = 1.8727*x* + 1100.16 R² = 0.2599 **Figure 3.** *Cont*.

1960 1980 2000 2020

years

1000

1100

1200

annual potential

evapotranspiration/mm

1300





cumulative anomaly


1960 1970 1980 1990 2000 2010 2020

1994

years

(**e**) (**f**)

annual temperature/℃

4

5

6

7

8

9

200

1960 1980 2000 2020

*y* = −0.5865*x* + 507.21 R² = 0.0102

years

*y* = 0.031*x* + 5.78 R² = 0.5657

years

300

400

500

annual precipitation/mm

600

700

800

**Figure 3.** Analysis of trends and abrupt changes of hydrological and meteorological elements in the basin. Note: *x* = year−1964. (**a**) Annual precipitation trend, (**b**) Analysis of abrupt change in annual precipitation, (**c**) Annual temperature trend, (**d**) Analysis of abrupt change in annual temperature, (**e**) Annual potential evapotranspiration trend, (**f**) Analysis of abrupt change in annual potential evapotranspiration, (**g**) Annual runoff depth trend, (**h**) Analysis of abrupt change in annual runoff depth. **Figure 3.** Analysis of trends and abrupt changes of hydrological and meteorological elements in the basin. Note: *x* = year−1964. (**a**) Annual precipitation trend, (**b**) Analysis of abrupt change in annual precipitation, (**c**) Annual temperature trend, (**d**) Analysis of abrupt change in annual temperature, (**e**) Annual potential evapotranspiration trend, (**f**) Analysis of abrupt change in annual potential evapotranspiration, (**g**) Annual runoff depth trend, (**h**) Analysis of abrupt change in annual runoff depth.

drological Station was 26.20 mm, and the annual runoff depth reached the maximum in

2013, about 49.01 mm, and was the minimum in 2002.

(**a**) (**b**)




cumulative anomaly



1960 1980 2000 2020

years

1960 1980 2000 2020

years

1997

100

300

cumulative anomaly

500

700

(**c**) (**d**)

The cumulative anomaly method was used to analyze the mutation of the annual precipitation, temperature, potential evapotranspiration and runoff depth series in the study area from 1965 to 2017. The results are shown in Figure 3. It was identified that the annual precipitation and runoff depth series did not have abrupt changes, the annual mean temperature series had a significant jump around 1997, and the annual potential evapotranspiration series had a significant jump around 1994. The simulation period is 2001–2017 and after the abrupt changes. The cumulative anomaly method was used to analyze the mutation of the annual precipitation, temperature, potential evapotranspiration and runoff depth series in the study area from 1965 to 2017. The results are shown in Figure 3. It was identified that the annual precipitation and runoff depth series did not have abrupt changes, the annual mean temperature series had a significant jump around 1997, and the annual potential evapotranspiration series had a significant jump around 1994. The simulation period is 2001–2017 and after the abrupt changes.

Figure 4 shows the changing process of the annual average NDVI series in the study area from 2001 to 2017. It shows that although the average annual NDVI value of the watershed fluctuated during the 17 years, the overall trend is significant increase, which indicates that the vegetation coverage in the study area during the 17 years was increasing continuously. Especially after 2011, the vegetation condition significantly improved. The annual NDVI value of the watershed was between 0.30 and 0.38 from 2001 to 2017. The average annual NDVI value of the watershed increased by 19.7%, and reached the maximum value of about 0.38 in 2013. Figure 5 shows the seasonality of NDVI on the monthly scale in the study area from 2001 to 2017. It shows that the vegetation change in the study area has a significant seasonality across the year. The maximum monthly NDVI in the study area generally occurs from June to August, while the minimum monthly Figure 4 shows the changing process of the annual average NDVI series in the study area from 2001 to 2017. It shows that although the average annual NDVI value of the watershed fluctuated during the 17 years, the overall trend is significant increase, which indicates that the vegetation coverage in the study area during the 17 years was increasing continuously. Especially after 2011, the vegetation condition significantly improved. The annual NDVI value of the watershed was between 0.30 and 0.38 from 2001 to 2017. The average annual NDVI value of the watershed increased by 19.7%, and reached the maximum value of about 0.38 in 2013. Figure 5 shows the seasonality of NDVI on the monthly scale in the study area from 2001 to 2017. It shows that the vegetation change in the study area has a significant seasonality across the year. The maximum monthly NDVI in the study area

NDVI generally occurs from January to March.

generally occurs from June to August, while the minimum monthly NDVI generally occurs from January to March. *Water* **2022**, *14*, x FOR PEER REVIEW 10 of 17

*Water* **2022**, *14*, x FOR PEER REVIEW 10 of 17

*y* = 0.0034*x* + 0.31

0.4

**Figure 4.** Annual average NDVI in the basin. Note: *x* = year−2000. **Figure 4.** Annual average NDVI in the basin. Note: *x* = year−2000.

month **Figure 5.** Seasonality of Monthly NDVI. **Figure 5.** Seasonality of Monthly NDVI.

#### **Figure 5.** Seasonality of Monthly NDVI. *4.3. Simulation of the Hydrological Process 4.3. Simulation of the Hydrological Process*

*4.3. Simulation of the Hydrological Process* According to the mutation analysis results of the hydrometeorological data, 2001– 2017 was selected as the simulation period of the model, 2001 was set as the warm-up period of the model, 2002–2013 was the calibration period, and 2014–2017 was the verification period. Sensitivity analysis was carried out on the model parameters, among which there are six parameters closely related to runoff, namely the number of runoff curves (CN2), the effective hydraulic conductivity of the main channel bed (CH\_K2), the According to the mutation analysis results of the hydrometeorological data, 2001– 2017 was selected as the simulation period of the model, 2001 was set as the warm-up period of the model, 2002–2013 was the calibration period, and 2014–2017 was the verification period. Sensitivity analysis was carried out on the model parameters, among which there are six parameters closely related to runoff, namely the number of runoff curves (CN2), the effective hydraulic conductivity of the main channel bed (CH\_K2), the soil evaporation compensation coefficient (ESCO), and the baseflow *α* factor of riparian According to the mutation analysis results of the hydrometeorological data, 2001–2017 was selected as the simulation period of the model, 2001 was set as the warm-up period of the model, 2002–2013 was the calibration period, and 2014–2017 was the verification period. Sensitivity analysis was carried out on the model parameters, among which there are six parameters closely related to runoff, namely the number of runoff curves (CN2), the effective hydraulic conductivity of the main channel bed (CH\_K2), the soil evaporation compensation coefficient (ESCO), and the baseflow *α* factor of riparian regulation and

soil evaporation compensation coefficient (ESCO), and the baseflow *α* factor of riparian

storage (ALPHA\_BNK), the Manning coefficient of the main channel (CH\_N2) and the saturated hydraulic conductivity of the soil (SOL\_K). (CH\_N2) and the saturated hydraulic conductivity of the soil (SOL\_K). The daily runoff at the Beidao Hydrological Station is simulated, and the measured and simulated flow processes in the calibration and verification periods are shown in

regulation and storage (ALPHA\_BNK), the Manning coefficient of the main channel

The daily runoff at the Beidao Hydrological Station is simulated, and the measured and simulated flow processes in the calibration and verification periods are shown in Figure 6, and the evaluation results of the daily flow simulation are shown in Table 1. It shows that the simulated daily flow process in the calibration period of the Beidao hydrological station is quite good, and the three evaluation indicators are all above 0.54, of which the KGE index reaches 0.70, and the simulation effect is poor in the verification period except for 2014. In general, the SWAT model has been successfully constructed in the study area, and has good applicability in the watershed. Further research is carried out based on this model. Figure 6, and the evaluation results of the daily flow simulation are shown in Table 1. It shows that the simulated daily flow process in the calibration period of the Beidao hydrological station is quite good, and the three evaluation indicators are all above 0.54, of which the KGE index reaches 0.70, and the simulation effect is poor in the verification period except for 2014. In general, the SWAT model has been successfully constructed in the study area, and has good applicability in the watershed. Further research is carried out based on this model.

*Water* **2022**, *14*, x FOR PEER REVIEW 11 of 17

**Figure 6.** Results of simulated and observed streamflow in calibration period and validation period. (**a**) calibration period. (**b**) verification period. **Figure 6.** Results of simulated and observed streamflow in calibration period and validation period. (**a**) calibration period. (**b**) verification period.

**Table 1.** Evaluation result of daily flow simulation value. **Table 1.** Evaluation result of daily flow simulation value.


#### *4.4. Time-varying Parameters of the Hydrological Model 4.4. Time-Varying Parameters of the Hydrological Model*

2009, 2016 and 2017.

In order to obtain the time-varying parameters of the SWAT model, the SUFI-2 algorithm was used to calibrate the model parameters in each year from 2001 to 2017, and the previous year of each calibration period was used as the model warm-up period. Then, the daily runoff process of the study area in each year was simulated by the SWAT model with one parameter set of this year. The Nash–Sutcliffe efficiency coefficient, *R<sup>2</sup>* and KGE were used to evaluate the simulation accuracy of the runoff in each period. The results are shown in Table 2. It shows that the daily runoff simulation accuracy in the other 10 years met the requirements, except for the poor simulation in 2002, 2004, 2008, In order to obtain the time-varying parameters of the SWAT model, the SUFI-2 algorithm was used to calibrate the model parameters in each year from 2001 to 2017, and the previous year of each calibration period was used as the model warm-up period. Then, the daily runoff process of the study area in each year was simulated by the SWAT model with one parameter set of this year. The Nash–Sutcliffe efficiency coefficient, *R* <sup>2</sup> and KGE were used to evaluate the simulation accuracy of the runoff in each period. The results are shown in Table 2. It shows that the daily runoff simulation accuracy in the other 10 years met the requirements, except for the poor simulation in 2002, 2004, 2008, 2009, 2016 and 2017.

**Year R<sup>2</sup> NSE KGE** 2002 0.33 0.26 0.53

**Table 2.** Evaluation result of simulated daily flow of each year.


**Table 2.** Evaluation result of simulated daily flow of each year. *Water* **2022**, *14*, x FOR PEER REVIEW 12 of 17

The hydrological model was calibrated year by year in 1-year increments. Four timedependent model parameter series, CN2, ESCO, ALPHA\_BNK and CH\_N2, were obtained for the study area. The relationship between model parameters and environmental characterization factors, i.e., precipitation, potential evapotranspiration and NDVI in the basin is plotted, as shown in Figures 7–10. It shows that the four parameters are fluctuating in different degrees during the study period. There was no obvious trend in CN2 and CH\_N2. The soil evaporation compensation coefficient (ESCO) and ALPHA\_BNK showed a downward trend, but the magnitude of change in ALPHA\_BNK was small and negligible. The smaller the value, the more evapotranspiration water can be obtained from the lower soil layer, which is consistent with the conclusion that the watershed evaporation is increasing year by year, as shown in Figure 8. The hydrological model was calibrated year by year in 1-year increments. Four time-dependent model parameter series, CN2, ESCO, ALPHA\_BNK and CH\_N2, were obtained for the study area. The relationship between model parameters and environmental characterization factors, i.e., precipitation, potential evapotranspiration and NDVI in the basin is plotted, as shown in Figures 7–10. It shows that the four parameters are fluctuating in different degrees during the study period. There was no obvious trend in CN2 and CH\_N2. The soil evaporation compensation coefficient (ESCO) and AL-PHA\_BNK showed a downward trend, but the magnitude of change in ALPHA\_BNK was small and negligible. The smaller the value, the more evapotranspiration water can be obtained from the lower soil layer, which is consistent with the conclusion that the watershed evaporation is increasing year by year, as shown in Figure 8. The hydrological model was calibrated year by year in 1-year increments. Four time-dependent model parameter series, CN2, ESCO, ALPHA\_BNK and CH\_N2, were obtained for the study area. The relationship between model parameters and environmental characterization factors, i.e., precipitation, potential evapotranspiration and NDVI in the basin is plotted, as shown in Figures 7–10. It shows that the four parameters are fluctuating in different degrees during the study period. There was no obvious trend in CN2 and CH\_N2. The soil evaporation compensation coefficient (ESCO) and AL-PHA\_BNK showed a downward trend, but the magnitude of change in ALPHA\_BNK was small and negligible. The smaller the value, the more evapotranspiration water can be obtained from the lower soil layer, which is consistent with the conclusion that the watershed evaporation is increasing year by year, as shown in Figure 8.

2017 0.29 0.12 0.53

2017 0.29 0.12 0.53

**Figure 7.** Time-varying process of CN2 to precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 7.** Time-varying process of CN2 to precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 7.** Time-varying process of CN2 to precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

**Figure 8.** Time-varying process of ESCO to basin precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

*Water* **2022**, *14*, x FOR PEER REVIEW 13 of 17

**Figure 9.** Time-varying process of CH\_N2 to basin precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 9.** Time-varying process of CH\_N2 to basin precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 9.** Time-varying process of CH\_N2 to basin precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

**Figure 8.** Time-varying process of ESCO to basin precipitation, potential evapotranspiration and

**Figure 8.** Time-varying process of ESCO to basin precipitation, potential evapotranspiration and

tion and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. *4.5. Analysis of the Relationship between Model Parameters and Environmental Factors*  **Figure 10.** Time-varying process of ALPHA\_BNK to basin precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 10.** Time-varying process of ALPHA\_BNK to basin precipitation, potential evapotranspiration and NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

#### The scatter plots of the hydrological model parameters and environmental repre-*4.5. Analysis of the Relationship between Model Parameters and Environmental Factors 4.5. Analysis of the Relationship between Model Parameters and Environmental Factors*

sentation factors of precipitation, potential evapotranspiration and NDVI were shown to identify the relationship between them, and the results are shown in Table 3, Figures 11– 14. Figure 11 shows that there is a weak negative correlation between CN2 and potential evapotranspiration. The number of runoff curves (CN2) changes inversely with the change of potential evapotranspiration, and has no obvious correlation with precipitation and NDVI. In Figure 12, ESCO showed a certain negative correlation with potential evapotranspiration and NDVI, and the correlation coefficient passed the 0.05 significance level test. The soil evaporation compensation coefficient (ESCO) decreased with the increase in potential evapotranspiration and NDVI, and had no significant correlation with precipitation. The smaller the value of ESCO, the more evaporative water the model can obtain from the lower soil layer, resulting in the reduction of surface runoff, interflow and subsurface runoff, especially the most obvious reduction in surface runoff. In Figure 13, CH\_N2 has a certain correlation with precipitation and potential evapotranspiration. It increases with the increase in precipitation and decreases with the increase in potential evapotranspiration. There is no obvious correlation with NDVI. The correlation coefficient between CH\_N2 and precipitation passes the 0.05 significance level test. In Figure 14, ALPHA\_BNK showed a certain negative correlation with potential evapotranspiration and NDVI, and the correlation coefficient passed the 0.05 significance level test. The baseflow *α* factor of riparian regulation and storage (ALPHA\_BNK) decreased with the increase in potential evapotranspiration and NDVI, and had a weak positive correlation with precipitation. The baseflow *α* factor of riparian regulation and storage is a response The scatter plots of the hydrological model parameters and environmental representation factors of precipitation, potential evapotranspiration and NDVI were shown to identify the relationship between them, and the results are shown in Table 3, Figures 11– 14. Figure 11 shows that there is a weak negative correlation between CN2 and potential evapotranspiration. The number of runoff curves (CN2) changes inversely with the change of potential evapotranspiration, and has no obvious correlation with precipitation and NDVI. In Figure 12, ESCO showed a certain negative correlation with potential evapotranspiration and NDVI, and the correlation coefficient passed the 0.05 significance level test. The soil evaporation compensation coefficient (ESCO) decreased with the increase in potential evapotranspiration and NDVI, and had no significant correlation with precipitation. The smaller the value of ESCO, the more evaporative water the model can obtain from the lower soil layer, resulting in the reduction of surface runoff, interflow and subsurface runoff, especially the most obvious reduction in surface runoff. In Figure 13, CH\_N2 has a certain correlation with precipitation and potential evapotranspiration. It increases with the increase in precipitation and decreases with the increase in potential evapotranspiration. There is no obvious correlation with NDVI. The correlation coefficient between CH\_N2 and precipitation passes the 0.05 significance level test. In Figure 14, ALPHA\_BNK showed a certain negative correlation with potential evapotranspiration and NDVI, and the correlation coefficient passed the 0.05 significance level test. The baseflow *α* factor of riparian regulation and storage (ALPHA\_BNK) decreased with the increase in potential evapotranspiration and NDVI, and had a weak positive correlation with precipitation. The baseflow *α* factor of riparian regulation and storage is a response The scatter plots of the hydrological model parameters and environmental representation factors of precipitation, potential evapotranspiration and NDVI were shown to identify the relationship between them, and the results are shown in Table 3, Figures 11–14. Figure 11 shows that there is a weak negative correlation between CN2 and potential evapotranspiration. The number of runoff curves (CN2) changes inversely with the change of potential evapotranspiration, and has no obvious correlation with precipitation and NDVI. In Figure 12, ESCO showed a certain negative correlation with potential evapotranspiration and NDVI, and the correlation coefficient passed the 0.05 significance level test. The soil evaporation compensation coefficient (ESCO) decreased with the increase in potential evapotranspiration and NDVI, and had no significant correlation with precipitation. The smaller the value of ESCO, the more evaporative water the model can obtain from the lower soil layer, resulting in the reduction of surface runoff, interflow and subsurface runoff, especially the most obvious reduction in surface runoff. In Figure 13, CH\_N2 has a certain correlation with precipitation and potential evapotranspiration. It increases with the increase in precipitation and decreases with the increase in potential evapotranspiration. There is no obvious correlation with NDVI. The correlation coefficient between CH\_N2 and precipitation passes the 0.05 significance level test. In Figure 14, ALPHA\_BNK showed a certain negative correlation with potential evapotranspiration and NDVI, and the correlation coefficient passed the 0.05 significance level test. The baseflow *α* factor of riparian regulation and storage (ALPHA\_BNK) decreased with the increase in potential evapotranspiration and NDVI, and had a weak positive correlation with precipitation. The baseflow *α* factor of riparian regulation and storage is a response index reflecting the discharge rate of subsurface runoff to the river discharge at the river bank, and the river discharge is positively correlated with it.


**Table 3.** The *R* <sup>2</sup> of each parameter with precipitation, potential evapotranspiration and NDVI. **Table 3.** The R<sup>2</sup> of each parameter with precipitation, potential evapotranspiration and NDVI. **R<sup>2</sup> Precipitation Potential Evapotranspiration NDVI Table 3.** The R<sup>2</sup> of each parameter with precipitation, potential evapotranspiration and NDVI. **R<sup>2</sup> Precipitation Potential Evapotranspiration NDVI**

index reflecting the discharge rate of subsurface runoff to the river discharge at the river

index reflecting the discharge rate of subsurface runoff to the river discharge at the river

Note: \* Indicates passing the 0.05 significance level test. Note :\* Indicates passing the 0.05 significance level test. Note :\* Indicates passing the 0.05 significance level test.

*Water* **2022**, *14*, x FOR PEER REVIEW 14 of 17

*Water* **2022**, *14*, x FOR PEER REVIEW 14 of 17

bank, and the river discharge is positively correlated with it.

bank, and the river discharge is positively correlated with it.

**Figure 11.** Correlation of CN2 to watershed precipitation, potential evapotranspiration, NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 11.** Correlation of CN2 to watershed precipitation, potential evapotranspiration, NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 11.** Correlation of CN2 to watershed precipitation, potential evapotranspiration, NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

**Figure 13.** Correlation of CH\_N2 to watershed precipitation, potential evapotranspiration, NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 13.** Correlation of CH\_N2 to watershed precipitation, potential evapotranspiration, NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

*y* = −8.6677*x* + 0.3521 R² = 0.2563

0 0.002 0.004 0.006

ALPHA\_BNK

**Figure 14.** Correlation of ALPHA\_BNK to watershed precipitation, potential evapotranspiration,

In summary, all four parameters of the SWAT model in the upper reach of the Wei River Basin show a certain correlation with potential evapotranspiration, and some parameters have a correlation with precipitation or NDVI, which indicates that the changes of the hydrological model parameters in the upper reach of the Wei River Basin are closely related to the dynamic changes of potential evapotranspiration. Therefore, potential evapotranspiration can be used as a representative factor of the dynamic change of

0.26

0.31

0.36

NDVI

0.41

In this study, the upper reach of the Wei River is selected as the study area, and the temporal and spatial evolution of the ecological hydrometeorological elements in the study area is revealed using the linear regression method and the cumulative anomaly method. The time-varying parameter series of the hydrological model is obtained based on the SWAT model in each year. The time-varying characteristics of the model parameters and environmental indicators such as precipitation, potential evapotranspiration and NDVI were analyzed, and the response relationship between them was identified. In addition to the decreasing trend of annual precipitation in the basin, the hydrometeorological elements such as annual potential evapotranspiration, annual average temperature, annual runoff depth and annual average NDVI all showed an increasing trend. The annual precipitation and annual runoff depth series in the basin did not change abruptly, but the annual mean temperature series and the annual potential evapotranspiration series jumped significantly around 1997 and 1994, respectively.

Except for the six years of 2002, 2004, 2008, 2009, 2016 and 2017 in 2002–2017, the daily runoff simulation accuracy in the other 10 years met the requirements. The four model parameters of CN2, ESCO, CH\_N2 and ALPHA\_BNK all fluctuated with time during the study period. The number of runoff curves (CN2) and CH\_N2 had no obvious

potential

evapotranspiration/mm

*y* = −27222*x* + 1232.6 R² = 0.5453

0 0.002 0.004 0.006

ALPHA\_BNK

(**a**) (**b**) (**c**)

NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

model parameters.

**5. Conclusions**

*y* = 24720*x* + 460.54 R² = 0.1036

0 0.002 0.004 0.006

ALPHA\_BNK

precipitation/mm

precipitation/mm

*y* = 1030.8*x* + 264.16 R² = 0.3364

0.1 0.2 0.3

CH\_N2

(**a**) (**b**) (**c**)

*y* = −387.29*x* + 1281.3 R² = 0.2062

0.1 0.2 0.3

CH\_N2

(**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

potential

evapotranspiration/mm

**Figure 14.** Correlation of ALPHA\_BNK to watershed precipitation, potential evapotranspiration, NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI. **Figure 14.** Correlation of ALPHA\_BNK to watershed precipitation, potential evapotranspiration, NDVI. (**a**) Precipitation, (**b**) potential evapotranspiration, (**c**) NDVI.

**Figure 13.** Correlation of CH\_N2 to watershed precipitation, potential evapotranspiration, NDVI.

*y* = −0.0119*x* + 0.3427

0.1 0.2 0.3

CH\_N2

R² = 0.0009 0.26

0.31

0.36

NDVI

0.41

In summary, all four parameters of the SWAT model in the upper reach of the Wei River Basin show a certain correlation with potential evapotranspiration, and some parameters have a correlation with precipitation or NDVI, which indicates that the changes of the hydrological model parameters in the upper reach of the Wei River Basin are closely related to the dynamic changes of potential evapotranspiration. Therefore, potential evapotranspiration can be used as a representative factor of the dynamic change of In summary, all four parameters of the SWAT model in the upper reach of the Wei River Basin show a certain correlation with potential evapotranspiration, and some parameters have a correlation with precipitation or NDVI, which indicates that the changes of the hydrological model parameters in the upper reach of the Wei River Basin are closely related to the dynamic changes of potential evapotranspiration. Therefore, potential evapotranspiration can be used as a representative factor of the dynamic change of model parameters.

#### model parameters. **5. Conclusions**

**5. Conclusions** In this study, the upper reach of the Wei River is selected as the study area, and the temporal and spatial evolution of the ecological hydrometeorological elements in the study area is revealed using the linear regression method and the cumulative anomaly method. The time-varying parameter series of the hydrological model is obtained based on the SWAT model in each year. The time-varying characteristics of the model parameters and environmental indicators such as precipitation, potential evapotranspiration In this study, the upper reach of the Wei River is selected as the study area, and the temporal and spatial evolution of the ecological hydrometeorological elements in the study area is revealed using the linear regression method and the cumulative anomaly method. The time-varying parameter series of the hydrological model is obtained based on the SWAT model in each year. The time-varying characteristics of the model parameters and environmental indicators such as precipitation, potential evapotranspiration and NDVI were analyzed, and the response relationship between them was identified.

and NDVI were analyzed, and the response relationship between them was identified. In addition to the decreasing trend of annual precipitation in the basin, the hydrometeorological elements such as annual potential evapotranspiration, annual average temperature, annual runoff depth and annual average NDVI all showed an increasing trend. The annual precipitation and annual runoff depth series in the basin did not change abruptly, but the annual mean temperature series and the annual potential In addition to the decreasing trend of annual precipitation in the basin, the hydrometeorological elements such as annual potential evapotranspiration, annual average temperature, annual runoff depth and annual average NDVI all showed an increasing trend. The annual precipitation and annual runoff depth series in the basin did not change abruptly, but the annual mean temperature series and the annual potential evapotranspiration series jumped significantly around 1997 and 1994, respectively.

evapotranspiration series jumped significantly around 1997 and 1994, respectively. Except for the six years of 2002, 2004, 2008, 2009, 2016 and 2017 in 2002–2017, the daily runoff simulation accuracy in the other 10 years met the requirements. The four model parameters of CN2, ESCO, CH\_N2 and ALPHA\_BNK all fluctuated with time during the study period. The number of runoff curves (CN2) and CH\_N2 had no obvious Except for the six years of 2002, 2004, 2008, 2009, 2016 and 2017 in 2002–2017, the daily runoff simulation accuracy in the other 10 years met the requirements. The four model parameters of CN2, ESCO, CH\_N2 and ALPHA\_BNK all fluctuated with time during the study period. The number of runoff curves (CN2) and CH\_N2 had no obvious trend changes, while ESCO and ALPHA\_BNK showed a downward trend. Because the magnitude of change in ALPHA\_BNK was small, it can be neglected, and the change of ESCO is consistent with the conclusion that the evaporation in the basin increases year by year.

All model parameters show a certain correlation with potential evapotranspiration in the basin, which indicates that the changes of the hydrological model parameters in the upper reach of the Wei River Basin are closely related to the dynamic changes of potential evapotranspiration. Therefore, potential evapotranspiration can be used as an environmental factor to characterize the dynamic changes of the model parameters.

**Author Contributions:** Conceptualization, H.W., D.L., M.H. and R.L.; Methodology, H.W., D.L. and M.H.; Software, H.W. and D.L.; Validation, H.W. and D.L.; Formal analysis, H.W., D.L. and M.H.; Investigation, H.W. and D.L.; Resources, H.W. and D.L; Data Curation, H.W. and D.L.; Writing—Original Draft, H.W., D.L., M.H. and R.L.; Writing—Review and Editing, H.W., D.L., M.H., R.L., Q.Y., G.M. and H.L.; Visualization, H.W. and D.L.; Supervision, D.L. and M.H.; project administration, D.L. and M.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is supported by the National Natural Science Foundation of China (Grant Nos. 42071335, 52109031 and 41874017), the science and technology fund of the Second Monitoring and Application Center, China Earthquake Administration (KJ20220212) and the Natural Science Basic Research Program of Shaanxi-Han to Wei Water Transfer Joint Foundation (2022JC-LHJJ-13).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-7114-0