The Cryosphere Surge dynamics on Bering Glacier , Alaska , in 2008 – 2011

A surge cycle of the Bering Glacier system, Alaska, is examined using observations of surface velocity obtained using synthetic aperture radar (SAR) offset tracking, and elevation data obtained from the University of Alaska Fairbanks LiDAR altimetry program. After 13 yr of quiescence, the Bering Glacier system began to surge in May 2008 and had two stages of accelerated flow. During the first stage, flow accelerated progressively for at least 10 months and reached peak observed velocities of ∼ 7 m d−1. The second stage likely began in 2010. By 2011 velocities exceeded 9 m d −1 or ∼ 18 times quiescent velocities. Fast flow continued into July 2011. Surface morphology indicated slowing by fall 2011; however, it is not entirely clear if the surge is yet over. The quiescent phase was characterized by small-scale acceleration events that increased driving stresses up to 70 %. When the surge initiated, synchronous acceleration occurred throughout much of the glacier length. Results suggest that downstream propagation of the surge is closely linked to the evolution of the driving stress during the surge, because driving stress appears to be tied to the amount of resistive stress provided by the bed. In contrast, upstream acceleration and upstream surge propagation is not dependent on driving stress evolution.


Introduction
Glacier surging is a unique glacier dynamic behavior, in which glacier flow speeds oscillate between two phases: a quiescent phase, characterized by slow flow that steepens glacier geometry, and a surge phase, characterized by extremely fast flow that flattens the geometry.Surge events exhibit flow speeds 10-100 times quiescent flow; they are relatively short, lasting from months to years, and can initiate and terminate rapidly (Cuffey and Paterson, 2010).The quiescent phase lasts decades, over which the glacier develops a steeper geometry that triggers another surge event.The time required for the glacier geometry to steepen during quiescence dictates the duration of the quiescent phase and the surge cycle overall (Meier and Post, 1969;Raymond, 1987;Harrison and Post, 2003).Though there have been many studies focused on glacier surge dynamics, most of what is known about surging comes from observations on a few small glaciers over only 1-2 surge cycles (Meier and Post, 1969;Raymond, 1987;Raymond and Harrison, 1988;Bindschadler, 1982;Heinrichs et al., 1996).
Quiescent flow velocities are, on average, slower than balance velocities; thus surface mass balance causes thickening in the accumulation zone, thinning in the ablation zone, and a steepening glacier geometry overall.Such an evolution has been closely observed on Variegated Glacier during quiescent phase (Raymond and Harrison, 1988).However, during quiescence, some surge-type glaciers have small acceleration events that redistribute thickening and thinning along the glacier profile (Meier and Post, 1969;Raymond, 1987;Harrison and Post, 2003;Raymond and Harrison, 1988;Heinrichs et al., 1996).Ultimately, the quiescent phase develops a region of thickening called the reservoir zone, which forms upstream of a thinned receiving zone.The surge phase theoretically reverses this process, thickening the receiving zone and thinning the reservoir zone and returns the glacier geometry to where it was at the end of the previous surge.The boundary between the reservoir and receiving zones is termed the dynamic balance line (DBL) (Raymond, 1987).

E. W. Burgess et al.: Surge dynamics on Bering Glacier
The location of the DBL is no coincidence.The DBL is a surface expression of key changes in bed conditions that ultimately control the surge onset and surge progression.At the DBL, upstream thickening and downstream thinning steepen the glacier, which increases local driving stresses.This increase in local driving stress creates a trigger point for surge initiation by changing fundamental dynamics in the basal hydrologic system.
For most glaciers in Alaska, warm summer temperatures supply the bed with vast amounts of water.During quiescent phases, the basal hydrologic system has the ability to generate a channelized drainage system that can evacuate water inputs efficiently and prevent hydrologic pressurization of the bed (Röthlisberger, 1972).If water inputs increase rapidly, such as in spring or during rain events, the bed will pressurize, but only temporarily, while the channels expand to accommodate the additional flux.
The surge phase is thought to begin when the increased driving stress at the DBL prevents the glacier from attaining an efficient channelized drainage system.Rather, the high driving stresses promote and maintain a distributed subglacial hydrologic system that can only increase flux by increasing basal hydrologic pressure (Kamb, 1987).As a result, basal water pressures remain extremely high, sometimes within a few bars of ice overburden pressure (Raymond, 1987).The high water pressure reduces the amount of shear stress the bed can support, and the glacier accelerates until it is able to attain force balance.The stability of the distributed drainage system is thus a key parameter.In a distributed system model, Kamb (1987) defines such a variable called the melt-stability parameter that is dependent on the basal shear stress and the smoothness of the bed.
The outstanding questions in understanding surge processes involve understanding exactly how the increased driving stress causes a discrete shift from a basal hydrologic system capable of channelization (quiescent phase), to a resilient distributed system (surge phase) before abruptly reverting back.The problem is we know little about actual conditions at the bed as they are very difficult to observe.Most surgetype glaciers are believed to have deformable sediment beds (Raymond, 1987) that can deform rapidly once the driving stress exceeds the yield stress of the sediment.Once the bed begins to deform, faster sliding velocities have little effect on the basal shear stress; therefore, the glacier can accelerate until lateral or longitudinal resistive stresses are able to attain force balance (Cuffey and Paterson, 2010).A key remaining question is how distributed/channelized systems evolve in sub-glacial systems with deformable or mixed beds.
Another outstanding problem is surges have been found to propagate downstream and upstream from their trigger point (Raymond, 1987).This propagation presumably is a consequence of reduced basal shear stress requiring increases in longitudinal or transverse stresses to maintain force balance.But exactly how this propagation evolves at the bed is not well understood.
The primary reason why these questions still stand is glacier beds are nearly impossible to observe over areas larger than a borehole.Thus, understanding how these processes evolve over large areas will require us to use readily observable surface expressions of bed dynamics to infer bed conditions.Here we examine the most recent surge cycle on the Bering Glacier system that initiated in May 2008 and is believed to have terminated in Summer 2011.We use altimetry data acquired between 1995 and 2011 and surface velocity data from 2007 to 2011 to examine the evolution of flow speed and driving stress during quiescent and surge phases.

The Bering Glacier system (BGS)
The BGS is the largest surging glacier outside of the ice sheets (Molnia, 2008), covering 4373 km 2 (Beedle et al., 2008) and accounting for 4 % of the ice area in Alaska (Beedle et al., 2008;Berthier et al., 2010) and 6 % of the mass loss (Arendt et al., 2002).It extends from ∼ 100 m to 3000 m elevation with an equilibrium line at approximately 1000 m (Molnia, 2008).The BGS has a broad piedmont lobe that calves into Vitus Lake to the south and abuts the Stellar lobe to the west.It is situated in coastal southern Alaska, in a maritime climate, and has a high rate of mass turnover.
The geography and nomenclature of the BGS is rather complex (Fig. 1).The BGS drains most of the westwardflowing Bagley Ice Valley (BIV) and eastward-flowing West Bagley (WB) (Fig. 1).Ice from the BIV and WB converges at 30 km on the WB profile in Fig. 1.There, the majority of the ice diverges south to form Bering Glacier.Ice on the northern edge of the WB and BIV is diverted north to form the smaller Tana Glacier -this ice is not considered part of the BGS.The Jefferies and Quintino Sella glaciers are tributaries of the BIV.The majority of the Jefferies ice flows north into the Tana and thus is not part of the BGS, while all of the Quintino Sella ice flows into the Bering and is part of the BGS.We herein refer to Bering Glacier as the portion of the BGS profile in Fig. 1 from 85 km to the terminus and the BIV as the portion from 0-85 km.Otherwise, references to locations will refer to kilometer distances along the BGS or WB profiles as indicated in Fig. 1.
The 1993-1995 event is the only BGS surge with quantified flow velocities; these results are worthy of summary to place the current surge in context.Surge onset occurred at approximately BGS-135 (Roush, 2003) (Fig. 1).The fastest velocity recorded was 59 m d −1 near the terminus (Roush, 2003).More generally, flow velocities were 10-20 m d −1 near the terminus (Roush, 2003) and up to 5 m d −1 in the lower BIV (Fatland and Lingle, 2002).Both of these areas had velocities < 1 m d −1 during quiescence (Post, 1972;Fatland and Lingle, 2002).Velocity data were unable to confirm the existence of an organized surge front propagating upstream or downstream from the onset location.A leading surface undulation propagated downstream from BGS-135 to the terminus at ∼ 100 m d −1 (Roush et al., 2003).Upstream propagation was less clear.Fatland and Lingle (1998) found a delay between surge onset on the Bering (∼ BGS-135 km) and a subsequent acceleration of the WB (Fig. 1), which lies ∼ 60 km upstream.Assuming that the delay was due an upstream propagation front originating from the surge onset location, Fatland and Lingle (1998) estimated the upstream propagation velocity was 200-500 m d −1 .Fatland and Lingle (1998) also considered the possibility that the WB accelerated because of a linked sub-glacial hydrologic system and may not have been a classic propagating surge front as observed on the Variegated Glacier (Kamb et al., 1985;Raymond, 1987).
During the 1993-1995 surge, Fatland and Lingle (2002) found a distinct longitudinal change in ice dynamical behavior at a point just downstream of the Jefferies Glacier confluence, at BGS-45 (Fig. 1).Upstream this point in 1994, acceleration rates were temporally uniform and velocities were relatively similar to quiescent velocities.Downstream of this point, acceleration rates fluctuated rapidly (on 3day timescales) and velocities increased dramatically in the downstream direction where uniform velocities existed prior to the surge.Short-term surface elevation changes, possibly indicative of transient sub-glacial water, were also more abundant downstream of this point during the 1993-1995 surge (Fatland and Lingle, 2002).

SAR offset tracking
We use L-, C-and X-band SAR platforms to generate ice displacement fields via offset/speckle tracking methods (Strozzi et al., 2002;Gray et al., 1998;Michel and Rignot, 1999).We use ALOS PALSAR Fine-Beam (46-day repeat), RADARSAT Fine-Beam (24-day repeat) ERS Ice Phases (3-day repeat) and TerraSAR-X StripMap (11-day repeat) data.Acquisitions are screened to obtain pairs with temporal baselines of 1 or 2 orbit intervals and perpendicular baselines < 400 m for RADARSAT and ERS, < 1000 m for PAL-SAR and < 20 m for TerraSAR-X.We obtained a total of 77 frames, providing 40 pairs acquired between 2006 and 2010.For the final analysis, we use 21 PALSAR pairs, one TerraSAR-X pair acquired in 2011, and 2 ERS pairs acquired during the 1993-1995 surge (same data as used by Fatland andLingle, 1998, 2002).The dates for these pairs are indicated on the timeline in Fig. 2. Pairs could not be obtained at regular intervals due to poor data availability.Summer pairs were generally not utilized because SAR offset tracking was found to be ineffective during extensive melt, though two summer PALSAR pairs did yield usable displacement fields.The single July TerraSAR-X pair, with an 11-day repeat, also produced good results.
Single-look complex (SLC) pairs are co-registered using spatial domain, normalized cross-correlation optimization offset tracking within GAMMA ® software ( , 2002).This method quantifies the displacement (offset) of pixel patterns between the two images.For each image pair, we derived offsets on stable ground adjacent to glaciers and assumed the offset vector field represents image co-registration.A polynomial function was fit to this offset field and then subtracted from the offsets observed on glaciers in the same image pair.The BGS and surrounding ice is extensive enough to cover most of a single SAR frame (∼ 70 km × 80 km); in these cases, glaciated ridgelines were allowed into the image co-registration in order to obtain an accurate offset polynomial.Related errors are discussed in Sect.3.3.
Ice displacements were generated from the SLCs using the same intensity cross-correlation offset tracking method used for the initial co-registration, but with smaller window sizes.No terrain correction was applied, as resultant uncertainties were small (Sect.3.3).SLCs were oversampled by a factor of two, and all offsets with signal-to-noise ratios (Strozzi et al., 2002) below 6.0 were eliminated.We identified optimum search window sizes by testing throughout the parameter space; optimum sizes were 96 × 156 pixels (nominally, 721 m × 726 m in range and azimuth) and 100 × 200 pixels (nominally, 748 m × 629 m), for RADARSAT and PALSAR respectively.For TerraSAR-X, we use a window size of 256 × 256 pixels (232 m × 474 m).
In a few cases, slightly larger windows were used to improve results in areas of poor correlation.Results were not dependent on window size as long as the window was not large enough to extend through shear zones or onto stable ground.Shear zones on Bering Glacier are generally wider than 1.5 km (Fatland and Lingle, 2002).
Erroneous offsets were eliminated with a highly effective culling routine.The routine first uses a preliminary version of the Randolph Glacier Inventory 1.0 (nearly identical to RGI1.0) (Arendt et al., 2012) to remove all off-ice offset vectors.Next, offset vectors v are filtered by orientation.We define a median vector ṽw i,j which has the value of the median range and azimuth components of all vectors within a moving window with a center at location i, j and width w.The angle between v i,j and ṽw i,j can be defined as .
If θ exceeds a threshold value, v i,j is removed.The orientation filter is run iteratively with threshold values of 24, 18 and 12 degrees.Finally, remaining vectors are filtered by direction and magnitude where v i,j − ṽw i,j must not exceed a threshold value.This routine removes more than 99 % of erroneous vectors and removes very few false negatives.Finally, each scene is manually inspected to remove any remaining erroneous vectors and inspected for other problems.
Displacements were geocoded using the ASTER GDEM (METI and NASA, 2011) and SAR imaging geometry.The GDEM is generated through automated processing of many stereo pairs acquired between 1999 and 2011.The GDEM is the most up-to-date DEM available but is of poor quality, having many artifacts on the BGS.Since the geocoding was performed after offset derivation, the DEM artifacts do not affect the calculated velocities and only induce negligible errors in the geolocation process.All geocoded displacements were gridded onto the same 30 mm UTM grid for comparison.Strain rates were computed following Nye (1959) and Bindschadler et al. (1996).Ice displacements were extracted along longitudinal profiles shown in Fig. 1 for the BGS, WB and Jefferies/Tana glaciers (profile not shown).For purposes of clarity, we divert the BGS profile off the centerline from BGS-115 to BGS-140 to circumvent data voids (Fig. 1).Centerline velocities were higher, but extremely patchy results prohibit us from knowing by how much.
We estimate our velocity uncertainty for each image pair by using the same offset tracking method (Strozzi et al., 2002) on stable-ground where offsets are assumed to be zero and any measured offset represents an error.Suitable areas for uncertainty estimates were manually delineated in each pair.Areas with steep topography or visible geometric distortion (foreshortening, layover, or shadowing) in the SAR images were excluded.These uncertainty estimates address errors associated with image co-registration, the lack of terrain correction, and the signal-to-noise ratio accepted in the offset tracking routine.
The nadir pointing laser system results in a single track of surveyed points along the flight path, spaced roughly 1.0 to 1.5 m apart, whereas the LiDAR system results in a 500-m-wide swath with roughly one surveyed point per square meter.Vertical accuracy of the surveyed points ranges from ±0.1 m to ±0.3 m and is dependent largely upon the quality of the trajectory solution for the aircraft (position and orientation from GPS and inertial measurement unit onboard the aircraft).To calculate changes in elevation between nadir pointing laser surveys, we followed techniques and error estimations as described by Arendt et al. (2008).Calculating elevation changes between nadir laser surveys and LiDAR surveys, and between LiDAR to LiDAR surveys is done by interpolating the later survey onto a 5 m gridded surface, then sampling that surface at the earlier survey's points.In both of these cases, the inherent ambiguity between local slope (from the old point to the new point) and elevation changes is greatly mitigated relative to comparing two single tracks of points.To minimize seasonal effects on elevation changes, each interval was flown within 8 days of the date of the previous survey, with the exception of 39 days (late) in 2008.

Efficacy of PALSAR in maritime climates
Not all SAR pairs produced useable offset results; the climate in Southeast Alaska rarely leaves a glacier surface unaffected by heavy snowfall or significant melt over a ∼ 1 month interval.We find that, in this environment, temporal de-correlation is event (storm) based, rather than a gradual process; thus avoiding orbit intervals that span extreme storms or melt events can increase chances of obtaining good velocity results.In this climate, C-band and X-band sensors generally require surface definition such as crevassing and cannot track speckle alone (ERS ice-phases with shorter repeats excluded), which limits spatial coverage significantly.
In contrast, L-band was extremely robust and able to track speckle reliably in this climate despite PALSAR's longer 46day-orbit interval.This strength is due to the longer wavelength providing deeper penetration and more stable scattering from subsurface snow and firn.TerraSAR-X was also effective because of the extremely short orbit interval (11 days) and high spatial resolution.However, it may be limited in observing slower moving ice.X-band decorrelates relatively quickly, thus requiring shorter intervals, which reduces accuracy.

Offset tracking uncertainties
We calculated uncertainty offset fields for 28 PALSAR pairs and 13 RADARSAT pairs.An average of 5030 offset values was calculated per image pair.The stable-ground offsets/errors are not normally distributed, primarily because of extreme outliers.Thus, we use robust statistics including the median, interquartile range (IQR) and normalized median absolute deviation (MADn, a robust equivalent to standard deviation) (Maronna, 2006) to describe stable-ground offset fields (Table 1).Dispersion metrics quantify stochastic uncertainty, while the mean absolute median provides an estimate of bias.
Both platforms have mean absolute median values < 0.01 m d −1 in both range and azimuth directions; random errors are larger (Table 1).PALSAR has larger stochastic errors in the range direction than in azimuth, whereas RADARSAT has isotropic stochastic error.PAL-SAR also has significantly larger perpendicular baselines than RADARSAT; thus it is likely that the larger errors in the range direction are a consequence of our lack of topographic corrections.This uncertainty is still too small to affect results.
During the surge, there is extreme velocity variability that has the appearance of noise in longitudinal profiles (visible in Fig. 2).Visual inspection of the velocity maps used to generate the profiles suggests that much of the spatial variability is indeed real and not noisy results.However, rapid changes in surface conditions during the surge lead to a few errant velocity vectors that appear as noise in the profiles.Surge phase image pairs do not have unusually high calculated uncertainties.We assume that these additional errors are stochastic and should not affect our conclusions.

Velocity results
Between November 2007 and early March 2008, BGS velocities were consistently ∼ 1 m d −1 (Fig. 2) from BGS-30-BGS-110.In late March-April 2008, flow velocities accelerated 20 %, to 1.2 m d −1 (at least between BGS-100-BGS-130).This acceleration is confirmed by a concurrent increase in the number of ice quakes observed by a seismic array placed on the BGS at 110 km (LeBlanc, 2009)  Unfortunately no SAR data were obtainable between May 2010 and July 2011.During this period, aerial observations found decreasing crevassing until January-February 2011 when classic surge morphology (Herzfeld and Mayer, 1997) appeared again on the BGS (Larsen, 2012;Molnia and Angeli, 2011).These observations would suggest slower velocities throughout the summer and fall 2010 and a second high-velocity stage beginning in January-February 2011 (Larsen, 2012;Molnia and Angeli, 2011).
The only velocity data obtained during the second highvelocity stage were a single TerraSAR-X pair, obtained on 5-16 July 2011, which captured an 11-day interval of velocities from BGS-80-BGS-123 (Fig. 3).The high resolution of TerraSAR-X reveals multiple arcuate propagation fronts on the glacier surface (Fig. 3) with velocities exceeding 9 m d −1 .Calculated longitudinal strain rates across the fronts exceed 6 a −1 but are likely much higher in reality due to smoothing effects of the offset tracking routine.
These propagation fronts appear as extreme variability in the longitudinal profiles (Fig. 2).While these smaller fronts are quite evident, the velocity data do not show evidence of a uniform propagation front moving across the length of the BGS.Such a front would show a temporal delay in acceleration as one moves upstream.If we examine the period of acceleration between January 2008 and February 2009 at two locations (the lower BIV, 80-90 km, and the BGS piedmont, 110-130 km), we see that the BIV reaches 45 % and 86 % of its peak velocity at the same time that the BGS piedmont reaches 33 % and 70 % of its peak velocity, respectively (Fig. 2).Therefore, the BIV reached a higher percentage of its total acceleration earlier than the BGS piedmont.Given the variability within the velocity data, we interpret this to imply that the BIV accelerated at least in unison with the Bering Glacier, if not before.
We cannot make the same comparison with the WB due to lack of data.However, the WB accelerated in unison with the BGS in February 2008 -measurable acceleration (up to 26 %) extended 20 km upstream of the confluence (Figs. 2,  4).The WB was back to quiescent velocities in November 2008, at the same time the BGS was accelerating rapidly.It accelerated again in January 2009 when the BGS was flowing 7 m d −1 .By 2010, BGS velocities had slowed but the WB continued to accelerate.
Surge morphology was visible on Bering Glacier through summer 2011.A time-lapse camera placed in the Grindle Hills (Burgess, unpublished data) looking in a northwesterly direction (location in Fig. 1) between 21 July and 15 September captured little ice movement (exact velocities have not been calculated).By fall 2011, continued aerial observations found the glacier surface had smoothed out significantly, but no velocity data are available to confirm the surge termination.

Changes in surface elevation and driving stress
Altimetry data extend back to the end of the previous surge in 1995, with intervals of 5, 3 and 4 yr during the quiescent phase and one-year intervals during the surge.Figure 5a and b present the rate of surface elevation change and cumulative elevation change (from 1995) over each observation interval during the quiescent phase.Figure 6a and b present the same for the surge phase.
In addition to surface elevation profiles, we provide a crude approximation of driving stress at the end of each interval in Figs.5c and 6c (quiescent and surge phases, respectively).For a rough approximation, thicknesses were assumed to be the surface elevation from the terminus to the equilibrium line (bed surface is close to sea level, Conway et al., 2009).Above the equilibrium line, we assume a linear decrease in thickness to zero at the divide.Slope angles were derived from the ASTER GDEM (METI and NASA, 2011), using 6-km boxcar smooth to remove small-scale variation and data artifacts.We assume that this DEM represents the glacier surface in 1995, and sequentially adjust the geometry with the surface elevation change data provided by altimetry.Driving stress is then calculated simply as ρgh sin α where ρ is the density of ice, h is the ice thickness and α is the slope angle.Our crude assumptions on ice thickness and slope from the DEM will not affect our conclusions.This is because temporal changes in driving stress are sensitive primarily to the changes geometry, which are derived from precise and high-resolution altimetry data.

Quiescent phase
During the quiescent phase, a series of acceleration events set up the BGS geometry with a dynamic balance line and a trigger point at the approximate location of the 1993-1995 surge initiation.Between 1995 and 2000, a small acceleration event caused ∼ 10 m of drawdown near BGS-80 and 5 m of thickening near BGS-97 (Fig. 5a and b).Between 2003 and 2007, another larger acceleration event occurred, which drew down the reservoir zone by up to 20 m at BGS-60 and thickened Bering Glacier by up to 35 m near BGS-110.These two events were discrete in time, but both had a DBL at about the same location (BGS-90).Furthermore, both events steepened the glacier geometry and, consequently, increased driving stresses at BGS-120-BGS-130 (Fig. 5b, c).
Over the quiescent phase, the DBL migrated downstream as seen on Variegated and Medvezhiy glaciers (Raymond, 1987), and eventually set up at BGS-123 in 2007.However, all new/added mass from the accumulation zone was redistributed to a confined reservoir zone between BGS-85 and the DBL.Driving stresses throughout most of the BGS were unchanged from the end of the previous surge to the start of the recent surge.
However, driving stresses near the DBL (BGS-120 to BGS-130) increased up to 70 % (or ∼ 50 kPa) during the quiescent phase.Upon surge onset, this region had the strongest acceleration and fastest velocities, which caused a clear switch in longitudinal stress (Figs. 2, 5d).Thus, we can consider this region, between BGS-120 to BGS-130, to be the trigger zone for the first stage of the surge (grey highlight in Figs. 2, 5, 6).The DBL during the first phase of the surge (2007)(2008)(2009) was at BGS-122 -only 1 km from the DBL during quiescence (Fig. 5d).Again, this was also the approximate location that the 1993-1995 surge is thought to have initiated (Roush, 2003).

Surge phase
The altimetry data show a distinct difference between the first and second stages of the surge.Note the location of the DBL during the surge (Fig. 6a).The two intervals from 2007 to 2009 cover the first stage and have a DBL at BGS-122 -precisely where it was during quiescence.Then, during the intervals from 2009-2011, the DBL makes a discrete shift downstream to BGS-147.An explanation for this downstream shift is readily available, but requires a closer look at the first stage.
During the initial acceleration of the first stage (2007)(2008), a peak in thickening formed at BGS-130.Over the following year, this peak reoccurred, but a much higher peak also formed downstream at BGS-142.These two thickening peaks have important consequences with respect to the driving stress.The upstream side of the peak at BGS-130 www.the-cryosphere.net/6/1251/2012/The Cryosphere, 6, 1251-1262, 2012  interval can be eased through comparison with the 2010 velocity profiles in Fig. 2 (green lines).Note the two velocity peaks at BGS-90 and BGS-150 fit well with the altimetry.The peak at BGS-90 is very broad, extends well up into the BIV and shows little month-to-month velocity change in 2010.The peak at BGS-150 is sharp and accelerates 0.4 m d −1 in less than one month.The morphological differences between the upstream and downstream velocities lead us to speculate that the broad accelerated velocities and drawdown in the BIV are a remnant of the first stage, while the sharp acceleration in April 2010 and thickening over the 2009-2010 interval represent the onset of the second stage.
Over the 2010-2011 altimetry interval, we see the main part of the second stage.During this period, the surge reached the terminus and advanced the terminus 2-4 km (Turrin et al., 2011).Like the first stage, the DBL did not move throughout the entire stage, but thickening moved downstream.Drawdown during the first and second stages was remarkably similar in extent, magnitude and shape (compare 2008-2009 and 2010-2011 intervals in Fig. 6a) with the exception that the reservoir zone extended further downstream during the second stage.
The second stage did little to subdue the large undulations in driving stress and topography created by the first stage (between BGS-125 and BGS-145).Rather, the drawdown was more or less uniform from BGS-80-BGS-140.However, upon close examination of the 2010-2011 interval in Fig. 6a, one will see an extremely subdued undulation in the drawdown that matches the 2008-2009 profile.This key point suggests that the mechanisms that created the twopeaked thickening during the first stage still existed in the second stage but were massively diminished.During 2008During -2010During and 1993During -1995, there was a key shift in ice dynamics at BGS-45 (Fig. 2).Upstream of BGS-45, surge velocities remained close to quiescent velocities.Downstream of BGS-45, velocities increased rapidly in a step-like fashion.Strain rates (not shown) indicate the peak at BGS-80 (Fig. 2) is caused by longitudinal compression as the BIV joins with the WB.The cause of the deceleration at BGS-65 is less clear but is probably due to bed topography.At this point, there is extensional lateral strain as ice spreads northward and compresses ice on the north side of the BIV.

A confined active surge zone in the BIV
The north side of the BIV receives its ice from the Jefferies Glacier, which eventually flows into the Tana (Fig. 1).This ice accelerated very little -from 0.5 to 0.7 m d −1 -during the BGS surge.A wider shear zone on the north side of the BIV is visible in Fig. 1.Also notable, the longitudinal step-like accelerations seen on the BGS portion of the BIV did not occur on ice originating from the Jefferies Glacier (Figs. 1, 2).Rather, velocities were relatively uniform from the Jefferies to the Tana.The Bagley fault runs along the BIV and possibly forms a longitudinal ridge structure that separates BGS ice and the Tana Glacier ice.This structure likely extends towards the north edge of the West Bagley and diverts the majority of the ice southward into Bering Glacier (Plafker, 1987;Bruhn et al., 2004Bruhn et al., , 2012)).Flow velocities on the Tana Glacier changed by only ∼ 0.1 m d −1 throughout the 2008-2011 BGS surge; the highest velocities were actually prior to the surge, in 2007.

Discussion
After a 13-yr quiescent phase, the BGS began a full-scale surge in May 2008 that appears to have ended in summer 2011 (uncertainty on termination will be addressed later in this section).While we do not have any data on bed conditions during the surge, close examination of the velocity and altimetry data allow us to make inferences about the relative amount of drag provided by the bed.Most importantly, we can reach conclusions about the basal hydrology by examining the persistency of basal drag features during surge evolution.
During the quiescent phase, small-scale acceleration events occurred that relocated any new mass in the accumulation zone to a small reservoir zone just downstream of the BIV/WB confluence at BGS-85-BGS-123.The fact that large acceleration events did not occur between 2000 and 2003 suggests that these acceleration events are not purely a consequence of increases in deformational velocity; rather, basal sliding likely plays a role.Such events are similar to events observed on Medvezhiy and Black Rapids Glaciers, acknowledging key differences in geometry and size (Dolgushin and Osipova, 1978;Raymond, 1987;Heinrichs et al., 1996).One consequence of the quiescent acceleration events was, when the surge began, the entire BIV was at the exact same level as it was at the end of the previous surge -driving stresses were unchanged as well.Thickening from both quiescent phase acceleration events stopped abruptly at the DBL; consequently, driving stresses increased by 70 %.This area represents a key transition in longitudinal stress and can be considered the trigger zone for the first stage.
During the surge, the DBL remained stationary during each stage but made a discrete downstream shift between the two stages (Fig. 6a).We hypothesize that this downstream shift is a consequence of driving stress evolution during the first stage.During the first stage, thickening occurred at two peaks at BGS-130 and BGS-142, both of which are upstream of the DBL for the second stage.Downstream of the second peak, thickening and steepening around BGS-145 increased local driving stresses and created a new trigger point and DBL for the second stage.
The two thickening peaks are the result of compressive longitudinal strain as the high flow speeds decline in the downstream direction.This reduction in flow speed and thickening is likely due to an area of relatively higher basal drag that resists rapid flow from upstream.The second peak had provided drag to stop the entire surge and prevent any thickening downstream of BGS-150.
During the second stage, the two points of relatively high drag largely disappeared.Thinning and extensional flow persisted past these points and allowed the surge to extend to and advance the terminus in 2011.This progression is very similar to that of the 1982 surge of Variegated Glacier (Raymond, 1987;Kamb et al., 1985).Thickening and compressional flow during the second stage occurred downstream of a new DBL at BGS-142 -immediately downstream of where driving stresses were elevated by the first stage.Note that the second stage did little to subdue the two thickening peaks from the first stage; rather it drew down the BGS uniformly, from BGS-65 to BGS-142.Without knowledge of the bed conditions, we can only speculate as to what reduced the basal drag at these two sticky points during the second stage.
One speculative explanation lies in the driving stress and is qualitatively supported by conclusions by Kamb (1987).Throughout the quiescent and surge phases, we see a reoccurring process.Areas with low driving stresses provide drag that causes ice to pile up behind the drag feature and E. W. Burgess et al.: Surge dynamics on Bering Glacier consequently increase driving stresses at the drag feature.In winter, the channelized system is destroyed and, if local driving stresses are elevated to a threshold, the elevated driving stress helps to promote and maintain a distributed hydrologic system and prevent channelization (Kamb, 1987).Low effective pressures can thus be maintained for extended periods of time.The shear stress on the bed is reduced, and the ice is able to accelerate further downstream until it reaches another high drag location.
The existence of these high drag spots does not appear to be entirely a function of driving stress.The small undulations in thinning rate between BGS-120-BGS-150 on the 2010-2011 altimetry interval suggest that there are features on the bed that continue to provide additional drag despite the collapse of the channelized system.Thus these areas at BGS-130 and BGS-142 could have some combination of higher surface roughness features and/or more hard beds, though this is highly speculative.
There is an obvious hitch with this hypothesis.The driving stress in 2011, after what is believed to be surge termination (Fig. 6c), is still elevated at sticky spot locations.This would imply that either (1) our explanation is missing a key component, or (2) the BGS is still not in a stable geometry and the surge is not yet over.We have no velocity data after our TerraSAR-X pair was acquired in July 2011.Aerial observations have seen moderate crevassing in 2012 and previous surges terminated with a large outburst flood (Molnia, 2008), which has not yet been observed with the recent surge.
The issue of surge propagation was only loosely observed during the previous surge, and our observations for the recent surge indicate a different picture.At the onset of the first stage, velocities did not only accelerate near the trigger zone; rather the majority of the BGS between BGS-50 and BGS-140, and the WB, all accelerated at the same time.Thus, the initial acceleration cannot be due to kinematic wave propagation and must be to hydrologic pressurization of most of the BGS in May 2008.If one looks at the four altimetry intervals during the surge, drawdown follows a consistent curve and eventually stops at about BGS-75, except for the 2009-2010 interval.Thus for the majority of the first stage and the second stage, we see no large-scale upstream propagation at all.In 2009-2010, drawdown did extend further up into the BIV but this drawdown and acceleration occurred with a slight decrease in local driving stress.Throughout the entire surge, driving stresses upstream of the DBL barely changed at all.Therefore, we conclude that the extensive thinning from BGS-45 to the DBL is purely a consequence of elevated hydrologic pressures at the bed, not due to evolving driving stress from kinematic wave propagation.Areas closer to the trigger zone accelerated more, perhaps simply due to higher basal water pressures, thus creating extensional flow upstream of the trigger zone and thinning.Downstream propagation, however, appears to be closely linked to changes in the driving stress.At persistent sticky points, elevated driving stress appears to reduce basal drag and allows high flow velocities to continue further downstream.
Why there does not appear to be high drag areas above BGS-123 is an interesting question.One simple explanation could be that the bed conditions are simply different above BGS-123 and are more amenable to maintaining pressurized distributed systems and low basal shear stresses.Another explanation could be that the bed upstream of BGS-123 is not significantly different from the bed downstream, and instead, the reason why there are not sticky spots above BGS-123 is because of higher driving stresses.If we consider the distribution of driving stress throughout the BGS length, the region between BGS-80 to the BGS-123 has significantly higher driving stresses than elsewhere -near 100 kPa.This area is also the reservoir area for the surge, and the area where rapid acceleration is seen at surge onset.Upstream of BGS-45, where driving stresses are low, ice may accelerate slightly due to kinematics, but surge velocities (1993-1995 and 2008-2010) and quiescent velocities (2007) vary by little (< 0.35 m d −1 ).During the previous surge, Fatland and Lingle (2002) noticed little short-term acceleration in velocity upstream of BGS-45.Thus, these data suggest that much of the BGS system requires a critical basal shear stress (80-100 kPa) to facilitate rapid sliding.Given our crude driving stress model, this interpretation is qualitative and the magnitude of the driving stress should not be taken exactly.
We have little evidence to address the question of why the surge evolves in two stages.The 1993-1995 surge had two phases that both began in winter and terminated in summer (Molnia, 2008).In the recent surge, the first stage began in spring and the second began in winter.Termination dates are unknown with certainty but could be in summer, in particular for the second stage.One explanation (Kamb, 1987;Eisen et al., 2005) is high ablation rates during summer can force the glacier to return to a channelized drainage system, increasing the effective pressure and basal shear stress.Over the course of fall and winter, the drainage system is able to close and revert back to a high-pressure distributed system, thus allowing the second phase to initiate.
In the case of the recent BGS surge, it is interesting that the first acceleration period lasted at least ten months (May 2008-February 2009) and took place primarily in fall but also in summer.The 2008 summer however was exceptionally cold (Alaska Climate Research Center, 2012), thus likely produced relatively little meltwater input to the bed and could have allowed the distributed system to persist into fall.The first stage terminated sometime during 2009, and velocity increased again the following winter.Thus, the BGS two-stage cycle could be explained by this theory as well.Surface velocity data, altimetry and calculated driving stresses reveal glacier dynamics throughout a complete surge cycle of the Bering Glacier system.Dynamics throughout this cycle show distinct similarities to surge cycles observed on much smaller glacier systems, which suggests consistency in mechanism despite a huge difference in size.We find that areas capable of rapid acceleration are confined to areas of relatively high driving stresses.Since driving stresses decline closer to the terminus due to thinning ice, periodic acceleration in the area of high driving stress causes longitudinal compression and thickening/steepening geometry.This thickening/steepening effectively expands the area of high driving stress and allows rapid flow to advance further downstream.However, as the driving stress evolves during each stage, the basal drag does not dynamically adjust to the evolving driving stress.Rather, we find that the rapid flow must shut down and reset for another stage before the DBL can move downstream.We suggest that high driving stresses are maintaining an existing distributed system that forms over a winter season but cannot initiate a transition from a channelized system to a distributed system given additional driving stresses.These conclusions are still highly speculative as no data on basal water pressures are available.Further analysis and modeling could provide more insight into the bed mechanisms responsible for these observed dynamics.

Fig. 1 .
Fig. 1.Composite velocity map of the BGS and Tana Glacier in Winter 2010.The BGS and WB longitudinal profiles apply to Figs. 2 and 4, respectively.Profiles begin at ice divides and distances increase in the direction of flow.White glacier outline provided by Armstrong et al. (2005).

Fig. 2 .
Fig. 2. Bering Glacier system longitudinal velocity profile (location in Fig. 1).Dates of image pairs are denoted by the colored boxes in the timeline.The width of each box indicates the length of the interval; heights of the boxes are for visual clarity only.Colors gradually darken as time moves forward each winter.Thin colored boxes below the timeline show dates of altimetry intervals in Figs. 5 and 6.BGS velocities in January 1994, during previous surge, are shown as grey data line.The location of the trigger zone for first stage is indicated by the grey box.
tracking results were unavailable for late spring and summer 2008, but the same seismic array(LeBlanc, 2009) found the number of ice quakes increased by an order of magnitude in the first two weeks of May, indicating an onset of activity much greater than the previous year's spring speedup.Between September 2008 and February 2009, the BGS accelerated progressively from BGS-80-BGS-135.Maximum observed velocities for the first stage were ∼ 7 m d −1 .Actual peak velocities were likely higher, as this maximum velocity was observed in a side shear zone where data were unavailable at the glacier centerline (BGS-130).In January-April 2010 surface velocities from BGS-110-BGS-140 slowed to quiescent speeds and velocities in the lower BIV (near BGS-90) were ∼ 2 m d −1 .At the terminus, a narrow velocity peak at BGS-150 accelerated ∼ 0.4 m d −1 between March and April 2010 (Fig.2).

Fig. 3 .
Fig. 3. High-resolution velocity field during second phase of the surge over 11 day interval (5-16 July 2011) from TerraSAR-X.Location of propagation fronts marked in inset map.

Fig. 4 .
Fig. 4. West Bagley longitudinal velocity profile (location in Fig. 1).Dates of each image pair are denoted in the timeline below as in Fig. 2. Colors gradually darken as time moves forward each winter.The grey data line shows WB velocities in January 1994.Note scales are different than Fig. 2.

Fig. 5 .
Fig. 5. Results from airborne altimetry on the BGS profile.Flight intervals are indicated on the timeline at the bottom.Blue and red hues correspond to quiescent and surge phases, respectively.Colors lighten though time.Grey box and vertical line mark the trigger zone and DBL for the first stage, respectively.(a) Surface elevation change rate during quiescence.(b) Cumulative elevation change from 1995.Colors correspond to the surface profile at the end of each interval.(c) Calculated driving stresses.Black line represents driving stress in 1995; colors correspond to the driving stress at the end of each interval.(d) Surface elevation change rate during the surge phase (same as Fig. 6a).

Fig. 6 .
Fig. 6. Results from airborne altimetry on the BGS profile cont'd.Colors, intervals and notation represented as in Fig. 5. Dotted line represents DBL for the second stage.(a) Surface elevation change rate during the surge phase (same as Fig. 5d, provided for comparison).(b) Cumulative elevation change from 1995.Colors correspond to the surface profile at the end of each interval.(c) Calculated driving stresses.Black line represents driving stress in 1995; colors correspond to the driving stress at the end of each interval.

Table 1 .
Statistics describing the 2-D distribution of measured ground displacements in ice-free terrain (units in m d −1 ).The STDEV, MADn IQR, and absolute median are calculated for all stable-ground offsets within each image pair.Subsequently, the mean STDEV, MADn, IQR and mean are calculated for all image pairs.