A Review on Snowmelt Models: Progress and Prospect

Abstract

The frequency and intensity of flood events have been increasing recently under the warming climate, with snowmelt floods being a significant part. As an effective manner of simulating snowmelt flood, snowmelt models have attracted more and more attention. Through comprehensive analysis of the literature, this paper reviewed the characteristics and current status of different types of snowmelt models, as well as the different coupling methods of models for runoff generation and confluence. We then discussed key issues in snowmelt modelling, including blowing snow model, frozen ground model, and rain-on-snow model. Finally, we give some perspectives from four aspects: data, model structure, forecast and early warning, and forecast and estimation. At present, most of the snowmelt models do not have blowing snow or frozen ground modules. Explicit consideration of blowing snow and soil freezing/thawing processes can improve the accuracy of snowmelt runoff simulations. With climate warming, rain-on-snow events have increased, but the mechanism of enhanced rain and snow mixed flooding is still unclear, particularly for the mechanism of rain-snow-ice mixed runoff generation. The observation and simulation of rain and snow processes urgently need further study. A distributed physical snowmelt model based on energy balance is an advanced tool for snowmelt simulation, but the model structure and parameter schemes still need further improvements. Moreover, the integration of satellite-based snow products, isotopes, and terrestrial water storage change, monitored by gravity satellites, can help improve the calibration and validation of snowmelt models.

1. Introduction

Floods are some of the most frequent and severe natural disasters in the world [1]. Since the 20th century, river floods have caused approximately 7 million deaths and direct annual losses of 10.4 billion U.S. dollars [2]. About two-thirds of the area in China faces the threat of floods [3]. Snowmelt floods often take place in northwestern China, which causes damages to road traffic, downstream reservoirs, channels, and other engineering facilities, and threatens the safety of people’s lives and properties. In the context of global warming, the speed of ice and snow melting has accelerated, causing advanced ablation periods and earlier snowmelt flooding peaks [4]. In addition, due to increasing rain-on-snow (ROS) events, the frequency and intensity of rain and snow mixed floods have increased significantly, and the losses caused by flood disasters are increasing accordingly [5].

From 1950 to 1993, snowmelt floods destroyed 142 small- and medium-sized reservoirs and dams in Xinjiang Province, China. In recent decades, serious snowmelt floods have impacted Wusu, Manas, Hutubi and many other places in Xinjiang. In a snowmelt flood event in March 2005, 110,000 people in the Yili River Valley were affected, with 12,000 houses being destroyed, 3,736 livestock drowned, and 796 hm² of farmland submerged [6]. In recent years, rain and snow mixed floods have occurred frequently in China. For instance, a rain and snow mixed flood was observed in the upper reaches of the Heihe River in the Qilian Mountains on 17 June 2013. As rainfall strongly eroded the snow surface in high mountainous areas, rapid rain and snow mixed runoff was formed. The total economic loss due to rain and snow mixed flooding exceeded 5.5 million yuan [5]. Therefore, accurate simulation of the snowmelt flood process is of great importance for flood disaster prevention and mitigation.

Hydrological models and hydraulic models are important tools for hydrology studies. The former focus on simulating the process of flow generation and convergence from the basin to the main river channel, while the latter focus on simulating the evolution process of water flow in the main river channel. Figure 1 shows the formation process of a snowmelt flood. Massive snow accumulation in mountainous areas may easily cause flooding once there is a rapid heating or rain event. In the meantime, mountainous terrain is conducive to the rapid concentration of snowmelt water. Therefore, when modelling snowmelt flood processes, it is better to couple two different types of models. The results of Hydrological modelling for different river sections can be used as input conditions for the hydraulic model, together with previous observations, so as to realize a rolling iterative flooding forecast [7].

The snowmelt model is a hydrological model that simulates the processes of snow accumulation, melt, and confluence. Since the snowmelt process involves the complex energy transfer and water transfer processes between the atmosphere and snow, between different snow layers, and between snow layers and soil, it is more complicated than the formation process of rainfall runoff [8]. If there is a rainfall event on snow, it is very difficult to simulate the runoff generation process due to the mixture of rain and snow. Snowmelt simulation is the basis for the simulation of snowmelt flood formation processes, and it plays a crucial role in accurately modelling snowmelt floods.

One remaining issue of snowmelt flood forecasting is the low accuracy of early warnings [5] for rain-snow-ice mixed floods, particularly under the background of rapid warming. A comprehensive review of the related progresses and challenges of snowmelt process models is critical to further development of the related field. Therefore, this study reviews the current status of different snowmelt models by extensively consulting relevant literature, summarizes the trends and challenges, and offers perspectives on the future of snowmelt flood simulation.

2. Snowmelt Modelling

2.1. Development History of Snowmelt Models

The development of snowmelt models has a long history. At the earliest stage, river flow is estimated by establishing a statistical relationship between measured variables (e.g., snow cover area, snow water equivalent, etc.) and snowmelt runoff [9]. Furthermore, empirical equations are also frequently used in this period. The principle of empirical equations, otherwise called the degree-day factor [10], is to assume that there is a linear relationship between air temperature and snowmelt. The most classic degree-day snowmelt model, Snow Runoff Model (SRM), was proposed by Martinec in 1975 [10]. Since then, it has been revised and improved constantly for easy generalization and wide applicability [11,12]. So far, the SRM model, recommended by the International Meteorological Organization [12], has been widely used in over 100 watersheds around the world. Because the degree-day model is easy to use and usually has high accuracy, it has been widely used. However, its disadvantage is that it cannot simulate the physical process of snow melting.

In the 1950s, the US Army Corps of Engineers first calculated snowmelt based on the energy exchange process between snow and the environment [13]. The basic principles of the energy balance model are energy balance and water balance, which is of strong physical meanings. After a series of applications and improvements [14,15,16], the energy balance equation has becoming more and more optimized. However, the energy balance model in this period is still a simple point-scale model.

In the 1970s, the rapid development of computer technology made the in-depth development of distributed physical hydrological models possible. In 1986, the first representative distributed physical hydrological model, the European Hydrological System (SHE), was developed [17]. In this model, the watershed plane is divided into grids in order to deal with the spatial distribution of model parameters, precipitation input and hydrological responses. The vertical plane is divided into several horizontal layers in order to deal with the problem of soil-water movement in different layers, and the energy balance is then used to calculate the snowmelt in each grid individually.

Since the 1990s, a series of distributed snowmelt models have been developed, such as SNOBAL [18], the Utah Energy Balance model (UEB) [19], Hydrological Simulation Program—FORTRAN(HSPF) [20], etc. Some distributed physical hydrological models also introduced snowmelt modules (including degree-day snowmelt algorithm and energy balance algorithm), for example the Distributed Hydrology Soil Vegetation Model (DHSVM) [21], Variable Infiltration Capacity (VIC) [22], and the Soil and Water Assessment Tool (SWAT) [23]. Moreover, based on energy balance models, two-layer or multi-layer snowmelt models have been developed, such as Snow Thermal Model (SNTHERM89) [16] and SNOWPACK [24]. These models, now widely used in snow engineering, take into consideration more detailed vertical distribution of snow.

In recent years, machine learning has become a revolutionary and versatile tool that has changed industrial applications and provided new and improved capabilities for scientific discovery and model development. It is increasingly used for hydrological modeling and hydrological prediction, becoming a new tool for hydrological research. The rapid development of machine learning is bringing more opportunities for snowmelt modeling, or so-called data-driven modeling.

2.2. Categories of Snowmelt Models

According to the different ablation algorithms used, snowmelt models can usually be classified into four categories: statistical, conceptual, physical and data-driven (green categories in Figure 2). Statistical snowmelt models use statistical methods or black box modules to establish the relationship between a certain snow hydrological characteristic parameter (e.g., snow area) and runoff to predict runoff. Conceptual snowmelt models usually establish an empirical relationship between snowmelt and temperature. These models, such as SRM and HBV [25], have physical meanings and are mature and widely used methods. The physical snowmelt models, such as SNTHERM and SNOWPACK, calculate snowmelt based on the energy balance of snow cover. They have strict physical meaning and have been widely used in distributed modelling of snowmelt. The emergence of data-driven models has benefited from the growth of massive data and the rapid increase in computational power. These models, Such as ANN [26] and LSTM [27], simulate the changes in snowmelt runoff using machine learning algorithms to select appropriate parameters (e.g., daily rainfall, temperature, solar radiation, snow area, snow water equivalent) from different data sources.

According to the spatial distribution characteristics of models, the snowmelt model can be divided into lumped, semi-distributed and distributed models (blue categories in Figure 2). Lumped models, representing by SRM, assume that snow cover and melting rate are consistent across the watershed or sub-basin. Semi-distributed models, including HBV, SWAT and the Precipitation-Runoff Modeling System (PRMS) [28], usually divide the watershed into sub-watersheds or hydrological response units according to a certain characteristic of the watershed (e.g., elevation, vegetation coverage, land use types, and topographic factors). Distributed models, such as SHE, DHSVM, and VIC, divide the watershed into grids, and assign different parameter values to each grid and then include wind, snow, and soil freezing processes in each individual grid. In recent years, with the rapid improvement of computing power, multi-layer hybrid and nested models are appearing, such as the Water and Energy Budget-based Distributed Hydrological Model (WEB-DHM) [29]. Sub-catchments, grids, and slopes are used for nested calculations, thereby improving the ability to adequately describe the snowmelt process.

From the perspective of modelling complexity, snowmelt models vary greatly. Therefore, when selecting a snowmelt runoff simulation model, it is necessary to comprehensively consider data availability, watershed characteristics, and time scale of simulation.

2.3. Statistical Snowmelt Models

The most widely used statistical method is to establish a relationship between the maximum snow water equivalent during the spring snowmelt period and the total amount of runoff, as in Equation (1),

Q = a + b SWE

(1)

where Q is the total amount of runoff, and SWE is the snow water equivalent. SWE can be obtained by manually measuring snow depth along a certain survey line or by automatic snow pillow measurement. The coefficients a and b are obtained by regression analysis of observed data. For example, DeWalle et al. [30] used early winter SNOTEL snowpack-water-equivalent data to predict snowmelt-runoff volumes of the Upper Rio Grande basin in Colorado. Identifying accurate snow sampling points that can represent the hydrological characteristics of snow cover in the watershed is key for statistical snowmelt models. Principal component regression analysis is usually used to find the optimal combination of survey lines or snow pillow measurement locations.

As the peak flow often has a good correlation with the total spring flow, Equation (1) can also be extended to predict spring peak flow and other hydrological parameters. There are other independent variables that may affect this relationship, such as accumulated rainfall in the previous days of snow cover and river base flow before snow melting. This model can also be used to provide indicators for soil water content and groundwater conservation before snowmelt.

Another statistical method for predicting snowmelt runoff is to use runoff characteristic curves. This method uses regression analysis to fit the ratio of snow cover in the spring melting period to the total runoff in the basin to obtain a series of characteristic curves. These characteristic curves should be consistent in different years in each basin, but different basins have different characteristics curves. Therefore, as long as there is an estimate of the snow cover area, these curves can be used to estimate the remaining snowmelt runoff. In addition, the characteristic curve can also provide information about the time when the characteristic flow occurs. As it becomes ever more convenient to obtain snow cover data regularly using remote sensing, this method is also able to provide a forecast of the runoff during the melting period.

Under average or typical situations, as long as there are enough observations, statistical models can be of high accuracy, but it is not appropriate for special situations that have not occurred before [9]. Besides, the generation mechanism of snowmelt runoff is still not clear yet. Therefore, for snowmelt flood process simulation, statistical snowmelt models established for one watershed are difficult to directly apply to other watersheds.

2.4. Conceptual Snowmelt Models

The conceptual snowmelt model is usually established based on the degree-day factor, so it can also be referred to as the degree-day factor model. The concept of degree-day factor was first proposed by Finsterwalder and Schunk [31] in 1887 for studying glacier changes, and was subsequently improved by many researchers. It is now widely used for ice and snow melting studies in the Alps [32], northern Europe [33], the Qinghai-Tibet Plateau [34], and the Greenland ice sheet [35]. Its basic assumption is that there is a linear relationship between the change of snowmelt rate and the daily average temperature above a certain temperature. The general form of the degree-day factor ablation algorithm is as Equation (2),

M = DDF (Ta−Tb)

(2)

where M(cm·day^‒1) is the amount of snow melting, Ta (°C) is the daily average temperature, and Tb (°C) is the threshold temperature. The threshold temperature Tb (°C) is generally set to be 0 °C, but in special cases, other values may be used. DDF (cm·°C^‒1·day^‒1) is the degree-day factor of snowmelt. It was found that DDF of snow varied greatly across different regions, ranging from 0.1 to 1 cm·°C^‒1·day^‒1. DDF can be determined based on measured data or empirical formula (DDF = 1.1·Ps/Pw. Ps is the density of snow and Pw is the density of liquid water). DeWalle et al. [36] suggested that SNOTEL data of the western United States can be used to calculate an accurate DDF value for a specific watershed. Zhang et al. [37] analyzed the spatial distribution of the degree-day factors in western China based on the observations of ice and snow melting and temperature in different glacier regions in High Asia, and then applied the results to calculate large-scale glacier meltwater.

There are many improvements to the conceptual snowmelt model, which can be summarized into three categories (

Table 1. Improvement of conceptual snowmelt model.

Improvement Category	Targeting Problem	Improvement Strategy	References
Improvement of degree-day factor	The value of the degree-day factor can be affected by the state of snow, forest cover, terrain, snow pollution, etc.	Divide the sub-regions according to a certain characteristic of the watershed, and enter the degree-day factors with temporal and spatial differences.	[11,38,39,40,41,42]
Increase degree items	The degree-day factor based on the temperature index cannot represent all the energy sources in the process of snow melting.	Further consider the energy source of snow melting through addition of parameters including solar radiation, wind speed, water vapor pressure, etc.	[10,43,44,45,46,47]
Change the model structure	Temporal and spatial resolution of the degree-day factor model is insufficient.	Establish a distributed-degree-time model based on digital elevation data (DEM) and spatial data of temperature and radiation.	[48,49,50,51,52]

). The first category is to improve the core parameters of the snowmelt algorithm, taking into account various factors that affect the degree-day factor, including the snow state, density, forest coverage, snow surface pollutants, topography, etc. For example, Martinec [11] used the empirical formula of snow density and DDF to calculate the change of DDF in different seasons. Anderson [38] simulated the seasonal variation between the minimum and maximum values of DDF in a year using a sine function. Federer and Lash [39] used a linear increment model to simulate DDF changes from January 1 to June 23.

Considering that forest cover may intercept a part of the solar radiation and slow down the wind speed, DDF in a forest cover area is lower than that in the open area. Marince et al. [11] suggested that the DDF of a vegetation coverage area should be adjusted appropriately according to the vegetation coverage status. Federer et al. [40] gave the range of DDF based on average daily temperature in North America as follows: 0.45–0.75 (No vegetation), 0.27–0.45 (fallen leaves), 0.14–0.27 (needles). Kuusisto [41] established a canopy coverage empirical calculation formula with DDF. Furthermore, the presence of dust and debris on the snow surface will reduce the snow surface albedo and increase the DDF. Singh et al. [42] considered the difference in DDF between clean snow and contaminated snow when estimating snowmelt in the Himalayan area. Several studies found that the difference in input shortwave radiation was caused by terrain, and thus tried to improve the accuracy of spatial distribution of degree-day factors by combining terrain factors such as altitude, slope, and aspect [53,54]. Compared with treating DDF as a constant, these methods greatly improved the simulation accuracy of the degree-day-factor-based conceptual snowmelt models.

The second category of improvement is to increase the energy input items when the snow melts. Zuzel et al. [43] found that adding wind, vapor pressure, and net radiation variables in the snowmelt prediction equation can significantly improve the simulation accuracy. Brubaker et al. [44] introduced radiation variables in the degree-day factor ablation algorithm, and proposed Equation (3),

M = m_QR_d + a_rT_d

(3)

where M (cm·day^‒1) is the amount of snow ablation, m_Q (cm·W⁻¹·m²·day⁻¹) is the physical constant of the depth at which energy is converted into water (generally taken as 0.026 cm·W⁻¹·m²·day⁻¹), R_d (W·m⁻²) is net radiation index, and a_r (cm·°C^‒1·day^‒1) is the corrected degree-day factor, not the same as DDF in Equation (2). Martine et al. [10] found the value to be 0.20–0.25 cm·°C^‒1·day^‒1 in the Swiss Alps. Td(°C) is the degree-day index that is equal to (Ta-Tb) in Equation (2). Morid et al. [46] found the value to be 0.20 cm·°C^‒1·day^‒1 at Ammameh station through observational studies. Adding radiant energy into the degree-day factor supplements the snowmelt calculation with energy exchange, which has a certain physical meaning and can improve the simulation accuracy. Kustas et al. [47] improved the simulation accuracy of the SRM model by nearly 40% by introducing the radiation factor.

The third category of improvement is changes to the model structure, such as the emergence of the spatially distributed degree-day model. In order to improve the temporal and spatial resolution of the snowmelt estimation of degree-day factor, many scholars have established distributed degree-day models based on digital elevation data (DEM) and considered the spatial difference of air temperature and radiation. Cazori [48] divided the basin into 20 m × 20 m grids and established a distributed hourly snowmelt model. The model uses air temperature and terrain-corrected clear sky shortwave radiation to improve snowmelt estimation. Hock et al. [49] establish a distributed hourly snowmelt model at 30 m resolution, using degree-day factors and calculations of clear sky shortwave radiation that took into account local slope and aspect. Pellicciotti et al. [50] proposed a distributed degree-day modelling approach that combines shortwave radiation and albedo. Moreover, they compared the difference in modelling performance between lumped and distributed models, and highlighted the robustness of the distributed model when appropriate extrapolation schemes for the redistribution of the input variables (e.g., temperature, albedo, and radiation) are used. Fang et al. [51] designed a distributed snowmelt model for the Juntang Lake basin in Xinjiang province, China. The concept of the unit time period was introduced in their snowmelt simulation, which resulted in a distributed degree-hour snowmelt model with high accuracy. Jost et al. [52] established an improved distributed degree-day model in the forested watershed, and proved its superiority to the classic degree-day model and Distributed Hydrology Soil Vegetation Model (DHSVM). There are also some widely used distributed hydrological models that also introduce the degree-day factor algorithm to simulate snow melting with spatial and temporal characteristics and have high accuracy, such as the Soil and Water Assessment Tool (SWAT) [23], Distributed Water-Heat Coupling Model (DWHC) [55], Geomorphology-Based Hydrological Model (GBHM) [56], etc.

Traditional conceptual snowmelt models based on degree-day factors assume that ablation occurs uniformly across the entire watershed. Most models are lumped or semi-distributed. Therefore, snow melting simulation has large uncertainties in their spatial distribution. In addition, the degree-day factor is constant on a daily scale. It is thus difficult to capture the intra-day change in snowmelt rate, especially the difference of the temperature and snow energy supply between day and night. The improved spatial distribution degree model fully considers the time and space differences and introduces energy exchange terms. It is therefore able to simulate the changes in snowmelt runoff on a daily scale and achieve fair results.

2.5. Physical Snowmelt Models

In 1956, the US Army Corps of Engineers first calculated the amount of snowmelt based on the energy exchange between snow cover and the environment. Anderson et al. [57] adopted the energy balance equations into the process of snow accumulation and ablation, and coupled them with the Stanford Model for snowmelt runoff simulation. Their results suggest that the physical process of snow ablation has a good simulation effect. Subsequently, they proposed a complete point-scale snow energy and mass balance equation. The general form of the energy balance snowmelt model is as Equation (4),

λM + Q_I = Q_NR + Q_S + Q_L + Q_G + Q_P

(4)

where λ is the heat of dissolution of snow (generally taken as 3.34 × 10⁵ J/kg), M is the amount of snowmelt, Q_I (J/m²) is the change in heat storage of snowdrifts, Q_NR (J/m²) is the net radiant flux, Q_S (J/m²) is the sensible heat flux, Q_L (J/m²) is the latent heat flux, Q_G (J/m²) is the heat flux from the soil, and Q_p (J/m²) is the heat transferred into the snow pile by rainfall. Figure 3 shows the components in the energy balance of snow surface. The energy exchange of the snow surface can be estimated from the relevant data of the atmosphere, and the energy consumed by snow melting can be converted into the amount of snow melting water.

The internal process of the snow can be better estimated by dividing it into many layers when the snow is thick enough. Anderson [14] uses a finite difference to calculate a multi-layer snow model, which considers fresh snow, compaction, change and retention, and liquid water transmission across different layers. Jordan [16] developed a one-dimensional numerical layered snow thermal model (SNTHERM.89), which divides the snow into multiple horizontal layers with each layer being controlled by the energy and mass balance equations. It can better simulate the features on the snow profile, such as snow depth, snow density, snow temperature, etc. Marks [8] took into account the temperature difference between the snow surface and snow layer, and proposed a two-layer energy balance snow melting model (SNOBAL). The critical freezing depth was calculated by the temperature gradient between the melting layer and the non-melting layer. The snow model in the Precipitation-Runoff Modeling System (PRMS) [28] is also a two-layer energy balance model. Bartelt et al. [24] developed a one-dimensional physical snow pack model (SNOWPACK). In this model, snow layers are defined by their size as well as their macroscopic properties (i.e., bulk density, mean stress, water content or temperature) and microscopic properties (i.e., ice grain size and shape, bond size or coordination number). This model can predict snow settlement, stratification, surface energy exchange and mass balance, and can be used to predict the occurrence of avalanches.

As the physical process of energy exchange between snow and air interface can be measured well, the snow ablation simulation based on energy balance can achieve high simulation accuracy at a point or in a small field [14]. However, there are many difficulties for large watersheds. For example, how to extrapolate the measurement data of a point scale or a small area to a large area, and how to determine the representativeness of the observation data on the point. In particular, in the process of turbulent exchange, how sensible heat and latent heat distribute over a large scale is not clear. In addition, the application of a complete energy balance to calculate the snow melt requires extensive inputs, such as wind speed, humidity, air temperature, etc., but the snow area usually has few monitoring sites, which makes it difficult to parameterize for large areas. Therefore, various empirical extrapolation methods have been applied in parameterization. Several studies [58,59,60,61] have examined the snow energy balance process in mountainous, grassland, and forest environments, and proposed and evaluated different empirical methods for energy flux estimation. Steven [62] compared the amount of snowmelt estimated by the energy balance between forest areas and bare land during rainfall, and found that the sensible heat and latent heat flux of bare land is almost three times that of forest areas because of its higher wind speed. Lu et al. [63] calculated the energy budget of the open field and seasonal snow under the forest canopy during the snowmelt period in the western Tianshan Mountains of China. Marks et al. [64] simulated and measured the sensible heat flux and latent heat flux between the snow under the pine forest and the atmosphere. Pomeroy et al. [65] conducted a simulation experiment on the turbulent flux in the process of blowing snow. However, research on the sensitivity of radiation flux and turbulence exchange for altitude, latitude, and vegetation density still needs to be further investigated [66].

The physical distributed snowmelt model has become increasingly popular. For example, the European Hydrological System (SHE) [17] uses a distributed single-layer energy balance snowmelt model to simulate the snow melting process. Wigmosta et al. [21] developed a distributed hydrological soil vegetation model (DHSVM), which includes a single-layer canopy snow interception and a double-layer energy balance snowmelt model under the canopy. The model clearly combines the effects of terrain and vegetation cover on the energy exchange of snow surfaces. Tarboton et al. [19] used a single-layer energy balance snowmelt model in the Utah Energy Balance model (UEB). Liang et al. [22] also used a single-layer energy balance snowmelt model in the Variable Infiltration Capacity (VIC) large-scale hydrological model that considers the freezing and thawing process of the soil.

The physical distributed snowmelt model usually has high temporal and spatial resolutions, and it has strict physical meanings. Compared with the conceptual snowmelt model, the calculation of snowmelt runoff in the study of complex terrain and vegetation coverage has more advantages. However, there are still many difficulties in the accurate spatial parameterization of various model parameters. Some recent studies [67] compared the VIC-CAS model simulation with the snow distribution and snow depth obtained by remote sensing, and found that the model has low accuracy for low snow years and thin snow simulation. Some studies [29] have successfully used remote sensing-derived time series snow products to directly drive the hydrological model.

With the rapid development of big data and high-performance compute capacity, more and more physical parameters can be accurately retrieved by remote sensing. The energy and hydrological parameters of the physical snowmelt model and the research on the parameter spatial discretization have made continuous progress. Distributed snowmelt models have developed to become key approaches in the study of snowmelt runoff forecasts.

Table 2. Comparison of main snowmelt models.

Model	Category	Spatiotemporal Resolution	Snow Melting Algorithm	Frozen Ground Consideration	Blowing Snow Consideration	Pros and Cons	References
SRM	Lumped, Conceptual	Daily, Elevation zone	Degree-Day	N	N	Less input data, less model parameters, flexible and easy to use, but it has high simulation accuracy, which is suitable for the snowmelt runoff forecast in mountainous areas with a lack of data, the forecast period is not precise enough	[11]
HBV	Semi-distributed, Conceptual	Daily, Elevation zone	Degree-Day	N	N	Less input data, less model parameters, flexible and easy to use. Used for hydrological forecasting of glaciers and snow-covered fields; lack of physical process research	[25]
SWAT	Semi-distributed, Physical	Daily, HRU	Degree-Day	N	N	Is widely used and becoming mature, and the simulation accuracy is high. Facing the long-term scale of large and medium watersheds; however, it has high requirements for input data and parameters, and it is difficult to apply in areas where monitoring data is lacking	[23]
GBHM	Distributed, Physical	Daily, Hillslope-based	Degree-Day	N	N	The application of the slope flow hydrological unit can better simulate the natural runoff process of the watershed. The model parameters have clear physical meanings and greatly reduce the amount of model calculations. It is suitable for the simulation calculation of the complete hydrological process of large and medium-scale watersheds in mountainous areas	[56]
DWHC	Distributed, Physical	Daily, Grid	Degree-Day	Y	N	The model has few adjustable parameters and input variables. It only needs conventional meteorological data and basic soil and vegetation parameters to continuously calculate the hydrological cycle elements such as the temperature and liquid water content of each layer of soil. Suitable for inland river areas in alpine mountainous areas; the model adopts a simplified method for many processes, which reduce the rigor of the physical process of the model	[55]
SNTHERM	point-scale, Physical	Hourly, point-scale	Energy Budget	Y	N	Is based on a more sophisticated physical basis and can well simulate various characteristic parameters on the snow profile, including snow depth, density, particle size and temperature. These parameters happen to be the input parameters required by the microwave remote sensing radiation transmission model. Can better simulate the temperature layer of the snow layer and frozen soil layer	[16]
SNOWPACK	point-scale, Physical	Hourly, point-scale	Energy Budget	N	Y	Only four micro-structure parameters can be used to describe the complex snow structure and simulate the change characteristics of the internal structure of the snow, but the snow water equivalent data cannot be obtained, and it is difficult to simulate the snowmelt runoff. Suitable for mountain avalanche prediction	[24]
PRMS	Semi-distributed, Physical	Hourly, HRU	Energy Budget	N	N	With a variety of simulation functions, it can simulate changes, in water balance, flood peaks and peak discharges, soil water, etc., in the processes of general precipitation, extreme precipitation and snow melting. The physical processes involved are also more comprehensive	[28]
DHSVM	Distributed, Physical	Hourly, Grid	Energy Budget	N	N	The model fully reflects the interaction and feedback mechanisms of climate, vegetation, snow cover, soil, and hydrology. However, it has high requirements for input data and parameters, and it is difficult to apply in areas where monitoring data are lacking.	[21]
SHE	Distributed, Physical	Hourly, Grid	Energy Budget	N	N	The model is very powerful, and software development is very mature. It can simulate all major hydrological processes in the terrestrial water cycle, including water movement, water quality, and sediment transportation; however, the model requires high input data, too many parameters, and a large amount of calculation.	[17]
VIC	Distributed, Physical	Hourly, Grid	Energy Budget	Y	N	It is widely used in soil moisture simulation, runoff forecasting, climate change impact analysis, and land use change impact analysis. It is suitable for large-scale land surface process simulation; however the assumptions made by the model conform to the characteristics of plain watersheds, not mountain watersheds.	[22]
UEB	Distributed, Physical	Hourly, Grid	Energy Budget	N	Y	No calibration is required, and fewer variables and parameters are required. However, there is a big difference between the simulated snow internal energy and the monitored data, and there are also problems in the parameterization of the turbulent flux.	[19]

compares the advantages and disadvantages of commonly used snowmelt models.

2.6. Data-Driven Model

The rise of data-driven models is closely related to the availability of data and the development of machine learning algorithms. Compared with statistical models, data-driven models mainly focus on improving the accuracy of prediction, usually using machine learning methods. They usually deal with more complex problems with higher dimensions and more variables [68]. The data-driven approach has the ability to determine which model inputs are critical [69]. At present, data-driven models and statistical models are almost being developed independently [70].

Unlike the process orientation of the physical model, the data-driven model is data-oriented. By learning about the complex relationship between input and output, it has high simulation accuracy even without prior knowledge of the underlying process. This is its biggest advantage. Over the past few decades, the field of hydrology has witnessed the flow and experience of several generations of machine learning methods, from regularized linear regression to support vector regression [71,72], from genetic programming to artificial neural networks (ANNs) [73,74], from classification and regression trees to random forests [75], from Gaussian processes to radial basis function networks [76,77]. All these methods provide useful solutions to different problems, but each method faces its own limitations. Here we will not introduce their principles and algorithms. Instead, we will focus on their application in hydrological studies, especially in the area of snowmelt modeling.

There have been many studies on snowmelt runoff prediction using machine learning. Vafakhah et al. [26] used artificial neural network (ANN) and Adaptive Neuro-Fuzzy Inference System (ANFIS) to simulate the snowmelt runoff in the Taleghan, Alborz. Acar et al. [78] utilized the Multilayer Perceptron (MLP) neural network to simulate the snowmelt runoff in Turkey’s semi-arid climate. Thapa et al. [27] developed a model based on a deep learning long short-term memory (LSTM) network for snowmelt-driven discharge modeling in a Himalayan basin that uses remote sensing snow products as the input, and compared it with Nonlinear autoregressive exogenous (NARX), Gaussian process regression (GPR) and support vector regression (SVR) models. Frederik et al. [79] used the LSTM model to simulate the snowmelt runoff in the New England area, and compared it with the Sacramento model (SAC) and the SNOW-17 model.

Machine learning has been applied to various aspects of snowmelt modeling, such as filling missing observation data [80,81], generating high accuracy rainfall observation data [82], merging satellite precipitation products and measured precipitation [83,84], forecasting real-time rainfall from radar [85], estimating snow cover from satellite data [86], estimating snow water equivalent [87,88], and calculating soil moisture and soil saturated hydraulic conductivity [89,90,91]. In addition, there are also many applications of machine learning in the field of hydrology, such as river water regime simulation [92,93,94], flood hydrograph prediction [95], groundwater dynamics simulation [96], parameter identification of hydrological models [97], and establishment of hydrological evaluation tools [98].

As a black box form of model, the data-driven model has a poor description of the modelling process. It is therefore impossible to evaluate the optimality of the results, and difficult to make meaningful comparisons between the performance of different models. This has caused the application of all data-driven models to be criticized [70,99], though the modelling process was not necessarily incorrect, just not clear enough.

The data-driven model should not act as a competitor to the physical process model, but as a cooperator. The “process-oriented” and “data-oriented” modeling methods are not complementary in themselves, but in many cases, the combination of the two for hydrological modeling is promising. For example, after a simple coupling, we can ensure that the energy and water balance of the entire model remains unchanged, and then determine model parameters in a simple manner through a data-driven method. Meanwhile, physical process models can also add monitoring data to data-driven models [100].

At present, new methods have been developed, such as theoretically guided data science (TGDS) and physical information machine learning. The motivation behind these methods is to improve the physical meaning of machine learning models by combining them with existing scientific knowledge. For example, Herath et al. [101] established a new model induction framework on Genetic Programming (GP), namely the Machine Learning Rainfall-Runoff Model Induction Toolkit (ML-RR-MI). ML-RR-MI is capable of developing fully-fledged lumped conceptual rainfall-runoff models for a watershed using the building blocks of two flexible rainfall-runoff modelling frameworks (FUSE and SUPERFLEX). They extend ML-RR-MI towards inducing semi-distributed rainfall-runoff models. Khandelwal et al. [102] proposed an LSTM-based deep learning architecture that is coupled with Soil and Water Assessment Tool (SWAT). The key idea of their approach is to model auxiliary intermediate processes that connect weather drivers to streamflow, rather than directly mapping runoff from weather variables. The goal of their work is to incorporate our understanding of physical processes and constraints in hydrology into machine learning algorithms. Okkan et al. [103] embedded artificial neural network (ANN) and support vector regression (SVR) into the monthly lumped CRR model. The dynamic water balance model was preferred as the CRR model. Then, all free parameters within these nested hybrid models were calibrated simultaneously. The ML parts within the nested schemes manipulate various output variants derived with three conceptual parameters for monthly runoff simulation. Kalpattern et al. [104] compiled a theoretical method framework that they collectively referred to as “theory-guided data science”.

In addition, observations on hydrology and water science are usually unbalanced in terms of region and time. There is usually a lack of comprehensive training data sets that can be used as machine learning algorithms. Moreover, there is a problem of how to predict situations beyond the scope of the training dataset [99,105]. Global changes are changing hydrology and related cycles, which makes the ability to predict changes beyond previous observations important [106]. Nevertheless, we still believe that machine learning may contribute greatly to hydrological modeling. The integration of hydrological knowledge, process-based models, and machine learning for snowmelt modeling can improve the simulation capabilities.

3. Key Issues in the Snowmelt Model

There are some key issues in snowmelt model that need attention, which include blowing snow, snow on frozen ground, and rain-on-snow.

3.1. Blowing Snow

The uneven distribution of snow caused by blowing snow makes the snow melting process vary considerably in space. The migration and sublimation of snow caused by blowing snow account for a large proportion of the snow hydrological cycle in some areas. MacDonald et al. [107] suggested that the sublimation of snow blowing in the Rocky Mountains caused 17–19% of its annual precipitation. Zhou et al. [108] found that the sublimation of blowing snow in the mountainous areas of western China induced 24% of annual snowfall there. The migration and sublimation of the snow blowing process also bring difficulties to snowmelt modeling. If these factors are not considered in the snowmelt runoff modeling, it will lead to a decrease in the accuracy of watershed runoff simulation.

A series of blowing snow models has been developed and applied to hydrological studies. These models can be roughly divided into two categories: physical blowing snow models based on the wind field (

Table 3. Features of existing physical blowing snow models.

Model	Model Input	Processes Considered	Applicability	References
PBSM	Wind speed, wind direction, air temperature, humidity, land use type	Blowing snow transportation, snow accumulation, snow erosion, blowing snow sublimation	Blowing snow in flat area	[109]
Snowtran-3D	Vegetation type, terrain, air temperature, humidity, precipitation, wind speed, wind direction	Blowing snow transportation, vegetation snow catching, wind speed and terrain changes, snow shear strength, wind-induced surface shear stress, suspended snow, snow erosion, blowing snow sublimation	Blowing snow in Alpine	[112]
SnowDrift-3D	Wind speed, wind direction, terrain, air temperature, humidity, precipitation	Snow transportation, snow accumulation, snow shear strength, wind-induced surface shear stress, snow erosion	Blowing snow in the mountains	[114]
Alpine3D	Wind speed, wind direction, precipitation, solar radiation, soil moisture, vegetation information	Blowing snow transportation, snow accumulation, snow erosion, blowing snow sublimation	Snow process and avalanche on steep mountain surface	[113]

), and conceptual blowing snow models based on terrain drift factors. Pomeroy et al. [109] developed a Prairie Blowing Snow Model (PBSM), which uses meteorological and land-use data to estimate the rotation, suspension and sublimation rate of blowing snow particles. Due to the complex algorithms and data requirements, PBSM is difficult to use. Therefore, Essery et al. [110] simplified the PBSM and used the Mason and Sykes three-Dimensional extension of the Jackson and Hunt theory/version 3R (MS3DJH/3R) model [111] to generate the wind field required for simulation. Their approach reduced the amount of calculation and maximized the reproduction of PBSM’s results. Linston et al. [112] constructed a physics-based numerical snow transport model (Snowtran-3D), which was then used to simulate snow depth evolution in grasslands or high mountain areas with sparse trees, and generate distributed snow depth or snow water equivalent maps. Lehning et al. [113] developed a high-resolution model (Alpine3D) that simulates the snow process on mountain surfaces. This model is composed of a radiation balance model and a drift snow model. The floating snow diffusion equation and the rotating transport equation are used to describe the snow drift process. This model is applicable on very steep terrain. Schneiderbauer et al. [114] developed a high-resolution three-dimensional model of atmospheric snow migration (SnowDrift-3D), which uses a passive transport equation to describe the rotation and suspension of particles. This model is directly connected with the numerical weather prediction model (ALADIN), with an unsteady wind field at a spatial resolution of up to 2 m as driving data.

The physical wind-blown snow model based on the wind field mainly faces three key problems. First, how to define the critical conditions for the occurrence of blowing snow. Studies have shown that blowing snow occurs when wind speed exceeds a certain transportation threshold. This threshold depends on the physical characteristics of the snow surface and environmental meteorological conditions [115]. It is obviously not accurate enough to simply use the critical wind speed to determine the occurrence of blowing snow. Li et al. [116] analyzed the meteorological data on the prairie in Western Canada, and found that the probability of blowing snow is related to wind speed, temperature, and snow age. By combining these three factors, the probability of wind and snow can be predicted. Bowling et al. [117] developed a blowing snow module for the VIC model, which relied on the normal probability distribution of wind speed, snow age, exposed vegetation roughness and temperature to calculate the probability of blowing snow. The second problem is that in complex terrain, especially mountainous areas, there is not enough accurate wind field data to drive the model. Studies have shown that the most sensitive modeling parameter in blowing snow modeling is driving wind field [118,119,120]. It is hard to generate an appropriate wind field simply by spatial interpolation of in-site monitoring data.

Some current solutions are to couple the wind field to an atmospheric model. Bernhardt et al. [121] used the Pennsylvania State University-National Center for Atmospheric Research MM5 model (The fifth-generation Penn State/NCAR Mesoscale Model) [122] to generate a 200 m resolution wind field to drive the SnowTran-3D model, and compared the simulation results to the wind field generated by interpolation. They found this method to be more reasonable. Mott et al. [123] used the ARPS weather model to generate wind fields with different resolutions (50 m, 25 m, 10 m, and 5 m) to drive the Alpine3D model and simulate the snow deposition pattern on the ridges of the Swiss Alps. Vionnet et al. [124] used the 50 m resolution wind field produced by the atmospheric model Meso-NH to drive a blowing snow model.

The third problem is that of data distributed parameterization for large-area blowing snow simulation. The parameterization of wind and snow on the grid requires extensive data of underlying surface, such as surface roughness, which is not easy to be obtained.

One of the problems in the blowing snow model is terrain. Some scholars have recognized that terrain parameters are important predictors of snow depth [125]. Different from the physical blowing snow models based on wind field, some researchers tried to parameterize the redistribution of wind and snow based on terrain analysis. Bloschl et al. [126] reckoned wind drift to be a function of terrain slope and curvature, and used it to interpolate the measured snow water equivalent to initializing the snowmelt model. Purve et al. [127] simulated snow accumulation patterns of a blowing snow model based on the relationship between wind and terrain. Hartrman et al. [128] used TOPMODEL to derive a terrain similarity index, also called humidity index, to redistribute snow from steep slopes to flat terrain. The UEB model also uses terrain parameters to simply describe the effect of wind and snow and simulates the redistribution process of snow by preprocessing the watershed terrain to give the drift factor value of snow. Winstral et al. [129] used digital terrain analysis to quantify the upwind terrain related to wind shielding and exposure, and provided parameters that determine the terrain exposure of each grid unit and the development potential of drift relative to the observed wind. An energy balance snowmelt model (SNOBAL) was then driven, which accurately simulates the observed snow distribution and snowmelt runoff. Winstral et al. [130] added wind speed and vegetation factors to further improve the simulation of snow redistribution. It has been successfully applied to three different watersheds, accurately describing the snow distribution from ridges to valleys. Through comparison, it is found that the snow distribution simulated by the SNOBAL model is more accurate than the model without the redistribution effect of blowing snow. And the simulated surface water input after adding the blowing snow module is consistent with the observed spring runoff pattern.

The physical blowing snow model generally describes the snow transportation process clearly, but the explicit modeling of these complex processes requires a larger amount of observation data and complex computation. The applications of the above-mentioned physical models are largely limited to small study areas. Applications in large areas require simplification. The physical blowing snow model can be applied to simulate a single wind-blown snow event and can also be coupled with the snowmelt model to improve the temporal and spatial resolutions of snowmelt simulation. The core idea of using terrain parameters to simulate snow redistribution is to parameterize the terrain to generate a spatial model in connection with the observed snow pattern. The calculation process is much simpler than that of the physical blowing snow model. It focuses on simulating the redistribution result of snow and does not consider the transportation process of blowing snow. A blowing snow model based on terrain parameters is more suitable for simulating the redistribution of blowing snow at a seasonal scale, but it is difficult to simulate a single blowing snow event. However, a blowing snow model that considers wind speed can also update the snow redistribution in a short time. At present, most snowmelt models do not consider the factor of blowing snow. Therefore, coupling the blowing snow model with the snowmelt model is a scientific issue that needs further research.

3.2. Snow on Frozen Ground

The phase change of water in frozen ground will change the hydraulic and thermodynamic properties of the soil, and is accompanied by absorptions or releases of huge phase change heat which affects the water and heat balance near the surface [131,132,133,134]. When the burial depth of frozen ground is small, its water barrier effect effectively hinders the infiltration of rainwater and snowmelt water, which may promote the formation of surface runoff and soil flow in the surface layer of melting soil and accelerate the response of runoff to rainfall or snowmelt. When there is no frozen ground, the soil has strong infiltration and water storage capacity. The surface runoff is significantly reduced, and the path of soil flow is significantly increased, resulting in a low runoff coefficient and a low response speed of runoff to rainfall or snowmelt water [135,136].

There have been some successful models in the field of water-heat coupling between frozen ground and snow. In 1973, Harlan et al. [137] proposed a water-heat coupling process model of frozen soil and snow. Fleechinger et al. [138] proposed a one-dimensional model of frozen soil water-heat coupling (Simultaneous Heat and Water Model, SHAW), and systematically simulated the land surface process of frozen ground and snow. The model was successfully used to simulate frost depth, snow depth, and soil temperature in an experimental watershed. Zhao et al. [139] used a numerical model (HAWTS) to establish a general parameter correlation equation for estimating snowmelt infiltration in frozen ground to estimate the amount of snowmelt infiltration in the northern forest and grassland environments. Keith et al. [140] used VIC to analyze the freezing and thawing process in the Mississippi River basin. They found that the infiltration of precipitation or snowmelt water would decrease, and the runoff of the watershed would increase with permafrost coverage. Gelfan et al. [141] used one-dimensional partial differential equations to describe the coupled heat transfer process of uniform unsaturated frozen soil based on in situ soil observations, and adopted the influence of phase change on water flow. Daniel Stadler et al. [142] used the SOIL hydrological model [143] to study the infiltration of snow water in the frozen ground layer, and found that the existence of permafrost had a greater impact on the thermal conductivity of the soil. Based on the improvement and fusion of the SOIL model, Jansson et al. [144] formed a Coup model containing multiple sub-modules. This model can simulate water and heat transfer processes, as well as carbon, nitrogen, and other nutrients in the soil-vegetation-atmosphere system. It is also used to simulate soil water and heat process in the tundra of the Heihe alpine meadow and the black soil area of the Songnen Plain [145,146]. Chen et al. [55] combined the SHAW model and the Coup model to establish a distributed hydrothermal coupling model DWHC for the inland river basin in the cold region. This model fully considers the freezing and thawing process of frozen soil and interacts with the MM5 climate model. Zhang et al. [147] coupled the SHAW model to the hydrological model GBHM, and constructed a model, SHAW-DHM, suitable for simulating the processes of water and heat transfer, frozen soil, snow, runoff and confluence in mountainous basins, and applied it to the upper reaches of the Heihe River Basin and Babao River Basin. Cao et al. [148] used a HYDRUS-1D software freeze-thaw module to analyze the influence of freeze-thaw on the hydrological process of active soil in the headwater of the Yellow River.

Studies have shown that explicit consideration of soil freezing and thawing can enhance the simulation capabilities of land surface process and hydrological models [149,150,151]. In the distributed hydrological model, the adjustment of soil freezing and thawing on hydrological processes is mainly reflected in the dynamic adjustment of soil liquid water content and hydraulic conductivity, which in turn affects soil evaporation, infiltration, and runoff processes. A hydrological model that does not include the frozen ground module cannot accurately simulate the snowmelt runoff in the basin with strong frozen soil, as the lack of the frozen ground module caused the model to significantly underestimate the snowmelt flood peak [152]. However, parametrically expressing and quantitatively simulating the influence of frozen ground is still a difficult issue. For different time scales and different regional environments, the relationship between changes in frozen ground and snow melt is still an issue that is worthy of attention. It is necessary to add a frozen ground module to the snowmelt model to study the runoff dynamics in cold regions.

3.3. Rain-on-Snow

Rain-on-Snow (ROS) is a widespread phenomenon in alpine regions around the world, and it plays an important role in generating high flow peaks. Such events are prone to produce mixed rain and snow floods that are more intense than snowmelt floods [153]. Harr [154] found that many peak flows of the Willamette River in Salem, Oregon are related to ROS events. Sui and Koehler [155] attributed the increase in the peak flow of the northern tributary of the Danube in Germany to the increase in ROS events. Merz et al. [156] found that 20% of the flood events in Austria between 1971 and 1997 were driven by ROS. Li et al. [157] found that about 70% of extreme runoff events in the west coast of the United States, the main mountains of the western inland and the Appalachian Mountains had a certain ROS contribution. Furthermore, the 2017 Oroville Dam crisis in California was exacerbated by the rapid snowmelt associated with extreme rainfall [158].

Studies have used observational datasets to study the relationship between ROS events and long-term changes in temperature rise [159,160,161], as well as global and regional climate models to investigate the impact of climate change on ROS events [162,163,164]. These studies demonstrated the correlation between the frequency of ROS events and changes in temperature. ROS-induced peak flow was reduced under climate warming in low and medium altitude areas, while it increased in high altitude areas. With the progress of climate change, ROS may become more widespread and frequent. However, for most areas in the world facing such risks, the runoff-generation mechanism of rain and snow mixed floods based on ROS is still unclear. What is certain is that if snow accumulation and ablation under ROS conditions are not considered, runoff and peak flow may be underestimated [165]. Therefore, a rain and snow mixed runoff model based on ROS is of great significance to the simulation and forecast of the rain and snow mixed flood formation process.

The runoff of ROS events is affected by many pre-conditions and physical processes, such as rainfall intensity and quantity, snow heat, snow type, humidity conditions, water movement of wet snow, interaction of melt water and underlying soil, or land flow of the snow base layer [166,167]. Therefore, a physics-based model is needed to continuously simulate the development of snow cover, as well as transpiration, sublimation and internal processes, in order to estimate the actual runoff generation [160,166]. At present, there are few studies on the manner of runoff under mixed rain and snow. A major difficulty is defining the intensity of ROS. Previous studies varied in the definition of ROS events. MacCabe and Surfleet et al. [160,168] stipulated that when snow is present, rainfall occurs with maximum daily temperature over 0 °C, and snow reduction can be observed. Ye et al. [169] stipulated that it only occurred when at least one of the four precipitation measurements per day is liquid and the ground is covered by snow for ≥1 cm. However, these definitions allow the identification of all possible ROS events, some of which may cause flood events eventually. Sui and Koehler [155] found that most ROS events in southern Germany occurred when the depth of snowmelt is greater than the total rainfall. Based on the reanalysis of the German ROS event in 2011, Kohn et al. [170] set the rainfall threshold to 3 mm, the snow cover of 10 mm, and the snowmelt content of 20% as indicators to identify an ROS event. Masamich et al. [171] defined an ROS day as a day when the daily rainfall exceeds 10 mm and the snow depth exceeds 10 cm. However, there are still problems with the applicability of different indicators, and how to define the intensity of ROS events is still a difficult problem. Moreover, in mountainous areas with complex terrain and sparse monitoring, how to separate rain and snow to determine snow surface rain events and how to parameterize snow surface rain events are also tough tasks.

The enhancement mechanism of ROS on snowmelt runoff remains unclear. It is generally agreed that snow rain has two effects on snowmelt and runoff. On the one hand, ROS increases the water content of the snow, reduces the albedo, increases the energy of the snow, and changes the energy balance of the snow layer [172]. At the same time, ROS enhances the metamorphic effect of snow cover, which leads to the coarsening of snow particle size [173], snow sinking and densification [174], increasing thermal conductivity [172] and thereby increasing the rate at which the snow melts. On the other hand, ROS will quickly saturate the water content of the snow cover, thereby generating maximum and rapid runoff.

A small number of studies conducted field simulation experiments of artificial rainfall on the snow surface. Three artificial rainfall experiments with different intensities were carried out in the alpine region of Austria and found that snowmelt runoff often flowed quickly, following the infiltration path of rainfall. Although the thick layer of snow has a strong water storage capacity, the icy layer in the snow has a greater water barrier effect, and the runoff movement rate caused by heavy rainfall is much greater than that of conventional snowmelt runoff [167,175]. Juras et al. [176] conducted artificial rainfall ROS experiments in the Krykonos Mountains in the Czech Republic, using isotope labeling to explore the contribution of ROS events to runoff, and found that rainwater in ROS events accounted for half of the total flow (52.7%), with the snow cover reserving approximately one-third of the rainwater input (33.6%). Eiriksson et al. [177] used dye tracers to measure ROS runoff on slopes in an artificial rainfall test in Boise National Forest, Idaho, USA. They found that the ROS event produced the most significant slope flow. During an ROS event, the direct contribution of the lateral flow in the snow accounted for about 12% of the total runoff. The lateral flow in the snow may be part of the reason for the rapid water transport during the ROS event and the high runoff flow during the ROS event. If the model uses lateral flow in the snow as a water path mechanism, the hydrological prediction of ROS events and snowmelt events may be significantly improved. However, the mechanism of ROS-based rain and snow mixed runoff generation is still indistinct, which requires further investigation.

In the future, using multi-source remote sensing data, isotopes, changes in surface water reserves and other data to conduct multi-element testing of simulated runoff results may be an important means to improve the accuracy of rain-snow-ice mixed flood simulation.

3.4. Challenges of the Snowmelt Models

Starting from the requirements of snowmelt flood forecasting, we sorted the existing challenges facing runoff generation models for snowmelt flood process simulation into four aspects: data, model, forecast and early warning, and prediction and estimation (

Table 4. Challenges of the snowmelt models.

Aspect	Challenges
1. Data	1.1 High-resolution data (air temperature, precipitation, snow, vegetation, soil moisture, etc.) 1.2 Accurate snow related estimation (snow cover, snow depth, snow water equivalent, etc.) 1.3 Fusion of different data sets and accuracy evaluation
2. Model	2.1 Data assimilation 2.2 The reliability test of runoff simulation 2.3 Process modeling of rain on snow and river ice 2.4 Coupling blowing snow, frozen ground, rain-on-snow 2.5 Coupling with atmospheric model 2.6 Deep learning algorithms
3. Forecast and early warning	3.1 High-resolution, accurate weather forecast 3.2 Uncertainty of weather forecast 3.3 Uncertainty of Snow algorithm 3.4 Robustness of the snowmelt model to floods of different scales
4. Prediction and estimation	4.1 Uncertainty of climate change scenarios 4.2 Uncertainty of downscaling method 4.3 New phenomena such as ROS under rapid heating Carbon black, aerosol, etc.

).

Snowmelt occurs mostly in mountainous areas in high latitudes, and the acquisition of various monitoring data is still a big problem. In particular, the acquisition of snow-related data with high temporal and spatial resolutions has serious limitations. Current automated measurements include sonic depth measurement, snow pillows, Gamma radiation measurement [178], and cosmic ray neutron measurement [179]. These automated methods have higher temporal resolution but lack important spatial variability due to their limited spatial footprints. LiDAR measurement [180] and GPS interferometric reflectometry (GPS-IR) [181] can measure snow cover in a small area, but cannot provide regularly repeated observations. Microwave remote sensing, InSAR techniques in particular, provides effective means to monitor snow cover in a large area, but how to balance the spatial and temporal resolution, and how to retrieve snow parameters in a complex terrain environment still require continuous exploration and research. At present, downscaling technology, data fusion methods, and deep learning algorithms provide us with a powerful way to obtain high-resolution data. In the future, we will need to establish three-dimensional air-space-ground monitoring networks to better monitor and simulate the snow melting process.

In the construction of the snowmelt model, we still need to consider more snow hydrological processes, such as blowing snow, frozen soil, and rain-on-snow mentioned above. In addition, the model needs to be designed to be more flexible and easier to couple with other models, such as climate models, land surface models, and hydraulic models, so that the model can simulate the snow melting process more reasonably and accurately.

The main purpose of snowmelt modeling is to provide more accurate snowmelt prediction simulation for early flood warning. Using the results of numerical weather forecasting as the input of the snowmelt model can help us predict the snowmelt runoff in advance [182]. At present, the most advanced numerical weather prediction models in the world include the MM5 model jointly developed by Pennsylvania State University and the U.S. Center for Atmospheric Research [122], the NCEP and WRF models developed by the National Environmental Center of the United States [183,184], the European mid-term ECMWF model developed by the Weather Forecast Center [185], and GRAPES developed by the China Meteorological Administration [186], etc. Different models usually provide several modes, with different forecast time-space resolution and forecast time effect. For example, the ECMWF model has a forecast time effect of up to 15 days. Its high-resolution TL639L62 mode aims for a 0–10 day forecast and has a horizontal resolution of 30 km, while the lower resolution TL319L62 mode aims for a 10–15 day forecast with a horizontal resolution of 65 km. The reliability and uncertainty of numerical weather forecasting are major problems faced by flood forecasting. We should further study the construction of high-resolution numerical weather prediction to reduce the uncertainty and improve the accuracy of flood forecasting. At the same time, we are required to test the robustness of the snowmelt model to deal with floods of different scales, and better deal with the uncertainty of future snowmelt flood modelling.

Snow cover responds quickly to climate change. Under the background of global warming, snowmelt flood risk rises accordingly [187]. Assessing the risk of snowmelt floods under future climate change scenarios is an urgent task. A major method of understanding the impact of climate change on floods is to couple climate models with hydrological models, and use future meteorological data to drive hydrological models to predict future flood risks [188]. Future climate scenarios are a reasonable description of the temporal and spatial distribution of future climate states based on certain driving forces and scientific assumptions. This can be roughly divided into two categories: emission scenarios and radiative forcing (socio-economic scenarios). Many global climate models have been developed so far. They can give the overall trend of future climate change, however, there are considerable differences among the models. Uncertainties broadly exist in climate simulation, for example the uncertainty of the model itself, the uncertainty of the emission scenario, and the uncertainty of the regionalization technology [189]. In addition to the uncertainties caused by the climate model, in the context of future climate change, new phenomena such as rain erosion and snow caused by rapid warming, carbon black, aerosols, etc., will also induce uncertainties to snowmelt models. The enhancement has also brought greater challenges to the future prediction of snowmelt floods.

4. Summary

Aiming at the current development of runoff generation models in snowmelt flood process simulation, this paper reviewed the characteristics and improvements of different snowmelt models and their application in snowmelt flood simulation. We analyzed the key issues faced by the simulation of runoff generation that are embedded in the current snowmelt flood formation process simulation. The main conclusions are as follows:

(1) The physical distributed snowmelt model has high temporal and spatial resolutions, and strict physical meaning. The spatial scheme of the energy balance component and the hydrological parameters of the physical snowmelt model, and the research on the spatial scheme of forcing data, require further progress. This type of model will be an important tool and a key development object for the simulation of the snowmelt flood process as computing power increases in the future.

(2) Machine learning adds new tools to hydrological modeling. The integration of hydrological knowledge, process-based models and machine learning for snowmelt modeling can further improve the simulation capabilities.

(3) At present, most snowmelt models still lack wind-blown snow modules and frozen ground modules. An explicit consideration of the effects of blowing snow and soil freezing and thawing in the snowmelt modeling process is expected to enhance the capabilities of snowmelt runoff simulation in land surface models and hydrological models.

(4) In the context of climate warming, global ROS events have become more and more extensive and intense, bringing great challenges to the snowmelt runoff simulation. At present, the runoff-generation mechanism for rain and snow mixed floods remains unclear. Therefore, the observation and simulation of ROS processes are an important direction in this field.