Synthesizing multiple remote-sensing techniques for subglacial hydrologic mapping: application to a lake system beneath MacAyeal Ice Stream, West Antarctica

Abstract

We present an analysis of the active hydrologic system of MacAyeal Ice Stream (MacIS), West Antarctica, from a synthesis of multiple remote-sensing techniques: satellite laser altimetry; satellite image differencing; and hydrologic potential mapping (using a satellite-derived DEM and a bedrock DEM from airborne radio-echo sounding). Combining these techniques augments the information provided by each one individually, and allows us to develop a protocol for studying subglacial hydrologic systems in a holistic manner. Our study reveals five large active subglacial lakes under MacIS, the largest of which undergoes volume changes of at least 1.0 km³. We discuss the hydrologic properties of this system and present evidence for links between the lakes. At least three of the lakes are co-located with sticky spots, i.e. regions of high local basal shear stress. We also find evidence for surface elevation changes due to ice-dynamic effects (not just water movement) caused by changes in basal resistance. Lastly, we show that satellite radar altimetry is of limited use for monitoring lake activity on fast-flowing ice streams with surfaces that undulate on ∼10 km length scales.

1. Introduction

The recent realization that subglacial activity causes significant surface displacements that are detectable by satellite has led to rapid developments in this field. Over the past 5 years, satellite observations (interferometric synthetic aperture radar (InSAR), radar altimetry (RA) and laser altimetry) have radically altered our picture of subglacial water systems under the Antarctic ice sheet (Gray and others, 2005; Wingham and others, 2006; Fricker and others 2007; Smith and others 2009). As recently as 2005, all known subglacial lakes were located near ice divides, and were thought to be isolated and relatively inactive (Kapitsa and others, 1996; Siegert and others, 2005). By inference from rapid (on the order of months), localized surface elevation changes, subglacial water has been observed to move rapidly in and out of reservoir areas below the ice. Refill and drainage can occur via steady flow or episodic floods, and in several cases (e.g. Wingham and others, 2006; Fricker and Scambos, 2009) coordinated changes imply that subglacial water is moving from one reservoir to another. Active lakes have now been detected in many glacial catchment basins in Antarctica (Smith and others, 2009). Subglacial lakes not only exist under the ice divides, but are also present under some of the upper catchment areas and the fast-flowing trunk portions of some ice streams. Under ice streams, hydrologic systems are an interconnected series of 0.1–1 km³ water bodies that fill and drain on monthly to annual scales. The impact of these subglacial water systems on ice dynamics and ice-sheet mass balance is yet to be determined, although in at least one important case (Byrd Glacier) there has been an observed change in local ice velocity linked to a subglacial flooding event beneath a large outlet glacier (Stearns and others, 2008).

Subglacial water movement beneath large ice sheets is a process about which we know very little, but which we must understand in order to assess its importance for ice-sheet dynamics and ice-sheet mass-balance determination and prediction. With this aim in mind, we need to know where subglacial reservoirs exist, how often water is exchanged between them, what effect water redistribution has on basal shear stress, and how the subglacial conduits are configured beneath the ice streams. Datasets for studying subglacial hydrologic networks are limited, however, due to the inaccessibility of the environment. Therefore, mapping by satellite methods is critical to advancing our understanding in this field.

Mapping subglacial water movement requires high spatial resolution combined with good temporal sampling. However, no single instrument has sufficient spatial and temporal resolution to provide a complete assessment of Antarctica’s active subglacial water system. InSAR data, used by Gray and others (2005), have good spatial coverage where and when they are available, but there are few datasets at high southern latitudes (>82° S). Although RA satellite data have excellent temporal sampling, their coverage extends only to 81.5° S, and the data have severe limitations in higher-relief ice-stream trunk areas because of the width of the RA beam and its interaction with undulating topography. The NASA Ice, Cloud and land Elevation Satellite (ICESat) data extend further south (to 86° S), but in the more northern parts of Antarctica the track spacing is quite sparse, so some active lakes may be missed (Smith and others, 2009). ICESat has provided laser altimeter data for several periods from 2003 to 2009 (Table 1), but at the time of publication (early 2010) laser operations for ICESat have ended, and there will be at least a 5 year gap before another polar-orbiting laser altimeter is launched. Because of the limitations of each dataset, a thorough assessment of Antarctica’s subglacial water system will require a combination of satellite and possibly airborne techniques.

Table 1. Acquisition dates for the 16 ∼33 day ICESat campaigns using the 9 day repeat orbit acquired up to March 2009. Just the last 33 days of the longer laser 2a campaign were repeated in subsequent campaigns. Note that laser 3k and laser 2d are from approximately the same time; laser 3 failed on 19 October 2008 and laser 2 was turned on to complete the set of ground tracks (with some overlap, i.e. tracks 96–145 were repeated). Data from the laser 1 (February–March 2003) and laser 2f campaigns (October 2009) are not analysed here (laser 1 was not in the 91 day repeat orbit, and laser 2f was short and not fully calibrated at the time of writing)

In this paper, we present a multi-technique study of the lake system under the lower MacAyeal Ice Stream (MacIS; formerly Ice Stream E), West Antarctica. The results ofourfirst analysis of MacIS using ICESat alone were published by Fricker and others (2007). Since then, additional ICESat data have been acquired during six further campaigns (making a total of 16 ICESat data-acquisition campaigns; Table 1). In addition, a new, 250 m digital elevation model (DEM) of the region has been generated using ICESat data enhanced with photoclinometry of Moderate Resolution Imaging Spectroradiometer (MODIS) imagery; we have used this DEM in combination with ice-thickness information to produce a hydrologic potential map for the region. We have also generated several MODIS image-differencing scenes for different epochs during the ICESat mission. These ancillary data have provided additional spatial coverage, which we use to better map and understand the active regions detected by ICESat. This is especially relevant for MacIS (∼80° S), where the ICESat track spacing is farther apart than in other regions we have studied (cf. Whillans Ice Stream at 84° S). Its more northerly latitude means that MacIS lies within the limits of European Remote-sensing Satellite (ERS)/Envisat RA coverage, which spans the period 1991 to the present. This opens up the possibility that RA data could be used to extend the time series back to times before ICESat, and forwards to fill in the gap between ICESat and its follow-on mission. By combining these multiple remote-sensing techniques, we have compiled a complete picture of the activity of the MacIS lake system for the period 2003–09. We report on the suitability of each technique for subglacial water mapping, and establish a protocol for studying other systems of active subglacial lakes by using a combination of techniques.

2. Study Region and Glaciological Setting

Our study region is the lower ∼180 km of MacIS (Fig. 1). The MODIS Mosaic of Antarctic (MOA) image over this ice stream (background image in Fig. 1) shows an undulating surface, linear flow stripes, and patches of intense crevas-sing, bounded by two curvilinear shear zones. The undulations are the surface manifestation of rigid bedrock features; the long, oblique feature labeled ‘Mac3’ in Figure 1 is an escarpment in the subglacial geology. In some areas, these bedrock features produce sticky spots, i.e. high local basal shear stresses beneath the ice stream (Stokes and others, 2007). Bindschadler and Scambos (1991) showed that the oblique undulations across MacIS had lower ice-flow velocity than the surrounding ice stream. They proposed that the undulations were the surface expressions of a large subglacial ridge and the associated basal drag. These bed undulations likely are a significant source of basal resistance, because they correspond to where the lubricating sediment has thinned or been totally removed (Stephenson and Bindschadler, 1988).

Fig. 1. (a) ICESat tracks across MacIS showing five regions with significant vertical motion detected by ICESat repeat-track analysis (2003–09), combined with information from a hydrologic potential map and image-difference images. ICESat tracks are numbered and are colour-coded by the total range (i.e. maximum minus minimum) in ICESat elevation from 2003 through 2009. Ice flow is from left to right. Background image is the MODIS MOA (Scambos and others, 2007). White outlines show inferred lake extents; black dotted line for southern limit of Mac3 is due to ambiguous signals there (section 4). The break-in-slope from MOA, which approximates the ice-shelf grounding line (Brunt and others, 2010), is the solid black curve to the lower right of the image. The black dashed squares show the extent of the image-difference images shown in Figure 2. (b) Map of estimated hydrologic potential (kPa, relative to the subglacial pressure at the grounding line) for the same region. The blue dashed curves are the estimated flow paths for water through the system. The yellow ellipses show the approximate starting (large ellipse) and ending (smaller ellipse) locations of the tremors detected by Winberry and others (2009), as derived from their figure 1.

Modeling of MacIS by inversion of the surface velocity field, surface slope and the ice thickness provided estimates of the basal shear stress variations in the trunk area (MacAyeal, 1992; MacAyeal and others, 1995; Joughin and others, 2004). Where vertical shear occurs in fast-flowing regions, it can contribute an order of magnitude more heat than conduction from the crust alone. The englacial heat production, the geothermal heat flux, the advection due to ice flow, and the ice surface temperatures all contribute to the englacial temperature distribution and the basal melt rate. Whereas Mercer, Whillans, Kamb and Bindschadler ice streams (formerly ice streams A, B, C and D respectively) produce the majority of their basal melt in the upstream/tributary sections of their catchments, MacIS produces the majority of its melt in the trunk region, experiencing high basal melt rates throughout much of its downstream area (Joughin and others, 2004). The high melt rate is due to the relatively steep slope and large number of basal shear sticky spots on the MacIS trunk compared to the other ice streams. Joughin and others (2004) showed that around 15% of the MacIS bed is composed of sticky spots, consisting of small (<100 km²) regions where basal shear is up to 50 times greater than in the surrounding areas. These areas and the surrounding areas of weak deformable basal till control the ice flow in regions away from the marginal shear zones. Consequently, the sticky spots convert a substantial quantity of gravitational potential energy into heat energy.

Winberry and others (2009) installed a network of 12 passive seismometers on MacIS in the 2005/06 austral summer. Over the 6 week observation period, they observed two harmonic tremors, each of ∼10 min duration, and interpreted them to be repeated small outburst-flood events,draining from one small lake to another. While the location of these observations does not appear to correlate with any of our identified lakes (the closest was ∼20 km away; Fricker and others, 2007), or match up with periods of drainage from any of them, it does confirm that there is water moving beneath MacIS. As we discuss below, we infer from hydrologic potential mapping that the Winberry and others (2009) events may have been located in an adjacent drainage system within a multi-basin regime beneath MacIS.

3. Methods

The primary technique we used to identify surface deformation indicative of subglacial activity uses satellite laser altimeter data from the Geoscience Laser Altimeter System (GLAS) on ICESat. We used a repeat-track technique in which several repeat elevation profiles along a segment of an ICESat ground track are analysed to locate elevation-change anomalies (Fricker and others, 2007; Fricker and Scambos, 2009). The technique is fundamentally limited, however, by the ICESat inter-track spacing, which is at a minimum near 86° latitude (north or south) and increases with decreasing latitude. MacIS is located at approximately 80° S, where the average ICESat track spacing is 20 km, giving sparser coverage than the more southerly Mercer and Whillans Ice Streams (where the average track spacing is ∼10 km). Therefore alternative techniques are required to fill in the gaps between ICESat tracks.

In this analysis, we combine ICESat data with ancillary remote-sensing datasets (satellite difference images and a hydrologic potential map) to investigate the subglacial hydrology of a system. While ancillary data have previously been used to confirm the validity of elevation-change observations from ICESat (Fricker and others, 2007; Fricker and Scambos, 2009), we show here that they can extend the ICESat data beyond a simple profile, and can also provide contextual information about the subglacial lake. Ancillary data can also be used to remove ‘false positives’ in detecting lakes from ICESat data. Since our study region falls within the coverage of the ERS and Envisat RA satellites, we also examine Envisat RA data over the lakes, to assess the suitability of these data for monitoring lake activity outside the 2003–09 window of the ICESat mission.

3.1. Satellite laser altimetry and repeat-track analysis

ICESat was launched in January 2003 and has since operated primarily in campaigns of ∼33 days duration occurring two or three times per year (see Table 1 for campaign dates). The same ground tracks were repeated during each campaign, so elevation changes can be analysed between repeated tracks. We used an interactive repeat-track analysis method to search for anomalous elevation-change signals in the repeat-track data (Fricker and Padman, 2006; Fricker and others, 2007; Fricker and Scambos, 2009). First, we made a track plot, showing which tracks cross the region. Next, we analysed each ICESat track sequentially to look for elevation-change anomalies. As a first-order cloud filter, we analyse the gain and energy values for all available repeats for each track. This is a simple and approximate method for filtering out cloud-affected surface returns, but we find that it is effective for our purposes (Fricker and Padman, 2006). As an update to earlier techniques (Fricker and Padman, 2006; Fricker and others, 2007), we added a ‘cross-track slope’ module in our processing to correct for the cross-track slope error introduced by the fact that ICESat tracks do not repeat exactly like those employed by Smith and others (2009). We do not, however, apply the cross-track correction here since it assumes (1) a spatially constant cross-track surface slope, and (2) a constant elevation change with time. Neither of these conditions is satisfied, generally, over active subglacial lakes. Furthermore, for elevation-change signals on this order (metre-level) the cross-track slope correction is not as critical as it would be for smaller-scale elevation-change detection.

Using all cloud-free repeats, we interpolated each elevation transect onto a common centre-line track of evenly spaced elevations. Next, we calculated the mean elevation profile of all track data to determine an overall mean elevation profile. We defined the ‘elevation anomaly’ h _i(X) for each repeat as the difference of each elevation profile from the mean profile. This is the same as referencing the profiles to a zero mean. Once a consistent elevation-change event spanning multiple campaigns in a coherent manner was detected along a particular track, we searched intersecting tracks, and then adjacent tracks, to determine the extent of the anomaly. When all events were identified, we calculated the average elevation time series for each lake, using the following procedure: (1) for each track across the lake, we calculated the elevation-change anomalies for each campaign, (2) we identified the limits of the anomaly, i.e. the two points on each track where the anomaly tends to zero (the edges of the lake), (3) we calculated the mean elevation for each campaign between those limits, and (4) we derived the weighted-mean and standard deviation of the ‘on-lake’ elevations for each campaign from all tracks (where the weights were equal to the length of that on-lake track segment divided by the sum of all on-lake track segment lengths present in that campaign). The time we attributed to each campaign epoch was the midpoint of the ∼33 day period. We did not attempt to do any interpolation if campaigns were missing on individual tracks (we just averaged over fewer tracks when fewer are present). This procedure led to a more complete time series than if we had just selected one track to represent the whole lake, since individual tracks have epochs missing due to clouds. The weighted average time series for each inferred lake was then combined with inferred lake areas from satellite images (or satellite image differences) to yield an estimated surface volume-change time series, with a value for each campaign. For each lake, the average elevation change (from all tracks that show a coherent signal, i.e. the elevation anomalies have the same temporal pattern) was calculated and combined with the area estimate for the lake to estimate volume change (from satellite image differencing if available, or from the MOA image if not; Scambos and others, 2007).

3.2. Satellite radar altimetry

MacIS falls within ERS/Envisat coverage (to 81.5° S). ERS RA elevation data at crossovers (points where ascending and descending orbits intersect) were used to infer subglacial lake drainage at Adventure Trench, East Antarctica (Wingham and others, 2006). This was the first time surface deformation related to subglacial lake drainage had been noticed in any satellite altimeter data. We used Envisat RA data (both from crossovers and along-track data) coinciding with the time of the ICESat mission to assess the applicability of this dataset to retrieve elevation time series over the MacIS lakes. If successful, then in principle RA data could be used to extend the history of the lake activity back as far as the start of the ERS 35 day repeat phase (April 1992), which would be helpful in establishing estimates of lake periodicities.

For our crossover analysis, we used elevation time series at Envisat orbital crossovers (Li and Davis, 2006) to examine the variability of elevations on the ice stream for some of the ICESat observation period (January 2003–December 2006) over the elevation-change regions identified by ICESat repeat-track analysis.

For the along-track analysis, we applied a technique similar to that used for ICESat repeat-track analysis. Since the Envisat orbit only repeats the same reference tracks to within +1.5 km, and since it is a 35 day repeat, we divided the ‘swath’ of parallel repeat tracks into thinner ‘strips’, such that all repeats in the same strip were close together. Since there are so many repeat passes, this approach leaves us with sufficient data in each sub-swath strip to examine temporal changes. For the repeat-track analysis, we chose passes that closely correspond in time with the times of the ICESat campaigns. Because the repeats are close together, we also did not account for the slope-induced error since we are considering relative elevations between tracks. For ERS/Envisat altimeter satellite orbits, the precision of the satellite position knowledge is ∼3 cm in the radial (vertical) and ∼10 cm in the lateral (along-track and cross-track) (personal communication from R. Scharoo, 2008). Over ice, the diameter of the RA pulse-limited footprint (PLF) is on the order of 5 km and the diameter of the beam-limited footprint (BLF) is on the order of 20 km.

3.3. Satellite image differencing

Image differencing has been successfully used on other ice streams to confirm the sign of ICESat-derived change and to more accurately delineate the area of the elevation change (Fricker and others, 2007; Fricker and Scambos, 2009; Bindschadler and others, in press). To evaluate the technique on MacIS we began by forming cumulated MODIS images (i.e. an average of several accurately geolocated (50 m accuracy; Wolfe and others, 2002) with similar sun and viewing geometry) from December for each of the years 2003–08. We then subtracted the images in pairs, selecting the pairs to bracket events noted in the ICESat data. The difference images record brightness changes, some of which correspond to changes in slope in regions where the surface elevation is changing (other changes in brightness can represent different frost patterns, residual clouds in the near-cloudless scenes selected, or image noise). We produced image difference series for both the upstream (lakes Mac4 and Mac5; Fig. 1) and downstream regions of MacIS (subglacial lakes Mac1, Mac2 and Mac3; Fig. 1) (Fig. 2 shows an example image difference for each region). We performed a detailed image-differencing analysis for subglacial lake Mac1, the largest lake of this system. For 2007, no clear-sky images at 0800–1000 GMT (the acquisition time selected for the other years) were available for early December, so an additional image set was acquired at ∼1500 GMT for 2006, 2007 and 2008, and a separate series of difference pairs was generated for the later years. The result is a detailed sequence of difference images over the lake, that provides a proxy observation for elevation change over various intervals.

Fig. 2. Difference images for the two regions of MacIS where subglacial lakes were detected by ICESat. (a) The upper region that covers subglacial lakes Mac4 and Mac5 and (b) the lower region that covers subglacial lakes Mac1, Mac2 and Mac3. Illumination in both the difference images is from the upper right. The region south of the outline for Mac1 is a possible elevation-change feature that is not well supported by the ICESat data. Boxes showing the extent of these two regions are shown in Figure 1a.

To determine lake area and extent, we considered where the end points (limits) of the ICESat elevation anomalies lay on the image-difference map. In general, the rules for defining lake boundaries are that (1) the feature(s) must be seen in the multiple individual scenes for a given year (i.e. not clouds or hoar frost, which are ephemeral on a scale of hours to days), (2) the feature(s) must follow some aspect of the non-difference surface topography (i.e. we assume that the lake boundary has some manifestation in the surface morphology), and (3) both the sense of the change and its scale must be consistent with ICESat elevation changes where they overlap. We traced the outline of the feature associated with the edge of the deformation on the difference scene. We estimated the error in this lake area measurement procedure to be 10% (determined through repeated measurements). Finally, we combined the area estimate with the average elevation change during the image-difference interval to estimate volume change.

3.4. Hydrologic potential maps

Water movement at the base of the ice sheet is controlled by the gradient of water-pressure potential, assuming a connected water system. The water-pressure potential is the sum of the gravitational potential at the ice base and the water pressure. The formula for calculating the water-pressure potential is derived from the Bernoulli equation and takes the sum of the bed elevation times the density of water and gravity and the overburden pressure, or the ice thickness times the density of ice times gravity minus the portion of the overburden pressure supported by till and bedrock, also known as the effective pressure (Paterson, 1994). However, it has proven useful to neglect effective pressure and rework the water-pressure potential equation to correspond to what we actually measure: ice surface elevation and ice thickness (Equation (1)). Because of the large effect of ice overburden pressure, the effect of surface topography is approximately nine times that of bed topography and tends to control regional flow, with bedrock topography dominating local flow where bedrock gradients exceed surface gradients by a factor of nine or greater. We calculated the hypothetical hydrologic potential field from maps of surface and bedrock topography and standard assumptions regarding density of the overburden ice and interface water, as follows:

(1)

where P is the water-pressure potential, p _i and p _w are the densities of overburden ice and water, respectively, g is the gravitational acceleration, h is ice surface elevation and H is bed elevation. This assumes that the water pressure is equivalent to the full overburden pressure, and neglects effective pressure (which is on the order of 1–10 m or 10–100 kPa; Engelhardt and others, 1990).

Subglacial water tends to pool in local minima in the water-pressure potential (e.g. Fricker and others, 2007) and flows from one minimum to another along inferred pathways that follow the water potential gradients (e.g. Fricker and Scambos, 2009). This allows us to identify which lakes lie up-/downstream of one another, and provides an independent means of validating the mass transport budget of the subglacial water system (i.e. water volume received by lake A should approximately equal water volume released from lake B, if lake A is the only lake downstream of lake B, unless lake A overflows before lake B finishes discharging). We determined a hydrologic potential gridded at 1 km using the thickness and surface-elevation datasets described below.

3.5. Surface DEM

We constructed a new DEM for MacIS using photoclinometry to enhance detail in an existing 1 km DEM produced from radar and laser altimetry (Bamber and others, 2009). Bindschadler and Vornberger (1994), Scambos and Haran (2002) and Haran and Scambos (2007) discuss development of the photoclinometry method. For this DEM we combined approximately 25 MODIS images acquired during the 2003/04 austral summer to create a quantitative surface slope field at 250 m resolution. The accuracy of the DEM is estimated to be 1–2 m at each 250 m grid node. This captures the ice-stream undulation field (varying at several tens of metres over a few km spatial scale) much better than the more accurate but spatially smoothed original 1 km DEM. We generated the DEM on the World Geodetic System 1984 (WGS84) ellipsoid, and combined the DEM with the 2009 EIGEN-04C geoid (http://icgem.gfz-potsdam.de/ICGEM/ICGEM.html), calculated for the WGS84 ellipsoid, to generate elevations above the geoid.

3.6. Bedrock topography

For the bedrock topography, we used ice thicknesses from radio-echo sounding (RES) surveys where possible, as part of three separate campaigns by the University of Texas Institute for Geophysics (UTIG) Support Office for Aerogeophysical Research (SOAR) from 1995 to 1999 (Blankenship and others, 2001). Additional data were obtained from the Scott Polar Research Institute (SPRI)/US National Science Foundation (NSF)/Technical Univerity of Denmark (TUD) radar sounding campaigns of 1971–75 (Drewry, 1975). Where data were absent, we initially used the BEDMAP-Plus ice-thickness dataset to interpolate between grids (personal communication from J. Johnson, 2009; original BEDMAP dataset is described by Lythe and others, 2001). Although this was effective in the western portion of our survey area, in the eastern regions the misfit between the RES data and BEDMAP-Plus interpolation was >500 m. To obtain ice thickness, we gridded the existing RES data to 5 km; we obtained bed elevation by subtracting the ice thickness from the surface elevation. Where no data were available, we performed a linear interpolation between the nearest measured data.

4. Results and Interpretation

We used the water-pressure potential map (Fig. 1b) combined with image differencing (Fig. 2) to guide our interpretation of the ICESat elevation-change events. The spatial distribution of ICESat-derived elevation-change events is shown in Figure 1a (updated from fig. S2 of Fricker and others, 2007). The distribution shows 14 independent elevation anomalies (i.e. on different tracks) that cluster into five distinct regions. The new water-pressure potential map for this region is shown in Figure 1b, covering the same region shown in Figure 1a. We note that the five regions where the ICESat elevation-change signals cluster correspond closely to local minima in the hydrologic potential.

The 2006–03 difference image for the lower part of MacIS captured elevation changes for the three lower regions (Mac1, Mac2 and Mac3 in Fig. 2b). During this period, the ice surfaces above Mac1 and Mac3 both subsided, and the ice surface above Mac2 uplifted slightly. A detailed imagedifferencing time series for subglacial lake Mac1 is presented below, and we discuss that lake separately. The two subglacial lakes on the upper part of MacIS (Mac4 and Mac5 in Fig. 2a) were captured best in the 2006–05 image. The ice surface above Mac4 drained significantly, while Mac5 filled only slightly during this time period: Mac4 shows a significantly stronger signal in the difference image than Mac5.

We interpret the five elevation-change regions as the surface expression of subglacial lake activity and refer to them hereafter as lakes Mac1–Mac5, numbered sequentially from the grounding line upstream (following the convention of Smith and others, 2009). We note, however, that Mac3 was not included in their analysis, so the numbers are not exactly the same (i.e. their Mac3 is our Mac4). White outlines of lake extent (Figs 1 and 2) are inferred by considering the ICESat signals, the water-pressure potential map and the difference images. We describe each of the lakes in turn, starting downstream and moving upstream.

4.1. Subglacial lake Mac1 (surface elevation 110–130 m)

Lake Mac1, the largest subglacial lake in our study area, produced the largest elevation-change signals (Fig. 3a), which we interpret as corresponding to filling and draining of this lake. Along track 275, the surface elevation decreased by ∼9 m at the point of maximum anomaly (indicated by the central vertical dashed line in Fig. 3a) in early 2004, which we interpret as subsidence of the ice surface as the lake drained substantially. The time series of ICESat elevations averaged over the three tracks across the lake (Fig. 3b) shows the drainage cycle from 2003 through 2009. The lake was ‘full’ at the start of the ICESat mission, then drained, and then filled again in late 2007.

Fig. 3. Subglacial lake Mac1. (a) Repeat ICESat profiles and elevation anomaly for track 275. (b) Time series of elevation/volume changes derived from all three ICESat tracks. Dashed lines indicate missing ICESat campaigns. (c–h) Selection of image-differencing results verifying various stages of the ICESat-derived fill/drain cycle. Red vertical lines in (b) show the approximate timings of the images for each of the years 2003–08.

Drainage of the lake between December 2003 and 2006 is also evident in the image difference for 2003–06, which shows a surface subsidence for this time period. The difference images also provide more spatial information within the lake itself. The drawdown region is characterized by a series of sub-parallel elongate troughs; initially, these appear to refill unevenly, as the 2004–06 image pair shows. Later difference images show a smooth refilling of the depressed area. The area of subsidence from the 2006–03 difference image includes a portion to the south of the main, deep part of the lake. This is a possible elevation-change feature that is not supported by the ICESat data. Two ICESat tracks cross this southern portion (tracks 102 and 275). Track 275 has a small elevation anomaly (∼0.75 m), suggesting subsidence between the two passes closest in time to the image-differencing epochs (laser 2a in November 2003 and laser 3g in November 2006, which are almost exact repeats); however, the laser 3g track is cloud-affected and this elevation anomaly could be entirely due to forward scattering. There are no cloud-free passes for track 102 until laser 3a (November 2004), and we see no evidence of lake activity coinciding with the southern portion of Mac1 on this track. Subglacial lake Mac1 is poorly sampled by ICESat in its southeast section. If there were no other information, this track sampling would lead to significant underestimation of the lake area. We used the image differencing and the hydrologic potential map to obtain a better area for the lake. From the 2006–03 difference image, we estimate the lake area (inside the solid curve) to be 150 km². The MODIS images and the ICESat elevation profiles (Fig. 3a) show the region is in a deep surface low on the ice stream, and the hydrologic potential map shows that the whole lake (including the southern part) sits over a large hydrologic potential minimum (Fig. 1b). Mac1 is just ∼30 km from the grounding line with the Ross Ice Shelf, and the hydrologic potential map indicates that the water from this lake discharges directly into the ocean.

4.2. Subglacial lake Mac2 (surface elevation 187 m)

The elevation-change signal attributed to Mac2 was detected on only two ICESat tracks, 22 and 102 (Fig. 4a for repeat elevation profiles across track 22). The next track to the west (track 394; Fig. 1a for location) showed no evidence of elevation change. The 2006–03 difference image over Mac2 confirms a broad, smooth uplift in this period (cf. the Mac2 region and the uneven and high-amplitude drawdown from Fig. 3c of Mac1, derived from the same images). The difference image also suggests that the area of uplift extends further upstream and downstream than a simple circular assumption from the ICESat anomaly area would imply. We use the difference image to estimate the area of the lake (140 km²). Local water-pressure potential minimum supports the shape and scale of the mapped extent. The average time series derived from these tracks is shown in Figure 4b; it can be seen that many campaigns are missing for track 102. It is possible, however, to discern from this time series that the lake filled between May and November 2005 (laser 3c and 3d campaigns).

Fig. 4. Subglacial lake Mac2. (a) Repeat ICESat profiles and elevation anomaly for track 22. (b) Time series of elevation/volume changes derived from tracks 22 and 102. Dashed lines indicate missing ICESat campaigns. Red vertical lines show the approximate timings of the images comprising the 2006–03 image difference.

4.3. Subglacial lake Mac3 (surface elevation 200 m)

Subglacial lake Mac3 is located under the oblique, stripe-like feature that crosses MacIS (visible in MOA; Fig. 1; discussed in section 3). We noticed that the time series for the six ICESat tracks that were thought to cross lake Mac3 were not all in phase with each other (cf. Fig. 5b and Fig. 5c and d). The four downstream signals on this feature (tracks 1322, 141, 340 and 221 in the northwest) were distinctly different from the two upstream signals (tracks 87 and 260 in the southeast). In addition, the two upstream signals were associated with ice surfaces that were ∼60 m higher than the average elevation of the four downstream signals, suggesting that these signals did not all correspond to the same lake.

Fig. 5. Subglacial lake Mac3. (a) Difference image for 2008–04. White outline shows inferred subglacial lake extents. Black dotted line for the southern limit is due to ambiguous signals there. (b) ICESat elevation profiles and elevation anomalies for track 141, interpreted as a true lake signal. (c, d) ICESat elevation profiles, elevation anomalies and differences from first repeat (laser 2a) for tracks 260 and 87 across the downstream region of the diagonal sticky spot shown in Figure 1. (e) Averaged time series derived from four ICESat tracks. Dashed lines indicate missing ICESat campaigns. Red vertical lines show the approximate timings of the images comprising the 2006–03 and 2008–04 image differences.

The 2006–03 difference image shows a significant drawdown/subsidence signal that is confined to the northwest (downstream) part of the stripe feature associated with subglacial lake Mac3 (Fig. 2). However, from the ICESat time series (Fig. 5e) it appears there is little or no activity between the ICESat campaigns closest to the image-acquisition dates. We believe this discrepancy is due to the relative timings of the images versus the ICESat campaigns, i.e. some activity was detected by the difference images (with image acquisitions in early to mid-December) that was not captured in the ICESat data (campaigns in October–November). The 2008–04 difference image shows a lower-magnitude signal over a wider area; this is the area we use to estimate the lake volume (142.6 km²). The southern portion of this lake is delineated by a dotted line, as explained below. The hydrologic potential map (Fig. 1b) indicates that the local hydrologic potential minimum coincides with approximately the same area as is outlined in the 2008–04 difference image. This is a powerful example of how different datasets can be used to guide ICESat data interpretation.

A time series constructed from the four lake signals falling within the lake outline defined by the 2008–04 image difference and the hydrologic potential map is shown in Figure 5e. This shows that the lake drained in early 2004 and again in late 2007/early 2008. This lake appears to have a ∼3 year cycle, although this is based on limited data.

We left the southern end of lake Mac3 ‘open’ as a dashed black line in the outline (e.g. Fig. 1a), because of the ambiguous nature of the signals from our various datasets in this complex region. In the southern part, the stripe-like feature is complicated by steep topography and intense flow dynamics, and has been identified as an unusually intense basal shear feature for some time (Bindschadler and Scambos, 1991). In this region, the vertical movements and the difference image features are not easily explained by simple filling and draining of water beneath the ice sheet. We believe that other events (e.g. sediment movement, erosion, or ice thinning and thickening associated with ungrounding and regrounding) are also occurring. Sergienko and others (2007) propose and model one plausible mechanism (thinning at the upstream end of an ungrounding event; thickening at the downstream end) that may partially explain what is happening in this area. Figure 5c and d show the two ICESat signals across this part of the feature (tracks 260 and 87 respectively). For these tracks, we have plotted the elevation-change anomaly and also the difference in elevation from the first repeat (laser 2a), for direct comparison with figure 4c of Sergienko and others (2007). For track 260 (Fig. 5c) we note that there is thinning upstream of the event and thickening downstream, and that the elevation-change signal has a dipole-like structure. This may be an example of the effect identified by Sergienko and others (2007). At this location, the ice crosses a sticky spot and experiences an increase in basal resistance. It is possible that some of the resulting elevation signal is a response to this effect, and not due to subglacial water movement. Analysis of several difference-image pairs reveals a series of subtle, narrow (∼750 m) ridges moving from upstream to downstream at rates similar to the ice-flow speed (∼350 m a⁻¹). One hypothesis for these features is berms of sediment being pushed off the sticky spot and carried downstream beneath the ice stream as the ice flows over the feature; another is variations in ice thickness. It is possible that there is some drainage and refilling in this region as well, but the ICESat signal may be modified by these additional processes. An alternative hypothesis is that these are wind depositional features, since the steep topography will affect the local winds. Since we do not have enough information at this time to test or characterize these hypotheses, we have dotted the southern end of the outline for Mac3, as we cannot exclude the possibility that there is subglacial water activity occurring beyond there. This will be the subject of future study.

4.4. Subglacial lakes Mac4 and Mac5 (surface elevations 375 m and 409 m)

The most upstream lakes on MacIS trunk (lakes Mac4 and Mac5) lie close together and at similar surface elevation. ICESat track 364 captures both lakes (Fig. 6a), and lake Mac5 is also crossed by track 325. Comparing the locations of the ICESat elevation-change signals with the hydrologic potential plot (Fig. 1b) shows that they both correspond to local minima in the estimated hydrologic potential. The areas of lakes Mac4 and Mac5 digitized from a 2006–05 difference image are 57 and 72 km² respectively. The area of deformation for Mac4 was not obvious in the difference image, but we extended it based on the ICESat elevation anomalies for track 325. The volume time series for these lakes from tracks 325 and 364 are shown in Figure 6b. These time series suggest that lakes Mac4 and Mac5 sometimes act as one lake and sometimes as individual lakes. There are some instances where draining of Mac5 appears to lead to filling in Mac4.

Fig. 6. Subglacial lakes Mac4 and Mac5. (a) ICESat elevation profiles and elevation anomalies for track 364 across both lakes. (b) Averaged volume time series derived from ICESat tracks 325 and 364. Dashed lines indicate missing ICESat campaigns (note there are no data for campaigns later than laser 3g (November 2006)). Red vertical lines show the approximate timings of the images comprising the 2006–05 image difference shown in Figure 2a.

4.5. Envisat RA-2 radar altimetry

We examined data from Envisat’s RA-2 instrument (Roca and others, 2009) to see whether it could provide additional information on the history of lake activity. Figure 7a shows the total range in elevation between 2003 and 2006 at each Envisat crossover, with our derived active lake locations overlaid. Note that, in general, crossovers with the largest range in elevation align with topographic features evident in the background MOA image. There are two 35 day repeat Envisat crossovers that show a large (>5 m) elevation range on lake Mac1 (labelled X1 and X2 in Fig. 7a). Unfortunately, there are no crossovers from the central part of the lake. Initial examination of the individual crossover time series at X1 and X2 suggests some evidence of a cyclical pattern. When we compared these with elevations from ICESat over the same lake and time period, however, we found no correlation between the two types of altimetric time series (Fig. 7b).

Fig. 7. Envisat RA data analysis over the MacIS subglacial lakes. (a) Locations of Envisat orbital crossovers, colour-coded by the total range in elevation, 2003–07; blue and red circles on track Env-197 illustrate the approximate diameters of the RA BLF and PLF respectively. (b) Elevation time series for two crossovers on subglacial lake Mac1. (c) Near repeats of Envisat track 197 across subglacial lake Mac1, and elevations derived from the photoclinometric/ICESat DEM.

We also analysed along-track Envisat RA-2 data for relative orbit (the European Space Agency’s terminology for ‘track’) 197 across lake Mac1 (Fig. 7a for track location). We interpolated the image-enhanced DEM of the region along this RA track, and compared the interpolated elevations with the Envisat elevations, and found significant topographic information to be missing from the RA-derived elevations (Fig. 7c).

The lack of correlation between the RA time series and the ICESat time series is the result of how the RA instrument operates over undulating topography (with undulations on the scale of the BLF, which is ∼20 km in diameter). Lake Mac1 lies within a topographic dip of 60 m relief, and is ∼10–15 km in length along the RA tracks. Over surface topographic features with length scales of the same order as the RA BLF, the relationship between the altimeter waveform and the surface elevation is not unique (Wingham and others, 1993). Topography on this scale interacts in a complex way with the RA waveform, which affects the elevation-change measurement. Indeed, Figure 7c shows that the RA profiles do not follow the surface topography of Mac1. These elevation profiles have not been corrected for the slope-induced error, but it is clear that the problem of missing topographic information would not be reconciled by performing this correction.

This example demonstrates that care must be taken when using RA data to look for elevation changes due to subglacial water movement, and suggests that the surface topography of the region needs to be assessed prior to undertaking such a study. In particular, it shows that RA data from ERS and Envisat are unlikely to be useful for studying active lakes in regions of rapid streaming ice flow where undulation topography varies by tens of metres over a few km. RA data can be useful for studying large lakes (Ridley and others, 1993) and for smoother, high-plateau regions of the ice sheet (e.g. Adventure Trench, in East Antarctica (Wingham and others, 2006)), but the surrounding topography must be smooth, on the length scale of the BLF. In these regions of the ice sheets, large changes in elevation from other sources are unlikely, so the active subglacial lake signals are unambiguous. Examples of other regions where RA is likely not appropriate include the other Siple Coast ice streams, the Ronne–Filchner ice streams, and the ice streams flowing into the Amery Ice Shelf. The next RA satellite (CryoSat-2, to be launched February 2010) is expected to perform better than RA-2 over rough surfaces, and therefore may provide useful data over ice-stream lakes.

5. Discussion

5.1. Flow of water through the MacIS subglacial lakes system

We estimated flow paths from the hydrologic potential map by following minima in the hydrologic potential (blue dashed curves in Fig. 1b). For increased accuracy we also identified the local minima in the hydrologic potential along the RES line data that lie closest to the intersection of our preliminary flow path and the flight path of the RES data; this allows us to verify that the hydrologic potential is monotonically decreasing along the flow path. When these flow paths are considered alongside the average volume time series 2003–09 for the MacIS lakes (Fig. 8) we can begin to make some conclusions about the hydrology of the system and how the lakes are connected. The flow paths suggest that, in general, flow of water through the system is from lake Mac5 to lake Mac4 to lake Mac3 to lake Mac1. Lake Mac2 appears generally separate, although it may have some connection to Mac1. The time series show that when lake Mac5 drained in late 2005/early 2006, both lakes Mac3 and Mac4 filled. Similarly, lake Mac1 filled between early and late 2007 as lake Mac3 drained. In each case, the volumes of water are approximately consistent (within the error bars). There was another drainage event for Mac3 in early 2004, but at that time Mac1 was also draining, so it is not possible to tell whether the floodwater from Mac3 was received by Mac1.

Fig. 8. Average volume time series for the five subglacial lakes in the MacIS system.

The southern hydrological flow path on MacIS appears to coincide with the locations of the seismic events observed in 2005/06 by Winberry and others (2009). This suggests that the Winberry and others (2009) lakes are relatively small, so we propose that they probably drain frequently. These lakes are hydrologically separate from the ICESat-detected lakes in the MacIS system.

5.2. Location of subglacial lakes with respect to basal shear stress

We have updated Joughin and others’ (2004) determination of basal shear stress using the new DEM, and the result of this inversion is shown in Figure 9, overlaid with our derived lake outlines. Three of the MacIS lakes (Mac1, Mac2 and Mac3) coincide with locations just downstream of local maxima in basal shear stress, whereas lakes Mac4 and Mac5 do not. High basal shear stress corresponds to sites of melting, as well as regions that lead to locally thinner ice, which lowers the basal water-pressure potential. Lake Mac3 is clearly related to high basal shear stress. The sticky spot (possibly due to lack of till, and rock protruding into the ice) may be the source of a large volume of water that influences the localized ice flow by virtue of the basal friction changes. As the ice flows over the sticky spot, high-energy dissipation and melting occurs. The subglacial topography is such that the water generated is funnelled towards the base of the subglacial escarpment, which essentially coincides with Mac3. In the ‘wake’ of the subglacial escarpment, subglacial lithostatic pressure may be lowered, because the ice must deform over a short downstream distance (a few ice thicknesses) to re-engage the bed. If subglacial water is produced by strain heating at the sticky spot, it is likely to accumulate in the subglacial pressure low downstream of the spot (Stokes and others, 2007). Note that this dynamically induced low pressure is a superimposed pressure reduction relative to the hydrologic potential map shown in Figure 1. Lakes Mac1 and Mac2 also seem to coexist with sticky spots, but perhaps not in such an obvious manner.

Fig. 9. Basal shear stress derived for MacIS using the method of Joughin and others (2004) with image-enhanced DEM; outlines of MacIS subglacial lakes are overlaid.

The basal traction contrast at a sticky spot causes the ice to thin downstream, reducing the pressure at the ice base. The reduction in overburden pressure downstream of a sticky spot creates a local minimum in the hydrologic potential in which water collects. The additional water in this depression further enhances the traction contrast. Downstream of the hydrologic potential depression, an additional traction contrast occurs where the ice regrounds. When ice regrounds on the downstream side, the corresponding increase in basal stress causes the ice to thicken, ultimately deepening the well in the hydrologic potential and often creating a seal (Fowler, 1999). Modelling and observational work indicate that sub-ice-sheet floods drain only the top few metres of a lake when this seal does eventually break (Evatt and others 2006; Carter and others, 2009), and thus at least a portion of the lake is expected to persist between floods, especially when at least a portion of the meltwater is sourced locally. Therefore a lake forming downstream of a sticky spot will not only have a constant source of water but also create glaciological conditions which support its continued existence. It therefore appears that some of the MacIS subglacial lakes may be a product of local melting, i.e. that there will be a ‘co-evolution’ of lake and sticky spot. We expect the lakes and the basal resistance in between the lakes to co-evolve and that this co-evolution will produce important influences on the velocity of ice flow. A more detailed exploration of this phenomenon is, however, beyond the scope of this paper.

6. Conclusions

We have examined the subglacial water system under MacIS using multiple techniques (ICESat laser altimetry, MODIS image differencing and a hydrologic potential map). Analysis of the ICESat-derived signals along with the image differencing and hydrologic potential map has allowed us to determine which signals are likely to be associated with lake activity, and which are due to the surface response of the ice to other subglacial conditions. MODIS image differencing over the lakes confirms the locations of the ICESat-detected events and the sign of the ICESat signals, and also provides more spatial detail for mapping the lake. This technique has its limitations, however, and requires a signal of around 2 m of net elevation change over a region of several km² for it to be discernible above the typical levels of noise and spurious features. The hydrologic potential map confirms that the lake positions and derived areas correspond to local hydrologic potential minima. Our study demonstrates how several datasets can be merged together to study a subglacial system in a holistic approach, to learn more information about the system than using any technique alone. These datasets complement each other and, when used interactively, can extend the ability of each to detect, map and track the evolution of an active subglacial lake.

Our analysis detected five active lakes in the MacIS system. We estimated the hydrological flow paths of water through the system, which suggest some connectivity between some of the lakes. This is generally consistent with the temporal patterns in the ICESat observations of surface elevation change over the lakes. We also reassessed the basal shear stress for the ice stream using the new DEM, applying the technique of Joughin and others (2004). Three of the active lakes appear to coexist with maxima in the basal shear stress, or sticky spots. This correlation suggests that not all of the water in the lakes is from the wide hydrological system, but that a significant amount is actually produced locally through basal melting.

Since ICESat-1 has reached the end of laser operation, and its follow-on (ICESat-2) will not be launched until at least 2015, there will be no more satellite laser altimetry data collected over Antarctica in the next 5 years. Some airborne laser altimeter data will be acquired by NASA’s ‘IceBridge’ mission, which is designed to fill the gap between ICESat-1 and its follow-on. We anticipate that such targeted airborne acquisitions over this and other active lake regions will contribute to continued monitoring of their activity. Several lakes (four on Whillans and Mercer Ice Streams and one on Recovery Ice Stream) have also been instrumented with continuous GPS that should operate throughout most of this gap. However, to obtain broader spatial coverage we need to consider how far our other potential methods can continue to monitor subglacial activity. Our study suggests that image differencing will be able to confirm whether lakes have drained or filled, provided the net elevation-change signal is sufficiently large, but can provide only approximate information on the magnitude of the signal, i.e. they provide qualitative and not quantitative information. Furthermore, image differencing can be ambiguous unless the elevation-change signal is large. Hydrologic potential mapping can help identify the most likely sites for subglacial water accumulations, but provides no information on activity. We have also demonstrated that the performance of the satellite RA is challenged over the ice streams, so that RA data are limited for identifying active lakes or extending the current ICESat time series backwards and forwards in time. It is anticipated that CryoSat-2 (to be launched in February 2010) will perform better than RA-2 over rough surfaces, but until we receive data the exact performance is unknown. We believe that the laser altimeter is currently the most suitable instrument for routinely monitoring lake activity on ice streams.