1. INTRODUCTION
The Canadian Arctic Archipelago (CAA) hosts the largest volume of ice outside the ice sheets (Pfeffer and others, Reference Pfeffer2014) and occupies latitudes that are currently experiencing some of the greatest rates of climate warming (Sharp and others, Reference Sharp2015, Reference Sharp2011), a tendency that is predicted to continue well into the future (Kirtman and others, Reference Kirtman, Stocker, Qin, Plattner, Tignor, Allen, Boschung, Nauels, Xia, Bex and Midgley2013; Lenaerts and others, Reference Lenaerts2013). A recent analysis of the CAA in situ surface mass-balance records indicates that the average glacier mass balance between 2005 and 2009 was five times more negative than the average from 1963 to 2004 (Sharp and others, Reference Sharp2011). Through modelling and remote sensing, Gardner and others (Reference Gardner2011) showed that 92% of the mass loss from the CAA can be explained by increased melt, while only 8% is attributable to frontal ablation (i.e. calving) between 2004 and 2009.
In terms of measuring glacier mass balance in alpine basins, there are two primary methods:
-
(1) The glaciological method (also often referred to as the direct method) involves interpolation between in situ measurements of accumulation and ablation at stakes drilled into the glacier surface. These stakes are typically located along the glacier centreline, and interpolation across the glacier basin is undertaken either by assuming that mass balance varies only with elevation (the profile method) or by mapping accumulation and ablation patterns (the contour method) (Østrem and Brugman, Reference Østrem and Brugman1991).
-
(2) The geodetic method, which differences elevation models to determine changes in ice volume over time. Estimates of the density of the ice/snowpack are then used to convert this volume change to mass loss or gain. In the present day, these elevation models are often derived from airborne or satellite remote sensing.
Laser altimetry data from the ICESat satellite have provided a valuable resource to determine recent ice volume changes in the CAA (Gardner and others, Reference Gardner2011). However, such data do not work well for relatively small alpine glaciers. Glacier mass changes can also be measured using the gravimetric method over large ice caps and ice sheets (e.g. using data from the Gravity Recovery and Climate Experiment; Sharp and others, Reference Sharp2011), but these measurements lack sufficient spatial resolution to detect changes at the scale of alpine glacier basins.
Remote sensing methods are promising for future glacier mass-balance monitoring, but to detect and understand climatic trends in context we require earlier datasets for comparison. Surface mass-balance measurements using the glaciological method began in 1959 across the Queen Elizabeth Islands (QEI) of the northern CAA and primarily focused on two smaller ice caps (Melville and Meighen), two outlet basins of larger ice caps (Devon and Agassiz) and two mountain glaciers (White and Baby, Axel Heiberg Island). Intermittent observations were also conducted for many years on Ward Hunt Ice Rise and Ice Shelf on northern Ellesmere Island and Prince of Wales Icefield on SE Ellesmere Island (Hattersley-Smith and Serson, Reference Hattersley-Smith and Serson1970; Braun and others, Reference Braun, Hardy, Bradley and Sahanatien2004; Koerner, Reference Koerner2005; Mair and others, Reference Mair2009). Those records, submitted to the World Glacier Monitoring Service (WGMS) and exceeding 10 year in length, are shown in Figure 1a. Due to the large area of ice in the QEI (~104 000 km2; Sharp and others, Reference Sharp, Kargel, Leonard, Bishop, Kääb and Raup2014) and logistical barriers, the CAA is comparatively under-sampled in contrast to some other Arctic regions (e.g. Iceland, Svalbard, northern Scandinavia; Sharp and others, Reference Sharp2011). As such, it is particularly important to periodically assess the accuracy and correctness of these key measurements, which have served as in situ validation for numerous other studies, through periodic reanalysis of the glaciological mass-balance datasets (Zemp and others, Reference Zemp2013).
Fig. 1. (a) Locations of glaciers in the CAA with published mass-balance records exceeding 10 years in duration. ‘Obs.’ refers to the number of years of observations. (b) SPOT 5 composite image of White Glacier, 29 August 2008. Mass-balance stakes measured in 2014 are indicated in red and the white box notes the location of the nunatak time-lapse camera.
The comparison of mass-balance data series measured by the glaciological method with independently derived measurements of mass change offers insight into the completeness and accuracy of the measurements (Østrem and Haakensen, Reference Østrem and Haakensen1999; Cogley, Reference Cogley2009; Zemp and others, Reference Zemp2013). While the glaciological method excels in capturing the large temporal and spatial variability of climate/glacier processes, the method is also subject to biases that might be small on an annual basis but can compound to cause significant errors over multi-decadal records. Whereas random errors rise proportionally to the square root of the number of years, systematic biases sum linearly from year to year. The biases are associated with the inherent limitations of the glaciological method, primarily the challenge of measuring superimposed ice and internal accumulation, and are often undetectable within the random error margin that is approximated to be 200 mm w.e. a−1 for most glaciological balance measurements (Cogley and Adams, Reference Cogley and Adams1998). As the biases build over many years, however, they can be detected by multi-year geodetic balance observations. For example, the Norwegian maritime glacier, Ålfotbreen, was observed to have a positive glaciological balance of +3.4 m w.e. over a 20 year period, while the geodetic balance determined from map comparison was found to be –5.8 m w.e. (Østrem and Haakensen, Reference Østrem and Haakensen1999). The discrepancy was primarily assigned to uncertainties in the glaciological method associated with high accumulation rates (sinking stakes), and overestimation of accumulation by snow-probe measurements. As such, a reanalysis of mass-balance data should not only reveal how well a glaciological measurement programme is performing, but may also indicate the magnitude and extent to which specific processes have impacted the cumulative data series and the corrections that can be imposed to reconcile the geodetic and glaciological mass-balance measurements (Zemp and others, Reference Zemp2013; Andreassen and others, Reference Andreassen, Elvehøy, Kjøllmoen and Engeset2015).
This study presents the first reanalysis of a surface mass-balance record in the CAA for White Glacier, Axel Heiberg Island. The cumulative record of glaciological mass-balance observations is compared with a geodetic balance calculated from elevation data spanning 1960 to 2014 using an approach that strives to ensure consistency between the datasets before comparison (Zemp and others, Reference Zemp2013).
2. STUDY SITE AND PREVIOUS RESEARCH
White Glacier (79.50°N, 90.84°W) is a 14 km long alpine glacier extending from 80 to 1782 m a.s.l. in the region of Expedition Fiord on Axel Heiberg Island, Nunavut, Canada (Fig. 1a). The glacier has a 5 km wide accumulation area and flows southeast into a narrow 0.8–1.1 km wide valley, terminating at a junction with Thompson Glacier to the east (Fig. 1b). The glacier area in 1960 was 41.07 km2, which decreased to 38.54 km2 by 2014 (Thomson and Copland, Reference Thomson and Copland2016). The region experiences mean annual temperatures of ~–20°C and annual precipitation ranging from 58 mm a−1 at sea level (as measured at Eureka, 100 km to the east) to 370 mm a−1 at 2120 m a.s.l. as measured in a 41 year snowpit record of annual accumulation on the Müller Ice Cap (Cogley and others, Reference Cogley, Adams, Ecclestone, Jung-Rothenhäusler and Ommanney1996). Over the period of observation (1960–2015), the average equilibrium line altitude (ELA) was 1075 m a.s.l. and the mean accumulation area ratio (AAR, accumulation area divided by the total area) was 0.55.
The mass-balance programme at White Glacier was initiated in 1959 and brought about several years of intensive glaciological studies, all of which were based at the nearby McGill Arctic Research Station (79.42°N, 90.75°W). Studies that included snow accumulation processes (Andrews, Reference Andrews1964; Havens and others, Reference Havens, Müller and Wilmot1965), glacier facies (Müller, Reference Müller1962), ice velocity fluctuations (Iken, Reference Iken1974) and polythermal regime (Blatter, Reference Blatter1987a) were led by Dr. Fritz Müller of McGill University from 1960 to 1980. Dr. Müller, for whom Müller Ice Cap on central Axel Heiberg Island was named, passed away in July 1980 resulting in a 3 year gap in the mass-balance record until Trent University re-established the mass-balance monitoring programme in 1983. Since this time, studies at White Glacier have focused on mass-balance and remote sensing techniques (Cogley and others, Reference Cogley, Adams, Ecclestone, Jung-Rothenhäusler and Ommanney1996, Reference Cogley, Adams and Ecclestone2011a; Cogley and Adams, Reference Cogley and Adams2000). In 2014, University of Ottawa undertook responsibility for the White Glacier mass-balance monitoring programme. As a result of these combined efforts, today White Glacier has the longest mass-balance record for an alpine glacier in the Canadian Arctic (55 years) and is one of 37 reference glaciers within the Global Terrestrial Network for Glaciers through the United Nations Framework Convention on Climate Change. These observations, submitted annually to WGMS (e.g. WGMS, 2015 and previous issues), are used with observations from 40 other glaciers to calculate a worldwide glacier mass-balance index (WGMS, 2008).
3. METHODS AND DATA
In the following sections we present the measurement methods used to calculate the White Glacier glaciological and geodetic mass balances and describe the approaches taken to homogenize the data. Data homogenization refers to corrections for “artefacts and biases that are not natural variations of the signal itself but originate from changes in instrumentation or changes in observational or analytical practice” (Cogley and others, Reference Cogley2011b). We discuss the systematic and random errors associated with these datasets and address the differences between the mass-balance processes captured by each of the methods. In this case, while the geodetic method is considered to measure surface, internal and basal mass balance (Zemp and others, Reference Zemp2013), the glaciological method only captures surface processes (i.e. accumulation by snowfall and ablation by surface melt and sublimation). We use estimates of the internal and basal mass balance in an attempt to correct for these generic differences, allowing for comparison of results obtained with the two methods. Finally, the homogenized and corrected balances over the period of record (PoR, 1960–2014) are compared and the discrepancy between them is assessed statistically to determine whether the glaciological series requires calibration.
3.1. Glaciological balance measurements and calculation
Owing to logistical and financial constraints, field measurements at White Glacier are conducted once per year in the spring (typically late April to early May) and as such the measurements do not allow for differentiation between winter and summer balances. The stratigraphic year generally starts and ends in September (at the end of the melt season) and is identified by the year in which the ablation season ends (Cogley and others, Reference Cogley, Adams, Ecclestone, Jung-Rothenhäusler and Ommanney1996); for example the 2014 balance year combines accumulation from winter 2013/14 and ablation from summer 2014. Winter accumulation of the current balance year is measured in snow pits at stakes located above the lower limit of the equilibrium line (~900 m a.s.l.) and the previous summer surface is commonly identifiable as a compact and slightly dirty layer overlain by well-developed depth hoar and cup-shaped ice crystals (Müller, Reference Müller1961). Stake readings below the equilibrium line (identifiable by glacier ice at the base of snow pits) indicate the summer ablation of the previous balance year by the change in height of the stake above the ice surface.
Annual surface mass balance (B) is calculated using a modified version of the profile method (Østrem and Brugman, Reference Østrem and Brugman1991; Cogley and others, Reference Cogley, Adams, Ecclestone, Jung-Rothenhäusler and Ommanney1996). A polynomial regression (order k = 3) through b, the annual stake balances, at measured elevations along the glacier centreline is used to describe the mass-balance gradient, which is then extrapolated across the glacier hypsometry in 25 m elevation bands where s(h,t) is the area of the elevation band with midpoint elevation h in year t, and t is 0 in 1960. The reasoning behind using a third-order polynomial relates to the desire to better represent the differing slope of the mass-balance gradient above and below the ELA, which becomes less steep at higher elevations. Statistically, the third-order polynomial fit returns a better R 2 value (0.93) than a linear (n = 1) fit where R 2 = 0.88, however the overall impact of using a third-order polynomial versus a linear regression is ultimately negligible in the scope of this reanalysis.
The elevation band balances b(h,t) are calculated from the regression equation at the elevation band midpoints and the annual specific balance (i.e. normalized by S(t), the glacier area) is calculated as follows:
$$B_{\rm t} = \; \displaystyle{{\sum s(h,t)b(h,t)} \over {S(t)}}$$
From 1960 to 2003, B t was calculated using the 1960 glacier hypsometry s(h,0), which was derived from the topographic map of Haumann and Honegger (Reference Haumann and Honegger1964). In 2003, a mapping campaign by Hember and others (Reference Hember, Cogley and Ecclestone2003) allowed for an update to the hypsometry below 400 m a.s.l. and thereafter calculations of B have used the 1960 hypsometry above 400 m, and the updated 2003 hypsometry from 400 m down to the terminus. In 2014 a new DEM of White Glacier, described below and in Thomson and Copland (Reference Thomson and Copland2016), was produced and will serve as the new hypsometry for mass-balance calculations until the next mapping campaign and reanalysis. Records of the original and updated mass-balance series for White Glacier, including matrices of the hypsometries used in the calculations, are available in the supplementary information document.
A detailed history of mass-balance measurements at White Glacier is provided by Cogley and others (Reference Cogley, Adams, Ecclestone, Jung-Rothenhäusler and Ommanney1996), along with a discussion of the mass-balance record from 1960 to 1991 and an assessment of errors. The two most negative pentadal averages of mass balance occurred within the last decade (2005–14) and this has been accompanied by increasing ELAs and decreasing accumulation area ratios (Fig. 2a–c). The stake network at White Glacier has evolved from a dense network of >100 stakes in the 1960s, to 30–40 stakes in recent years (Figs 1b, 2d), resulting in an approximate coverage of 1 stake km−2. The 39 stakes measured in the last reported balance year (2014) ranged from 122 to 1519 m a.s.l. (~97% of the area of the glacier), and stake balances (b) from –2884 to +138 mm w.e.
Fig. 2. (a) The 55 year record of glaciological mass balance (B) with pentadal averages and standard errors; (b) ELA; (c) AAR; and (d) number of mass-balance stakes in network.
3.2. Conventional balances
An advantage of calculating glacier-wide mass balance over a historic reference surface is that it results in mass-balance values that exclusively reflect changes in glacier/climate interactions (Elsberg and others, Reference Elsberg, Harrison, Echelmeyer and Krimmel2001; Cogley and others, Reference Cogley2011b). However, this reference-surface balance does not reflect the actual mass change in a given year but rather the mass change that would have resulted from using a constant reference topography for which the glacier hypsometry (i.e. the s(h,t) and S(t) terms) remained unchanged through time. Therefore, a cumulative glaciological balance data series that would be equivalent to the geodetic balance would use conventional balances (Elsberg and others, Reference Elsberg, Harrison, Echelmeyer and Krimmel2001) that are calculated from glacier hypsometries that are concurrent in time with the mass-balance measurements. However, annual measurements of glacier hypsometry are impractical and nonexistent for most glaciers.
Conventional balances were calculated for the White Glacier mass-balance series using a linear regression that updates the glacier hypsometry annually between 1960 and 2014, the two years with complete topographic coverage based on the 1 : 10 000 map by Haumann and Honegger (Reference Haumann and Honegger1964) and the new 2014 DEM of White Glacier (Thomson and Copland, Reference Thomson and Copland2016). The time series of 25 m elevation band areas was calculated from the linear relationship,
$$s(h,t) = s(h,0) + t\; \times \left( {\displaystyle{{s(h,n) - s(h,0)} \over n}} \right)$$
where s(h,0) is the band area in the first year (1960) and t increases from 0 in unit steps over the n years in the record. These interpolated s(h,t) values were used to calculate the conventional balance using Eqn (1). The reference and conventional balances, with associated equilibrium line elevations and AAR, are provided in Table 1.
Table 1. Surface mass balance values and attributes calculated using the reference (ref) and conventional (con) hypsometry matrices, which are recommended as an update to the previous WGMS published versions (see supplementary information document)
Mass balance (B) is provided in mm w.e., ELA is the equilibrium line altitude (m) and AAR refers to the accumulation area ratio.
3.3. Geodetic balance
A 54 year geodetic balance of White Glacier was calculated by differencing the earliest large-scale map of White Glacier (1 : 10 000) based on 1960 aerial photography (Haumann and Honegger, Reference Haumann and Honegger1964), with a new DEM produced using Structure from Motion techniques from >400 oblique air photographs collected in July 2014 (Thomson and Copland, Reference Thomson and Copland2016). The 1960 White Glacier 1 : 10 000 map is available through the WGMS ‘Fluctuations of Glaciers Maps’ (http://wgms.ch/products_fog_maps/), along with supporting documentation that summarizes the field surveys, photogrammetry and plotting techniques used to produce it (Haumann, Reference Haumann and Muller1961, Reference Haumann1963; Blachut, Reference Blachut and Müller1961, Reference Blachut1963; Cogley and Jung-Rothenhäusler, Reference Cogley and Jung-Rothenhäusler2002). The 10 m contours of the 1960 map were manually digitized at a working scale of 1 : 2000 and converted into a 5 m × 5 m DEM using the Esri ArcGIS 10.1 ‘Topo to Raster’ tool, which is based on the ANUDEM software developed by Hutchinson (Reference Hutchinson1988). The 2014 DEM was originally produced at 1 m resolution in NAD27. The surface was then down-sampled to 5 m resolution to allow direct comparison with the 1960 DEM and reprojected into a local planar coordinate system used to define the original 1960 map. A coregistration of the two DEMs was finally performed to minimize the error in the observed volume change. To determine the glacier-wide average change in height, the geodetic balance was calculated as the sum of observed changes in elevation, pixel by pixel, within the 1960 White Glacier extent, divided by the mean area between 1960 and 2014 to arrive at the glacier-wide average change in height (Zemp and others, Reference Zemp2013). Specific details of the 2014 DEM production, projection systems, the application of ground control, coregistration techniques and error analysis can be found in Thomson and Copland (Reference Thomson and Copland2016).
3.4. DEM coregistration
Coregistration, guided by minimization of the elevation differences between DEMs over areas of stable terrain near the glacier, can be considered as a means of homogenizing the geodetic datasets (Zemp and others, Reference Zemp2013). The universal co-registration method (Nuth and Kääb, Reference Nuth and Kääb2011) exploits the relationship between elevation differences (dh), terrain slope and aspect, and was applied here to identify the translational offsets between the geodetic DEMs. The full details of the coregistration process, including parameters and error analysis are provided in Thomson and Copland (Reference Thomson and Copland2016). This was assessed over three regions of stable terrain, together covering 4.5 km2 and with a slope range of 5°–30° and elevation range of 50–1410 m a.s.l (Fig. 3a). Introducing a small rotation correction before using the Nuth and Kääb (Reference Nuth and Kääb2011) method significantly improved the co-registration results (Table 2). The rotation correction, about the z-axis, was determined using Horn's absolute orientation solution (Horn, Reference Horn1987), which produces a least-squares fit between the coordinates of common points between the 2014 DEM (source points) and 1960 DEM (target points).
Fig. 3. (a) Average change in elevation (dh) between 1960 and 2014 observed in each 25 m hypsometric band (1960 hypsometry and extent shown). The upper right legend refers to elevation changes over ice, while the lower left legend indicates elevation differences over stable terrain following the DEM coregistration. The position of the Thompson Glacier terminus as of 2014 is illustrated in the lower right. (b) White Glacier hypsometry (1960; grey bars) and observed elevation changes (blue curve) over PoR for 100 m elevation bands.
Table 2. Coregistration parameters for homogenization of 1960 and 2014 DEM areas over stable terrain using rotation and translation, where ccw indicates counter clockwise rotation
Final parameters applied to the full DEMs are in bold.
Rotation and translation parameters that resulted in the minimum average dh and standard deviation σ dh over stable terrain were applied to the 2014 DEM, from which the 1960 DEM was then subtracted to produce a dh map (Fig. 3a). Small voids in the 2014 DEM, consisting of <3% of the basin area, were treated by assigning the void space a dh value equal to the average dh of its 25 m elevation band. For example, if a void were judged to lie within the 1000–1025 m a.s.l. elevation band, it would be assigned the average of all dh values found within that band. The same method was applied in the area where Thompson Glacier has advanced since 1960 into the area previously occupied by White Glacier (Fig. 3a). The distribution of dh across the glacier hypsometry is shown in Figure 3b.
4. RANDOM AND SYSTEMATIC ERRORS
Uncertainties of the glaciological method are associated with point measurements in the field and the extrapolation of these measurements across elevation bands (Cogley and Adams, Reference Cogley and Adams1998). Regarding the latter, the longitudinal design of the stake network limits the extent to which we can observe and assess variability of mass balance in the glacier basin. For the 1960–62 balance years, Young (Reference Young1972) conducted a parameterization of surface mass balance as a function of elevation, slope, aspect and curvature across the White Glacier basin and found that the parameterization resulted in slightly higher estimates of melt, which differed from an elevation-only parameterization by –0.10 ± 0.04 m w.e. a−1. At the same time, studies of spatial variability in the ablation area at three cross-profiles ranging from 250 to 850 m a.s.l. showed a standard deviation in stake balances of 0.196 m w.e. a−1 (Adams, Reference Adams1966). The variability of the stake balances with elevation (Fig. 4) has not changed measurably over the PoR and the polynomial regression (k = 3) used to estimate band balances each year has an average R 2 = 0.93. As has been predicted (e.g. Dyurgerov and Dwyer, Reference Dyurgerov and Dwyer2000) we do observe a slight (though not statistically significant) steepening of the linear mass-balance gradient towards the present day, as well as a shift towards more negative mass-balance values, which corresponds with increasing ELAs over the PoR (Fig. 4b). From the above studies and previous literature concerning the errors associated with the glaciological method (Adams, Reference Adams1966; Cogley and others, Reference Cogley, Adams, Ecclestone, Jung-Rothenhäusler and Ommanney1996), we take 0.250 m w.e. a−1 to be the random error associated with the glaciological record (σ glac.a), which incorporates 0.20 m w.e. a−1 as the random sampling error and an estimated 0.050 m w.e. a−1 to account for unquantified systematic errors.
Fig. 4. (a) Stake mass balances for 1960–2014 tending from blue to orange colours over the PoR. The overlying black line indicates the average linear mass-balance gradient. (b) Progression of the decadal mass-balance gradients (m w.e. km−1) from 1960 to 2014, ending with the 2010–14 pentadal gradient.
For the geodetic balance, the coregistration results over stable terrain are taken to define errors associated with the DEM comparison. The bias in the geodetic balance over the PoR, εgeod.PoR = 0.663 ± 0.047 m w.e., is derived from the average dh over stable terrain (0.78 m) and an assumed density of 850 ± 60 kg m−3, as described below. The random error in the geodetic balance, σ geod.PoR = 0.866 m w.e., is taken from the DEM standard error and the uncertainty in the density conversion. The DEM standard error is calculated to be twice the standard deviation of dh over stable terrain (2 × 5.08 m) divided by the root of the sample size. Given the spatial autocorrelation between pixels of a DEM derived from contours (as is the case for the 1960 DEM), we approximate the sample size to be 319, based on the number of individual contours over glacier ice in the 1960 map. This is likely a conservative estimate for the degrees of freedom considering that it would also be reasonable to use the number of photographs used to produce the 2014 DEM (507) as a proxy for the independent sample size, which would reduce σ geod.PoR to 0.814 m w.e.
In both the glaciological and geodetic balance calculations, we consider random uncertainties to have a normal error distribution and to be uncorrelated when combined (i.e. in the calculation of geodetic balance errors).
5. CORRECTIONS OF GENERIC DIFFERENCES
5.1. Density conversion
One of the generic differences between glaciological and geodetic mass-balance measurements is that the former are derived from mass change measurements that are derived from well constrained in situ density measurements, whereas the latter are seated in volume change measurements, which require conversion to mass by an estimated density for the material lost or gained. While density conversion is a significant issue (e.g. Moholdt and others, Reference Moholdt, Nuth, Hagen and Kohler2010), this study follows the recommended approximation for density conversion of 850 ± 60 kg m−3 developed by Huss (Reference Huss2013) which, through modelling of synthetic and monitored glaciers, incorporates elements of densification, internal refreezing and changing volume of the firn. The approximation is suitable for White Glacier, which meets the requirements of Huss (Reference Huss2013) that the geodetic observations span >5 years, and that the glacier demonstrates relatively stable mass-balance gradients, non-negligible changes in volume and an established and stable region of firn. It also enables consistent comparison with other mass-balance reanalysis studies (Zemp and others, Reference Zemp2013).
5.2. Survey differences
Temporal differences between the geodetic and glaciological observations need to be addressed because the geodetic observations span the air photo survey dates of 2 August 1960–10 July 2014, rather than the end-of-summer dates in each of these years. To align the geodetic dates with the stratigraphic year used for the glaciological mass-balance measurements, a simple degree-day model (Eqn (3) of Zemp and others, Reference Zemp2010) was used to estimate melt between the observation dates in 1960 and 2014 and the ends of their respective stratigraphic years, which were taken to be 1 September 1960 and 8 September 2014, the dates on which temperatures near the terminus dropped below 0°C. Temperature data for August and September 1960 were acquired from Andrews (Reference Andrews1964) at ‘Lower Ice Station’, and from a Campbell Scientific automatic weather station (AWS) at the White Glacier end moraine (500 m southwest of the terminus) for the same months in 2014. Following Braithwaite and Raper (Reference Braithwaite and Raper2007) we assumed a degree-day factor for ice of 7.0 mm w.e. d−1 K−1 (their ‘Medium’ assumption) to estimate mass-balance corrections of +0.095 m w.e. for 1960 and –0.197 m w.e. for 2014 that were applied as corrections to the geodetic balance, as shown in Figure 5.
Fig. 5. Uncorrected and corrected balances from (a) the glaciological method and (b) the geodetic method. A breakdown of the corrections applied to each method to homogenize the balances is shown as inset graphs.
The gap in the mass-balance record (1980–82) was filled by Glenday (Reference Glenday1989) using a multivariate regression of glaciological balance values from 1970–80 with temperature and shortwave radiation data available from the nearest weather station, at Eureka, Nunavut, ~100 km east of White Glacier. The parameterization returned balances of –0.178, –0.175 and –0.092 m w.e. for 1980–82. Unfortunately, information on the uncertainty of this parameterization is not available, so we estimate the uncertainty to be the standard deviation of the 51 years with observations (0.320 m w.e. a−1).
5.3. Internal and basal mass balance
The glaciological method measures surface mass balance, and is unable to capture internal and basal accumulation and ablation. Here, simple estimates of these processes are considered and applied to the glaciological balances to allow comparison with the geodetic balance. In temperate glaciers, potential energy loss by the transfer of mass downslope leads to internal strain heating and frictional warming at the bed. Considering that White Glacier comprises cold ice (–15 to –5°C for the majority of the ice volume; Blatter, Reference Blatter1987a), we assume that any strain or frictional energy will go towards only slightly warming this cold ice, rather than promoting melt. Internal ablation due to the loss of potential energy (A PE ) of englacial and subglacial runoff was estimated using
$$A_{PE} = \; \displaystyle{{mhg} \over {L_{\rm f} \rho S}}$$
where m is the mass of flowing water, g is acceleration due to gravity, h is the elevation loss, L f is the latent heat of fusion, ρ is the density of water, and S is the glacier area (average between 1960 and 2014). The cumulative measured surface ablation over the 54 year record was used to estimate the mean annual discharge, and subsequently m, to be 1.4 × 1010 kg a−1. The elevation loss was approximated by assuming that on average, meltwater originates at half the vertical distance from the terminus to the average ELA, h = (1070–50)/2 = 510 m. Together, these resulted in a calculation of A PE = 0.274 m w.e. over the PoR (i.e. 5.1 mm w.e. a−1). We have assumed no melting by turbulent transfer of energy from the flowing water to the ice, but have also assumed that all the meltwater loses all its potential energy to internal ablation and none to warming of the bedrock. With these assumptions we are neglecting biases of the same order as our estimate of A PE , and they will tend to cancel.
Internal accumulation by surface melt and refreezing in the firn does not lead to a bias in glaciological balance measurements at White Glacier because the snow accumulation measurements, made in the spring, are presumed to be representative of winter accumulation and therefore are likely not to have been affected by melting. However, summer accumulation is not necessarily negligible and may constitute a bias such that reported annual balances are too negative. At other monitored glaciers in the CAA, summer accumulation (excluding that by rainfall that freezes) gained between April–May mass-balance measurements and the end of the melt season has been determined by the use of records from AWS (Koerner, Reference Koerner2005). Currently the only snow-depth sounder installed at White Glacier is situated in the ablation area to monitor the onset and duration of surface melt. Instead, photography from a time-lapse camera located on a nunatak of White Glacier (white square in Fig. 1b) with a field of view covering elevations between ~900 and 1400 m a.s.l. was used to detect the occurrence of snowfall in the accumulation area between April 2013 and July 2015. These observations indicate that between 1 and 3 snowfall events have occurred each year in the period following spring mass-balance measurements and before the melt season begins. Together, these events amount to <20 cm of snow accumulation in each year on average (<22 mm w.e. using ρfresh snow = 108 kg m−3; Metcalfe and others, Reference Metcalfe, Ishida and Goodison1994). This contribution to the total mass balance of the glacier is considered negligible because it comprises the first layers of snow that are exposed at the onset of summer melt and therefore the first mass to melt and run off. During the melt season only one snowfall event per year was observed by the time-lapse camera in 2013 and 2014, and these events resulted in snowfall accumulations of 5 cm (<6 mm w.e.), which rapidly melted at lower elevations and likely upper elevations as well (given observations of ‘glazing’ of the snowpack surface within a day or two of snowfall).
Basal ablation from geothermal heating was estimated for the area of White Glacier where the bed is known to be temperate along the centreline of the glacier trunk, which totals ~5 km2, or 12.5% of the total glacier area (Blatter, Reference Blatter1987a). Basal ablation by geothermal melt (A GT ) was estimated to be 0.015 m w.e. over the PoR when calculated using
$$A_{GT} = \displaystyle{{qt} \over {L_{\rm f} \rho}} $$
where the geothermal heat flux q is taken to be 0.023 W m−2 (Blatter, Reference Blatter1987b), and t is the time elapsed (s) over the PoR. Basal accumulation can be considered negligible for non-tidewater glaciers (AMAP, 2011).
6. RESULTS AND DISCUSSION
In the glaciological record, a +0.445 m w.e. correction was applied to convert the cumulative reference-surface balance (–11.214 m w.e.) to a conventional balance (–10.769 m w.e.) using a linear interpolation of changing hypsometry over the 51 years with observations. Figure 6 illustrates the difference in balances calculated from the reference surface (1960) and conventional surfaces (interpolated). The conversion from reference to conventional balances incorporates two scenarios that present contrasting effects: (1) as the glacier terminus retreats, the glacier loses area in the region of most strongly negative balance resulting in a net positive influence on the glacier-average balance, and (2) as the glacier thins, its surface lowers, which exerts a net negative influence on the mass balance because ablation is stronger at lower altitudes. The positive correction applied to arrive at our reference-surface balance would suggest that the former process dominates over the latter. Investigating the relative roles of these processes requires calculating and comparing mass change related to a fixed elevation with varying area (B HdS for scenario 1), and the mass change resulting from elevation change (generally ice thinning) assuming a constant area (B SdH for scenario 2). For this approach it is necessary to work with total w.e. mass change rather than specific balances to remove the second area term (denominator in Eqn (1)). We find that B ref = –486 321 kt and B con = –447 206 kt, requiring that a correction C of +39 115 kt be applied to B ref to arrive at B con . This correction is the sum of the positive B HdS and the negative B SdH . From the geodetic balance calculations, we calculate B SdH using
$$B_{SdH} = \rho \sum \left( {\overline {dh} \; s(h_0 )} \right)$$
where dh is the average elevation change observed over the PoR across a given elevation band, s is the band area according to 1960 hypsometry, and ρ is the assumed density (850 ± 60 kg m−3). Consequently, if C = B SdH + BHdS we find that with B SdH = –386 423 kt from Eqn (5), B HdS therefore equals 425 538 kt. The greater absolute magnitude of B SdH therefore indicates that glacier retreat has a greater impact than glacier thinning on the mass balance of White Glacier.
Fig. 6. (a) Difference (dB) between the reference (measured) and conventional (interpolated assuming a linear change in hypsometry) glaciological balances from 1960 to 2014. (b) Cumulative balances over the 54 year period with the cumulative difference (dB) plotted above.
The remaining corrections to the glaciological record, including estimates of mass balance during the gap years (1980–82) and melt induced by the loss of potential energy and by geothermal heating, together lead to a more negative corrected glaciological balance by –0.445, –0.274 and –0.015 m w.e., respectively, over the PoR (Fig. 5a). The final estimate of the glaciological cumulative balance is –11.503 ± 2.037 m w.e. between 1960 and 2014.
For the geodetic record, the average dh across the 1960 hypsometry shows that glacier thinning (maximum –30 m) was prominent from the terminus up to elevations of ~1400 m, and that moderate thickening (maximum 15 m) was observed at the highest elevations (Fig. 3). The average elevation change across the mean glacier extent (39.81 km2) was –11.97 m over the 1960–2014 PoR. When the elevation change is multiplied by density, the resulting initial geodetic estimate of mass balance is –10.17 ± 0.866 m w.e., with an estimated error of 0.72 m w.e. Corrections to the geodetic balance (Fig. 5b) include accounting for the observed non-zero elevation change over stable terrain after coregistration, which increased the geodetic balance by 0.663 ± 0.047 m w.e. over the PoR. The geodetic balance also required small adjustments of +0.095 and –0.197 m w.e. to account for melt between the geodetic survey dates and the end of the respective 1960 and 2014 stratigraphic years. We cannot estimate the uncertainty of these modelled adjustments but we assume that they are small.
The random errors (σ) in the two methods are significantly different due to the way in which they propagate over time. Over n years,
$$\sigma _{glac.a} = \displaystyle{{\sigma _{glac.PoR}} \over {\sqrt n}} $$
and
$$\sigma _{geod.a} = \displaystyle{{\sigma _{geod.PoR}} \over n}$$
such that σ glac.a = 0.277 m w.e. a−1 (2.037 m w.e. over PoR) determined from the combined 0.250 m w.e. a−1 error for the 51 years with observations in the glaciological record and 0.320 m w.e. a−1 for the parameterized gap years. For the geodetic balance, σ geod.a = 0.016 m w.e. a−1 (0.866 m w.e. over PoR), which incorporates errors associated with the density conversion (0.718 m w.e. over PoR) and the standard error in the elevation change observation (0.484 m w.e. over PoR). In this study the geodetic balance is therefore taken to be the standard against which the glaciological balance is compared, owing to its lower random error.
Figure 5 illustrates the original glaciological and geodetic balances and the corrections applied to homogenize the data. The resulting glaciological and geodetic balances over the PoR, determined over the average of the 1960 and 2014 glacier extents, are –11.50 m w.e. and –9.61 m w.e., respectively. Note that the corrections of biases did not reduce the ~1 m w.e. difference between the uncorrected balances, but rather they almost doubled the difference, from 1.04 to 1.89 m w.e.
The statistical significance of the difference between balances is determined from the reduced discrepancy:
$$\delta = \; \; \displaystyle{{\Delta _{PoR}} \over {\sigma _{common\; PoR} \;}} $$
where ∆ PoR is the difference between the corrected glaciological and geodetic balances (–1.89 m w.e.) and σ common.PoR is the common variance (2.21 m w.e.), which describes the combined random errors of the glaciological and geodetic balances (Zemp and others, Reference Zemp2013). As the reduced discrepancy (δ = –0.85) falls within the 95% confidence interval (–1.96 < δ < 1.96), the balances can be deemed statistically equivalent. Therefore, the glaciological method can be considered to have accurately captured the long-term mass changes at White Glacier when compared against the geodetic balance. However, it is not possible to statistically identify any biases that might explain the observed –1.89 m w.e. discrepancy over the PoR. From the common variance it is also possible to calculate the smallest bias that could theoretically be detected in the glaciological record within this 54 year period (Zemp and others, Reference Zemp2013); the bias calculated at the 5% risk limit, which can be of either sign, is 0.148 m w.e. a−1 and well below the random error of 0.250 m w.e. a−1 in the glaciological balance measurements.
In determining the geodetic mass balance for this study, DEM coregistration combined with rotation was an important step considering the challenges of working with the 1960 map, for which there is limited information about the map projection and datum (Cogley and Jung-Rothenhäusler, Reference Cogley and Jung-Rothenhäusler2004). The coregistration method of Nuth and Kääb (Reference Nuth and Kääb2011) enabled alignment of the DEMs independent of georeferencing difficulties and reduced the random error and biases between the DEMs. In fact, DEM alignment without coregistration would have resulted in a positive geodetic balance of 1.63 m w.e. over the PoR. Coregistration by translation alone was found to decrease the random error by 14% and, when combined with rotation, to reduce it by 39%.
While the observed discrepancy (–1.89 m w.e.) between the glaciological and geodetic balances is not statistically significant, it is possible to explore the parameter space in which certain parameters vary. Of particular note is the uncertainty associated with density applied in the calculation of the geodetic balance (Moholdt and others, Reference Moholdt, Nuth, Hagen and Kohler2010; Huss, Reference Huss2013). Holding all other parameters fixed, we find that allowing the density to vary within 850 ± 60 kg m−3 results in a maximum possible discrepancy of –2.27 m w.e. over the PoR when density equals 790 kg m−3. The resulting delta value of this discrepancy is –1.11 m w.e., which still falls within the 95% confidence interval (–1.96 < δ < 1.96). A similar test was conducted by holding all parameters fixed and allowing the modelled mass balance for each of the three gap years to fluctuate within the observed variance of the mass-balance record. The resulting possible glaciological balances ranged from –10.86 to –11.57 m w.e., the latter of which produced a discrepancy of –2.34 m w.e. and a delta value also within the 95% confidence interval. If both density and the modelled balances over the gap years are combined, the maximum discrepancy is –2.34 m w.e. with a delta value of –1.29, which still remains within –1.96 < δ < 1.96.
It is possible that the more negative balance observed from the glaciological record is attributable to the assumption that summer snowfall is negligible; however the lack of observational data limits our ability to estimate the impact of this assumption. Similarly, while internal refreezing below the summer snowpack has been shown to impact glaciological (i.e. summer mass balance) and geodetic measurements in other studies (e.g. Bolch and others, Reference Bolch2013; Huss, Reference Huss2013) we are unable to treat it quantitatively. Another possible explanation for the more negative glaciological balance is the missed contribution of mass from superimposed ice originating from snowmelt off nunataks. In theory superimposed ice, recognized as such, does not impact mass-balance measurements because it is a redistribution of mass within the system. However, snow originating on nunataks, which is not included in the mass-balance calculation, regularly melts and runs off during the summer months and could refreeze in the accumulation area if the nunatak is above the ELA, or may freeze as superimposed ice. At present, nunataks comprise 6% of the White Glacier basin above 1000 m a.s.l. (the approximate average ELA) and including their area in the glacier hypsometry and subsequent mass-balance calculations would only lead to a small increase (<1%) in the calculated balance. However, as nunataks expand with increasing melt, as we observed by an increase in the overall nunatak area by 2% from comparison of the White Glacier 1 : 10 000 maps by Haumann and Honegger (Reference Haumann and Honegger1964) and Thomson and Copland (Reference Thomson and Copland2016), this proposed effect may impact mass-balance measurements in future years.
7. CONCLUSION
Reanalysis of the White Glacier mass-balance record from 1960 to 2014 indicates that, in comparison with the observed geodetic balance, the glaciological method successfully captures the long-term mass change and does not require correction or calibration. This finding suggests that the assumptions made in the glaciological method, including that summer accumulation is negligible at White Glacier, are reasonable for the PoR. One of the more significant findings of this study was the importance of employing the coregistration techniques of Nuth and Kääb (Reference Nuth and Kääb2011), and in this case DEM rotation as well, before calculating the geodetic balance. Indeed, failure to do so would have resulted in a positive geodetic balance over the 54 year PoR, as opposed to the measured −9.61 m w.e. loss. Homogenization of the glaciological record by conversion to conventional balances reduced the discrepancy with the geodetic balance by ~0.5 m w.e. More detailed analysis of the individual roles of glacier retreat and thinning suggests that glacier thinning impacts calculated mass balance to a lesser extent than glacier retreat, such that loss of significantly negative mass-balance regions presents a buffering effect under the negative mass-balance conditions. Corrections applied to address processes that define the generic difference between the glaciological and geodetic balances (i.e. internal and basal ablation) were found to increase the difference. This may suggest that the surface mass balance is a reasonable approximation to the total (climatic–basal) mass balance for land-terminating glaciers in the CAA, which is perhaps not surprising considering the short summer season at high latitudes during which internal and basal melt might occur.
SUPPLEMENTARY MATERIAL
To view supplementary material for this article, please visit https://doi.org/10.1017/jog.2016.112
ACKNOWLEDGEMENTS
This research was supported by the Natural Sciences and Engineering Research Council of Canada, Canada Foundation for Innovation, Ontario Research Fund, ArcticNet, University of Ottawa, the W. Garfield Weston Foundation and the Canadian Polar Commission's Northern Studies Training Program. The authors would like to acknowledge the continued support of the Canadian Polar Continental Shelf Program for providing field logistics at White Glacier. Particular thanks go to Wayne Pollard, current Director of the McGill Arctic Research Station, where the glacier research programme has operated for 55 years. The authors are grateful for thoughtful reviews by E. Thibert and an anonymous reviewer.
Open access
