Climatology of the Canadian Cordillera

The Canadian Cordillera comprises the Rocky, Columbia and Coast Mountain ranges extending from the US border into northern British Columbia. These mountains have a glaciated area of approximately 38 613 km² (Reference OmmanneyOmmanney, 2002) and experience a long snow-accumulation season, generally lasting from October to June. The climate of this region is influenced by maritime polar (mP) air masses that originate in the north Pacific Ocean, maritime tropical (mT) air masses that originate in the subtropical Pacific, and continental polar (cP) air masses that form over northern Canada (Fig. 1a).

Fig. 1.

(a) Opabin Plateau and Opabin Glacier showing snow-pit sites OG (Opabin Glacier) and OP (Opabin Plateau). The inset shows the position of the study area in North America along with the source regions of maritime polar (mP), maritime tropical (mT) and continental polar (cP) air-masses. (b) Haig Glacier showing snow-pit sites HG (Haig Glacier) and HF (Haig Forefield).

The precipitation regime of southwestern Canada is dominated by mid-latitude cyclones that cross southern British Columbia in the prevailing westerly flow. These weather systems bring mP and mT air and meet a series of steep north–south-oriented mountain ranges on their path from the Pacific coast to the Rocky Mountains, with altitudes ranging from ∼1500 to ∼3500 m a.s.l. These barriers generally increase in elevation towards the continental divide and induce orographic uplift, cooling and precipitation formation.

The most commonly used wintertime synoptic classification for the Rocky Mountains was developed by Reference Changnon, McKee and DoeskenChangnon and others (1993), who describe seven 500 hPa synoptic patterns associated with precipitation over the northern US Rocky Mountains. This classification scheme was adopted by Reference Byrne, Berg and TownshendByrne and others (1999) and Reference Lapp, Byrne, Kienzle and TownshendLapp and others (2002) to forecast river runoff and winter precipitation for western North America, respectively, and was shown by Reference Moore and McKendryMoore and McKendry (1996) to be consistent with patterns identified using Kirchhofer climate typing for spring snowpack anomalies in British Columbia.

Reference Moran, Marshall, Evans and SinclairMoran and others (2007) identified three circulation types for winter storms bringing snow to Robertson Glacier (adjacent to Haig Glacier): northwesterly, westerly and upslope. Northwesterly and westerly circulation is often associated with intense accumulation events in this region, particularly when it interacts with outbreaks of cold air from the Arctic (Reference StewartStewart and others, 1995). Surface high-pressure systems in northern Canada push cold cP air southwards on the eastern side of the continental divide. This alone does not always produce a significant amount of snowfall in the Rocky Mountains, but, when it is combined with an impulse of moist Pacific air from the west, the region can receive heavy snowfall over several days. Upslope storms refer to snow that is derived from the easterly sector, i.e. storms travel in opposition to the prevailing westerly flow in this region. These weather systems originate from either lee-slope cyclogenesis or the passage of a mid-latitude cyclone to the south of Alberta (Reference Chung, Hage and ReineltChung and others, 1976). They typically occur from March until May and are associated with heavy snowfall that accompanies the southeasterly advection of mP or mT air to the eastern slopes of the Rocky Mountains.

Field sites and analysis

The two field sites selected for this study are in catchments that are the focus of long-term studies of glacier mass balance, surface energy-balance processes and alpine hydrology. Opabin Glacier is a small (∼0.4 km²) niche glacier located at 2300–2450 m a.s.l. in Yoho National Park, British Columbia, on the west-facing (windward) slopes of the Rocky Mountains (Fig. 1a). Haig Glacier is in an eastfacing (lee-slope) setting in Peter Lougheed Provincial Park, Alberta (Fig. 1b). It is the largest arm of the southernmost icefield in this region and has an area of approximately 2.9 km². The upper glacier (2800 m a.s.l.) is a broad gently sloping plateau that flows southeast from the North American continental divide to approximately 2430 m a.s.l.

These field sites were visited during the 2004/05, 2005/06 and 2006/07 snow-accumulation seasons and into the subsequent summer melt seasons. In 2004/05, work was focused on Haig Glacier, and a late-season snow pit was sampled at Opabin Glacier. In 2005/06, snow pits were sampled intensively at both field sites through the accumulation season, and in 2006/07 mid- and late-winter pits were sampled at both sites.

At Haig Glacier, snow pits were dug in the vicinity of two automatic weather stations (AWS). One AWS is located mid-glacier near the equilibrium line (HG) and the other is approximately 100 m from the glacier terminus on the forefield moraine (HF). Similarly, at Opabin Glacier, snow pits were dug near an ultrasonic depth gauge (UDG) at a site near the centre of the glacier (OG) and on the deglaciated Opabin Plateau near an AWS (OP). The lower HF and OP sites are referred to as forefield sites, while the upper HG and OG sites are referred to as glacier sites. On each site visit, snow pits were dug to expose the cumulative seasonal precipitation. Each snow-pit wall was sampled for density, stable-isotope analysis and the temperature profile. At the glacier sites snow pits were dug back to the summer-melt surface, and at the forefield sites snow pits were dug to the moraine surface. Samples were taken upward from the glacier or moraine surface, using a 100 cm³ steel sampler at 10 cm vertical resolution in 2004/05 and 5 cm vertical resolution in 2005/06 and 2006/07.

Table 1 summarizes the dates, names, total snow water equivalent (SWE) and number of samples from all pre-melt snow pits. The mean April snowpack was similar on both glaciers over the study period, averaging 1146 mm w.e. On average, April SWE was greater by a factor of 1.9 on the glaciers, relative to the forefield settings.

Table 1.

Snow-pit dates, names, total SWE and number of samples (n) for all pre-melt snow pits. Snow pits are numbered consecutively followed by the abbreviated year (e.g. HF105 is the first snow pit of 2004/05 at the Haig Forefield)

Snow samples were bagged and transported to Calgary frozen, where they were melted at room temperature and bottled for stable-isotope analysis. Analysis was conducted at the University of Calgary Stable Isotopes Laboratory using the CO₂–H₂O equilibration technique for δ¹⁸O and chromium reduction for δD. In both cases a dual-inlet isotope-ratio mass spectrometer was used to determine the δ values, which are reported in ‰ relative to the standard, Vienna Standard Mean Ocean Water (V-SMOW). Accuracies are better than ±0.2‰ for δ¹⁸O and ±2.0‰ for δD.

Locating events in isotope stratigraphies

The snow-surface height at each snow-pit site was recorded at half-hour intervals using an SR-50 ultrasonic depth gauge suspended on a horizontal bar above the snow surface. A surface-height record was produced at each site by correcting for air-temperature effects and removing spurious data that occur when blowing snow passes the sensor. We only recovered complete snow-surface height data from the forefield sites due to data loss from the HG AWS and excessive drift in the record from the OG site. Analysis of snowfall events is therefore based on the forefield snow-height records at both field sites.

Accumulation events, defined from the time of minimum to maximum snow height, were identified in the records from both sites. These are characterized by rapid increases in snow height followed by a slower, usually exponential, decrease in snow height during settling. They were numbered as either single events or as multi-stage events if there was minimal settling between periods of accumulation. Only those events that resulted in >5 cm total accumulation were included in the analysis to separate snowdrift from surface-height changes generated by synoptically driven accumulation.

While the SR-50 records present a seasonal chronology of accumulation events, isotope stratigraphies are derived from subsurface snow that is subject to progressive densification over the accumulation season. This densification is spatially variable, both between sites and at different depths within the snowpack, and occurs at varying rates throughout the course of the winter. For this reason, the datasets need to be processed to identify unique accumulation events in the snow-pit isotope stratigraphies.

Snow-pit density data were used to calculate the total SWE for each snow column:

where t _l is the thickness of the sampled layer (mm), ρ ₁ is the density of each layer (kg m⁻³) and ρ _w is the density of water (1000 kg m⁻³). An estimate of the mean fresh-snow density from each snow pit, , was obtained from the accumulated SWE_f and the total accumulation from the SR-50 record, d _f , between snow-pit sampling dates,

Mean fresh-snow densities derived using this method ranged from 97.3 to 131.7 kg m⁻³, which is consistent with field measurements of fresh-snow density at these sites. The mean fresh-snow density for each snow pit was then used to estimate the SWE of each accumulation event at the OP and HF AWS sites. The same fresh-snow densities were applied to the HG and OG sites, but the ratio of snow-pit SWE between the forefield and glacier sites was used to scale the accumulation events to account for the precipitation gradient between the lower and upper sampling locations.

Once snow-accumulation records were converted to SWE, accumulation events were fitted to the snowpack δ¹⁸O stratigraphies. For the first snow pit of each season, the stratigraphies of the forefield and glacier sites were compared to identify features from the higher-elevation sites that were established prior to snow accumulation at the forefield (AWS) sites. This early season snow was removed from the analysis because there are no corresponding SR-50 records from the glacier sites from which we could identify early-season accumulation events. In subsequent snow pits, the structure of the δ¹⁸O record was used to find the upper edge of the previously sampled snow, and accumulation events were distributed from this point upward within the snowpack. This was repeated for all snow pits until each storm was located in an isotope stratigraphy.

In some cases, accumulation event boundaries were adjusted so that they encompassed the same features of the δ¹⁸O stratigraphy at the forefield and glacier sites. The storm boundaries were also cross-checked between the two field sites so that they occurred at consistent locations within the δ¹⁸O stratigraphies. An assumed relationship between stable isotopes and temperature was also used to guide the final position of storm boundaries in the snowpack. Minor adjustments were made so that very cold events, for example, were associated with isotopically depleted layers. While it is expected that stable water isotopes will have a stronger relationship with condensation temperature than near-surface temperatures, sampling sites lie at or above the cloud base during winter precipitation, making surface AWS data a good surrogate for condensation temperature in this case.

Figure 2 illustrates the process of fitting accumulation events to snow-pit stratigraphies. The OG406 and OP406 δ¹⁸O stratigraphies are shown along with the snow-accumulation and temperature records from the OP AWS from 9 January (when previous snow pits were sampled) up to the 30 January sampling date. Four major accumulation events occurred over this time period and were recorded in both the OP and OG snow-pit stratigraphies. These events were distributed upwards from the top of the 9 January OP and OG δ¹⁸O stratigraphies. The SWE of each event in the OP snowpack was based on the estimated mean fresh-snow density, and these events were scaled to fit the OG snowpack.

Fig. 2.

Snow-accumulation and temperature records from the OP AWS from 9 to 31 January 2006 showing accumulation events 14–17 and the mean temperature for each event (in grey bars). Also shown are snow pits OG406 and OP406, which were both sampled on 30 January 2006. The solid grey lines mark the final accumulation event boundaries placed using an estimated fresh-snow density.

Trajectory classification procedure

In order to investigate the effect of vapour pathway on isotopic stratigraphies, a trajectory classification was produced to group each accumulation event into a trajectory class. There are two main approaches to developing a synoptic classification scheme. The first is subjective and involves manually classifying the atmospheric continuum into categories. While manual classifications are conceptually straightforward, this process is labour-intensive, difficult to replicate and best suited to smaller datasets (Reference YarnalYarnal, 1993; Reference Frakes and YarnalFrakes and Yarnal, 1997). The second approach is multivariate classification analysis, such as cluster analysis, which avoids preconceived groupings but requires a larger dataset, usually over a longer time period (e.g. Reference Vuille, Hardy, Braun, Keimig and BradleyVuille and others, 1998). A manual approach is used in this study, due to the relatively small number of storms (n = 73) along with the regional nature and short duration of the study. The suite of tools used to produce the classification includes the following:

1. 1. 500 hPa geopotential height maps: 500 hPa is generally considered representative of mid-tropospheric circulation and is a useful source of information about atmospheric flow paths (Reference Changnon, McKee and DoeskenChangnon and others, 1993). Daily 500 hPa maps are available from both Environment Canada and the US National Centers for Environmental Prediction (NCEP)/US National Center for Atmospheric Research (NCAR) reanalysis dataset over a 2.5° × 2.5° global grid (Reference KalnayKalnay and others, 1996; http://www.cdc.noaa.gov). NCEP/NCAR geopotential height data were extracted for an area covering western Canada (40–65° N, 145–105° W), and daily 500 hPa maps were produced to cover the duration of each storm, from 1 day before until 1 day after each event (based on GMT start and end times).

2. 2. Back trajectories: Daily 72 hour back trajectories from the US National Oceanic and Atmospheric Administration’s HYSPLIT (HYbrid Single-Particle Lagrangian Integrated Trajectory) model (www.arl.noaa.gov/ready/hysplit4.html) were produced using NCEP/NCAR reanalysis data. These are a useful tool to help assess the vapour source and pathway of an advected air mass (Reference Davidson, Jaffrezo, Small, Summers, Olson and BorysDavidson and others, 1993; Reference Kahl, Martinez, Kuhns, Davidson, Jaffrezo and HarrisKahl and others, 1997), and, in this case, each trajectory describes the most probable upwind path taken by a parcel of air arriving at the study sites at 2500 m a.s.l. This height corresponds to around 750 hPa, and was chosen because it is the approximate altitude of vapour advection. The 72 hour time frame was selected to enable trajectories to extend into the Pacific Ocean (for storms from the westerly sector), so that the entire trajectory over land could be visualized.

A sequencing technique using daily geopotential height maps was used to classify accumulation events (Reference Yarnal and FrakesYarnal and Frakes, 1997). This reduces the model error in manual synoptic classifications by grouping storm days so that the source region of a mid-latitude storm system can be defined and the trajectory tracked over the course of the event. This process was aided by compiling an ensemble of daily back trajectories into a geographical information system for each accumulation event. The trajectory lines were assigned different symbols, depending on the total daily accumulation, so that high-accumulation days could be easily identified in the ensemble. Finally, composite geopotential height maps were used to define the average flow conditions over the course of each accumulation event. The use of composites is well established in synoptic climatology and can help to resolve some of the within-storm variability that emerges from sequencing weather maps (Reference Yarnal and FrakesYarnal and Frakes, 1997). For this study, the composites were produced for each storm from NCEP 500 hPa height data, which are available four times daily (00Z, 06Z, 12Z and 18Z). The time frame of each composite was based on the closest 6 hour interval to the storm start and end times at both study sites.

The meteorological characteristics of each accumulation event were investigated using data from each field site and from an array of weather stations in the region (Fig. 3). For all Environment Canada stations, data were downloaded at http://climate.weatheroffice.ec.gc.ca/climateData/hourly-data_e.htm. Where available, the average relative humidity, wind direction, wind speed, pressure and temperature along with the maximum change in temperature, dT _max, were calculated for each storm. At Environment Canada stations, precipitation data were derived from daily totals and the occurrence of precipitation was used to guide the storm onset and end times. These times were refined by examination of the hourly weather observations for each storm in the case of Calgary Airport (YYC) and Golden. For Banff and Pocaterra, where there are no recorded manual weather observations, storm durations were derived from daily precipitation data only.

Fig. 3.

Map showing the locations of the glacier study sites and weather stations used in the storm trajectory classification. All sites are Environment Canada weather stations, aside from Haig and Opabin Glaciers and Kananaskis Field Station (K.F.S.).

Upper-air conditions, including zonal and meridional flow, total flow strength and flow direction, were determined for each accumulation event from daily NCEP 500 hPa data for a point over Haig Glacier (50° N, 115° W) (Reference Shea and MarshallShea and Marshall, 2007). These flow indices were also used to guide the classification process described above. Mixing-ratio data from upper-air soundings at Kelowna International Airport (YLW) were downloaded from the University of Wyoming website at http://weather.uwyo.edu/upperair. Kelowna is upwind of the field sites under westerly flow, but is more than 400 km from the field sites so it is not strictly representative of lower-tropospheric conditions in the Canadian Rockies. Nevertheless, the percentage change in mixing ratio from the start to the end of each storm does give an indication of the degree of drying associated with each event.

Trajectory classification

In order to effectively relate storm trajectories to isotopic values in winter snowpacks, it was considered advantageous to produce a broad classification scheme. This included the westerly (W), northwesterly (NW) and upslope (UPS) classes defined by Reference Moran, Marshall, Evans and SinclairMoran and others (2007), along with one additional category, southwesterly (SW). We were unable to classify 6 of 73 storms, which are grouped together as unclassified (UNC).

Figure 4 shows the percentage of total accumulation represented by each class. Accumulation amounts from each event at the OP and HF AWS sites were averaged to derive totals for each storm in 2005/06 and 2006/07, while accumulation totals for 2004/05 are derived from the HF AWS (as snow pits were not sampled at Opabin Glacier over that season). Evident in Figure 4 is the dominance of southwesterly storms in 2004/05. These brought 46.5% of the total accumulation, while upslope storms were responsible for 24.5% of the seasonal moisture supply. In the following two seasons, westerly and southwesterly events brought more than 75% of the seasonal precipitation. Westerly storms were relatively more important to snow accumulation in 2006/07 (41.4%) than in 2005/06 (35.6%). For southwesterly storms, this was reversed, with these trajectories being relatively more important to accumulation in 2005/06 (48.6%) compared with 2006/07 (34.4%).

Fig. 4.

Percentage of total seasonal accumulation represented by each storm trajectory class.

Meteorological characteristics of trajectory classes

Figure 5 shows the 500 hPa geopotential height composite maps and an ensemble of daily back trajectories for a sample accumulation event from each trajectory class. The key meteorological variables associated with each trajectory class for the alpine field sites, Calgary and Banff are presented in Table 2. These values are used in conjunction with the synoptic classification tools described above to characterize each trajectory class.

Fig. 5.

500 hPa geopotential height composite maps (left) and daily back-trajectory ensembles (right) for one accumulation event from each trajectory class. Black boxes encompass alpine field sites, and trajectories are calculated backward from Haig Glacier.

Table 2.

Mean accumulation and surface temperature, and maximum change in temperature (dt _max, for each trajectory class from the Haig Forefield, the Opabin Plateau, Calgary and Banff. Also given are the mean 500 hPa wind direction and wind speed from flow indices, and the mean mixing ratio and change in mixing ratio (dq) from Kelowna radiosonde data

The isotopic characteristics of Rocky Mountain snowpacks

The δ¹⁸O and δD values of precipitation are controlled by a combination of continental, altitude and temperature effects as identified by Reference DansgaardDansgaard (1964). Continentality, or distance inland, influences the degree of distillation of heavy isotopes from an air mass during rain-out as the vapour travels away from its source region (Reference Friedman, Smith, Gleason, Warden and HarrisFriedman and others, 1992). The altitude effect is intrinsically linked to the continental effect because, as air masses travel inland and meet steep topography, the induced orographic rain-out causes depletion in heavy isotopes with elevation (Reference Rozanski and SonntagRozanski and Sonntag, 1982). This was shown by Reference Moran, Marshall, Evans and SinclairMoran and others (2007) to be complicated in the Rocky Mountains by trajectory. Eastward-moving weather systems produce a conventional depletion in stable isotopes on windward slopes, while the continued rain-out on the lee slopes of the continental divide produces ongoing depletion of heavy isotopes as elevation decreases on eastern slopes. Reference Peng, Mayer, Harris and KrousePeng and others (2004) hypothesize that this is the reason why Calgary (Fig. 3) has lower annual average δ¹⁸O and δD values than expected for its latitude.

Figure 6a shows the density-weighted mean δ¹⁸O values for all end-of-season snow pits. The Opabin Glacier sites, OG/OP, had lower mean δ values than the Haig sites, HG/ HF, in all three accumulation seasons. The average δ¹⁸O for all seasons at Haig Glacier was −19.5‰ compared with −21.0‰ at Opabin Glacier. The finding by Reference Moran, Marshall, Evans and SinclairMoran and others (2007) that continued distillation of heavy isotopes beyond the continental divide leads to more negative δ values downwind suggests that Haig Glacier should receive snow that is more depleted in heavy isotopes than Opabin Glacier due to its easterly aspect. The latitude of these field sites along the continental divide may, however, provide an explanation for the discrepancy. Opabin Glacier lies downwind of Glacier National Park, British Columbia, which contains steep north–south-oriented topography and numerous peaks over 3000 m a.s.l. It is likely that this topography induces a significant amount of rain-out of heavy isotopic species before precipitation-bearing weather systems reach the continental divide.

Fig. 6.

(a) Mean δ¹⁸O at each field site in each accumulation season. (b) Mean deuterium excess (d) at each field site in each accumulation season.

Figure 6a also reveals that snow at lower-elevation forefield locations (HF and OP) was consistently more depleted in heavy isotopes than snow at the glacier sites. The average δ¹⁸O at the HF site for all samples was −20.9‰ compared with −19.9‰ at the HG site. This trend was also apparent at Opabin Glacier, where the average δ¹⁸O at the OP site was −22.5‰ compared with −21.0‰ at the OG site. This inversion of the usual altitude effect may be due to the bias introduced by early-season (isotopically heavier) snow that is retained in the glacier snowpacks and does not accumulate at the forefield sites. For example, the mean late-season δ¹⁸O at the HF site in 2004/05 (HF305) was −21.6‰, while the mean late-season δ¹⁸O of HG305 was −19.5‰ for the entire snow pit and −20.5‰ with early-season snow removed. The same pattern was observed in the subsequent two seasons at both field sites: removal of early-season snow from the glacier sites brings mean glacier δ¹⁸O closer to the forefield values, but does not fully account for the discrepancy between sites. At Haig Glacier, there could also be ongoing distillation of heavy isotopes as water vapour moves across the continental divide from the head of the glacier to the terminus.

The global relationship between δD and δ¹⁸O, known as the global meteoric waterline (GMWL), was defined by Reference CraigCraig (1961) as δD = 8δ¹⁸O + 10. The slope and intercept of the line of best fit for all δD–δ¹⁸O data from glacier sites are 7.9 and 8.8 respectively. Regression lines were also produced for HG and OG data, and for each of these sites in each accumulation season. The slope of the regression lines was comparable across all seasons, ranging from 7.6 to 8.1. The intercepts were similar in most cases, ranging from 7.7 to 12.7, with the notable exception of the HG data in 2004/05, which had an intercept of −0.9.

The deuterium excess (d) of alpine snowpacks in the Canadian Cordillera has the potential to be a useful diagnostic tool in identifying vapour source regions for individual accumulation events. Defined as d = δD – 8δ¹⁸O, this parameter is thought to be controlled primarily by the relative humidity of the atmosphere above the oceanic source region and secondarily by both the air temperature and wind speed during evaporation (Reference Merlivat and JouzelMerlivat and Jouzel, 1979). Source regions with lower relative humidities generally produce higher d values in vapour due to enhanced kinetic fractionation during evaporation, and this is augmented by in-cloud kinetic fractionation during the formation of snow (Reference Jouzel and MerlivatJouzel and Merlivat, 1984). Deuterium excess is also modified when additional moisture is added during vapour transport (Reference Ciais, White, Jouzel and PetitCiais and others, 1995; Reference Araguás-Araguás, Froelich and RozanskiAraguás-Araguás and others, 1998). In some regions, the effect of relative humidity on d results in a distinct seasonality due to source region shifts (Reference Schlosser, Oerter, Masson-Delmotte and ReijmerSchlosser and others, 2008). Reference Peng, Mayer, Harris and KrousePeng and others (2004) found no discernible seasonality in 10 years of deuterium excess data from Calgary, although we see evidence of lower d in spring, with a mean d from both field sites of 8.2‰ in March–April vs 11.1‰ for October–February.

The mean d for all snow-pit data in this study is 10.0‰. While this is identical to the global average determined by Reference CraigCraig (1961), d values differed significantly at adjacent glacier and forefield sites (Fig. 6b). As with mean δ values, these differences can be partly attributed to the accumulation of early-season snow at the glacier sites and not at the forefield sites. We have also observed post-depositional modification of isotope stratigraphies at the forefield sites not seen at the glacier sites (Reference Sinclair and MarshallSinclair and Marshall, 2008). This may be due to melting conditions at the base of the forefield snow pits, which can lead to the presence of liquid water (hence evaporative losses) and mass exchange with the underlying ground. The isotopic effects of post-depositional modification are compounded in a derived variable like deuterium excess, so it is probably unwise to use forefield d values for vapour trajectory analysis.

The glacier snowpacks show no sign of post-depositional modification prior to the onset of spring melt, so we believe that these records can be reliably interpreted (Reference Sinclair and MarshallSinclair and Marshall, 2008). The mean end-of-season d was 11.1‰ and 11.4‰ in 2005/06 and 10.2‰ and 10.3‰ in 2006/07 at the HG and OG sites, respectively (Fig. 6b). The broadscale similarity in d between the glacier sites suggests that Opabin and Haig Glaciers received moisture from the same source regions over these accumulation seasons. The mean end-of-season d at Haig Glacier in 2004/05 was significantly lower (8.4‰), and, as shown in Figure 4, the dominant storm trajectories in 2004/05 differed from the subsequent two seasons. While snow derived from southwesterly trajectories was important in all seasons, it was comparatively more dominant over westerly trajectories in 2004/05. In addition, many of the southwesterly events in 2004/05 were associated with Aleutian coastal low-pressure systems. These draw moisture from the north Pacific, rather than the subtropics, and are expected to originate from higher relative-humidity source regions. This is supported by the fact that the mean d from southwesterly events in 2004/05 was 7.1‰.

Temperature and trajectory controls

The mean storm δ¹⁸O values derived using SR-50 records and snow-pit isotope stratigraphies were used to assess the relative strength of temperature and trajectory controls on δ¹⁸O. There are several potential sources of error in the methodology used to locate individual storm systems within winter snowpacks. These include: non-continuous sampling of snow-pit walls, which means that we may have missed some features of the δ¹⁸O stratigraphies; the spatial variability in both the δ¹⁸O records and rates of snowpack densification; and the fact that the factors used to upscale forefield SR-50 records to the glacier sites are different for every storm. Because of these variable sources of uncertainty, it is difficult to assign a quantitative error to the mean storm δ¹⁸O values. Sampling and scaling errors associated with the determination of discrete accumulation-event layers are likely to be large relative to the analytical error (0.2‰). However, these sources of uncertainty are expected to be random and unbiased, making mean storm δ¹⁸O values amenable to regression analyses. Since a strong temperature control on δ¹⁸O is expected, this analysis provides a means of assessing the effectiveness of this novel methodology to locate storms within snowpacks based on δ¹⁸O stratigraphies.

Figure 7 shows the relationship between mean local storm temperature and mean δ¹⁸O for each accumulation event in each season, and for all data combined. The relationship between these datasets is significant in all cases at 95% confidence, with adjusted R ² values ranging from 0.23 to 0.76. The higher R ² values and regression line slopes at Opabin Glacier indicate that temperature is a better predictor of snowpack δ¹⁸O at this site (Fig. 7c and e) than at Haig Glacier, where the δ¹⁸O–temperature relationship varied by a factor of ∼2 over the 3 years of study (Fig. 7b, d and f).

Fig. 7.

Correlation between the mean temperature and mean snow-pit δ¹⁸O for each accumulation event for (a) all HG and OG data (2004/07) and (b–e) each site in each accumulation season. Adjusted R ² values are given at 95% confidence.

The δ¹⁸O–temperature relationship is very consistent at Opabin Glacier, though there are only 2 years of data so this is a tentative conclusion. It is nonetheless in line with our expectations that continental influences such as air-mass mixing and upslope precipitation events, felt more strongly at Haig Glacier, confound isotope signals. This is supported by the results of regressing δ¹⁸O and temperature for each storm class (Table 3). Temperature is a significant predictor of δ¹⁸O for westerly, southwesterly and northwesterly classes, with adjusted R ² values ranging from 0.34 to 0.55. Upslope storms show no significant relationship between δ¹⁸O and temperature, and the regression equation has a low slope compared with other classes. Given that there were no upslope events in 2005/06 (Fig. 4), it follows that the δ¹⁸O–temperature slope at Haig Glacier was close to that of Opabin Glacier that season, but it was weakened by the occurrence of upslope events in the other two seasons (Fig. 7).

Table 3.

Results from regressing δ¹⁸O and temperature for each storm class. Bold p values indicate a significant correlation at 99% confidence β and values give the unstandardized slope of each regression equation. Also shown are the mean and standard deviation of δ¹⁸O for each storm class

Despite the strong relationships between δ¹⁸O and temperature shown in Figure 7, there are a number of outlying storms at both sites in all years of the study. Assuming that these accumulation events are properly placed in the isotope stratigraphies, it is apparent that local temperatures have little bearing on the δ¹⁸O of some accumulation events and that the δ¹⁸O values of these outlying storms must be governed by factors other than temperature, such as vapour pathway.

The trajectory classification was used to attempt to resolve some of the variability in δ¹⁸O not attributable to temperature. Figure 8 shows the mean δ¹⁸O of accumulation events from both the HG and OG sites grouped by trajectory. The white box plots show the mean and upper and lower percentiles for these data, and the grey box plots show the same statistics for temperature in each trajectory class. Figure 8 again illustrates the distinct temperature control on δ¹⁸O, with colder storms producing, overall, the most negative mean δ¹⁸O values. While this is a reflection of the results presented in Figure 7, grouping events by trajectory also shows a distinct split in the mean and variation in δ¹⁸O between westerly/southwesterly events and northwesterly/upslope events (Table 3). Westerly and southwesterly events have less variability and higher mean δ¹⁸O values than northwesterly and upslope classes. Table 4 shows the results of paired t tests performed for each combination of storm classes. Only westerly/southwesterly events and northwesterly/upslope events have statistically similar mean δ¹⁸O values at 95% confidence, supporting the conclusion that there is a significant difference between all other groupings.

Fig. 8.

Mean δ¹⁸O for each accumulation event grouped by trajectory class and season. Box plots show the mean and 25th and 75th percentiles, and whiskers show the 90th and 10th percentiles, for all data from each class. The white boxes are δ¹⁸O and the grey boxes are temperature. Each square represents one accumulation event and these are grouped by class and season.

Table 4.

Results of paired t tests for mean δ¹⁸O values and residual δ¹⁸O from the temperature–δ¹⁸O regression for each combination of storm class. Bold p values indicate a significant difference at 95% confidence for a two-tailed test

It is not clear from these results whether the statistical difference between classes results from vapour transport per se or from the unique temperature characteristics of each class. Figure 8 shows that the separation in mean δ¹⁸O between different storm classes is at least partly driven by temperature. Mean accumulation-event temperatures were higher for westerly and southwesterly events than for northwesterly and upslope classes. However, while mean temperatures for westerly and southwesterly events were within 1°C at both field sites, there was a greater separation in mean temperatures between upslope and northwesterly classes (Fig. 8; Table 2). Northwesterly storms exhibited higher condensation temperatures than upslope events; but were associated with snow that had similar mean δ¹⁸O values to colder upslope events. The comparatively low δ¹⁸O of snow derived from northwesterly trajectories may be due to the long travel distance and extensive rain-out experienced by water vapour during these weather systems.

To attempt to extract a vapour-trajectory effect from the temperature effect for each class, we analyzed the residuals of the δ¹⁸O–temperature regression for all storms. These are considered to reflect the degree of variation in δ¹⁸O not controlled by temperature. Paired t tests were again performed for each combination of storm class using the residual δ¹⁸O data to reveal any statistical differences between storm classes. In all cases, the classes were found to be statistically similar, indicating that none of the residual δ¹⁸O variability could be differentiated in terms of vapour trajectory alone (Table 4). Finally, a multivariate regression was performed with δ¹⁸O as the dependent variable. Mean temperature and storm class were entered as independent variables. Although westerly, northwesterly and upslope classes were entered into the model, temperature was the only significant predictor, and the model returned an adjusted R ² of 0.46, which is no higher than the result obtained from simply regressing δ¹⁸O against temperature only (adjusted R ² = 0.47).