Introduction
Stable isotopes of oxygen and hydrogen in precipitation are strongly correlated with local air temperature (Reference DansgaardDansgaard, 1964), providing the basis for temperature reconstructions from both polar and lower-latitude ice cores (e.g. Reference DansgaardDansgaard and others, 1993; Reference Thompson, Mosley-Thompson, Davis, Lin, Henderson and MashiottaThompson and others, 2003). However, one of the key barriers to reconstruction of paleotemperatures from cryospheric records is the effect of changing vapour source regions and air-mass trajectories on stable water isotopes, δ18O and δD (e.g. Reference Charles, Rind, Jouzel, Koster and FairbanksCharles and others, 1994; Reference JouzelJouzel and others, 1997). Air-mass trajectories control the distance traversed by atmospheric water vapour and the nature of the topography it encounters prior to the eventual precipitation of this moisture over land. These factors in turn control the degree of distillation, or rain-out, of heavy isotopes from the vapour mass and have a marked influence on the isotopic composition of the resulting precipitation (Reference DansgaardDansgaard, 1964). In order to reconstruct a temperature history from isotopic records, it is important to extract the effect of changing vapour trajectories from the δ–temperature relationship. This is particularly relevant for the interpretation of low- and mid-latitude ice cores, as these show a high degree of spatial variability in δ–temperature relationships on annual to millennial timescales (Reference Thompson, Cecil, Green and ThompsonThompson, 2004).
One approach to deciphering the relationships between air-mass trajectories and their isotopic fingerprint in cryospheric records is to consider their interactions on a seasonal timescale. Previous work in Antarctica (Reference Isaksson and MelvoldIsaksson and Melvold, 2002; Reference McMorrow, Curran, van Ommen, Morgan and AllisonMcMorrow and others, 2002) compares glaciochemical records in polar snow and firn with contemporaneous meteorological observations to locate regional accumulation events in snowpack stratigraphies. Reference Helsen, van de Wal, van den Broeke, van As, Meijer and ReijmerHelsen and others (2005) and Reference Shuman, Bromwich, Kipfstuhl and SchwagerShuman and others (2001) related seasonal snow-pit isotope data to temperature records in the Antarctic and Arctic, respectively, while Reference Schlosser, Reijmer, Oerter and GrafSchlosser and others (2004) considered the influence of precipitation origin on the δ18O–temperature relationship at Neumayer station, Antarctica, finding that the quality of the relationship depended on trajectory class. However, the linkages between vapour trajectories, snow accumulation and precipitation isotopes have not been explored to a great extent in mountainous mid-latitude regions, where vapour trajectories and source regions are highly variable (Reference Vuille, Hardy, Braun, Keimig and BradleyVuille and others, 1998).
The main objective of this research is to extract the effect of vapour trajectory from the dominant influence of condensation temperature on stable isotopes in winter snowpacks in the Canadian Rocky Mountains. We approached this using alpine snow-pit data from three consecutive accumulation seasons, based on the following steps: (1) determine the isotopic characteristics of these snowpacks; (2) locate each accumulation event from each season in an isotope stratigraphy to obtain mean δ18O, δD and deuterium excess values; (3) produce a trajectory classification for all events; (4) investigate the meteorological characteristics of each trajectory class; and (5) examine the effect of temperature vs trajectory class on stable-isotope ratios.
It has been shown that changing atmospheric circulation patterns in central North America from the Pleistocene to the present have influenced the isotopic composition of precipitation (Reference Amundson, Chadwick, Kendall, Wang and DeNiroAmundson and others, 1996; Reference Edwards, Wolfe and MacDonaldEdwards and others, 1996). However, there has been surprisingly little investigation into the effect of intra-seasonal changes in storm trajectories on the isotopic composition of precipitation in this region. Given that each dominant storm trajectory is likely to have a different isotopic signature, the number of storms of a given trajectory and the accumulation associated with these events will have a marked influence on mean δ values in the end-of-season snowpacks. Understanding these relationships can enable the interpretation of the complete winter snowpack stratigraphy in terms of weather systems that dominate the seasonal moisture supply, offering a new tool for glacial mass-balance and water-resource forecasting. Furthermore, understanding these relationships will be valuable in terms of modelling historical vapour trajectories and developing robust temperature reconstructions for future paleoclimatic investigations in western Canada.
Study Area and Data Collection
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 km2 (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 km2) 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 km2. 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 cm3 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 CO2–H2O equilibration technique for δ18O 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 δ18O and ±2.0‰ for δD.
Methodologies
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), ρ
1 is the density of each layer (kg m−3) and ρ
w is the density of water (1000 kg m−3). An estimate of the mean fresh-snow density from each snow pit, 
, was obtained from the accumulated SWEf 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−3, 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 δ18O 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 δ18O 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 δ18O 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 δ18O 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 δ18O 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 δ18O 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. 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. 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.
Results and Discussion
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
Westerly
Westerly storm systems are characterized by strong zonal flow across western Canada (Fig. 5a) and a high precipitation gradient between the glacier sites and stations further to the east (Table 2). It was common for Calgary to receive only trace precipitation during westerly storm events, while Haig Glacier received in excess of 20 cm of accumulation. A lag in storm onset was observed for most storms, with Opabin Glacier often receiving snowfall 12–24 hours earlier than Haig Glacier, which is consistent with the position of these field sites along westerly storm trajectories (Fig. 3). The mixing ratio at Kelowna was highest during westerly storms, averaging 2.9 g kg−1. There also tended to be little change in mixing ratio in Kelowna during westerly storms, reflecting the strong influx of Pacific moisture over the duration of these events.
Temperatures were mild at both field sites during westerly events, averaging −4.9°C and −5.6°C at Haig and Opabin Glaciers, respectively. Conditions in Calgary and Banff were generally warm and dry during westerly flow, as expected when chinook (föhn) winds prevail in the lee of the Rocky Mountains. However, in some cases this situation was reversed and the continental divide sites would be warmer than Calgary during westerly storms, despite being ∼1300 m higher. This reflects the climatic divide that frequently sets up in the region, separating mP and cP air masses.
Southwesterly
The broad-scale synoptic pattern associated with southwesterly flow over western Canada is a trough descending over the northern Pacific and a ridge of high pressure to the south of the Canadian Rockies (Fig. 5b). This synoptic set-up can draw large amounts of tropical or subtropical (mT) moisture towards western Canada on what is known as the ‘Pineapple Express’. The other synoptic mechanism that produces southwesterly trajectories is coastal low-pressure systems that develop near the Aleutian Islands in the Gulf of Alaska and travel south down the Pacific coastline. When these are situated in the vicinity of Vancouver Island, strong southwesterly flow and heavy precipitation develop over southwestern Canada.
This trajectory class is associated with the highest average precipitation at both field sites, and heavy precipitation to sites east of the continental divide (Table 2). Moisture advection over Kelowna was correspondingly high during these storms, and the mixing ratio increased by 14% over the storm durations, on average, in contrast with the drying that accompanied all other trajectory classes. The average temperature during southwesterly storms was mild at both field sites, but these events were associated with strong temperature decreases: −6.9°C and −4.5°C at Haig and Opabin Glaciers, respectively. This was also the case in Calgary and Banff, where temperatures dropped by >9°C during these storms.
Northwesterly
Northwesterly trajectories are associated with meridional flow over western Canada. They generally have an inverse synoptic structure compared with southwesterly storms, with an upper-level ridge over the Pacific and an inland trough of low pressure, associated with Arctic (cP) air (Fig. 5c). These storms have a strong northerly component and moderate amounts of accumulation (Table 2), depending on the relative strength of the ridge vs the trough. If the trough gains strength during the course of an accumulation event, upper-level flow swings further north, effectively cutting off the supply of moisture from the west. The storm may then develop an upslope (easterly) component, either for the remainder of the event, or until the trough weakens and moisture can again reach the Rocky Mountains.
Haig and Opabin Glaciers received considerably less snow from northwesterly trajectories than from other classes. At Haig Glacier, for example, these events brought, on average, 40% and 52% less snow than westerly and southwesterly events, respectively. Calgary, in contrast, received more snow during northwesterly events than during westerly and southwesterly events. This is likely related to upslope snow that accumulates on the eastern side of the Rockies but does not reach the continental divide. Northwesterly storms were also associated with the lowest average mixing ratio (2.3 g kg−1), and the strongest drying (−25.4%) of all classes.
Temperatures were lower at Haig and Opabin Glaciers during northwesterly storms, relative to westerly and southwesterly events. Temperatures at Calgary and Banff were also lower than for westerly events. Temperature decreases seen at eastern slope sites during these snow events were not observed at the Opabin Plateau, indicating that the Arctic air descending over Alberta during these events generally does not cross the continental divide.
Upslope
Winter upslope events in western Canada are commonly associated with the southward penetration of an upper-level trough of low pressure. This cold air pools across Alberta and often spills into British Columbia, bringing frigid temperatures and widespread accumulation (Fig. 5d). These air masses do not bring significant amounts of water vapour with them, but interact with Pacific moisture to bring sustained snowfall to the Rocky Mountains.
Accumulation totals in Calgary and Banff were the highest for upslope events. These systems were also associated with high accumulation (22.6 cm) at Haig Glacier, but Opabin Glacier received less snow from upslope storms than from westerly and southwesterly events, which indicates that these storms weaken considerably on the western side of the continental divide. The influence of cold Artic air during these storms is also reflected in the very low average temperatures at all sites and large temperature decreases, particularly at sites east of the continental divide.
The isotopic characteristics of Rocky Mountain snowpacks
The δ18O 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 δ18O and δD values than expected for its latitude.
Figure 6a shows the density-weighted mean δ18O 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 δ18O 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 δ18O 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 δ18O 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 δ18O 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 δ18O at the HF site in 2004/05 (HF305) was −21.6‰, while the mean late-season δ18O 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 δ18O 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 δ18O, known as the global meteoric waterline (GMWL), was defined by Reference CraigCraig (1961) as δD = 8δ18O + 10. The slope and intercept of the line of best fit for all δD–δ18O 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δ18O, 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 δ18O values derived using SR-50 records and snow-pit isotope stratigraphies were used to assess the relative strength of temperature and trajectory controls on δ18O. 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 δ18O stratigraphies; the spatial variability in both the δ18O 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 δ18O 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 δ18O values amenable to regression analyses. Since a strong temperature control on δ18O is expected, this analysis provides a means of assessing the effectiveness of this novel methodology to locate storms within snowpacks based on δ18O stratigraphies.
Figure 7 shows the relationship between mean local storm temperature and mean δ18O 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 2 values ranging from 0.23 to 0.76. The higher R 2 values and regression line slopes at Opabin Glacier indicate that temperature is a better predictor of snowpack δ18O at this site (Fig. 7c and e) than at Haig Glacier, where the δ18O–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 δ18O for each accumulation event for (a) all HG and OG data (2004/07) and (b–e) each site in each accumulation season. Adjusted R 2 values are given at 95% confidence.
The δ18O–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 δ18O and temperature for each storm class (Table 3). Temperature is a significant predictor of δ18O for westerly, southwesterly and northwesterly classes, with adjusted R 2 values ranging from 0.34 to 0.55. Upslope storms show no significant relationship between δ18O 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 δ18O–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 δ18O 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 δ18O for each storm class
Despite the strong relationships between δ18O 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 δ18O of some accumulation events and that the δ18O 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 δ18O not attributable to temperature. Figure 8 shows the mean δ18O 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 δ18O, with colder storms producing, overall, the most negative mean δ18O 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 δ18O between westerly/southwesterly events and northwesterly/upslope events (Table 3). Westerly and southwesterly events have less variability and higher mean δ18O 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 δ18O values at 95% confidence, supporting the conclusion that there is a significant difference between all other groupings.
Fig. 8. Mean δ18O 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 δ18O 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 δ18O values and residual δ18O from the temperature–δ18O 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 δ18O 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 δ18O values to colder upslope events. The comparatively low δ18O 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 δ18O–temperature regression for all storms. These are considered to reflect the degree of variation in δ18O not controlled by temperature. Paired t tests were again performed for each combination of storm class using the residual δ18O 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 δ18O variability could be differentiated in terms of vapour trajectory alone (Table 4). Finally, a multivariate regression was performed with δ18O 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 2 of 0.46, which is no higher than the result obtained from simply regressing δ18O against temperature only (adjusted R 2 = 0.47).
Conclusions and Future Directions
The identification of linkages between isotopic stratigraphies in alpine snowpacks and contemporaneous meteorological conditions has the potential to be used in two key ways. Firstly, if the relationships are well understood at a given study site, the isotope stratigraphy of the end-of-winter snowpack can be used to determine which vapour trajectories were dominant in a given season. Knowing which trajectories are most important for snow accumulation can, in turn, help determine how changing atmospheric circulation patterns will impact the spatial and temporal distribution of the winter snowpacks and hence water resources in mountainous regions. The second main application of this type of analysis is to help interpret longer records of climate change. Understanding the transport history of material in firn and ice cores is essential to accurate reconstruction of temperature histories from these records.
This study used data from winter snowpacks over three consecutive accumulation seasons at two alpine field sites to assess the effect of temperature and vapour trajectory on δ18O. The isotope stratigraphies preserve a record of the isotopic fingerprint of individual storm events throughout each season and show a large degree of variability through the winter at each site. This variability reflects the rain-out and temperature history of the succession of winter storms that bring snow to the region. On average, Opabin Glacier snowpacks have consistently lower δ values than Haig Glacier snowpacks, which we suspect to be due to Opabin’s position downwind of the orographic obstacle imposed by Glacier National Park, British Columbia. The lower-elevation OP and HF sites had consistently lower δ values than the glacier sites, which is partially explained by the accumulation of early-season (isotopically heavy) snow at higher sites. Some part of this difference may also be due to post-depositional modification of the snowpack at the forefield sites.
Accumulation events were distributed in snow-pit isotope stratigraphies based on the relative amounts of SWE from each event, and a storm trajectory classification enabled quantification of the relative importance of different winter-storm trajectories for snow accumulation over the study period. The meteorological characteristics of each storm class were defined and the classification appears robust in light of the parameters considered.
The classification also enabled an analysis of the relative importance of temperature and vapour trajectory on the mean δ18O of accumulation events. Temperature was the dominant control on δ values, and described 47% of the variability in δ18O at both field sites, although the strength of this relationship differed on western and eastern sides of the continental divide. At Opabin Glacier, temperature explained about 75% of the variability in δ18O, but at Haig Glacier the influence of upslope snow weakened this relationship. An attempt was made to isolate the effects of vapour trajectory by grouping events by trajectory class and considering the mean and range in δ18O for each class. While there is a high degree of overlap in the range of δ18O between trajectory classes, there is a distinct separation between westerly/southwesterly events and northwesterly/ upslope classes. It appears, however, that the differences in δ18O are strongly related to temperature differences between classes rather than storm trajectories alone.
It is unclear at this stage whether the lack of a vapour-trajectory signal is a result of physical processes, and it is possible that this is masked by the broad classification scheme employed in this study. To further investigate the effect of temperature and storm trajectory on δ values, condensation and rain-out along different air-mass trajectories will be modelled using a single-stage Rayleigh distillation model coupled to a simple model of orographic rain-out along specific pathways in western North America. This approach will be utilized to further test the hypothesis that these pathways influence the isotopic character of the mountain snowpack. Furthermore, the impact of vapour source region on d in winter snowpacks in this region is currently not well understood, and modelling individual storm systems may elucidate the relationships between d and vapour source regions for western Canada.
Acknowledgements
We thank the many people who have assisted us with fieldwork over the course of this study, particularly J.L. Hood, T. Moran and B. Horton. We also acknowledge ongoing assistance and collaboration from M. Hayashi and his research group who provided us with weather station data from the Opabin Plateau and supported our fieldwork at this site. R. Smith at the University of Calgary weather research station and J. Buchanan-Mappin at the Kananaskis Field Station also provided us with meteorological data. Two anonymous reviewers provided many thoughtful comments on the draft manuscript, and we also thank the scientific editor, D. Peel. This study was supported by grants from the Natural Sciences and Engineering Research Council (NSERC) of Canada and the Canadian Foundation for Climate and Atmospheric Sciences.
