Multivariate Validation at Multistation of Distributed Watershed Hydrological Modeling Based on Multisource Data on Chinese Loess Plateau

Abstract

Earlier, hydrological simulation calibration and validation relied on flow observations at hydrological stations, but multisource observations changed the basin hydrological simulation from single-flow validation to multivariate validation, including evaporation, soil water, and runoff. This study used the Soil and Water Assessment Tool (SWAT) distributed hydrological model to simulate and investigate hydrological processes in the Jinghe River Basin in China. After a single-station, single-variable calibration using flow observation data at the Zhangjiashan Hydrological Station, multisource data were used to validate actual evaporation, soil water, and runoff. Using the flow station data from Zhangjiashan station for parameter calibration and validation, the simulated values of R², NSE, and KGE were all above 0.64, the PBIAS was within 20%, and the values of all the metrics in the calibration period were better than those in the validation period. The results show that the model performed satisfactorily, proving its regional applicability. Qingyang, Yangjiaping, and Zhangjiazhan stations had R², NSE, and KGE values above 0.57 and PBIAS within 25% during regional calibration, considering spatial variability. Additionally, simulation accuracy downstream increased. R², NSE, and KGE were above 0.50, and PBIAS was within 25% throughout validation, except for Qingyang, where the validation period was better than the calibration period. The Zhangjiashan station monthly runoff simulation improved after regional calibration. Runoff validation performed highest in the multivariate validation of evaporation–soil water–runoff, followed by actual evaporation and soil water content in China. The evaluation results for each hydrological variable improved after additional manual calibration. Multivariate verification based on multisource data improved the hydrological simulation at the basin scale.

1. Introduction

The calibration and verification of previous hydrological models relied heavily on the flow observation data provided by hydrological stations. As a typical distributed watershed hydrological model, the SWAT model was formulated by the United States Department of Agriculture (USDA). By integrating spatial data derived from geographic information systems (GIS) and remote sensing (RS), it is capable of simulating surface runoff, evapotranspiration, and water quality processes, and analyzing the impact of agricultural management methods and land use change on regional water supply [1]. It has been applied all over the world. To explore the evolution and development mechanism of runoff in the middle reaches of the Yellow River, Zhang and Zhi [2] simulated and discussed the runoff change process of the Yanhe River Basin under the condition of greening of the Chinese Loess Plateau by constructing SWAT-MODFLOW surface water–groundwater coupling model. Meanwhile, Peng et al. [3] conducted a case study in the Jinghe River Basin located on the Loess Plateau, utilizing the distributed eco-hydrological model RHESSys to simulate ecological and hydrological processes. The findings demonstrated that the RHESSys model is suitable for conducting eco-hydrological simulations in the large watershed of the Loess Plateau, and the model accurately captured the impact of vegetation changes on runoff, evapotranspiration, and soil water dynamics within the basin. Usually, the observed discharge is the main data to calibrate and verify the parameters of the hydrological models.

Nevertheless, there is still room for improvement in the SWAT model’s simulation of the flood peak process during flood season. The SWAT model is highly applicable in the upper sections of the Minjiang River of China, known for its complex topography and significant altitude variations. The observed and modeled daily runoff patterns exhibit overall consistency, particularly during the main flood period in June and July [4]. Hao et al. [5] constructed a land–atmosphere coupling simulation platform for the Weihe River Basin on the Chinese Loess Plateau to simulate the rainfall and runoff process in one month. They compared it with the measured rainfall and runoff daily process. Zhao et al. [6] employed three distinct approaches to quantify the sensitivity and uncertainty of SWAT model parameters in the Jingchuan River Basin on the Loess Plateau of China. For the model prediction uncertainty in the same parameter range, SUFI-2 (Sequential Uncertainty Fitting) is slightly better than GLUE (Generalized Likelihood Uncertainty Estimation). Based on the above considerations, how to more comprehensively calibrate the model parameters and verify the model effect is still an essential issue for hydrological simulation.

Global hydrometeorology data products have evolved from direct observation products to fusion products that combine data from multisource in recent years. The use of multisource observation data has advanced the hydrological simulation of the basin, allowing for the verification of evaporation, soil water, and runoff using multiple variables instead of only a single flow. Szeles et al. [7] initially calibrated the HBV model in a small watershed using only flow data. Subsequently, they progressively incorporated additional variables such as precipitation, snowfall, snow thickness, evapotranspiration observations, soil moisture, slope runoff, and groundwater level. The assessment findings at the yearly, seasonal, and daily levels indicated an enhancement in the overall consistency of the process. Chao et al. [8] thoroughly assessed five evapotranspiration data packages using ground observations and GRACE (Gravity Recovery and Climate Experiment) satellite observation data. They validated the hydrological simulation by calibrating it with evaporation. Lu et al. [9] utilized the SWAT distributed hydrological model to create a series of novel calibration techniques for suitability models. These methods considered various watershed characteristics, multiple time scales, multivariate data, and multiple stations. The application of these methods yielded favorable outcomes. Hydrological simulations tend to calibrate multivariate hydrological systems simultaneously utilizing multisource observational data.

On the Chinese Loess Plateau, considerable progress has been made in the observation and analysis of precipitation, potential evaporation, actual evaporation, soil moisture, and runoff [10,11,12,13,14,15]. Evaporation plays a crucial role in hydrological cycles. Wu et al. [16] employed the Penman–Monteith approach and they utilized potential evapotranspiration data (PET601) from the E-601 pan to operate the SWAT model in the upper reach of a hydrological station in the Weihe River’s northern region. The findings indicate that the PET values derived using various approaches exhibit significant disparities; however, they share comparable patterns yearly. The PET calculated using the PET601 approach was less than that calculated using the Penman–Monteith method. Xie et al. [17] utilized GRACE data with observed hydrological and meteorological data, including precipitation and runoff, to estimate the monthly evapotranspiration time series in the source region of the Yellow River. This estimation was achieved by using the water balance equation of the basin. The findings indicate that the predicted results align with the evapotranspiration simulation results of the Noah model of GLDAS-2.0 (Global Land Data Assimilation System). The effective storage of soil water is determined by soil infiltration, which is strongly associated with vegetation, atmospheric water, surface water, and groundwater [18,19]. The drawbacks of conventional techniques can be overcome by using RS technology, which has the advantages of speed, efficiency, continuity, and comprehensive coverage inverting soil moisture. There is considerable practical utility in conducting research on soil water RS inversion techniques [20]. Lu and Shi [21] examined the patterns and trends in the spatial and temporal distribution of surface soil water in China throughout the first part of the twenty-first century using AMSR-E surface soil water inversion products. It was found that the ensemble analysis data could overcome the systematic bias in a single inversion product to a certain extent. To compensate for the limitations of a single inversion product, they obtained spatial and temporal distribution information for remote sensing surface soil water in China. Huang et al. [22] examined the application of several soil moisture products (ECMWF-ERA5, GLDAS-Noah v2.1, and GLDAS-CLSM v2.2) in typical northern regions of China. The findings showed that the three model products’ soil moisture content was higher in both summer and autumn and that their geographical distribution was comparable.

Runoff describes complex hydrological processes due to several variables, including climate and the underlying surface. The impact of a changing environment on the water cycle is readily observed through alterations in runoff [23]. According to Ni and Lü [24], the Loess Plateau exhibits a high sensitivity of runoff changes to alterations in the underlying surface resulting from human activities. Conversely, increases in PET have the most negligible impact on runoff sensitivity in this region.

This study aimed to enhance hydrological modeling and utilize multisource observational data. It employed multisource datasets for verification purposes and used the distributed hydrological model SWAT to conduct a multivariate verification study on evaporation, soil water, and runoff in the Jinghe River Basin of the Chinese Loess Plateau. As a typical watershed on the Loess Plateau in the middle reaches of the Yellow River, the study of hydrological processes in the Jinghe River Basin is not only of great significance for understanding the hydrological cycle mechanism of the entire Loess Plateau but also provides a scientific basis for the planning and implementation of future work on soil and water conservation, so as to promote the continuous improvement of the ecological environment of the Loess Plateau.

2. Materials and Methods

2.1. Study Area

The Jinghe River is situated inside the central region of the Loess Plateau in China. Originating from Liupan Mountain, located at the coordinates of 106°14′–108°42′ E and 34°46′–37°19′ N in Jingyuan County, Ningxia Hui Autonomous Region, the basin spans an area of 45,421 km² and is characterized by a continental monsoon climate. The mean annual precipitation is 516.7 mm, primarily concentrated during summer, accounting for up to 60% of the total yearly rainfall with high intensity. The catchment area of the upstream of Zhangjiashan Hydrological Station is 45,180 km². The basin contains a total of three hydrological stations and eight meteorological stations. The study area is shown in Figure 1.

2.2. Data

The data needed to construct the model in this study included Digital Elevation Model (DEM), land use, soil data, meteorological, and runoff data. The DEM (Digital Elevation Model) data were procured from the SRTM-DEM (Shuttle Radar Topography Mission) with 90 m resolution, courtesy of the geospatial data cloud. Additionally, land use data with a 1 km spatial resolution were retrieved from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences. Moreover, soil data with a 1 km spatial resolution were acquired from the World Soil Database (WSD). Meteorological data from 1979 to 2018 were derived from CMADS-LV1.0 (The China Meteorological Assimilation Driving Datasets for the SWAT model Version 1.0). The runoff data for the three hydrological stations were derived from the Yellow River Basin Hydrological Yearbook.

The actual evaporation validation data were selected from China Terrestrial Evapotranspiration Dataset (TEDAC) (1982–2017) [25] and the evaporation dataset of the Global Land Data Assimilation System (GLDAS-NoahV2.1). In contrast, the validation data of soil water were selected from the NNsm global soil moisture data [26] and GLDAS-NoahV2.1 soil water content dataset.

2.3. Research Method

2.3.1. Distributed Watershed Hydrological Model SWAT

Watershed hydrological modeling can be divided into two stages based on the hydrological formation process. The first stage was the confluence of the land slope. Runoff, sediment, nutrients, and pesticides in each sub-watershed are transported to and transferred into rivers during this process. The second part is the confluence and calculation stage of the river channel, that is, the process of water flow, sediment, and non-point source pollutants in the river channel migrating to the outlet of the basin.

In the first stage of the terrestrial water cycle, the model is based on the following water balance equation:

$$S{W}_t=S{W}_0+\sum _{i=1}^t(P-Q_S-ET-W-Q_{gw})$$

(1)

where $S{W}_0$ and $S{W}_t$ are the initial and final water contents, respectively, $P$ is the precipitation on day $i$, $Q_s$ is the surface runoff on day $i$, $ET$ is the evaporation on day $i$, $W$ is the amount of water infiltrated into the deep groundwater on day $i$, and $Q_{gw}$ is the groundwater outflow on day $i$.

In the second stage of confluence and calculation, the Manning formula was used to calculate the flow discharge and velocity:

$$Q=\frac{A\cdot R^{2/3}\cdot l^{1/2}}{n}$$

(2)

$$v=\frac{R^{2/3}\cdot l^{1/2}}{n}$$

(3)

where Q is the cross-sectional flow of the river; A is the wetted area; R is the hydraulic radius; l is the bottom slope of the river; n is the Manning coefficient of the river; and v is the flow velocity.

Water flow is simulated by the Muskingum formula in river channel calculus:

$$V=K\cdot [C{Q}_{in}+(1-C)\cdot Q_{out}]$$

(4)

$$K=\frac{600\cdot L}{v}$$

(5)

where V is the river water storage capacity, K is the impounding time, C is the weight factor, Q_in and Q_out are the inlet and outlet flows, respectively, and L is the river length.

2.3.2. Hydrological Evaluation Index

To make the simulation results more reasonable, we selected the R², NSE, PBIAS, and KGE four hydrological evaluation indicators to evaluate the applicability of the model to the Jinghe River Basin.

The index calculation formula is as follows:

$$R^2={\left[\frac{{\sum }_{i=1}^n(Q_1-\overline{Q_1})(Q_2-\overline{Q_2})}{\sqrt{{\sum }_{i=1}^n(Q_1-\overline{Q_1}{)}^2{\sum }_{i=1}^n(Q_2-\overline{Q_2}{)}^2}}\right]}^2$$

(6)

$$NSE=1-\frac{{\sum }_{i=1}^n(Q_1-Q_2{)}^2}{{\sum }_{i=1}^n(Q_1-{\overline{Q}}_1{)}^2}$$

(7)

$$PBIAS=\frac{{\sum }_{i=1}^n(Q_1-Q_2)\times 100}{{\sum }_{i=1}^n{Q}_1}$$

(8)

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

(9)

$$\beta =\frac{{\mu }_s}{{\mu }_0}, \gamma =\frac{{\sigma }_s/{\mu }_s}{{\sigma }_0/{\mu }_0}$$

where Q₁ and Q₂ are the runoff observation and simulation values, respectively; ${\overline{Q}}_1$ and ${\overline{Q}}_2$ are the observed average runoff and simulated average runoff; n is the length of the time series; ${\mu }_s$ and ${\mu }_0$ are the mean values of the simulated and measured flows, respectively; ${\sigma }_s$ and ${\sigma }_0$ are the mean square deviations of the simulated and measured values, respectively; and r is the linear correlation coefficient between the simulated and measured values.

Different scholars may adopt different evaluation criteria for the indicators. The evaluation criteria adopted in this study are as follows:$R^2\ge 0.6$, $NSE\ge 0.5$, and $PBIAS$ within $\pm 25\%$, it is considered that the simulation effect is satisfactory.

2.3.3. Parameter Sensitivity Analysis

The SWAT-CUP contains five methods: SUFI-2, the generalized likelihood uncertainty algorithm, the particle swarm optimization algorithm, the parameter-solving method, and the Monte Carlo Markov chain. The SUFI-2 algorithm was used in SWAT-CUP (Calibration and Uncertainty Programs for SWAT) for model calibration and verification, parameter sensitivity, and model uncertainty analysis.

The SUFI-2 algorithm uses two methods for parameter sensitivity analysis: global sensitivity analysis and one-time single-parameter sensitivity analysis. This study adopts a global sensitivity analysis. The global sensitivity analysis approach aims to assess the combined effect of multiple parameters on the model output when they are varied simultaneously. This approach takes into account the interactions and possible non-linear relationships between the parameters. Its limitation is that the parameter range and number of runs affect the relative sensitivity of the parameters, with the advantage of producing more reliable results. Sensitivity ranking was defined by evaluating two coefficients: T-stat and p-value. The T-stat value represents the parameter sensitivity. If the absolute value of T-stat is larger, the parameter sensitivity is higher. The p-value represents the significance of the parameter sensitivity. The closer the value is to 0, the stronger the significance of the parameter sensitivity.

2.4. Model Evaluation

The SWAT model (v. 2016) was constructed for the area above the Zhangjiashan Hydrological Station. The warm-up period was 1997, the calibration period was 1998–2007, and the validation period was 2008–2017. As the SWAT model requires a consistent projection coordinate system for spatial data, the spatial data were uniformly projected, and the projection coordinate system was Beijing_1954_3_Degree_GK_CM_108E.

In GIS software (v. 10.2), the downloaded DEM data is merged and projected, and then the vector shape file of the Jinghe River Basin is used to obtain the DEM data of the Jinghe River Basin. The minimum catchment area was set to 500 km², and the study area was divided into 51 sub-watersheds, as shown in Figure 2a. To ensure the accuracy and time cost of calculation, less than 10% of land use and 15% of soil distribution were split and merged into other types. The study area was divided into 425 hydrological response units (HRUs).

Land use is a crucial factor in determining the effectiveness of the SWAT model. It directly influences how precipitation flows across the land surface and ultimately affects the simulation outcomes. The land use of the investigated area was categorized into six types: grassland (PAST), cultivated land (RICE), forest land (FRST), water area (WATR), urban and rural residential land (URBN), and unused land (BARR). Figure 2b shows the spatial arrangement of land use in the basin. Out of them, grassland areas constituted the greatest share, amounting to 46.05%. A total of 40.89% of the land is designated as agricultural land. At its minimum, the percentage of unutilized land was a mere 0.06%.

The physical properties of soil determine the movement of water and gas in the soil profile and play an essential role in the water cycle in the HRU. Establishing a database of soil physical properties requires the calculation of several parameters. Some parameters were obtained directly from the HWSD (Harmonized World Soil Database) dataset. Some parameters must be calculated using the SPAW 6.02 software. Specific parameters must be calculated using this equation. The other datasets had default values. The soil in the study area was classified into seven types: calcareous soil (CMc), calcareous alluvial soil (FLe), artificial soil (ATf), simple calcareous soil (CLh), transitional red sand soil (ARb), saturated loose lithological soil (RGe), and simple, highly active leaching soil (LVh). Calcareous and alluvial soils were the predominant soil types. The distribution of the soil types is shown in Figure 2c.

3. Results

3.1. Parameter Sensitivity Analysis

In the SWAT-CUP, SUFI-2 was used for the global sensitivity analysis. A total of 28 parameters related to hydrological processes were selected for the iterative operation, 500 times in one iteration [27]. To calibrate the model, we selectively used the top 18 parameters in the sensitivity ranking, as shown in Figure 3.

From the overall perspective shown in Figure 3, the sensitivity of the surface runoff-related parameters was the strongest, followed by that of the underground runoff-related parameters. The more sensitive parameters in the Jinghe River Basin are CN2 (the number of SCS runoff curves under wet condition II), HRU_SLP (average slope), SOL_BD (soil saturated bulk density), SLSUBBSN (average slope length), BIOMIX (biological mixing efficiency), CH_K2 (effective hydraulic conductivity of the main river bed), etc.

3.2. Model Calibration and Verification of the Discharge at Zhangjiashan Station

The year 1997 was set as the warming-up period, 1998–2007 was set as the warm-up calibration period of the model, and 2008–2017 was set as the verification period for the model, respectively. Based on the monthly flow measured at the Zhangjiashan Hydrological Station, the applicability of the model was evaluated using R², NSE, PBIAS, and KGE. An evaluation of the simulation results is presented in

Table 1. Evaluation of monthly runoff simulation result at Zhangjiashan station in the Jinghe River Basin.

Period	R2	NSE	PBIAS (%)	KGE
Calibration period (1998–2007)	0.78	0.67	−2.1	0.75
Verification period (2008–2017)	0.77	0.64	18.2	0.72
Simulation period (1998–2017)	0.76	0.62	4.8	0.72

. As shown in Table 1, the R², NSE, and KGE of the model were 0.64 and higher during the calibration and verification periods, and the PBIAS was within 20%. The accuracy of each index was higher during the calibration period than during the verification period. Thus, the model applies to the Jinghe River Basin.

Figure 4 displays the monthly flow simulation results. The simulated flow consistently matched the measured flow process. The performance of the simulation during the non-flood season fell below that of the flood season. The peak flow often occurs during the flood season, which spans from July to October. The peak timings of the simulated and measured flow processes are unified. Throughout the verification period, there was a substantial increase in the runoff peak in 2013, leading to considerable improvements in the findings of the SWAT simulation. Simultaneously, the high level of simulation value will result in continuity and impact the simulation effect in subsequent years. Peak overestimations were observed in 2014 and 2017. The simulation effect of the model was further analyzed through a correlation diagram of the measured and simulated monthly flows during the calibration and verification periods. In Figure 5, the black line represents a 1:1 proportional straight line, and the red line represents the regression line between the measured and simulated values. The scatter distribution is relatively concentrated in the calibration and verification periods, and only a small number of scatter points deviate from the linear relationship, indicating that the simulated and measured flows have good consistency in the calibration and verification periods.

3.3. Multistation Calibration and Verification of the Discharge

In a large-scale basin, the calibration and verification of the hydrological parameters only at the outlet stations of the entire basin may ignore the spatial heterogeneity of the basin. They may not meet the simulation accuracy requirements of each sub-catchment area. Therefore, flow observation data from three hydrological stations (Qingyang, Yangjiaping and Zhangjiashan) in the Jinghe River Basin were selected for further multistation calibration. The Qingyang station is located on the Malian River, a tributary of the Jinghe River, and the Yangjiaping and Zhangjiashan hydrological stations are located on the main stream of the Jinghe River. The sub-regions of the three stations are shown in Figure 6.

Among these are 11 sub-basins controlled by the Qingyang hydrological station, numbered 1–10 and 16, respectively, with a total control area of 10,603 km². There are 14 sub-basins controlled by the Yangjiaping hydrological station, numbered 13–14, 19–24, 31–35, and 38, with a total control area of 14,114 km². The Zhangjiashan Hydrological Station controls 26 subbasins, numbered 11–12, 15, 17–18, 25–30, 36–37, and 39–51, with a total control area of 20,463 km².

The Qingyang, Yangjiaping, and Zhangjiashan stations were calibrated and verified according to the order of first upstream and downstream, the first tributary, and then the mainstream.

Firstly, parameter sensitivity analyses need to be carried out for each region. The results of the parameter sensitivity ranking and final values for the three stations are shown in

Table 2. Ranking of parameter sensitivities for each station and final parameter values.

Sensitivity Ranking	Parameter at Qingyang	Parameter Value at Qingyang	Parameter at Yangjiaping	Parameter Value at Yangjiaping	Parameter at Zhangjiashan	Parameter Value at Zhangjiashan
1	CN2	82.01	CN2	61.64	CN2	50.08
2	HRU_SLP	0.54	SOL_AWC	−0.11	HRU_SLP	0.16
3	SOL_BD	−0.17	BIOMIX	0.41	SOL_BD	0.19
4	SOL_AWC	0.72	ESCO	0.95	SLSUBBSN	−0.42
5	BIOMIX	0.53	ALPHA_BF	0.63	SOL_AWC	−0.28
6	ESCO	0.08	HRU_SLP	0.35	SOL_K	−0.49
7	SOL_K	0.44	EPCO	0.95	CH_K2	3.01
8	SLSUBBSN	0.15	GW_DELAY	159.04	BIOMIX	−0.01
9	ALPHA_BF	0.04	CH_K2	120.52	ALPHA_BNK	0.86
10	CH_K2	58.61	SLSUBBSN	−0.14	EPCO	0.03

. The parameters are listed in

Table 3. Definition of parameters.

Parameters	Definition
CN2	The number of initial SCS runoff curves under wet condition II
HRU_SLP	Average gradient
SOL_BD	Soil-saturated bulk density
SOL_AWC	Effective water capacity of soil layer
BIOMIX	Biomixing efficiency
ESCO	Soil evaporation compensation factor
SOL_K	Soil-saturated hydraulic conductivity
SLSUBBSN	Average slope length
ALPHA_BF	Base flow alpha factor
CH_K2	Effective hydraulic conductivity of main channel riverbed
GW_DELAY	Groundwater hysteresis time
ALPHA_BNK	Baseflow α factor of riparian water storage
EPCO	Crop consumption compensation factor

.

The most sensitive parameter at the three stations was CN2 (number of SCS runoff curves under wet condition II); however, the values differed. The CN2 values at Qingyang, Yangjiaping, and Zhangjiashan stations were 82.01, 61.64, and 50.08, respectively. This parameter represents the runoff capacity of land and soil types in the region. There was an exponential relationship between CN2 and runoff. The larger the value, the stronger the runoff capacity and this parameter plays a decisive role in the simulation of the flood peak flow. Among the three catchments, the runoff capacity was greatest in the catchment under the control of Qingyang Station, followed by Yangjiaping. Conversely, the catchment governed by Zhangjiashan Station exhibited the weakest runoff capacity. The sensitivity of HRU_SLP (average slope) at Qingyang and Zhangjiashan stations was second and sixth at Yangjiaping station, respectively, indicating that the sensitivity of this parameter was high at all three stations. The values at Qingyang, Yangjiaping, and Zhangjiashan stations were 0.54, 0.35, and 0.16, respectively. This parameter mainly affects lateral flow, and the change in slope is positively correlated with the change in runoff yield.

Similarly, 1997 was used as the warm-up period, 1998–2007 was the calibration period, and 2008–2017 was the verification period. Based on the monthly flows measured at the three hydrological stations, the applicability of the model was evaluated using R², NSE, PBIAS, and KGE. An evaluation of the simulation results is presented in

Table 4. Evaluation of simulated monthly runoff at three stations.

	Calibration				Verification				Simulation Period
Stations	R2	NSE	PBIAS (%)	KGE	R2	NSE	PBIAS (%)	KGE	R2	NSE	PBIAS (%)	KGE
Qingyang	0.62	0.57	20.8	0.69	0.75	0.59	13.5	0.68	0.65	0.56	17.9	0.73
Yangjiaping	0.75	0.70	22.7	0.72	0.64	0.50	24.8	0.65	0.68	0.58	25.0	0.69
Zhangjiashan	0.76	0.70	6.4	0.82	0.73	0.69	15.5	0.74	0.74	0.64	11.5	0.77

.

Table 4 shows that the R² of the three stations is greater than 0.6, NSE is greater than 0.5, PBIAS is within ±25%, and PBIAS is less than 0 in the calibration period and the verification period, indicating that the simulated monthly flow is less than the observed monthly flow in both the calibration period and the verification period. In addition to the simulation effect of the Qingyang station, the verification period is better than the calibration period, and the simulation effect of the other two stations is better than the verification period. In general, the more downstream the model, the higher the simulation accuracy. Comparing the hydrological evaluation indicators of the Zhangjiazhan station rate period and validation period in Table 1 and Table 4, the results are shown in Figure 7.

Figure 7 shows that the blue area is larger than the red area, which means that the overall effect of the multi-station setting is better than that of the traditional single-station setting. The calibration period effect is better than that of the validation period for both the single-station and multi-station settings.

The monthly flow simulation results for each station are shown in Figure 8. In Figure 8a, the monthly flow trend at Qingyang station is consistent with the measured monthly flow process. The simulation results for 2001 and 2003 were the best. During the verification period, there were three years of flood peak overestimation: 2010, 2013, and 2014. In Figure 8b, the simulation effect of the Yangjiaping station in the calibration period is better than that in the validation period, and the simulation effect is better in the summer and autumn of the calibration period. The best simulation results were obtained in 1998 and 2003, and the simulation flow during the non-flood season of the validation period was generally small. The simulated values in the dry season are significantly lower than the measured values, and this difference mainly stems from the fact that the runoff in the dry season mainly relies on the generation of groundwater, and the simulation process of groundwater flow and flow in the Jinghe River Basin is quite complicated, which leads to the deviation of the simulated flow in the dry season.

Figure 9 shows the correlation between the measured and predicted monthly flows for calibration and verification periods for the three stations. By comparing the 1:1 scale straight line and the regression line in Figure 9, most of the regression lines are at the lower right of the scale straight line, indicating a trend of the simulated monthly runoff value being less than the measured monthly runoff value at the three stations throughout the simulation period. The periodic R² of Qingyang station was the smallest, at only 0.614. As shown in Figure 9a, the distribution of the points is scattered, particularly in the middle-value and high-value areas. The Zhangjiashan calibration period R² was the largest, reaching 0.76, indicating that the station had a good consistency between the measured and simulated monthly flows.

3.4. Multivariate Verification

The use of multisource observation data has advanced the basin’s hydrological simulation, expanding from single-flow verification to multivariate verification of evaporation, soil water, and runoff. Thus, after verifying the runoff data from the aforementioned stations, two actual evaporation datasets and two soil water content datasets were selected to evaluate the simulation effect of the model on the actual evaporation and soil water.

3.4.1. Actual Evaporation Validation

As a key part of the global water and energy cycle, surface evapotranspiration affects climate through a wide range of feedback acting on temperature, humidity, and precipitation. Two validation datasets for actual evaporation were chosen: the China Terrestrial Evapotranspiration Dataset (TEDAC) [25] and the GLDAS-Noah V2.1 (Global Land Data Assimilation System) evaporation dataset. China’s Terrestrial Evapotranspiration Dataset is a product that calculates the amount of water that evaporates from the land surface and is released through transpiration by plants using the complementary technique. The dataset had a monthly temporal resolution and a spatial resolution of 0.1°. NASA’s GLDAS aims to use advanced land surface modeling and data assimilation techniques to ingest satellite and ground-based observational data products, generating an optimal field of land surface state and flux. At present, GLDAS operates with four land surface models, namely Noah, Catchment, the Community Land Model, and the Variable Infiltration Capacity model. This study used the evapotranspiration dataset obtained from the GLDAS dataset, generated explicitly by the Noah land surface model. The spatial resolution was 0.25° × 0.25° and the temporal resolution was monthly.

The two actual evaporation datasets were processed in batches using MATLAB R2016b. The NC files were converted into TIFF files, and the actual evapotranspiration data of the Jinghe River Basin were obtained. In addition, the parameter values of each station were returned to the SWAT model, and the output consisted of monthly actual evapotranspiration data from 1998 to 2017.

The time series used for verification by TEDAC covered the period from 1998 to 2017, whereas the time series used for verification by GLDAS covered the period from 2001 to 2017. The applicability of the model was evaluated using four indicators: R², NSE, PBIAS, and KGE.

Table 5. Evaluation table of actual evapotranspiration verification results of each station.

	TEDAC (1998–2017)				GLDAS (2001–2017)
Sub-Region of the Stations	R2	NSE	PBIAS (%)	KGE	R2	NSE	PBIAS (%)	KGE
Qingyang	0.62	0.59	15.80	0.74	0.79	0.73	1.6	0.81
Yangjiaping	0.65	0.63	11.94	0.77	0.65	0.56	9.59	0.72
Zhangjiashan	0.67	0.63	16.69	0.74	0.64	0.54	18.55	0.69

contains the evaluation of the simulation results. Comparison was made between the average actual evaporation in each subregion of the three stations, Qingyang, Yangjiaping, and Zhangjiashan, and the data from TEDAC and GLDAS.

Table 5 shows that the R² of TEDAC and the simulation at each station is greater than 0.62, NSE is greater than 0.59, PBIAS is within ±17%, and the index accuracy is the highest at Yangjiaping station, followed by Zhangjiashan. Qingyang is the lowest, indicating that the dataset can be used for the actual evaporation verification of the SWAT model in the Jinghe River Basin. The R² of GLDAS at each station is greater than 0.66, NSE is greater than 0.54, and PBIAS is within ±19%. The index accuracy of Qingyang station is the highest, R² reaches 0.79, and PBIAS is only 1.6%. The index accuracy of the Yangjiaping and Zhangjiashan stations is not as good as that of the Qingyang station. Still, mostly satisfactory results indicate that the dataset is also suitable for the actual evaporation verification of the SWAT model in the Jinghe River Basin.

A comparison between the actual evaporation simulation value of each station and the hydrological process of the actual evaporation data from TEDAC and GLDAS is shown in Figure 10. Overall, the actual evaporation simulation values of each station were highly consistent with the change trends of the two datasets. The SWAT model can better simulate the evaporation hydrological process; however, the SWAT simulation value is generally underestimated in summer, and the actual evaporation simulation effect is better in spring, autumn, and winter. As shown in Figure 10d, at the Qingyang station, the SWAT simulation value had the highest consistency with the changing trend of the GLDAS evaporation data. From Table 5, the KGE coefficient reaches 0.81.

3.4.2. Soil Water Content Validation

Soil water is closely linked to atmospheric water, surface water, groundwater, and vegetation, and soil infiltration determines the extent to which soil water can be effectively stored. The validation data of the two soil waters were the soil moisture dataset of GLDAS-Noah V2.1 and the NNsm dataset [26] The soil moisture dataset provided by the Noah land surface model in the GLDAS model had a spatial resolution of 0.25° × 0.25° and a monthly time resolution. The NNsm is a 36 km global soil moisture dataset with similar accuracy to SMAP soil moisture on a long-term daily scale. A run slot was used to derive month-by-month and day-by-day soil moisture data from the model simulation results for 2001–2017.

The soil water data used for verification by GLDAS were monthly data from 2001 to 2007, and the soil water data used for verification by NNsm were daily data from 2009 to 2017. The model’s applicability was evaluated using four indicators: R², NSE, PBIAS, and KGE. An evaluation of the simulation results is shown in

Table 6. Evaluation of soil moisture verification results at each station.

	GLDAS				NNsm
Stations	R2	NSE	PBIAS (%)	KGE	R2	NSE	PBIAS (%)	KGE
Qingyang	0.46	0.43	19.39	0.52	0.58	0.54	20.11	0.60
Yangjiaping	0.45	0.44	23.94	0.50	0.40	0.39	22.10	0.44
Zhangjiashan	0.52	0.52	11.36	0.59	0.32	0.30	10.56	0.52

.

Based on the data presented in Table 6, the R² and NSE values for each station are higher than 0.43 when compared to the GLDAS dataset. The PBIAS falls within the range of ±24%, indicating acceptable levels of accuracy. Additionally, there is not a significant difference in the accuracy of any station index. The accuracy of each station index was relatively low compared to the NNsm dataset. Overall, the simulation of soil water in the Jinghe River Basin at each station was lower than the simulations of runoff and actual evapotranspiration.

3.5. Further Calibration and Evaluation

The optimization algorithm uses evaluation metrics as the basis for parameter optimization, which does not make full use of the physical meaning of the parameters, and the optimization algorithm has been run several times and can no longer be effective in improving the simulation. Since the actual evapotranspiration and soil moisture content data cannot now be used as target variables for optimization in the model, the parameters were manually adjusted to optimize the multivariate simulated values, and based on the validation results of runoff, actual evapotranspiration, and soil moisture, the parameters were manually adjusted based on the physical properties of each parameter in Swat-cup, with the expectation of obtaining a better evaluation result. The adjusted parameters are listed in

Table 7. Adjusted values for each parameter at each station.

Parameter at Qingyang	Parameter Value at Qingyang	Parameter at Yangjiaping	Parameter Value at Yangjiaping	Parameter at Zhangjiashan	Parameter Value at Zhangjiashan
v_CN2	82.01	v_CN2	57.65	v_CN2	40.08
v_HRU_SLP	0.54	r_SOL_AWC	−0.15	v_HRU_SLP	0.18
r_SOL_BD	−0.15	r_BIOMIX	0.41	r_SOL_BD	0.17
r_SOL_AWC	0.48	v_ESCO	0.85	r_SLSUBBSN	−0.42
r_BIOMIX	0.36	v_ALPHA_BF	0.63	r_SOL_AWC	−0.28
v_ESCO	0.45	v_HRU_SLP	0.30	r_SOL_K	−0.50
r_SOL_K	0.24	v_EPCO	0.95	v_CH_K2	3.01
r_SLSUBBSN	0.05	v_GW_DELAY	159.04	r_BIOMIX	0.00
v_ALPHA_BF	0.81	v_CH_K2	120.52	v_ALPHA_BNK	0.86
v_CH_K2	58.62	r_SLSUBBSN	−0.14	v_EPCO	0.03

.

The catchments controlled by Yangjaping and Zhangjiashan stations are similar in terms of vegetation endowment, which enhances the number of common parameter values in the simulated effects at these two stations to some extent. Through in-depth analyses, it is found that changes in the values of HRU_SLP and SOL_BD have a significant effect on runoff, especially at Yangjaping and Zhangjiashan stations, where an increase in the value of HRU_SLP directly leads to an enhancement of runoff volume. In addition to these two parameters, the main parameters of CN2 (SCS runoff curve number under wet condition II), ESCO (soil evaporation compensation factor), and SOL_AWC (effective water capacity of soil layer) were also adjusted. The value of CN2 is closely related to surface runoff, and the change in its value directly affects the accuracy of the runoff simulation, and the value of ESCO is closely related to evaporation, and the model-simulated runoff is significantly affected by the decrease in ESCO. When the value of ESCO decreases, the maximum evapotranspiration simulated by the model increases accordingly. In addition, SOL_K (soil-saturated hydraulic conductivity) is also closely related to runoff, and as the value of SOL_K increases, the runoff volume also shows an increasing trend. The fine-tuning of these parameters is essential to improve the accuracy of hydrological simulation. The adjusted parameters were brought back to SWAT and re-run to obtain the simulated monthly flows from 1997 to 2017, which were compared with the measured monthly flows of the three sites, respectively. The evaluation results obtained from the calculations are shown in

Table 8. Evaluation of monthly runoff simulation results of three stations.

Stations	R2	NSE	PBIAS (%)	KGE
Qingyang	0.65	0.58	12.40	0.76
Yangjiaping	0.69	0.57	17.00	0.73
Zhangjiashan	0.77	0.67	13.51	0.76

.

A comparison of Table 4 and Table 8 reveals that the precision of each level of the evaluation index improved after the parameter adjustment. Figure 11 shows the correlation between the measured and simulated monthly flows during the simulation period for the three stations. In the Jing River Basin, the simulation accuracy gradually increased from upstream to downstream.

The actual evaporation results of the three stations in the model were derived and compared with those of TEDAC and GLDAS, and the evaluation results are shown in

Table 9. Evaluation table of actual evapotranspiration verification results of each station after further calibration.

	TEDAC (1998–2017)				GLDAS (2001–2017)
Sub-Region of the Stations	R2	NSE	PBIAS (%)	KGE	R2	NSE	PBIAS (%)	KGE
Qingyang	0.67	0.63	17.15	0.75	0.81	0.77	3.30	0.84
Yangjiaping	0.70	0.66	14.25	0.77	0.69	0.62	12.19	0.74
Zhangjiashan	0.71	0.70	9.67	0.77	0.68	0.63	12.73	0.77

. A comparison of Table 5 and Table 9 shows that at the Qingyang and Yangjiaping stations, except for the accuracy of PBIAS, the accuracy of the other three indicators improved. The accuracy of the four hydrological evaluation indices at the Zhangjiashan station was improved.

A comparison between the actual evaporation simulation value of each station and the hydrological process of the actual evaporation data of TEDAC and GLDAS is shown in Figure 12.

The soil moisture data obtained from the analysis were compared with the soil moisture data from GLDAS and NNsm. The results of this comparison are presented in

Table 10. Evaluation of soil moisture verification results at each station after further calibration.

	GLDAS				NNsm
Stations	R2	NSE	PBIAS (%)	KGE	R2	NSE	PBIAS (%)	KGE
Qingyang	0.51	0.51	11.19	0.56	0.54	0.53	18.02	0.62
Yangjiaping	0.52	0.50	10.82	0.61	0.44	0.40	21.45	0.46
Zhangjiashan	0.54	0.52	3.27	0.63	0.32	0.31	11.30	0.55

. Comparing Table 6 and Table 10, it is evident that the evaluation indices of the three stations improved after adjusting of parameters and the performance of the simulation is improved.

The correlations between the simulated values of soil moisture with GLDAS and NNsm at the three stations are shown in Figure 13.

Figure 13 compares the 1:1 proportional line and regression line, showing that most of the regression lines are at the lower right of the 1:1 line. This indicates that the simulated soil water content was lower than the soil water content in the two datasets in the three regions during the validation period. From the regression relationship, the simulation performance of the Zhangjiashan station in Figure 13f is the best, and its linear coefficient value is 1.011, which is the closest to 1.

4. Conclusions

This study evaluates multivariate validation at numerous stations in the Jinghe River Basin on the Chinese Loess Plateau. Single-station and single-variable calibrations were performed at Zhangjiashan Hydrological Station based on the flow observation data from 1997 to 2017, followed by multistation calibration. The multisource data were then used for multivariate verification. The model’s applicability was evaluated using four hydrological simulation evaluation indices (R², NSE, PBIAS, and KGE).

The Jinghe River Basin SWAT model was constructed utilizing spatial and attribute data. A total of 51 sub-basins and 425 hydrological response units were identified. The monthly discharge simulation at Zhangjiashan station had R², NSE, PBIAS, and KGE values of 0.75, −2.1%, 0.67, and 0.78, respectively, during calibration while 0.77, 0.64, 18.2%, and 0.72, respectively, during validation. The results were all satisfactory, showing that the model is adequate.

Multistation calibration utilized spatial differences. Qingyang, Yangjiaping, and Zhangjiashan stations had R², NSE, and KGE above 0.57 and PBIAS within 25%. During verification, R², NSE, and KGE were above 0.50, and PBIAS was within 25%. The Qingyang station outperformed the other two stations during verification, whereas the latter performed better during calibration. Typically, the impact of the downstream simulation was superior. The monthly runoff simulation accuracy of the Zhangjiashan Station was improved after partition calibration.

In the multivariate validation, runoff performed best, followed by actual evaporation, whereas the soil water content was slightly worse. In the two datasets of actual evaporation data, TEDAC performed better in validation at the Yangjiaping and Zhangjiashan stations, whereas GLDAS performed better at the Qingyang station. The evaluation results for each hydrological variable improved after further calibrating the parameters [28]. Multivariate validation using multisource data improved distributed hydrological simulation at river basin scale. This study provides a reference for evaluating water resources and parameter calibration of hydrological models at the river basin scale.