Using surface velocities to calculate ice thickness and bed topography: a case study at Columbia Glacier, Alaska, USA

Abstract

Information about glacier volume and ice thickness distribution is essential for many glaciological applications, but direct measurements of ice thickness can be difficult and costly. We present a new method that calculates ice thickness via an estimate of ice flux. We solve the familiar continuity equation between adjacent flowlines, which decreases the computational time required compared to a solution on the whole grid. We test the method on Columbia Glacier, a large tidewater glacier in Alaska, USA, and compare calculated and measured ice thicknesses, with favorable results. This shows the potential of this method for estimating ice thickness distribution of glaciers for which only surface data are available. We find that both the mean thickness and volume of Columbia Glacier were approximately halved over the period 1957–2007, from 281 m to 143 m, and from 294 km³ to 134 km³, respectively. Using bedrock slope and considering how waves of thickness change propagate through the glacier, we conduct a brief analysis of the instability of Columbia Glacier, which leads us to conclude that the rapid portion of the retreat may be nearing an end.

1. Introduction

Knowledge of glacier volume and ice thickness distribution is essential for hydrological applications, ice flow modeling, assessing the impact of climate change on glaciers, and sea-level rise predictions, among other applications. Direct measurements of ice thickness at a point (e.g. from a borehole) or along a track (e.g. using radio-echo sounding or seismic methods) are time-consuming and expensive. Direct measurements of total ice volume generally require many such points or tracks. In addition, errors in interpolation or extrapolation from along-track measurements of ice thicknesses can introduce large anomalies in calculated ice-flux divergence if the interpolation technique used does not conserve mass or ice flux (Seroussi and others, 2011). A mass-conserving method for interpolating ice thickness has been developed by Morlighem and others (2011). Techniques using radar tomography and interferometry have recently been developed, but they have mostly been applied to portions of the Greenland ice sheet with flat surfaces, and their performance and applicability to outlet and tidewater glaciers is not yet known (Paden and others, 2010).

Previous studies have focused on calculating ice volume and ice thickness distribution of glaciers without direct measurement. One of the simplest methods is volumearea scaling (e.g. Bahr and others, 1997; Radić and others, 2008). These methods are easily implemented, requiring information about glacier area and parameterization of the shape of the bed topography. In general, the shape of the bed of individual glaciers is not well known, but the method has proven useful for characterizing regional ice volumes (Radić and Hock, 2010).

Other studies have inferred ice thickness distribution through direct evaluation of the mass continuity equation (e.g. Rasmussen, 1988; Morlighem and others, 2011), while others use a simplification of the equation to overcome data gaps (e.g. Fastook and others, 1995; Warner and Budd, 2000; Farinotti and others, 2009a), or through applications of the shallow-ice approximation (e.g. Li and others, 2011). More recently, a method has been proposed which uses neural networks along with simplifications of the mass continuity equation to estimate the bed topography and ice volumes of entire regions, where little more than surface topography might be known (Clarke and others, 2009).

Here we propose a new method for calculating ice thickness, by evaluating the mass continuity equation between adjacent flowlines, rather than through a local solution or on a large grid. With velocity fields that cover some portion of a glacier, a digital elevation model (DEM) of the glacier surface, rates of surface mass balance and thinning, and knowledge of the ice thickness at the boundary of the domain of interest, it is possible to calculate the ice thickness distribution of the glacier over the region covered by the velocity field used. Because no assumption is made about the ice flux through the terminus of the glacier, this method is directly applicable to both land and marine/lake-terminating glaciers.

We apply this method to Columbia Glacier, a large tidewater glacier in Alaska, USA. We investigate the sensitivity of the calculated ice thicknesses to flowline separation and input variables, and produce a new, highresolution gridded bed topography map for Columbia Glacier. We then use this bed topography map to calculate ice thickness, evaluating the success of this computation through direct comparison with radar profiles of ice thickness. Using this, we calculate total ice volume in 2007 and 1957. Finally, we examine spatial patterns of volume change at Columbia Glacier, the stability of the current tidewater extent, and potential analogues to Columbia Glacier with both Glacier Bay and Icy Bay, Alaska.

2. Columbia Glacier

Columbia Glacier (Fig. 1) is a large tidewater glacier, ~30 km to the west of Valdez, Alaska. At present (2011), the glacier has a surface area of ~910km², ranging from sea level to 3700ma.s.l. It reached its most recent extended position around 1850, terminating near the northern edge of Heather Island (Calkin and others, 2001). Around 1980, it began a rapid retreat that continues through the present (e.g. Meier and Post, 1987; O’Neel and others, 2005; Pfeffer, 2007; Walter and others, 2010; Rasmussen and others, 2011). Since then, it has retreated >23 km. During the course of the retreat, ice speeds near the terminus have exceeded 25 m d^-1 (O’Neel and others, 2005). Since 1976, the glacier has been monitored by the United States Geological Survey (USGS), as well as the University of Colorado at Boulder (CU-Boulder), with aerial photogrammetry, and since 2004 with time-lapse photogrammetry (e.g. Krimmel, 2001; O’Neel and others, 2005; Walter and others, 2010). It is the single largest contributor of Alaska glaciers to sea-level rise, accounting for ~6% of the total regional contribution over the period 1962–2006 (Berthier and others, 2010). In the periods 1957–2007 and 2007–11, the glacier lost ~132km² (~12%) and ~20km² of area, respectively. Columbia is the best-studied example of tidewater glacier retreat in the world. Ice flow, ice discharge and Columbia’s tidewater retreat are all extensively documented, providing rich insight into the underlying processes that modulate tidewater glacier behavior and stability.

Fig. 1. Columbia Glacier, showing 1957 glacier extent in gray. Contours indicate 1957 surface elevations. Open circles indicate distance, ξ, from the head of the glacier following the central flowline defined by Meier and others (1985). Extent of measured bathymetry is shown, as well as extent of bed topography map by Engel (2008). Thick black lines indicate location of radar tracks, numbered 1–5. Location of terminus in June 2011 is shown as a dashed line.

3. Datasets

For the proposed method, the required input data are glacier surface topography in the form of a DEM, surface velocities, surface mass-balance rates, rates of ice surface elevation change, and ice thickness at the boundary of the domain of interest. The method is validated using known ice thicknesses and bed elevations. A summary of available datasets is shown in Table 1. A center-line coordinate system ξ, where ξ = 0 km at the head of the glacier, was proposed by Meier and others (1985), and is used here to present center-line data.

Table 1. Overview of datasets, available for Columbia Glacier, used in this study. Numbers in parentheses indicate number of datasets available during given time period. Datasets without citation are unpublished

3.1. Map data and digital elevation models

Full glacier coverage DEMs are available for 1957 and 2007. The 1957 DEM is a digitized United States Geological Survey (USGS) 1:63 360 topographic map, with an associated glacier outline. The 2007 DEM is generated using a SPOT (Satellite Pour l’Observation de la Terre) panchromatic image acquired on 22 September 2007, and has a spatial resolution of 40 m (Korona and others, 2009).

The 2007 Columbia Glacier outline was manually digitized from the same SPOT panchromatic image as the 2007 DEM. Additionally, two Landsat images acquired on 13 September 2006 and 3 July 2009 were used to digitize the glacier outline where it was obscured by cloud cover or snow in the SPOT image. A watershed algorithm was used to find the Columbia Glacier flow basin using the 1957 DEM (Kienholz, 2010). This was then manually corrected to represent 2007 flow divides, using the 2007 DEM. Terminus position for 2007 was mapped to correspond to the 2007 DEM.

Smaller spatial coverage DEMs, determined using photogrammetric techniques, are associated with each surface velocity dataset, and are available yearly for the periods 1976–2001 (Krimmel, 2001) and 2004–10. The majority of these DEMs are limited to the eastern trunk of the glacier and the former terminus area (the current proglacial fjord), covering ~265 km², or 30% of the glacierized area in 1957.

3.2. Fjord bathymetry and glacier bed elevation

High-resolution bathymetric data are available in the large proglacial fjord that has opened up since the onset of retreat, up to the location of the 2004 terminus (Noll, 2005). Bathymetric data were collected in 2005 by the US National Oceanic and Atmospheric Administration (NOAA) ship Rainier using multibeam sonar. These data cover 6% of the 1957 glacier extent (Fig. 1) and have a reported accuracy of 7% of water depth.

A partial bed topography map, covering 16% of the 1957 glacier extent, is also available above the location of the 2004 terminus, but has limited accuracy (Engel, 2008). In addition, the coverage of this map is limited to the lower part of the eastern trunk of the glacier, a substantial portion of which has since calved off. The upstream limit of this map is ~35 km from the head of the glacier, and 13 km from the June 2011 terminus. This bed topography map was produced using the mass continuity equation, surface velocity fields, and surface elevation change data derived from center-line altitude profiles. Engel (2008) also assumes that surface mass balance is zero over the period of retreat. Details of the method used to produce this topography map are provided by O’Neel and others (2005) and Engel (2008).

3.3. Ice thickness

Ice thickness data were collected on 22 April 2010, from a de Havilland Otter airplane towing an impulse radar with 2 MHz center-frequency antennas, in a configuration similar to that used in 2006 to sound Bering and Malaspina Glaciers, Alaska (Conway and others, 2009). Waveforms were geolocated using GPS on board the airplane. Ice thickness at each location was estimated from the difference in two-way travel time from the surface to the bed, assuming a wave speed in ice of 170 m μs^-1. The resolution (1/4 wavelength) when using a 2 MHz antenna is ~20 m in ice. Uncertainty in the estimate of ice thickness comes from uncertainty in the wave speed (~2 m μs^-1, which corresponds to 1.2% of the ice thickness, or 7.2 m for 600 m thick ice), and from picking the travel time to the surface and to the bed. This uncertainty when picking the two-way travel time of the maximum reflection amplitude with our 2 MHz system is about 0.1 μs^-1, which corresponds to 8.5 m. The uncertainty would increase in the presence of side reflectors, which can introduce ambiguity in bed pick positions. Assuming these uncertainties are not correlated, the combined uncertainty for 600 m thick ice is ~14 m. This uncertainty is comparable in magnitude to differences in estimates of thickness at crossover points between tracks, which typically agree to better than 10 m. Ice thickness data are available through the Advanced Cooperative Arctic Data and Information Service portal (http://www.aoncadis.org/home.htm) under ‘Collaborative Research IPY: Dynamic Controls on Tidewater Glacier Retreat’.

3.4. Surface mass balance

Estimates of the vertical profile of annual surface mass balance during 1948–2007 were made by Rasmussen and others (2011) using upper-air temperatures and winds. Their model was calibrated with 67 mass-balance measurements made mainly in 1977 and 1978. It assumes that both precipitation and the positive degree-day factor are constant with altitude. Accumulation increases with altitude, however, because both precipitation and the fraction of precipitation falling as snow increase. Ablation decreases with altitude because the number of positive degree-days decreases. The model does not explicitly consider radiation or redistribution of fallen snow. Reported root-mean-square error (RMSE) is 1.0 m w.e. a^-1, with coefficient of determination r² = 0.88 (Bevington, 1969). Mean values over the period range from –6.7mw.e. a^-1 at 100ma.s.l. to 5.4mw.e. a^-1 at 3700ma.s.l., and the mean equilibrium-line altitude is 940ma.s.l. over the period. Over the period 1948–81, the modeled surface mass balance of the entire glacier is positive (~0.8 km³ w.e. a^-1). The modeled surface mass balance remains positive through 1995 (~1.2 km³ w.e. a^-1), and then turns negative after 1995 (~-0.4 km³ w.e. a^-1).

3.5. Surface velocities

Surface velocity fields covering part of the glacier have been calculated over the periods 1976–2001 (Krimmel, 2001; O’Neel and others, 2005) and 2004–10 using aerial photogrammetry. Krimmel (2001) presents a record of 121 flights made between 1957 and 2001 that are unevenly spaced in time. Manual feature tracking is used to determine average velocity fields between successive pairs of images, as well as mean surface topography between successive pairs of images. Spacing between image pairs ranges from 10 to 138 days. Time of year for flights typically falls into two ranges, January-March and July-October. To maintain consistency with the full-coverage DEMs, only velocity sets from mid-July onward are considered, resulting in a set of 29 flight pairs.

For flights made since 2004, photogrammetric analysis is done using automated feature tracking, allowing for a much denser spatial coverage (Ahn and Howat, 2011). Densities of velocity measurements range from approximately one measurement for every 20 gridcells, to nearly one measurement per 100 m x 100 m gridcell.

Full glacier coverage velocity fields for 2011 were determined with standard feature- and speckle-tracking techniques (Joughin, 2002; Strozzi and others, 2002) applied to TerraSAR-X synthetic aperture radar data. Formal errors, based on the statistics of the matches, are computed and are generally small (1–10 m a^-1). There may be velocity errors of less than ~3% due to error in the DEM-derived slopes used to correct for vertical motion (Joughin and others, 1996).

4. Method

The method used here is based on principles first described by Rasmussen (1988), and also described by O’Neel and others (2005). We solve for ice thickness between adjacent flowlines, in order to avoid the consideration of cross-flow gradients. We choose this method of evaluating the continuity equation because its implementation is relatively simple.

First, consider a cell (of area S), through the lateral boundaries of which there is no flow (Fig. 2). In practice, this could be two adjacent flowlines, an ice catchment basin such as those employed by Farinotti and others (2009a) or an entire glacier. If the latter two examples are used, care must be taken with the spatially averaged surface velocities and ice thicknesses, as discussed below.

Fig. 2. Schematic illustration of glacier surface. Black arrows indicate flow vectors, thick black lines indicate flowlines, and black lines transverse to flow are boundaries of cells. Inset: Map view of a cell, with area S, through the lateral boundaries of which there is no flow. R and P designate the upstream and downstream boundaries of the cell, respectively. v_in and v_out indicate the ice velocity at the upstream and downstream boundaries of the cell, respectively.

Assuming a constant ice density, the mass conservation equation is given by

(1)

where H is the ice thickness, ∙ is the ice flux divergence and is the climatic-basal mass-balance rate (i.e. the sum of the surface, internal and basal mass balances; Cogley and others, 2011). In general, and at Columbia Glacier in particular, the climatic-basal balance rate is dominated by the surface mass-balance rate , so we replace with from here.

Using Gauss’s theorem, rearranging and then integrating Eqn (1) over the surface area S of any cell, we obtain

(2)

where q _in and q _out are the ice fluxes through the upstream and downstream boundaries (R and P; Fig. 2) of the cell, respectively. Next, we consider the fact that the integral of the velocity profile over the ice thickness H is

(3)

where γ ϵ [0.8,1] is a factor relating the surface velocity v _sfc with the depth-averaged velocity (Cuffey and Paterson, 2010).

Along R, we have the following (the result for P is similar; here, x is defined along the boundary P or R):

(4)

where v_sfc is measured normal to R, and W_R is the length of the boundary R. For simplicity, we assume γ is spatially and temporally constant. This is not necessarily the case, but sensitivity experiments (see Section 6.1) indicate that changing values of γ have little impact on the calculated ice thickness in this study.

We now desire to solve Eqn (4) for . From the mean value theorem, we know that there exists x₀ € [0, W _R] such that the following is true:

(5)

What we then seek is the spatial scale on which both H and v_sfc are approximately constant, so that any variation in H or v_sfc is small, and therefore . Clearly, this is not the case when the cell covers most or all of the glacier width, so care must be taken if the cell covers a substantial portion of the glacier width. An analysis of the sensitivity of the calculated ice thickness to initial separation distance of adjacent flowlines follows.

On a scale where H and v_sfc are approximately constant, is approximated by Hv_{s fc}, and we can combine Eqns (2) and (4) to solve for H at the upstream boundary, R, given the ice flux q_out at the upstream boundary P:

(6)

We can also solve for H at the downstream boundary, starting at the upstream boundary of the glacier. The direction of computation (up-or down-glacier) will depend on the data availability. This results in a calculation of the ice thickness at many cross sections between a pair of adjacent flowlines; with many such pairs of flowlines, we can extend the calculation to cover the entire region of interest.

To solve Eqn (6), we need the following information: values of surface mass balance values of surface elevation change ∂H/∂t, a surface velocity field, and an estimate of the ice thickness at the boundary of the domain of interest. Values of can be interpolated from measurements, or obtained from model outputs. Considerable thought has been put into the differences between rates of surface elevation change derived from profiles versus those derived from DEM differencing (e.g. Arendt and others, 2002; Larsen and others, 2007). Thinning rates derived from profiles can be used, but may not be representative of the margins of the glacier (Berthier and others, 2010). The surface velocity field can be derived from various sources such as photogrammetry or radar-based measurements, but it must cover the region of the glacier over which we wish to calculate ice thickness. If the available velocity fields cover the entire glacier, the fact that H → 0 at the boundary of the glacier is a sufficient initial estimate of the ice thickness at the boundary of the domain of interest; otherwise, some prior knowledge of the ice thickness is necessary. The next section details the datasets used for this study, and how they are used to solve Eqn (6).

5. Application to Columbia Glacier

With the fjord bathymetry and surface DEMs, we have ice thicknesses for years when Columbia Glacier extended into the proglacial fjord, and thus have ice thicknesses at the downstream boundary of the flowlines. Because of the glacier’s continued retreat since bathymetry measurements were last taken in 2005, the full-coverage velocity maps in 2011 do not overlap with the region of measured bed topography. To calculate the bed topography using the 2011 velocity fields, we must first calculate the bed topography in the region between the proglacial fjord and the current glacier (‘Confluence’, Fig. 1). The velocity maps from the mid-2000s are then used to calculate the bed topography well upstream of the region of the 2011 terminus, and the remainder of the bed topography is calculated using the 2011 velocity fields.

We selected velocity fields from 1984, 1985, 2005, 2006 and 2009, based on measurement density and time of year (Table 1). Each velocity field is associated with a DEM and a modeled profile of surface mass balance. Thinning rates for each of these datasets are calculated by differencing two subsequent DEMs.

We calculated ice thicknesses upstream of the region where bathymetry is known (Fig. 1) in the following manner. To avoid interpolating velocities in the region near the calving front, we selected upstream and downstream boundaries for the flowlines; these boundaries are never within 3 km of the calving front for any particular velocity field. Flowlines are calculated using MATLAB’s built-in stream2 function (MATLAB ©1984–2012 Mathworks, Inc.) and are recalculated for each velocity field. The resulting flowline vertices define cell boundaries for the thickness calculation. It should be noted that these cell boundaries do not have constant area or dimension: their size changes according to the flow field (Fig. 2). At the downstream edge of this flowband, we calculated ice thickness using the known surface elevation and the bathymetry We then calculated the ice thickness at the upstream boundary of each cell, according to Eqn (6), using both rates of surface mass balance and rates of surface elevation change.

Before retreat began at Columbia Glacier, the majority (>90%) of the ice motion was calculated to be sliding (Rasmussen, 1988), and this has held throughout the retreat. Assuming Glen’s flow law with exponent n = 3 and 90% of the motion due to sliding gives γ = 0.98 (Cuffey and Paterson, 2010). Based on this, and the results of the sensitivity test in the next section, we set γ = 1, which corresponds to the case where vertical shear is negligible or absent. We also set the initial flowline spacing to 250 m, as this allows for faster calculation.

We repeated the calculation of H over each of the available partial velocity field sets to provide redundancy in the calculation, and calculated bed elevation by subtracting the ice thickness from the surface elevation. We then gridded the calculated bed elevations to a 100 m grid by taking the arithmetic mean of calculated values that fall within each gridcell. On average, the spread between individual values for a gridcell is ~13 m, which is well within the uncertainties associated with the input data. This gridded topography provides a map of bed elevation that covers the region of the confluence, as well as part of the west branch, and the eastern trunk (Fig. 1).

This incomplete bed topography map is then used as input with the 2011 velocity fields, and the bed topography calculation is thereby extended to the entire glacier extent. Because a full-coverage surface DEM of Columbia Glacier that matches with the 2011 TerraSAR-X velocity fields is not currently available, we use previous full-coverage DEMs to estimate surface elevations in areas not covered by the 2010 partial DEM as follows. First, thinning rates between 1957 and 2007 are calculated using the two full-coverage DEMs. A full-coverage map of relative thinning rates, A, is produced by dividing the map of these thinning rates by their mean value over the same area covered by the 2010 DEM. Next, annual thinning rates between 2007 and 2010 are calculated, using the 2010 DEM, and are assumed to be valid for 2011. The mean of these values is calculated, and the final surface elevation map is calculated:

(7)

where Z is the surface elevation for the subscripted year, is the mean value of the 2007–10 thinning rates evaluated over the 2010 domain, and A is the ratio between the 1957–2007 thinning rates and the spatial mean. This estimation can add considerable uncertainty into the ice thickness calculation, and this is discussed in Section 6.2.

6. Results

General information about the calculated topography map is given in Table 2, and the final topography map itself is shown in Figure 3. Bed topography data (netCDF format) are available as supplemental material to this paper (http://www.igsoc.org/hyperlink/11j249.html).

Table 2. General results for the new bed topography map. S indicates glacier map area, V total ice volume, ‘V below s.I.’ is volume of ice that is below sea level, H _max is maximum ice thickness, is mean ice thickness (± one std dev.), H_med is median ice thickness and z_bed is bed elevation

Fig. 3. Calculated bed topography map for Columbia Glacier, 1957 extent. May 2011 terminus location shown as a red dashed line. Location of thickest ice in both 1957 and 2007 is shown as a red star.

6.1. Sensitivity analysis

We conduct analysis of the sensitivity of calculated ice thickness to the main components of the model-initial flowline separation (i.e. the starting flowband width), γ, surface mass-balance rates, rates of surface elevation change and the velocity direction – using the four velocity fields from 1984 and 1985. The sensitivity analysis is limited to the area downstream of the 2004 terminus (i.e. where the bathymetry is known), as the density of bed topography measurements there allows for a more complete comparison of ice thicknesses. The soft-rock geology that characterizes the Chugach Range (e.g. Wilson and others, 1998) garners expectations of rapid erosion, with the implication that bathymetric measurements do not perfectly represent the glacier bed when ice was present. Whether these expected high rates of erosion directly result in sedimentation in the fjord is not known; for simplicity, we assume that the measured bathymetry represents the glacier bed.

First, we investigate the effect of changing initial separation of flowlines on calculated ice thickness. Because the most dense velocity fields (derived from TerraSAR-X) are available in a 100 m grid, this is taken as a lower limit for the flowline separation. Initial flowline separation is then varied, in 50 m increments, up to 500 m. Results are shown in Figure 4a. Overall, flowline separation shows very little effect on either mean absolute difference or RMSE. As flowline separation increases beyond 300 m, however, maximum error increases significantly, suggesting that the limit where is flowline separations of ~300m.

Fig. 4. Results of analysis of calculated ice thickness sensitivity to changes in (a) initial flowline offset, and (b) γ (Eqns (3–6)). Error is defined as the difference between measured and calculated ice thicknesses.

Second, we investigate the sensitivity of the calculated bed topography to values of γ by varying γ from 0.8 to 1 (Fig. 4b). RMSE for the given values of γ are all within a 1.5 m range, indicating that values of γ are not as important as the initial flowline separation. For a typical gridcell in this region (covered by the bathymetry data; Fig. 1), typical values are about —0.5 m w.e. a^-1 for , 4.5 m w.e. a^-1 for ∂H/∂t and 2 x 10⁸m³a^-1 for q_in and q_out. Given this, calculated ice thicknesses are on the order of 500 m; for this value, varying γ from 0.8 to 1.0 changes the calculated ice thickness by less than 1 m. Further analysis of Equation (6) shows that even awayfrom the fast-flowing portion of the glacier, changes in γ introduce changes in ice thickness that are much smaller than the uncertainties associated with the components of Eqn (6).

Third and fourth, we investigate the effect of varying surface mass-balance rates by increasing and decreasing by 0.5 and 1.0ma^-1, and by increasing and decreasing values of ∂H/∂t by 1.0 and 2.0 m a^-1. Results are summarized in Table 3. Very little change in error is observed. Again, given typical values of components of Eqn (6), values of and ∂H/∂t are typically one to four orders of magnitude less than ice velocities, and ice fluxes are typically several orders of magnitude larger than the terms on the right-hand side of Eqn (2). This is a reflection of the fact that for tidewater glaciers in rapid retreat like Columbia Glacier, climatic balances are dwarfed by dynamic changes (Meier and Post, 1987; Pfeffer, 2007), at least on the lower region of the glacier.

Table 3. Results of analysis of calculated ice thickness sensitivity to changes in surface mass balance (), surface elevation change (∂H/∂t) and velocity vector direction (a), over the domain covered by measured bathymetry (Fig. 1). Here error is defined as the difference between measured and calculated ice thickness. Mean error is reported as arithmetic mean ± one std dev.

Finally, we investigate the effect of perturbing the input velocity fields. To this end, we perturb the angle of the velocity vectors using a uniform random distribution on the interval [-Δa, Δa], where Δa is measured in degrees. We use values of 0.5,1,2,5 and 10 for Δa. Results are summarized in Table 3. For large (>5°) values of Δa, maximum error increases significantly (by a factor of 2 over unperturbed values), but mean error and RMSE remain relatively unchanged. Most likely, this is due to the random nature of the perturbations; with equal numbers of positive and negative changes, we would expect the overall effect on the mean to be small.

6.2. Uncertainty and error analysis

Because full-coverage DEMs are not available for the same time as full-coverage surface velocity maps, surface elevations must be estimated using the existing full- and partial-coverage DEMs, as shown by Eqn (7). To estimate the uncertainty introduced by this method, we calculate a full-coverage DEM for 2010, and compare the estimated to the measured surface elevations. The resulting RMSE in surface elevations is 26.6m.

For the full-and partial-coverage DEMs, we estimate uncertainties by comparing elevations in non-glacierized areas, resulting in an estimated vertical position error of 2m. With a mean time separation of 1 year between successive DEMs, this results in an uncertainty of 2.8 m a^-1 in calculated rates of surface elevation change.

For velocity datasets derived from non-digital photogrammetry (dates before 2004), Krimmel (2001) estimates an uncertainty in displacement of 4 m between successive images. Given a typical time separation between successive images of 0.115 years, this translates to an uncertainty in velocity of ~35 m a^-1.

For velocity datasets derived from digital photogrammetry (2004–10), uncertainties arise from point identification on photographs and interpolation of irregularly spaced data to grid nodes. We use the following equation to estimate the uncertainty (ma^-1) in velocities derived from aerial photographs:

(8)

where C is uncertainty in image registration and feature tracking in pixels (p), Δx is the image resolution (m p^-1), and Δt is the time separation between successive images (days). Using typical values of 1–2 p for C, 2mp^-1 for Δx and 40 days for Δt, we estimate an uncertainty of 18–36ma^-1 in velocity values. Given this, we use the upper bound on the uncertainty (36 m a^-1), to encompass both digital and non-digital datasets.

We estimate the uncertainty introduced by interpolating irregularly spaced velocities to grid nodes by calculating the standard deviation of the difference between interpolated and uninterpolated values at the same point. The resulting interpolation error is ~13 m a^-1.

As reported by Rasmussen and others (2011), the RMSE in values of used is 1.0 m w.e. a^-1. The total uncertainty in calculated ice thickness, assuming that each of the estimated uncertainties (DEM, thinning rates, velocity and ) is independent and quadratically additive, is 46.7 m.

This treatment is for random uncertainty only, and does not address systematic error. To ensure that systematic errors are not present in values of and ∂H/∂t, we check the value of — ∂H/∂t, integrated over the 2007 glacier domain, obtaining a value of 6.5 km³ a^-1. This is comparable with estimates of calving fluxes for Columbia Glacier (O’Neel and others, 2005; O’Neel, 2012), lending confidence that the errors in measurement are not systematic.

6.3. Ice thickness and bed topography

We calculate ice thickness distribution for the glacier for the two years in which we have full-coverage DEMs, 1957 and 2007 (Fig. 5). Over that period, the mean thickness (± one std dev.) of the glacier is nearly halved, from 280 ± 215 m in 1957 to 145 ± 140 m in 2007. The median thickness of the glacier is also more than halved, from 225 m in 1957 to 100 m in 2007.

Fig. 5. Calculated ice thickness map (a) 1957 and (b) 2007, and (c) thickness decrease 1957–2007 for Columbia Glacier.

The maximum thickness of the ice, however, is nearly the same in the two years: 1040 m in 1957, and 1005 m in 2007. The location of thickest ice for both years is indicated in Figure 3, and is very high on the glacier, where very little thickness change has occurred.

The calculated bed topography, over the 1957 glacier extent, ranges from 525 m below sea level at its deepest point (ξ ≈ 52 km; Fig. 1), to 3670 m a.s.l. at its highest point (ξ ≈ 0 km). The deepest part of the bed covered by the current (2011) glacier extent is, by contrast, 280 m below sea level (ξ ≈ 41.5). At the June 2011 location of the terminus (ξ ≈ 48.2 km), the maximum calculated depth is 154 m below sea level.

The portions of the bedrock topography lying below sea level are 11%, 6% and 4% for the glacier extents in 1957, 2007 and 2011, respectively. In the center-line coordinate system, the calculated bed rises above sea level at ξ ≈ 36.4 km, which is in good agreement with previous estimates (Mayo and others, 1979). This is also 11.7 km from the current (2011) terminus.

We also compare calculated ice thicknesses to ice thicknesses derived from radar measurements. A comparison of individual radar tracks is shown in Figure 6, and comparison of all measured and calculated ice thicknesses is shown in Figure 7. Based on comparison to radar, the method appears to work best in areas with many overlapping velocity datasets, like the eastern trunk (radar tracks 3–5, Figs 1 and 6). Only one radar track (track 3, Figs 1 and 6) is crossflow, at least for part of its extent. For the part of the track that is cross-flow (distance along track 6–12 km; track 3, Fig. 6), the difference between the calculated and measured ice thickness is as small as it is for the along-flow portion of the track, lending confidence that the method succeeds for both cross-flow and along-flow directions.

Fig. 6. Comparison of calculated thickness (dash-dot) to measured thickness (solid line) along radar profiles, oriented west to east (Fig. 1). Track number is indicated in upper left corner of each panel. Difference between measured and calculated ice thickness is expressed as mean value ± one std dev., indicated in upper right corner of each panel.

Fig. 7. Comparison of calculated and measured ice thickness for all radar points (Fig. 1). Legend indicates to which radar track each point belongs. The statistics at the bottom right refer to the ensemble of points (n: number of points; avg dev: average deviation; RMSE: root-mean-square error).

6.4. Volume change

We calculate the total ice volume for Columbia Glacier in 1957 and 2007 using the ice thickness maps (Fig. 5). For 1957, we calculate a total ice volume of 294 km³, and 134 km³ for 2007. The calculated volume loss over the period 1957–2007 is then 160 km³. The resulting percentage volume loss of 54% indicates that Columbia Glacier has lost over half of its volume since 1957, most of which occurred following the onset of retreat, around 1982. Assuming a density of 900 kg m^-3, the calculated mass loss is 144 Gt.

For the 1957 glacier extent, the calculated ice volume below sea level is 18 km³ (6% of the total; cf. 11% of the glacier area below sea level), and for the 2007 glacier extent the calculated ice volume below sea level was 7 km³ (5%; cf. 6%).

Figure 8 illustrates the existence of a clear separation between the upper (ξ < 26.5 km) and lower (ξ > 26.5 km) glacier. This separation is also demonstrated in Figure 9a. A majority (91%) of the volume change over the period 1957–2007 has occurred at 1957 surface elevations below 1400 ma.s.l. Despite making up a substantial portion of the total 1957 area (34%), elevations above 1400 m contributed little (9%) to the total volume change over the period.

Fig. 8. (a) Calculated center-line bed topography and surface elevation in 1957 and 2007, with location of apparent separation between upper and lower glacier regions indicated. (b) Calculated ice thickness along center line in 1957 and 2007, along with 2011 center-line surface velocity, showing jump in velocity near location of apparent hinge point between upper and lower glacier regions.

Fig. 9. (a) Percent of total volume change 1957–2007, along with percent of total area (1957 hypsometry), and (b) thickness change 1957–2007 for each 100 m elevation band (from 0 to 3700 m; 1957 hypsometry).

7. Discussion

Our method to compute spatially distributed glacier thickness is based on evaluation of the mass continuity equation between adjacent flowlines that are interpolated from surface velocity fields. The main advantage this method has over other methods for calculating ice thickness is its assimilation of available data, thereby giving sound constraints on the calculated ice thickness. Because the method does not calculate the ice thickness on the whole glacier domain at once, it is likely faster than other methods that solve the mass continuity equation on a grid.

The requirement of accurate and dense surface velocity fields is both an advantage and a disadvantage of this method. The method can only be applied in regions where ice is flowing and surface velocity fields are available. It is also limited in application to one glacier at a time, and significant work is required to estimate the total ice volume of an entire region. For individual glacier applications, however, it is both accurate and easily implemented. At present, there are no other tidewater glaciers in the world with the wealth of data available at Columbia Glacier, but increased availability and application of remote-sensing data, such as TerraSAR-X and TanDEM-X imagery, will serve to reduce the input data barrier.

We tested sensitivity of calculated ice thickness to different components of Eqn (6) in the early phase of Columbia Glacier’s retreat (1984 and 1985). Because typical values of each component are small compared with typical values of ice flux at that time, we find that calculated ice thickness is not sensitive to relatively large changes in the value of each component. It should be stressed that this is most likely not typical for alpine glaciers. It is also likely not true for the upper region of the glacier, where dynamic thinning has been largely non-existent, and so errors in the other inputs (, velocity) are large relative to the ice thickness. Despite this, the method performs well on the upper region of the glacier, indicating its ability to calculate ice thickness in both fast-and slow-flowing regions of a glacier.

We have tested the sensitivity of the calculated ice thickness to uniform values of γ only, but the possibility remains that values of γ may not be uniform at the upstream and downstream boundaries of each cell: that is, Eqn (6) can be reduced to (with γR, γp the values of γ at the upstream and downstream boundaries, respectively)

(9)

where H_R and H_P are the ice thicknesses at the upstream and downstream boundaries, respectively. If the ratio of γ_P to γ_R reaches 0.8, then the error introduced by assuming a uniform γ could reach 20%. In practice, however, the existence of such a large change in the value of γ over such a small distance (<200–300m) seems unlikely. If larger cell sizes are used, however, it may become necessary to consider nonuniform values of γ.

A comparison between the newly presented bed topography map and the previous (Engel, 2008) bed map shows that the previous map overestimated the ice thickness by as much as 70%, a fact that is acknowledged by the author Results of this comparison are summarized in Table 4. There may be several reasons for this discrepancy between the two calculated bed topographies. First, Engel (2008) noted that spurious overdeepenings in calculated bed topography may be introduced when flowlines originating near the terminus are used. To avoid this problem, we do not calculate flowlines in the near-terminus region (within ~3 km of the terminus for any particular velocity field). Second, we use many overlapping velocity fields that were unavailable at the time of the previous study Third, the previous study assumed that surface mass-balance rates were zero and used thinning rates calculated from center-line laser altimetry profiles in the mid-1990s to estimate thinning rates in both the 1980s and in 2004. This assumption would most likely result in an overestimation of ice thickness, but a quantification of the extent is not presently available.

Table 4. Comparison of new bed topography map (Fig. 3) and topography map produced by Engel (2008). The new map is evaluated only over the domain covered by the previous map (Fig. 1). ‘% below s.l.’ indicates the percentage of the bed that is below sea level

We also compare the bed topography calculated using our method with the estimate by Rasmussen (1988). The estimate computed by Rasmussen covers only the lowest 15 km of the 1957 glacier extent, which is approximately the limit of the bathymetry measurements. Comparing Rasmussen’s estimate with the measured bathymetry yields a mean difference of 42 ± 130 m, while the mean difference between our calculated bed topography and the measured bathymetry is 8 ± 43 m.

With the full bed topography map, we calculated total ice volume for the years when full-coverage surface DEMs are available. As far as we know, this is the first estimate of total ice volume for Columbia Glacier. It is not, however, the first estimate of volume change at Columbia Glacier. Using thinning rates derived from laser altimetry, Arendt and others (2006) calculate a total mass loss over the period 1957–2004 of 141 Gt (3.01 Gta^-1). Arendt and others do not take into account ice below sea level. Over the period 1962–2006, Berthier and others (2010) estimated the rate of volume loss above sea level at Columbia Glacier to be 2.43 Gta^-1. Our estimate of 144Gt (2.88 Gta^-1; 2.68Gt a^-1 above sea level) during the period 1957–2007 is between these two reported estimates.

Figures 8 and 9 illustrate a decoupling between the upper and lower regions of the glacier at ξ ≈ 26.5; this separation is apparent in both the surface elevation and ice thickness profiles as a ‘hinge point’. Above this hinge point, the glacier surface is largely the same as it was in 1957; below this hinge point, the glacier has thinned in excess of 500 m in places. The calculated bed topography in the region of the hinge point shows several prominent bedrock humps along the center line, where the ice thickness drops from 300–400 m to 100–200 m, approximately a 50% reduction in ice thickness over a length scale of 1–2 km. These rapid drops in ice thickness likely serve to impede stress transmission between the upper and lower regions of the glacier, and retard the drawdown of ice to the lower glacier region.

There are abundant studies that calculate glacier volume change over recent time periods, but very few that calculate initial volumes, and that therefore can report relative volume changes. In a study estimating the glacier ice volume throughout the Swiss Alps, Farinotti and others (2009b) found that the ice volume of the Swiss Alps decreased by 12% over the period 1999–2008. In a second, smaller study of 20 glaciers in the southeastern Swiss Alps, Huss and others (2010) find a regional decrease of 47% of ice volume over the period 1900–2008. The authors determine total ice volume using the method developed by Farinotti and others (2009a). For individual glaciers in the region, percentage volume loss ranged from 30% to 75%. For the largest glaciers in the region (with surface areas of 7–1 7 km²), volume loss ranged from 38% to 62%. These mostly land-terminating glaciers exhibit the same relative changes in volume as Columbia Glacier, but over a time period that is over twice as long. Considering that volume change at Columbia Glacier was very nearly zero over the period 1957–81 (Meier and Post, 1987; Rasmussen and others, 2011), Columbia Glacier has lost a comparable percentage of volume in one-quarter of the time taken by these glaciers (26 years). This discrepancy is not surprising, given that most of the volume change from Columbia Glacier is due to dynamic causes such as calving, whereas volume loss from glaciers in the Swiss Alps is due to more negative surface mass balances alone.

Based on our findings at Columbia Glacier, we suggest that the potential for a rapid acceleration of glacier discharge does not require much of the glacier to be grounded below sea level; indeed, Columbia Glacier lost 50% of its ice volume between 1957 and 2007, despite having only 11% of the bed topography over the 1957 glacier extent below sea level. The Greenland ice sheet has a relatively large portion of its ice grounded below sea level (e.g. Bamber and others, 2001) and is drained by many outlet glaciers, some of which have deep troughs that extend into the interior of the ice sheet (e.g. Pfeffer and others, 2008). Recent studies of the outlet glaciers of the Greenland ice sheet (e.g. Howat and others, 2005; Rignot and Kanagaratnam, 2006; Moon and Joughin, 2008) have shown an increase in both ice velocity and ice discharge, in a relatively short period of time. In light of the larger portion of the Greenland ice sheet being grounded below sea level, we might expect the recent, rapid adjustments in outlet glacier geometry to continue.

In addition to the potential implications for the Greenland ice sheet, our finding of a large, rapid, tidewater response with a small portion of the glacier grounded below sea level has implications for other mountain glaciers and ice caps. Tidewater glaciers are prevalent in many mountain glacier and ice-cap systems around the world, in particular in the Canadian and Russian Arctic, Svalbard and the periphery of Greenland and Antarctica. Despite their importance for sea-level rise (e.g. Meier and others, 2007; Hock and others, 2009), relatively little is known about the dynamics and stability of these systems.

Pfeffer (2007) proposed a simple mechanism to initiate irreversible tidewater retreat. By examining how waves of thickening and thinning propagate through a glacier, he found that instability arises when changes in glacier geometry serve to reduce resistive stresses more than driving stresses, and that stabilization can occur where there is a ratio of ice thickness to water depth h/hw above 1.49, a negative longitudinal flux gradient (∂q/∂x) and high rates of diffusive thickening D. Applying this analysis to the 2011 data for Columbia Glacier, and averaging values over the 2.5 km nearest the terminus, we find an average value of h/h_w = 1.7, ∂q/∂x that is positive but nearly zero, and a rate of diffusion (D/W = 2.64x 10⁸ m² a^-1, expressed per unit width) in the same order of magnitude as Pfeffer’s values for Columbia Glacier (D/W = 8.33 x 10⁸ m² a^-1). In addition, the non-diffusive propagation speed is positive, indicating that perturbations propagate downstream, which in turn indicates stability.

We can also examine the slope of the bedrock along the center line at Columbia Glacier. Weertman (1974) calculated that an ice sheet whose bed slopes inward toward the center is unstable, a finding supported by numerical experiments of tidewater glacier evolution (Vieli and others, 2001). At present, the bed slope at the terminus of Columbia Glacier slopes inward; it stops sloping inward at ξ = 41.5 km, 6 km from the present terminus location. Based on this, and the application of Pfeffer’s (2007) analysis, it is possible that the rapid part of Columbia Glacier’s retreat is nearing an end; at the very least, it does not seem likely that the retreat rate will match that seen in the 1980s and 1990s.

Because of its high calving rates, thinning rates at Columbia Glacier are not presently representative of most glaciers in Alaska. However, the well-studied retreat of Columbia Glacier could serve as an analogue to the opening of both Icy Bay and Glacier Bay, Alaska. The glaciers of Icy Bay retreated ~40 km between the early 20th century and the present (Barclay and others, 2006), while the glaciers of Glacier Bay retreated >100 km between the late 18th and early 20th centuries (e.g. Larsen and others, 2005; Molnia, 2008; Barclay and others, 2009). Because Columbia Glacier has been studied since before its retreat, it can provide some insight into how the retreats of the glaciers of Icy Bay and Glacier Bay evolved. Because the glaciers of Icy Bay and and Glacier Bay have largely slowed or finished their retreats (e.g. Porter, 1989; Molnia, 2008; Barclay and others, 2009), they can provide insight into how Columbia Glacier will continue to evolve, until its retreat eventually ends.

8. Conclusion

We have presented a method, based on conservation of mass, for estimating spatially distributed glacier ice thickness and bed topography. The method requires velocity data of sufficient density to ensure that over the grid spacing used, ice thickness and surface velocity can be accurately interpolated. The other data requirements (multiple full-coverage DEMs, and surface mass-balance rates) are large, but not necessarily prohibitive. Increased availability of remote-sensing data products, such as velocity fields from TerraSAR-X, and surface DEMs from TanDEM-X, will likely serve to further reduce this barrier. Because the method is applicable to both fast and slow-moving ice, it can potentially be applied to any glacier for which sufficient data are available.

Using this method, we have presented a high-resolution bed topography map for Columbia Glacier. The new topography map covers the entire glacier extent, both current and former, and is in excellent agreement with radar ice thickness measurements, with a mean difference between measured and calculated ice thickness of –5 m, and RMSE of 44 m. The majority of this error is likely due to our approximation of the current (2011) glacier surface elevation using DEMs extrapolated from previous years (2007 and 2010).

With this bed topography map, we have calculated the total ice volume at Columbia Glacier in 1957 and 2007, finding that the glacier has lost > 50% of its volume during this period. The majority of the calculated volume change is from areas where the 1957 surface is below 1400 m a.s.l. This is most likely due to the difficulty of transferring stress across several prominent bedrock humps that occur along the center line, ~26km from the head of the glacier.

The small fraction of the 1957 bed extent that is below sea level shows that having a large percentage of the bed below sea level is not necessary to initiate a rapid tidewater response. Because little is known about the stability and dynamics of tidewater glacier systems around the world, the implications of this finding for these other systems are unknown. Application of an analysis of tidewater stability to the present state of Columbia Glacier, as well as examination of its bedrock slope, indicate that the rapid portion of Columbia Glacier’s retreat is perhaps nearing an end.

Author Contribution Statement

R.McN. developed the method, performed all the calculations and wrote the paper. R.H. provided guidance and feedback throughout the work and contributed to the writing. S.O’N. initiated the study, M.T. suggested the general path of the approach, and both together with L.A.R., W.T.P., M.T. and B.E.S. provided comments on the method and the paper. L.A.R. provided the surface mass-balance data, Y.A. the photogrammetric velocity data, M.B. and I.J. the satellite- derived velocity data, B.E.S. and H.C. the radar data, and S.H. digitized the glacier outline.