Data collection at Haig Glacier took place during the 2017 melt season, consisting of a snow survey in May and four field visits, each 2 weeks apart, in July and August (see Table 1 for dates). The May snow survey included probing of snow depth at 33 sites along the glacier centerline. At three of these sites, snow density measurements were collected from snow pits dug through to the underlying glacier ice. Snow densities were measured at 10-cm intervals using a 100-cm³ box cutter. Average values of snowpack density from the snow pits were used to interpolate density to probed sites for calculation of snow water equivalent (SWE).

The winter snow distribution is used to seed a distributed energy balance model of glacier melt (Ebrahimi and Marshall, 2016), using meteorological forcing from the forefield AWS, which is at an elevation of 2,340 m. The AWS measures incoming and outgoing longwave and shortwave radiation, temperature, wind speed and direction, humidity, air pressure, rainfall, and snow depth. From 2001–2015, an additional AWS was maintained at an elevation of 2,665 m on the glacier surface. This period of overlap informs transfer functions to map off-glacier AWS data onto the glacier surface (Marshall, 2014), providing the meteorological fields that are needed to drive the energy balance model. The data are further distributed as a function of elevation on the glacier surface using lapse rates for temperature, humidity, and incoming longwave radiation (Ebrahimi and Marshall, 2016). The distributed energy balance model runs on 30-min time steps over grid cells of ∼30 m, with outgoing longwave and shortwave radiation calculated locally as a function of surface temperature and albedo (Marshall, 2014). Surface albedo evolves in space and time as a function of surface type (snow vs ice), with melt-season reductions in snow albedo calculated as a function of cumulative positive degree days (Marshall and Miller, 2020). The energy balance model provides an estimate of daily meltwater production and runoff, as well as the transition from seasonal snow cover to exposed glacier ice across the glacier.

During each July and August field visit, samples of Haig Glacier runoff water were collected from the stream site every 2 hours into plastic 20-mL sampling bottles over a 24-h period starting around 8 a.m. each morning and continuing until the same time the following day. The water samples were filtered through a 0.20-µm filter to remove particulates.

Glacier surface snow/ice samples were collected at eleven of the sampling points (Figure 1) on the same days as stream samples, approximately between 9 a.m. and 3 p.m. Snow/ice samples were collected using a trowel into freezer bags. The surface layer (approximately 1 cm) was scraped clear and samples represent depths of 1–5 cm in the snow and ice. Samples were melted upon return to the base camp, and transferred into 20-mL plastic bottles.

Additional water samples were collected on 14 October 2017 from the proglacial stream. The stream was still flowing at this time, though it was covered by a thin layer of ice. There was no supraglacial meltwater generation or runoff at this time, since the surface of Haig Glacier was snow-covered, the waterfall at the glacier terminus was frozen (though with water running behind it), and temperatures had been below 0°C for several days. The ice cover on the stream will have eliminated atmospheric snow/rainfall inputs as well as evaporation. The sample collected from underneath the stream ice cover is therefore interpreted to be subglacial outflow, which may have also passed through the shallow groundwater system. The meltwater likely originated in summer 2017, but spent a period of many days to weeks in transit through the glacier drainage system, passing through the subglacial environment.

Data from 2004–2007 (Sinclair and Marshall, 2008) and are used in this study to place 2017 summer samples in context with winter and spring samples. Data from 2004–2007 are divided into dry and wet snow samples. Samples considered “dry snow” were collected in the months of January through early May, corresponding to measured snow temperature below 0°C. Samples considered “wet snow” were collected from May through August, when either snow temperatures of 0°C were measured or liquid water was observed in the snowpack.

Laboratory analyses of major ions and stable isotopes of the proglacial stream and glacier surface samples were conducted at the University of Calgary. Stable isotopes of oxygen (δ¹⁸O) and hydrogen (δ²H) were analyzed at the Isotope Science Laboratory. Since snow/ice samples from the glacier surface were not filtered in the field, filtration of these samples took place before each sample was pipetted into a 1.5-mL glass test vial. The samples were run on a Los Gatos Research Isotopic Water Analyzer (model IWA-45EP), and the delta value was reported in per mil (‰) with a precision of ±0.1‰ for δ¹⁸O and ±1.0 for δ²H relative to Vienna–Standard Mean Ocean Water (V-SMOW).

Major ions, including sulphate (SO₄ ^2-), were analyzed on a Metrohm 930 Compact Iron Chromatography (IC) Flex system connected to a Metrohm 858 Professional Sample Processor (an autosampler). The concentration was reported in mg L^-1 with a precision of 5% and a detection limit of 0.10 mg L^-1. The bedrock underlying Haig Glacier is mainly dominated by sandstone. Pyrite (FeS₂) and gypsum (CaSO₄) are often found in sandstone. Sulphate (SO₄ ^2-) is produced by pyrite oxidation and gypsum dissolution making SO₄ ^2- the best ionic fingerprint for water that has been in contact with the bedrock.

Snow depths ranged from 2.84 m on the lower glacier to 4.70 m at the continental divide site in the upper accumulation area in May 2017, with average snowpack densities from 395 to 440 kg m⁻³. Peak annual accumulation typically occurs in late May. Averaged over the glacier, winter specific mass balance was B _w = 1.49 m w. e., about 15% above normal (Adhikari and Marshall, 2013). The summer of 2017 was warm and dry, with sustained temperatures above 0°C from June through August and few precipitation events. Warm conditions drove extensive melting, with the snowline retreating to the upper accumulation area by mid-August. Modelled melt season (May through September) runoff was B _s = 3.31 m w. e., giving a net annual mass balance of B _n = −1.82 m w. e. Geodetic (LiDAR) mass balance measurements from 2017 are in good agreement with the glaciological (field- and model-based) estimates, giving B _n = −1.91 ± 0.21 m w. e. (Pelto et al., 2019).

Figure 4 shows the modelled daily surface energy balance and meltwater production from May 1 through 16 September 2017, at which point fresh snow began accumulating on the glacier and the summer melt season abruptly ended. Within the model, meltwater runoff can be partitioned into seasonal snow vs the underlying glacier ice and firn. The firn area of the glacier is small, roughly 5% of the total glacier area, so the combined meltwater from firn and ice is dominated by the component from exposed glacier ice. In 2017, snow melt accounted for 1.25 m w. e. (38%) of the summer runoff, while ice and firn melt made up 2.06 m w. e. (62%). Based on model results, the transition from snow-to ice-dominated runoff took place around July 13 (Figure 4D). Snowmelt from the upper glacier continued to contribute to runoff through the rest of the summer, but at low levels after the first week of August.

Table 1 reports the estimated meltwater runoff derived from seasonal snow vs glacier ice and firn for the four sample collection periods. All four periods were warm (average daily temperatures of 8°C–10°C on the glacier), with meltwater production of 19–62 mm w. e. d⁻¹, averaged over the glacier. There was extensive seasonal snow cover on the glacier during the first sampling visit (July 13–14), with roughly equal amounts of snow and ice melt on these days. The glacier albedo was higher during this visit and conditions were overcast for much of the period, so net energy and meltwater production were close to the average daily values on Haig Glacier for summer 2017, despite the warm conditions. Snow cover was largely depleted on the remaining three visits. Combined with clear-sky conditions and warm temperatures, net energy levels and melt rates were high (roughly twice the mean daily values for summer 2017).

In earlier work at Haig Glacier, Shea et al. (2005) report time lags of ∼3 h for meltwater runoff to reach the stream site during the late summer, when glacier ice is exposed and crevasses are open to funnel the water efficiently through the glacier drainage system. Daily meltwater production therefore provides a reasonable estimate of total daily runoff in the proglacial stream, but meltwater storage and delays in runoff can occur, particularly from the upper glacier, and base flow may also contain older water. The surface meltwater values (Table 1) therefore do not precisely reflect the amount of meltwater flowing through the stream site on a given day.

The δ¹⁸O and δ²H values of snow, ice, and meltwater runoff measured at Haig Glacier are averaged over multiple years (Figure 8). AWS temperature and snow accumulation (winter mass balance) data indicate interannual variability but no statistically significant trends in temperature or precipitation at the site from 2001–2017. Insofar as these are the main expected controls on isotopic ratios within the snowpack, we argue that snow isotopic values from different years can be compared and combined. While there is no expectation or evidence of changes in the mean annual precipitation isotopes with time, we observe a systematic shift in isotopic ratios associated with post-depositional processes on the glacier. Specifically, the average δ¹⁸O value increases from −20.1‰ for dry snow to −19.2‰ for wet snow and −18.6‰ for glacier ice, the same trend is observed for δ²H values.

The δ¹⁸O in the Haig proglacial stream varies diurnally (Figure 7). Runoff near the end of the melt season (August 22) is most enriched compared to the earliest sampled runoff (July 13), which is most depleted in ¹⁸O (the same trend is observed in deuterium ratios, see Figure 5). Meltwater runoff had a distinct peak in δ¹⁸O in the mid-afternoon of July 26, indicating enrichment in ¹⁸O at this time of high meltwater production. Stable water isotope ratios in the proglacial stream on each sampling day are summarized in Figure 5.

In the stream, average values of δ¹⁸O and δ²H increased as the melt season progressed, indicating a systematic enrichment in heavy isotopes through the melt season. Shapiro-Wilk normality test showed that the data for each stream sampling interval are normally distributed. Analysis of Variance statistical tests were used to compare the mean δ¹⁸O values for all four stream-sampling intervals. In all cases, statistical tests have values of p < 0.001, which indicate that the river samples from each day are from different populations, and confirming the statistically-significant change in isotopic values of the glacier runoff over the melt season.

The δ¹⁸O- δD relationship is examined for Haig Glacier through a series of meteoric water lines (MWL). Dry snow samples (sampled in winter and early spring, with snow temperatures below 0°C) are assumed to represent local Haig Glacier precipitation, and have been used to construct a local meteoric water line (LMWL)for the study area (Figure 6). Meteoric water lines have also been constructed for wet snow samples (snow temperatures of 0°C) from 2004–2007 and 2017 (Figure 6B), 2017 ice (Figure 6C), and the 2017 proglacial stream (Figure 6D).

Table 2 summarizes the MWL parameters at Haig Glacier. From the dry-snow samples, the slope of the LMWL is 7.74. The wet snow from 2004–2007 and the 2017 field seasons has a steeper slope of 7.47. MWL values for the glacier ice sampled in 2017 are significantly lower, 6.98 and the MWL for the proglacial stream samples was equal to 7.83, similar to the wet-snow value.

We also calculated d-excess values for the proglacial stream and the sampled water sources. As summarized in Figure 5, in the proglacial stream the d-excess varies between 8.4‰ and 8.6‰ for the first three sample periods and increases to 9.3‰ for 22 August. The average summer d-excess in the proglacial stream was 8.7‰, similar to the value for wet snow (8.8‰). Glacier ice (6.4‰) and dry snow (10.2‰) bracket the wet-snow value (Table 2), indicating a reduction in d-excess through the transition from dry snow to wet snow to glacier ice.

The diurnal variation of stream δ¹⁸O values (Figure 7) indicates the role of different source water in contributing to stream runoff. This is further supported by ion concentrations measured in the stream (Miller, 2018), in particular sulphate (SO₄ ^2-), which is present in the glacier runoff but undetectable in the glacier snow and ice. SO₄ ^2- is in the bedrock underlying the glacier and any water in contact with the bedrock acquires a sulphate signal. A subglacial water sample collected on 14 October 2017 from the proglacial stream has the highest SO₄ ^2- (26.75 mg L⁻¹) concentration of all the stream samples. The subglacial water is rich in numerous other dissolved ions, such as calcium and magnesium, but sulphate is the most pronounced in the suite of measured ions, making it the best tracer for subglacial water at this site. This single sample does not capture temporal variability of SO₄ ^2- concentration in the Haig Stream. SO₄ ^2- concentration was measured in surface ice samples and the proglacial stream in the summer, with average concentration of <0.1 mg L^-1 (n = 44) and 2.54 ±1.3 mg L^-1 (n = 57), respectively. The concentration of SO₄ ^2- in surface ice and Haig Stream when it is diluted with meltwater is very low compared to the autumn sample. Therefore, we assume that the measured concentration of SO₄ ^2- is a fingerprint of subglacial water.

The results of the mixing analysis indicate that water from the subglacial component in the stream increases from 5.5% in early July to 10.5% in late August. This signals an evolution of water routing through the glacier. As the melt season progresses, pathways open up for meltwater generated at the surface to connect to the subglacial environment. While the contribution of subglacial water to the proglacial stream increases, the stream volume is still dominated by supraglacial water during the late melt season, or water that passes quickly through the englacial and subglacial systems, with limited contact with the bed.

We consider the potential enrichment of surface snow associated with evaporation of liquid water in the pore space or ponded on the glacier surface. Assuming an initial isotopic value equal to the dry-snow value, δ 0 18 O = −20‰, Equation 4 gives a fractionation factor of 0.92. This indicates that 8% of the water volume must evaporate to change the δ¹⁸O value from −20‰ to −19‰ on the glacier surface. We assume that evaporation occurs within the near-surface snow, with well-mixed liquid water and water vapour in the upper 10 cm of the snow cover. The snow water volume in the July snowpack at Haig Glacier is estimated at 4%–8% by volume (Samimi and Marshall, 2017), so the top 10 cm of wet snow holds 4–8 mm of water. For a liquid water content of 6% (6 mm), 0.48 mm must evaporate to shift the δ¹⁸O value by 1‰.

The average summer (JJA) net energy, Q, available for evaporation or melting on Haig Glacier in summer (JJA) 2017 was calculated to be 112 W m⁻². If all of this energy was directed to evaporation, the rate of evaporation E) would be 4.5 x 10⁻⁸ m s⁻¹. Based on this, 3 h would be required to evaporate 0.48 mm of water. Not all available net energy is used for evaporation; most is used for melting. If 10% of the available energy is consumed by evaporation (i.e., for an evaporative energy flux of −11 W/m², typical of early summer), this would require ∼30 h of active evaporation (a few days) to produce enough evaporation to explain a 1‰ enrichment of the residual water.

Dry snow spans a large range of δ¹⁸O (Figure 8), as a function of the weather system that conveys moisture to the study site (Sinclair and Marshall, 2009). Wet snow δ¹⁸O values are less varied, and this pattern continues for glacier ice and stream samples, which have low variability compared to dry snow. Homogenization of the snowpack at the onset of the melt season at Haig Glacier was previously reported by Sinclair and Marshall (2008), consistent with the observations of Raben and Theakstone (1998). The progressive enrichment in heavy isotopes in meltwater runoff over the course of the melt season echoes the results of Taylor et al. (2001) and Williams et al. (2006). This systematic seasonal enrichment is also recorded in the supraglacial snow and ice at Haig Glacier, with an increase in δ¹⁸O of ∼1.0‰ from dry to wet snow and an additional increase of ∼0.5‰ from wet snow to glacier ice (Figure 8). The latter shift is identical to that reported by Theakstone and Knudsen (1996), who document δ¹⁸O values of −11.6‰ for glacier ice, compared with −12.1‰ for July snow.

The glacier ice that we sampled is not contemporaneous with the seasonal snowpack. It takes several years for the snow to transition from firn to ice at this site, and glacier ice is then transported downslope over timescales of decades. Mechanical deformation of the firn and ice will not change its bulk isotopic composition, but multi-year firn and ice have had more time and exposure for isotopic exchange in freeze-on processes that we describe in Section 5.3. This may explain the enrichment in heavy isotopes (see below). It is also possible that the original isotopic ratios of the glacier ice were different, since it developed from snow that was deposited on the glacier at higher elevations and in prior decades. The glacier is low-sloping and slow-moving (velocities of about 5 m yr⁻¹ reported by Adhikari and Marshall, 2013), with a total elevation range of about 325 m, so the altitude effects will be minor and surface ice that we sampled may not have travelled far from its source. Moreover, ice that originated higher up on the glacier and in earlier (cooler) decades is expected to be more depleted in heavy isotopes; therefore, the inferred ∼0.5‰ enrichment in heavy isotopes through post-depositional modification is a minimum, and the enrichment processes are enough to overcome the potential source effects that act in the opposite direction.

Proglacial stream δ¹⁸O values are intermediate between those of snow and ice, but the progressive isotopic enrichment over the course of the melt season (Figure 5) is consistent with the seasonal evolution of the meltwater source. At the beginning of the sampling period, most runoff from the glacier is from melting snow. As the melt season progresses and glacier ice becomes exposed, the more enriched ice contributed an increasingly larger fraction to the meltwater, giving isotopically heavier stream water at the end of the sampling period (late August).

The enrichment of ice in δ¹⁸O compared to supraglacial snow helps to explain the seasonal evolution of δ¹⁸O in the meltwater runoff, consistent with previous studies, but the reason for the isotopic enrichment in wet snow and ice is uncertain. We suggest two potential mechanisms for enrichment of ice compared to supraglacial snow after the onset of melt.

The enrichment of oxygen isotopes in ice relative to snow is consistent with evaporative processes that take place from the liquid water fraction as snow melts. Wet snow contains liquid water on the surfaces of snow grains and in the pore spaces of the snowpack, and liquid water can also pool at the snow-ice interface. The wet snowpack eventually melts and drains away, leaving the ice surface exposed, but liquid water commonly pools on the glacier surface where surface slopes are low. Because lighter isotopes are favoured during evaporation, residual water remains enriched with heavier isotopes of oxygen and hydrogen as evaporative processes take place (Gonfiantini et al., 2018). Wet snow and meltwater that remain on the glacier past the summer melt season turn to firn and eventually to glacier ice, leaving ice enriched relative to the meteoric snowfall.

The analysis of local meteoric water lines shows a reduction in the MWL slope for the glacier ice, relative to the 2017 snow samples (Table 2), which may be indicative of post-depositional kinetic fractionation processes during evaporation. Zhou et al. (2013) documented a decrease in slope from precipitation to glacier firn to glacier runoff on the Tibetan Plateau, and attributed this in part to kinetic evaporative processes. Where the residual water eventually contributes to the glacier ice, evaporation leads to enrichment of δ¹⁸O isotopes in glacier ice compared to glacier snow, as well as a lower MWL slope (Williams et al., 2006; Zhou et al., 2014; Wang et al., 2015).

We numerically illustrate a scenario in which evaporation leads to a 1‰ change in δ¹⁸O values of liquid water on the glacier surface over the course of 30 h. This amount of evaporation is readily achieved in the spring and early summer, when near-surface water in the snow aquifer can have a residence time of days to weeks. Evaporation is likely reduced once the seasonal snow cover is depleted in late summer, since supraglacial drainage is more efficient on the exposed ice and the residence time of water declines.

Evaporation can explain the enrichment in heavy isotopes seen in the meltwater runoff, but it does not necessarily account for the observed enrichment in the snow and ice itself. Enriched meltwater that refreezes or is retained within the surface snow and ice will cause an enrichment of the snow and ice matrix, but much greater amounts of evaporation would be needed to shift the bulk matrix fraction of the snowpack by 1‰ through this mechanism.

Freeze-thaw processes offer an alternative mechanism to explain the isotopic enrichment of wet snow and glacier ice. When there is meltwater refreezing and drainage, equilibrium fractionation processes cause preferential refreezing of the heavier isotopes, leaving lighter isotopes in the residual pore water of the snowpack (the opposite effect to evaporation). If the isotopically lighter residual water drains, this results in an overall enrichment of heavy isotopes in the remaining bulk snow (Taylor et al., 2001; Ohlanders et al., 2013). The initial meltwater runoff should therefore be isotopically depleted, but as the melt season progresses the enriched snow and ice is melted, and runoff should become progressively more enriched in heavy isotopes.

Williams et al. (2006) initially proposed this idea to explain the seasonal evolution of δ¹⁸O observed in the runoff from Green Lake 5 rock glacier in Colorado, but specific to the englacial environment of the rock glacier. When liquid water infiltrates to the subsurface, any freeze-on during contact with englacial (or buried) ice leaves the ice enriched in heavier isotopes compared to the meltwater runoff and the initial surface snow. Williams et al. (2006) also report higher values of deuterium excess in late-season runoff relative to the initial conditions. Taken together with the enrichment in heavy isotopes, they interpret this as an indication that meltwater from the buried glacier ice has undergone multiple freeze-thaw episodes.

Consistent with the observations of Williams et al. (2006), d-excess in the Haig Glacier proglacial stream increased over the course of the 2017 melt season, from 8.6‰ to 9.3‰ (Figure 5). However, there was no progressive increase in d-excess for the supraglacial snow samples collected at Haig Glacier in July and August 2017, and glacier ice values had significantly lower d-excess (6.4‰) than the snow samples (8.8‰) (Figure 8). This could reflect the effects of evaporation in wet snow and during the multi-year transition from firn to ice; kinetic fractionation processes during evaporation favour ²H over ¹⁸O (e.g., Benetti et al., 2014), which increases d-excess in water vapour while leaving the residual water in the snow/ice reservoir with lower values of d-excess. Additionally, this argues against the Williams et al. (2006) freeze-thaw process as an explanation for the enrichment in heavy isotopes in glacier ice, relative to the supraglacial snowpack. D-excess may be a potential tracer of glacial melt, and it may be worth considering river water deviations from the local meteoric water line, a method known as line-conditioned excess, lc-excess, after Landwehr and Coplen (2006). Future work to trace these signals could be rewarding.

This process is most effective if snow, firn, or ice are below 0°C, to promote refreezing. This is not the case at Haig Glacier, so englacial and subglacial refreezing are unlikely, but diurnal freeze-thaw cycles occur at the glacier surface through much of the summer, with overnight refreezing of ponded surface water and near-surface pore water (Samimi and Marshall, 2017). These freeze-thaw cycles may lead to preferential freeze-on of heavy isotopes and removal of lighter isotopes via drainage. Fractionation by this process could also occur through the depth of the snowpack during spring, as the snow warms to an isothermal state and surface meltwater infiltrates and refreezes. Glacier energy balance modeling from summer 2017 indicates that there were an estimated 108°days with freeze-thaw cycles on the glacier from 1 May through 30 September, with 89°days on the upper glacier and 119°days at the lowest elevations. This is based on days with melting where the net energy Q _N was negative for 1 hour or more, most commonly associated with overnight refreezing. The supraglacial freeze-thaw cycles could explain the progressive seasonal enrichment of the residual snow and glacier ice, assuming that some of the isotopically depleted liquid meltwater drains (leaves the system) through percolation and runoff.

This process could also be occurring in conjunction with daytime liquid water evaporation, which would explain the mixed signals within the supraglacial water system and meltwater runoff. Freeze-on processes deplete this water in heavy isotopes, while evaporation enriches it in heavy isotopes, with these processes partially offsetting each other. The seasonal evolution of the stable isotopes in the meltwater runoff largely reflects the melt source (snow vs ice), accounting for the progressive seasonal enrichment.