Glaciological method

Snow pits and probing are used to directly measure accumulation, whereas stakes measure surface ablation (e.g. Østremand Brugman, 1991; Reference Kaser, Fountain and JanssonKaser and others, 2003). We use the stratigraphic system to define the annual mass balance, whereby measurements are made between successive annual minima, typically in early September at Castle Creek Glacier. Conversion to w.e. is made by assuming the density of ice to be 900 kg m^_3, and through direct measurement of snow density from snow pits (2−4). Measured snow density at the end of the ablation season is less spatially variable than snow depth, with standard deviations of <1% and 33% respectively. We measured annual mass balance at a point (b _a) at 21, 12 and 18 sites during the 2009, 2010 and 2011 balance years, equating to sampling densities of 2.2, 1.3 and 1.9 km^_2 respectively. Glaciological measurements are absent for a portion of the middle of the glacier where an icefall impedes safe travel.

We apply the balance-gradient (or regression) method to extrapolate from b _a measurements to glacier-wide annual mass balance (B _a), using the surface area within 50 m elevation bins (e.g. Reference Fountain and VecchiaFountain and Vecchia, 1999). For the glaciological method, we use a three-part linear spline to represent the variation of b _a with elevation. This spline is derived from the b _a measurements, with the intercept set to the observed elevation of the annual snowline (b _a = 0). Our measurements do not reach the highest elevation bins of Castle Creek Glacier, so we use measurements from our highest observations for these uppermost elevation zones of the glacier (e.g. Reference Cogley, Adams, Ecclestone, Jung-Rothenhäusler and OmmanneyCogley and others, 1996). Use of a linear spline to interpolate from mass-balance observations has been found to be similar to a quadratic interpolation and superior to the contour method (Reference Fountain and VecchiaFountain and Vecchia, 1999).

Photogrammetric method

We use aerial photographs taken in 2008, 2009 and 2011 to derive B _a for 2009, and cumulative balances for the periods 2008−11 and 2009−11 (Table 1). Ground sampling distance of the 2008 and 2009 images is 0.25, and 0.53 m for the 2011 images. Ground control points (GCPs) were obtained from stereo models of 2005 aerial triangulation scans made available by the Province of British Columbia.

Table 1. Aerial photography and stereo model details (GCP: ground control point)

We created stereo models from the three years of photography using the Vr Mapping photogrammetry software suite (Cardinal Systems LLC). A set of 18 common GCPs, consisting of bedrock features or stable boulders distributed around the glacier at various elevations, and 50−70 tie points were used for exterior orientation and generation of stereo models (e.g. Reference Schiefer and GilbertSchiefer and Gilbert, 2007; Reference Barrand, Murray, James, Barr and MillsBarrand and others, 2009). The use of the same GCPs for all years ensured that positional errors were randomly distributed (Ka_a_b and Vollmer, 2000; Reference Schiefer and GilbertSchiefer and Gilbert, 2007; Reference Schiefer and GilbertSchiefer and others, 2007). Analysis of 11 check patches allows us to measure the relative accuracy and assess systematic bias between stereo models. These check patches consist of 25 individual check points in a 5 m grid, located on stable bedrock near the glacier(Figs 1 and 2). We estimated systematic bias among stereo models and derived trend surfaces based on the mean residuals of the 11 check patches; these trend surfaces were then used to apply a correction for elevation points on the glacier. We also compared the average slope angle of each check patch with the mean residual of the 25 individual points to detect any horizontal bias among models.

Fig. 2. Elevation residuals of check patches collected from three different stereo models. Each box plot displays the distribution of surface-elevation residuals of 25 check points, defined as the difference in surface elevation of the same horizontal coordinates from the first stereo model to the second in a given epoch. The three box plots for each check patch show residuals for the epochs 2008−09, 2009−11 and 2008−11 respectively from left to right. See Figure 1 for check-patch locations.

To map glacier surface elevation, we manually digitized mass points (series of x,y,z data points collected on a predetermined grid) on a 100 m grid within the glacier extent as manual measurements yield better results than automated extraction methods (Reference McGlone, Mikhail and BethelMcGlone and others, 2004). Good contrast in the 2008 and 2009 photography enabled us tomeasure every gridpoint (n = 937), but poor photographic contrast due to fresh snow cover in the2011 photography reduced measurements by 26% (Fig. 3).

Fig. 3. (a) Grids of photogrammetric mass points for 2011 (black dots) and 2008 and 2009 (black dots and open circles); (b) subsets of mass points for _long profile_ (black dots) and _grid_ points (open triangles); and (c) _walkable routes_ (open circles) and _arrays_ (black triangles).

The largest source of error in our photogrammetric methodology is in the operator_s ability to perceive and measure the glacier surface, a methodological shortcoming noted by others (e.g. Reference Rasmussen and KrimmelRasmussen and Krimmel, 1999). This bias decreased, however, with repeated measurements (Fig. 4). We define blunders as points that fall outside 68% (±σ1?) of measured elevation change for a given 50 m elevation bin, and we remeasured blunders five times for the 2008 and 2009 models, and three times for the 2011 models. To correct persistent blunders (<5%), we performed ordinary kriging interpolation from the non-blunder mass points, and extracted elevations from these trend surfaces.

Fig. 4. Thickness change for balance year 2009 from manual aerial photogrammetry. The two sets of mass points are from the first (v1) and sixth (v6) measurements of mass points, illustrating improved performance of the photogrammetrist.

As dates of photography did not exactly coincide with our in situ measurements, we made corrections using a temperature-index model (Reference HockHock, 2003). We derived glacier surface temperature for the model from observations at an automated weather station near the east margin of the glacier (2105 m a.s.l.; Fig. 1), which is outside the influence of katabatic winds (Déry and others, 2010), and assume a lapse rate of 0.0068C m^_1. We apply melt factors (mm w.e.°C^_1 d^_1) for snow (k _s) and ice (k _i) derived from mass-balance data (Reference Shea, Moore and StahlShea and others, 2009) for Place Glacier, British Columbia (k _s = 2.76, k _i = 4.67), and compare our findings to results using melt factors from Peyto Glacier, Alberta, Canada (k _s = 2.34, k _i = 5.64), and our more limited observations at Castle Creek Glacier (k _s = 3.45, k _i = 4.33).

Geodetic mass-balance methods necessitate assumptions of the density of ice and snow lost or gained from the surface of a glacier (e.g. Reference HussHuss, 2013). We tested three scenarios to assess our density assumptions: The first scenario (A) assumes that the density profile of the glacier remains unchanged with time (Sorge?s law; Reference BaderBader, 1954), and 900 kg m^_3 is used to convert from ice eq. to w.e. Our second scenario (B) uses 900 kg m^_3 for all points below the equilibrium-line altitude (ELA), and 750 kg m^_3 for all points above the ELA, assuming that the loss or gain of material in the accumulation area is not entirely composed of ice, but at least partially of firn and snow (e.g. Reference ZempZemp and others, 2010). The third scenario (C) uses 900 kg m^_3 for all points below the ELA, and 600 kg m^_3 for all points above the ELA, assuming that the loss or gain of material in the accumulation area has a density equivalent to end-of-season densities measured in our snow pits. These three scenarios use maps of accumulation and ablation areas based on glacier extent and observations of the ELA in the latter year of each period. Densities for Castle Creek Glacier elevation change vary from 800 to 850 kg m^_3 for B, and from 699 to 726 kg m^_3 for C.

To achieve B _a via the photogrammetric method, we use two spatial extrapolation methods. For balance year 2009, when the 100 m grid is measured in its entirety, we apply the arithmetic mean of all 937 points. This assumes that all vertical velocities sum to zero (Eqn (2)). The second method sums the product of the average dh (z) from each 50 m bin and the bin_s corresponding surface area. In this second method, the integral of the emergence velocities from the subset of points might be zero (e.g. along a flowline; Reference CogleyCogley, 2005) or nonzero, introducing an error or bias in the estimate of B _a.

GPS method

To perform RTK GPS surveys (hereafter denoted as GPS) we used Topcon GB-1000 dual-frequency receivers (measurement precision of ±0.03 m). Reoccupation of previously measured points indicated a horizontal accuracy of ±0.03 m, and repeat measurements of three check points on stable bedrock respectively indicated horizontal and vertical accuracy of 0.02 ±0.01 and 0.03 ± 0.02 m. We attached the rover antenna to a short antenna pole that was in turn attached to the backpack of the operator. The distance from the rover antenna to the glacier surface was measured on a flat surface and assumed to be constant (e.g. Reference Nolan, Arendt, Rabus and HinzmanNolan and others, 2005).

We made measurements of dh via GPS along longitudinal profiles and at gridpoints (Figs 1 and 3). Forinitial measurement of the longitudinal profiles, we collected points in kinematic mode at 5 s intervals and differentially corrected these points. From these kinematic profiles we selected points every 10 m in elevation for GPS measurementin successive years. We established grids of points comprising cross-glacier profiles at common elevations, yielding 20 points in the ablation area and 16 points in the accumulation area (Fig. 1).

Measured dh was converted to w.e. using the three density assumption scenarios discussed above. GPS measurements were made at the same time as the glaciological measurements, so no melt correction was needed. Unfortunately, poor line of sight with the base station and insufficient base radio power resulted in few measurements in the accumulation area, so we adopt the balance-gradient method to achieve B _a using a linear spline and the period-specific hypsometry. For GPS measurements, we use a two-part linear spline fit to our 2011 observations (Fig. 5). This spline is then shifted to match 2009 and 2010 ablation-area b _a measurements. With this shift we assume the profile shape does not change from year to year, and that few lower-elevation observations can be used to adequately define b (z) (e.g. Reference Rasmussen and KrimmelRasmussen and Krimmel, 1999).

Fig. 5. Measurements of at-a-point mass balance and thickness change, and associated profiles with elevation, for balance years 2009, 2010 and 2011. Spatial extrapolation uses the hypsometry displayed at left. Error bars in the top panel indicate 1 of the photogrammetric measurements within each 50 m elevation bin.

Additionally, we tested the viability of using fewer, spatially limited GPS measurements to estimate B _a via four subsets of the 937 photogrammetric measurements for balance year 2009 (Fig. 3). These subsets include: (1) a ′long profile′ of 41 points along the glacier center line, (2) ′walkable routes′ consisting of 56 points, which are the safely navigable routes of our GPS longitudinal profiles, (3) ′arrays′ consisting of the 36 ablation and accumulation area points measured in situ, and (4) a ′grid′ consisting of 61 points, regardless of safe travel.

Estimation of vertical velocity

We employed an Eulerian frame of reference to estimate vertical velocity at fixed coordinates using Eqn (7) with GPS measurements of ? and glaciological measurements of ?. From August to September (2008−10), we estimated vertical velocity for point arrays in the ablation and accumulation areas of Castle Creek Glacier (Fig. 1). We placed 20 ablation stakes (ablation array) in four across-glacier profiles. At each stake, which typically travelled 5−20 m a^_1 down-glacier, we measured surfaceablation; dh was measured with GPS at fixed coordinates, where stakes were initially placed. We assume that ablation measured at a transient stake is representative of ablation at the site where it was initially placed.

To estimate vertical velocity in the accumulation area, we made four probing observations within 3 m of the location where surface elevation was measured with GPS. The average of these multiple observations was used to minimize errors stemming from a non-vertical probe, the observer_s ability to accurately probe the previous summer surface, and the effects of decimeter-scale variability in the summer surface (e.g. from suncups, meltwater channels and differential ablation). We assume that the difference between probing observations at successive times is representative of b, even though horizontal flow in the period between the two observations (2−10 m) results in a different snowpack and surface being probed. Complications with GPS radio transmission reduced our initial 16-point accumulation array to seven observations in one August_September period, (2009).

We compare these estimates of vertical velocity made at fixed coordinates (Eulerian frame of reference) with those from the often-used geometric relation at the surface (Eqn (4)) made at a transient marker (Lagrangian frame of reference). We employ Eqn (4) with GPS measurements of ? made at the transient ablation stake, and ? derived from a common DEM.

Error analysis

Our error analysis assumes all compounded error terms are uncorrelated. Error estimates in our measurement of cumulative mass balance from glaciological and GPS methods likewise assume annual measurements are uncorrelated. We thus estimate the uncertainty as the sum of the measurement variance.

Glaciological method

Many studies have reported a random error of _ ±0.20 m w.e. a^_1 for b _a, and this estimate is often taken to be a reliable estimate of B _a uncertainty given the spatial autocorrelation of mass-balance measurements (Reference Cogley and AdamsCogley and Adams, 1998). Recent re-analysis of glaciological measurements finds that this error is ±0.34 m w.e. a^_1 (Reference ZempZemp and others, 2013). A major sourceof error in glaciological mass-balance measurements is the spatial variability of b _a not captured by the sampling network (e.g. Reference Kaser, Cogley, Dyurgerov, Meier and OhmuraKaser and others, 2006). We quantify random error in our glaciological measurements from uncertainties in the measurements and their extrapolation. Random errors from stake measurements arise from the determination of the surface due to surface roughness and ablation caused by the stake and average ±0.10 m w.e. a^_1 (e.g. Reference Huss, Bauder and FunkHuss and others, 2009); we adopt this error term as it accords with our observations.

Accumulation measurements rely on depth and density measurements in pits and depth measurements by probing. Measurement errors of snow depth include misidentification of the previous year_s surface and determination of the undulating present-year surface; pit measurements, however, are less problematic than soundings with a probe. Penetration by a probe into underlying firn would overestimate mass balance (Reference Thibert, Blanc, Vincent and EckertThibert and others, 2008). We found that mis-identification of an overlying ice lens as the previous year_s surface was as common as probing into the underlying firn, and treated it as a random error term instead of a bias. Deviation of a probe from vertical, however, overestimates snow depth and introduces a positive bias (Østrem and Haakensen, 1999). Over three balance years we probed snow depth at 98 sites; we probed snow depth in the four cardinal directions around each site. The average standard deviation of these 392 probing observations was 0.07 m ice eq. Many studies estimate the compounded random error of accumulation measurements (typically ±0.30 m w.e. a^_1) is greater than those from the ablation area (e.g. Reference Huss, Bauder and FunkHuss and others, 2009). However, our average standard deviation of 0.07 m ice eq. from four probing measurements at 98 locations suggests a reduced error. We usean error estimate of ±0.10 m ice eq. a^_1 for all stake, pit and probing measurements of ablation or accumulation depth.

Errors in measurements of b _a also arise from an assumption of the density of ice (900 kg m^_3) and the measurement of snow density in pits. We used a large 500 cm³ tube core for the purpose of sampling snow within pits. Snow cutters have a typical measurement error of 11%, and the larger tube cutters are ofhigher precision (Reference Conger and McClungConger and McClung, 2009). Our field scale has a measurement error of ±3.3%, or ±10 g for an average sample weight of 300 g. We thus conservatively assume an error in our accumulation-area density measurements of 10%, or ±60 kg m^_3 for average conditions in our three years of study.

We therefore estimate error in b _a as

(8)

where σ? _l is the estimated error in our measurements of the length of ablation and accumulation (±0.10 m ice eq. a^_1), l is an area-weighted average ablation and accumulation measurement (2.0 m ice eq. a^_1), ?_? is an area-weighted average of error in density assumptions and measurements expressed as a conversion factor (±0.04), and ? is an area-weighted density expressed as a conversion factor (0.72).

We calculate sampling error (extrapolation from b (z) to B _a via the hypsometry) by the standard deviation of the residuals between our observations and the linear spline. These residuals are 0.36, 0.37 and 0.25 m w.e. for balance years 2009, 2010 and 2011 respectively.

Error in planimetric area (Reference Granshaw and FountainGranshaw and Fountain, 2006; Reference Bolch, Menounos and WheateBolch and others, 2010), defined as the sum of squared horizontal error in stereo model registration and digitizing error (±5 pixels), yields ±0.3% for balance years 2009 and 2010, and ±0.6% for balance year 2011.

Errors for each balance year are thus

(9)

where ? _Glac is the estimated error for glaciological B _a, ?_{b a} is measurement error and ? _Ext is extrapolation error.

Photogrammetric method

We quantify the uncertainty in dh using the standard deviation of elevation residuals of 275 points in 11 check patches (Figs 1 and 3). These residuals reveal the combined error in our stereo models (x, y, z) and the operator_s precision and accuracy in manually digitizing points (z). We follow Reference Rolstad, Haug and DenbyRolstad and others, (2009) to assess uncertainty in sequential DEM analysis when the correlation range (the extent of spatial autocorrelation) is less than the averaging area:

(10)

where σ² _A is the variance of the spatially averaged elevation difference, σ²Δ _z σis the variance of the elevation difference, A _cor is the correlation area and A is the glacier surface area. We calculate A _cor (Reference Rolstad, Haug and DenbyRolstad and others, 2009) as

(11)

where a ₁ is the correlation range determined using the Incremental Spatial Autocorrelation tool in ArcGIS 10.1, which gives the Global Moran_s I statistic over a series of increasing distances (e.g. Reference Getis and OrdGetis and Ord, 1992).

For all geodetic measurements of B _a presented below, we convert to w.e. using our B scenario that assumes 900 kg m^_3for all points below the ELA, and 750 kg m^_3 for all points above the ELA. To quantify error in this density assumption we use a range of possible density values, with the maximum density error defined by the A scenario, which assumes a density of 900 kg m^_3 for all points. The minimum density is defined by the C scenario, which assumes 900 kg m^_3 below the ELA and 600 kg m^_3 above the ELA. We consider these extrema as the plausible range of densities (±3?) and use these as a conservative estimate of density-assumption error.

We propagate the DEM error and density-assumption error as

(12)

where σ?_{B p} is the error in our photogrammetric B _a measurements (m w.e. a^_1), ? _A is the σDEM error (0.30 m ice eq.), A is area-weighted surface thickness change (1.00 m ice eq. a^_1),

?_? σis our estimate of density-assumption error from a range of possible values expressed as a conversionfactor (±0.09), and ? is an area-weighted density expressed as a conversion factor (0.81).

To estimate error in our correction for differing observation dates, we use the standard deviation of the residuals of measured and modeled ablation from our midsummer array measurements. The total number of array points is 71 from 2008,2009 and 2010, but most of these points (64) are in the ablation area. Residuals from all points have a mean of −0.01 m ice eq. and a standard deviation of 0.20 m ice eq. However, the seven points from the accumulation area have a residual mean of −0.24 m ice eq., which indicates there is a potential systematic bias in our date correction.

We estimate error in our date-corrected photogrammetric B _a measurements (? _Phtgrm) as

(13)

where σ _Corr is the estimated error in our correction for differing observation dates.

GPS method

The random error in dh arises from movement of the antenna attached to the operator_s pack, the stance of the operator on an uneven surface, and foot or ski penetration into a firn or snow surface. We thus adopt a conservative measurement error of ±0.10 m (e.g. Reference Nolan, Arendt, Rabus and HinzmanNolan and others, 2005) which is three times greater than the measurement error (±0.03 m) we observed by resurveying benchmarks.

For GPS measurements, we employ the same density assumptions and error estimates as our photogrammetric approach. As with our glaciological measurements, we estimate error in specific GPS measurements (b _GPS) as

(14)

where σ _h is the estimated error in our measurements of thickness change (±0.10 m ice eq. a^_1), h is an area-weighted average measurement of thickness change (−0.70 m ice eq. a^_1), σ_ρ is an area-weighted average of error in density assumptions and measurements expressed as a conversion factor (±0.09), and ? is an area-weighted density expressed as a conversion factor (0.81).

GPS interpolation error is approximated by the standard deviation of the residuals between our observations and the linear spline. We use ±0.57 m to estimate sampling error for all three balance years; this value arises from measurements made in 2011, which included the accumulation area.

The estimated error for the GPS method is

(15)

where σ _GPS is the estimated error for glacier-wide GPS mass balance (B _GPS) and σ _Ext is extrapolation error.

Results

The glaciological method yields two years of negative mass balance, followed by a third year of mass gain, resulting in a cumulative balance for the period 2009−11 of 0.10 ± 0.63 m w.e. (Table 2). The photogrammetric resultfor balance year 2009 (−0.15 ± 0.36 m w.e.) overlaps the glaciological method (−0.12 ± 0.39 m w.e.). Photogrammetrically based mass change for the periods 2009−11 and 2010−11, however, is more negaive than that derived by the glaciological method (Table 2)

Table 2. Estimates of glacier-wide mass balance (B _a) by the three different methods

Our GPS-derived mass-balance estimates for 2010 (−0.34 ± 0.57 m w.e.) and the period 2010−11 (_0.44 ± 0.81 m w.e.) overlap with the glaciological method (−0.31 ± 0.40 and −0.44 ± 0.26 m w.e. respectively), but the GPS results for balance years 2009, 2011 and cumulative period 2009−11 are more negative than either the glaciological or photogrammetric results (Table 2). For the period 4 August to 14 September 2009 (an additional period for which we made both glaciological and GPS measurements, with points distributed in both the accumulation and ablation areas) we find similar results of −1.02 ± 0.28 and −0.98 ± 0.28 m w.e.

We tested the reliability of using the reduced subset of photogrammetric points measured in 2011 (n = 698, 74%) to represent glacier-wide elevation change by comparing B _a for 2009 derived from 100% of the points with a synthetically thinned subset of points matching the point locations of those measured in 2011. This comparison yielded B _a of −0.15 m for both 100% and 74% coverage, suggesting that incomplete coverage of the glacier in 2011 is likely not the source of differences in the methods.

Using a spatially limited subset of geodetic point measurements (e.g. via GPS) yielded results comparable to those from measurements across the entire glacier surface (e.g. photogrammetry) if measurements were made in all elevation bins.We measured B _a as −0.08 (n = 41), −0.15 (n = 56), 0.02 (n = 37) and −0.10 m w.e. (n = 61) in four subsets, compared to −0.15 ± 0.36 m w.e. (n = 937) for the 100 m grid covering the entire glacier (Fig. 6). The subset with a slightly positive B _a (0.02, n = 37) is the only subset that does not include measurements in all elevation bins, with a data gap over the mid-glacier icefall (Fig. 3).

Fig. 6. Four subsets of the 2009 photogrammetric mass points: (a) points on a longitudinal profile along the center of the glacier; (b) points along safely walkable longitudinal profiles; (c) 37 array-point locations; and (d) points from an evenly spaced grid. Values in the upper right corner of each panel indicate B _a for the associated subset. Spatial locations of each subset are displayed in Figure 3.

Random errors

Estimated error (±0.35 m w.e.) in our glaciological measurements is dominated by the spatial variability of b _a. The precision of our accumulation-area measurements is _0.10 m, three times more precise than Reference Huss, Bauder and FunkHuss and others, (2009). The coefficient of variation (ratio of standard deviation to mean) of accumulation-area b _a, however, is twice that of the ablation area. The large error term in our glaciological measurements thus arises from the high spatial variability of b _a in the accumulation area.

Errors in our photogrammetric methodology are dominated by the precision and accuracy of our manual measurements of photogrammetric mass points (DEM uncertainty), and the number of independent samples in our sampling grid, which depends on the spatial autocorrelation of the surface elevation data. Photogrammetric measurements of dh are five times less precise than either the glaciological or GPS data. Reference Rolstad, Haug and DenbyRolstad and others, (2009) found spatial correlation within DEMs generated by automated photogrammetry at three spatial scales: hundreds of meters, a few kilometers and at tensof kilometers. They hypothesize that correlation at hundreds of meters was the result of matching errors in their automated methodology, whereas correlation at tens of kilometers was due to inaccurate georeferencing. Others (e.g. Reference Nuth, Kohler, Aas, Brandt and HagenNuth and others, 2007; Reference Barrand, James and MurrayBarrand and others, 2010) have relied on this study to justify an assumption of a correlation area of 1 km². These studies use automated techniques to investigate decadal change, in contrast to our manual methods to assess annual change. Based on our assessment of correlation range, we determined values of 1600, 700 and 1100 m in balance year 2009, and periods 2010−11 and 2009−11 respectively, resulting in correlation areas of 8.0, 1.5 and 3.8 km², which yield more conservative estimates of DEM error than if we had assumed a correlation area of 1 km².

Our errors in B _a using a GPS arise from few upper-elevation measurements, and the lack of glacier-wide elevation measurements may bias dh (z) (Figs 5 and 7). We estimate interpolation error of ±0.57 m, based on the residuals between a two-piece linear spline and observations for 2011; a shorter monitoring period (August_ September 2009) where we could make measurements in both the ablation and accumulation areas yielded a smaller interpolation error (±0.17 m). However, we cannot reliably determine dh (z) for the GPS method given the lack of suitable accumulation-area measurements. Further investigation should be made using the GPS method with measurements well distributed across the entirety of a glacier_s surface, and with the rover antenna on a stadia rod to maximize measurement precision.

Fig. 7. Same as Figure 5, but for two multi-year periods: 2009_11 and 2010_11.

To convert our geodetic measurements of dh to B _a (m w.e.) we used a bulk density of 810 ± 90 kg m^_3, which is similar to the density of 850 ± 60 kg m^_3 recommended by Reference HussHuss, (2013). Our estimated error due to assumed density (±0.05 to ±0.15 m w.e.) is lower than other error sources used in our geodetic estimates of mass change. Density errors annually vary due to the changes in extent of ablation and accumulation areas, and the magnitude of dh. However, we find error from density assumptions to be minimal on an annual basis.

We do not attempt to quantify errors due to advection of topography, but recognize that this process may inflate geodetic errors of B _a, especially in cases when the photogrammetric method does not sample the entire glacier surface. We observed topographic advection on a number of occasions as we navigated back to a point for remeasurement via GPS and found a crevasse where a relatively level surface previously existed. This same advection of surface features also likely plays a role in the apparent _blunders_ in our photogrammetric methodology, particularly over the minor, mid-glacier icefall.

Systematic differences

All three methods used in this study differ in their sources of systematic bias. Our glaciological measurements may include a positive bias from the omission of internal and basal mass balance. The photogrammetric measurements include a negative bias introduced by the manual operator and our temperature index model. In contrast, the GPS data do not include these biases.

The glaciological method suffers from potential biases associated with probing the previous summer surface and no measurement of internal and basal mass balance. Our profiles of cumulative mass balance are more positive than cumulative values for our uppermost pit alone (Fig. 7), indicating a potential positive biasof the probing measurements. We use the average difference (0.07 m w.e.) between B _a derived from all points and only from pits as an approximation of this potential bias.

Sinking (self-drilling) of ablation stakes may produce a negative bias in the glaciological method (Reference Riedel, Wenger and BowermanRiedel and others, 2010). We tested this possibility by measuring ablation at adjacent stakes, one with an insulated cap at the base of the stake and one without. No differences were noted, so we assume that self-drilling did not occur.

The glaciological measurement does not capture internal and basal mass balance. These mass changes are glacier-dependent, with internal and basal accumulation playing a more significant role in cold, continental climates (e.g. Storglacia_ren, Sweden; Reference ZempZemp and others, 2010), whereas internal and basal ablation dominates in warm, maritime climates (e.g. Franz Josef Glacier, New Zealand; Reference Alexander, Shulmeister and DaviesAlexander and others, 2011). Castle Creek Glacier is a temperate, continental glacier anddoes not experience a high geothermal flux. Estimates of internal ablation for temperate glaciers vary dramatically from1 cm w.e. a^_1 (e.g. Reference Thibert, Blanc, Vincent and EckertThibert and others, 2008) to 2.5 m w.e. a^_1 at lower elevations (Reference Alexander, Shulmeister and DaviesAlexander andothers, 2011). We conservatively estimate internal ablation and accumulation to respectively be 10% and 4% of B _a (Reference ZempZemp and others, 2010; Reference Alexander, Shulmeister and DaviesAlexander and others, 2011). Combining these values with our probing bias, we estimate positive biases of 0.08, 0.09 and 0.10 m w.e. for 2009, 2010 and 2011 respectively, and 0.27 and 0.19 m w.e.for the periods 2009−11 and 2010−11 respectively.

We corrected for stereo model bias in our photogram-metric measurements of B _a, but neglect densification and subglacial erosion, processes that we assume to be negligible in our study. To assess bias in the temperature index model, we differenced modeled from observed ablation at array points for August 2008, 2009 and 2010. Themean of the residuals for 64 observations in the ablation area is 0.02 m w.e., whereas it is − 0.24 m w.e. for seven observations in the accumulation area, suggesting a possible negative bias. However, this potential bias is perhaps related to error in our probing measurements at the seven accumulation-area locations, and the limited number of observations. When comparing observed and modeled ablation measurements, we found melt factors derived from mass-balance measurements (Reference Shea, Moore and StahlShea and others, 2009) yielded better agreement when using Place Glacier melt factors (−0.01 ± 0.20 m w.e.) thanthose from Peyto Glacier (+0.23 ± 0.28 m w.e.). We are hesitant to rely on Castle Creek Glacier melt factors given our short period of observation and few measurements for snow surfaces (n = 7). The use of melt factors derived from either Place or Peyto Glacier does not substantially affect our estimated photogrammetric mass change, however (Table 3).

Table 3. Varying photogrammetric B _a results when using different melt factors (mm w.e.°C^_1 d^_1) in our temperature index model for correcting for dates of photography. Percentages in parentheses indicate percent differences from results using Place Glacier melt factors.

After correcting for stereo model bias, manual digitization can still introduce significant bias in photogrammetric measurements of dh (Reference McGlone, Mikhail and BethelMcGlone and others, 2004). Our repeated measurements from the same stereo models of 275 check points assess this bias (Fig. 8), and yield a potential bias of −0.26 ± 0.46, −0.34 ± 0.58 and 0.09 ± 0.74 m for stereo models from 2008, 2009 and 2011 respectively. The mean of all 825 points yields a bias of −0.15 m. These results indicate a potential manual-operator bias, but one that is inconclusive in terms of both sign and magnitude. Additionally, measurement uncertainty varies according to glacier surface characteristics. Repeat measurements of surface elevation from 2009 models yielded residualsof 0.01 ± 0.44, 0.16 ± 0.66 and −0.15 ± 1.44 m ice eq. for ice, firn and snow respectively. This error, however, did decrease with repeated measurements, revealing the importance of operator experience (Fig. 4). We found GPS and photo-grammetric point measurements of dh accord in the ablation area (Figs 5a and 7).

Fig. 8. Frequency distributions of repeat-measurement residuals of 275 check points for stereo models from 2008, 2009 and 2011 aerial photography used to constrain potential error and bias of the analyst in manual photogrammetry. Each distribution displays surface-elevation residuals, defined as the difference in surface elevation of the same horizontal coordinates from the same stereo model and by the same operator, but from initial measurements and a repeat measurement at a later date.

We are unaware of other studies assessing the potential biases of GPS B _a measurements, and we do not attempt to quantify GPS bias, which would require a complete understanding of a glacier_s spatially varying vertical velocity.

The two geodetic methods yield results that are comparable to, or significantly more negative than, those of the glaciological method. These results thus suggest a greater mass loss for the glacier using geodetic measurements than was measured using the glaciological method. Our findings accord with work at South Cascade Glacier (Reference KrimmelKrimmel, 1999), whereas an analysis of 29 glaciers by Reference CogleyCogley, (2009) found only a minor, insignificant systematic difference (−0.07 m w.e. a^_1).

Estimation of vertical velocity

Our estimates of ablation, dh, and estimated vertical velocity for the ablation area vary for 2008, 2009 and 2010 (Fig. 9). Our method of estimating vertical velocity (Eqn (7)) yields values that average 0.30 m ice eq. more emergence than estimates from Eqn (4). GPS measurements at seven points in the accumulation area in 2009 produce a submergence estimate of −0.38 ± 0.15 m ice eq.

Fig. 9. Spatial distribution of mass balance, thickness change and vertical velocity (ice eq.) from the 20-point ablation-area array for August in 2008, 2009 and 2010. We derived mass balance from stake measurements, thickness change from RTK GPS, and vertical velocity as the difference of the two. Locations of stakes and GPS measurements are indicated by black dots.

We do not quantify errors in our estimate of vertical velocity. Using our approach, this error term will be dominatedby advection of surface topography, but we lack a detailed estimate of the glacier_s surface roughness. Advection of an irregular surface (0.5−2 m roughness) will greatly exceed our estimated precision in measuring height change from ablation stakes or with a GPS (±0.10 m ice eq.). We suspect that some of the variability in vertical velocity (Fig. 9) may arise from advection of topography.

Discussion and Recommendations for Future Mass-Balance Monitoring

Limitations and advantages are inherent in each method to assess glacier mass change. Surface mass balance estimated by the glaciological method should be continued for index glaciers and potentially expanded for under-represented mountainous regions. A key benefit of the glaciological method is its ability to record b (z), a prerequisite for modeling the response of a glacier to changes in climate (e.g. Radi? and Hock, 2011). Unfortunately, the glaciological method also suffers from a number of logistical shortcomings, which include the time- and energy-intensive nature of the method, and its limitation to glaciers that are relatively small and safe for travel.

We recommend increased use of the photogrammetric method to monitor B _a. When GCPs alreadyexist, fieldwork is not necessary, which can significantly reduce measurement time and cost. However, a lack of density measurements remains a source of error. Remote-sensing geodetic methods afford the best possibility to monitor representative glaciers, including those that are large, complex and difficult to visit. Additionally, such methods enable the monitoring of more glaciers in a region, avoiding expensive and time-intensive field studies.

Poor contrast, particularly for glaciers covered by fresh snow cover, limits the use of manual or automated feature extraction from aerial photography (e.g. Reference KrimmelKrimmel, 1999; Reference Bamber and RiveraBamber and Rivera, 2007). Multispectral aerial photography or high-resolution, multispectral satellite imagery, however, improves contrast in the accumulation area of glaciers. The inclusion of the near-infrared band in our 2011 stereo models, for example, enabled us to measure elevation for 74% of the glacier_s surface, which was covered by fresh snow.

We advocate the adoption of our GPS technique to measure B _a. Previous studies concluded that geodetic measurements alone cannot be used to measure B _a in the absence of a knowledge of dynamics (e.g. Reference Hagen, Eiken, Kohler and MelvoldHagen and others, 2005), but we find that a well-distributed subset of point measurements can adequately mitigate the confounding role of vertical velocity to yield reliable estimates of B _a. In the case of Castle Creek Glacier, we find a sample density of 4 km^_2 can be used to derive B _a using geodetic-grade GPS receivers (Fig. 6). If completed for a well-distributed subset of points across the glacier surface, the GPS method may circumvent some of the limitations inherent in the other two approaches, but more studies are required to quantify errors and bias inherent in the GPS approach.

Use of the in situ GPS method is restricted to those glaciers that are accessible, and safe for travel. Additionally, measurements yield dh, which is modulated by vertical velocity, making the GPS method ill-suited for determination of b _a and b (z) (cf. Figs 5 and 7). GPS is advantageous on larger glaciers where the glaciological method is impractical due to necessary commitments of time and energy. Use of the GPS method at the end of the accumulation season, when glacier surfaces are covered with snow, may enable travel and measurement on surfaces inaccessible at the end of the balance year. We recommend future efforts to refine the GPS method be undertaken at an established field station with a source of power adequate for the high-powered base station radio (35 W).

Combining at-a-point glaciological and GPS measurements provides one method to assess the spatial and temporal changes of vertical velocity for a glacier. Our methodology to estimate vertical velocity (Eqn (7)) compares well with that ofthe geometric relation at the surface (Eqn (4)). However, a major shortcoming is that it is field-intensive, whereas it is possible to employ the kinematic boundary condition remotely (e.g. Reference Gudmundsson and BauderGudmundsson and Bauder, 1999). Furthermore, the use of an Eulerian frame of reference in this method imparts a potential error due to advection of topography, an error that is likely to be on the order of 0.5 m but may be in excess of 5 m in areas of complex surface topography. Changes in ice flux have been found to be significant in determining recent thinning (Reference Berthier and VincentBerthier and Vincent, 2012), and combining GPS and glaciological measurements allows insight into the fine-scale structure of seasonal ice dynamics. Additionally, our in situ methodology enables the estimation of submergence in the accumulation area. Remote-sensing studies, which derive submergence from measurements of horizontal motion and elevation change (Eqn (4)), often fail in the accumulation area due to a lack of surface features to track to determine surface velocity. However, a study of the potential error imparted by advection of topography is required.

Our future efforts at Castle Creek Glacier will include continued monitoring of annual length change (Reference Beedle, Menounos, Luckman and WheateBeedle and others, 2009), and glaciological measurements of annual mass balance. Additionally, we plan to acquire aerial photographs ofCastle Creek Glacier for geodetic measurements of annual and cumulative mass balance.

There is a need to continue long-term mass-balance measurements, resume interrupted series, expand to important regions and more representative glaciers, and improve error analysis (Reference Fountain and VecchiaFountain and others, 1999; Reference Zemp, Hoelzle and HaeberliZemp and others, 2009). Geodetic methods provide a valuable measure of B _a that is complementary to those of the glaciological method. Photogrammetric and GPS methods provide means to improve our understanding of glacier change, and to help understand and predict the fate of mountain glaciers.