Enhancing Soil and Water Assessment Tool Snow Prediction Reliability with Remote-Sensing-Based Snow Water Equivalent Reconstruction Product for Upland Watersheds in a Multi-Objective Calibration Process

Abstract

Precipitation occurs in two basic forms defined as liquid state and solid state. Different from rain-fed watershed, modeling snow processes is of vital importance in snow-dominated watersheds. The seasonal snowpack is a natural water reservoir, which stores snow water in winter and releases it in spring and summer. The warmer climate in recent decades has led to earlier snowmelt, a decline in snowpack, and change in the seasonality of river flows. The Soil and Water Assessment Tool (SWAT) could be applied in the snow-influenced watershed because of its ability to simultaneously predict the streamflow generated from rainfall and from the melting of snow. The choice of parameters, reference data, and calibration strategy could significantly affect the SWAT model calibration outcome and further affect the prediction accuracy. In this study, SWAT models are implemented in four upland watersheds in the Tulare Lake Basin (TLB) located across the Southern Sierra Nevada Mountains. Three calibration scenarios considering different calibration parameters and reference datasets are applied to investigate the impact of the Parallel Energy Balance Model (ParBal) snow reconstruction data and snow parameters on the streamflow and snow water-equivalent (SWE) prediction accuracy. In addition, the watershed parameters and lapse rate parameters-led equifinality is also evaluated. The results indicate that calibration of the SWAT model with respect to both streamflow and SWE reference data could improve the model SWE prediction reliability in general. Comparatively, the streamflow predictions are not significantly affected by differently lumped calibration schemes. The default snow parameter values capture the extreme high flows better than the other two calibration scenarios, whereas there is no remarkable difference among the three calibration schemes for capturing the extreme low flows. The watershed and lapse rate parameters-induced equifinality affects the flow prediction more (Nash-Sutcliffe Efficiency (NSE) varies between 0.2–0.3) than the SWE prediction (NSE varies less than 0.1). This study points out the remote-sensing-based SWE reconstruction product as a promising alternative choice for model calibration in ungauged snow-influenced watersheds. The streamflow-reconstructed SWE bi-objective calibrated model could improve the prediction reliability of surface water supply change for the downstream agricultural region under the changing climate.

Highlights:

(1) Calibration of snow parameters with reference streamflow might have a negative impact on the streamflow prediction.

(2) Multi-objective calibration behaves more reliably than single-objective calibration.

(3) Snow water-equivalent predictions mainly deviate at low elevation bands among different calibration strategies.

(4) Equifinality issues affect streamflow predictions more than snow water-equivalent predictions.

1. Introduction

Upland watersheds are of vital importance for freshwater supplies, since snowpack accumulated in the regional mountains acts as a seasonal water reservoir that provides major water storage for urban and agricultural users [1,2,3,4,5,6,7]. Mountain snowpacks store water during winter and discharge it, starting from spring or summer. As the global temperature is projected to rise, it becomes increasingly important to evaluate the vulnerability of mountain snow to climate change, since temperature variability could influence the spatiotemporal distribution and persistence of mountain snow, which further affects the seasonal streamflow characteristics, hydropower production, and surface water supply downstream [8,9,10,11,12,13]. A warmer climate is also expected to exaggerate the frequency and intensity of extreme events, such as droughts and floods [14,15,16,17,18,19,20,21]. Despite these global trends, regional differences exist, which, in some cases such as California, show exacerbated responses in mountain snowpacks to changes in the climate. California has experienced progressively higher winter air temperatures and more summer heat waves in recent years, which have resulted in a reduced mountain snowpack during the winter and, in part, contributed to the multi-year drought from 2012 to 2016 [22,23,24,25]. Agriculture in the Central Valley is heavily reliant on surface water that originates from the Sierra Nevada snowmelt [26]. Hence, earlier snowmelt caused by a warming climate and a diminishing snowpack would have serious effects on the Central Valley agriculture production and ecosystem function [27,28,29].

Hydrologic models are commonly applied to evaluate the impact of climate change on hydrologic systems and water resources [30,31]. In snow-dominated regions, the performance of a watershed-scale hydrologic model is significantly affected by the model’s ability to correctly capture the partitioning of precipitation into its solid state (snow) and liquid state (rain) in space and time [32,33,34]. The Soil and Water Assessment Tool (SWAT) is a widely used semi-distributed watershed model, which has a snow module that adopts the temperature-index approach to generate snow output for different elevation bands as snow water-equivalent (SWE), defined as the water contents of the snowpack in units of precipitation [35]. The calibration parameters controlling hydrological processes in SWAT can be broadly classified into three categories: (1) watershed parameters, which describe the general physiological characteristics/properties of the watershed, such as the curve number and the baseflow separation factor (alpha); (2) snow parameters, which include seven parameters that directly control snow accumulation and snowmelt (SFTMP, SMTMP, SMFMX, SMFMN, TIMP, SNOCOVMX, and SNO50COV); and (3) elevation lapse rate parameters, which control the temperature lapse rate and precipitation lapse rate. The descriptions of each parameter are shown in Table 1. Typically, SWAT is calibrated with streamflow observation at the watershed outlet by finding the proper parameter values from their common ranges that fit the streamflow observation best [36,37,38]. However, the single-objective calibration of all the parameters with only one reference streamflow might lead to significant prediction bias when using the model in forecast mode to predict future streamflows in response to climate scenarios or when using the calibrated model to predict variables other than flow. The streamflow integrates the contributions of different hydrological components in natural watersheds, and it is difficult to calibrate models such that they accurately capture these specific contributions from individual components, which can easily lead to model equifinality. Equifinality indicates a given model state can be reached by many potential means or parameter combinations [39]. Comparatively, multi-objective model calibration has shown its benefits for prediction reliability compared with single-objective calibration and reducing the chance of equifinality [40,41,42,43]. However, only very few studies have included observed/reference snow information into the model calibration process due to a lack of historical SWE reference data available for mountainous watersheds. The Snow Telemetry (SNOTEL) network, which was initially installed in the mid-1960s, is composed of over 800 automated data collection sites located in remote, high-elevation mountain watersheds in the Western United States to monitor snowpack and record snow water-equivalent and snow depth information [44]. The SNOTEL ground-based SWE observations can be applied for the model calibration process for those watersheds where the in-situ observation site density is high and equally distributed. However, in a general view, the current SNOTEL sites are still not enough, and SWE records for many snow-dominated watersheds are scarce, thus calling the academic communities to find new alternatives that have wider applications. Estimating accurate watershed-wide snow water-equivalents (SWE) is one of the most important unsolved problems in mountain hydrology [45,46]. Direct sensing of SWE is not reliable in most cases due to the rugged topography and deep snow in mountain regions [47]. With the rapid development of measurement devices (i.e., LiDAR-Light Detection and Ranging), data-processing tools, statistical applications, and remote-sensing techniques), snow reconstruction products are becoming increasingly available for major mountain ranges in the United States [48,49], thus making the calibration of SWAT models with both streamflow and SWE reference datasets feasible [50,51,52].

Based on the aforementioned discussions, this study compares applying ParBal-reconstructed SWE as the second calibration objective, in addition to streamflow with single-objective calibration (with streamflow only), and investigates the way of assigning snow parameter values (i.e., using default values, automatic calibration values, or manual calibration values) in affecting the model prediction reliability. Specifically, the objectives are to: (1) investigate the influence of snow parameters on streamflow and SWE predictions, (2) evaluate the influence of applying SWE reconstruction data as the second reference dataset in addition to streamflow in the SWAT model calibration process, and (3) assess if and to what extent equifinality occurs as a result of calibrating the watershed and elevation lapse rates parameters. Practically, this study evaluates the Parallel Balance Model SWE reconstruction product as a potential alternative for ground-based SWE records in the model calibration process, and the remote-sensing-based SWE reconstruction product has wider applications for ungauged watersheds. The streamflow-reconstructed SWE-calibrated hydrological model has high prediction reliability for investigating the snowmelt-driven surface water supply change for the agricultural area in Southern Central Valley.

2. Study Area

The study area is located in the Southern Sierra Nevada Mountain range in California, USA. The study area encompasses four watersheds, the Kings watershed (area: 4410 km², elevation: 140 m–4299 m), the Kern watershed (area: 5360 km², elevation: 760 m–4405 m), the Tule watershed (area: 1400 km², elevation: 198 m–3112 m), and the Kaweah watershed (area: 960 km², elevation: 151 m–3825 m), that are draining into the Southern Central Valley, one of the most important agricultural regions in the world [53,54] (Figure 1). The lower elevation of the Southern Central Valley region has a Mediterranean climate with a pronounced cool and moist late-fall and winter seasons and warm and dry spring and summer seasons. The higher elevation on the east side of the area is affected by both rain and snow. The mean annual precipitation varies between 150 mm in the valley to 1000–2000 mm in the mountains [55,56]. The primary sources of precipitation during the rainy season are low-pressure disturbances that move in from the northwest Pacific. In addition, storms from the southwest containing abundant subtropical moisture, so-called atmospheric rivers, can generate heavy localized precipitation and high runoff from the surrounding mountains. Most of the annual precipitation occurs from November to March. In the Sierra Nevada high-elevation ranges, snowfall provides the majority of the annual precipitation, and seasonal snowmelt is a major water supply for agriculture and domestic usage in the downstream area [53,57]. The snow-covered area in each watershed is approximately 79%, 69%, 70%, and 49%, respectively, for the Kings, Kern, Tule, and Kaweah watersheds. Land use in these upland watersheds ranges from grass rangeland and oak woodland in the low elevations to mixed-conifer forest at the middle elevations and alpine vegetation with dwarf pine and hemlock trees at high elevations.

3. Methods

3.1. Data

ArcSWAT (an ArcGIS add-in), was used to develop the SWAT models for each watershed with topography, land use, and soils data. In ArcSWAT, each watershed was discretized into sub-watersheds based on the topography using a United States Geological Survey (USGS) 30-m resolution, 1-arc second digital elevation model (DEM, http://ned.usgs.gov). A DEM was also used to derive the terrain slope, another variable that is used, together with soil type and land use, to divide each sub-watershed further into multiple hydrological response units (HRUs). The topography, land use, and soil data used in this study were obtained from the National Elevation Dataset 30-m Digital Elevation Model (NED DEM, http://ned.usgs.gov), the National Land Cover Dataset 2001 land use data (NLCD, https://www.mrlc.gov/data/nlcd-2001-land-cover-conus), and the State Soil Geographic dataset (STATSGO, https://www.nrcs.usda.gov/) respectively. When the SWAT models were built, Daymet (https://daymet.ornl.gov/) historical daily precipitation and temperature data (1-km by 1-km resolution, 1981–2013) were used as the climate drivers for the SWAT model to calibrate the streamflow and SWE. Monthly referenced streamflow data recorded at the outlets of the four watersheds were used to calibrate the model. Monthly reference streamflow data was calculated based on the observed reservoir operation data and was obtained from the California Department of Water Resource (DWR). The Parallel Energy Balance Model (ParBal) snow reconstruction dataset was used as the snow reference data for SWAT model calibration and validation. The daily ParBal SWE data was obtained from the University of California, Santa Barbara, with 500-m spatial resolution for the years 2001 to 2017 [58]. Reconstruction of the snow water-equivalent from remotely sensed data is a technique where the snowpack is built up in reverse from melt-out to peak using downscaled satellite-based energy balance forcings and a fractional snow-covered area, albedo, and snow grain size data. The Parallel Balance Model (ParBal) SWE data is a promising reconstructed SWE product derived from Moderate Resolution Imaging Spectroradiometer (MODIS) snow-covered area and grain size data. Although the ParBal data is not the direct observations, considering the limited coverage of ground-based SWE observations in the upland watersheds, the ParBal dataset becomes an important alternative for in-situ SWE observations.

3.2. SWAT Model Implementation and Calibration

3.2.1. SWAT Model Implementation

The SWAT Calibration and Uncertainty Program (SWAT-CUP, version 2012) was applied to calibrate the SWAT models. The Sequential Uncertainty Fitting (SUFI-2) algorithm was selected as the calibration technique in SWAT-CUP. SUFI-2 finds an optimal set of parameter values that best fit the reference data by narrowing down their predefined initial ranges over multiple iterations [59,60]. A total of 17 parameters were selected for the calibration of SWAT for each of the four upland watersheds (

Table 1. Parameters included in the SWAT model calibration, their initial ranges, and the final calibrated parameter values.

Parameter		Description	Initial Range	Kings	Kern	Tule	Kaweah
Snow Parameters	V_SFTMP	Snowfall temperature (°C)	0–5	4	4.8	4.58	4
	V_SMTMP	Snowmelt base temperature (°C)	0–5	5	5	5	4.5
	V_SMFMX	Melt factor for snow on June 21 (mm H2O/°C-day)	0–10	5	8	4.59	5
	V_SMFMN	Melt factor for snow on December 21 (mm H2O/°C-day)	0–10	3	0.1	0.62	4
	V_TIMP	Snow pack temperature lag factor	0.8–1	0.88	0.9	0.86	0.8
	V_SNOCOVMX	Minimum snow water content that corresponds to 100% snow cover, SNO100 (mm H2O)	0–500	400	320	360.79	490
	V_SNO50COV	Fraction of snow volume represented by SNOCOVMX that corresponds to 50% snow cover	0.01–0.99	0.4	0.1	0.66	0.03
Watershed Parameters	r_HRU_SLP	Average slope steepness	−0.84–0.84	0.21	−0.83	−0.79	−0.81
	r_CN2	SCS curve number	−0.5–0.5	0.48	−0.11	−0.49	−0.34
	V_ALPHA_BF	Baseflow alpha factor (days)	0–0.99	0.92	0.06	0.53	0.58
	V_GW_DELAY	Groundwater Delay (days)	0–500	97.5	169	95.1	104.32
	r_SOL_AWC	Available water capacity of the soil layer (mm H2O/mm soil)	−0.5–0.5	−0.36	−0.013	0.22	−0.044
	V_GWQMN	shallow aquifer required for return flow to occur (mm)	0–5000	4175	3469	931.16	2654.74
	V_GW_REVAP	Groudwater “revap” coefficient	0.02–0.2	0.18	0.096	0.073	0.085
	V_RCHRG_DP	Deep aquifer percolation factor	0–0.9	0.2	0.44	0.069	0.4
Elevation Lapse rate Parameters	V_TLAPSE	Temperature lapse rate (°C/km)	−10–0	−7	−8.5	−9.12	−7
	V_PLAPSE	Precipitation lapse rate (mm H2O/km)	0–100	62.61	33.04	27.01	99.16

SWAT refers to the Soil and Water Assessment Tool. The descriptions of each parameter are shown in Table 2.

). These parameters were selected following the SWAT documentation and sensitivity analyses conducted on the SWAT model parameters by the California Department of Water Resource [36,61]. The 17 selected parameters were divided into three groups comprising snow parameters, watershed parameters, and elevation lapse rate parameters, as mentioned before (Table 1). It needs to be pointed out that the elevation lapse rate parameters (i.e., temperature and precipitation lapse rates, Tlapse and Plapse) could affect the partitioning between rain and snow from precipitation, as well as the distribution of snow in the high elevations, although they are not direct snow parameters. To account for orographic effects in mountainous watersheds, which both influence temperature and precipitation, SWAT allows the definition of up to 10 elevation bands in each sub-watershed. For each elevation band, SWAT simulates the accumulation, sublimation, and melting of snow. Outputs from each elevation band are averaged for each sub-watershed and reported in output files. Based on previous studies, it has been shown that, typically, using 3–5 elevation bands satisfies the accuracy of most model predictions while keeping the computational demand (e.g., running time) low [62]. In this study, five elevation bands were defined for each of the four upland watersheds (Kings, Kern, Tule, and Kaweah watersheds) [63,64,65]. The developed SWAT models were used to test and compare the different calibration strategies.

3.2.2. SWAT Model Calibration

Three calibration scenarios were applied in this study. Scenario S1 and S2 applied the single-objective calibration, where the simulated streamflow was evaluated with the reference streamflow at the watershed outlet. The S3 scenario involved a multi-objective calibration, where, in addition to the streamflow, a second calibration objective was set for each watershed, consisting of the calibration of SWE at the sub-watershed scale using ParBal-reconstructed SWE reference data. The following describes in more detail how SWAT was calibrated for each calibration scenario. In the S1 calibration scenario, the SWAT models were calibrated with monthly streamflow reference data at the watershed outlet using the watershed parameters and elevation lapse rate parameters only. The snow parameters were left at their default values (

Table 2. Description of snow parameters and their default values in SWAT.

	Description	Default Value in SWAT-CUP
V_SFTMP	Snowfall temperature (°C)	1.0
V_SMTMP	Snowmelt base temperature (°C)	0.5
V_SMFMX	Melt factor for snow on 21 June (mm H2O/°C-day)	4.5
V_SMFMN	Melt factor for snow on 21 December (mm H2O/°C-day)	4.5
V_TIMP	Snow pack temperature lag factor	1.0
V_SNOCOVMX	Minimum snow water content that corresponds to 100% snow cover, SNO100 (mm H2O)	1.0
V_SNO50COV	Fraction of snow volume represented by SNOCOVMX that corresponds to 50% snow cover	0.50

The default values are found in the SWAT Calibration and Uncertainty Program (SWAT-CUP) 2012 version input/output documentation, chapter 4 [63]. SWAT refers to the Soil and Water Assessment Tool.

). In the S2 scenario, instead of applying the default snow parameter values, the snow parameters were also calibrated with reference streamflow data at the watershed outlet, in addition to the watershed parameters and elevation lapse rate parameters. In the S3 calibration scenario, a new two-stage lumped calibration strategy was developed, and all three groups of parameters were calibrated with both the monthly streamflow data at the watershed outlet and daily SWE data in the snow-influenced sub-watersheds. The total calibration period of the SWAT models was from 1981–2007, while the first three years were used as the warm-up period. Initially, the SWAT model conducted one iteration that consisted of 500 simulations. The simulation outcomes were evaluated using standard statistical measures of agreement between the observed and simulated streamflow, including the Coefficient of Determination (R², Equation (1)), Nash-Sutcliffe Efficiency (NSE, Equation (2)), and p-value and r value. These performance measures are frequently used by many hydrological modeling studies. Further descriptions of these methods can be found below. Unless all the performance measures are satisfactory, a new iteration will be initiated again with another 500 simulations.

The Coefficient of Determination (R²) is calculated using Equation (1). Where Obs_t is the observation data at time step t, Sim_t is the simulation data at time step t, and $\overline{\mathrm{Obs}}$ is the average observation value.

$${\mathrm{R}}^2=\frac{{\left[{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}\left({\mathrm{Obs}}_{\mathrm{t}}-\overline{\mathrm{Obs}}\right)\left({\mathrm{Sim}}_{\mathrm{t}}-\overline{\mathrm{Sim}}\right)\right]}^2}{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Obs}}_{\mathrm{t}}-\overline{\mathrm{Obs}}\right)}^2{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Sim}}_{\mathrm{t}}-\overline{\mathrm{Sim}}\right)}^2}$$

(1)

The Nash-Sutcliffe Efficiency (NSE) is a normalized value determining the relative magnitude of the residual variance compared to the observed data variance. The NSE ranges between −1 and 1, with NSE = 1 being the optimal value. The mathematical equation for calculating NSE is shown in Equation (2). Where Obs_t is the observation data at time step t, Sim_t is the simulation data at time step t, and $\overline{\mathrm{Obs}}$ and $\overline{\mathrm{Sim}}$ are the average observation and average simulation values between time step 1~T, respectively.

$$\mathrm{NSE}=1-\frac{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Obs}}_{\mathrm{t}}-{\mathrm{Sim}}_{\mathrm{t}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Obs}}_{\mathrm{t}}-\overline{\mathrm{Obs}}\right)}^2}$$

(2)

The p-value and r value are used in the streamflow calibration to evaluate the prediction ensemble relative to the observations. The p-value is used to assess the percentage of reference data bracketed by the 95% prediction uncertainty calculated at the 2.5 and 97.5 percentiles of the cumulative distribution of the simulated variables (95PPU), and the r factor is the thickness of the 95PPU envelop.

The multi-objective lumped calibration strategy applied in the S3 scenario included two calibration stages. In the first step, the SWAT model was automatically calibrated with monthly reference streamflows at the watershed outlet using the calibration period from 1981–2000 (the first three years were used as a warm-up period, and all three groups of parameters were calibrated). In the second step, the SWAT model parameters obtained in the previous step were further manually calibrated with referenced SWE data at the sub-watershed scale for 2001–2007. Here, the best parameter set derived from the first step was used as the initial values for the SWE calibration in the second step. In each calibration round, the snow parameters and elevation lapse rate parameters were modified within their respective value ranges before the SWAT-CUP simulation was conducted. This process was repeated until the simulated SWE reached a satisfactory level of agreement with the reference SWE data. During the SWE calibration process, it was ensured that the quality of the simulated streamflow with the modified parameters stayed at acceptable fit levels.

For all three scenarios, the SWAT model performance and calibrated parameters were evaluated during the validation period (2008–2013). Both the monthly streamflow at the watershed outlet and daily SWE prediction of each sub-watershed were evaluated with reference data to assess the performance of the calibrated models. The choice of periods for model calibration or validation was mainly determined by the availability of reference data. Since we had longer streamflow reference data (1981–2013), whereas the ParBal reconstructed SWE data was relatively shorter (available starting from 2001), the period before 2001 (1981–2000, the first three years were used as the warm-up period) was applied as the streamflow calibration period (stage 1 calibration). The years after 2001 were split into two generally equal time spans (2001–2007 was used as the SWE calibration period, and 2008–2013 was used as the validation period) to ensure each period contained both dry and wet years. A longer validation period (i.e., validation starting earlier) was not applied, since it would squeeze the SWE calibration period and affect model calibration reliability, and the current way of splitting the periods was the ideal choice. The details for all three calibration scenarios are summarized in

Table 3. Description of the three SWAT model calibration scenarios.

Scenarios		S1	S2	S3
Streamflow calibration	Include?	YES	YES	YES
	Period	1981–2007	1981–2007	1981–2000
	Reference data	Outlet streamflow	Outlet streamflow	Outlet streamflow
Snow water equivalent calibration	Include?	No (with default snow parameter values applied)	Yes (calibration of snow parameters)	Yes (calibration of snow parameters)
	Period	N/A	1981–2007	2001–2007
	Reference data	N/A	Outlet streamflow	Sub-watershed SWE
Validation	Include?	YES	YES	YES
	Period	2008–2013	2008–2013	2008–2013
	Validation objectives	Outlet Streamflow and sub-watershed SWE	Outlet Streamflow and sub-watershed SWE	Outlet Streamflow and sub-watershed SWE

S1–S3 refer to Scenario 1, Scenario 2, and Scenario 3, respectively. SWE: snow water-equivalent.

and Figure 2.

3.2.3. Evaluation of Calibration Scenarios

The predicted streamflow from the three calibration scenarios were evaluated with a reference streamflow for the calibration period and validation period, respectively. The NSE and R² values were compared among the three scenarios to assess the impact of using a single vs. multi-objective calibration on the predicted flow and SWE. In addition, extreme flows are calculated and compared for the three modeling scenarios. Extreme flow conditions are strongly related to natural hazards such as droughts (an extreme low-flow condition) and floods (an extreme high-flow condition), thus needing more attention. Especially in recent years, continuous severe drought conditions (2012–2016) have significantly affected California and caused huge amounts of economic loss and property damage [66,67]. For this reason, we calculated the 10%, 20%, 80%, and 90% quantile flows based on the streamflow predictions for 2001–2013 for all three scenarios. The calculated low flows and high flows were compared and evaluated among the three modeling scenarios.

In addition to the streamflow, the SWE predictions were also evaluated for the calibration and validation periods at both the sub-watershed scale and elevation band scale. For the SWE evaluation, the NSE values were calculated for all snow-dominated sub-watersheds for all three calibration scenarios, and their corresponding distributions (expressed as box-whisker plots) were evaluated. Since investigating the predicted SWE at the sub-watershed scale alone could not provide a clear picture about how the three modeling scenarios performed at different elevations, we also assessed the model performance under each scenario at the elevation band scale. For this, we took one sub-watershed from each watershed and calculated the HRU area weighted daily SWE for Elevation Band 1 to Elevation Band 5 for comparison.

3.3. Evaluating the Impact of Parameter Equifinality on Model Predictions

Different combinations of parameter values might lead to a similar level of fit between simulation and observation, a modeling phenomenon known as parameter equifinality [68,69]. In the S3 scenario, although the two-stage calibration strategy took both streamflow and SWE references into account, the second stage of calibration only started with the best set of parameter values derived from the first stage. Therefore, there might exist alternative parameter combinations that could result in similar or even better model predictions. Since it is difficult to manually test all the possible parameter combinations in the parameter space, the parameter equifinality could be investigated by evaluating the first few sets of best parameter combinations from the stage 1 calibration in S3. In this study, the values of the watershed parameters and elevation lapse rate parameters from the first 10 sets of the best parameters (listed in

Table A1. The 1st to 10th best parameter sets (rank1–rank10 correspond to the 10 sets of highest NSEs) from the last 500 simulations of the stage 1 calibrations in S3 for each watershed.

Kings Watershed
Parameter	1	2	3	4	5	6	7	8	9	10
V_SFTMP	1.57	3.23	4.46	2.48	4.51	3.94	2.03	3.19	0.60	2.15
V_SMTMP	2.76	1.72	2.74	1.21	4.18	3.24	3.85	4.85	4.78	4.26
V_SMFMX	2.79	2.41	1.33	2.97	2.21	4.19	2.01	3.43	5.89	4.85
V_SMFMN	0.63	1.57	1.03	0.03	2.75	1.95	0.89	2.71	1.51	5.01
V_TIMP	0.97	0.82	0.96	0.82	0.86	0.90	0.85	0.89	0.81	0.89
V_SNOCOVMX	444.5	498.5	129.5	344.5	132.5	36.5	345.5	25.5	29.5	56.5
V_SNO50COV	0.22	0.03	0.71	0.19	0.36	0.68	0.34	0.67	0.93	0.59
r_HRU_SLP	0.21	0.82	0.32	−0.41	0.51	0.53	−0.31	0.74	0.47	0.70
r_CN2	0.48	−0.39	0.39	0.26	0.02	−0.42	0.18	−0.16	−0.16	0.04
V_ALPHA_BF	0.92	0.98	0.95	0.05	0.68	0.21	0.79	0.02	0.76	0.94
V_GW_DELAY	97.5	437.5	320.5	478.5	26.5	259.5	81.5	54.5	330.5	486.5
r_SOL_AWC	−0.36	−0.42	0.08	−0.12	−0.34	0.39	−0.30	0.29	−0.36	0.33
V_GWQMN	4175	3575	3555	3065	4365	2715	1295	1145	4045	385
V_RCHRG_DP	0.20	0.68	0.82	0.57	0.65	0.63	0.61	0.48	0.68	0.28
V_TLAPSE	−7.91	−8.87	−1.57	−7.79	−2.99	−8.51	−2.21	−5.11	−8.79	−9.91
V_PLAPSE	16.50	88.70	58.10	48.5	17.30	44.5	98.30	95.90	5.3	39.5
NSE	0.886	0.866	0.864	0.839	0.833	0.831	0.828	0.826	0.817	0.812
Kern Watershed
Parameter	1	2	3	4	5	6	7	8	9	10
V_SFTMP	4.88	4.19	4.86	4.55	3.64	4.87	4.73	4.20	4.74	3.77
V_SMTMP	1.14	3.67	3.94	4.49	3.47	3.74	4.04	3.51	3.59	3.76
V_SMFMX	7.75	6.97	7.91	7.84	7.57	7.63	7.68	7.66	7.63	7.74
V_SMFMN	1.14	0.69	0.93	1.46	1.28	0.74	0.99	0.93	1.49	1.54
V_TIMP	0.92	0.95	0.93	0.91	0.92	0.90	0.95	0.95	0.91	0.91
V_SNOCOVMX	479.89	399.11	383.90	359.90	357.20	386.26	440.01	450.82	381.53	482.93
V_SNO50COV	0.45	0.45	0.39	0.45	0.40	0.46	0.33	0.45	0.29	0.41
r_HRU_SLP	−0.83	−0.69	−0.74	−0.82	−0.80	−0.84	−0.78	−0.71	−0.77	−0.81
r_CN2	−0.11	−0.18	−0.08	−0.01	−0.05	−0.08	−0.15	−0.17	−0.15	0.00
V_ALPHA_BF	0.06	0.05	0.05	0.05	0.05	0.15	0.09	0.05	0.08	0.05
V_GW_DELAY	168.98	254.04	313.07	278.35	319.57	245.79	145.98	118.20	222.36	241.89
r_SOL_AWC	−0.01	−0.11	−0.18	−0.11	0.05	−0.10	0.06	−0.14	0.04	−0.08
V_GWQMN	3469	3710.12	4295.67	4673.77	3937.65	4948.14	4493.08	3974.45	3338.71	3994.53
V_GW_REVAP	0.10	0.10	0.08	0.09	0.12	0.08	0.10	0.11	0.08	0.08
V_RCHRG_DP	0.44	0.59	0.56	0.58	0.44	0.52	0.54	0.37	0.60	0.43
V_TLAPSE	−6.54	−6.52	−8.28	−6.37	−7.59	−7.42	−7.30	−6.86	−6.74	−7.76
V_PLAPSE	75.60	73.30	70.93	54.10	64.53	73.02	56.88	78.31	49.92	75.46
NSE	0.815	0.800	0.798	0.795	0.789	0.788	0.787	0.784	0.779	0.778
Tule Watershed
Parameter	1	2	3	4	5	6	7	8	9	10
V_SFTMP	4.52	3.91	4.37	4.36	4.23	4.28	3.64	3.77	3.80	4.06
V_SMTMP	1.98	1.85	1.21	3.67	3.11	2.90	3.70	1.15	2.98	3.39
V_SMFMX	5.09	6.30	5.33	4.27	5.15	5.40	2.70	5.70	6.39	4.44
V_SMFMN	5.69	5.04	4.72	5.51	4.40	5.22	3.00	4.68	4.49	6.38
V_TIMP	0.93	0.87	0.91	0.88	0.86	0.92	0.92	0.84	0.85	0.85
V_SNOCOVMX	336.35	354.82	417.06	342.85	391.42	368.16	375.00	408.17	448.87	340.46
V_SNO50COV	0.65	0.82	0.68	0.79	0.55	0.66	0.64	0.49	0.42	0.67
r_HRU_SLP	−0.79	−0.84	−0.81	−0.80	−0.78	−0.79	−0.82	−0.83	−0.84	−0.83
r_CN2	−0.49	−0.33	−0.39	−0.48	−0.43	−0.48	−0.48	−0.46	−0.48	−0.24
V_ALPHA_BF	0.53	0.33	0.43	0.44	0.38	0.27	0.45	0.51	0.55	0.56
V_GW_DELAY	95.10	71.88	181.54	103	117.82	61.50	186.98	108.93	216.62	56.07
r_SOL_AWC	0.22	0.13	0.36	0.09	0.29	0.13	0.18	0.17	0.11	0.34
V_GWQMN	931.16	1431.35	261.26	288.06	1150	305.92	372.91	1918.15	792.72	1525.14
V_GW_REVAP	0.07	0.07	0.07	0.06	0.06	0.08	0.08	0.07	0.06	0.05
V_RCHRG_DP	0.07	0.11	0.03	0.27	0.079	0.27	0.04	0.13	0.002	0.14
V_TLAPSE	−6.94	−6.30	−4.69	−2.81	−6.06	−5	−3.89	−4.37	−4.92	−6.45
V_PLAPSE	27.32	32.91	29.75	27.02	18.75	23.06	15.52	28.17	19.35	17.35
NSE	0.808	0.785	0.774	0.773	0.771	0.759	0.758	0.754	0.749	0.747
Kaweah Watershed
Parameter	1	2	3	4	5	6	7	8	9	10
V_SFTMP	3.66	4.37	4.27642	3.52	3.46	2.74	4.12	4.70	3.23	3.45
V_SMTMP	2.26	2.45	0.85	1.07	0.89	0.97	1.02	1.61	1.12	2.35
V_SMFMX	3.81	4.82	3.01	3.87	5.27	4.66	3.18	2.05	1.95	5.61
V_SMFMN	0.70	1.93	0.22	0.79	0.46	0.50	1.09	1.10	1.08	0.11
V_TIMP	0.96	0.97	0.97	0.89	0.93	0.90	0.99	0.93	0.89	0.92
V_SNOCOVMX	330.56	498.60	483.98	419.35	259.75	267.05	272.67	323.25	392.94	378.33
V_SNO50COV	0.66	0.41	0.78	0.48	0.78	0.55	0.32	0.82	0.59	0.76
r_HRU_SLP	−0.81	−0.75	−0.70	−0.79	−0.77	−0.83	−0.82	−0.69	−0.80	−0.70
r_CN2	−0.34	−0.47	−0.19	−0.22	−0.27	−0.48	−0.36	−0.26	−0.28	−0.30
V_ALPHA_BF	0.58	0.86	0.65	0.56	0.83	0.73	0.48	0.87	0.93	0.58
V_GW_DELAY	104.32	291.52	134.88	42.57	13.33	181.67	173.88	73.13	20.47	147.88
r_SOL_AWC	−0.04	0.19	−0.02	−0.15	−0.30	−0.15	−0.17	0.15	0.19	−0.14
V_GWQMN	2655	2391.41	2561.47	2259.74	4174.36	3746.45	4103.04	4733.93	4212.76	2402.38
V_GW_REVAP	0.09	0.10	0.09	0.17	0.10	0.11	0.16	0.15	0.17	0.16
V_RCHRG_DP	0.40	0.58	0.58	0.47	0.42	0.42	0.36	0.44	0.22	0.33
V_TLAPSE	−2.36	−6.13	−3.393	−1.79	−3.74	−5.30	−4.55	−3.59	−2.73	−4.13
V_PLAPSE	55.90	6.06	37.41	38.16	15.05	32.42	14.80	31.92	5.18	23.30
NSE	0.697	0.672	0.657	0.655	0.648	0.645	0.645	0.625	0.621	0.619

in the Appendix A for the four upland watersheds) were used to substitute the corresponding values in the final calibrated parameter values from the S3 scenario. The modified 10 sets of parameters were applied for the SWAT simulation, respectively, for the calibration period 2001–2007 to test the equifinality effect on the streamflow and SWE, respectively.

4. Results

4.1. Impact of Three Model Calibration Scenarios on the Streamflow Predictions

Figure 3 and Figure 4 show the comparison between the simulated and observed streamflows for the three modeling scenarios (S1–S3) for the calibration and validation periods, respectively. The model performance is evaluated with NSE and R² values. The calibrated parameter values are shown in Table 1. The results indicate that all three calibration strategies generally reached a satisfactory level of streamflow prediction accuracy. For all three calibration scenarios, the NSE ranged from 0.4–0.9 and R² values ranged from 0.5–0.9 for all four watersheds. When comparing the two single-objective calibration scenarios, S1 and S2, the streamflow predictions performed slightly worse under the S2 than S1 scenario (Figure 3). This indicates that including snow parameters in the streamflow calibration process increases the model complexity, which might have a negative effect on the streamflow prediction accuracy and might reduce the prediction reliability. The results further show that the NSEs from S3 scenario for the Kern and Tule watersheds are 0.81 and 0.78, respectively, which are slightly higher than the values obtained for the S1 (0.80 and 0.77) and S2 scenarios (0.78 and 0.56) for both watersheds. However, for the Kings and Kaweah watersheds, S3 resulted in a slightly worse model performance, with NSEs of 0.85 and 0.67 compared to S1 (NSEs were 0.93 and 0.76) and S2 (NSEs were 0.89 and 0.74). This might be because the Kings and Kaweah watersheds have higher elevations and are located at the relatively northern side of the study region, and more snow is accumulated there. Therefore, a longer calibration period (1981–2007, from S1 and S2) would better capture the hydrological processes. Moreover, no matter which calibration strategy is used, when applying the calibrated model to the validation period (2008–2013), the streamflow prediction accuracy slightly decreased compared to the calibration period. For example, for the Kern watershed, the sub-watershed NSEs ranged between 0.78–0.81 in the calibration period and between 0.60–0.65 in the validation period (Figure 3). The decrease in model performance can, in part, be explained by the nonstationary nature of the streamflow under variable climate forcing and the fact that model calibration using a specific period from the historic record might not always reflect the same hydro-climatic conditions as observed during the validation period, as shown in Figure 4. Figure 4 further indicates that the calibration of snow parameters or involving the snow reference dataset could not improve the streamflow prediction accuracy. In general, the predicted streamflow performed better in the S1 modeling scenario in the Kings (R² = 0.93) and Kaweah (R² = 0.76) watersheds or reached a similar level of fit in all three calibration scenarios in the Kern (R² = 0.80) and Tule (R² = 0.78) watersheds. The insets in Figure 3(A3–D3) show the predicted streamflow extremes (10%, 20%, 80%, and 90% quantiles) versus observed data (black dashed line) for the 1983–2013 modeling period. The results indicate that, in general, the S1 scenario performs the best for predicting high-magnitude events among the three calibration scenarios (i.e., Kings, Tule, and Kaweah). However, for the low flow conditions (10% and 20% quantiles), there are no clear differences between the three calibration scenarios and watersheds. In contrast, Figure 3 points out that, in California, a place frequently affected by hydrological extremes, applying the default SWAT snow parameter values and calibrating the model with only streamflow reference data provides the best agreement between the observed and predicted extreme streamflow values. However, when predicting a low flow, any of the three calibration strategies could be applied, resulting in very similar low-flow predictions.

4.2. Impact of the Three Model Calibration Scenarios on the SWE Predictions

The simulated SWE from the SWAT output is in the form of a snow water-equivalent value at the HRU level for each elevation band. The HRU area weighted SWE for each sub-watershed is further calculated for all four watersheds, respectively. The sub-watershed scale-calculated SWE from S1–S3 are compared with ParBal-referenced SWE data. In order to better represent the findings, we separated the snow-dominated sub-watersheds into two groups: calibrated sub-watersheds and poor-performing sub-watersheds. The calibrated sub-watersheds (marked as blue in Figure 5) are those NSEs of the SWE predictions that have reached or are close to the acceptable level, and they will be included in the following analysis and comparison, whereas the poor-performing sub-watersheds (marked as red in Figure 5) are those SWE predictions from all S1–S3 modeling scenarios that are very poor, even after the two-stage calibration (i.e., NSE<0) and far away from acceptable NSE values. The poor-performing sub-watersheds are excluded from the following analysis. The sub-watersheds need to be separated in this way, because the lumped calibration strategy was applied in S1–S3, and the parameter values found for the entire watershed might not be suitable for each individual sub-watershed simultaneously. Figure 5 indicates that the S3 calibration scenario could improve the overall SWE predictions of the calibrated sub-watersheds. Although there exist cases that the maximum values of NSE are not the highest after the snow calibrations of S3 (i.e., the Tule watershed; see Figure 5A), the variations of NSEs from the calibrated sub-watersheds are smaller (i.e., the range of NSEs is about 0–0.62 after snow calibration for the Tule watershed, whereas the ranges are −0.18–0.63 before snow calibration). Figure 5B shows the average NSE values for the snow-calibrated sub-watershed (shown in Figure 5C); the results indicate that the average NSE improved for S3 compared with the other two scenarios. Therefore, from the overall perspective, the SWE prediction is more reliable when applying the S3 calibration scenario compared with the other two calibration strategies. Figure 6 compares the SWE predictions at the elevation band scale among the three calibration scenarios for one typical snow-dominated sub-watershed in each watershed with reference snow data. The typical sub-watersheds selected from each watershed are those having high SWE accumulation during 2001–2013 based on the observation data. The Kings watershed and Kaweah watershed have relatively higher SWE, and the peak SWE values are more than 1000 mm, while Kern and Tule watersheds have smaller SWE values, the peaks of which are around 500 mm. This might be caused by the relative locations of the four watersheds (Kings and Kaweah are located relatively north and have higher elevations). The NSE of the sub-watershed-level snow NSEs are summarized in

Table 4. SWE Nash-Sutcliff Efficiencies for all three scenarios for the typical sub-watersheds in each watershed.

NSE	S1	S2	S3
Kings	0.5708	0.5300	0.6488
Kern	0.6900	0.6799	0.7204
Tule	0.6729	0.6554	0.6896
Kaweah	0.6932	0.6590	0.6478

S1–S3 refer to Scenario 1-Scenario 3.

. It indicates that the S3 SWE prediction performance is better and more reliable among the three modeling scenarios for the Kings, Kern, and Tule watersheds. Additionally, S2 performs slightly worse than S1 for all four watersheds. This points out that involving the snow parameter alone without the SWE reference data in the SWAT model calibration would only increase the model parameter complexity, which would not improve the SWE prediction accuracy. Moreover, the second row to sixth rows show the SWE predictions for elevation bands 1–5, respectively, for each typical sub-watershed. The results point out that SWE predictions mainly deviate from each other at low elevations among the three calibration scenarios, whereas, at high elevations, the three calibration strategies show similar levels of SWE predications. Figure 6 indicates that the S3 scenario is better at capturing SWE at low elevations, which makes its predictions more accurate and reliable than S1 and S2.

4.3. Influence of Parameter Equifinality on the Prediction Reliability

The calibrated snow parameters from S3 are applied, and the model equifinality was evaluated through the changing of the watershed and elevation lapse rate parameters in this study. Figure 7 and Figure 8 show the simulated streamflow and SWE based on the 1st–10th best sets of parameter values from calibration stage 1 in the S3 scenario for 2001–2007, respectively (previously, the S3 modeling scenario only considered the 1st best parameter set from 500 simulations from stage 1). The final snow parameter values derived from S3 scenario stage 2 were applied unchanged for all ten simulations. Figure 7 and Figure 8 indicate that both the streamflow and SWE predictions based on the best parameter set from the stage 1 calibrations were not eventually the best among the ten sets. The streamflow predictions were more significantly affected by the equifinality issue, since the NSEs changed 0.2–0.4 from their smallest value to the largest value (i.e., the NSE values of Kings varied between 0.55–0.9; the change was about 0.35), given the ten sets that yielded almost similar levels of prediction in the stage 1 streamflow calibration period. Comparatively, the SWE predictions were also affected by the equifinality issue; however, the influence was smaller than the streamflow predictions. The medium of the NSE values varied less than 0.1 in most cases (i.e., it ranged 0.6–0.69 for the Kings watershed).

5. Discussions

5.1. Applying ParBal-Reconstructed SWE Data in the Calibration Process

In general, this study shows that the SWAT calibration is strongly influenced by snow processes in the upland watersheds by matching the flow peak and smoothing the hydrograph related to snowmelt [70,71]. As expected, the results from this study prove that incorporating SWE reference data in the model calibration process would improve the model SWE prediction accuracy in general and enhance the overall prediction reliability. Despite that some previous studies applied remote-sensing (i.e., MODIS)/in-situ-based snow cover reference datasets to validate the SWAT SWE predictions, they seldom directly included SWE reference information into the model calibration process as a means to further improve the results [72,73]. Although the multi-objective calibration of the SWAT model was not firstly introduced in this study (previously, researches applied second-reference datasets, such as soil moisture, evapotranspiration, or sediment loads, in addition to streamflows, to demonstrate the advantage of multi-objective calibration in improving the model performance [40,42,74]), SWE reference data is seldom considered, since in-situ snow records are lacking for many upland watersheds. Ground-based snow observation stations such as SNOTEL or other in-situ snow records could be applied to calibrate and validate the snow predictions. One recent study uses the SWE information derived from in-situ snow depth observations in the Upper Adige River watershed (6875 km², a similar size to the Kern watershed) in the Northeastern Italian Alps for SWAT simulations [75]. Nevertheless, that study required sub-watershed levels and elevation band level snow depth observation stations densely distributed in the study area. However, to date, not many upland watersheds offer such highly dense in-situ SWE observation networks. In addition, they used a short SWE calibration and validation period (two years for each), which may not be able to capture both the wet/dry cycles. Comparatively, in our study, seven years were applied for the SWE calibrations, and another six years were used for the validation. Both wet and dry years were in the two periods. The remote-sensing-based SWE reconstruction dataset (such as ParBal) has wider applications as reconstructed SWE products become available for more mountainous regions. The two-stage calibration approach used in this study could be easily applied to other ungauged watersheds.

5.2. The Role of Snow Parameters in SWAT Model Predictions

The findings from this study demonstrated the importance of snow parameters in the hydrological model prediction accuracy. Comparing S1 and S2 in assigning snow parameter values to SWAT, the results showed that, in general, the default snow parameter values resulted in a similar performance compared with the calibrated snow parameters when only the streamflow reference data was applied. However, when comparing these two calibration strategies with S3, it pointed out that calibrating the snow parameters with the snow reference data would best benefit the SWE prediction accuracy. This indicates that involving snow parameter calibrations (instead of applying the default parameter values) would only introduce model parameter complexity, and it would not be useful unless the corresponding SWE reference dataset was used in the calibration process. One former study pointed out that the major application of default snow parameter values is to assess the applicability of the SWAT in ungauged watersheds, and there is little observed data available for calibrations [76]. This partly aligned with the findings from our study; however, we also found that default parameters might also be useful for calculating extreme flows, especially the high flow conditions for the upland watersheds. S1 behaves the best when simulating the 80% and 90% quantile high flows, and this broadened the application of the default snow parameters of the SWAT in identifying hydrological responses under extreme high-flow situations (such as flooding).

5.3. The Impact of Parameter Equifinality on the Model Prediction

The identification of parameter values is a major concern in calibrating the model for a specific watershed [36,77]. The results from this study point out that the calibration of snow parameters with corresponding SWE reference datasets definitely improves the SWE prediction reliability in general and reduces the equifinality issue to a certain level. However, this does not indicate that the final calibrated parameters got rid of the equifinality issues. To translate this in another way, the calibrated parameters just control the equifinality issue within our satisfactory range for the two dimensions (flow and SWE) in this study. The results indicate that, on the one hand, the more types of parameters and reference datasets (i.e., evapotranspiration and soil moisture) involved in the calibration process, the higher the chance of a calibrated model closer to reality (less of an equifinality issue). On the other hand, the model parameters calibration process also plays the role of compensating the overall uncertainty from different components in the modeling process to the reference dataset. For instance, the uncertainty comes from the simplification of the model structure and governing equations, as well as input climate data uncertainty, among others [16,78,79]. The calibrated parameters should satisfy the goal for the targeted variables predicted by the model. Therefore, we suggest the SWAT users apply the appropriate number of calibration parameters and reference datasets regarding the prediction variables that they are interested in instead of just simply involving as many parameter as possible and aiming to make it the “true model”.

6. Summary and Conclusions

In this study, SWAT models were implemented for four uplands watersheds (Kings, Kern, Tule, and Kaweah) in Tulare Lake Basin. Three model calibration scenarios were applied in four upland watersheds located in the southern part of the Sierra Nevada Mountains using single (monthly streamflow at the watershed outlet) or multi-objective (monthly streamflow and daily ParBal-reconstructed SWE at the sub-watershed scale) calibrations under the lumped calibration scheme. The extreme flow conditions were also evaluated. Additionally, under the multi-objective calibration scheme, given the fixed calibrated snow parameters, the influence of the watershed parameters and elevation lapse rate parameters on the model parameter equifinality was investigated. Specifically, the following conclusions were drawn from this study:

• Although all three calibration strategies provide satisfactory streamflow predictions in general, slight differences exist among the three modeling scenarios, and the S3 calibration scenario does not always perform the best for streamflow prediction accuracy. In general, S2 performs slightly worse than the S1 scenario, which indicates that involving model parameter complexity might have a negative impact on the streamflow predictions if the model is calibrated only with a streamflow.
• The S3 calibration scenario shows its advantages for SWE predictions compared with S1 and S2. Even though, in some cases, the maximum NSE values of the sub-watersheds are not the highest after snow calibration, the variations are smaller. Therefore, the two-stage calibration strategy is recommended under the lumped calibration scheme for the upland watersheds, where snow plays a significant role.
• The S1 modeling scenario with default parameter values provides the best prediction for high flow conditions (related to floods), whereas no obvious differences exist among the three scenarios in simulating low flow conditions (droughts).
• All three scenarios work well for capturing large amounts of snow in the high elevations; however, the S3 scenario performs better than the other two scenarios in capturing small amounts of SWE in low elevations.
• The streamflow predictions are affected by the equifinality issue, and the changes of the NSE values vary between 0.2–0.3 for the ten parameter sets. Comparatively, the SWE prediction is also affected by the equifinality issue; however, the influence is slighter; the variation of the NSE is less than 0.1 for most cases.
• SWAT users are suggested to apply the appropriate number of calibration parameters and reference datasets for the goal (i.e., SWE and streamflow, in this study) that they are interested in. In most cases, there is no need to involve extra parameter complexity that is not related to the goal, and it might have a negative impact on the model predictions.

There exist certain limitations for this study. Since a lumped calibration strategy is applied for the SWE calibrations, it is difficult to make all the sub-watersheds reach a good fit simultaneously, which is caused by the heterogeneity natures of the different sub-watersheds. Additionally, although the model parameter equifinality issue is investigated, it is only partially tested for the first 10 best parameter sets from calibration stage 1 because of the large amounts of simulation work. In the future, the two-stage calibration strategy could be applied to more snow-dominated watersheds to strengthen the findings from this study.