1. Introduction

Jakobshavn Isbræ (Sermeq Kujalleq; Ilulissat Ice Stream) calves into Jakobshavn Isfjord (Ilulissat Isfjord; Kangia) at ~69° N, 50° W in West Greenland. Jakobshavn Isbræ drains 6% of the Greenland ice sheet and has been known as the continuously fastest-moving glacier on Earth, moving at speeds of 20.6 m d^-1 (7.5 km a^-1) (Reference Pelto, Hughes and BrecherPelto and others, 1989) and 8.360 km a^-1 according to Reference Echelmeyer and HarrisonEchelmeyer and Harrison (1990) at the calving front. The glacier entered a phase of rapid retreat in 1999, which came to a halt in 2007 and continued in 2009 and 2010. The retreat has been accompanied by surface lowering and acceleration to about twice the previous velocity, modulated by a seasonal variability of slower winter speeds and thickening and faster summer speeds and increased thinning (Reference WeidickWeidick, 1984; Reference Pelto, Hughes and BrecherPelto and others, 1989; Reference Echelmeyer and HarrisonEchelmeyer and Harrison, 1990; Reference Echelmeyer, Clarke and HarrisonEchelmeyer and others, 1991, Reference Echelmeyer, Harrison, Clarke and Benson1992; Reference Iken, Echelmeyer, Harrison and FunkIken and others, 1993; Reference KrabillKrabill and others, 1999; Reference Mayer and HerzfeldMayer and Herzfeld, 2000, Reference Mayer and Herzfeld2001; Reference Johnson, Prescott and HughesJohnson and others, 2004; Reference Joughin, Abdalati and FahnestockJoughin and others, 2004; Reference Podlech and WeidickPodlech and Weidick, 2004; Reference ThomasThomas, 2004; Reference ThomasThomas and others, 2004; Reference Luckman and MurrayLuckman and Murray, 2005; Reference JoughinJoughin and others, 2008, Reference Joughin, Smith, Howat, Scambos and Moon2010, Reference Joughin2012; Reference Mayer and HerzfeldMayer and Herzfeld, 2008; Reference Nick, Vieli, Howat and JoughinNick and others, 2009; Reference Moon, Joughin, Smith and HowatMoon and others, 2012).

The objective of this paper is to investigate dynamic changes in Jakobshavn Isbræ, using an approach that combines structural glaciology and remote-sensing data analysis and is facilitated by mathematical characterization of generalized spatial surface roughness. Generalized spatial surface roughness is a mathematical concept that encompasses surface roughness and spatial structural properties of the ice surface, morphology and micro-morphology (short: surface roughness). This surface-roughness characterization serves to provide parameters that are related to ice flow, ice dynamics and interaction of the ice with bed topography.

The idea is to provide a method that facilitates mapping of geophysical parameters that relate to ice-dynamic change, using satellite and airborne altimeter data, and to apply this method to Jakobshavn Isbræ. First, we will demonstrate that the roughness approach can be applied to map dynamic provinces in fast and slow-moving ice and that several dynamic types can be distinguished. The focus of the paper is on ICESat (Ice Cloud and land Elevation Satellite) GLAS (Geoscience Laser Altimeter System) data analysis, and ATM (Airborne Topographic Mapper) data are included for high-resolution analysis and to validate GLAS data analysis. In addition, ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) data of the NASA Terra satellite are used for validation.

Second, the approach will then be applied to investigate elevation and roughness changes in Jakobshavn Isbræ south ice stream during rapid retreat, which occurred between 1999 and 2007, to contribute to solving the Jakobshavn problem: what is the cause of this rapid retreat? Possible causes of the observed changes include: (1) a fjord-glacier cycle (Reference Meier and PostMeier and Post, 1987; Reference PfefferPfeffer, 2007; Reference Vieli and NickVieli and Nick, 2011), (2) calving of the front and acceleration started in consequence of climatic warming in the Jakobshavn glacier region (Reference Box and SteffenBox and Steffen, 2001; Reference Steffen, Box and StraneoSteffen and Box, 2001; Reference Box, Yang and BromwichBox and others, 2009), (3) acceleration due to an internal change in the glacier’s dynamic system, and (4) increased calving and front retreat due to warming of the water in Ilulissat Isfjord, which then led to upward propagation of acceleration due to lack of back-holding force (Reference Rignot, Fenty, Menemenlis and XuRignot and others, 2012; Straneo and others, 2012; cf. Reference Truffer and EchelmeyerTruffer and Echelmeyer, 2003). Jakobshavn Isbræ is arguably the best-studied tidewater glacier on Earth, yet its recent retreat has not been fully attributed to any one of these possible causes, or to a combination.

2. Approach

The approach used here combines structural glaciology and analysis of generalized surface roughness. 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. 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 Liu and LohoefenerLiu, 2002; Reference GreveGreve, 2003). Crevassing is a form of brittle deformation that occurs when local forces exceed a threshold. Characterization of generalized spatial surface roughness is a mathematical approach that utilizes parameters derived from spatial statistical functions to capture spatial properties of a surface (Reference HerzfeldHerzfeld, 2008). This applies to crevassed and non-crevassed ice surfaces. Combining structural geology and mathematical roughness analysis allows us to derive deformation characteristics in fast-moving glaciers, as described in theory in Reference Herzfeld, Mayer and GreveHerzfeld and others (2004), and to map deformation provinces in surging and continuously fast-moving glaciers (Reference Herzfeld and MayerHerzfeld and Mayer, 1997; Reference HerzfeldHerzfeld 1998, Reference Herzfeld2008; Reference Herzfeld, Mayer, Feller and MimlerHerzfeld and others, 2000; Reference Mayer and HerzfeldMayer and Herzfeld, 2000, Reference Mayer and Herzfeld2001; Reference Herzfeld and ZahnerHerzfeld and Zahner, 2001; Reference HerzfeldHerzfeld, 2008; cf. Reference Vornberger and WhillansVornberger and Whillans, 1990; Reference MarmoMarmo and Wilson, 1998; Reference RistRist and others, 1999). Hence crevasse patterns and surface roughness can be used as a source of geophysical or glaciological information. For mapping and monitoring of dynamic changes in fast-flowing glaciers and ice streams in Greenland, Antarctica and other remote areas, it will be useful to be able to perform this from space. In this paper, we investigate to what extent ice surface roughness affects satellite and airborne laser altimeter data and, in turn, to what extent roughness characteristics derived from laser altimeter data can be employed to map dynamic provinces. The area of Jakobshavn Isbræ and its neighborhood provides a good test area for this approach, because many different crevasse types and roughness types are encountered in the fast-moving glacier and the adjacent slow-moving ice.

Two types of laser altimeter data are used: GLAS data collected during the ICESat mission (2003–09) (Reference Schutz, Zwally, Shuman, Hancock and MarzioSchutz and others, 2005; Reference ZwallyZwally and others, 2005) and ATM data collected from NASA aircraft (Reference Krabill, Thomas, Martin, Swift and FrederickKrabill and others, 1995). The ICESat ground tracks cross the fast-moving Jakobshavn Isbræ south ice stream at angles that sample the ice stream’s dynamical provinces almost perfectly, as will be derived below. ATM data that were simultaneously collected underflying ICESat in March 2006 and March 2007 will be employed to validate the GLAS data analysis and provide higher-resolution information.

Applications are (1) mapping ice-dynamic provinces or types in Jakobshavn Isbræ and its environs, (2) analysis of changes in elevation and surface roughness and (3) glacio-logical interpretation.

3. Dynamic Provinces of Jakobshavn IsbrÆ Determined Using Structural Glaciology

An overview of the dynamic provinces of Jakobshavn Isbræ is given in Figure 1, based on synthetic aperture radar (SAR) data collected with the European Remote-sensing Satellite ERS-2 on 22 and 25 May 1997 (Reference Mayer and HerzfeldMayer and Herzfeld, 2001). The analysis shown in Figure 1 uses a composition of two ERS-2 datasets: ERS-2 orbit 10958, descending frame 2205–25May1997; and ERS-2 orbit 10915, descending frame 2205–22May1997.

Fig. 1. Structural glaciological provinces of Jakobshavn Isbræ based on 1995 ERS-2 SAR data composition of (1) ERS-2 orbit 10958, descending frame 2205–25May1997 and (2) ERS-2 orbit 10915, descending frame 2205–22May1997 (data collected by European Space Agency). (From Reference Mayer and HerzfeldMayer and Herzfeld 2001.) The main ice stream, marked here by double lines, is Jakobshavn south ice stream; the north ice stream joins the south ice stream north of the easternmost land promontory seen in this image. The center line of Jakobshavn south ice stream is indicated in the lower part of the ice stream. The symbols indicate directions of extension.

A close-up of the fast-moving south ice stream allows a finer classification into five different crevasse classes, which are oriented parallel to the center line and occur symmetrical to the center line: (1) closed compressional and conjugate shear faults, (2) very smooth surface with pervasive shear deformation, (3) closed conjugate shear faults with one dominant direction curving into flow direction towards the center, (4) closed conjugate shear faults and (5) slightly open extensional crevasses in multiple directions (Fig. 2a). Zone (1) in Figure 2 contains the center line (marked in Fig. 1), and zone 5 is the shear margin, located outside the double lines that mark the fast-moving Jakobshavn south ice stream. A comparison with an ASTER image of the NASA Terra satellite, taken in May 2003, shows that the longitudinal zones of Jakobshavn south ice stream had remained structurally unchanged (Fig. 2b). Aerial photographic observations collected in 2005 indicated that the longitudinal zones of the south ice stream continued to stay the same (Reference Mayer and HerzfeldMayer and Herzfeld, 2008). Large melt ponds remained throughout the decade 1995–2005, but the number of small melt ponds increased.

Fig. 2. Longitudinal zonation of the focused core of Jakobshavn Isbræ south ice stream, 1996 and 2003. (a) Based on 1996 aerial photograph. From center (left) to margin (right), view in flow direction: (1) closed compressional and conjugate shear faults, (2) very smooth surface with pervasive shear deformation, (3) closed conjugate shear faults with one dominant direction curving into flow direction towards center, (4) closed conjugate shear faults, (5) slightly open extensional crevasses in multiple directions. Zones repeat symmetrically on the other side of the center line. (After Reference Mayer and HerzfeldMayer and Herzfeld, 2001, Fig. 5.) Photograph taken slightly downstream of the turn in the south ice stream. (b) Same zones indicated on top of 2003 Terra ASTER image (see Fig. 3); crossing of 1996 flight line is just to the left of the numbers. The zone names ‘shear margin’ and ‘center of ice stream’ have been added for general reference. Comparison of (a) and (b) shows that dynamic provinces have not changed.

The structural glaciological provinces of Jakobshavn Isbræ south ice stream are longitudinally aligned centered around the flowline, because Jakobshavn Isbræ follows a subglacial trough, which runs generally west–east (see Fig. 1 and compare with Figs 3 and 4 and analysis in Section 6). The trough causes the acceleration from slow-moving inland ice to fast-moving ice stream. Therefore, the first and seventh authors flew the 1996 observation tracks in directions normal to the ice stream’s flowline to sample surface properties of the structural provinces. Figure 2 documents these. ICESat GLAS tracks happen to be oriented and located almost perfectly to analyze spatial roughness characteristics of these longitudinally oriented provinces, because the tracks run approximately north–south, with the ascending and descending tracks forming a narrow angle. The track selected in Figure 4 crosses the south ice stream just down-glacier of the bend and thus is located close to one of our 1996 flight tracks and was underflown with the ATM instrument in March 2005 and March 2006. The coincidence of these three datasets facilitates the analysis in this paper.

Fig. 3. Generalized roughness measures for Jakobshavn Isbræ south derived from ICESat GLAS data and ATM data. (a) Geostatistical pond _var and pond _res parameters allow us to map ice-dynamic roughness provinces of Jakobshavn Isbræ’s south ice stream. GLAS L3I data, November 2007. From left to right: geostatistical classification parameters pond _var, pond _res. Tracks offset to show different parameters; left track is in true location. UTM coordinates. (b) Ice-surface roughness mapped from full-resolution ATM data using geostatistical pond _var parameter indicates ice-dynamic provinces of Jakobshavn Isbræ’s south ice stream. 2009 IceBridge ATM data. (a, b) Background ASTER data May 2003, NASA Terra satellite.

Fig. 4. Ice-surface roughness from GLAS data, subglacial topography from CReSIS MCoRDS data and crevassing and surface image from Terra ASTER data for Jakobshavn Isbræ region. GLAS L3I data from November 2007; roughness measure pond _res. Subglacial topography regridded and contoured from CReSIS MCoRDS data. Background ASTER data May 2003. UTM coordinates. Value of pond _res parameter plotted normally to flight track (red) demonstrates relationship to crevasse provinces and subglacial topography.

4. Mathematical Background: Geostatistical Classification Based On ICESat Glas And Atm Altimeter Data

The analysis in this paper applies the geostatistical characterization and classification method as introduced in Reference HerzfeldHerzfeld (2002, Reference Herzfeld2008). In the interest of brevity, we limit the mathematical description to the functions and parameters employed here. A mathematical treatment and comparative analysis of generalized roughness analysis for laser altimeter data, especially ICESat GLAS data, using waveform parameters and geostatistical classification parameters is deferred to a later paper. To mathematically summarize spatial surface roughness in a given area, a vario function is formed. In this paper we use the first-order vario function and the residual vario function, motivated and defined as follows (Eqn (1)). If we are interested in characterizing the surface structure in a given area, then the actual value of surface elevation at each location is not as relevant as the parameters that describe the morphology more generally for the whole study area. Therefore, we form differences of measurement values and average over all points that have the same common distance (or distance and direction), creating vario functions, defined by

for pairs of points (x_i, z (x_i )), (x_i ) h, z (x_i ) h)) D, where D is a region in R² (case of survey profiles) or R³ (case of survey areas) and n is the number of pairs separated by h; the distance value h is also termed lag. The function v ₁(h) is called the first-order vario function. This function exists always and has a finite value, because only finitely many data points enter the calculation. In situations where a regional trend or a local drift underlies the data, the residual vario function is often more useful to analyze roughness. With

the residual vario function is defined as:

If average surface elevation is constant for distance lag h, then m (h) = 0 and res ₁(h) = v ₁(h), and if the surface is stationary in the statistical sense, then res ₁ ≡ v ₁. The latter is heuristically not the case for glaciers and ice streams, therefore use of the residual vario function is often more appropriate.

The pond parameter is defined as the maximal value of the vario function. The higher the spatial variability, the higher is the pond value. In our application, the vario function is evaluated for a moving window along GLAS or ATM tracks, using a maximum lag distance of 0.8 times the window length, then the pond parameter is the maximum of the vario-function values in that window. The terms pond _var and pond _res designate pond determined from the first-order vario function v ₁ and the residual vario function res₁, respectively. For ice regions with approximately constant average height (e.g. some sub-areas of the Greenland inland ice), the pond _var parameter may be better suited, whereas for areas with sloping topography (e.g. drainage basins of outlet glaciers), the pond _res may be more appropriate. We examine this further in Sections 6 and 9.

5. Instrumentation and Data: Icesat Glas and Atm

The GLAS instrument aboard ICESat is a pulse-limited laser altimeter operated at 40 Hz, 532 nm and 1064 nm using a doubler crystal. The 1064 nm (infrared) laser channel is designed for measurements of surface altimetry and heights of dense clouds, and the 532 nm (green) lidar channel for measurements of the vertical distribution of clouds and aerosols. Altimeter data are derived using waveform fitting of the returned signal, geolocation and inversion of two-way travel time for range determination. With the help of precise orbit determination, laser-pointing correction and several other corrections, elevation of the ice surface is calculated (Reference Schutz, Zwally, Shuman, Hancock and MarzioSchutz and others, 2005; Reference Herzfeld, Bride, Zwally and MarzioHerzfeld and others, 2008). All satellite altimeter data follow ground tracks with points (65 m diameter footprints for GLAS) spaced along-track (172 m for GLAS) and large gaps in between, with spacing dependent on latitude. For most of the ICESat mission, a 33 day subcycle of a 91 day repeat orbit was used. Following failure of the first of the lasers (labeled L1, L2, L3) early in the mission, a 33 day observation cycle was carried out three times per year (February or March, May, October–November), and later twice per year (February or March and October–November). Track spacing at 70° latitude is 10.1 km for the designed 91 day orbit and 29.7 km for the actual 33 day subcycle of the 91 day repeat orbit. Note that, other than a radar altimeter, the laser altimeter signal does not penetrate optically dense clouds (but thin clouds); this results in significant losses of surface returns, especially in the Arctic (Reference Spinhirne, Palm and HartSpinhirne and others, 2005).

Small variabilities exist in the ground-track locations for the repeat tracks, and along-track locations of points also vary between the 33 day missions. Elevation and pond difference values for data from two different tracks are only formed if pairs of points are found in a small window with a small tolerance. This reduces the effects of slopes (in any direction) especially on the elevation differences and differences in pond _var; parameter pond _res is not affected by surface slope.

The ATM, in contrast, uses a scanning laser technique, which results in ~10 000 points per near-circular ellipsoidal scan (Reference Krabill, Thomas, Martin, Swift and FrederickKrabill and others, 1995). Because the aircraft carrying the laser moves forward, the ground pattern of observation is a spiral. Points are spaced 3–4 m apart at near-nadir locations and 5–6 m apart along the outer margins of the swath. Our analysis uses full-resolution data. Application of the classification method requires prior aggregation of data into location-referenced bins.

6. Roughness and Dynamic Provinces Of Jakobshavn IsbrÆ From Glas and Atm Data

Geostatistical classification parameters pond _var and pond _res are calculated for ICESat GLAS data of all missions, using along-track moving windows, and for the much finer ATM data. The result is shown in Figure 3a for GLAS L3I data (November 2007) and in Figure 3b for ATM data (2009), superimposed on a color composite of the high-resolution (15 m pixel) bands 1,2,3 of ASTER data from the Terra satellite. The ASTER data show the dynamic provinces in the Jakobshavn region, such as the central fast-moving south ice stream, the wide shear margins north and south of that, undulating ice surface and smooth ice with melt ponds in the southern part of the area in Figure 3a, and large crevasses that cause disintegration of the lower Jakobshavn Isbræ. The spatial distribution of values of the pond _var parameter matches the location of the ice-dynamic provinces. The same is true for the pond _res parameter, but values are generally lower, because the residual with respect to the drift is formed (Eqn (3)). Highest values occur over the heavily crevassed shear zones north and south of the south ice stream. A comparison of Figures 3a and 2 shows that of the five longitudinal crevasse zones identified in photographs, the central one (zone 1) and the neighboring, smoothest zone 2 with closed crevasses are resolved in the ASTER data near the bend in the south ice stream (564000 W, 7671000 N). Zones 3 and 4 appear similar, and the shear margin is clearly visible. At the location where the GLAS track crosses, the focused core does not show zones 1 and 2 separately (cf. Fig. 2b). The pond _var and pond _res values from GLAS data do not resolve the longitudinal bands in the central core, but the transverse increase in roughness from the center out is well resolved. The analysis is of course limited by the resolution of the data.

Application of the classification parameters to the much finer ATM data (Fig. 3b) reflects that the fast-moving ice stream is smoothest in longitudinal bands to the side of the center and then increases in several steps towards the shear margin, where highest pond _var values are derived. The region south of the southern shear margin shows the overall lowest roughness values; this is an area of uncrevassed ice with melt ponds.

7. Ice-Surface Roughness and Subglacial Topography

Attribution of acceleration in ice streams and resultant crevassing generally falls into two categories, dynamically and topographically induced. Many outlet glaciers of the Greenland ice sheet follow troughs in the bedrock, evident in surface roughness as mappable by the above-mentioned parameters from GLAS altimeter data.

The fast glacier flow in Jakobshavn Isbræ south ice stream is caused by a subglacial trough. The combined effects of convergent flow of ice into the trough and increased basal sliding (due to thicker ice and hence higher pressure at the base of the trough) lead to spatially accelerated flow of the ice stream, and the faster flow causes heavy crevassing. In order to analyze the spatial relationships between subglacial bed topography, ice-surface roughness and crevassing, we reanalyzed bed topographic data and created a high-resolution map, co-referenced with the ASTER and ICESat satellite data (Fig. 4). The bed-topographic data were collected by the Center for Remote Sensing of Ice Sheets (CReSIS), University of Kansas, USA, using the Multichannel Coherent Radar Depth Sounder (MCoRDS) (Reference GogineniGogineni and others, 2001; Lohoefener, 2006; Li, 2009). Figure 4 shows bed topography, crevassing from 2003 ASTER data and surface roughness from GLAS data (2007), quantified by the pond _res parameter. Figure 4 demonstrates that the fast ice movement of Jakobshavn Isbræ south ice stream, which is indicated by heavy crevassing that extends >80 km into the ice sheet, follows the location of a subglacial trough. The composite of radar data, image data and roughness parameters provides evidence that the spatial acceleration of Jakobshavn Isbræ is caused by the presence of the trough, which is not the case for the north ice stream. Spatial surface roughness properties calculated from GLAS data allow us to more closely relate spatial acceleration of the ice to dynamic provinces, as analyzed in more detail in Section 6. While highest velocities are known to occur in the center of the ice stream, the roughest terrain is found in the shear margins of the south ice stream and in the area where the north ice stream accelerates. Overall, the fast-moving ice stream has higher surface roughness than the surrounding inland ice, which moves at 0.3 m a^-1. Highest roughness values occurring over the shear margins are related to the spatial velocity gradient. The fastest-moving center of the ice stream lines up with the center of the trough and has an intermediate roughness. The inner edge of the shear margins coincides in location with the outer edge of the subglacial trough in the CReSIS data.

Surface processes such as wind and melting also influence generalized spatial ice surface roughness, but occur on a smaller spatial scale than crevassing, to ~10 m, and are morphologically less significant, resulting in much smaller pond _var and pond _res parameters outside crevassed areas (Figs 3 and 4; Reference Herzfeld, Mayer, Feller and MimlerHerzfeld and others, 1999, Reference Herzfeld, Mayer, Feller and Mimler2000). In crevassed regions such as the lower reaches of Jakobshavn Isbræ and other fast-flowing Greenland outlet glaciers, ice dynamics is the main source of surface roughness.

8. Change In Elevation In Jakobshavn South Ice Stream

Since 1999 Jakobshavn Isbræ has exhibited a rapid retreat of the ice front and surface lowering. A goal of this paper is to contribute to analyzing and understanding these changes. We apply the methods described in the previous sections and investigate the variables (1) elevation change, (2) change in surface roughness as quantified by the pond parameters, and (3) subglacial topography. For clarity, the analysis is first carried out for two points in time (March 2006 and March 2007; for exact mission intervals see Table 1) during rapid retreat and surface lowering of Jakobshavn Isbræ. The analysis of change detection is based on ICESat GLAS data from March 2006 (GLAS L3E) and March 2007 (GLAS L3H). These two mission time frames were selected because data are available for the same season and the same track locations and data density along-track is very good, which allows a comparative analysis. In Section 10 the analysis will be extended to time series.

Differences between March 2007 and November 2006 altimeter data indicate a surface lowering of 10–20 m in 1 year over the south ice stream (Fig. 5a). The highest values occur at the center of the south ice stream. Elevation differences gradually decrease over the northern shear zone to lower values of ~4 m in the slow-moving ice area north of Jakobshavn Isbræ south ice stream. Elevation difference values also decrease to the south of the south ice stream, with a high variability over the southern shear margin. Local variability may be attributed to a forward-scattering effect or multiple returns over crevassed areas. However, all elevation difference values are negative. Considering that spatial variability in elevation change is likely due to crevassing, the amount of surface lowering generally increases toward the center of the south ice stream.

Fig. 5. Jakobshavn Isbræ south surface lowering and surface roughness changes between March 2006 (GLAS L3E) and March 2007 (GLAS L3H). (a) Surface lowering rate was ~15 m a^-1, varying locally to 20 m a^-1. (b) Surface roughness changes, calculated as differences in parameters pond _var (left) and pond _res (right) between March 2006 (GLAS L3E) and March 2007 (GLAS L3H). Results for pond _res offset to the right for visualization purposes; correct location of the track is that of the left track. (a, b) Background ASTER data May 2003, NASA Terra satellite.

9. Changes in Roughness in Jakobshavn South Ice Stream

Changes in surface roughness in the same time frame of March 2006 (GLAS L3E) to March 2007 (GLAS L3H) are shown in Figure 5b. Surface-roughness change in time is calculated as

for two points in time t ₁, t ₂. Surface-roughness changes as indicated by the pond _var parameter are small throughout the study area. However, pond _var shows variability especially in areas of crevasse fields that may be associated with undulations in surface topography and with the sloping margins of the south ice stream (the ice stream is lower than the surrounding inland ice; see Section 10). As expected from Eqns (2) and (3), the values of pond _var are affected by the local trend of the underlying topography. The fact that the local trend is nonzero is also visible in the fact that pond _var and pond _res are locally different (Fig. 3). Hence the pond _res parameter is used in the following analysis of roughness changes. Spatial variability is very low for pond _res. In conclusion, during the time of rapid elevation drop, the spatial and temporal distribution of the roughness provinces did not change. The implications of this observation are discussed in Section 11.

10. Time Series of Elevation Changes and Roughness Changes

The objective of this section is to investigate elevation and roughness changes over several years. The ICESat mission lasted from 2003 to 2009, but it is difficult to find time frames with data that are useful for comparison. Limiting factors are the following: Data were collected only for two to three 33 day intervals per year, termed missions (Table 1).

While the basic ground track pattern was repeated in each of those operation times, off-pointing manoeuvres in the area of lower Jakobshavn Isbræ limit the coincidence of tracks in the same locations. Furthermore, lidar data are generally affected by dense cloud cover, which further limits the number of comparisons between tracks that can be used for change detection by along-track data evaluation (Reference Spinhirne, Palm and HartSpinhirne and others, 2005). For the purpose of this paper, we have analyzed all GLAS data collected in the study area during the entire ICESat time frame 2003–09 and formed elevation and roughness differences for all possible pairs of missions (Table 1). Pairs of GLAS data that yielded sufficiently many points are investigated in the sequel.

Data from missions L3H and L3E were selected as primaries, because footprints for these two tracks match well in location, and hence elevation and roughness differences can be calculated for an entire across-flow profile. To observe the retreat of Jakobshavn Isbræ, off-pointing manoeuvres of tracks whose original position was near the front of the glacier to the front of the ice stream were frequently undertaken; this yields the differences in location between the tracks that cross the lower ice stream (track 1 is in location 1a or 1b) in Figure 6 and is also the reason that the track seen at (E 563000/N 7660000) to (E 553000/N 7700000) in Figure 6b appears to be missing in Figure 5. We use the names track 1A for the track at (E 555000/N 7660000) to (E 545000/N 7700000) in Figure 6a, track 1B for the track at (E 563000/N 7660000) to (E 553000/N 7700000) in Figure 6b and track 2 for the track that is common to L3H and L3E and crosses the center of Jakobshavn south ice stream near its southernmost point (Figs 5, 6a and b and 7a and b).

Fig. 6. Jakobshavn Isbræ surface elevation change per year. Elevation change per year is calculated by subtracting elevations from two overlapping tracks and dividing by time (months/12), for all mission pairs possible that include L3H or L3E. (a) Relative to March 2006 (GLAS L3E): differences per year are formed for L3E – {L3B, L3C, L3D, L3H, L2E}. Track 1A near calving front (difference pairs are L3E – {L3B, L3C}), track 2 crosses near bend in south Isbræ. (b) Relative to March 2007 (GLAS L3H): differences per year are formed for L3H – {L2B, L3D, L3E, L3I, L3J}. Track 1B and track 2 cross near bend in south Isbræ. (a, b) Results are offset to the right for visualization purposes. Actual location of GLAS ground tracks coincides with the location of the leftmost track in each set, and with the tracks in Figure 5 (note that not all five differences exist in each location; for details, see Figs 8–11). Background ASTER data May 2003, NASA Terra satellite.

Figure 6 shows surface elevation change rates in map view, given as elevation change per year with rates accurate to the number of months. Figure 6a shows surface elevation change per year relative to March 2006 (L3E), and Figure 6b the same relative to March 2007 (GLAS L3H). Elevations for L3B, L3C, L3D, L3H, L2E were subtracted from elevations for L3E, and elevations for L2B, L3D, L3E, L3I, L3J were subtracted from elevations for L3H. The elevation differences were divided by the time (months/12) between the two missions in the pair. Elevation differences were only formed if matching elevation points were found in two tracks, with a small tolerance in location of the points. This avoids matching points across large differences in location, where slope effects in any direction might bias an elevation difference. Surface lowering is observed not only during the year March 2006–March 2007 (Fig. 5), but occurred monotonously throughout the observation time frame at rates of 10–15 m a^-1. For the tracks that cross the south ice stream near its bend (tracks 1B and 2), elevation losses are largest over the ice stream itself and the proximal regions of the northern and southern shear margins. For track 1A (Fig. 6a), this statement holds in principle, but spatial variability is much higher than farther up-glacier. Track 1A pointed at locations near the ice front (note that the correct location is that of the left plotted track) in an area affected by heavy calving in 2005. The track differences seen are from L3B (February 2005) to L3E (February 2006) and L3C (May 2005) to L3E (February 2006), during which the north ice stream retreated. The time interval between L3I (October 2007) and L3H (March 2007) appears to be an outlier of exceptionally high change rates in the interior of the ice stream. This is further investigated below.

The corresponding change rates have been calculated for generalized ice surface roughness measured by the parameter pond _res, shown in map view in Figure 7. As seen in Figure 7, surface roughness change rates are close to zero for any pair of tracks in the interior of the south ice stream (tracks 1B and 2). Exceptions are a few isolated locations at the northern shear margin in track 1B (Fig. 7b), which may be attributable to changes in the location of individual large crevasses, possibly caused by the drawdown of the central ice stream relative to the margin. These results are in principle the same as those derived for the 2006–2007 comparison in Section 9. The track that crosses near the calving front (track 1A, Fig. 7a) shows roughness changes locally in the area where the north ice stream retreated in 2005, according to field observations of the first author. In this area, Jakobshavn Isbræ was not a connected, crevassed ice stream any more, but a semi-connected melange of crevassed ice pieces where the calving front had become a transitional area between ice stream and fjord with large icebergs.

Fig. 7. Jakobshavn Isbræ surface roughness change per year. Surface roughness change per year is calculated by subtracting surface roughness values (pond _res) from two overlapping tracks and dividing by time (months/12), for all mission pairs possible that include L3E or L3H. (a) Relative to March 2006 (GLAS L3E): differences per year are formed for L3E – {L3B, L3C, L3D, L3H, L2E}. Track 1A near calving front (difference pairs are L3E – {L3B, L3C}), track 2 crosses near bend in south Isbræ. (b) Relative to March 2007 (GLAS L3H): differences per year are formed for L3H – {L2B, L3D, L3E, L3I, L3J}. Track 1B and track 2 cross near the bend in south Isbræ. (a, b) Results are offset to the right for visualization purposes. Actual location of GLAS ground tracks coincides with the location of the leftmost track in each set, also with the tracks in Figure 5 (note that not all five differences exist in each location; for details, see Figs 8–11). Background ASTER data May 2003, NASA Terra satellite.

Elevation and roughness changes are investigated in more detail using values plotted along tracks in Figures 8–11. Elevation changes (Figs 8a, 9a, 10a and 11a) are plotted such that a continued surface lowering through time is pictured by a drop in elevation from one panel to the next. The top panel for each figure shows the elevation profile that follows the tracks and hence runs approximately across-flow for the ice stream. Note that there are fewer roughness-difference values than elevation-difference values, because roughness is calculated for a moving window. Running means of elevation and roughness values over 1000 m along-track windows are calculated for every 100 m along-track (red in Figs 8–11).

From front up-glacier:

Track 1A crosses near the front of both the north and south ice streams. Only the earliest pairs have data differences (the front retreated later beyond this area). The inner margin, where north and south ice streams converge, is near N 7677000. Elevation loss was largest over the south ice stream. The elevation loss in the north ice stream extends up into the entire crevassed region and along the margin of the north ice stream (to N 7690000). The pond _res difference plots show that roughness variability is highest along the shear margins and the inner confluence of the north and south ice streams, which is explained by crevassing. The laser signal can hit the bottoms of crevasses, which change their location rapidly. Notably, the area of elevation loss extends across the entire north and south ice streams, except for a thickening at N 7680000 in February 2005– February 2006 that is not present in the May 2005– February 2006 difference. Roughness changes do not follow the trend of elevation changes, and overall roughness changes are very low.

Track 1B crosses just upstream of the so-called ‘rumple’ (Reference Echelmeyer, Clarke and HarrisonEchelmeyer and others, 1991), an area of heavy crevassing upstream of the confluence of the north and south ice streams. Elevation is lowest (for L3H in March 2007) immediately north of the south ice stream. Figure 10a shows that surface lowering at this crossing is much higher over the south ice stream than over the north ice stream. Surface lowering extends into the southern shear margin of the south ice stream, but is not as strong as over the main fast-moving provinces (boundary at ~N 7670000). Elevation loss was 40–60 m from L2B (February 2004) to L3H (March 2007), 20–30 m between L3H (March 2007) and L3J (October 2007), and only 15–20 m between L3H (March 2007) and L3J (March 2008). Hence a small elevation increase must have occurred between October 2007 and March 2008 in the south ice stream, which indicates thickening in winter.

North of N 7680000, track 1B crosses only several shear zones east of the north ice stream, not the north ice stream itself. Elevation losses are much lower than for the south and north ice streams (Fig. 8) and vary locally.

Fig. 8. Jakobshavn Isbræ surface elevation and roughness changes. Surface roughness changes are calculated as differences in parameter pond _res. (a) Surface elevation profile and elevation differences for L3E – track 1A. (b) Surface roughness profile and roughness differences for L3E – track 1A. In (a, b), middle two panels, black dots show differences from pairs of measurements, and red dots show differences of running means of measurements (where means are calculated over 1000 m windows for every 100 m along-track). Bottom panels show elevation and pond _res differences, formed by subtracting 1000 m running means (L3C – L3B).

The most interesting track is track 2, which traverses the upper Jakobshavn Isbræ south ice stream near the turn in direction (Figs 9 and 11) (note that the leftmost changerate track indicates the correct location of a GLAS track). The following landmarks or boundaries of glaciological provinces are: N 7663000 – southern boundary of southern shear margin of south ice stream; N 7667000 – southern shear margin of south ice stream/south ice stream; N 7668000–7669000 – center of south ice stream; N 7673000 – south ice stream/northern shear margin of south ice stream; N 7678000 – northern boundary of northern shear margin of south ice stream.

Fig. 9. Jakobshavn Isbræ surface elevation and roughness changes. Surface roughness changes are calculated as differences in parameter pond _res. (a) Surface elevation profile and elevation differences for L3E – track 2. (b) Surface roughness profile and roughness differences for L3E – track 2. In (a, b), second to second-to-last panels, black dots show differences from pairs of measurements, and red dots show differences of running means of measurements (where means are calculated over 1000 m windows for every 100 m along-track). Bottom panels show elevation and pond _res differences through time 2005–2009, formed by subtracting 1000 m running means (L3C, L3D, L3H, L2E – L3B).

Fig. 10. Jakobshavn Isbræ surface elevation and roughness changes. Surface roughness changes are calculated as differences in parameter pond _res. (a) Surface elevation profile and elevation differences for L3H – track 1B. (b) Surface roughness profile and roughness differences for L3H – track 1B. In (a, b), second to second-to-last panels, black dots show differences from pairs of measurements, and red dots show differences of running means of measurements (where means are calculated over 1000 m windows for every 100 m along-track). Bottom panels show elevation and pond _res differences through time 2004–2008, formed by subtracting 1000 m running means (L3I, L3J – L2B).

Fig. 11. Jakobshavn Isbræ surface elevation and roughness changes. Surface roughness changes are calculated as differences in parameter pond _res. (a) Surface elevation profile and elevation differences for L3H – track 2. (b) Surface roughness profile and roughness differences for L3H – track 2. In (a, b), second to fourth panels, black dots show differences from pairs of measurements, and red dots show differences of running means of measurements (where means are calculated over 1000 m windows for every 100 m along-track). Bottom panels show elevation and pond _res differences through time 2004–2008, formed by subtracting 1000 m running means (L3D, L3E, L3J – L2B).

The series of plots of difference values demonstrates that L3H–L3E is the best data pair (used in Sections 8 and 9). The values over the south ice stream are intermittent, but indicate that surface lowering continues throughout the observation time frame, as indicated in the bottom panels of Figures 9 and 11. Between L3B (February 2005) and L3E (February 2006), surface lowering over the central ice stream was 15–20 m, ~30 m between L3B and L3H (March 2007) and 40–60 m between L3B (February 2005) and L2E (March 2009), an average of 10–15 m a^-1, varying locally. The area north of the ice stream also lost elevation, ~10 m in 4 years (2.5 m a^-1; Fig. 9). An exception to the general surface lowering over the ice stream is a thickening in winter between October 2005 and February 2006 (Fig. 11), which repeated in 2007/08, as noted above, and corresponds to a slowing of the flow velocity reported in Reference JoughinJoughin and others (2012). Thickening is limited to the ice stream itself, while the area north of the ice stream lost elevation during the same winter.

Roughness change over the same time frames is again insignificant and does not show a trend over time.

Largest variabilities in roughness occur over the shear margins. Figure 11 indicates essentially the same process of surface lowering and a lack of roughness change, for the years 2004–08. Surface lowering is 40–60 m in 4 years, and 10–15 m a^-1.

In summary, the detailed analysis has provided more accurate values on spatial and temporal variability of elevation change, which is 10–15 m a^-1 over the south ice stream, and indicated a winter thickening. Roughness shows a variability in the shear margins, but did not change during the period of drastic surface lowering. A comparison of roughness and elevation change comes qualitatively to the same result for the analysis of time series as for the analysis of the pair of L3H–L3E datasets.

11. Summary of Results and Interpretation: Attribution of Changes During Rapid Retreat of Jakobshavn IsbrÆ

The analyses in the previous sections led to the following results:

1. For glaciers and ice streams generalized spatial surface roughness (at scales larger than 10 m) is an indicator of ice dynamics.

2. Surface roughness may be quantified by the pond _var and pond _res parameters as derived from ICESat GLAS and ATM data. The pond _res parameter is better suited for analysis of ice-stream regions because of local slopes in surface elevation.

3. Surface lowering is largest over the south ice stream, above the center of the trough.

4. While surface elevation decreased by ~10–15 m a^-1 between 2004 and 2009, surface roughness remained almost unchanged.

5. Spatial acceleration is caused by the existence of a subglacial trough. The trough in the bedrock morphology leads to a convergence of flow, concentration of ice in the trough and faster flow due to surface slope in the longitudinal direction of the trough. These processes lead to a spatial sequence of crevasse provinces arranged in the transverse direction.

We now use the above results to investigate what may have caused the observed acceleration in Jakobshavn Isbræ. Possible causes include: (i) an internal dynamic change in the glacier system, such as a kinematic wave, (ii) a climatically caused acceleration in the upstream glacier system (faster flow following warming of ice), (iii) front retreat as part of a fjord-glacier cycle (Clarke, 1987; Reference Meier and PostMeier and Post, 1987), (iv) front retreat because of fjord-water warming (Holland and others, 2008; Reference Rignot, Fenty, Menemenlis and XuRignot and others, 2012; Straneo and others, 2012), and (v) increased basal sliding due to increased water at the base of the glacier (cf. Reference GreveGreve, 2005; Reference Greve and BlatterGreve and Blatter, 2009; Reference Herzfeld, Fastook, Greve, Donald, Wallin and ChenHerzfeld and others, 2012).

Both (i) and (ii) involve dynamic changes in the interior of the glacier, whereas (iii) and (iv) are initiated at the calving front. We analyze which of the possible causes are consistent with our results.

The detailed spatio-temporal analysis, based on all pairs of ICESat GLAS datasets with matching points, has provided accurate values on spatial and temporal variability of elevation change and roughness change. Overall, elevation change was ~10–15 m a^-1 surface lowering over the southern branch of Jakobshavn Isbræ, with an indication of some winter thickening. The time series analysis showed that this surface change pattern persisted between 2004 and 2009 over the south ice stream.

Roughness parameters such as pond _var and pond _res provide a proxy mapping of structural glaciological provinces, due to the effects of crevasse types or crevassing intensity, and hence dynamic provinces. Since roughness did not change, we conclude that the distribution of the dynamic provinces remained essentially the same. The dynamic provinces are largely caused by the several glaciological effects on ice flow associated with movement of ice into and down the subglacial trough (in the area around the location of the GLAS track, which is ~15 km from the 2003 calving front). Given the combined observations of acceleration of glacier flow and significant surface lowering, but no detected changes in roughness as measured by our parametric approach, an internal dynamic change in the glacier system is not likely a cause for the surface lowering, front retreat or acceleration, as this would result in a change in roughness provinces in the central ice stream.

At this point we take an excursion to ascertain the validity of the conclusion relating change in roughness, measured by the pond parameter, and change in crevasse types, as a means to rule out an internal dynamic change in the glacier system. An internal dynamic change would result in new gradients in the kinematic forces that cause crevassing. The simplest, first case is that crevassing affects previously uncrevassed ice; in this case, spatial variability in an image or altimeter dataset increases and hence the pond parameter increases. A second case is that a new force affects previously crevassed ice; the new force field is not linearly dependent on the existing force field; new crevasses form, overprinting the existing crevasse fields; spatial variability increases and hence the pond parameter increases. The remaining, third case is that a new force affects previously crevassed ice and the new force field is linearly dependent on the existing force field; this is the case of a simple acceleration in the same direction as the existing flow, which for a previously crevassed fast-flowing glacier may (or may not) change the spacing of crevasses; in this case the pond parameter remains the same (because the pond parameter is defined as the maximal value for all distance classes in the vario function or residual vario function). The third case is the situation for Jakobshavn Isbræ. For examples illustrating changes in surface crevasse types caused by changes in the kinematic forces during a glacier surge, the reader is referred to Reference Mayer and HerzfeldMayer and Herzfeld (2000) and Reference HerzfeldHerzfeld (1998). The argument in this section holds more generally than for surges, namely for any sequences of kinematic forces that cause or change crevassing.

The interpretation – that an internal dynamic change in the glacier system is not (or not likely) a cause for the surface lowering, front retreat or acceleration – does not contradict the changes observed near the calving front of Jakobshavn Isbræ during rapid retreat (for a description of the structural changes during retreat see Reference Mayer and HerzfeldMayer and Herzfeld, 2008). Instead, the results of the analyses support the hypothesis of a climatically induced surface lowering and are consistent with a front retreat as part of a fjord-glacier cycle or following warming of fjord water, as described by Holland and others (2008), Reference Rignot, Fenty, Menemenlis and XuRignot and others (2012) and Straneo and others (2012). Increasing temperatures over the Greenland ice sheet as observed by Reference Box and SteffenBox and Steffen (2001), Reference Steffen, Box and StraneoSteffen and Box (2001) and Reference Box, Yang and BromwichBox and others (2009) have also been considered a likely cause for surface lowering in the Jakobshavn region, but this signal may not be strong enough to cause the observed acceleration. Results of our analysis of the tracks that cross the interior of the south ice stream are consistent with both increasing temperatures and warming fjord water as initiators of acceleration. The fact that surface roughness of the track near the calving front changes over time suggests an initiation of the acceleration processes near the front and therefore favors the ice–ocean-interaction line of thought.

The observations and analyses are also consistent with the hypothesis of a front retreat that triggers a disintegration of the glacier near the front, which propagates up-glacier due to lack of back-holding force to the area of the confluence of the north and south ice streams and causes further acceleration in the ice stream. The classic glaciological explanation of the sudden retreat of a fjord glacier lies in the fjord-glacier cycle, according to which a fjord glacier retreats over a long time, until it reaches a minimal extent and then advances again (cf. Clarke, 1987; Reference Meier and PostMeier and Post, 1987; Reference PfefferPfeffer, 2007; Reference Vieli and NickVieli and Nick, 2011). The duration of a fjord cycle is typically on the order of 1 ka. This possibility is supported by the observations that Jakovshavn Isbræ has been retreating slowly throughout recent decades, experienced rapid retreat between 1999 and 2007 and then came to a halt in 2007, as the front reached the confluence of the north and south ice streams, which was considered a possible natural limit for retreat. However, the glacier has meanwhile retreated beyond this line. Recent research indicates that warming of the fjord water that interacts with the floating tongue of the glacier may cause rapid retreat as observed at Jakobshavn Isbræ (Holland and others, 2008; Reference Rignot, Fenty, Menemenlis and XuRignot and others, 2012; Straneo and others, 2012). A number of explanations of physical processes that may have caused the retreat and acceleration have been examined in Reference JoughinJoughin and others (2012), using models. This paper demonstrates how characterization of spatial surface roughness can be derived from satellite and airborne laser altimeter data and applied to narrow down possible physical processes.

12. Outlook To Icesat-2

Time-dependent analysis of the roughness characteristics in the along-flow direction of the ice stream is expected to reveal changes; however, the limited spatial coverage by GLAS data does not allow such investigations. As a consequence of the higher resolution of the ATM data compared to GLAS data, application of the pond parameter resolves the across-track sequence of dynamic provinces much better. Future satellite ICESat-2, to be launched by NASA in 2016, will carry a micro-pulse photon-counting altimeter with six beams in the across-track direction, whose expected resolution is 0.7 m in the along-track direction (Reference AbdalatiAbdalati and others, 2008, Reference Abdalati2010; Reference HerzfeldHerzfeld and others, 2014). Spatially and temporally well-resolved observation of ice-surface elevation and elevation change, especially in the highly dynamic regions of the outlet glaciers, has been identified as critical for understanding the changes in the Greenland and Antarctic ice sheets, the expected volume loss from the ice sheets, and the impending sea-level rise as a consequence of mass loss. The analyses undertaken here suggest that ICESat-2 data can be used to map glacio-dynamic provinces in detail. From the order-of-magnitude better spatial and temporal resolution of future ICESat-2 data, compared to GLAS data, more information can be expected in two ways: (1) simply due to improved mapping capabilities of the sensor, using the same analysis methods as for ICESat data, and (2) through application of multiple parameters (for examples, see Reference HerzfeldHerzfeld, 2008) that capture spacing of crevasses, melt pond distribution, deformation characteristics and other glaciologic variables.