Improvement of the SWAT Model for Snowmelt Runoff Simulation in Seasonal Snowmelt Area Using Remote Sensing Data

Abstract

The SWAT model has been widely used to simulate snowmelt runoff in cold regions thanks to its ability of representing the effects of snowmelt and permafrost on runoff generation and confluence. However, a core method used in the SWAT model, the temperature index method, assumes both the dates for maximum and minimum snowmelt factors and the snowmelt temperature threshold, which leads to inaccuracies in simulating snowmelt runoff in seasonal snowmelt regions. In this paper, we present the development and application of an improved temperature index method for SWAT (SWAT+) in simulating the daily snowmelt runoff in a seasonal snowmelt area of Northeast China. The improvements include the introduction of total radiation to the temperature index method, modification of the snowmelt factor seasonal variation formula, and changing the snowmelt temperature threshold according to the snow depth derived from passive microwave remote sensing data and temperature in the seasonal snowmelt area. Further, the SWAT+ model is applied to study climate change impact on future snowmelt runoff (2025–2054) under the climate change scenarios including SSP2.6, SSP4.5, and SSP8.5. Much improved snowmelt runoff simulation is obtained as a result, supported by several metrics, such as MAE, RE, RMSE, R², and NSE for both the calibration and validation. Compared with the baseline period (1980–2019), the March–April ensemble average snowmelt runoff is shown to decrease under the SSP2.6, SSP4.5, and SSP8.5 scenario during 2025–2054. This study provides a valuable insight into the efficient development and utilization of spring water resources in seasonal snowmelt areas.

1. Introduction

Snowmelt runoff is an essential component of water resources in many parts of the world [1]. At the same time, it may become a disaster in the form of a snowmelt flood under certain conditions [2]. Accurate simulation and prediction of snowmelt runoff are crucial in water resources management and snowmelt flood disaster prevention [3,4]. Naturally, snowmelt runoff is affected by a variety of meteorological factors and involves many complex physical processes [5]; as such, modeling snowmelt runoff is more challenging than that for conventional rainfall-runoff processes [6,7]. These challenges become even more exacerbated when it comes to representing the climate change impact on hydrological processes [8,9].

Many researchers have attempted to develop snowmelt runoff models directly or to add snowmelt modules to existing hydrological models in order to properly simulate snowmelt runoff, including, for example, HBV [10,11], SRM [12,13], UBC [14,15], CRHM [16,17], ARHYTHM [18], and SWAT model [19,20], to name just a few. Among them, the soil and water assessment tool (SWAT) model has been widely used for both simulating and predicting snowmelt runoff, as it considers the impact of snowmelt and frozen soil on runoff generation and concentration [21,22,23].

The SWAT model has a dedicated snowmelt module, which considers only the impact of snowpack temperature and air temperature on snowmelt [24,25]. This often leads to underestimating snowmelt runoff [26]. Some attempts have been made to address this problem, for example, by integrating an energy balance snowmelt model into the SWAT model [27,28]. A clear drawback of this approach is that it often requires extra large amounts of data for parameterization and, when such requirement cannot be met, results from the energy balance model can become worse than other simpler methods [29,30]. Meanwhile, others tried to improve the temperature index method, for example, by including the radiation factor [31,32,33] on the ground that radiation is one of the main energy sources of snow melting [34,35,36]. It is further argued that solar radiation is fundamentally more influential compared with the temperature, as part of the solar radiation energy is used to warm the snow up and the rest to melt the snow [37]. In other words, the snowmelt process is under joint influences of both temperature and solar radiation [35]. From this perspective, research added direct solar radiation to the temperature index method to simulate snowmelt runoff and achieved good results [38]. Apart from direct solar radiation, scattering radiation is also important for snowmelt [39]. Accordingly, adding total radiation on the snow surface into the SWAT model is worth trying to mitigate the snowmelt runoff underestimation problem.

Another perspective from which the snowmelt module of the SWAT model can be improved is related to an important parameter, i.e., the snowmelt factor, which varies seasonally [40]. In the existing setting, the seasonal variation in the snowmelt factor is approximated using a sinusoidal function, which assumes the minima and maxima to be on the 21 December and 21 June, respectively [41]. Clearly, this is not an even close approximation for snowmelt, as snowmelt mainly occurs in spring in many seasonal snowmelt areas [42,43]. In addition, the snowmelt factor decreases as the snow depth reduces [44] and, consequently, the maximum snowmelt rate in the seasonal snowmelt area in the northern hemisphere occurs in late spring [45]. In other words, the assumption with respect to the snowmelt factor is not well supported by the actual processes, which may be a key factor for the low accuracy of snowmelt runoff simulations in the SWAT model. Some researchers proposed using discharge for snowmelt factor determination [46]. However, there is a large time difference between the snowmelt occurrence and snowmelt runoff occurrence in mountainous areas with long confluence time [47,48]. Therefore, in this study, we chose to use the snowmelt occurrence time directly so as to avoid the interference caused by long confluence time.

In addition to the snowmelt factor, another important parameter in the snowmelt module is the snowmelt temperature threshold whose default value is set as 0 °C [49]. Some studies suggested that only 0 °C is not necessarily a best threshold value to determine whether the snow melts. In addition, they advised using the so-called accumulated temperature, which, instead of using a single observed temperature value, determines whether snowmelt occurs by using the sum of all temperature readings in a day that are greater than 0 °C [38,46,50,51]. However, the study ignored the possibility that snowmelt could also occur below 0 °C [52]. The threshold of snowmelt temperature is also closely related to the geographical environment. For example, the snowmelt temperature threshold is found to be below 0 °C in Altay Mountains [53] but well above 0 °C in central Switzerland [54]. Therefore, accurately determining the snowmelt temperature threshold by considering the characteristics of the regional geographical environment needs to be addressed.

Although previous studies have noticed the deficiencies in the existing snowmelt module of SWAT in terms of no consideration of solar radiation, inaccurate snowmelt temperature threshold determination, and no representation of the regional variation in snowmelt factor [38,46,50,51], no attempts have been made to address them in a comprehensive manner. In the study presented in this paper, we systematically addressed these issues by (i) integrating both direct solar radiation and scattering radiation into the temperature index; (ii) using local observations to determine the dates of the maxima and minima of snowmelt factors and modify accordingly the snowmelt factor seasonal variation formula; (iii) determining the temperature threshold of snow melting by identifying the relationship between snow depth and temperature and changing the original temperature threshold; and, (iv) finally, applying the modified model in predicting snowmelt runoff in the study area—the Baishan basin. To facilitate discussion, the modified SWAT model is hereafter referred to as SWAT+.

This paper is structured as follows: Section 2 describes the study area and the data used; Section 3 presents the overall methodology and discusses the details of the improvements made towards the snowmelt module and the bias correction measures in utilizing climate projection data. The results are explained and discussed in Section 4 and Section 5. Finally, the conclusions drawn from this study are presented in Section 6.

2. Materials

2.1. Study Area

Located in the northeast of China (126°28′–128°51′ E, 41°42′–43°21′ N), the Baishan basin is a typical seasonal snowmelt area (Figure 1). The Baishan basin is the main source area of the Second Songhua River, which originates from the Changbaishan Tianchi and consists of two major tributaries: the Toudao Songhua river and the Erdao Songhua river, with a total drainage area of 18,645.4 km². Snowmelt runoff is a major contribution to the spring water resources in the Baishan basin, lasting from March to April annually [55]. Spring floods often occur when there is a continuing and rapid rise in temperature in spring. According to the long time series of daily snow depth derived from passive microwave remote sensing data provided by the National Qinghai Tibet Plateau data center [56,57], the timing of the first stable snowfalls in the basin ranges from 5 November to 11 November over the 40-year period of 1980–2019, with an averaged starting time of 9 November. The snow cover can last up to 6 months, with the maximum daily snowfall depth up to 24.64 mm. The main snowmelt start time is from 2 March to 14 March.

2.2. Data

Several types of data are used in this study, including spatial data, historical meteorological and hydrological observations, and climate projections. The spatial data consist of digital elevation model (DEM) (obtained from http://www.gscloud.cn website, accessed on 22 December 2021), soil map, and land use. The soil map and land use data come with a spatial resolution of 1 km and are collected from the harmonized world soil database (version 1.1) (http://www.fao.org/soils-portal/data-hub/soil-maps-and-databases/harmonized-world-soil-database-v12, accessed on 25 December 2021) and resource and the environment science data center website (https://www.resdc.cn/, accessed on 15 January 2022), respectively. The historical meteorological and hydrological observations used in the study include daily precipitation, air temperature, wind speed, relative humidity, total solar radiation, the snow depth at three gauged meteorological stations (Jingyu, Donggang, and Songjiang), and discharge from the Baishan reservoir from 1980 to 2019. They were obtained from the National Meteorological Science Data Center (https://data.cma.cn, accessed on 15 January 2022), the National Qinghai Tibet Plateau data center (http://data.tpdc.ac.cn, accessed on 5 April 2022), and Baishan hydropower plant of State Grid Xinyuan company, respectively. Among them, the snow depth is inverted using passive microwave remote sensing, which has also been used in much climate and hydrological research in China [58,59,60] and proved well applicable in Northeast China [61].

Climate projection data from the Coupled Model Intercomparison Project (CMIP6) (http://cmip-pcmdi.llnl.gov/cmip6/, accessed on 3 July 2022), i.e., the BCC-CSM2-MR (Beijing Climate Center/China) model, are selected for future snowmelt runoff changes analysis. The baseline period is set as 1980–2019 and the future period is 2025–2054. There are three different climate change scenarios considered, i.e., SSP2.6, SSP4.5, and SSP8.5, representing sustainable development, moderate development, and conventional development, with the artificial radiation force reaching 2.6, 4.5, and 8.5 W/m², respectively [62,63].

3. Methods

3.1. Modification and Improvement of the SWAT Model

This section presents the modification made to the original SWAT model, aiming to improve its snowmelt runoff simulation performance. The changes are based on the revision 682 of the SWAT source code (https://bitbucket.org/blacklandgrasslandmodels/swat_development/src/master, accessed on 3 March 2022). The modifications focus on the three components of the snowmelt module in SWAT, i.e., snowmelt calculation taking into account the contribution from radiation, snowmelt factor with more accurate seasonal representation, and snowmelt threshold temperature determined by the local environment.

3.1.1. Snowmelt Calculation Considering Total Solar Radiation

In the original SWAT snowmelt module, the change in snowpack depends on sublimation, melting, and snowfall, which is described as the snow water equivalent (SWE), representing water stored in the snowpack. The SWE is calculated as follows:

$${\mathrm{SWE}}_{\mathrm{i}}={\mathrm{SWE}}_{\mathrm{i}-1}+\mathrm{P}-\mathrm{E}-\mathrm{SM}$$

(1)

where ${\mathrm{SWE}}_{\mathrm{i}}$ and ${\mathrm{SWE}}_{\mathrm{i}-1}$ are the water content of the snowpack (mm) on the day and the day before; $\mathrm{P}$, $\mathrm{E}$, and $\mathrm{SM}$ are the snowfall (mm), sublimation (mm), and snowmelt (mm), respectively. Snowmelt ($\mathrm{SM}$) in the original SWAT model is calculated using the temperature index method, i.e.:

$$\mathrm{SM}=\mathrm{SMF}\times {\mathrm{sno}}_{\mathrm{cov}}\times \left[\frac{{\mathrm{T}}_{\mathrm{snow}}+{\mathrm{T}}_{\mathrm{mx}}}{2}-{\mathrm{T}}_{\mathrm{mlt}}\right]$$

(2)

where $\mathrm{SMF}$ is the snowmelt factor (mm H₂O d⁻¹·°C⁻¹). ${\mathrm{sno}}_{\mathrm{cov}}$ is the snow cover fraction. Some studies simulated and predicted snow cover fraction through evolutionary algorithm [64]. In SWAT model, ${\mathrm{sno}}_{\mathrm{cov}}$ is calculated based on water content of the snow pack using Equation (3). ${\mathrm{T}}_{\mathrm{snow}}$, ${\mathrm{T}}_{\mathrm{mx}}$, and ${\mathrm{T}}_{\mathrm{mlt}}$ is the snowpack temperature, maximum air temperature, and snowmelt temperature threshold (°C), respectively.

$${\mathrm{sno}}_{\mathrm{cov}}=\frac{\mathrm{SNO}}{{\mathrm{SNO}}_{100}}\times {\left[\frac{\mathrm{SNO}}{{\mathrm{SNO}}_{100}}+\exp \left({\mathrm{cov}}_1-{\mathrm{cov}}_2\times \frac{\mathrm{SNO}}{{\mathrm{SNO}}_{100}}\right)\right]}^{-1}$$

(3)

where $\mathrm{SNO}$ is the water content of the snow pack on a given day (mm H₂O), ${\mathrm{SNO}}_{100}$ is the threshold depth of snow at 100% coverage (mm H₂O), and ${\mathrm{cov}}_1$ and ${\mathrm{cov}}_2$ are coefficients that define the shape of the curve. The values used for ${\mathrm{cov}}_1$ and ${\mathrm{cov}}_2$ are determined by solving Equation (3) using two known points: 95% coverage at 95% ${\mathrm{SNO}}_{100}$ and 50% coverage at a user-specified fraction of ${\mathrm{SNO}}_{100}$.

Given the important role of solar radiation in the snowmelt process, some researchers have added solar radiation to the classical temperature index method [65,66] and reported improved accuracy in snowmelt simulations. We derived a similar approach in this study and introduced the total solar radiation to the original snowmelt module of the SWAT model. Now the snowmelt is calculated using a modified equation as shown in Equation (4):

$$\mathrm{SM}={\mathrm{sno}}_{\mathrm{cov}}\times \left(\mathrm{SMF}\times \left(\frac{{\mathrm{T}}_{\mathrm{snow}}{+\mathrm{T}}_{\mathrm{mx}}}{2}-{\mathrm{T}}_{\mathrm{mlt}}\right)+\mathrm{SR}\times \left(1-\mathsf{\alpha }\right)\times {\mathrm{m}}_{\mathrm{Q}}\right)$$

(4)

where SR is the incoming total solar radiation (W/m²), α is snow albedo, and ${\mathrm{m}}_{\mathrm{Q}}$ is the coefficient of energy to water depth conversion. Following a previous study [67], we chose to use α = 0.679, ${\mathrm{m}}_{\mathrm{Q}}$= 0.26 mm·W⁻¹·m²·day⁻¹.

3.1.2. Snowmelt Factor with an Improved Representation of Seasonal Variation

The snowmelt factor (SMF) in the snowmelt module of the SWAT model already considers the seasonal variation, as shown in Equation (5):

$$\mathrm{SMF}=\frac{{(\mathrm{SMF}}_6{+\mathrm{SMF}}_{12})}{2}+\frac{{(\mathrm{SMF}}_6-{\mathrm{SMF}}_{12})}{2}\times \sin (\frac{2\mathsf{\pi }}{365}\times (\mathrm{n}-81\left)\right)$$

(5)

where SMF₆ and SMF₁₂ represent the snowmelt factor on the 21 June (172nd day) and the snowmelt factor on 21 December (the 355th day in a year), respectively, that is, when n = 172 and n = 355, 2π/365 × (n − 81) reaches π/2 and 3π/2 phases, respectively. The sin(2π/365 × (n − 81)) formula accordingly reaches maximum values of 1 and minimum values of −1. In addition, n is the order position of the day in a year, i.e., the nth day. In other words, the maximum and the minimum snowmelt factors are assumed to happen on the two fixed dates, which is clearly unfit for seasonal snowmelt basins where solar radiation is an important element affecting the snowmelt rate [68].

In this study, instead of using the two fixed dates, we attempted to determine a more accurate representation of the seasonal variation in the snowmelt factor by, firstly, investigating the relationship between the 5-day moving average of snow depth and temperature from November to April. As shown in Figure 2, the snow depth decreases as the radiation gradually rises. The black arrows indicate the time when the minimum snowmelt factor occurs [69]. Conversely, when the snow depth gradually decreases, the time corresponding to the maximum radiation will be the time of the maximum snowmelt factor, as indicated by the green arrows in Figure 2.

Then, the multi-year average occurrence time for the maxima and minima of the snowmelt factor is obtained by examining the temporal distributions of the two dates, i.e., using a box plot. As shown in Figure 3A, the median date of the snowmelt factor minima is 8 March (the 67th day in a year) and that for the maximum snowmelt factor is 16 April (the 107th day) (Figure 3B). In this study, these two medians are used to represent the maximum and minimum snowmelt factors. Following the sinusoidal variation principle, the formula of the snowmelt factor is modified as Equation (6). The minimum sinusoidal value and maximum sinusoidal value occur when the values of 67 and 107 are substituted into the formula, respectively. That is, when n = 67, the sin function should be in the $3\mathsf{\pi }/2$ phase and the sinusoidal value should be −1. Meanwhile, the sin function is $\mathsf{\pi }/2$ phase and sinusoidal value is 1 for n = 107. Based on this, the formula of the snowmelt factor is modified as follows:

$$\mathrm{SMF}=\frac{\left({\mathrm{SMF}}_3+{\mathrm{SMF}}_4\right)}{2}+\frac{\left({\mathrm{SMF}}_4-{\mathrm{SMF}}_3\right)}{2}\times \sin (\frac{9\mathsf{\pi }}{365}\times \left(128-\mathrm{n}\right))$$

(6)

where ${\mathrm{SMF}}_{3 }$ is the snowmelt factor on 8 March and ${\mathrm{SMF}}_4$ is the snowmelt factor on 16 April, that is, the largest and minimum snowmelt factors in the seasonal snowmelt areas.

3.1.3. Snowmelt Temperature Threshold

The snowmelt temperature threshold is another critical parameter in the snowmelt calculation method described by Equation (2), which determines the time when snowmelt starts. In this study, we determine the value of this parameter by investigating the relationship between the 5-day moving averages of the snow depth and temperature from November to April in the next year. As shown in Figure 4, there is a clear negative correlation between the two quantities. In addition, the temperature from which the snow depth starts to fall sharply is determined as the snowmelt temperature threshold, as indicated by the green arrows.

The distribution of this parameter is plotted in Figure 5A, which shows that the snowmelt temperature threshold mainly concentrates around −2.38–0.40 °C from 1980 to 2019, with the annual average being −0.54 °C. Following this finding, we then modified the snowmelt temperature threshold in the model as −0.54 °C. The corresponding snow depth at the beginning of the snowmelt is 10.30–15.87cm (Figure 5B).

3.2. Setup and Validation of the Modified Model SWAT+

The setup of the modified model largely follows the same procedure of that for the original SWAT model, as the modifications are made only to the internal calculation components. Firstly, the DEM data (Figure 6A) were used to delineate the watersheds where the study area is divided into 31 watersheds (Figure 6B). Secondly, the hydrological response unit (HRU) is separated according to the terrain, land use, and soil types. Overall, there are 6 land use types (Figure 6C) and 17 soil types (Figure 6D). Before generating HRU, the spatial resolution of the land use, soil type, and DEM data is unified to 30 m × 30 m. In addition, the study area is classified into 227 HRUs. Finally, the meteorological data were imported and the initial model parameters were set.

The parameters of the model were then calibrated using the SUFI-2 algorithm provided by the standard Soil and Water Assessment Tool Calibration and Uncertainty Procedures (SWAT-CUP) component [70]. The SUFI-2 algorithm generates a group of parameters randomly by Latin-Hypercube simulations for seeking the best parameters combination [71]. In total, 28 parameters were chosen to be calibrated. The original range and fitted value of these parameters for both SWAT and SWAT+ are listed in

Table 1. Parameters for SWAT and SWAT+.

Parameter Name	Meaning	Max	Min	SWAT Fitted Value	SWAT+ Fitted Value
SFTMP	Snowfall temperature (°C)	5	−5	0.82	−0.07
ESCO	Soil evaporation compensation factor	1	0.01	0.30	0.65
EPCO	Plant uptake compensation factor	1	0.01	0.98	0.09
CN2	SCS runoff curve number for moisture condition II	0.5	−0.5	0.29	0.52
BIOMIX	Biological mixing efficient	0.5	−0.5	0.92	−0.21
SOL_Z	Depth from soil surface to bottom of layer (mm)	0.5	−0.5	0.22	−0.13
SOL_AWC	Available water capacity of the soil layer (mm H2O/mm soil)	0.5	−0.5	0.68	0.64
SOL_K	Saturated hydraulic conductivity (mm/h)	0.8	−0.8	1.13	0.09
SOL_ALB	Moist soil albedo	0.5	−0.5	0.01	−0.09
GW_DELAY	Groundwater delay time (days)	500	0	431.71	111.34
ALPHA_BF	Baseflow Alpha Factor	1	0	1.69	0.47
SMTMP	Snow melt base temperature (°C)	5	−5	0.00	−1.46
GWQMN	Threshold depth of water in the shallow aquifer required for return flow to occur (mm H2O)	5000	0	1380.81	7057.22
GW_REVAP	Groundwater “revap” coefficient	0.2	0.02	0.19	0.12
REVAPMN	Threshold depth of water in the shallow aquifer for “revap” or percolation to the deep aquifer to occur (mm H2O)	500	0	106.88	−21.94
CH_N2	Manning’s “n” value for the main channel	0.3	0	0.15	0.35
CH_K2	Effective hydraulic conductivity in main channel alluvium (mm/h)	150	0	145.12	−46.50
RCHRG_DP	Deep aquifer percolation fraction	1	0	0.95	0.42
SMFMX	Melt factor for snow on 21 June (mm H2O/°C-day)	10	0	0.36	7.73
SMFMN	Melt factor for snow on 21 Dec (mm H2O/°C-day)	10	0	0.58	6.31
TIMP	Snow pack temperature lag factor	1	0.01	0.50	0.13
SURLAG	Surface runoff lag coefficient	24	1	4.93	5.88
TLAPS	Temperature lapse rate (°C/km)	10	−10	4.80	−1.95
SLSUBBSN	Average slope length (m)	150	10	−7.74	15.78
CANMX	Maximum canopy storage (mm H2O)	100	0	−2.37	−11.74

.

3.3. Evaluation of the Modified Model SWAT+

We employed several indexes to compare the performance of SWAT+ against the original SWAT for snowmelt runoff simulation, including the root mean square error (RMSE), mean absolute error (MAE), relative error (RE), the coefficient of determination (R²), and Nash Sutcliffe coefficients (NSE) [72], as shown in Equations (7)–(11).

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{R}}_{\mathrm{obs},\mathrm{i}}-{\mathrm{R}}_{\mathrm{mod},\mathrm{i}})}^2}{\mathrm{n}}}$$

(7)

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|\left({\mathrm{R}}_{\mathrm{obs},\mathrm{i}}-{\mathrm{R}}_{\mathrm{mod},\mathrm{i}}\right)\right|$$

(8)

$$\mathrm{RE}=\left({\mathrm{R}}_{\mathrm{mod},\mathrm{i}}-{\mathrm{R}}_{\mathrm{obs},\mathrm{i}}\right)/{\mathrm{R}}_{\mathrm{obs},\mathrm{i}}\times 100\%$$

(9)

$${\mathrm{R}}^2=\frac{{\left[{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{R}}_{\mathrm{sim},\mathrm{i}}-\overline{{\mathrm{R}}_{\mathrm{sim}}}\right)\left({\mathrm{R}}_{\mathrm{obs},\mathrm{i}}-\overline{{\mathrm{R}}_{\mathrm{obs}}}\right)\right]}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{R}}_{\mathrm{sim},\mathrm{i}}-\overline{{\mathrm{R}}_{\mathrm{sim}}}\right)}^2{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{R}}_{\mathrm{obs},\mathrm{i}}-\overline{{\mathrm{R}}_{\mathrm{obs}}}\right)}^2}$$

(10)

$$\mathrm{NSE}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{R}}_{\mathrm{obs},\mathrm{i}}-{\mathrm{R}}_{\mathrm{mod},\mathrm{i}})}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{R}}_{\mathrm{obs},\mathrm{i}}-\overline{{\mathrm{R}}_{\mathrm{obs},\mathrm{i}}})}^2}$$

(11)

where ${\mathrm{R}}_{\mathrm{obs},\mathrm{i}}$ and ${\mathrm{R}}_{\mathrm{mod},\mathrm{i}}$ are the observed snowmelt runoff and simulated snowmelt runoff and $\mathrm{n}$ is the total number of time steps (days) of the simulations.

3.4. Bias Correction Method for Climate Projection Data

The modified model SWAT+ was also used to simulate the climate change impact on future snowmelt runoff in the basin by coupling it with climate projections. Due to the limited spatial resolution, simplified physical processes, and incomplete understanding of the climate system, there are systematic deviations in the climate model projection data. As such, many bias correction methods have been developed and applied, e.g., delta method, multiple linear regressions, and quantile mapping [73,74,75]. In this study, we chose the delta method, which is widely used for climate projection data bias correction for its easy understanding and operation [76,77].

The delta method can eliminate the bias between climate scenario data and observations effectively [78]. We applied bias correction to two types of projection data inputs: precipitation and temperature. For precipitation, the correction is made using the variation ratio between the modeled precipitation and observed precipitation [79]. Similarly, temperature bias correction is based on the absolute variation between climate model temperature and observed temperature. The two correction methods are shown in Equations (12) and (13):

$${\mathrm{P}}_{\mathrm{adj}}={\mathrm{P}}_{\mathrm{ori}}\times \frac{{\mathrm{P}}_{\mathrm{obs}}}{{\mathrm{P}}_{\mathrm{ori}}}$$

(12)

$${\mathrm{T}}_{\mathrm{adj}}={\mathrm{T}}_{\mathrm{ori}}+\left({\mathrm{T}}_{\mathrm{obs}}-{\mathrm{T}}_{\mathrm{ori}}\right)$$

(13)

where ${\mathrm{P}}_{\mathrm{adj}}$ and ${\mathrm{P}}_{\mathrm{ori}}$ are the corrected and original precipitation of climate model data and ${\mathrm{P}}_{\mathrm{obs}}$ is observed precipitation;${ \mathrm{T}}_{\mathrm{adj}}$ and ${\mathrm{T}}_{\mathrm{ori}}$ are the corrected and original temperature of climate model data. ${\mathrm{T}}_{\mathrm{obs}}$ is observed temperature.

4. Results

4.1. Model Performance Assessment

This section compares the SWAT+ model and the original SWAT model in terms of their performance for daily snowmelt runoff simulation in the study area. The results of the daily snowmelt runoff simulation (March–April) for the calibration period (1980–1999) and the validation period (2000–2019) are shown in Figure 7. It can be seen that, for low flow regimes, the snowmelt runoff is largely underestimated by the original SWAT model setting for both the calibration and validation periods, whereas the SWAT+ model performs better in this case.

The performances from the two models are further evaluated using the selected metrics and the results are shown in

Table 2. Snowmelt runoff simulation accuracy of SWAT model and SWAT+ model.

Metrics	Calibration			Validation
	SWAT	SWAT+	Change	SWAT	SWAT+	Change
MAE	98.66	77.27	−21.68%	106.98	84.99	−20.56%
RE	−32.09	−20.85	−35.03%	−343.90	−23.67	−93.12%
RMSE	155.76	117.41	−24.62%	144.39	121.05	−16.16%
NSE	0.55	0.74	0.19	0.47	0.63	0.16
R2	0.63	0.78	0.15	0.63	0.70	0.07

, which indicates a better performance of the SWAT+ model compared with the original SWAT model, showing improvements of MAE, RE, and RMSE by 20.56%, 93.12%, and 16.16%, respectively. In short, the snowmelt runoff simulation performance of the SWAT+ model is improved in comparison with the SWAT model, especially when the snowmelt runoff is small. This improvement in snowmelt runoff simulation underlines the benefits of the modification we made to the original SWAT model, i.e., adding total solar radiation in the SWAT model and localization of snowmelt temperature thresholds and snowmelt factors.

4.2. Projected Change in Snowmelt Runoff under Climate Change

Table 3. Precipitation (P) and temperature (T) during the snowmelt runoff period from 2025 to 2054.

Baseline		SSP2.6		SSP4.6		SSP8.5
		Projections	Change	Projections	Change	Projections	Change
T	1.05 °C	2.53 °C	1.48 °C	3.22 °C	2.17 °C	4.07 °C	3.02 °C
P	66.14 mm	75.82 mm	14.64%	72.29 mm	9.29%	78.89 mm	19.28%

shows the projected precipitation and temperature during the snowmelt runoff period (March–April) from 2025 to 2054 and their corresponding changes relative to the baseline period (1980–2019). It appears that the Baishan basin will be warmer during 2025–2054 compared with the baseline period. Annual average temperature in 2025–2054 will increase by 2.53, 3.22, and 4.07 °C under the SSP2.6, SSP4.5, and SSP8.5 scenarios, respectively, compared with the baseline annual average temperature. In addition, the projected annual average precipitation will increase by 14.64% under SSP2.6, 9.29% under SSP4.5, and approximately 19.28% under SSP8.5.

The SWAT+ model was used to simulate the snowmelt runoff for both the baseline period (1980–2019) and the future (2025–2054) to assess the impact of climate change. For the baseline period, the simulation is driven by historical observations, whereas the future case is conducted using the projected forcing data, i.e., precipitation and temperature. As shown in Figure 8, the overall snowmelt runoff will reduce from March to April under the future climate scenarios compared with the baseline period. The annual average snowmelt runoff under SSP2.6, SSP4.5, and SSP8.5 is 148.3, 130.76, and 117.76 m³/s, with decreases of 31%, 39%, and 45%, respectively, from the baseline period (214.14 m³/s). The ensemble average of the snowmelt runoff presents clear and significant increase in March and a decrease in April, with peak flows shifting from April to March. This is due to the fact that snowmelt peak shifts to earlier spring as the temperature increases, causing significant increase in snowmelt runoff in March and a considerable decrease in April.

5. Discussion

5.1. Effect of the Modified Temperature Index Method

The traditional temperature index method in the SWAT model cannot represent the real snowmelt physical processes and leads to low simulation accuracy of snowmelt runoff. Therefore, some studies proposed to couple an energy balance model to the SWAT model to improve the low accuracy of snowmelt runoff simulation and reported good results [80,81]. However, other studies compared the performance of the improved temperature index model and the energy balance model and concluded that the improved temperature index model performed better than the energy balance model, with an extra benefit of requiring much less data [33,82]. In this study, we added the total solar radiation, which includes direct solar radiation and scattering solar radiation, to the traditional temperature index method, and further modifies the seasonal variation formula and determines the snowmelt temperature threshold. Compared with the SWAT model simulation results, the SWAT + decreased by 21.68% to 35.03% in MAE, RE, and RMSE metrics during calibration period. The NSE and the R² indicators increased by 0.19 and 0.15, respectively. During the validation period, the MAE, RE, and RMSE indicators increased by 16.16% to 93.12%, while the NSE indicators and R² indicators increased by 0.16 and 0.07.

5.2. Snowmelt Temperature Threshold

The 0 °C threshold used by the original SWAT model is only a rough indicator of snowmelt [83]. Some scholars noticed this problem and chose to use the accumulated temperature to determine whether the snowmelt has melted [84,85,86]. However, when calculating the accumulated temperature, snow melting below 0 °C is also ignored. This study determines the snowmelt temperature threshold by using the snow depth and temperature data from remote sensing. It helps mitigate the problem of poor availability of spatial and temporal coverage in snow observation data due to sparse distribution of snow measurement stations, especially in mountainous areas. The result shows that a better estimate of the temperature threshold of snowmelt in the Baishan watershed is around −0.54 °C, which is substantially lower than the prescribed 0 °C. This also confirms the results from some previous studies [87,88] and supports the argument that snow may melt when the temperature is well below 0 °C [89]. The reason is that there are many factors affecting snowmelt temperature threshold, such as impurities in the snow [90]. As such, snowmelt does not merely depend on the air temperature, it also depends on the balance of the energy fluxes, e.g., it can occur through solar heating when the air temperature is below 0 °C.

5.3. SWAT+ Model Uncertainty

Even though our results show the benefits of the SWAT+ model for snowmelt simulation in a seasonal snowmelt basin, we also recognize that the introduction of new variables will bring in extra uncertainty to the model. Studies on SWAT model uncertainty have been abundant in recent years; for example, the SUFI-2 algorithm and GLUE algorithm are widely used in the uncertainty analysis of the SWAT model [91,92]. This paper uses the SUFI-2 algorithm to judge the uncertainty of the SWAT+ model based on the criteria of the 95% prediction uncertainty (95 PPU) [93]. The 95 PPU has two characteristic indicators: P-factor and R-factor. P-factor is the percentage of observed data that are enveloped by 95 PPU and R-factor is the average thickness of 95 PPU. Larger P-factor and smaller R-factor indicate less model uncertainty. The simulation uncertainty is generally considered acceptable for P-factor > 0.7 and R-factor < 1 [94]. The P-factor of the SWAT + during the validation period was 0.86 and 0.67, respectively, and the R-factor was 0.92 and 0.77, respectively. It shows that the SWAT+ model uncertainty is acceptable.

5.4. Limitations and Further Studies

Our results show that the SWAT+ model uncertainty is acceptable. Further, there is also a certain amount of uncertainty in the process of determining the snowmelt threshold, snowmelt factor, as well as in the data derived from remote sensing, which needs to be fully examined and studied. The snow depth data used in this paper are derived from a passive microwave remote sensing data product. Although this product has been confirmed to have good applicability in Northeast China in general [61], its accuracy over the study area has not been extensively evaluated. At present, the accuracy evaluation of remote sensing data can be divided into two categories. One is the mutual verification conducted using observed data [95,96] and the other one can be mutually verified according to other types of remote sensing data [97,98]. Due to the lack of measured data, the future study will evaluate accuracy according to other types of remote sensing data. There are two ways to obtain multiple remote sensing snow depth data. One is to extract information about snow depth based on satellite remote sensing data by ourselves, for example, using a new algorithm to extract snow depth based on passive microwave measurements [99] or using relatively high spatial resolution active microwave measurements to capture the spatial variability of snow depth in a mountain area [100]. The other is to use ready-made remoting sensing snow depth data for accuracy evaluation, such as America National Snow and Ice Data Center AMSR-E data (http://nsidc.org/data/docs/daac/ae_swe_ease-grids.gd.html, accessed on 5 September 2022) and Japan Aerospace Exploration Agency AMSR2 data (https://suzaku.eorc.jaxa.jp/, accessed on 8 September 2022). Due to the limit of the scope of this study, the uncertainty quantification and representation of remote sensing data used in this paper are not included in this study. However, further work is planned to investigate the uncertainty associated with remote sensing data uncertainty [101,102].

6. Conclusions

SWAT model is a distributed hydrological model with a physical foundation and is widely used in hydrological simulations. The snowmelt module in the SWAT model adopts a relatively simple temperature index method together with assuming both the dates for maximum and minimum snowmelt factors and the snowmelt temperature threshold, which leads to underestimating snowmelt runoff obviously. This paper presents a study of improving snowmelt runoff simulation using SWAT by systematically addressing the issues mentioned above. Firstly, this study integrates the total solar radiation into the temperature index. Secondly, accurate dates of the maxima and minima of snowmelt factors are determined and snowmelt factor seasonal variation formulas are modified. Thirdly, the temperature threshold of snowmelt is determined by identifying the relationship between snow depth and temperature. Finally, the Baishan basin in Northeast China was selected to test SWAT+ model performance. The results show that SWAT+ performs better than the SWAT model in seasonal snowmelt zones, especially when the snowmelt runoff is in its low to medium range. In comparison with the SWAT model, the SWAT+ model has remarkable improvements in all selected performance metrics, i.e., MAE, RE, RMSE, R², and NSE. Further, the effects of climate change on snowmelt runoff in the Baishan basin were investigated using the SWAT+ model. Compared with the baseline period (1980–2019), overall snowmelt runoff will reduce under SSP2.6, SSP4.5, and SSP8.5 scenarios. The annual average snowmelt runoff under SSP2.6, SSP4.5, and SSP8.5 are 148.3, 130.76, and 117.76 m³/s, a decrease of 31%, 39%, and 45%, respectively, from the baseline period (214.14 m³/s).