Abstract

Shallow (<1 m deep) snowpacks on agricultural areas are an important hydrological component in many countries, which determines how much meltwater is potentially available for overland flow, causing soil erosion and flooding at the end of winter. Therefore, it is important to understand the development of shallow snowpacks in a spatially distributed manner. This study combined field observations with spatially distributed snow modelling using the UEBGrid model, for three consecutive winters (2013–2015) in southern Norway. Model performance was evaluated by comparing the spatially distributed snow water equivalent (SWE) measurements over time with the simulated SWE. UEBGrid replicated SWE development at catchment scale with satisfactory accuracy for the three winters. The different calibration approaches which were necessary for winters 2013 and 2015 showed the delicacy of modelling the change in shallow snowpacks. Especially the refreezing of meltwater and prevention of runoff and infiltration of meltwater by frozen soils and ice layers can make simulations of shallow snowpacks challenging.

INTRODUCTION

In many countries, at low altitude, snowpacks in agricultural areas are typically shallow (i.e., <1 m deep). Shallow snowpacks are much more sensitive to changes in air temperature and fluctuations between precipitation as rain or snow than deep snowpacks and have a reduced ability to dampen peak outflows compared with deep snowpacks (Wever et al. 2014). Furthermore, snowpacks less than 20 cm deep cannot insulate the soil from low air temperatures, but can amplify deep-freezing of soils due to an increase in ground albedo caused by the snow (Zhao et al. 2013). This deep-freezing of the soil profile can change the soil hydraulic properties dramatically (Stähli et al. 2001; Al-Houri et al. 2009). In southern Norway this causes severe soil erosion through impeded infiltration of snowmelt and rainwater at the end of winter (Øygarden 2003). The amount of snow at the end of winter determines how much meltwater is available for creating overland flow and soil erosion on agricultural soils (Lundberg et al. 2016). Besides damaging the nutrient rich top layer of agricultural soils, sediment transfer from fields to streams during winter and spring represents a major part of the annual loss of phosphorus and nitrogen from agricultural catchments (Su et al. 2011). Therefore, it is important to better understand the development of shallow snowpacks during winter, including the spatial distribution of snow water equivalent (SWE).

Snow is one of the most changeable elements in the hydrological cycle and poses a great range of challenges for monitoring and measuring (Doesken & Robinson 2009). Below, the main factors influencing snow processes of interest for this study are summarised. A more detailed description can be found in Gray & Male (1981) and Lundberg et al. (2016).

Catchment properties such as slope, aspect and vegetation cover have a major influence on the condition and morphology of snowpack. Slope, curvature and aspect determine exposure to wind (Lapen & Martz 1996), precipitation and solar radiation (Gray & Male 1981). Wind interacting with terrain and vegetation heavily influences the deposition and redistribution of snow, enhancing or reducing snow accumulation (Gascoin et al. 2013). More snow is usually accumulated on leeward slopes and forest edges than on windward sides, due to reduced exposure to higher wind speeds. Snow depth is often greater in depressions due to reduced wind speeds and increased turbulence. Furthermore, vegetation (e.g., crop residues) enhances the roughness of the surface during the exposed period before total snow cover, resulting in increased turbulent wind flow and causing complex snow structures and patterns. Another important factor is the air temperature at the time of snowfall, which controls the dryness, hardness and crystalline form of new snow, and thereby its erodibility by wind. For instance, heavy wet snow falling at temperatures close to 0 °C is less susceptible to wind transport and is common for low altitude areas in the temperate zone (Richard & Brun 2008), like southern Norway, where this study was conducted.

The main process responsible for differential melting of snow in a catchment is energy exchange. Important energy inputs and outputs are solar radiation, rain, incoming long-wave radiation, outgoing long-wave radiation and sensible and latent heat fluxes such as evaporation, condensation and sublimation of ice (Tarboton et al. 2000).

The complexity of snow accumulation, distribution and melt processes and the speed of change, which is typical for shallow snowpacks, make spatial investigation of snow properties challenging and costly (Gascoin et al. 2013). Therefore, the spatial resolution of snow-related data is often poor, or no data on spatial distribution are available (Tarboton et al. 2000). Distributed snowpack models, using physically based one-dimensional mass and energy balance, have been successfully used for describing spatial variability in snowpack properties (Luce & Tarboton 2001). However, most studies have focused on snowmelt modelling in mountainous, high altitude areas with usually deep snowpacks (e.g., Blöschl et al. 1991; Cline et al. 1998; Hiemstra et al. 2006; Bernhard et al. 2010; Zwaaftink et al. 2011) and forested areas (Lundberg et al. 2016). Little research has been done on modelling shallow snowpack processes with high spatial and time resolution for agricultural areas. These areas are often neglected due to the focus of most studies on water recharge (for hydropower and drinking water) by snow and ice (Lundberg et al. 2016). The effect of snowpacks on freeze-thaw and water transport in agricultural soils has received great attention (e.g., Flerchinger 2000; Parkin et al. 2013); however, most of these studies were conducted at point scale. Some studies have investigated distributed soil erosion induced by snowmelt in agricultural areas, but have either focused on a single snowmelt event (Ollesch et al. 2005, 2006; Weigert & Schmidt 2005) or could not produce satisfactory results for the whole winter period due to lack of incorporation of snowpack development in the model (Grønsten & Lundekvam 2006).

The need to better mitigate soil erosion and peak discharge during winter and spring from agricultural areas requires a better understanding of the temporal and spatial changes of shallow snowpacks and that spatial distributed snow models can simulate these changes, to provide continuous datasets, e.g., for soil erosion models. Therefore, the aims of this study were to: (1) quantify and visualise the temporal and spatial changes of shallow snowpacks (during three winter periods) on arable land; and (2) test a spatially distributed snow model on observed snow cover patterns. For this aim, the UEBGrid model (Sen Gupta & Tarboton 2013) was applied.

METHODS

Study area and weather data

The study area, Gryteland (Figure 1), is located in the Ski municipality, approximately 30 km south of Oslo, Norway. The catchment covers an area of 0.29 km2. In this area, snow cover processes have been monitored since 2008, with the main focus on arable land. A monitoring station was installed at the outlet of the sub-catchment in 2008 (Kramer & Stolte 2009). This station measures precipitation, runoff and the suspended solids content by means of automated sampling. The catchment is characterised by undulating landscape (altitude: 106–141 m a.s.l., slope 2–10%) covered by approximately 60% arable land and 40% coniferous forest. All aspects are represented except east-facing slopes.

Figure 1

Slope map for the arable land of the Gryteland catchment, showing the locations of the snow poles and the SWE sampling point transect. The six SWE sampling points used for the calibration of the model are named in the map.

Figure 1

Slope map for the arable land of the Gryteland catchment, showing the locations of the snow poles and the SWE sampling point transect. The six SWE sampling points used for the calibration of the model are named in the map.

Mean annual temperature is 5.3 °C, with minimum average daily temperatures of −4.8 °C in January/February and maximum average daily temperatures of 16.1 °C in July. Mean annual precipitation is 785 mm, with a minimum monthly amount of 35 mm in February and a maximum of 100 mm in October (Thue-Hansen & Grimenes 2015). Winter is usually relatively unstable, with alternating periods of freezing and thawing and several snowmelt events (Kværnø & Øygarden 2006). Snow data have been collected in the area since 2009 and show maximum snow depth of around 55–60 cm in the years 2009–2011 (Stolte & Kværnø 2013).

As input for modelling and for understanding the snow accumulation and melting processes, detailed measured weather data were necessary. A weather station (106 m a.s.l) was installed in the catchment at the end of 2013, providing data for the winters 2014 and 2015. Data from a station (96 m a.s.l.) 6 km from the catchment were used for winter 2013 (Thue-Hansen & Grimenes 2015).

Field measurements

During the three winters, SWE was sampled after weather changes were expected to result in changes in SWE. Three to five samples were collected at each measuring point, depending on the snow depth. These samples were combined to measure the mass of the snow and to calculate its density (average mass divided by average snow volume). SWE was calculated by dividing the product of average snow density and average depth of the snowpack by the density of water (1,000 kg m−3). The samples were taken along transects that covered all aspects and slope angles occurring in the catchment (Figure 1). To get a better understanding of the spatial development of the snow cover, snow poles (Figure 1) marked with a measuring scale were installed for winter 2015. These poles allowed snow depth to be measured from a distance by reading the value on the poles with a spyglass, keeping the disturbance of the snowpack to a minimum. Due to the strong correlation of snow depth with SWE, it has been recommended to substitute the time-consuming SWE measurements with snow depth measurements (Jonas et al. 2009).

On each measuring occasion, the SWE at each pole was estimated by multiplying the measured snow depth (at the pole) by the average snow density measured along the snow sampling transect. These spatially distributed pole measurements were then compared with the SWE maps produced by the model. In addition, aerial pictures were taken with a drone to visualise snow patterns in the catchment.

UEBGrid

In southern Norway, the point version of the Utah Energy Balance snow model (UEB) (Tarboton & Luce 1996; You 2004; Luce & Tarboton 2010; Mahat & Tarboton 2012; Mahat et al. 2013) has been used successfully (Kramer & Stolte 2009). To extend its application in this area for spatial distributed modelling, the grid version of UEB (UEBGrid) (Sen Gupta & Tarboton 2013, 2016; Brown et al. 2014) was used in this study. It models the snow pack as a single layer to avoid over-parameterisation (Sen Gupta & Tarboton 2013). The snow pack is characterised by state variables, SWE, energy content (Es) and the age of the snow, which is only used for albedo calculations. Es is defined as the energy content of the snowpack plus a top layer of soil with depth (i.e., thermal active depth of the soil (De)). Soil temperature measurements at different depths in the catchment showed that major fluctuations usually happened between 10 and 40 cm depth (depending on the winter). Therefore, De, which interacts with the snow, was initially set to 10 cm (Table 1). Es represents energy content including this soil depth. High frequency oscillating ground heat fluxes above De are absorbed into Es. This procedure provides a simple approximation of the effects of frozen and unfrozen ground on the snowpack.

Table 1

Input parameters used for the three simulations (winter 2013, 2014, 2015)

Parameter Same for all three simulations 2013 2014 2015 
Tr [°C]    
Ts [°C] −1    
z0 [m] 0.01    
df [–]    
Ksat [m h−120    
Liquid holding capacity of snow (Lh) [−] 0.05    
Average atmospheric pressure [Pa] 101,325    
Canopy coverage fraction map values [−]    
Canopy height map value [m]    
Subtype map value [−] Ground (0)    
De [m] 0.1    
Bristow–Campbell B parameter for January [°C]  5.8 3.1 5.0 
Bristow–Campbell B parameter for February [°C]  7.1 3.1 4.9 
Bristow–Campbell B parameter for March [°C]  10.9 – – 
Bristow–Campbell B parameter for April [°C]  9.2 – – 
Longitude [degrees] 10.84    
Latitude [degrees] 59.67    
Parameter Same for all three simulations 2013 2014 2015 
Tr [°C]    
Ts [°C] −1    
z0 [m] 0.01    
df [–]    
Ksat [m h−120    
Liquid holding capacity of snow (Lh) [−] 0.05    
Average atmospheric pressure [Pa] 101,325    
Canopy coverage fraction map values [−]    
Canopy height map value [m]    
Subtype map value [−] Ground (0)    
De [m] 0.1    
Bristow–Campbell B parameter for January [°C]  5.8 3.1 5.0 
Bristow–Campbell B parameter for February [°C]  7.1 3.1 4.9 
Bristow–Campbell B parameter for March [°C]  10.9 – – 
Bristow–Campbell B parameter for April [°C]  9.2 – – 
Longitude [degrees] 10.84    
Latitude [degrees] 59.67    
Es of the snow pack is determined by incoming long-wave radiation, advected heat from precipitation, ground heat flux, sensible heat flux, latent heat flux due to sublimation/condensation, advected heat removed by meltwater and outgoing long-wave radiation. The change in SWE depends on the rainfall rate (Pr), snowfall rate (Ps), meltwater outflow (Mr) and sublimation from the snowpack (E) (Equation (1)):  
formula
(1)
Ps is determined according to Equation (2):  
formula
(2)
where, Ps = precipitation falling in the form of snow [m h−1]; P= total (measured) precipitation [m h−1]; Pr = precipitation falling in the form of rain [m h−1]; df= drift factor [−].
In Equation (2), the amount of Pr is determined by the threshold temperatures for which all precipitation is rain (Tr) or snow (Ts). According to the following relations (Equation (3)), precipitation is partitioned based on the measured air temperature (T):  
formula
(3)
To initiate snowmelt in UEBGrid, the temperature of the snowpack (Ta) has to reach 0 °C. By using Es as a state variable the model does not explicitly forecast the temperature of the snowpack. Ta is therefore obtained from the following equations (Equations (4)–(6)):  
formula
(4)
 
formula
(5)
 
formula
(6)
where, pw = density of water [1,000 kg m−3]; C = specific heat of ice [2.09 kJ kg−1°C−1]; pg =soil density [kg m−3]; Cg =specific heat of soil [1.4 kJ kg−1 °C−1]; hf =heat of fusion [333.5 kJ kg−1]; Cw = specific heat of water [4.18 kJ kg−1 °C−1].
According to the described equations, Es determines the liquid content of the snowpack. Together with Darcy's law for flow through porous media, Mr is determined (Equation (7)):  
formula
(7)
where, Ksat = snow saturated hydraulic conductivity [m h−1]; S= relative saturation in excess of water retained by capillary forces [–].

The model is driven by input of air temperature, precipitation, wind speed, humidity and solar radiation at time steps sufficient to resolve the diurnal cycle (six hours or less). Sen Gupta & Tarboton (2013) reorganised the input–output workflow of UEB to use network Common Data Form map files (netCDF) to allow the model to be run over a grid with input varying in space and time, resulting in the UEBGrid model used in this study. The model runs separately at each grid cell and the melt outputs are accumulated for predefined sub-watersheds.

Measured weather data (precipitation, wind speed, solar radiation, air temperature, humidity) at hourly time intervals were used here as input for the model (Figures 2 and 3). For the different years, the monthly Bristow–Campbell B parameter was given as an input (Table 1) for the model to calculate atmospheric transmissivity. Required maps for watershed delineation (Figure 1), slope and aspect were derived in ArcMap 10.2.2 from a digital elevation model (DEM), created from light detection and ranging data (LIDAR) (Starkloff & Stolte 2014). The calculations for the agricultural area were done with a raster resolution of 5 m × 5 m. The forest was not further considered in this study, because in all three winters the snow depth was negligible.

Figure 2

Measured air temperature and precipitation in winter 2013, 2014 and 2015. The grey rectangles, with start and end dates, indicate the periods with measurable snow cover.

Figure 2

Measured air temperature and precipitation in winter 2013, 2014 and 2015. The grey rectangles, with start and end dates, indicate the periods with measurable snow cover.

Figure 3

Diagrams of measured wind speed and precipitation per hour, with marked precipitation as snow (S) and rain (R) and wind rises of wind speed and direction in winter 2013, 2014 and 2015. The radial durations are in hours with the wind speeds [m s−1] in each cardinal direction covering 30 degrees. For illustration purposes, a horizontal line, marking wind speeds of 4 m s−1, was included in the diagrams.

Figure 3

Diagrams of measured wind speed and precipitation per hour, with marked precipitation as snow (S) and rain (R) and wind rises of wind speed and direction in winter 2013, 2014 and 2015. The radial durations are in hours with the wind speeds [m s−1] in each cardinal direction covering 30 degrees. For illustration purposes, a horizontal line, marking wind speeds of 4 m s−1, was included in the diagrams.

UEBGrid allows definition of grid cells for which detailed output can be produced. In this study, grid cells containing the location of a snow pole or SWE measuring point were selected for output, resulting in a total of 41 points. Six points were chosen for calibration and are discussed in detail in this paper. These six points represent the major characteristics of the catchment, with a north facing slope (‘North 1’), a south facing slope (‘South 1’), a northwest facing plain (‘Point 12’), two depressions (‘Station 1’ and ‘Station 2’) and a point influenced by the forest (‘Station 3’) (Figure 1).

The spatially distributed snow observations of 2015 allowed evaluation of the performance of the model for all locations in the study area. Therefore, the best calibration results for 2015 were spatially and temporally compared with the measured SWE measurements. To statistically evaluate the spatial performance of UEBGrid for 2015, the Nash–Sutcliffe model efficiency coefficient (NSE) was calculated for each measuring occasion.

For calibrating the model, the user guide of UEBGrid (Tarboton & Luce 1996) suggests adjusting the parameters surface aerodynamic roughness height (z0), surface head conductance (Ks), Ksat and new snow visible band reflectance (avo). A previous study carried out with the point version of UEB for the same area as used in this study (Kramer & Stolte 2009) showed that the parameters Ks and avo had no significant effect on the modelling results and suggested that correction of the threshold temperatures Tr and Ts might be necessary. In a later study, Tarboton et al. (2000) suggested adjusting the drift factor (df) to better account for snow transport by wind. Additionally, due to the importance of the thermal interaction of shallow snowpacks with the underlying soils (Norum et al. 1976), the parameter De was chosen as a calibration parameter. Considering this, the parameters Ksat, z0, df, Tr, Ts, and De were used for calibrating the model. The initial input parameters used for this study are presented in Table 1.

The parameter z0 gives the height below which the wind speed is zero in the logarithmic wind profile layer of the model. It effects the calculation of the heat conductance (Kh) between snow surface and air (Equation (8)) and the vapour conductance (Ke) at the snow surface (Equation (9)):  
formula
(8)
where, k = von Kármán's constant [0.4]; V = wind speed [m s−1]; z = measuring height of wind speed [m].  
formula
(9)

The df parameter is used in the model to assign different cells a multiplication factor (Equation (2)) to account for spatial differences in SWE due to snow drift by wind.

RESULTS AND DISCUSSION

Weather data

Measured air temperature, precipitation, and wind speed and direction are presented in Figures 2 and 3, respectively. Solar radiation and relative humidity are not presented, because no significant difference was observed between the years, which means that differences in snow cover properties were not caused by these factors.

As can be seen in Figure 2, the periods with measurable snow cover had different durations in different years.

Wind speed and wind direction are important factors for the distribution and shape of snowpacks, especially during snowfall and immediately after deposition. Already wind speed of about 3 m s−1 can redistribute light and dry snow crystals (Gray & Male 1981). In Figure 3, measured wind speed together with measured precipitation for the snow measuring periods in the three winters are presented. As can be seen, the highest (hourly average; 8 m s−1) measured wind speed occurred in winter 2015, possibly resulting in the redistribution and densification of deposited snow. For comparison, an average typical blizzard occurring in the USA has an average wind speed of 10 m s−1. Furthermore, during the main snowfall episodes in 2015, the wind speeds exceeded 4 m s−1, reaching 5–6 m s−1. In comparison, during winter 2013 and 2014, the wind speeds reached only 3.5 m s−1 during snowfall.

These measurements confirmed the observations in the field that wind strongly shaped the snowpack in winter 2015, which had to be considered in the modelling.

Field measurements

Measured snow depth and SWE for the three winter periods are presented in Figure 4, which shows results from the six locations used for model calibration and validation.

Figure 4

Measured snow water equivalent (SWE) and snow depth at the same six locations in winter 2013, 2014 and 2015. Circles mark measurements at Station 2, which were set to zero because running meltwater on snow and ice made correct measuring impossible. The last measuring points for 2014 and 2015 show the dates when snow had disappeared at all measuring locations. However, some locations were snow-free earlier, but no exact dates could be given due to rapid melting.

Figure 4

Measured snow water equivalent (SWE) and snow depth at the same six locations in winter 2013, 2014 and 2015. Circles mark measurements at Station 2, which were set to zero because running meltwater on snow and ice made correct measuring impossible. The last measuring points for 2014 and 2015 show the dates when snow had disappeared at all measuring locations. However, some locations were snow-free earlier, but no exact dates could be given due to rapid melting.

The smallest snow depths and SWEs measured at Station 2 in 2013 were set to zero because running meltwater on snow and ice made it impossible to measure during these dates (Figure 4). The measurements at Station 1 after the second measuring occasion in winter 2015 were not representative because of a thick (10–20 cm) ice layer on top of the soil, that was established after snowmelt followed by refreezing. Any SWE or snow depth measurements carried out at this location afterwards represent measurements on top of this ice layer until the end of winter.

The development of snow depth and SWE differed between the three winters (Figure 4). In winter 2013, there was an early decrease in snow depth and SWE due to a melting event at the end of January (Figure 2), followed by a sudden increase in snow depth and SWE, and a constant snow depth and slow increase of SWE over time until the spring melt. In total, 10 measuring occasions were possible in 2013.

The short winter of 2014 resulted in only four measuring occasions (note: the unusually high first snow depth value at Station 1 might be a measuring mistake). While the snow depth decreased already after the second measuring occasion, SWE continued to increase until all snow melted rapidly after the fourth measuring occasion. An increase in SWE followed by rapid snowmelt is typical behaviour of snowpack before each melting episode (Gray & Male 1981). As the snow depth decreases, the snow density increases due to retention of water (meltwater, often supplied by rain) in the snowpack, until the snowpack reaches a point where it is saturated with water (generally 3–5% water by weight) and the internal temperature of the pack reaches 0 °C (Gray & Male 1981). When the snowpack has reached this point, i.e., when it is ‘ripe’, then all additional water is quickly released and the snowpack collapses. In 2014 no major precipitation events (Figure 2) occurred before all the snow was gone. The air temperature (Figure 2) stayed below 0 °C for almost the whole period with snow, until melting started (measuring occasions 3 and 4).

During winter 2015, snow cover was established for a longer period than in 2014, resulting in seven measuring occasions. Similar to winter 2013, winter 2015 was more complex than winter 2014 with an extensive melting event before sampling was started, followed by an increase in snow depth and SWE and an additional snowmelt event on 24–26 January. Meltwater that had accumulated in the depression froze during the following cold period beginning on 26 January, which resulted in thick ice layers in the depressions (Figure 5).

Figure 5

Aerial picture of snowmelt pattern in the study area, 20 February 2015. The white rectangle shows the area covered by the small photos of ice layers, taken 28 January and 3 March 2015.

Figure 5

Aerial picture of snowmelt pattern in the study area, 20 February 2015. The white rectangle shows the area covered by the small photos of ice layers, taken 28 January and 3 March 2015.

Figure 5 shows an image of the arable part of the study area patched together from aerial pictures taken at the end of the melting period in 2015. In large parts of the area, all the snow had melted except in the south.

The aerial images (Figure 5) visualised spatial patterns of the snow cover. From these patterns, it was possible to obtain information about the processes determining the distribution of snow at catchment scale. Examples are: (1) ice development due to refreezing of accumulated meltwater in depressions, which were the main pathways for surface runoff; (2) the closed snow cover in the southern parts of the catchment, on arable land, showing the shadowing effect of the forest against solar radiation and deposition of snow along the forest edge due to reduced wind speeds.

Modelling

The modelling results for winter 2013 are presented in Figure 6. Changing the parameters df, lh, Ksat, Tr and Ts did not improve the modelling results. The best possible fit was achieved with a slight change of z0 from 0.01 to 0.006 m and the adjustment of De to 0.3 m (Figure 6, black line).

Figure 6

Measured and modelled snow water equivalent (SWE) for winter 2013. The arrows mark false zero measurements (3 April), as discussed in the text. Lines ‘SWE 0.3 m z0.006’, ‘− z0.007’ and ‘− z0.01’ present results of model runs with different values for z0.

Figure 6

Measured and modelled snow water equivalent (SWE) for winter 2013. The arrows mark false zero measurements (3 April), as discussed in the text. Lines ‘SWE 0.3 m z0.006’, ‘− z0.007’ and ‘− z0.01’ present results of model runs with different values for z0.

However, it should be noted that similar fits were achieved by setting the parameters lh, Ts and Tr to values of 0.8, 3 and 3.8 °C, respectively. Snowfall at temperatures around 3 °C has been observed in the study area; however, setting Ts to 3 °C alone did not change the modelling results, and the highest measured lh in soaked snow cannot greatly exceed 0.15 (Kinar & Pomeroy 2015).

That the adjustment of De was necessary to achieve better modelling results for 2013 indicated that the interaction of the snowpack with the underlying soil was a major process occurring in winter 2013. Observations carried out during winter 2013 supported this assumption. Winter 2013 was characterised by a major melting event during the measuring period (between 30 January and 6 February) (Figure 2). During this event an (approximately 10 cm thick) ice layer, which had developed earlier in January on the soils, prevented infiltration of meltwater (image 3 in Figure 7), resulting in high water contents in the bottom layer of the snowpack. During snow sampling on 30 January, it was observed that a slight disturbance of the snow resulted in a total collapse of the lower 5 cm snow into water (Figure 7). De directly influences Ta of the snowpack in UEBGrid, as can be seen from Equations (4) and (6). An increasing of De results in a decrease of Ta and therefore in a reduction of the liquid content of the snowpack and a slowdown of melting.

Figure 7

Photos 1, 2 and 4 taken at different measurement points on 3 April 2013. Photo 3 shows water standing on an ice layer on top of the soil during snow sampling on 30 January 2013.

Figure 7

Photos 1, 2 and 4 taken at different measurement points on 3 April 2013. Photo 3 shows water standing on an ice layer on top of the soil during snow sampling on 30 January 2013.

A slight overprediction for the North 1 and Station 1 locations can be seen, but not exceeding 2 cm difference between measured and modelled SWE. For these locations, the initial z0 value of 0.01 m resulted in a better fit (Figure 6). The model predicted a longer period of snow cover for Stations 1 and 3 and Point 12 compared with the measured values (Figure 6, arrow markings). However, these measurements require some explanation. As can be seen from Figure 7, there was still snow present on 3 April 2013.

There was no snow on the south facing slopes on 3 April (image 1 in Figure 7), where the model predicted the timing of total exposure correctly. The north facing slopes were still covered with snow (image 1), which the model predicted correctly. At Point 12 (image 4), Station 3 (image 2) and Station 1 (not shown), the snow had disappeared in areas where the stubble had broken the snow surface (therefore zero measurements at these locations), but not in the wheel tracks, which had no stubble (image 4 in Figure 7). In comparison, neighbouring areas (not shown) without any plant cover or covered with only poorly developed winter wheat were still completely covered with snow on 9 April. The model was unable to account for the change in albedo and increased heat conductance due to stubble. For bare soil, the model would probably have correctly predicted the timing of total disappearance of snow at all locations. The last measurement points in Figure 6 indicate when all snow had disappeared at all locations (3 April at South 1 and 15 April at the other locations), which fits well with the timing predicted by the model.

The modelling results for 2014 are presented in Figure 8. Calibration was done by adjusting z0 visually from 0.01 to 0.006. Increasing De from 0.1 to 0.3 m, as used for the simulation of winter 2013, did not improve the modelling results.

Figure 8

Measured and modelled snow water equivalent (SWE) for winter 2014. Line ‘SWE 0.1 m z0.01’ presents the model results with De set to 0.1 m and z0 set to 0.01 and line ‘SWE 0.1 m z0.006’ the results with z0 set to 0.006. Line ‘SWE 0.3 m z0.006’ presents modelling results with the same input as used for winter 2013.

Figure 8

Measured and modelled snow water equivalent (SWE) for winter 2014. Line ‘SWE 0.1 m z0.01’ presents the model results with De set to 0.1 m and z0 set to 0.01 and line ‘SWE 0.1 m z0.006’ the results with z0 set to 0.006. Line ‘SWE 0.3 m z0.006’ presents modelling results with the same input as used for winter 2013.

However, it should be noted that the limited number of measurements in 2014 made detailed evaluation of the model results difficult, but it can be seen that the modelled SWEs had the same magnitude as the measured values and the model simulated the measured increases and decreases in SWE (Figure 8). Winter 2014 was the least complex of the three observed winters, which was also mirrored by the simple calibration necessary.

For winter 2015, no satisfactory fit was achieved by adjusting De and z0 to values used for winter 2013 and 2014 (Figure 9, dotted line). This indicated that other processes, which did not occur or had only insignificant effects in 2014 and 2013, needed to be taken into account for the 2015 simulations. That wind played a greater role in winter 2015, compared with the other winters, was confirmed by the wind measurements (Figure 3). In addition, during the conducted field measurements it was observed that the following interacting processes occurred: (1) different snow patterns due to the redistribution of snow by wind, depending on the terrain (slope, curvature and aspect); (2) the forest acting as an obstacle to wind (reducing wind speeds, causing snow to be deposited along the forest edge); (3) the shading effect of the forest against solar radiation. More snow could accumulate in the southern part of the catchment due to process 2, because the forest edge protected the area against the major wind direction of southwest to southeast in 2015 (Figure 3), contrary to 2014 and 2013, where major wind directions were north to southeast (Figure 3). These interacting processes cannot be directly simulated with UEBGrid or, indeed, with most available snow models (Essery et al. 2013). However, in UEBGrid the drift factor (df) can be used to account for the effect of processes 1 and 2 (Tarboton et al. 2000; Luce & Tarboton 2001). Accordingly, the df (Figure 8) was calibrated for the six locations (Figure 8). As can be seen, depending on the location, different values for df were necessary to fit the modelled to the measured SWE.

Figure 9

Measured and modelled snow water equivalent (SWE) for winter 2015. Lines ‘SWE drift 1.2’, ‘− 1.4’ and ‘− 1.5’ present results of model runs with different values for df. Line ‘SWE 0.3 m z0.006’ presents the model run with the same input as used for 2013.

Figure 9

Measured and modelled snow water equivalent (SWE) for winter 2015. Lines ‘SWE drift 1.2’, ‘− 1.4’ and ‘− 1.5’ present results of model runs with different values for df. Line ‘SWE 0.3 m z0.006’ presents the model run with the same input as used for 2013.

With these adjustments, a satisfying fit was achieved. However, the complete disappearance of snow was two days late for the South 1 location, while for North 1, Station 3 and Point 12 this was one day too early. A slight overprediction occurred at all locations, especially after the fifth measuring occasion, with a maximum difference of approximately 0.03 m for the seventh measuring occasion at the North 1 location.

That the highest df value (1.5) had to be applied to the areas close to the forest edge suggested that process 2 played a major role in winter 2015. Also, the south facing slopes required a value of 1.5, suggesting increased deposition of snow on these slopes during winter 2015. Based on this calibration, a map with different values for df, depending on the location, was created (Figure 10), enabling an evaluation of the model performance for the whole area. The highest df value (1.5) was given to the areas close to the forest edge (similar to Station 3) and south facing slopes (similar to South 1). Areas with similar characteristics to Point 12 and Station 2 were given a value of 1.4, while north facing slopes with similar characteristics to the points North 1 and Station 1 were given a value of 1.2.

Figure 10

Aspect and drift factor maps used as an input for the 2015 modelling.

Figure 10

Aspect and drift factor maps used as an input for the 2015 modelling.

Figure 11 shows the modelled SWE maps for the seven measuring occasions in 2015. In addition, for each date the modelled SWE is compared with the measured SWE at the snow pole locations.

Figure 11

Modelled SWE maps for the seven measuring occasions in 2015 and the corresponding graphs with measured SWE vs. modelled SWE at the snow poles. In the graphs, the modelled points are connected for better visibility. In each graph the calculated NSE value is presented and snow poles with the largest differences between measured and simulated SWE (i.e., >0.005 m SWE) are marked in the graphs. Additionally, the location of these poles is shown in the simulated SWE maps, where black dots indicate overprediction and white dots underprediction of SWE by the model.

Figure 11

Modelled SWE maps for the seven measuring occasions in 2015 and the corresponding graphs with measured SWE vs. modelled SWE at the snow poles. In the graphs, the modelled points are connected for better visibility. In each graph the calculated NSE value is presented and snow poles with the largest differences between measured and simulated SWE (i.e., >0.005 m SWE) are marked in the graphs. Additionally, the location of these poles is shown in the simulated SWE maps, where black dots indicate overprediction and white dots underprediction of SWE by the model.

To better quantify the performance of the model, for each measuring occasion, the calculated NSE is presented in the graphs of Figure 11. Except for 27 January, the model performed well, with NSE ranging between 0.55 and 0.73. The NSE for 27 January was 0.06, indicating that the model was only slightly better than the average of the measured data. As can be seen from Figure 11, the model predicted lower SWE at many locations compared with the measurements. The measurements taken on 27 January were characterised by a decrease in snow depth, but an increase in SWE compared with 23 January (Figure 4). This suggests that the snow pack was actually maintaining more water from precipitation (snow and rain at 1.5 °C), which had occurred on 26 January (Figure 2), than UEBGrid anticipated. The parameter lh can be used in the model to account for an increased water storage of the snowpack. However, the given value for lh is constant for the whole modelling period, which would have affected all results.

Based on the modelling experiences gathered in this study, four suggestions on how the UEBGrid model could be improved were found. (1) To improve the prediction of different snow patterns due to the combination of wind speed, wind direction and terrain, the df could be linked to slope angle, aspect and input of measured wind directions. (2) To account for accumulation of meltwater in snowpacks at the bottom of slopes and in depressions, a flow routing routine could be implemented in the model. This would also help in the context of soil erosion prediction due to snowmelt. (3) To allow for ice layers within snowpacks and on soil surfaces, theoretically, using a multiple layered snow pack and (4) including the possibility of a soil profile below the snow in the model (as in the case of some non-distributed snow models) could improve the simulation results. However, suggestions 3 and 4 could result in a difficult to calibrate, overparameterised and measurement-intensive model, reducing its applicability. Therefore, when improving the model, a good balance between model complexity and input data requirements, which UEBGrid provides, has to be maintained (Avanzi et al. 2016).

Overall, the UEBGrid model proved capable of simulating the development and spatial distribution of SWE in the shallow snowpacks in the three winter periods. However, the different calibration approaches that were necessary for winter 2013 and 2015 showed the delicacy of modelling the change in shallow snowpacks. Due to the rapid changes occurring in such snowpacks, as shown by the observations and measurements, a certain simulation error has to be expected. Especially, the refreezing of meltwater and prevention of runoff and infiltration of meltwater by frozen soils and ice layers can make simulations of shallow snowpacks challenging. Simulation results therefore have to be treated with care, and correct results can probably only be acquired in combination with detailed field observations and measurements. Many additional tools for improving the observation and measurements of rapidly changing snowpacks already exist (e.g., Egli et al. 2009; Hopkinson et al. 2012; Garvelmann et al. 2013; Kinar & Pomeroy 2015; De Michele et al. 2016), which should be more widely utilized.

CONCLUSIONS

The three winter periods studied differed in length, amount of snow, and processes shaping the snow cover. These interacting processes, caused by different weather conditions and state of the soil, resulted in higher (2013, 2015) or lower (2014) complexity of the modelling. By considering weather data and quantitative (measurements) and qualitative (aerial and ground images) field observations, it was possible to correctly validate the model results. Especially, the importance of detailed photographic image and video documentation, as a substitute for quantitative measurements, for process understanding cannot be stressed enough.

UEBGrid proved capable of simulating the development and spatial distribution of SWE in the study area for 2013, 2014 and 2015. The acquired calibration values, especially for df, should be tested in coming winters and in other areas with similar characteristics, to evaluate how consistent they are. UEBGrid captured the essential physics of accumulation and melt processes of shallow and patchy snow covers. The measurements and calibrations underlined how sensitive such snowpacks are to wind, small changes in the composition of precipitation and soil temperature conditions, requiring a greater focus on automated, differentiating precipitation and spatially distributed wind and SWE measurements. Also, more research into processes related to shallow snowpacks is necessary, with the prospect of decreasing snow depths in areas where snowpacks now are usually deep, due to global warming (Barnett et al. 2005).

However, the fact that it was possible to satisfactorily simulate temporal and spatial changes of shallow snowpacks made us confident that we are a step closer to the overall goal: to better quantify and finally reduce soil erosion during winter and spring.

ACKNOWLEDGEMENTS

The research was funded by the SIS-Catchy project of NIBIO, Norwegian Institute of Bioeconomy Research.

REFERENCES

REFERENCES
Al-Houri
,
Z. M.
,
Barber
,
M. E.
,
Yonge
,
D. R.
,
Ullman
,
J. L.
&
Beutel
,
M. W.
2009
Impacts of frozen soils on the performance of infiltration treatment facilities
.
Cold Regions Science and Technology
59
,
51
57
.
Avanzi
,
F.
,
De Michele
,
C.
,
Morin
,
S.
,
Carmagnola
,
C. M.
,
Ghezzi
,
A.
&
Lejeune
,
Y.
2016
Model complexity and data requirements in snow hydrology seeking a balance in practical applications
.
Hydrological Processes
30
,
2106
2118
.
Bernhard
,
M.
,
Liston
,
G. E.
,
Strasser
,
U.
,
Zängl
,
G.
&
Schulz
,
K.
2010
High resolution modelling of snow transport in complex terrain using downscaled MM5 wind fields
.
The Cryosphere
4
,
99
113
.
Blöschl
,
G.
,
Kirnbauer
,
R.
&
Gutknecht
,
D.
1991
Distributed snowmelt simulations in an alpine catchment 1. Model evaluation on the basis of snow cover patterns
.
Water Resources Research
27
(
12
),
3171
3179
.
Brown
,
M. E.
,
Racoviteanu
,
A. E.
,
Tarboton
,
D. G.
,
Sen Gupta
,
A.
,
Nigro
,
J.
,
Policelli
,
F.
,
Habib
,
S.
,
Tokay
,
M.
,
Shrestha
,
M. S.
,
Bajracharya
,
S.
,
Hummel
,
P.
,
Gray
,
M.
,
Duda
,
P.
,
Zaitchik
,
B.
,
Mahat
,
V.
,
Artan
,
G.
&
Tokar
,
S.
2014
An integrated modelling system for estimating glacier and snow melt driven streamflow from remote sensing and earth system data products in the Himalaya
.
Journal of Hydrology
519
,
1859
1869
.
De Michele
,
C.
,
Avanzi
,
F.
,
Passoni
,
D.
,
Barzaghi
,
R.
,
Pinto
,
L.
,
Dosso
,
P.
,
Ghezzi
,
A.
,
Gianatti
,
R.
&
Vedova
,
G. D.
2016
Using a fixed-wing UAS to map snow depth distribution: an evaluation at peak accumulation
.
The Cryosphere
10
,
511
522
.
Doesken
,
N. J.
&
Robinson
,
D. A.
2009
The challenge of snow measurements
. In:
Historical Climate Variability and Impacts in North America
(L. A. Dupigny-Giroux & C. J. Mock, eds)
.
Springer
,
Heidelberg
.
Egli
,
L.
,
Jonas
,
T.
&
Meister
,
R.
2009
Comparison of different automatic methods for estimating snow water equivalent
.
Cold Regions and Technology
57
,
107
115
.
Essery
,
R.
,
Morin
,
S.
,
Lejeune
,
Y.
&
Menard
,
C. B.
2013
A comparison of 1701 snow models using observations from an alpine site
.
Advances in Water Resources
55
,
131
148
.
Flerchinger
,
G. N.
2000
The Simultaneous Heat and Water (SHAW) Model: Technical Documentation. Technical Report NWRC-2000-09
,
USDA-ARS Northwest Watershed Research Center
,
Boise, Idaho
.
Garvelmann
,
J.
,
Pohl
,
S.
&
Weiler
,
M.
2013
From observation to the quantification of snow processes with a time-lapse camera network
.
Hydrology and Earth System Sciences
17
,
1415
1429
.
Gascoin
,
S.
,
Lhermitte
,
S.
,
Kinnard
,
C.
,
Bortels
,
K.
&
Liston
,
G. E.
2013
Wind effects on snow cover in Pascua-Lama, Dry Andes of Chile
.
Advances in Water Resources
55
,
25
39
.
Gray
,
D. M.
&
Male
,
D. H.
1981
Handbook of Snow
.
Pergamon Press
,
Toronto
,
Canada
.
Hiemstra
,
C. A.
,
Liston
,
G. E.
&
Reiners
,
W. A.
2006
Observing, modelling and validating snow redistribution by wind in a Wyoming upper treeline landscape
.
Ecological Modelling
197
,
35
51
.
Hopkinson
,
C.
,
Collins
,
T.
,
Anderson
,
A.
,
Pomeroy
,
J.
&
Spooner
,
I.
2012
Spatial snow depth assessment using LiDAR transect samples and public GIS data layers in the Elbow River watershed, Alberta
.
Canadian Water Resource Journal
37
,
69
87
.
Jonas
,
T.
,
Marty
,
C.
&
Magnusson
,
J.
2009
Estimating the snow water equivalent from snow depth measurements in the Swiss Alps
.
Journal of Hydrology
378
,
161
167
.
Kinar
,
N. J.
&
Pomeroy
,
J. W.
2015
Measurement of the physical properties of the snowpack
.
Reviews of Geophysics
53
.
doi:10.1002/2015RG000481
.
Kramer
,
G. J.
&
Stolte
,
J.
2009
Cold-Season Hydrology Modeling in the Skuterud Catchment. Bioforsk Report
,
Norway
,
Vol. 4, No. 126
.
Luce
,
C. H.
&
Tarboton
,
D. G.
2001
Modeling Snowmelt over Area: Modeling Subgrid Scale Heterogeneity in Distributed Model Elements
.
Civil and Environmental Engineering, Utah State University
,
USA
.
Lundberg
,
A.
,
Ala-Aho
,
P.
,
Eklo
,
O.
,
Klöve
,
B.
,
Kværner
,
J.
&
Stumpp
,
C.
2016
Snow and frost: implications for spatiotemporal infiltration patterns – a review
.
Hydrological Processes
30
,
1230
1250
.
Mahat
,
V.
&
Tarboton
,
D. G.
2012
Canopy radiation transmission for an energy balance snowmelt model
.
Water Resources Research
48
,
W01534
.
Mahat
,
V.
,
Tarboton
,
D. G.
&
Molotch
,
N. P.
2013
Testing above- and below-canopy representations of turbulent fluxes in an energy snowmelt model
.
Water Resources Research
49
,
1107
1122
.
Norum
,
D. I.
,
Gray
,
D. M.
&
Male
,
D. H.
1976
Melt of shallow prairie snowpacks: basis for a physical model
.
Canadian Agricultural Engineering
18
,
2
6
.
Ollesch
,
G.
,
Sukhanovski
,
Y.
,
Kistner
,
I.
,
Rode
,
M.
&
Meissner
,
R.
2005
Characterization and modelling of the spatial heterogeneity of snowmelt erosion
.
Earth Surface Processes and Landforms
30
,
197
211
.
Ollesch
,
G.
,
Kistner
,
I.
,
Meissner
,
R.
&
Lindenschmidt
,
K.
2006
Modelling of snowmelt erosion and sediment yield in a small low-mountain catchment in Germany
.
Catena
68
,
161
176
.
Parkin
,
G.
,
von Bertoldi
,
A. P.
&
McCoy
,
A. J.
2013
Effect of tillage on soil water content and temperatures under freeze-thaw conditions
.
Vadose Zone Journal
12
,
1
9
.
Richard
,
A.
&
Brun
,
E.
2008
Snow and Climate Physical Processes, Surface Energy Exchange and Modeling
.
Cambridge University Press
,
New York
.
Sen Gupta
,
A.
&
Tarboton
,
D. G.
2013
Using the Utah energy balance snow melt model to quantify snow and glacier melt in the Himalayan Region
. In:
Proceedings of the Western Snow Conference 81st Annual Meeting
.
Adaptive Water Management in a Changing Climate
,
Jackson Hole, Wyoming
,
15–18 April
, pp.
103
114
.
Stähli
,
M.
,
Nyberg
,
L.
,
Mellander
,
P.
,
Jansson
,
P.
&
Bishop
,
K. H.
2001
Soil frost effects on soil water and runoff dynamics along a boreal transect: 2. Simulations
.
Hydrological Processes
15
,
927
941
.
Stolte
,
J.
,
Kværnø
,
S. H.
2013
Snowmelt and runoff in two small field-scale catchments
. In:
Agriculture and Environment – Long Term Monitoring in Norway
(
Bechmann
,
M.
&
Deelstra
,
J.
, eds).
Akademika Publishing
,
Trondheim
,
Norway
.
Su
,
J. J.
,
van Bochove
,
F.
,
Thériault
,
G.
,
Novotma
,
B.
,
Khaldoune
,
J.
,
Denault
,
J. T.
,
Zhou
,
J.
,
Nolin
,
M. C.
,
Hu
,
C. X.
&
Bernier
,
M.
2011
Effects of snowmelt on phosphorus and sediment losses from agricultural watersheds in Eastern Canada
.
Agricultural Water Management
98
,
867
876
.
Tarboton
,
D. G.
&
Luce
,
C. H.
1996
Utah Energy Balance Snow Accumulation and Melt Model (UEB), Computer Model Technical Description and Users Guide
.
Utah Water Research Laboratory and USDA Forest Service Intermountain Research Station
,
USA
.
Tarboton
,
D. G.
,
Blöschl
,
G.
,
Cooley
,
K.
,
Kirnbauer
,
R.
,
Luce
,
C.
2000
Spatial snow cover processes at Kühtai and Reynolds Creek
. In:
Spatial Patterns in Catchment Hydrology: Observations and Modelling
(
Grayson
,
R.
&
Blöschl
,
G.
, eds).
Cambridge University Press
,
Cambridge
,
UK
.
Thue-Hansen
,
V.
&
Grimenes
,
A. A.
2015
Meteorologiske Data for Ås 2014
,
Norwegian University of Life Science
,
Norway
.
Weigert
,
A.
&
Schmidt
,
J.
2005
Water transport under winter conditions
.
Catena
64
,
193
208
.
Wever
,
N.
,
Jonas
,
T.
,
Fierz
,
C.
&
Lehning
,
M.
2014
Model simulation of the modulating effect of the snow cover in a rain-on-snow event
.
Hydrology and Earth System Sciences
18
,
4657
4669
.
You
,
J.
2004
Snow Hydrology: The Parameterization of Subgrid Processes Within a Physically Based Snow Energy and Mass Balance Model
.
PhD Thesis
,
Civil and Environmental Engineering, Utah State University
,
Logan
,
USA
.
Zhao
,
Y.
,
Huang
,
M. B.
,
Horton
,
R.
,
Liu
,
F.
,
Peth
,
S.
&
Horn
,
R.
2013
Influence of winter grazing on water and heat flow in seasonal frozen soil of Inner Mongolia
.
Vadose Zone Journal
12
.
doi:10.2136/vzj2012.0059
.
Zwaaftink
,
C. D. G.
,
Löwe
,
H.
,
Mott
,
R.
,
Bavay
,
M.
&
Lehning
,
M.
2011
Drifting snow sublimation: a high-resolution 3-D model with temperature and moisture feedbacks
.
Journal of Geophysical Research
116
,
D16107
.