An evaluation of mass-balance methods applied to Castle creek Glacier, British Columbia, Canada

Abstract

We estimate the glacier mass balance of a 9.5 km² mountain glacier using three approaches for balance years 2009, 2010 and 2011. The photogrammetric, GPS and glaciological methods yielded sampling densities of 100, 5 and 2 points km^-2, with measurement precisions of ± 0.40, ± 0.10 and ± 0.10 m w.e. respectively. Our glaciological measurements likely include a positive bias, due to omission of internal and basal mass balance, and uncertainty in determining the interface between snow and firn with a probe (±0.10 m w.e.). Measurements from our photogrammetric method include a negative bias introduced by the manual operator and our temperature index model used to correct for different dates of imaging (0.15 m w.e.), whereas GPS measurements avoid these biases. The photogrammetric and GPS methods are suitable for estimating glacier-wide annual mass balance, and thus provide a valuable measure that complements the glaciological method. These approaches, however, cannot be used to estimate mass balance at a point or mass-balance profiles without a detailed understanding of the vertical component of ice velocity.

Introduction

Geodetic estimates of glacier mass balance quantify changes in surface elevation over a given time period and employ density assumptions for snow and ice. Most geodetic studies assess mass change over periods of a decade or more (e.g. Arendt and others, 2002; Schiefer and others, 2007). Some studies, however, use geodetic methods to determine mass balance on shorter temporal scales. For example, Meier and Tangborn, (1965) used aerial photography taken three years apart to estimate the annual mass balance and investigate the short-term ice dynamics of South Cascade Glacier, Washington, USA. Tangborn and others, (1975) made repeat survey measurements of a 112-point grid on South Cascade Glacier during one ablation season. They found that their geodetic measurements accorded with those of the glaciological method (Cogley and others, 2011). Rasmussen and Krimmel, (1999) used aerial photogrammetry to derive annual specific mass balance over a portion of South Cascade Glacier for balance years 1993 and 1994. They identified the utility of geodetic measurement, but their study also revealed the potential for systematic biases in both geodetic and glaciological mass-balance measurements. Krimmel, (1999) presented a comparison of South Cascade Glacier annual mass balance derived from both glaciological and photogram-metric methods for 12 balance years (1986−97). The South Cascade Glacier cumulative geodetic balance was more negative than glaciological balance by ˜0.25 m w.e. a^_1, indicative of bias in one or both of the methods, and possibly due to basal melt, density estimates and melting of stakes into the glacier surface. The study concluded that photogrammetry can be used to determine annual mass balance if geodetic control is consistent for stereo models.

In the late 1990s, GPS methods were tested to measure mass balance (Eiken and others, 1997; Gandolfi and others, 1997; Jacobsen and Theakstone, 1997). These studies demonstrated the potential of using GPS in kinematic mode, whereby pointmeasurements are made continuously at a preset time interval, to record glacier surface elevation. Without using real-time kinematic (RTK) GPS, where a real-time differential correction is received from the base station through radio transmission allowing centimeter accuracy, however, these early studies were only able to reoccupy points on a glacier to within tens of meters. Hagen and others, (1999) used GPS to measure the mass balance of Kongsvegen, Svalbard, for the period1991_95 and found that the results fit well with field measurements of a glacier with negligible vertical velocity. Previous work noted the rapidity and accuracy of kinematic GPS measurements, which might lead to an increase in the number ofmonitored glaciers (Hagen and others, 1999; Theakstone and others, 1999). Hagen and others, (2005) presented multi-year comparisons of GPS profiles along central flowlines of three Svalbard glaciers. They concluded that changes in glacier geometry cannot be used to assess mass balance without independent knowledge of vertical velocity. In that study, RTK GPS was not used and reoccupation of previously measured points was made only to within 30_90 m. Other researchers have usedRTK GPS to accurately reoccupy and remeasure survey points on a glacier to assess multi-year mass balance (e.g. Nolan and others, 2005).

Here we compare glaciological, photogrammetric and RTK GPS methods to estimate annual glacier mass balance. Our objectives are to: (1) test the RTK GPS method; (2) investigate potential error in the glaciological method by comparison with two geodetic methods; and (3) make recommendations for future mass-balance monitoring.

Theory: conservation of mass and vertical velocity

Geodetic measurements of glacier thickness change incorporate two dominant terms, surface mass balance and the vertical component of ice flow, necessitating consideration of flux divergence. Conservation of mass at a point on the surface of a glacier can be stated as

(1)

where is the rate of thickness change, is the specific surface mass-balance rate, ρ is density, and is a flux-divergence term (e.g. Rasmussen and Krimmel, 1999; Cuffey and Paterson, 2010). Implicit in Eqn (1) are the assumptions that densification, internal and basal mass balance, isostatic displacement, and erosion of the bed surface are negligible. When integrated across the entire glacier surface, and assuming a negligible influence from changing surface area and no flux across the glacier boundary (e.g. from avalanching or calving), flux divergence is zero, yielding

(2)

where and are glacier-wide integrations of thickness change and surface mass balance respectively.

Flux divergence, however, is not zero for a given point, and the vertical component of ice velocity at the surface (vertical velocity, w _s) plays a confounding role in deriving from measurements of . Previous studies use Eqn (1) to estimate from geodetic measurements and include a complete treatment of the flux divergence term in Eqn (3) , with consideration of the vertical profile of velocity and the componen of surface flow due to sliding (e.g. Gudmundsson and Bauder, 1999; Rasmussen and Krimmel, 1999).

In studies to quantify w _s, some efforts include consideration of subsurface glacier flow or assume steady-state conditions (e.g. Reeh and others, 1999, 2003), while others neglect the vertical profile of velocity (e.g. Meier and Tangborn, 1965; Holmlund, 1988; Pettersson and others, 2007), and rely on the kinematic boundary condition at the glacier surface:

(3)

where u _s and v _s are the horizontal components of ice velocity atthe glacier surface S (Cuffey and Paterson, 2010). Often, Eq (3) is used to estimate w _s in only the ablation area and is further reduced to a simple geometric expression (e.g. Meier and Tangborn, 1965; Holmlund, 1988; Pettersson and others, 2007):

(4)

where is thickness change measured at a marker (usually a stake) moving with the glacier ice, u _s is oriented along the flow, v _s is assumed to be zero, and ? is the slope of the glacier surface (Cuffey and Patterson, 2010).

Previous work that employs either Eqn (1) to estimate from geodetic measurements, or Eqns (3) and (4) to estimate w _s relies on a Lagrangian frame of reference, whereby horizontal flow and occasionally are measured at a marker (e.g. a stake) as it moves with the ice. In an Eulerian frame of reference, wheremeasurements are made at fixed coordinates, changes as the sum of ? and w _s (Cuffey and Paterson, 2010). However, horizontal ice flux at the surface (u _s and v _s in Eqn (3)) advects glacier surface topography through locations where geodetic repeat measurements are made. We discuss this as a source of error below, but otherwise neglect advection of topography in our geodetic measurement of mass balance and estimation of w _s. Omission of these horizontal velocity terms yields

(5)

which may be rearranged to solve for specific mass balance

(6)

or for vertical velocity

(7)

In addition to advection of topography, Eqns (5_7) also neglect densification, internal and basal mass balance, isostatic displacement, and erosion of the bed surface.

Vertical velocity at the surface (w _s) is typically negative or downward (submergence) in the accumulation area and positive or upward (emergence) in the ablation area. To avoid ambiguity in our discussion of these, we use the term emergence to refer to positive (upward) flow, submergence to refer to negative (downward) flow, and vertical velocity as a general term without a specified sign.

Methods

We measured the mass balance of Castle Creek Glacier, Cariboo Mountains, British Columbia, Canada (Beedle and others, 2009). This 9.5 km² mountain glacier flows north for 5.9 km, has an elevation range of 2827−1810 m a.s.l. and contributes meltwater to Castle Creek, a tributary of the Fraser River (Fig. 1).

Fig. 1. Location of Castle Greek Glacier and nearby glaciers with reference mass-balance series(inset).The inset also includes nearby towns,major cities and two stems of the Fraser River, into which Castle Creek meltwater flows.Check patches are labeled by their respective mean elevation.

Our terminology and notation follow the recommendations of Cogley and others, (2011). The glaciological method measures surface mass balance, whereas geodetic methods measure elevation change of the glacier surface, from which volume and mass changes are estimated. Geodetic measurements sum the effects of mass change at the surface, internal to the glacierand at the bottom of the glacier. For simplicity, however, we refer only to annual mass balance regardless of method. A mass-balance profile, b (z), is defined as the variation of mass balance with elevation. We use dh as notation for thickness change, and dh (z) to refer tothe profile of thickness change with elevation. All measurements of glacier-wide mass balance presented are conventionalbalances whereby values are averaged over a changing glacier surface (Elsberg and others, 2001), allowing direct comparison of glaciological and geodetic methods. The changing ice geometry is based on digital elevation models (DEMs) from 2008, 2009 and 2011 aerial photos. We estimate the 2010 glacier geometry by linear interpolation from the 2009 and 2011 DEMs. Castle Creek Glacier surface areas used to convert from volume change to specific (per unit area) mass balance are 9.56 km² (2009), 9.52 km² (2010) and 9.49 km² (2011).

Glaciological method

Snow pits and probing are used to directly measure accumulation, whereas stakes measure surface ablation (e.g. Østremand Brugman, 1991; Kaser 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. Fountain 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. Cogley 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 (Fountain 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. Schiefer and Gilbert, 2007; Barrand 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; Schiefer and Gilbert, 2007; Schiefer 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 (McGlone 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. Rasmussen 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 (Hock, 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 (Shea 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. Huss, 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; Bader, 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. Zemp 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; Cogley, 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. Nolan 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. Rasmussen 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 (Cogley and Adams, 1998). Recent re-analysis of glaciological measurements finds that this error is ±0.34 m w.e. a^_1 (Zemp 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. Kaser 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. Huss 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 (Thibert 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. Huss 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 (Conger 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 (Granshaw and Fountain, 2006; Bolch 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 Rolstad 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 (Rolstad 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. Getis 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. Nolan 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 Huss 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. Rolstad 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. Nuth and others, 2007; Barrand 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 Huss, (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 (Riedel 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; Zemp and others, 2010), whereas internal and basal ablation dominates in warm, maritime climates (e.g. Franz Josef Glacier, New Zealand; Alexander 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. Thibert and others, 2008) to 2.5 m w.e. a^_1 at lower elevations (Alexander andothers, 2011). We conservatively estimate internal ablation and accumulation to respectively be 10% and 4% of B _a (Zemp and others, 2010; Alexander 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 (Shea 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 (McGlone 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 (Krimmel, 1999), whereas an analysis of 29 glaciers by Cogley, (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. Krimmel, 1999; Bamber 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. Hagen 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. Gudmundsson 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 (Berthier 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 (Beedle 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 (Fountain and others, 1999; Zemp 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.