A strategy to represent impacts of subgrid-scale topography on snow evolution in the Canadian Land Surface Scheme

Abstract

This sensitivity study applies the offline Canadian Land Surface Scheme (CLASS) version 3.6 to simulate snowpack evolution in idealized topography using observations at Likely, British Columbia, Canada over 1 July 2008 to 30 June 2009. A strategy for a subgrid-scale snow (SSS) parameterization is developed to incorporate two key features: ten elevation bands at 100 m intervals to capture air temperature lapse rates, and five slope angles on four aspects to resolve solar radiation impacts on the evolution of snow depth and SWE. Simulations reveal strong elevational dependencies of snow depth and SWE when adjusting temperatures using a moist adiabatic lapse rate with elevation, with 26% peak SWE differences between that at the average elevation versus the mean of the remainder of the elevation bands. Differences in peak SWE on north- and south-facing slopes increase from 3.0 mm at 10° slope to 17.9 mm at 50° slope. When applied to elevation, slope and aspect combinations derived from a high-resolution digital elevation model, elevation dominates the control of peak SWE values. Inclusion of the range of SSS effects into a regional climate model will improve snowpack and hydrological simulations of western North America's snow-dominated, mountainous watersheds.

INTRODUCTION

Many global and regional climate models (GCMs/RCMs) are now available to study the evolution of the climate system including the impacts of climate change on the cryosphere. When coupled with state-of-the-art land surface models, GCMs and RCMs form the most comprehensive tools used to simulate and project the climate on decadal to century timescales. The recent generation of climate models now strives to resolve more complex land surface processes as computational power increases. To do so, their land surface schemes often parameterize subgrid-scale processes for more representative simulations of land surface processes including snowpack evolution and its contribution to runoff generation (Prein and others, 2015; Zhao and Li, 2015; Kauffeldt and others, 2016). This is especially important in western North America (among other regions worldwide) where snowpacks form the largest component of seasonal water storage (Cayan, 1996; Barnett and others, 2005; Hamlet and others, 2005; Lapp and others, 2005; Mote and others, 2005; Mote, 2006; Kang and others, 2014, 2016). Snow also influences the surface energy fluxes, radiative budgets and hydrologic conditions during the winter accumulation and spring ablation seasons (Cayan, 1996; Brown, 2000; Leung and Qian, 2003; Pan and others, 2003; Hernández-Henríquez and others, 2015).

Snow depth and SWE are two key variables quantifying the seasonal snowpack and its potential contribution to snowmelt runoff and flooding (Mote and others, 2003). Accurate simulations of snow depth, SWE and snowmelt runoff, among other hydrologic variables, rely on realistic estimates of the amount and timing of snowfall, the overall surface energy and radiative balances, and other snow-related processes. Due to the high temporal and spatial variability of snow, it is important to incorporate subgrid-scale representations of snow-related processes and properties in models (Liston, 2004). Advancements in simulating subgrid complexity in snow evolved through the parameterization of such factors as snow redistribution by wind (Déry and Yau, 2002; Gordon and others, 2006), snow–vegetation interactions (e.g. Bartlett and Verseghy, 2015), snow-specific surface area (Roy and others, 2013), snow-cover fraction and non-uniform snow distribution (Liston, 2004; Nitta and others, 2014), and topography (slope, aspect, elevation) (Déry and others, 2004). The validity of including one parameterization over another depends on the grid scale and location of analysis. For example, blowing snow processes only contribute significantly to snow distribution at small spatial scales in regions such as Arctic and alpine tundra, ice sheets and ice shelves, prairies and steppes (Déry and Yau, 2001; Gordon and others, 2006; Lenaerts and others, 2010, 2012). Despite recent subgrid-scale snow (SSS) developments, inaccuracies related to the simulation of snowpack evolution persist (e.g. Liston, 2004; Roy and others, 2013). Further analyses of subgrid variations on the topographic control of air temperatures and incoming solar radiation with impacts on snow evolution are thus warranted, particularly when considering the positive snow-albedo feedback on climate change in high latitude and altitude regions (Hernández-Henríquez and others, 2015; Wang and others, 2016).

Recent increases in air temperature have been observed in many of the world's mountainous regions (Pepin and Lundquist, 2008; Ohmura, 2012; Sharma and Déry, 2016). This induces modifications to winter and spring snowpack levels, in turn leading to shifts in seasonal streamflow and a reduction in summer flows (Dettinger and Cayan, 1995; Arnell, 2003; Déry and others, 2016; Kang and others, 2016). Thus, altitude emerges as an important variable to consider since air temperature trends are most sensitive at snow lines (Jefferson, 2011; Wi and others, 2012; Pepin and others, 2015). High elevations may thus respond more strongly to global and regional warming trends than their lower counterparts (Liu and Chen, 2000; Diaz and Eischeid, 2007; Liu and others, 2009; Pepin and others, 2015).

This paper's main objective is to develop a strategy for a subgrid-scale parameterization for a land surface model with consideration of the effects of elevation, slope and aspect on snow evolution in an idealized model setting. Observational data from a meteorological station over 1 year are utilized to force and evaluate the model. Air temperatures are first adjusted with the moist adiabatic lapse rate (MALR) for a range of elevations. Then incoming direct-beam shortwave radiation is modified according to varying slope angles and aspect quadrants to test model sensitivity of snowpack simulations for these new exposures. Finally, the spatio-temporal evolution of SWE at ~75 m resolution considering the combined effects of elevation, slope and aspect yields insights on the dominant topographic controls on snow evolution at the scale of an RCM grid cell. As such, this sensitivity study provides guidance in the development of an SSS parameterization for future implementation into the land surface scheme of the Canadian Regional Climate Model (CRCM) for improved climate simulations and projections of western North America's snow-dominated, mountainous watersheds (Poitras and others, 2011; Šeparović and others, 2013).

MODEL, DATA AND METHODS

Canadian Land Surface Scheme (CLASS)

This study uses the current operational version of the CLASS, version 3.6, which is the standard land surface model for the CRCM (Verseghy, 1991, 2012; Verseghy and others, 1993; Šeparović and others, 2013). CLASS 3.6, a second-generation land surface model (Pitman, 2003), provides physically based energy, momentum and water exchanges between the land surface, a vegetation canopy and the atmosphere. It represents explicitly four vegetation types: needleleaf trees, broadleaf trees, crops and grasses. The soil column is discretized into three standard layers: 0.10, 0.35 and 4.10 m. Additional deeper layers in the multi-layered configuration can also be specified by the user (Verseghy, 1991; Ganji and others, 2017). Soil temperatures vary over time following a one-dimensional heat conduction equation with heat capacities and thermal conductivities dependent on liquid and frozen soil moisture contents and soil texture (Langlois and others, 2014). Heat conduction transfers rely on the quadratic relationship between soil temperatures and depth.

CLASS incorporates snow as a fourth variable-depth bulk ‘soil’ layer with a thermal regime distinct from the underlying soil. It resolves snow density for fresh snow based on empirical functions involving air temperature (Hedstrom and Pomeroy, 1998). Following deposition, snow density increases exponentially with time to a maximum possible value depending on snowpack depth for compaction (Brown and others, 2006). Thermal conductivity is derived from an empirical function of snowpack density (Sturm and others, 1997). Snowpack temperature variation with depth is solved from a quadratic equation, while a linear equation relates the surface heat flux to the snowpack temperature gradient. As such, snow temperature constitutes a prognostic variable that is computed through the same quadratic equation for soils and is resolved iteratively from the surface energy balance. Snow accumulates when surface air temperatures are below a threshold value T _rs of 0 °C during periods of precipitation (although T _rs can be modified and include mixed (solid and liquid) precipitation at warmer temperatures (Brown and others, 2006)). Implementation of T _rs is usually not required when CLASS is coupled to an atmospheric model such as the CRCM that then supplies amounts for each phase of the precipitation. Snowmelt occurs when the snow surface or bulk snow layer temperature reaches 0 °C. In this case, the melting of the snowpack occurs partially or completely through the excess energy and snow temperature is reset to 0 °C. Snowpacks retain liquid water up to 4% (of snowpack weight) prior to drainage. While this approach is consistent with some land surface schemes (e.g. Noah; Niu and others, 2011), others (e.g. the Interaction Soil Biosphere Atmosphere (ISBA) model (Noilhan and Mahfouf, 1996), the European Centre for Medium-Range Weather Forecasts (ECMWF) land surface model (Viterbo and Beljaars, 1995) and the Versatile Integrator of Surface and Atmosphere (VISA) processes model (Yang and Niu, 2003)) parameterize the maximum liquid water retention capacity as a function of density, with increased water retention capacity as snow densifies. Snow meltwater can also refreeze (causing the snow density to increase). This process continues until the entire snowpack becomes isothermal. Thereafter meltwater may pond on the surface and eventually infiltrate or contribute to surface runoff.

Albedo, an important parameter for assessing slope and aspect effects on snow, is modeled in CLASS using empirical exponential decay functions. If the amount of new snowfall at a given time step ⩾0.1 mm, the snow albedo is set to the fresh snow value of 0.84. The background albedo value for old snow (if the melt rate is non-negligible or the snowpack temperature exceeds −0.01 °C) and for new snowfall on bare ground is 0.50. CLASS integrations typically operate on time steps of ⩽30 min as longer values may induce numerical instabilities (Verseghy, 2012). Verseghy (1991), Bartlett and others (2006), Music and others (2009), Marsh and others (2010), Bartlett and Verseghy (2015), Verseghy and others (2017) and Ganji and others (2017) provide further details on the CLASS model structure, its operation, evaluation and recent improvements.

Study area and forcing data

The University of Northern British Columbia's Quesnel River Research Centre (QRRC) is the site of a meteorological station and is part of the Cariboo Alpine Mesonet climate network (CAMnet; Déry and others, 2010). The QRRC site resides in Likely, British Columbia (BC), Canada on the western perimeter of the Cariboo Mountain chain at 744 m a.s.l. with coordinates 52°37′60″ N and 121°35′24″ W. The Cariboo Mountains span 44 150 km² between BC's Interior Plateau and the Rocky Mountain Trench, forming the northern extension of the Columbia Mountains (Sharma and Déry, 2016). The region observes strong orographic forcing owing to the natural obstruction of moisture transport from the Pacific Ocean. Its climate is wetter compared to the Rocky Mountains to the east and the Interior Plateau to the west with annual precipitation ranging from 500 mm in the town of Quesnel up to 2500 mm along the spine of the Cariboo Mountains (Burford and others, 2009; Sharma and Déry, 2016). More than half of the precipitation falls as snow leading to heavy snowpacks that reach peak accumulations of ~900 mm SWE at treeline (Déry and others, 2014). Above treeline (>1700 m), vegetation in the Cariboo Mountains comprises alpine meadows, mosses and lichens, while at mid altitudes, lodgepole pine (with most mature trees decimated by a recent mountain pine beetle outbreak), Engelmann spruce and subalpine fir abound. Valley bottoms are the setting of a unique inland temperate rainforest (ITR) dominated by old-growth red cedars and hemlocks (Sharma and Déry, 2016). While the QRRC resides in the ITR along the north-facing wall of the Quesnel River Valley, the weather station itself remains in an open, flat clearing with short grass as vegetation cover.

CLASS uses seven variables as meteorological forcing, all of which are available at 15 min intervals over 1 July 2008 to 30 June 2009 from QRRC meteorological station values: incoming long and shortwave radiation, 2 m air temperature and relative humidity (with respect to water), 10 m wind speed, atmospheric pressure and precipitation. A heated tipping bucket gauge protected by an Alter shield quantifies precipitation amounts but without discriminating snowfall from rainfall. These measurements are not modified for gauge undercatch given the phase of the precipitation remains unknown (near T _rs = 0 °C) and winds remain light at the QRRC (see Meteorological forcing for and preliminary tests with CLASS section). A standard anemometer records 10 m wind speed and direction, while a radiometer measures the incoming solar and infrared radiation and an ultrasonic sensor gauges snow depth, which is used only for validation purposes. Details of all sensors at this weather station are available in MacLeod and Déry (2007) and Déry and others (2010). CLASS requires the reference height of air temperature and wind speed measurements to be higher than the vegetation height. This is especially important where canopies exist with heights >20 m. For stations with bare soil or short grass vegetation types (e.g. QRRC), the air temperature and wind reference heights are 2 and 10 m above ground, respectively, sufficiently above the maximum vegetation height.

Given the high frequency of data availability, CLASS is forced at 15 min time intervals with simulation output at the same frequency in addition to daily averages or totals. While air temperature and incoming shortwave radiation are modified to account for subgrid-scale topographic effects, all other CLASS forcing variables remain unchanged (see CLASS simulations section). Analyses focus on CLASS-simulated snow depth and SWE owing to their importance to BC's hydrological cycle. Peak SWE results are analyzed as a measure of total seasonal accumulation, which typically occurs annually between mid-April to mid-May in the study area (Déry and others, 2014). Prior to comparisons with CLASS output, the observed snow depths recorded by the ultrasonic ranger are adjusted for a 20 cm positive bias that appears in its original time series. This bias correction is supported by the lack of subfreezing precipitation and low shortwave albedo values in early fall 2008. The source of this positive bias in the original snow depth data may arise from vegetation at the site (grasses) and/or an incorrect value for the distance to target in the data logger code.

Incorporation of elevation bands

Partitioning of the land surface into ten different elevation bands centers on the QRRC's altitude of 744 m, hereafter referred to as the ‘mean elevation’. The ten elevation bands, all with the same areal weighting of 0.1, are distributed into 100 m vertical increments from 244 to 1244 m excluding 744 m (Arola and Lettenmaier, 1996). The 15 min observed air temperature data from the QRRC meteorological station are then modified for each elevation band using the MALR, taken as 6.4 °C km⁻¹ (Whiteman, 2000; Marshall and others, 2006), with additional sensitivity tests to other lapse rate values (see CLASS simulations section).

Incorporation of slope angles and aspect ratios

The observed incoming shortwave radiation (K↓, W m⁻²) data are first partitioned into direct-beam (I _b, W m⁻²) and diffuse (D, W m⁻²) components using (DeWalle and Rango, 2008):

(1) $$\displaystyle{D \over {K\!\! \downarrow}} = 1 - 1.2\left( {\displaystyle{{K\!\!\downarrow} \over {I_{\rm q}}}} \right) + 0.13\left( {\displaystyle{{K \!\!\downarrow} \over {I_{\rm q}}}} \right)^2,$$

where I _q (W m⁻²) is the calculated potential incoming solar radiation. From Eqn (1), we then set I _b = K↓ − D, which is then modified to account for different slope angles and aspects. Five slope angles of 10°, 20°, 30°, 40° and 50° are prescribed to idealized surfaces, in combination with four different aspect quadrants: north (N), east (E), south (S) and west (W) for each of those slope angles, for a total of 20 different slope/aspect combinations. Inclusion of topographic effects on I _b follows Li and Lam (2007):

(2) $${\rm \;} I_{{\rm bs}} = \left( {\displaystyle{{I_{\rm b}} \over {\sin (\alpha )}}} \right)(\sin (\alpha )\cos (\beta ) + {\rm cos}(\alpha )\sin (\beta )\cos \left \vert {\gamma - \gamma _{\rm n}} \right \vert ).$$

Here, I _bs (W m⁻²) represents the adjusted incoming direct-beam shortwave radiation for a sloped surface, α denotes the solar elevation angle (radians), β is the slope angle (radians), γ denotes the solar azimuth (radians) and γ _n = 0, π/2, π and 3π/2 for north-, east-, south- and west-facing slopes, respectively (Li and Lam, 2007). Thereafter, I _bs and D are summed to provide K↓ _s, the total incoming solar radiation on a sloped surface as the forcing data in the CLASS simulations.

CLASS simulations

Four sets of CLASS offline simulations in idealized settings are performed to explore the impacts of topographic variations on snow evolution at the QRRC over 1 July 2008 to 30 June 2009. First, CLASS is run over the study period for a horizontal surface at the mean elevation without adjustment of the meteorological forcing data and this baseline simulation is evaluated with mean daily snow depth observations (as SWE observations are unavailable). Sensitivity tests to the T _rs parameter are also performed. Next, CLASS is operated for ten elevation bands with air temperatures adjusted by the MALR while all other forcing variables remain unchanged. Two additional simulations are run with the lapse rates set to 4.4 and 8.4 °C km⁻¹ to test the CLASS model sensitivity to this parameter with the focus on peak seasonal SWE accumulation. In the third set of simulations, CLASS incorporates 20 different slope and aspect combinations whereby observed incoming direct-beam solar radiation is modified following Eqn (2) while all other input variables remain unmodified. Then results for snow depth and SWE averaged either over all ten elevation bands (referred to as the ‘mean of all remaining elevations’ or MAREs) or 20 slopes/aspects (referred to as the ‘mean of all slopes’ or MASs) are compared with those from the baseline simulation.

In the final set of CLASS simulations, the effects of elevation as well as slope angle and orientation are both considered. Here the simulation domain spans an area of 7.8 km × 10.3 km encompassing the QRRC. The Shuttle Radar Topography Mission (SRTM) 3 arc second (56.2 m longitudinally and 92.7 m latitudinally) resolution digital elevation model (DEM; Jarvis and others, 2008) provides grid-scale elevation, from which slope and aspect ratios are retrieved for 15 429 points. Then we partition the DEM data according to seven 100 m elevation bins (centered on 644, 744, 844, 944, 1044, 1144 and 1244 m), five 10° slope angle ranges (centered on 5°, 15°, 25°, 35° and 45°) and four aspect quadrants (N, E, S and W) for a total of 140 topographic combinations. Cells with a 0° slope are assumed to represent flat, open water surfaces (e.g. Quesnel Lake and River; see Déry and others, 2010 or Petticrew and others, 2015) and are excluded from the simulations. We then run CLASS adjusting the observed air temperature and incoming solar radiation for each of the 140 topographic exposures and track the biweekly evolution of SWE for each grid cell in the 80.4 km² domain during spring 2009.

The time step for all model integrations is 15 min and a 1-year spin-up period is applied whereby the meteorological data are used recursively to ensure prognostic soil variables approach equilibrium prior to analysis. The surface at the QRRC is taken as 100% covered by grasses, and the corresponding grassland surface parameters for CLASS are prescribed in all simulations (Table 1).

Table 1. Surface parameters used in the CLASS simulations

RESULTS

Meteorological forcing for and preliminary tests with CLASS

Daily averages or totals of seven meteorological variables driving CLASS show the evolution of atmospheric conditions at the QRRC during 1 July 2008 to 30 June 2009 (Fig. 1). Meteorological conditions exhibit a seasonal pattern typical for western Canada, with high air temperatures in summer, followed by a cooling trend in autumn with the coldest conditions in winter (December and January). Air temperatures then rise in late winter and spring to exceed consistently 0 °C on 3 April onward. Precipitation shows no particular seasonal pattern although more intense events arise with the passage of synoptic storms in December and January. Incoming solar radiation follows the expected annual cycle with minimal values during winter. The 10 m mean daily wind speeds remain light (0.4 m s⁻¹ on average) and surpass 1 m s⁻¹ only on eight occasions, in contrast to conditions in nearby alpine terrain (Déry and others, 2010; Petticrew and others, 2015). These quiescent conditions suggest that gauge undercatch and blowing snow have minimal impact on the observed precipitation and snow depth accumulation, respectively. For the study period overall, the QRRC average air temperature and total precipitation are 3.2 °C and 652.3 mm, respectively, with more than half of the precipitation falling as snow (370 mm w.e.).

Fig. 1. Daily averages and precipitation totals of meteorological data from the QRRC station, 1 July 2008 to 30 June 2009. Shown are (a) incoming shortwave radiation partitioned as direct, diffuse and total, (b) incoming longwave radiation, (c) daily total precipitation, (d) 2 m air temperature, (e) 2 m specific humidity, (f) 10 m wind speed and (g) atmospheric pressure.

In response to the autumn cooling trend, snow begins to accumulate in November, with each storm passage yielding rapid rises in snow depth (Fig. 2). Following periods of accumulation, snow depth declines relatively slowly as the snow settles and compacts. Mean daily snow depth stabilizes at ~80 cm from early February until early April at the onset of melt. As air temperatures exceed the melting point, snowmelt proceeds quickly with complete ablation by early May.

Fig. 2. Daily average snow depth simulated by CLASS with different values of T _rs and observed at the QRRC meteorological station from 1 July 2008 to 30 June 2009.

While the simulated snow depth closely matches observations from November through January, thereafter CLASS underestimates snowpack depth compared to the QRRC observations, leading to early melting of snow. Thus, the CLASS threshold air temperature (T _rs) partitioning precipitation into rainfall or snowfall was increased from 0 °C to preset temperature flags for rain of 2 and 6 °C but the snow simulations showed little improvement (Fig. 2). For example, the QRRC observed maximum daily snow depth of 94 cm occurred on 15 March 2009, while CLASS simulated peak daily snow depths of 86, 91 and 92 cm for T _rs = 0, 2 and 6 °C, respectively, on 2 February 2009. Thus, as recommended by Brown and others (2006), the original (unadjusted) precipitation data with T _rs = 0 °C are used henceforth to run CLASS, and it is assumed that the model underestimates the snow depth and simulates an early melt at the QRRC (as reported for other locations; see Pomeroy and others, 1998; MacDonald and others, 2016).

Effects of elevation

Snow depth and SWE simulated by CLASS forced by the QRRC data between 1 July 2008 and 30 June 2009 exhibit strong elevational dependence (Fig. 3). Indeed, snow depth and SWE increase as a function of elevation with maximum values at the highest elevations. A comparison between the QRRC (mean elevation) simulated snow depth and SWE with the MAREs reveals notable differences (Fig. 4). For example, the mean elevation maximum SWE equals 197.3 mm, whereas the corresponding value for the MAREs achieves 151.6 mm. Similarly, the maximum snow depth for the mean elevation attains 86.4 cm, while for the MAREs it reaches 58.4 cm. These differences persist throughout the year, from the accumulation to the ablation period, due to the inclusion of air temperature change at elevation in MARE values compared to the horizontal surface at the QRRC.

Fig. 3. Elevational dependence of average daily (a) SWE and (b) snow depth, simulated by CLASS from 1 July 2008 to 30 June 2009 when applying the moist adiabatic lapse rate (MALR).

Fig. 4. Elevational dependence of average daily (a) SWE and (b) snow depth, comparing simulated CLASS results for the mean elevation (744 m) with the mean of all remaining elevations (MAREs), from 1 July 2008 to 30 June 2009.

Sensitivity tests to the prescribed lapse rates reveal stronger differences in the simulated maximum SWE at lower than mid or higher elevations (Fig. 5). The simulated snowpack is more sensitive to air temperatures at low elevations where wintertime precipitation phase and snowmelt remain important. At high elevations, the peak annual SWE shows relatively minor differences between the three lapse rates given that air temperatures remain ⩽0 °C throughout winter and the precipitation forcing unchanged from that at lower elevations. Thus, orographic enhancement of precipitation may be a more important factor controlling an elevation dependence in maximum SWE in regions where air temperatures remain subfreezing through the accumulation periods.

Fig. 5. Maximum SWE as a function of elevation for three different lapse rates simulated by CLASS during the winter of 2009.

Effects of slope and aspect

Table 2 and Fig. 6 report peak SWE accumulation values and their timing for different slopes and aspects relative to that for a horizontal surface at the QRRC. Peak SWE for the north-facing slope at 10° attains 198.9 mm on 2 February 2009, whereas for the south-facing slope, peak SWE reaches 195.9 mm on 2 February 2009. With slope angles increasing from 10° to 50°, peak SWE values on north-facing surfaces rise by 11.1 mm, while on south-facing surfaces, they decrease by 3.7 mm owing to altered incoming direct-beam solar radiation. For example, the daily mean shortwave radiation intercepted by a north-facing (south-facing) surface on 1 April at 10° and 50° slope angles attains 76.3 and 55.2 W m⁻² (87.4 and 97.1 W m⁻²), respectively. Table 2 also presents the timing of complete snowmelt on different slopes and aspects. The date for complete snowmelt for the north- and south-facing 10° slopes is 19 April 2009 and 15 April 2009, respectively. With increasing slope angle, the north-facing slope melts more slowly (timing of complete melt for a 50° slope is 7 May 2009), whereas the south-facing surface melts earlier with increasing slope (timing of complete melt is 11 April 2009 for the 50° slope angle). Next for each aspect, the simulated SWE for a horizontal slope angle is compared to the mean SWE of the MASs to assess how simulations differ when using a sole mean slope angle versus a range of subgrid slope angles (Fig. 7). Averaging SWE simulated from different slopes influences the magnitude of peak SWE and the timing of snow disappearance, depending on aspect. This effect is greatest on north-facing slopes where energy input is lowest (Fig. 7a) but lesser on other aspects where energy input is higher (Figs 7b–d). Comparing mean SWE of all aspects for different slope angles to that for a horizontal surface demonstrates the effects of slope on the timing of snowmelt (Fig. 8). Increasing slope angles delay snowmelt, reinforcing the influence of north-facing slopes on idealized snowmelt simulations. Given the scant incoming solar radiation recorded at the QRRC (Fig. 1a), slope and aspect controls on snow evolution become significant only during the ablation period.

Fig. 6. Average daily SWE plots for five slope angles (10°, 20°, 30°, 40° and 50°) with SWE values plotted for comparison between aspects (horizontal, N, E, S, W), simulated by CLASS from 1 July 2008 to 30 June 2009. A similar pattern is seen for snow depth (not shown).

Fig. 7. Average daily SWE simulated by CLASS for a 30° slope angle compared with the mean of all slopes (MASs) for different aspects, 1 July 2008 to 30 June 2009. A similar pattern is seen for snow depth (not shown).

Fig. 8. Average daily SWE for a horizontal surface compared with the mean SWE of all four aspect quadrants plotted for five different slope angles (10°, 20°, 30°, 40° and 50°). SWE is simulated by CLASS from 1 July 2008 to 30 June 2009. A similar pattern is seen for snow depth (not shown).

Table 2. Peak SWE accumulation values (PV, mm), dates of peak accumulation (PD) and complete melt (CM) for different slope angles and aspects in 2009 as simulated by CLASS

Combined effects of elevation, slope and aspect

The final set of simulations incorporates the effects of elevation on air temperatures plus slope and aspect on incoming direct-beam solar radiation over an area of 80.4 km², which approaches that in the current generation of RCMs. The area encompasses the QRRC in the valley of the Quesnel River, just to the west of its outlet from Quesnel Lake with higher elevations to its southwest (Fig. 9a). During spring 2009, topography has a dominant control on SWE with the most abundant accumulations at the highest elevations and the least amounts along the valley floor (Figs 9b–f). This pattern persists over time, but by April, the effects of slope and aspect emerge. For instance, on 15 April, areas on the north-facing slope along the Quesnel River Valley show higher SWE values than exposed areas just to the south at higher elevations (Fig. 9e). With the prolonged days in April, the snow ablates rapidly such that only bare ground or vegetation remains by the end of the month.

Fig. 9. (a) Topographic map of a 7.8 km × 10.3 km area encompassing the QRRC with elevations inferred from a 3 arc second resolution DEM. The CLASS simulated SWE considering elevation, slope and aspect on (b) 1 March, (c) 15 March, (d) 30 March, (e) 15 April and (f) 30 April 2009. Areas in white denote water surfaces or absence of snow.

SUMMARY AND CONCLUSIONS

This paper presents preliminary results on the development of a SSS parameterization for CLASS simulations to explore the effects of varying elevations, slopes and aspects on snow depth and SWE evolution in idealized settings. CLASS was first validated against observational data at the QRRC weather station in Likely, BC for the period 1 July 2008 to 30 June 2009, revealing negative biases in the simulated snowpack accumulation and an early onset of melt. In the next set of simulations, ten elevation bands 100 m apart vertically spanning 244–1244 m (excluding the mean elevation of 744 m) were implemented and observed air temperatures were modified using the MALR. Differences of 28 cm in peak snow depth and 45.7 mm in peak SWE emerged between the mean elevation versus the MAREs. In the subsequent CLASS simulations, the incoming direct-beam shortwave radiation was adjusted for slope angles ranging from 10° to 50° at 10° intervals on four different aspects (north, east, south and west) and results compared to those for a flat surface. Peak SWE values not only rose in magnitude (up to 11.1 mm increase) on north-facing slopes but their timing of occurrence was also delayed by nearly 3 weeks. Averaged simulated SWE for a combination of different tilt angles and aspects compared to values on a horizontal plane revealed the importance of subgrid-scale topography on peak accumulation and the timing of snowmelt, particularly for north-facing aspects. Finally, slope angles and aspect ratios were combined with elevation bands, based on topographic bins generated from a high-resolution DEM, to illustrate the spatio-temporal effects of complex topography on snowpack evolution at the subgrid scale of a typical RCM grid cell. Overall, these results revealed distinct trends in snow depth and SWE with elevation, indicating that elevation has more influence than slope and aspect angles when comparing spatial averages. Indeed, representing snowpacks using only mean (model grid cell) topographic characteristics masks the non-linear effects elevation, slope and aspect introduce in their evolution through time. As such, consideration of SSS effects in relatively coarse-scale climate models, particularly in regions of complex topography such as western North America, remains a priority for accurate simulations of basin-scale cryospheric and hydrologic processes.

Additional idealized simulations with the SSS parameterization in CLASS are being undertaken to further assess the robustness of the preliminary results reported here. CLASS experiments will be extended by several years to test the sensitivity of the simulated snowpack to interannual variations in climatic conditions. This will reveal whether subgrid topographic effects on snow evolution are more important in some years over others (or remain relatively constant over time) and the atmospheric controls on such variations (or lack thereof). Finally, CLASS model performance will be further explored by using SWE data from nearby snow pillow stations. This extended investigation will assist in demonstrating the importance of incorporating SSS effects within the coupled CRCM/CLASS framework, leading to an improved representation of snowpacks, snowmelt, and streamflow generation across the mountainous terrain of western North America.