*2.2. Model Description*

The SWAT model was developed by the Agricultural Research Service of the United States Department of Agriculture (USDA-ARS) to predict impacts of climate and land use/management changes on water, sediment, and agricultural chemical yields [1,38–40]. This physically based distributed model has been widely applied for the estimation of watershed issues under climate change and human activities at watershed scales. The main outputs of SWAT are streamflow, ET, soil water, groundwater, sediment load, and nitrogen and phosphorus loads [41–44]. More details can be found in the SWAT input/output documentation [45]. SWAT is considered to be one of the best tested models in the field of hydrological simulation and has been widely used globally. In particular, the model is open-source and can facilitate users' modifications to meet their specific needs. Therefore, SWAT was selected and modified technically, in this study, to adapt it for ET calibration, using SWAT Calibration and Uncertainty Programs (SWAT-CUP). SWAT-CUP is an integrated program for sensitivity analysis, calibration, validation, and uncertainty analysis of SWAT models.

#### *2.3. Model Input and Setup*

The SWAT model requires inputs on weather and topography data. A digital elevation model (DEM) was obtained from NASA's Shuttle Radar Topography Mission (SRTM V3.0) with 3 arc second (~90 m) resolution [46]. The land use (2015) and soil type maps with 1 km resolution were obtained from the Ecological and Environmental Science Data Center for West China (http://westdc.westgis.ac.cn). Daily meteorological data of five weather stations (1951–2016) were obtained from the National Meteorological Information Center (http://data.cma.cn), including precipitation (P), maximum and minimum air temperature, relative humidity (RH), wind speed (WS), and sunshine duration (SD). The solar radiation (SR) was calculated based on SD and latitude [41]. The Geographic Information System (GIS) interface, ArcSWAT (version 2012), was used to delineate the watershed, resulting in 86 sub-basins and 2141 HRUs (Hydrologic Response Units).

#### *2.4. Model Calibration*/*Validation and Uncertainty Analysis*

SWAT-CUP was employed in this study for model calibration and uncertainty quantification. We used two hydrological variables, namely, streamflow from two stations (Zhaoping and Jingnan as shown in Figure 1) and ET of the entire GRB for model calibration and validation. Monthly streamflow data at these two gaging stations were obtained from the local hydrological department, and the remotely sensed ET (MOIDS16A2 ET) was provided by NASA Land Processes Distributed Active Archive Center (LP DAAC, https://modis.gsfc.nasa.gov/data/dataprod/mod16.php), covering the period from 2000 to 2014.

In this study, the SWAT model was examined with streamflow and MODIS ET using a 13 year (1998–2010) record for calibration and another 5 year (2011–2015) period for validation (Figure 2). The influences of initial conditions were minimized using a 5 year (1993–1997) warm-up period [40]. Based on a review of literature related to SWAT calibration and our own experience [20,43,46–51], we selected seven important parameters for optimization (this step can be implemented using parameter sensitivity analysis for a customized project), and their adjustment modes and initial ranges are defined in Table 1. These seven parameters can affect the hydrologic responses of surface runoff subsurface lateral flow and channel routing procedures in SWAT simulations. As stated in Table 1, CN2, SOL\_AWC, and SOL\_K were adjusted using a relative adjustment approach (i.e., multiplying the initial value by a coefficient (1+ a given value)) to maintain their spatial variation, with adjustment ranges between −20% and 20% [20,33,52]. ALPHA\_BF, CH\_K2, ESCO, and SURLAG were replaced by a new given value, with adjustment ranges based on prior knowledge of the GRB and a review of existing publications [15,30,53]. Before the scenario design, we compared the results of 1000, 2000, 3000, and 5000 model runs and found that the results showed no significant difference after 2000 runs. Hence, we used SUFI-2, PSO, and combined SUFI-2–PSO approaches for uncertainty quantification

and reduction, with a total of 2000 model runs for each approach. Below is a brief introduction to these approaches. *Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 5 of 26

**Figure 2.** Availability of the data (streamflow and evapotranspiration (ET)) used in this study. **Figure 2.** Availability of the data (streamflow and evapotranspiration (ET)) used in this study.


**Table 1.** Calibrated parameters of the SWAT model for the Guijiang River Basin.

v\_SURLAG.bsn Surface runoff lag time [0.1, 4] Note: "v\_" means the existing parameter value is to be replaced by a given value, "r\_" means an Note: "v\_" means the existing parameter value is to be replaced by a given value, "r\_" means an existing parameter value is multiplied by (1+ a given value).

existing parameter value is multiplied by (1+ a given value).

#### 2.4.1. SUFI-2

2.4.1. SUFI-2 SUFI-2 is a Bayesian-based optimized algorithm [54] that quantifies the uncertainties using sequential and fitting processes. In the parameter identification process, a sensitivity matrix and Hessian matrix are calculated to update and narrow the parameter ranges. A short step-by-step SUFI-2 is a Bayesian-based optimized algorithm [54] that quantifies the uncertainties using sequential and fitting processes. In the parameter identification process, a sensitivity matrix and Hessian matrix are calculated to update and narrow the parameter ranges. A short step-by-stepdescription of SUFI-2 is as follows:

description of SUFI-2 is as follows: **Step 1:** Define an objective function. We selected the commonly used Nash–Sutcliffe coefficient (NSE) [55] (Equation (A1)) as the objective function in this study. R2 (Equation (A2)) and percentage of bias (PBIAS) (Equation (A3)) were also used for evaluating model performance. **Step 1:** Define an objective function. We selected the commonly used Nash–Sutcliffe coefficient (NSE) [55] (Equation (A1) in Appendix A.1) as the objective function in this study. R<sup>2</sup> (Equation (A2) in Appendix A.2) and percentage of bias (PBIAS) (Equation (A3) in Appendix A.3) were also used for evaluating model performance.

**Step 2:** Set the initial ranges for the selected parameters which are to be calibrated (Table 1). Several sets of parameters are generated by Latin hypercube sampling (LHS) that are used for model simulations. Then, the NSE for each variable is calculated to assess the model performance. **Step 2:** Set the initial ranges for the selected parameters which are to be calibrated (Table 1). Several sets of parameters are generated by Latin hypercube sampling (LHS) that are used for model simulations. Then, the NSE for each variable is calculated to assess the model performance.

**Step 3:** Evaluate each sampling round by a series of measures [22]. The sensitivity matrix J, the Hessian matrix H, and the parameter covariance matrix C are calculated as: **Step 3:** Evaluate each sampling round by a series of measures [22]. The sensitivity matrix J, the Hessian matrix H, and the parameter covariance matrix C are calculated as:

$$J\_{\rm ij} = \frac{\Delta \mathbf{g}\_i}{\Delta b\_j} i = 1, \dots, \mathbf{C}\_{2'}^n \ j = 1, \dots, m \tag{1}$$

ଶ(்)ିଵ (3)

is the number of all

<sup>ଶ</sup> is the variance of the

=

objective function values, i.e., NSEs, in this study.

$$H = \mathbf{J}^{T}\mathbf{J}\tag{2}$$

$$\mathbf{C} = \mathbf{S}\_{\mathcal{S}}^2 \mathbf{(J}^T \mathbf{J})^{-1} \tag{3}$$

where *b* is the parameter vector and *g* is the objective function, *C n* 2 is the number of all combinations of the two simulations, *m* is the number of parameters, *S* 2 *g* is the variance of the objective function values, i.e., NSEs, in this study.

**Step 4:** Quantify the uncertainty. The 95% prediction uncertainty (95PPU) [22], *R*-factor, and *P*-factor are used to illustrate the fitting degree and uncertainty.

The average width of 95PPU is calculated as:

$$\overline{d\_{\mathcal{Q}}} = \frac{1}{K} \sum\_{l=1}^{K} (Q\_{ll} - Q\_{L}) \tag{4}$$

As in Equation (1), *K* is the number of observations. The subscripts "U" and "L" mean upper (97.5th) and lower (2.5th) boundary of 95PPU, respectively.

*R*-factor is the relative width of 95PPU and *P*-factor is the percentage of observations bracketed by the 95PPU [9,56]. These are calculated by the following equations:

$$R - factor = \frac{d\_{\mathcal{Q}}}{\sigma\_{\mathcal{Q}}} \tag{5}$$

$$P - factor = \frac{k\_c}{K} \tag{6}$$

where σ*<sup>Q</sup>* is the standard deviation of variable *Q* and *k<sup>c</sup>* is the number of values covered by 95PPU.

**Step 5:** Update the parameter ranges for further iteration.

#### 2.4.2. PSO

PSO is a kind of self-adaptive random algorithm based on a group hunting strategy [57,58]. The particles are initially randomized in the given ranges followed by an iterative search for optima. Two optima, pbest and gbest, are defined for the update of particles. Pbest is the best solution for a particle, and gbest is the global optimum obtained by any particle in the population [59].

The PSO framework in SWAT-CUP is as follows:

**Step 1:** Initialize all particles. In SWAT-CUP, users define the number of simulations (particles) S and iterations I. Then, S sets of parameter values are generated randomly.

**Step 2:** Calculate the fitness value NSE, which is also the objective function in SUFI-2, obtain pbest and gbest, and compare them with previous values. If the present value is better, then reset the best value.

**Step 3:** Calculate the velocity and position of each particle for the next iteration.

**Step 4:** Implement Step 2 using the updated particles for further iterations until the number of iterations reaches I.

#### 2.4.3. SUFI-2 and PSO Combination and Scenario Setting

SUFI-2 enables precise position of parameter ranges that are usually narrower than the last iteration. However, the posterior parameter ranges need to be checked and manually adjusted by users for the next iteration. SWAT-CUP can automatically provide a set of parameter ranges after one iteration. Users can use the recommended ranges or adjust them to ensure their physical meaning for the next iteration. PSO is an effective automatic global optimization method in which the optimization efficiency decreases with the increase in simulation times. This means PSO can achieve a satisfactory simulation with less calibration time. To improve optimization efficiency and decrease modeling uncertainty, we combined SUFI-2 and PSO in different arrangements as listed in Table 2. The primary principle of the combination approach is to obtain posterior parameter ranges from SUFI-2, followed by optimizing parameters using PSO. Approaches No.1, No.2, and No.6 (using SUFI-2 or PSO only) were set as the reference to compare the effects of different combined approaches. Approaches No.3, No.4, and No.5 were the different combinations of SUFI-2 and PSO algorithms.


**Table 2.** Scenario setting of different approaches.

Note: "R" means calibration using streamflow data only, "RE' means calibration using streamflow and ET simultaneously, and the number of "\*" represents the number of manual operations.

Each approach mentioned above was used for streamflow calibration only (R) and streamflow plus ET calibration (RE), leading to twelve scenarios in total (scenario codes are listed in Table 2).
