Study site

The study site is located in a montane forest opening in the Henry’s Lake Mountains of southwestern Montana, approximately 15 km west of West Yellowstone, Montana, USA (44°42′ N, 111°17′ W). Slope angle ranges between 24° and 29° and slope aspect ranges between 34° and 49°. The surrounding forest is composed primarily of Subalpine fir (Abies lasiocarpa) and Douglas fir (Pseudotsuga menziesii), including two prominent trees at the base of the slope which protrude into the slope’s viewshed. The study site itself is sparsely vegetated by grasses, forbs, small shrubs and occasional small coniferous trees and sagebrush. Sheltered from westerly ridge-top winds and facing northeast, the study site frequently develops surface hoar in winter, the focus of several past studies (Reference Landry, Birkeland, Hansen, Borkowski, Brown and AspinallLandry and others, 2004; Reference Logan, Birkeland, Kronholm and HansenLogan and others, 2007). Reference LutzLutz (2009) provides a detailed site description. The study site is located on an unconfined slope, approximately 100 m below the ridge top and 100 m above the drainage floor, where a small seasonal (snowmelt-fed) creek exists in spring and summer. No obvious vapor source or confining topographic feature exists on or adjacent to the site.

Overview of field observations

Winter field observations documented the thickness and shear strength of a buried surface hoar layer that caused snowpack instabilities during the late winter season of 2005. The surface hoar formed during the third week of January 2005 and persisted on the surface for an additional 11 days before it was buried by two small snowfall events in early February. On 28 February 2005, 3 weeks after surface hoar burial, we sampled a cross-shaped plot spanning the entire site (Fig. 1, plot 1). This allowed the quantification of slope-scale trends of snow properties at the beginning of the study. Subsequently on 1, 8, 14 and 21 March we sampled the square corner plots (plots 2–5), allowing us to quantify spatial patterns and changes over time.

Fig. 1. Sample layout. Dashed lines separate plots 1–5 (indicated by circled numbers). SMP profiles are circles; shear frame tests are black diamonds.

During the five sampling days, we collected 824 SnowMicroPen (SMP) profiles and 352 shear frame tests, described in detail in the following subsections. The observations were located on a sample layout optimized for geostatistical analysis of SMP profiles and pit-based regression analysis. An effective sample layout captures the existing spatial structures of a natural phenomenon. This particular layout was developed by Reference LutzLutz (2009) and was found to be more effective than other gridded layouts tested by Reference Kronholm and BirkelandKronholm and Birkeland (2007). Although Reference Kronholm and BirkelandKronholm and Birkeland (2007) found randomly located observations produce the most reliable estimates, these prove impractical to map in the field and to coordinate spatially coincident observations (SMP and shear frame tests).

At each plot, we nested a majority of the observations in nine clusters, each including six shear frame tests and nine SMP profiles within a 1 m² surface area (Fig. 1). We refer to these clusters as ‘pits’, since observers prepared recessed benches for shear frame testing. Previous spatial studies utilizing the SMP centralized all the closely spaced observations in one group near the middle of a plot (e.g. Reference Birkeland, Kronholm, Logan and GanjuBirkeland and others, 2004; Reference Kronholm, Schneebeli and SchweizerKronholm and others, 2004; Reference Logan, Birkeland, Kronholm and HansenLogan and others, 2007). Alternatively, by distributing the closely spaced observations into nine decentralized clusters, the variance structure at short distances is representative of the entire plot (instead of just the central portion).

We established a slope-oriented coordinate system anchored to wooden posts and rebar that had been fixed on the site perimeter before snowfall. We located sample placements by pulling and offsetting survey tapes between the wooden posts. A Geographic Information System (GIS) that linked sample coordinates with high-resolution terrain and vegetation models allowed the spatial reconstruction of snowpack data.

SnowMicroPen (SMP) measurements and derived surface hoar thickness

We derived surface hoar thickness, h _wl, from penetrometer profiles recorded using the SMP (Reference Schneebeli and JohnsonSchneebeli and Johnson, 1998; Reference Johnson and SchneebeliJohnson and Schneebeli, 1999) Serial No. 9-9. Previous work identified an inverse relationship between surface hoar thickness and shear strength (Reference Davis, Jamieson and JohnstonDavis and others, 1998; Reference Jamieson and JohnstonJamieson and Johnston, 1999; Reference Jamieson and SchweizerJamieson and Schweizer, 2000). The sensor resolution and tip size enable thin layers such as surface hoar to be effectively detected (Reference Pielmeier and SchneebeliPielmeier and Schneebeli, 2003). Microstructural information pertinent to snowpack stability (Reference Lutz, Birkeland and MarshallLutz and others, 2009; Reference Pielmeier and MarshallPielmeier and Marshall, 2009) was derived using a micro-mechanical model (Reference Marshall and JohnsonMarshall and Johnson, 2009) and is presented by Reference LutzLutz (2009).

We assessed and classified the signal quality of all profiles, as defined by Reference LutzLutz (2009). Signal problems, such as micro-variance dampening and artificial drifting, influence signal interpretation and should be accounted for when analyzing SMP signals (Reference LutzLutz, 2009; Reference Pielmeier and MarshallPielmeier and Marshall, 2009). Two main signal criteria assisted in manually delineating h _wl in SMP profiles: (1) the load and fracture signature typical of large hoar crystals interacting with the sensor tip (Reference Johnson and SchneebeliJohnson and Schneebeli, 1999; Reference Schneebeli, Pielmeier and JohnsonSchneebeli and others, 1999) and (2) short (e.g. 0.1–3 mm) ‘flat line’ segments that exhibit no obvious autocorrelation and possess near-zero mean resistance values. Since the adjacent layers were composed of small facets that produced smaller and more frequent loads in the resistance signal, these criteria could be safely assumed.

Shear frame measurements

Shear frame measurements quantified the in situ shear strength of the targeted buried surface hoar. The shear frame test has been used extensively for avalanche research and forecasting (Schleiss and Schleiss, 1970). The presence of additional buried surface hoar layers in the soft slab above the targeted layer precluded the use of column-type stability tests.

Following standard operating procedure, a 250 cm² shear frame was positioned within 5 mm of the boundary between the surface hoar and the superstratum (Perla and Beck, 1983; Reference SommerfeldSommerfeld, 1984; Reference GreeneGreene and others, 2009). Glide-planeparallel force was applied in the downslope direction to a force gauge (attached to the shear frame) at a pull rate within the range of brittle fracture (e.g. 1s (Reference Jamieson and JohnstonJamieson and Johnston, 2001)). The shear strength, τ ₂₅₀ (Pa), is defined as the pull-generated shear force at failure F _fail divided by the tested area A _frame,

(1)

We approximated the shear component of force exerted by the experimental set-up, τ _expt., which is dependent on slope angle, α, and the mass of the shear frame, m _frame, and contained snow, m _expt.snow, the latter of which was approximated from slab density measurements such that

(2)

The total shear stress at failure, τ_∞ , was then estimated by multiplying the sum of τ ₂₅₀ and τ _expt by a size-correction parameter (Reference SommerfeldSommerfeld, 1984) such that

(3)

Terrain and vegetation observations

We surveyed topography, vegetation and plot corner markers with a Nikon DTM-500 total station during the snow-free season. Collecting the three-dimensional location of ∼2300 positions within the study site and several hundred surrounding it resulted in a point density greater than 2 points m⁻². We also surveyed all woody plants (sagebrush and coniferous trees) taller than 0.4 m on the site, as well as the height and crown radius and shape of all trees within 30 m of the study site perimeter (Fig. 2).

Fig. 2. Map of tree locations proximal to the study site (large square). Tree crowns are represented proportionally as circles. Solid black lines separate plots 1–5.

A high-resolution (0.5 m × 0.5 m) elevation model was generated from these points using universal kriging. Cross-validation calculated RMS (root-mean-square) prediction errors of 0.04–0.05 m (Reference LutzLutz, 2009). A tree-surface model was generated using inverse-distance weighted interpolations that approximated the size and shape of the observed trees (Reference LutzLutz, 2009). These surface models could then be integrated with the snow observations.

Weather station

Meteorological conditions during the surface hoar formation and persistence period were recorded continuously by a weather station located 14 m upslope of the study site and intermittently by field observers. A Campbell Scientific CR-10X data logger recorded snow depth, air temperature and wind speed and direction at ∼1.5 m above ground. Hourly snow-depth measurements using a Campbell Scientific Sonic Ranger 50M-45 allowed us to pinpoint when the surface hoar layer was buried by subsequent snowfall. A Campbell Scientific 107 Temperature Probe with radiation shield recorded hourly average air temperature, and a MetOne 034A wind sensor measured wind speed and direction. The logger saved basic statistics of wind speed and direction to memory every 10 min.

Sky visibility modeling

Sky visibility, v_% , was modeled using a modified viewshed analysis implemented in Arc Macro Language (AML). For each gridcell across the study site, the model calculated the percentage of the hemisphere obstructed by adjacent ground and trees. It builds conceptually on Reference Dozier and FrewDozier and Frew (1990) with the distinction of a higher-resolution hemisphere (5° altitudinal bands, 2° azimuthal bands) and a terrain model that incorporates proximate trees. Reference LutzLutz (2009) conducted independent optical measurements that produced spatial patterns of v_% similar to those produced by this GIS-based approach (r ² = 0.89; p < 0.0001).

Longwave radiation modeling

Surface hoar formation and growth depends on the development of a pronounced temperature gradient directly above the snow surface and adequate water vapor supplied through turbulent transfer, fostered by light winds (Reference NybergNyberg, 1938; Reference ColbeckColbeck, 1988; Reference Hachikubo and AkitayaHachikubo and Akitaya, 1997; Reference CoopersteinCooperstein, 2008). Determining the spatial patterns of the temperature gradient necessitates measuring the snow surface temperature, T _snow, across the site, which was not possible since field observations across the site would mean disturbing the site during the formation period. Remote-sensing technologies such as thermographic cameras and small deployable temperature and humidity sensors/loggers were not available at the time of this fieldwork.

Instead, we estimated the spatial pattern of incoming longwave radiation, L _↓, as a proxy of the night-time cooling potential at the snow surface during the surface hoar formation period. In the absence of shortwave (solar) radiation, the snow surface temperature, T _snow, is mainly controlled by the difference between incoming and outgoing longwave radiation fluxes, L _↓ and L _↑ respectively, referred to as the net longwave flux, L _net. Although the spatial pattern of L _↓ is nearly identical to that of v_% , modeling L _↓ enabled us to compare the scale of variation (in W m⁻²) with that of the modeled incoming shortwave radiation.

Based on the Stefan–Boltzmann equation, L _↓ can be described as a function of the percentage hemispheric viewshed occupied by sky (v_% ) and vegetation (1−v_% ) and their respective emissivity, ε, and temperature, T,

(4)

with the Stefan–Boltzmann constant of proportionality, σ, equal to 5.67 × 10⁻⁸ W m⁻² K⁻⁴. As proposed by Reference Bader and WeilenmannBader and Weilenmann (1992) and Reference HöllerHöller (2001) for snow-free canopies, T _tree was set equal to T _air, reducing Equation (4) to

(5)

This assumption seemed reasonable since under clear-sky conditions the tree surfaces cooled throughout the night. L _↓ was estimated at locations across the study site as defined in the previous equation, utilizing the GIS-derived estimates of v_% . T _air was set to −1.7°C to simulate conditions experienced in the early morning of 20 January, as recorded by the weather station. Following Reference Gubler and RychetnikGubler and Rychetnik (1991) and Reference HöllerHöller (2001), ε _tree = 0.94 and ε _air = 0.70. Because clear skies dominated at night-time we assumed that a relatively low emissivity could be applied. T _air was obtained from weather station measurements on 20 January and was treated as spatially constant across the study site.

Shortwave radiation modeling

After surface hoar is formed, its preservation depends largely on the snow surface remaining cold. Shortwave radiation heats the snow surface, which hinders crystal growth (by diminishing the air-temperature gradient) and can degrade surface hoar crystals (by metamorphism or melting). Hence, we modeled incoming global shortwave radiation at 15 min increments for the duration of a clear-sky day using the physically based Bird Clear Sky Model (Reference Bird and HulstromBird and Hulstrom, 1981 as implemented by D.R. Myers (http://rredc.nrel.gov/solar/models/clearsky/BIRD_02_05_2010.xls)). As with L _↓, we modeled shortwave radiation instead of solar indices because we wanted to assess the scale of spatial patterns in W m⁻².

We attempted to measure shortwave radiation with a LI-COR Li-200 pyronometer at the weather stations to verify model output. However, the measurements were not reliable for two reasons: (1) the pyronometer experienced tree shading intermittently throughout the day, including at noon, which made it impossible to compare values from the study site; (2) the sensor had been mounted parallel to the slope to record radiation values that correspond directly with the orientation of the study slope (personal communication from K. Reference KronholmKronholm, 2004). This, however, resulted in highly oblique solar measurements that were influenced by a cosine (i.e. Lambertian) response or effect (Reference IqbalIqbal, 1983).

We ran the Bird model for 20 January, a fairly warm, clear day during the surface hoar formation period for which we had concurrent weather observations. The model produced multiple estimates of solar irradiance, including diffuse and direct irradiance on horizontal surfaces, I _ah and I _dh respectively, as well as the direct irradiance on a surface oriented perpendicular to the solar rays, I _db.

The model utilizes several atmospheric parameters. Atmospheric extinction was parameterized using latitude, longitude, time of day and transmission information (atmospheric pressure, precipitable water, aerosol and ozone contents) extrapolated from surrounding weather stations in southwestern Montana and Idaho within 150 km of the study site (Reference LutzLutz, 2009). Detailed information about the sensor at these stations, data and data processing is available in NREL (2007). For the Bird model, the extraterrestrial solar constant was set to 1376 W m⁻² and the snow surface albedo was set to 0.9.

We integrated the Bird model results into an AML-based GIS model that spatially distributed the diffuse and direct shortwave radition across the study site incrementally over the modeled day.At 15 min intervals, a discrete hillshade analysis allocated an I _db value of 0 to cells in shaded areas. For non-shaded cells, the wattage of direct shortwave radiation as experienced on variable terrain, I _ds, was calculated by scaling I _db by the relative intensity, which is a function of the cosine of the angle of incidence, φ (the angle between an incoming solar ray, i, and the surface normal vector, s, at the point of incidence). To achieve this practically in grid-model space, I _ds was calculated using the three-dimensional vector components of i and s such that

(6)

where i_x , i_y and i_z are the x-, y- and z-components of the incoming solar direct radiation vector of strength, I _db, and t_x , t_y and t_z are the x-, y- and z-components of the slope normal vector, s. When calculating diffuse shortwave radiation, we followed Reference Wilson and GallantWilson and Gallant (2000) and accounted for terrain and trees partially obstructing sky visibility. The diffuse global shortwave radiation on the slope, I _as, was spatially distributed by multiplying I _ah by v_% :

(7)

At each time step, the sum of I _as and I _ds equaled the modeled global shortwave irradiance at each gridcell on the snow surface, I _s. We did not incorporate reflection from adjacent trees and slopes. Model output presented here includes the maximum global shortwave radiation I _max (W m⁻²) experienced at each cell location during the day, and the total global shortwave radiation experienced over the course of the day, ∑I (MJ m⁻²). Additional model outputs are described by Reference LutzLutz (2009).

Geostatistical analysis

Geostatistical analysis quantified spatial patterns of h _wl as trend surface and autocorrelation properties, which drove spatial interpolations. R statistical software packages (gstat and geoR) facilitated this analysis. Eleven linear and polynomial trend surfaces were tested to describe broad spatial features in the cross-slope, upslope and combined directions (Table 1). Reference KronholmKronholm (2004), Reference Logan, Birkeland, Kronholm and HansenLogan and others (2007) and Reference Lutz, Birkeland, Kronholm, Hansen and AspinallLutz and others (2007) similarly applied global spatial regression modeling in snow studies. In the absence of significant (p≤0.05) spatial trends, the plot mean was assumed the best-fit model. Otherwise, the model with highest significance (smallest p-value) was selected. Semi-variance analysis tested for positive spatial autocorrelation of residuals of the best-fit trend model. When spatial auto-correlation existed, a spherical semivariogram model was fitted using a weighted least-squares function (Reference CressieCressie,1993).

Table 1. Eleven tested spatial trend surface models. β _0−i are coefficients, x and y are sample coordinates, μ(z) is the predicted local mean and μ is the global mean

The trend surface and semivariance information drove spatial interpolations of h _wl. Depending on the trend surface and autocorrelation properties, local averages were interpolated using the plot mean, global regression, ordinary kriging or universal kriging, as noted below in the results. Reference LutzLutz (2009) also applied this information to random field interpolations, which can serve as drivers for propagation models (e.g. Reference Kronholm and BirkelandKronholm and Birkeland, 2005) and to illustrate how local variability obscures spatial patterns of local mean values.

Linear regression

Spatially implicit linear associations between radiation and h _wl and τ_∞ were quantified with weighted least-squares (WLS) regression. First, because both the shear frame and SMP measurements are destructive, they could not be paired spatially on the same unit of snow. Hence, it was not possible to determine whether a relationship existed between the variances of the two types of observations, based on individual pairs of observations. Given that the observations were clustered within tightly spaced local grids (Fig. 1), it was possible to estimate the local variance of both types of observations.

In addition, WLS is useful for small datasets and when the variance of a response variable is inconsistent across the range of explanatory values. The weights were inversely proportional to the variance of the predicted variable across the range of explanatory values. WLS also tested for associations between snow observations grouped by pit. Both snow observation techniques disturb the snowpack and are not spatially coincident, thereby warranting a grouped analysis of nested observations.

Surface hoar formation

The presence of hollow columnar surface hoar was observed on 20 January. Crystal height was typically 4–8 mm but in some instances was up to 15 mm. The crystals consisted of open-sided columns, ∼2–4 mm wide, with distinct scrolling on both columnar edges (Fig. 3). All slope aspects near the site were densely covered. Five days later (25 January 2005), tall plates (∼15 mm) had sprouted from the tops of the hollow columns, now partially decomposed (Fig. 3). This composite crystal structure was observed repeatedly before its burial on 5 February.

Fig. 3. (a) Enlarged illustration of initial hollow columnarcrystal form with scrolling edges on one side, observed on 20 January. (b) Photograph of surface hoar crystal, taken 3 weeks after crystal formation and 10 days after burial; stem and plate structures are evident.

Meteorological conditions

During the second half of January 2005, a diurnal cycle of relatively warm days with light cloud cover and cool clear nights developed. The surface hoar stem structure formed on 20 January. The meteorological conditions on 20 January and preceding days were mostly clear, with a night-time breeze averaging 0–2 m s⁻¹ and maximum gusts less than 10 m s⁻¹ (Fig. 4a). Winds were consistently out of the northeast (Fig. 4b), corresponding to an upslope movement, and were probably associated with a regional wind. Temperatures were mild but remained below 0°C (Fig. 4c).

Fig. 4. Meteorological conditions at the study site, 18 January–5 February 2005: (a) wind speed, (b) wind direction and (c) air temperature. Vertical dashed lines demark the first and second formation periods and the persistence period. Gap in coverage on 25 January is concurrent with field maintenance.

The sprouting plates formed sometime between 20 and 25 January. A diurnal pattern of T _air prevailed, with daily minima (−5°C) occurring in the early morning and maxima (1–3°C) occurring in the early afternoon (Fig. 4c). T _air varied by only 8°C during these 6 days. A diurnal pattern is also evident in wind speed: morning breezes grew steadily and reached peak values by late afternoon or early evening, with a maximum wind speed of 24 m s⁻¹ (Fig. 4a). Winds subsided to complete calm by the late evening or early morning. Winds were generally from the southwest (down-slope), although significant shifts occurred during daylight hours (Fig. 4b).

Between 25 January and 5 February, the surface hoar remained exposed on the surface. During this period, T _air reached 0°C twice, and on one occasion reached 3°C. Night-time lows varied between −5 and −10°C. Winds were slightly stronger than in previous days but maintained the diurnal pattern. Periods of sustained 5–10 m s⁻¹ wind occurred. Clear to partially cloudy conditions persisted throughout the period.

Radiation estimates

The site-scale spatial pattern of the GIS-distributed estimates of incoming longwave radiation, L _↓, was slightly concave, with high values at the slope base and along the right and upper edges of the site (Fig. 5a). L _↓ was greatest at the base of the slope where v_% was affected by downslope trees (v_% = 57%), and lowest near the middle of plot 4, where v_% was greatest (v_% = 75%). We do not include a map of v_% since it exhibits the same, but inverted, spatial pattern as L _↓.

Fig. 5. GIS-distributed radiation estimates: (a) incoming longwave flux L _↓ (W m⁻²) (T _a = −1.7°C), and Bird-model derived estimates of (b) maximum global shortwave radiation I _max (W m⁻²), (c) cumulative global shortwave radiation ∑I (MJ m⁻²) and (d) cumulative exposure time to direct shortwave radiation ∑t_I (hours). 1 m resolution. Black lines on scales show cumulative distribution of cell values.

We present three estimates of incoming solar radiation, including the maximum global shortwave radiation, I _max, the cumulative global shortwave radiation, ∑I, and the total exposure time to direct shortwave radiation, ∑t_I (Fig. 5). A cross-slope trend existed in I _max, with the right half of the site (plots 2 and 3) receiving about 85 W m⁻² more than the left half of the site (plots 4 and 5). Regarding cumulative solar energy, ∑I, plot 2 received the greatest amount of solar energy, and the outer corners of the remaining square plots received the smallest amounts. Lastly, a site-scale trend of solar exposure time ∑t_I ran diagonally downslope, with the lowest values occurring at plot 3 and the highest values occurring at plot 5.

Buried surface-hoar layer presence and thickness

Of 824 SMP profiles, 806 exhibited the surface hoar layer. Nine of the 18 profiles not exhibiting the weak layer were located directly where a rogue ski track had been cut through the outer corner of plot 3 in January (see ‘+’ in Fig. 6). Reference LutzLutz (2009) found no association between the locations of the remaining weak-layer absences and proximity to buried vegetation.

Fig. 6. Spatial interpolations of local mean h _wl (mm) using SMP data (a) from individual plots separately and (b) combined, based on best-fit trend surfaces and semivariance information (Table 2). ‘+’ indicates locations where surface hoar was absent in SMP profiles. Black curves on scales show cumulative distribution of cell values.

Trend surface analysis showed that weak-layer thickness, h _wl, increased in the upslope direction at plots 1, 2 and 5 independently (Table 2; Fig. 6a). This trend was also observed in surface hoar crystal size, observed by the shear frame operator. Individual h _wl values varied between 3 and 21 mm, with the minimum at plot 2 and maximum at plot 4. A subtle but significant cross-slope trend also existed. At common plot boundaries, the local interpolated mean values appear similar between adjacent plots (Fig. 6a). The semivariance analysis showed that at plots 1, 2 and 4, h _wl was locally correlated to distances of 2.6, 5.0 and 4.0 m respectively (Table 2). At plots 3 and 5, h _wl had no autocorrelation.

Table 2. Spatial patterns of h _wl. Estimates derived from SMP measurements (mm)

To determine whether h _wl observations from all five plots could be treated as a single site-scale dataset, we tested for significant differences of adjacent observations at plot 1 and the subsequent four square plots (Fig. 7), using the Wilcoxon rank-sum test (Reference WilcoxonWilcoxon, 1945). For each statistical test, only observations within 2 m of the common boundary between the compared plots were included. No thinning was evident between the adjacent subsets of these plots (Fig. 7). Weak-layer thinning may have been difficult to quantify because of the limited slab load which existed for most of the sampling period (Reference LutzLutz, 2009).

Fig. 7. Box plots compare h_wl at initial conditions (plot 1) with adjacent subsets of subsequent plots, namely (a) plot 2, (b) plot 3, (c) plot 4 and (d) plot 5, indicating that surface hoar layer thinning did not occur during the 3 week sampling period. Black square represents the median, whiskers represent the range, and box limits delineate the interquartile range. Triangles on boxes represent approximate 95% non-parametric confidence intervals of the mean (Reference McGill, Tukey and LarsenMcGill and others, 1978). The delta value is the difference in medians, accompanied by p-value when Wilcoxon rank sum was significant (p≤0.05). n is sample size.

The only significant (p = 0.006) difference was plot 2 possessing a thicker surface hoar layer than plot 1, which was probably a spatial phenomenon. Because no thinning was evident between the initial and subsequent observations along common plot boundaries, we considered it appropriate to also analyze all h _wl observations as a large site-scale dataset (Table 5; Fig. 6b). Both interpolations clearly show that h _wl was at a minimum at the base of the slope and was greatest in the upper left corner of the site in plot 4 (Fig. 6).

Snowpack associations with radiation estimates

Significant (p≤0.05) inverse linear relationships existed between h _wl and all radiation estimates at the site scale, i.e. when all observations were included and when only plot 1 observations were included (Table 3). Positive correlations between h _wl and v_% mirror those with L _↓ (recall that the spatial pattern of L _↓ is largely determined by v_% ). Significant correlations also existed at the smaller spatial extents for L _↓ at plots 2 and 5, for I _max at plot 2 and for ∑t_I at p lot 4. All significant fits indicate that h _wl decreased with increasing exposure to longwave and shortwave radiation.

Table 3. WLS regression fits between h _wl and sky visibility and radiation estimatesFootnote *

* WLS weights set to the inverse of h _wl pit variance. Significant (p ≤0.05) fits shown in bold.

Unlike surface hoar thickness h _wl, shear strength τ _∞ changed significantly over time (Reference LutzLutz, 2009). Thus, site-scale correlations between τ _∞ and radiation estimates could not utilize observations from all plots pooled together because that would incorporate variance associated with shear strengthening observed over the sampling period. At plot 1, which spanned the entire site, significant positive linear correlations were identified between pit estimates of τ _∞ and the shortwave radiation variables ∑I and ∑t_I (Table 4), indicating that higher shear strength was indirectly associated with greater exposure to solar radiation at the site scale. No significant correlations existed between τ _∞ and I _max or L _↓ at the plot scale.

Table 4. WLS regression fits between τ _∞ and radiation estimatesFootnote *

* WLS weights set to the inverse of h _wl pit variance. Significant (p ≤0.05) fits shown in bold.

Direct associations between surface hoar thickness and shear strength

Significant inverse linear relationships existed between h _wl and τ _∞ on the upper slope, at plots 3 and 4 in the WLS analysis (Table 5). At these plots, the regression fits indicate that shear strength increased across space at a rate of 52 and 88 Pa mm⁻¹, respectively. All other plots had no significant relationship between observations grouped by pit, although the signs of the fitted lines were generally also negative.

Table 5. WLS regression fits between τ _∞ and h _wl Footnote *

* WLS weights set to the inverse of h _wl pit variance; Significant (p ≤0.05) fits shown in bold. p-value is significance level, r ² is explanatory strength, b is y-intercept and m is coefficient.

Meteorological conditions during surface hoar formation and persistence

The composite stem-plate structure of the surface hoar crystals resulted from two consecutive crystal growth regimes (Reference Jamieson and SchweizerJamieson and Schweizer, 2000). General meteorological conditions for both formation periods can be drawn from the weather station data. The surface hoar stem crystals formed during light wind conditions, which are known to foster hoar formation through turbulent transfer of water vapor (Reference ColbeckColbeck, 1988; Reference Hachikubo and AkitayaHachikubo and Akitaya, 1997; Reference CoopersteinCooperstein, 2008). The absence of moderate or strong winds also allows the necessary temperature gradient at the snow surface to develop (Reference NybergNyberg, 1938; Reference ColbeckColbeck, 1988; Reference Hachikubo and AkitayaHachikubo and Akitaya, 1997).

The subsequent plate-like crystals formed during a period dominated by diurnal patterns associated with solar-induced warming. Winds varied considerably during this time, but were generally strongest in the afternoon or early evening, subsiding to complete calm in the late evening or early morning. Irregularity in wind speed, indicated by the range in 10 min minimum and maximum recorded wind speeds (Fig. 4a), may have also fostered crystal growth by periodically resupplying vapor to the boundary layer without causing a temperature gradient to be diminished, as would occur under sustained winds (Reference ColbeckColbeck, 1988).

During both surface hoar formation periods, T _air at 1.5 m above the snow surface never dropped below −5°C and often approached 0°C (Fig. 4c). Reference Chalmers and JamiesonChalmers and Jamieson (2001) observed similar air temperatures during surface hoar formation. The development of pronounced temperature gradients at the snow surface (Reference NybergNyberg, 1938; Reference KobayashiKobayashi, 1961; Reference Lang, Leo and BrownLang and others, 1984) could foster surface hoar development under these conditions.

Our work shows that on a northeast-facing slope at this latitude in the Northern Hemisphere, surface hoar can remain intact for many days despite high winds, relatively high air temperatures and direct shortwave radiation. While weather station data helped establish general meteorological conditions, the spatial patterns of wind and temperature were not observed or simulated. Such information would be valuable in future field studies examining the spatial distribution of surface hoar growth and persistence. Potential spatial patterns in winds may have resulted in trends of sensible heat flux and vapor transport across the study site, both of which would influence growth and persistence of surface hoar.

Associations between incoming radiation and snowpack properties

Within the 31 m × 31 m study site, h _wl varied by a factor of 7 (3–22 mm) and was inversely correlated with spatial patterns of incoming radiation estimates. Shea and Jameson (2010) have also observed extreme variations in surface hoar thickness across sparsely forested terrain.

The shear strength, τ _∞ , was positively correlated with exposure to shortwave radiation, ∑I and ∑t_I , across plot 1. Incoming longwave radiation inhibits surface cooling at night, a precondition for surface hoar formation, and shortwave radiation causes surface warming by day, which fosters crystal degradation. The inverse relationships between L _↓ and h _wl at plots 1, 2 and 5 (Table 3) support previous associations identified by Reference Gubler and RychetnikGubler and Rychetnik (1991) and Reference HöllerHöller (2001). Negative correlations between h _wl and I _max, ∑I and ∑I_t were also identified (Table 3). Our estimates of ∑I coincide with the minimum values in the surface hoar preservation regime observed by Reference Feick, Kronholm and SchweizerFeick and others (2007).

These findings support Gubler and Rychetnik’s (1991) assertion that surface hoar formation will be most productive in places where both shading and sky visibility spatially coincide. In the Northern Hemisphere these microclimatic factors are easily coupled on forest openings on northerly aspects, where sky visibility and daytime shading coincide near the center of the openings and together encourage low surface temperatures throughout the night and day (Fig. 8). In addition, trees located at the base of such openings can receive intense solar radiation on their upslope boughs during the day, increasing re-radiation on the lower portion of the site and further enhancing the existing pattern of L _↓. Such factors could be well modeled using distributed radiation modeling (e.g. Reference Adams, Gauer, McKittrick, Curran, Levenan and GanjuAdams and others, 2004).

Fig.8. Conceptual model of (a) night-time and (b) daytime patterns in surface temperature on a northeast-facing (right) slope in the Northern Hemisphere. Bold dashed lines represent approximate upslope extent of study site.

In this dataset, spatial colinearity of longwave and shortwave radiation estimates made it difficult to discern whether one was more responsible for the spatial patterns of h _wl. Based on the following three characteristics, it is likely that both radiation types influenced the spatial patterns of h _wl: (1) correlations with both radiation types were significant, (2) all significant correlations were logical in the sign of the coefficients and (3) when the albedo of global shortwave radiation is accounted for (by multiplying I _max values by 0.1, assuming an albedo of 0.9), the magnitude of the spatial variations of absorbed radiation is similar for long- and shortwave estimates.

Associations between surface hoar thickness and shear strength

Direct correlations between the spatial patterns of h _wl and τ _∞ were significant only at plots 3 and 4. This was surprising given that both snow properties were spatially correlated with incoming radiation patterns. The lack of statistical significance may have been the result of the small sample sizes (six shear frames per pit, nine pits per plot).

Although here we examine the spatial patterns of the surface hoar properties, it is worth noting that no weak-layer thinning coincided with the observed shear strengthening over the 3 week-long sample period (Reference LutzLutz, 2009). Reference Chalmers and JamiesonChalmers and Jamieson (2001) also documented a young buried surface hoar layer that strengthened over 4 weeks despite no evidence of thinning during that period. Surface hoar strengthening is typically associated with thinning caused by inclination or penetration into the substratum (Reference Davis, Jamieson and JohnstonDavis and others, 1998; Reference Jamieson and JohnstonJamieson and Johnston, 1999; Reference Jamieson and SchweizerJamieson and Schweizer, 2000; Reference Lutz, Birkeland, Kronholm, Hansen and AspinallLutz and others, 2007). Our field observations indicate that the strengthening occurred at the transition with the superstratum. At initial conditions the surface hoar remained attached to the substratum after each shear test. As the sampling period progressed, the surface hoar ripped away from the substratum and remained attached to the superstratum (in the shear frame). Presumably the substratum interface would have strengthened more under a heavier slab load (Reference LutzLutz, 2009).

Implications for stability evaluation

The observed spatial patterns of h _wl and τ _∞ highlight the conundrum that avalanche forecasters, ski patrollers and backcountry enthusiasts face when trying to choose and sample representative slopes or portions of slopes. Assessing stability at the edge of a forest opening, where the surface hoar is less developed than in the middle, can result in an overestimation of stability, or a false stable assessment. Yet sampling the middle of the slope, if steep enough, could be hazardous. On the other hand, if the observer finds any indication of unstable conditions at the edge of an opening, it is very likely that the buried surface hoar layer could be thicker and weaker in the middle to upper portion of the opening.

For researchers, these findings have ramifications for field-based snow studies. During the planning stage of field studies, scientists need to consider the potential for slope-scale trends to be present, even on relatively uniform terrain. This is especially critical for temporal studies attempting to monitor changes in snowpack characteristics over time. Using our, or similar, sampling strategies, temporal studies can quantify and account for spatial patterns at initial conditions only if complex spatial patterns can be discerned and accounted for. This methodology could be applied to extensive high-resolution terrain, vegetation and snowpack surface models (e.g. Deems and Painter, 2006) to develop sampling strategies that eliminate potential colinearity of environmental determinants. This could be done a priori by modeling spatial patterns of radiation and other environmental factors before conducting snow observations. The spatial patterns of incoming radiation, although generated statistically, illustrate the complexity that vegetation and terrain introduce to the local conditions and the importance of incorporating spatial heterogeneity into process-driven snowpack and boundary layer modeling of surface hoar and near-surface faceted crystal formation (Reference Lehning, Bartelt, Brown and FierzLehning and others, 2002; Reference Adams, Gauer, McKittrick, Curran, Levenan and GanjuAdams and others, 2004; Reference Stössel, Guala, Fierz, Manes and LehningStössel and others, 2010).