1. Introduction

2. Approach

Much of the research on surges is based on a few case studies of mostly smaller glaciers, such as Variegated Glacier, Alaska, observed during its 1982/83 surge (e.g. Reference PostPost, 1960; Reference KambKamb and others, 1985; Reference Humphrey, Raymond and HarrisonHumphrey and others, 1986; Reference RaymondRaymond, 1987) and Trapridge Glacier, Yukon (Reference Clarke, Collins and ThompsonClarke and others, 1984, and many later studies), which are easier to equip with geophysical instrumentation but behave differently than large glaciers. Reference Meier and PostMeier and Post (1969) distinguished three types of surging glaciers: (1) large-sized with high surge velocities, large displacement and large lowering of ice reservoir, (2) large-sized with low surge velocities and small vertical and horizontal ice displacement, and (3) steep, small glaciers with small reservoir lowering and ice displacement. Bering Glacier is a prototype of a large surge-type glacier of type 1, whose size is comparable to large fast-moving ice streams in Greenland, which individually play a role in understanding ice-sheet stability and discharge into the ocean. Observing and analyzing surge stages in a large glacier requires a different method than measurement of geophysical properties on the ice, which is impractical. The approach used here is a form of automated spatial analysis of airborne remote-sensing data, which includes automated image analysis, GPS data analysis and analysis of elevation data. The analysis progresses by application of a classification based on a set of spatial parameters that together capture components of the dynamics of surge progression.

Analysis of altimeter data allows elevation changes to be derived which document the mass transfer from the reservoir area to the receiving area. In addition to the large-scale elevation changes, high-resolution altimeter data also capture changes in small-scale topography associated with characteristic crevassing. The latter are needed as they contain the vertical component in deformation matrices or tensors describing deformation, crevassing and generally small-scale kinematics associated with the surge (Reference Herzfeld, Clarke, Mayer and GreveHerzfeld and others, 2004). The analysis of elevation changes and classification of crevassed ice surfaces forms a central part of this paper.

3. Observations

During an aerial observation campaign in September/October 2011, the following data were collected: (1) digital photographic data, (2) digital video data, (3) GPS data and (4) laser altimeter data (see Fig. 1). Flyable weather permitted geophysical observations on 25 September (flight 1) and 26 September (flight 2) 2011. Laser altimeter data were collected with a lightweight laser altimeter operating at 905 nm (near-infrared) wavelength with 400 Hz data collection rate. The instrument is a single-beam nadir-pointing Laser Technology, Inc., Universal Laser System (ULS) rangefinder that is used to measure the distance from the aircraft to the ground surface. The instrument also measures laser pulse return power. Beam divergence is 3 mrad (30 cm diameter footprint at 100 m range). The ULS collects surface topography measurements at 400 Hz (derived from 4000 Hz system-internal measurements ), corresponding to a 7.5 cm along-track sample spacing when flown at 30 m s–1, and has a 60 cm footprint when flown at 200 m above ground level. Accuracy of the instrument is 2 cm according to manufacturer’s specifications, and measurement error (observed standard deviation of measurements over a known, flat surface) is 0.6 cm (Reference Crocker, Maslanik, Palo, Adler, Herzfeld and EmeryCrocker and others, 2011; Table 1).

Fig. 1. Location and flight paths of data collection campaign over Bering Glacier, Steller Glacier and Bagley Ice Valley, September 2011. (a) Location and data collection: Bering Glacier is the glacier that calves into a series of proglacial lakes near 60.158 N, –143.58 E; Steller Glacier is west of Bering Glacier; Bagley Ice Valley is seen in the northern part of the map. Data collected over tracks marked ‘GPS’ are digital photographic data, digital video data and GPS data. Over tracks marked as laser runs, laser data were collected in addition to photographic, video and GPS data. Tracks indexed 1 were observed on 25 September 2011, tracks indexed 2 on 26 September 2011. (b) Location of Bering Glacier, Alaska (terminal lobe marked by red dot).

Table 1. Universal Laser System (ULS) specifications

To meet the rapid time constraints of observing the surge, altimeter data were collected by simply holding the altimeter out of the plane secured to an outstretched arm of the scientist holding it, pointing at approximately nadir. During operation, the ULS and a GPS unit were connected to a laptop computer (a MacBook), and laser measurements and GPS times were co-registered. In post-processing, altimeter measurements were corrected by GPS positions (aircraft positions). The same laser system was previously used as part of an integrated system on manned and unmanned aircraft over glaciers and sea ice (Reference Crocker, Maslanik, Palo, Adler, Herzfeld and EmeryCrocker and others, 2011). Because mounting an instrument outside a plane requires the approval of the US Federal Aviation Administration (FAA), in addition to longer-term coordination with the flight service, hand-held operation was the only option to collect such data during the rapidly evolving surge. However, this simple method worked surprisingly well, as the following analysis will demonstrate. Individual crevasses are visible and different crevasse provinces can be identified in the data. Maximal crevasse depth resolved is almost 60 m (Figs 1 and 2).

Fig. 2. Crevasse profiles and and ice-surface profiles observed during Bering Glacier surge, 25 September 2011. Laser altimeter data, corrected with GPS data. (a, b) Overview (a) and close-up (b) of location of examples shown in (c–e). (c–e) Elevation vs along-track distance; aspect ratio 5 : 1 . Blue: data points; red: piecewise linear interpolation. (c) Area 4a (note two different surface provinces separated by a 60 m deep rift; see Fig. 3) (flight 1 , laser run 1a, 25 September 2011); (d) area 3a (flight 1 , laser run 1b, 25 September 2011); and (e) area 1a2 (flight 1 , laser run 1b, 25 September 2011).

In flight 1 , laser data were collected in two along-flow profiles, from a cross-flow line connecting Khittrov Hills and Grindle Hills to a cross-flow line at Point 3687 in central Bering Glacier (laser run 1a), with a second line to validate the method (laser run 1b). In flight 2, longer tracks were flown in lower and central Bering Glacier, including Tashalich arm, the main lobus and two profiles approximately normal to the flowline (Fig. 1). The cross-flow profiles were repeated offset such that videos and photographs of the same locations could be obtained. A total of 141.48 km of laser altimeter data were collected (23.2 km in flight 1 , run 1a; 8.7 km in flight 1 , run 1b; 55.8 km in flight 2, run 2a; 53.78 km in flight 2, run 2b). Laser data used in the following analysis were corrected using GPS data collected at the same time.

Figure 2 gives examples of laser-altimeter profile sections, which indicate (1) that different structural provinces can be distinguished, (2) that the laser beam reaches into the crevasses, so that crevasse depth can be measured, and (3)0that laser points are sufficiently dense on the ground to capture the ice surface morphology, including roundedness of crevasse edges and small surface features between major crevasses. These principal properties of the data facilitate the analyses in this paper. Prior to the analyses, we introduce the originally geological concept of a structural province (Section 4) and the methodological concept of mathematical characterization and classification for discrimination of structural provinces (Section 5). The following elevation data analysis has three parts: First, we analyze the first-order variables that have long been of interest to glaciologists and mountaineers, namely average and maximal crevasse depth and average crevasse spacing (Section 6). Second, we derive characteristic parameters for crevasse provinces (Section 7). Third, absolute elevation changes are calculated from 2010 and 2011 data (Section 8).

4. Structural Glaciological Provinces

The principal concept is that the dynamics of fast-moving ice manifest themselves in crevassing, and consequently the deformation history of the ice can be reconstructed by analysis of crevasse patterns (Reference Herzfeld and MayerHerzfeld and Mayer, 1997; Reference HerzfeldHerzfeld, 1998). Crevassing is a form of brittle deformation that occurs when local forces exceed a threshold. Examples of deformation include extension, compression and shear. A structural province is an area that is homogeneous with respect to structural glaciologic properties (properties caused by deformation) and that is maximal (i.e. any larger area would contain a second structural glaciologic type and hence be inhomogeneous). Examples of structural provinces are seen in Figure 3: a large rift separates two structural provinces, a crevasse field of rhombic crevasses and a region that is almost uncrevassed, where older, closed crevasses dominate. The rift is most likely the result of shear. The field of rhombic crevasses to the north of the rift is resultant of shear and extension caused by surge kinematics, whereas the surge has only marginally affected the area south of the rift.

Fig. 3. Rift and adjacent crevasse provinces, Bering Glacier, 25 September 2011. Photograph by U.C. Herzfeld. The laser profile given in Figure 2c and analyzed in Figures 6 (classification parameters) and 7 (typical variograms of structural provinces) crosses the rift at a location shown just below the middle of the photo at a small angle, where the rift appeared deepest; flight direction is from the heavily crevassed region in the north to the almost uncrevassed region south of the rift (viewers near-field), where a few small crevasses exist near the rift (only). Approximate location of flight path shown in yellow.

Structural geologic principles provide links between dynamics, kinematics and deformation, which can be physically formalized and quantified using continuum mechanics (Reference MeansMeans, 1976; Reference SuppeSuppe, 1985; Reference Twiss and MooresTwiss and Moores, 1992; Reference Ramsay and LisleRamsay and Lisle, 2000; Reference LiuLiu, 2002; Reference GreveGreve, 2003). Use of crevasse patterns as a source of geophysical information in an objective and automated approach was first applied in analysis of image data (video data) from the 1993–95 surge (Reference HerzfeldHerzfeld, 1998; Reference Herzfeld and ZahnerHerzfeld and Zahner, 2001). This approach leads to classification of surge-crevasse types correctly in 95% of all cases and allows association of deformation characteristics to specific images (Reference Herzfeld and ZahnerHerzfeld and Zahner, 2001). Applications of the method facilitate mapping of deformation provinces in surging and continuously fast-moving glaciers (Reference Herzfeld and MayerHerzfeld and Mayer, 1997; Reference HerzfeldHerzfeld, 1998; Reference Herzfeld, Mayer, Feller and MimlerHerzfeld and others, 2000; Reference Mayer and HerzfeldMayer and Herzfeld, 2000, Reference Mayer and Herzfeld2001, Reference Mayer and Herzfeld2008; Reference Herzfeld and ZahnerHerzfeld and Zahner, 2001; cf. Reference Vornberger and WhillansVornberger and Whillans, 1990; Reference Marmo and WilsonMarmo and Wilson, 1998; Reference RistRist and others, 1999).

However, the analysis is based on two-dimensional data (imagery), which provide mathematical information for derivation of deformation matrices in two dimensions, as described in theory in Reference Herzfeld, Clarke, Mayer and GreveHerzfeld and others (2004). In this paper, we utilize high-resolution laser profilometer data to explore the third dimension. For example, laser altimeter data of the rift and its neighboring provinces are seen in Figure 2c.

The analysis in this paper focuses on crevasses induced by surge kinematics. As is well known, in any glacier, crevasses form induced by topography of the glacier bed and margins. The interaction between surge forces and preexisting topographically induced crevasses is discussed by Reference Herzfeld and MayerHerzfeld and Mayer (1997) and Reference HerzfeldHerzfeld (1998).

5. Methodological Concept: Mathematical Characterization and Classification for Discrimination of Structural Provinces

The analysis of laser altimeter data is based on an approach that quantitatively relates parameters derived from elevation profiles to structural glaciological provinces. In this section, we introduce the concepts connecting observations and attribution of glacial change. These mathematical principles are summarized in the geostatistical classification method (Reference HerzfeldHerzfeld, 2008). For the surge application we include the additional parameter of roundedness of crevasse edges, as surge crevasses are typically clear-cut because of the sudden occurrence in a single deformation event (Reference Lingle, Post, Herzfeld, Molnia, Krimmel and RoushLingle and others, 1993; Reference HerzfeldHerzfeld, 1998).

The objective of this paper is to present laser altimeter data collected during the current (2011) surge of Bering Glacier and to demonstrate how these laser altimeter data can be used in a systematic way to understand certain quantitative aspects of the surge process. The approach combines two methods, (1) elevation change determination and (2) spatial classification–characterization of crevasse properties, and hence covers large (in 1) and small (in 2) spatial scales, as well as a classic quantitative method (in 1) and (in 2) a new method that is being developed specifically for analysis of crevasse formation during a surge.

For a glacier, the classification has three components: (1) characterization, (2) classification and (3) segmentation. The objective of mathematical characterization is to determine a set of parameters that uniquely describe an object. The objective of classification is to associate a given object to one of a number of classes; so that the class association can be carried out for each case automatically, a rest class is required to collect all objects that do not belong to any of the characteristic classes. In each new application, the characterization problem must be mastered before the classification problem. By moving a classification operator over a large spatial dataset (e.g. an image), each window is associated to a class, and a segmentation may be achieved. A segmentation into reasonably coherent subsets or segments may require smoothing; however, in a good classification this is not necessary (for examples and details see Reference Herzfeld and HigginsonHerzfeld and Higginson, 1996). The characterization–classification approach stems from a family of methods derived to address the general problem of retrieving hidden, overprinted or secondary information from spatial data (Reference HerzfeldHerz-feld, 2008) and is adapted here to analyze laser altimeter data with the objective of understanding crevasse formation and surge progression. The mathematical principles are given in the Appendix.

Related to the surge phenomenon, the characterization is applied to uniquely describe (mathematically characterize) crevasse types that occur during typical surge stages. Applied to laser data, we expect to answer questions such as: Does crevasse depth increase during the surge? How does transport of crevasses affect spacing and depth of surge crevasses? Can crevasse depth and crevasse significance be used to discriminate two different generations of crevasse opening (as occurs when a new kinematic wave affects a previously crevassed area)? Can we distinguish surge crevassing that affects previously smooth ice from surge crevassing that affects previously crevassed ice? In this paper, we demonstrate the principles of this type of analysis and show that characteristic parameters can be derived. To this end, we establish links between crevasse types and associated parameters based on laser altimeter data.

The approach links the crevasse classes to their position in the surge. To demonstrate this, we select two topics: (a) the region of a large rift that separates two flow regimes and (b) the position of crevasse types relative to the bulge and the front of the surge (from synopsis of crevasse types and elevation change analysis).

Once the characterization is established, an automated classification can be computationally implemented and carried out. For example, an automated classification algorithm can be presented with a new section of laser data of a surge glacier and associate a crevasse type. The third step is the application of the automated classification to a dataset that covers the entire glacier (or a sufficiently dense grid of tracks), to map deformation provinces. (Such a dataset does not exist currently.) A map of deformation provinces gives the state of the surge at one point in time. Notably, only one set of data collected at one point in time is required for analysis (not two, as is the case with interferometry or elevation-change detection or derivation of velocity from image data (Reference Herzfeld and MayerHerzfeld and Mayer, 1997)). In this paper, the characterization is carried out. We present the data and a proof of concept for the overall mapping segmentation, along with applied results that follow directly from the characterization step. Then we combine that with results from the elevation-change analysis.

6. Crevasse Depth and Crevasse Spacing

Previously, glaciologists believed maximal crevasse depth in any glacier is 30 m, but much larger depths can be reached. The rift shown in Figure 3 has a depth of ^ 6 0 m (Fig. 2). This is the maximal depth observed in the collected lidar profiles, which cover only a small subset of the ice surface. However, the rift is the most prominent crevasse feature that existed during the time of survey and visually appeared to have the largest depth of any open crevasse.

Results of along-track calculation of (1) maximal crevasse depth, (2) average crevasse depth, (3) average crevasse spacing and (4) mean curvature of crevasse edges are shown in Figure 4; all calculations are taken over a 500m along-track window. Data analysis is based on laser altimeter data that have been GPS-corrected for geolocation (elevation, position corrected by plane altimetry and position).

Fig. 4. Depth and spacing of crevasses and curvature of crevasse edges, derived from laser profilometer data collected during Bering Glacier surge, 25 and 26 September 2011. Parameters calculated from GPS-corrected laser altimeter data in 100m along-track windows. White interpolation of flight track indicates that the crevasse parameters could not be determined, because no crevasses exist within a 100 m along-track window. Superimposed on Google image for approximate reference. Image collected before the surge. (a) Maximum crevasse depth (m); (b) average crevasse depth (m); (c) average crevasse spacing (m); and (d) average curvature (m–1) of crevasse edges, used as a measure of age of crevasses since formation (for definition see Section 7.1).

Average crevasse depth ranges from 2 to 54 m depending on the crevasse provinces (Fig. 4b). Average crevasse spacing is typically 20-30m and reaches 100m in some crevasse fields. Spacing of crevasses is evaluated along-track and hence provides a relative measure for near-straight segments; actual spacing is less than or equal to the value shown here. This value is best used in combination with image analysis.

Maximal crevasse depth varies across the surveyed area, but appears to be associated with crevasse provinces, i.e. depth values are similar for a given length of track, then change to another value. Average track depth suggests the same hypothesis, but overall shows more variability than maximal depth. The hypothesis is investigated in the following.

7. Characteristic Parameters of Crevasse Provinces

Two types of parameters are used: (1) curvature or roundedness of crevasse edges and (2) geostatistical classification parameters. Geostatistical classification parameters are defined in the Appendix.

7.1. Roundedness as an indicator of crevasse age and surge progression

The motivation of this part of the analysis is the observation that surge crevasses open suddenly, as the surge wave hits a given location, resulting in characteristic sharp edges. This is well documented in many observations of the 1993–95 surge (Reference Lingle, Post, Herzfeld, Molnia, Krimmel and RoushLingle and others, 1993; Reference Herzfeld and MayerHerzfeld and Mayer, 1997; Reference HerzfeldHerzfeld, 1998). In contrast, crevasse fields surveyed with laser altimetry in the lower and central Bering Glacier formed several weeks/months before our survey. Melting and additional snowfall cause rounding of the crevasse edges. The high vertical accuracy and spatial density of the crevasse observations permit quantification of roundedness, as performed in the following analysis.

The roundedness algorithm includes the following steps: (1) discrimination of crevasse interior (A) and ice surface outside of crevasses (B) as two domains, (2) identification of the bottom of crevasses, (3) identification of the start of the crevasse edge, (4) interpolation of surface outside of crevasses, (5) calculation of differences between (A) and (B) as crevasse depth, (6) identification of the convex region in the profile, starting at the crevasse edge identified in (3) (using the edge that is located after the crevasse in observation time), (7) fitting of a circle to the convex region and (8) calculation of the radius r of the circle. (9) Curvature is defined as 1/r, i.e. larger curvature values indicate fresher crevasses. The method is illustrated for area 3A in Figure 5.

Fig. 5. Calculation of roundedness of crevasse edges from laser altimeter profiles of Bering Glacier. Area 3A. Subsets of laser run 1a, flight 1 (see Fig. 2d). Location shown in Figure 2a.

This algorithm has been applied to all laser data collected during our campaign. The result is shown in Figure 4d. Notably, curvature values in crossover locations of tracks match. A concentration of crevasses with more rounded edges is observed at 60.3–60.358 N, 143.4–143.358 W (reservoir area), over the Khittrov crevasse field at approximately (60.278 N, 143.458W) and in Tashalich arm (near (60.28 N, 143.558 W)). The surge started in the reservoir area and reached the front in Tashalich arm soon after that. The elevation change analysis for these areas is given in Section 8. The Khittrov crevasse field is characterized by wide open crevasses and fallen seracs. The curvature analysis holds potential for systematic full surveys and repeat surveys.

7.2. Geostatistical classification parameters for discrimination of crevasse types

Geostatistical classification uses parameters proven to characterize spatial surface roughness and spatial surface properties in a summarizing manner (Reference HerzfeldHerzfeld, 2008). Parameters are derived from vario functions and residual vario functions, as described in the Appendix. Here we employ the parameters pond (overall surface roughness), p1 (significance, absolute size), p2 (significance, relative value) and mindist. In the crevasse analysis, the mindist parameter yields the average spacing of crevasses in a homogeneous area (a subarea of a structural province). To demonstrate the approach, laser data from the area of the rift and the two adjacent provinces are analyzed.

Comparison of the panels in Figure 6 indicates that differences between the vario function and the residual vario function do not affect the parameters much. This means that the elevation trend over the window sizes investigated (200m) is small compared to the crevasse structure. To simplify the text, we base the following interpretation on the parameters from the residual vario function, but use parameters from both the vario and the residual vario function in the characterization in Figure 6. The two parameters, pondres and p1res allow the classification of three structural provinces: (1) the heavily crevassed region to the left (north) of the rift, (2) the rift itself and (3) the area of low-relief surface microtopography to the right (south of the rift). The latter may be an area of closed old crevasses (non-surge) and surface melt features that is in the early stages of being affected by rhombic surge crevasses (see Fig. 3). A deterministic characterization of the three main provinces is possible using only two or three components. Based on first analysis of the parameter plots in Figure 6, the rift separates two different flow regimes. A simple statistical cluster analysis, based on geostatistical parameters pond, p1 and mindist and using the centroid method for clustering and automatically searching for two clusters, confirms that the provinces separated by the rift can be automatically discriminated using this approach.

Fig. 6. Geostatistical classification parameters calculated from laser altimeter data for the region of the rift. Top panel shows laser altimeter profile. Next two panels show classification parameters derived from vario functions: red – pond; blue – p1; green – p2; black – mindist. Bottom two panels show classification parameters derived from residual vario functions: red – pond; blue – p1; green – p2; black – mindist. Parameters calculated for windows of 200 m, maximum vario function lag 150 m, windows offset by 10 m, resulting in a feature vector every 10 m. Parameters are plotted as relative values (ratio of actual parameter value and maximum parameter value in the set). All values are location-referenced to along-track distance, geographic latitude and longitude and also to Universal Transverse Mercator (UTM) north and east. Vertical lines indicate structural provinces and sub-provinces, based on visually aided analysis of the feature vectors composed of classification parameters.

However, a visual interpretation of the parameter plots in Figure 6 suggests that province 1, the crevassed region north of the rift, can be subdivided into five sub-provinces (1a, 1b, 1c, 1d, 1e ordered in along-track direction), which is not obvious from visual interpretation of the photograph in Figure 3. To refine the classification, it is useful to introduce the concept of a feature vector, with vector elements v = (pond(var),p1(var), p2(var), mindist(var),pond(res), p1 (res), p2(res), mindist(res)) for each of the four parameters, derived for two types of vario functions. We consider the variograms (vario-function graphs) from which the parameter plots are calculated (Fig. 7) as a means for a more detailed analysis, as well as a demonstration of the functionality of the characterization-classification method.

Fig. 7. Variograms for sections of rift-traverse profile. A typical variogram is given for each province and sub-province of the rift-traverse profile, identified in Figure 6. Difference between vario-gram (blue) and residual variogram (green) quantifies amount of drift. Variograms are calculated for 200 m windows, offset by 10 m in along-track direction for the entire profile section shown in top panel of Figure 6, using a unit lag of 2 m and a maximum lag of 150m (75 points, 75% of window size). Every fifth variogram is numbered. Variogram labels (variogram number - along-track distance) allow the variograms shown here to be matched with along-track distance of the profile shown in top panel of Figure 6 and geostatistical classification parameters shown in Figure 6.

As a first step, the difference between the variograms and the residual variograms in Figure 7 is a means to check for second-order stationarity of surface. The variogram is always larger than or equal to the residual variogram (see Eqn (A3) in the Appendix). Figure 7g (variogram 46) has the typical shape of a transitional variogram, with low values indicative of good spatial correlation for short distances, increasing to higher values until a sill is reached (sill value 0.6 m2 for res and 0.75 for var), indicative of degrading correlation with increasing separation of points. The variogram can increase steeply after a given distance, while in this case the residual variogram gives the characteristic pond value. As seen in Figure 6, the pond parameter is low for the non-crevassed area (but there are some surface features, i.e. the surface is not entirely smooth). The low value of parameter p1 tells us that there are no significant crevasses. In region 2 of the rift, the highest pond values of the entire laser data section are encountered; a typical variogram is given in Figure 7f (variogram 34); p1 is also high, because a large crevasse (the rift) dominates (high pond, high p1).

Provinces 1, 2 and 3 can be discriminated based on parameters pond and p1: for province 3, use pond less than 10 (in the absolute values), p1 less than 0.01 (in the relative values in Fig. 6); for the rift (region 2), pondis larger than 200; for the crevassed region (1), pond is between 0.01 and 200.

To refine the classification of province 1 into sub-provinces, a variogram is analyzed for each sub-province. For example, sub-province 1e has low values of p1, pond and p2 (not many deep and significant crevasses). Section 1d (Fig. 7d, variogram 23) has a pond value of 100-180 (second-highest to the rift) but also a high p2 value, which indicates the presence of a mix of large crevasses with (a different generation of) less deep, but relatively significant crevasses in between (p2 is the relative-significance parameter; see Appendix). We conclude that sharp surge-type crevasses formed here. This is validated by a look at the photograph (Fig. 3) and the crevasse-edge curvature analysis (Fig. 4d). Variogram 18 (Fig. 7c) shows sets of three small crevasses in each group of large crevasses, but the small crevasses are not relatively significant (and hence probably older); this is characterized by a low p2 value and low pond (50) and p1. Figures 7a (variogram 2) and 8b (variogram 9) show that sub-provinces 1a and 1b can be discriminated based on parameters mindist (average crevasse spacing) and p2, while p1 and pond are similar for 1a and 1b.

In summary, the analysis demonstrates that the classes can be discriminated deterministically. To implement an automated classification, a powerful class association method needs to be selected; an example is the connec-tionist–geostatistical method (Reference Herzfeld and ZahnerHerzfeld and Zahner, 2001).

8. Elevation Change

In this section we compare altimeter data collected by us during the surge in September 2011 with altimeter data collected by C. Larsen (University of Alaska Fairbanks) as part of NASA’s Operation IceBridge. Our flight tracks were planned without knowledge of the existence of Larsen’s observations, but the flight tracks of both campaigns are sufficiently closely located for a comparison. Therefore, regions of crossovers shown in the track plots are used to derive elevation differences (Figs 8 and 9). Pairs of crossover points were identified as follows. First, crossover regions were determined from the track plots. Each region was divided into squares, such that 50 squares cover the x-direction; then laser measurement points were aggregated into the squares, averages within each square formed for each observation year, and finally differences formed according to

(1)

Fig. 8. Laser altimeter tracks over Bering Glacier. Location of laser altimeter data collected by C. Larsen (University of Alaska Fairbanks), 2010, under NASA Operation IceBridge, and by U.C Herzfeld, September 2011, as part of US National Science Foundation project.

Fig. 9. Elevation change from 2010 and 2011 laser altimeter dataover Bering Glacier. (a) Lower central Bering Glacier (uppermost area). (b) Lower central Bering Glacier (downstream of (a)). (c) Front of Bering Glacier in Tashalish arm. Green: data collected by C. Larsen, 2010. Other colors: data collected by U.C. Herzfeld, September 2011. Squares show elevation differences in crossover locations and near-crossover locations. For locations of areas see Figure 9.

Elevation differences shown in Figure 9 indicate a surface lowering of 40-70 m on the northwestern side of the lower central part of Bering Glacier and a thickening of ?s20-40 m in lower Bering Glacier near the ice front in Tashalish arm. While the analysis is limited by the coverage of the laser tracks, the area of largest surface lowering corresponds in location to the area where the surge may have started, based on visual interpretation of crevasse patterns. Hence this area is interpreted as the former reservoir area of the surge. (The first laser run was selected to include the region where the surge may have started, based on visual crevasse interpretation and knowledge of the 1993-95 surge process.)

The area of largest thickening coincides with the area where the surge first reached the ice front. The fact that surface increase is still observed indicates that the additional mass transferred from the reservoir area has not been entirely redistributed in the receiving area, by September 2011.

9. Summary of Results

The Bering Glacier–Bagley Icefield System is currently surging (2011). The dynamics of a surge-type glacier change over a quasi-cyclic pattern that includes a long quiescent phase of normal flow and a short surge phase of fast flow. A surge glacier builds up ice in a reservoir area during the quiescent phase of the surge; this ice is then rapidly released during the surge phase. In consequence, drastic elevation changes occur during the surge phase, including rapid thinning in the former reservoir area, drawdown along the margins of the glacier and thickening in the receiving area downstream of the reservoir area. The kinematic wave associated with the surge propagates upstream and downstream, causing heavy crevassing that is characteristic of the rapid, brittle deformation. In this paper, we have investigated the large-scale elevation changes and the small-scale surface elevation and elevation-change characteristics as a means to understand the surge. With this objective in mind, airborne laser altimeter data were collected along several profiles totaling 141.48km over Bering Glacier in September 2011.

Elevation data analysis was carried out in three parts: First, we analyzed the first-order variables average and maximal crevasse depth and average crevasse spacing. Maximal crevasse depth is 60 m, reached in a rift that separates two deformation domains, indicative of two different flow regimes. Surge crevasse depth generally reaches 20–30 m.

Second, we derived characteristic parameters of structural provinces by applying the geostatistical classification method. Parameters include significance and spacing of crevasses, surface roughness and crevasse-edge curvature (indicative of crevasse age). Average crevasse spacing is typically 20–30m and reaches 100m in some crevasse fields. Because freshly formed surge crevasses have a characteristic clear-cut edge, which rounds over time due to erosion and melting, the curvature parameter provides a means for relative dating of crevasse age and hence allows us to infer the progression of the surge wave through the glacier, based on a single-time observation.

Combined into a feature vector, the parameters characterize spatial surface roughness. A classification based on the feature vector serves to objectively discriminate structural provinces, indicative of surge progression down-glacier and up-glacier. Here this has been demonstrated for the region of a large rift in central Bering Glacier, where several crevasse provinces and sub-provinces can be distinguished.

The characterization of altimeter data from the region of the rift shows that the rift is a unique feature that separates different structural provinces. In conclusion, different flow regimes exist on either side of the rift: a field of rhombic crevasses resultant from shear and extension caused by surge dynamics (north of the rift) and an area south of the rift that has only at its margins been affected by the surge (at time of survey). The rift has the effect of impeding progression of surge kinematics. The rift itself may be attributed to shear between the two dynamic provinces. This result motivates the hypothesis that the rift may be associated with a large bed-topographic feature, such as a ridge between two large channels.

Third, absolute elevation changes were calculated by differencing our 2011 altimeter data and laser altimeter data collected in 2010 under NASA’s Operation IceBridge by University of Alaska scientists. Results show surface lowering of 40–70m in the lower central Bering Glacier on the northwest side (location of the bulge that marks the former reservoir area) and thickening of 20–40 m near the front in Tashalich arm (where the surge first reached the ice front).

Elevation differencing documents the mass transfer from the reservoir area to the receiving area. Results of the crevasse profilometry and crevasse pattern analysis allow us to mathematically–physically reconstruct the rapid progression of the surge to the time of observation. The geostatistical and structural–glaciologic analysis permits us to distinguish a small number of deformation events (one to three), but is limited once many events turn the crevasse patterns into a chaos of seracs and blocks, as typically happens near the ice front. Using parameters such as roundedness (relative age) and overprinting of crevasses, more than one deformation event or crevasse generation can be identified. Especially because the surge only started a few months earlier, such reconstruction is feasible for most of the glacier. At the time of observation in September 2011 the surge was still in progress, and future observations will be needed to capture the entire surge process.

Methodologically, the paper introduces an approach that utilizes the principles of structural glaciology and mathematical characterization and classification as a means to identify crevasse provinces (structural provinces in case of the brittle deformation that dominates during a surge) and then to derive characteristic parameters related to surge progression. This is the first work in which we apply the spatial characterization–classification approach to laser altimeter data of crevasses and thus derive quantitative parameters of the vertical component of surge deformation stages. When combined with image classification using the same method (Reference Herzfeld and ZahnerHerzfeld and Zahner, 2001) three-dimensional descriptors of surge deformation can be obtained.

There are three components to this type of analysis of ice dynamics: (1) characterization (establishment of a unique set of parameters for each class), (2) classification (association of new objects to a class) and (3) segmentation (mapping of structural provinces for an entire glacier or glacier system). In this paper, the characterization step has been carried out for typical surge-crevasse types. Notably, the characterization step alone is sufficient to provide a basis for forward modeling (assuming all crevasse types that occur during a surge have been sampled with laser-altimeter and image data). To provide a parameterization for model input is one of the goals of our study. Following collection and analysis of data covering the entire glacier system throughout the surge, the structural classification will be carried out systematically, aided by connectionist–geostatistical classification.

As a further application, the same approach generalizes to investigating kinematics and brittle deformation associated with ice acceleration in temperate glaciers in or near the Arctic, including accelerations of outlet glaciers of the Greenland ice sheet in a warming climate (surge-related and not surge-related).