1. Introduction
There are 379 known subglacial lakes under the Antarctic ice sheet, which have been inferred through a variety of remote-sensing techniques (Reference Wright and SiegertWright and Siegert, 2012). Before 2005, airborne radar sounding (e.g. Reference Oswald and RobinOswald and Robin, 1973; Reference Siegert, Carter, Tabacco, Popov and BlankenshipSiegert and others, 2005; Reference Carter, Blankenship, Peters, Young, Holt and MorseCarter and others, 2007; Reference Wright and SiegertWright and Siegert, 2012) was the primary tool used to identify subglacial lakes. Recently, satellite measurements of surface elevation change (e.g. Reference Gray, Joughin, Tulaczyk, Spikes, Bindschadler and JezekGray and others, 2005; Reference Wingham, Siegert, Shepherd and MuirWingham and others, 2006; Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference Smith, Fricker, Joughin and TulaczykSmith and others, 2009) and satellite image differencing (e.g. Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference Fricker and ScambosFricker and Scambos, 2009) have been used to show where the surface deforms in response to subglacial water movement. This technique has been used to identify 130 ‘active’ subglacial lakes. Analyses of repeat-track satellite laser altimetry from NASA’s Ice Cloud and land Elevation Satellite (ICESat) have shown that subglacial lakes underlie many of Antarctica’s ice streams and outlet glaciers, undergoing volume changes of several km3 on monthly to decadal timescales (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others 2007; Reference Fricker and ScambosFricker and Scambos, 2009; Reference Smith, Fricker, Joughin and TulaczykSmith and others, 2009). In one location (Byrd Glacier), a lake drainage event was concurrent with temporary acceleration of ice downstream (Reference Stearns, Smith and HamiltonStearns and others, 2008), implying that lakes may influence ice dynamics. More recently, studies have shown that where subglacial lakes exist, they capture nearly all the meltwater from the catchment above them (Reference CarterCarter and others, 2011; Reference Carter and FrickerCarter and Fricker, 2012), making it possible to infer water routing through a system from simultaneous observations of volume change for all lakes in the system.
Progress in Antarctic subglacial lake research over the past two decades has led to three separate, large international projects that directly access and explore a subglacial lake. One such project is the Whillans Ice Stream Subglacial Access Research Drilling (WISSARD) to explore Subglacial Lake Whillans (SLW) (Reference Fricker, Siegert, Kennicutt and BindschadlerFricker and others, 2011; Reference PriscuPriscu and others, 2013). SLW is located under a low-relief yet rapidly sliding region called the Whillans Ice Plain (WIP) (Fig. 1). WIP is the site of a substantial number of observations over the last 50 years (Reference Bindschadler, Vornberger and GrayBindschadler and others, 2005) and likely has undergone much larger-scale flow reorganization in the region over the last millennium (Reference Catania, Hulbe, Conway, Scambos and RaymondCatania and others, 2012). To complement the WISSARD fieldwork, we investigate the origin of water entering SLW, explore and characterize the hydrologic connections throughout the surrounding region, and evaluate how changes in the hydrology relate to changes in the flow of the overlying ice. We examine the entire WIP using repeat-track ICESat altimetry data (2003–09) with a water budget model to obtain a detailed picture of the water system to see how it responds to changes in the ice geometry and the filling and drainage of lakes within the system and further upstream over the ICESat observation period.
Map model domain with labels of primary areas of interest. Abbreviations: BIS: Bindschadler Ice Stream; KIS: Kamb Ice Stream; WIS: Whillans Ice Stream; WIP: Whillans Ice Plain; SLE: Subglacial Lake Engelhardt; SLW: Subglacial Lake Whillans; L78: Lake 78; SLM: Subglacial Lake Mercer; SLC: Subglacial Lake Conway; USLC: Upper Subglacial Lake Conway; W6: Lake Whillans 6; W7: Lake Whillans 7; W8: Lake Whillans 8; K1: Lake Kamb 1; K5: Lake Kamb 5. All maps in this paper use the Polar Stereographic projected with a standard latitude of 71S. Shading denotes ice surface velocities from Reference Rignot, Mouginot and ScheuchlRignot and others (2011).
2. Known Subglacial Hydrology of Whillans Ice Plain
WIP is at the confluence of Whillans Ice Stream (WIS) and Mercer Ice Stream (MIS). Here ice surface velocities range from 100 to 600 m a–1 (Reference Rignot, Mouginot and ScheuchlRignot and others, 2011) despite low surface and basal slopes (generally <0.001), low driving stresses (Reference Bindschadler, Vornberger and GrayBindschadler and others, 2005), and widespread freezing conditions at the ice base with accretion rates approaching 5–6 mm a–1 (Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004) (Fig. 2). Although small quantities of meltwater are generated by frictional heating over regions of enhanced basal traction or ‘sticky spots’ (Reference Stokes, Clark, Lian and TulaczykStokes and others, 2007; Reference Sergienko and HulbeSergienko and Hulbe, 2011), most of the basal lubrication responsible for sustaining such high flow rates comes from basal meltwater generated upstream (Reference Bougamont, Tulaczyk and JoughinBougamont and others, 2003; Reference Parizek, Alley and HulbeParizek and others, 2003). This volume of meltwater delivered to WIP may be enhanced by possible water piracy from neighboring Kamb Ice Stream (KIS) following its shutdown ∼165 years ago (Reference Anandakrishnan and AlleyAnandakrishnan and Alley, 1997). It is uncertain, however, whether this additional water offsets any of the ongoing basal freezing in lower WIS. Changing hydrological conditions have been implicated in the slowdown of WIP glaciers: lower WIS and MIS ice flow velocities decreased on average by 0.6% yearly since 1985 (Reference JoughinJoughin and others, 2005; Reference Scheuchl, Mouginot and RignotScheuchl and others, 2012; L.H. Beem and others, unpublished information), with greater deceleration downstream. This deceleration has led to a regionally variable dynamic thickening of the ice of up to 0.5 m a–1 (Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009).
Map of model domain of Kamb, Mercer and Whillans ice streams showing data used for this study. Background shading is from the MODIS Mosaic of Antarctica (MOA) image of Reference Haran, Bohlander , Scambos and FahnestockHaran and others (2005). Boxes are shown for the areas covered by Figures 3, 7 and 10. Abbreviations: SPRI: Scott Polar Research Institute, UK; RIGGS: Ross Ice Shelf Geophysical and Glaciological Survey; SOAR: US Support Office for Aerogeophysical Research; ASAID: Antarctic Surface Accumulation and Ice Discharge.
Movement of water through a subglacial system is governed by the hydropotential (θh ) at the base of the ice. In Antarctica, the hydropotential is typically calculated as the sum of the bed elevation and water pressure (assuming water pressure is equal to the overburden pressure) using models for the surface elevation and the ice thickness (Reference ShreveShreve, 1972). In this framework (Appendix), surface elevation gradients are 11 times more significant than bed gradients in determining the direction of water flow (Reference ShreveShreve, 1972; Reference Iken and BindschadlerIken and Bindschadler, 1986). Bedrock gradients can, however, locally exceed surface gradients by more than a factor of 11.
Earlier investigations into the subglacial hydrology of WIP suggest that three distinct hydrologic basins exist, defined by the regional hydropotential derived from a high-resolution Moderate Resolution Imaging Spectroradiometer (MODIS)- based photoclinometric digital elevation model (DEM) for the surface, an older database for ice thickness (Reference Shabtaie and BentleyShabtaie and Bentley, 1987; Reference Fricker and ScambosFricker and Scambos, 2009), and connected lake events inferred from ICESat data (Figs 2 and 3) (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference Fricker and ScambosFricker and Scambos, 2009). The three inferred hydrologic basins, named after their relative geographic location, are: (1) the Southern Basin following the southern margin of WIS trunk and WIP to the lower trunk of MIS, containing Upper Subglacial Lake Conway (USLC), Subglacial Lake Conway (SLC), Subglacial Lake Mercer (SLM) and ‘Lake 7 and Lake 8’ (originally thought to be two separate lakes but now considered a single lake, L78); (2) the Central Basin containing Subglacial Lake Whillans (SLW); and (3) the Northern Basin following the northern margin of WIS trunk and WIP, containing Subglacial Lake Engelhardt (SLE) (Fig. 1).
Steady-state water flux over model domain (see box in Fig. 2 for location) from ‘control’ model with contours of hydropotential (black), and difference between tuned hydropotential and original hydropotential (pink and green).
Reference Carter and FrickerCarter and Fricker (2012) identified several possible cascading lake drainage events that originated in the upper reaches of KIS and WIS and culminated in the drainage of lakes in WIP, indicating that all three basins are supplied with water via the WIS trunk. The point where the water divided downstream into three flow paths (hereafter referred to as the ‘flow trifurcation region’) is also a region that has no RES data coverage (Fig. 2). However, on a larger regional scale, the new Bedmap2 dataset shows that basal relief of the area is ∼200 m over tens of km (Reference FretwellFretwell and others, 2013), suggesting that surface topography, with undulations roughly 20 m or more over the same length scales, is the primary driver of water routing (Eqn (1)). We model the linkages between the lakes in upper KIS, northern and southern WIS and WIP to the subglacial hydrology of the region in greater detail, using the full ICESat observation period (2003–09), exploring how changes in the flow trifurcation region affected water distribution among the three drainage basins during this observation period. Our study focuses on the largest lakes in the system, primarily, K1, K5, W8, W7, W6, USLC, SLC, SLM, L78, SLW and SLE. Although Reference Fricker and ScambosFricker and others (2009) and Reference Smith, Fricker, Joughin and TulaczykSmith and others (2009) reported additional lakes in our study area, work by Reference Carter and FrickerCarter and Fricker (2012) suggested that the volume change undergone by these lakes is too small to play a significant role in the water budget, and they are not included in this study.
3. ICEsat Analysis of Lake Volume Time Series and Local Elevation Change
3.1. ICESat data and analysis
ICESat operated two or three times yearly between 2003 and 2009, primarily in 33 day ‘campaigns’. Each campaign used the same reference ground tracks (within 111 m RMS deviations; Reference Siegfried, Hawley and BurkhartSiegfried and others, 2011), allowing detection of along-track elevation changes. In earlier work, we identified the lakes using an interactive repeat-track analysis method to search for anomalous elevation-change signals in the repeat-track data (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference Fricker and ScambosFricker and Scambos, 2009). Previous studies on lake volume change (e.g. Reference Fricker and ScambosFricker and Scambos, 2009; Reference Smith, Fricker, Joughin and TulaczykSmith and others, 2009) and regional surface elevation change (e.g. Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012) were based on early releases of ICESat and did not include data from later campaigns (2008–09). In this study we extend the time series over all lakes and use a later data release (GLA12 product, release 633).
We used ICESat data for two applications, which required two different processing methods: (1) estimation of lake volume change (Section 3.1.1), and (2) estimation of changes in the surface elevation (Section 3.1.2). For (1), we limited our calculations of surface elevation anomalies to previously determined lake outlines from Reference Fricker and ScambosFricker and Scambos (2009), while for (2) we constrained our analysis to the trifurcation region where our understanding of the hydropotential surface from previous work (Reference Carter and FrickerCarter and Fricker, 2012) indicates the presence of a critical point for allocating water among the three basins of WIP.
To preprocess the ICESat data for both (1) and (2), we first used the instrument gain as a first-order cloud filter. This is a simple and approximate method for filtering out cloud-affected surface returns, but we found it to be adequate for our purposes (e.g. Reference Fricker and PadmanFricker and Padman, 2006). We also applied a correction for inter-campaign biases (personal communication from T. Urban, 2013, updated from Reference Urban and SchutzUrban and Schutz, 2005; Reference Siegfried, Hawley and BurkhartSiegfried and others, 2011). We note that since the release of R633, a data-processing issue has been identified that may lead to errors of up to 10 cm (Reference Borsa, Moholdt, Fricker and BruntBorsa and others, 2013). We do not apply this shot-by-shot correction, but rather apply the inter-campaign bias correction, which is constant over each campaign, as the differences between the corrections are an order of magnitude below the signal we observe. Future studies will examine the impact of the new correction on our findings as more is understood about this issue.
Additionally for (2), we applied a cross-track slope correction between non-repeating ground tracks. This correction assumes: (a) a stable cross-track surface slope on scales of 0.5 km or less, and (b) a constant elevation change with time – conditions that generally appear valid for the trifurcation region (Section 3.2; Fig. 4). Over active subglacial lakes, however, neither of these conditions is generally met, so we did not apply the correction for application (1). The expected impact of the cross-track slope is believed to be small relative to the total elevation change over an active subglacial lake and consequently less important than it would be for detecting regional surface change (Reference Moholdt, Hagen, Eiken and SchulerMoholdt and others, 2010).
Elevation anomaly from ICESat over time for selected gridcells in the flow trifurcation region.
After these preprocessing steps, for both (1) and (2) we calculated the surface change relative to the L2a campaign (November 2003) using repeat track analysis, as detailed in the following subsections.
3.1.1. Lake volume estimates
We estimated the volume change for all lakes in our study area using the following approach:
-
1. Identify all ICESat tracks that cross a given lake, using the lake outlines from Reference Fricker and ScambosFricker and Scambos (2009) and Reference Smith, Fricker, Joughin and TulaczykSmith and others (2009). Collect all available repeats along those tracks and identify the limits of the along-track elevation change anomalies as described by Reference Fricker, Scambos, Bindschadler and PadmanFricker and others (2007).
-
2. Within these limits, calculate the average elevation for each repeat (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007), applying the corrections identified in Section 3.1 and an additional correction for regional ice thickness change (Reference Smith, Fricker, Joughin and TulaczykSmith and others, 2009), producing a precise time series of average elevation for each track across the lake. Not all ICESat tracks contain valid data during every campaign; some tracks may be missing due to instrumental or environmental (e.g. clouds) issues. For tracks with missing data for a campaign, we interpolate the average elevation using simple linear interpolation (Fig. 5), using other tracks on the same lake that had data for that campaign.
As an improvement over previous studies, we include ‘off-pointed’ repeat tracks in the volume-change time series. During the ICESat mission, some repeats were purposefully pointed off-track for certain campaigns to sample specific locations known as ‘targets of opportunity’ (TOOs), resulting in a new set of repeat tracks parallel to the original track but offset by 0.5–5 km. In previous studies, these off-pointed repeats resulted in a temporal data gap for that campaign; in our study, when the off-pointed repeats land on a lake, we incorporate them into the volume-change time series, simultaneously improving both the spatial and temporal sampling.
-
3. The final average elevation change for each campaign is a length-weighted average of the elevation anomaly for each repeat:
(1)
where ΔV L is the lake volume change, n is the number of tracks crossing the lake, Δhk is the average elevation anomaly across the lake for track k, lk is the length of the track k crossing the lake and A L is the area of the lake. This formula assumes that lake area remains constant with lake level. When compared with the results of Reference Fricker, Scambos, Bindschadler and PadmanFricker and others (2007), Reference Fricker and ScambosFricker and Scambos (2009) or Reference Smith, Fricker, Joughin and TulaczykSmith and others (2009), whose algorithms did not include any interpolation or extrapolation to account for missing campaigns nor include off-pointed tracks, the degree of monotonicity for the lake volume time series increased substantially (Fig. 6).
Time series of estimated surface elevation changes for three major subglacial lakes in the system: (a, b) SLE; (c, d) SLM; and (e, f) K5. Colored lines show the mean elevation estimates from the individual ICESat tracks before (a, c, e) and after (b, d, f) interpolation for missing campaigns. The thicker black line is the weighted mean of all tracks. For K5, tracks 220 and 101 were off-pointed for several campaigns to survey a target of opportunity (TOO) and were therefore considered separate tracks (220b and 101b). Note also the elimination of outlier data points for tracks 1331 and 220.
Volume time series for all lakes in this study. Previous time series (Reference Smith, Fricker, Joughin and TulaczykSmith and others, 2009, for K1, K5, W6, W7 and W8; Reference Fricker and ScambosFricker and Scambos, 2009, for USLC, SLC, SLM, L78, SLW and SLE) are shown with a dashed line. Dotted lines denote linear extrapolation. Panels are divided among the various basins in our study area with (see Fig. 1) (a) upper WIS, (b) upper KIS, (c) Central Basin, (d) Southern Basin and (e) Northern Basin.
Our lake volume change calculation has two primary sources of error: errors in the estimate for the area of the lake and errors in the elevation differences (Reference Fricker and ScambosFricker and Scambos, 2009). For lakes K1, K5, W6, W7 and W8, for which no MODIS analysis was performed, the net sum of these errors was ±50% (Reference Smith, Fricker, Joughin and TulaczykSmith and others, 2009) while the errors for the remaining lakes were closer to ±20% as their shorelines were better constrained (Reference Fricker and ScambosFricker and Scambos, 2009).
3.1.2. Surface elevation change over flow trifurcation region
We used ICESat data to estimate the surface elevation change in the flow trifurcation region. We divided the region into 2.5 km × 2.5 km gridcells (the same resolution as the hydropotential grid (Appendix), but less than ICESat track spacing in this region). To accommodate the cells not crossed by ICESat tracks, the final model contained some degree of spatial interpolation, described later in this section. For each gridcell crossed by an ICESat track we: (1) average the measurements for each campaign relative to the L2a; (2) apply a 1 year moving-average filter to each bin’s time series, filling in for missing campaigns, producing a continuous time series; and (3) for each campaign, perform an inverse-distance weighted interpolation of the elevation anomaly to fill in the empty gridcells, obtaining a continuous grid of elevation anomalies across the trifurcation region (Fig. 7). After campaign Laser 3k any surface change in this area is calculated through 2008 and assumed to have remained constant thereafter.
Surface elevation anomaly and water flux across the flow trifurcation region of lower WIS trunk (see box in Fig. 2 for location), for selected time-steps of the ‘main model run’.
In the process of converting elevation change to fluctuations in hydropotential, we assume all elevation change results from ice thickness change rather than bed elevation change or changes in the bulk density of the firn layer. Although Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others (2011) suggest that inter-annual variability in snowfall may lead to apparent surface elevation changes on the order of tens of cm, such changes are expected to be sufficiently spatially homogeneous across our relatively small study area. We also note that the modeled elevation change rate due to firn effects over the adjacent Ross Ice Shelf was on the order of a few cm a–1 (Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012). From these assumptions the relation between surface elevation change and hydropotential change is
The importance of surface change over this region and its connection to lake activity are described in more detail in our model results and discussion.
3.2. ICESat results for lake time series and elevation-change detection
3.2.1. Lake volume time series
Overall the shape and magnitude of the volume time series for all lakes in the lower WIS/MIS changed only slightly from previous work of Reference Fricker and ScambosFricker and Scambos (2009). Our improved estimation of lake-volume change extended the time series and improved their degree of monotonicity in most cases. One exception is SLM, which is a challenging case from its high number of missing tracks and relatively short period between drainage cycles. An example of the improvement is shown at SLE (Figs 5a and 6e), where one short episode of apparent filling and draining was eliminated as an artifact of the previous track-averaging schemes. The new time series indicates that SLE drained consistently from November 2003 until reaching minimum volume in June 2006, at which time filling resumed and continued monotonically until at least October 2009. At SLW, our method resolved slightly more about how the lake level changed between November 2006 and October 2007 than Reference Fricker and ScambosFricker and Scambos (2009), although the difference was negligible (Fig. 6c).
More significant changes in the volume time series were found for upstream lakes of the system, such as K5, K1 and W6; K5 and W6 reach minimum volumes nearly a year later than in the previous time series of Reference Smith, Fricker, Joughin and TulaczykSmith and others (2009). This improvement mainly came from the inclusion of TOOs, which improve the spatial and temporal resolution (Fig. 5e and f).
Volume change estimates for many lakes, including USLC, K1, K5, W6, W7 and W8, are compromised by an absence of data from later campaigns, in particular USLC and K1, both of which have no data after the Laser 3k (October 2008). While USLC appears to have started filling following a low stand in March 2008 (Fig. 6d), behavior of K1 during this time frame is less certain. While it appears to have been filling up to March 2008, given the vertical ranges of other lakes in the region, it is possible that K1 began draining shortly thereafter. While initially we assume it to be filling after this time, we address the uncertainty with a separate model run, in which it drains after March 2008 (Section 5.5). Additional complications arise for volume time series for K1 and K5 associated with high rates of regional surface uplift (approaching 1 ma–1; Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009), which complicated efforts to separate out lake surface elevation change. In other areas (Section 3.2.2), regional surface change appears to be significant for the routing of water between the various lakes in our study.
3.2.2. ICESat-derived elevation change over the flow trifurcation region
Our dh/dt analysis revealed surface uplift in the trifurcation region that was widespread and heterogeneous. Most of the region thickened by at least 1 m between 2003 and 2008, with thickening of >2 m directly over the subglacial hydrologic pathway where flow diverted from northern WIS and KIS is apportioned between SLC and SLM in the Southern Basin and SLW in the Central Basin (Fig. 7). The trifurcation region is on the edge of a larger region of thickening associated with the slowdown of lower WIS and increased buttressing of the Ross Ice Shelf in this region (Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012; Reference Scheuchl, Mouginot and RignotScheuchl and others, 2012), and has been reported in previous studies (e.g. Reference Bindschadler, Vornberger and GrayBindschadler and others, 2005; Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009). None of these studies considered the effect of this change or its heterogeneities on the flow of subglacial water. In the following sections we explore how change in the trifurcation region may relate to the observed lake activity.
4. Modeling of the Whillans Ice Plain Subglacial System
The new time series for lake volume change for the region, along with our calculated surface change, provides the tools necessary to describe the evolving hydrology over the region and infer a source region for SLW. We show how the subtle observed changes over a small but critical area cause the system to reorganize.
4.1. Subglacial water model
We used the subglacial water model of Reference CarterCarter and others (2011) and Reference Carter and FrickerCarter and Fricker (2012) to simulate basal water distribution and subglacial lake filling rates in WIP and the hydrologic catchments that supply it. The model directs subglacial meltwater down the hydrologic potential and calculates the change in water distribution from the storage and release of water by previously documented subglacial lakes (Fig. 3). Our model does not address the physical mechanisms of flood initiation and termination (e.g. Reference Evatt, Fowler, Clark and HultonEvatt and others, 2006; Reference FowlerFowler, 2009), nor does it attempt to identify a particular flow mechanism (i.e. focused conduits versus distributed sheets or cavities; e.g. Reference Flowers and ClarkeFlowers and Clarke 2002; Reference Creyts and SchoofCreyts and Schoof, 2009). It simply tests whether the hydrologic potential surface, basal melt distribution and lake volume change are consistent with each other. Once this is established, we confirm hydrologic connections between the various subglacial lakes in our domain and infer water sources for each of the lakes. To improve the model’s water routing we have increased the spatial resolution from 5 km × 5 km to 2.5 km × 2.5 km.
Water is routed via a simple ‘steady-state D8’ (Reference Quinn, Ostendorf, Beven and TenhunenQuinn and others, 1998), in which flux out of a cell is equal to the sum of incoming flux and local melt:
where Q
in and Q
out are water flux in and out of a cell, respectively, 
is the basal melt rate (negative, if water is freezing to the base), and Δx and Δy are the cell’s horizontal dimensions. Q
out is apportioned among all ‘downstream’ cells using the formula
with k being the number of adjacent cells with lower hydrologic potential, and ds the distance to the adjacent cell 
such that more flow goes toward adjacent cells with steeper downward gradients.
For cells lying within subglacial lakes, we apply additional forcing depending on whether a lake was inferred to be filling or draining during each time interval (3–4 months, defined by the timing of the ICESat campaigns (Section 3)). For intervals when a lake was inferred to be filling, the corresponding cells are treated as hydrologic sinks, setting Q out to zero. We then compare the flux of water into the lake predicted by the model against the observed volume increase. When lake volumes are decreasing, we add the inferred volume loss to the melt rate for cells within this lake and allow water flowing in from upstream to pass through as described in Eqns (3) and (4). In summary, observations of lake volume increases are used to validate our model while observations of lake volume decreases are used to force it.
4.2. Model inputs
Our water model requires four inputs: (1) a basal melt rate to produce water for the system; (2) a time series of lake volume change to control the water supply and compare to model output; (3) a hydropotential grid to determine where water is routed; and (4) an estimate for the variation in the hydrologic potential over time. For (1), we use published modeled melt rates (Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004), which are based on an inversion for basal shear stress from observed topography and velocity of the Siple Coast ice streams having a basin-wide error of ±10%. The published error may also have some systematic biases due to its being based on older datasets for ice surface and bedrock topography (e.g. Reference Lythe and VaughanLythe and others, 2001 vs Reference FretwellFretwell and others, 2013), and its exclusion of longitudinal stress effects (e.g. Reference Beem, Jezek and Van der VeenBeem and others, 2010). Quantifying the effect of these differences on the actual melt rates is, however, outside the scope of this study. For (2) and (4), we use the two ICESat time series described in Section 3. We created (3) following the Appendix, from a published DEM of the ice surface z srf (Reference Haran and ScambosHaran and Scambos, 2007), ice thickness data obtained by airborne radio-echo sounding (RES) campaigns from 1971 to 1999 (Drewy, 1975; Reference Retzlaff, Lord and BentleyRetzlaff and Bentley, 1993; Reference BlankenshipBlankenship and others, 2001) and the recently published Bedmap2 ice thickness interpolation for h i (Reference FretwellFretwell and others, 2013). Although these measurements were used in the final Bedmap2 grid, the smallest resolvable feature in that project was on the order of 5 km while the along-track resolution of some of these datasets approached 30 m. Following work by Reference CarterCarter and others (2011) we have found that the addition of the original ice thickness measurements allows us to incorporate the effects of smaller-scale features in the ice and bedrock topography that likely is important for basal water routing.
From the initial hydropotential grid, we apply a stream-etching algorithm (Appendix) to incorporate the possible effects of narrow flow paths and optimize the high along-track resolution (30 m) of the RES ice thickness profile data. In the resulting hydropotential map, the combined uncertainties from ice surface elevation and ice thickness lead to an error of up to ±2 m in the calculation of the final hydrological potential surface where coverage is reasonably dense (Reference Bindschadler, Vornberger and GrayBindschadler and others, 2005), but may be up to ±10 m across large gaps in the sampling (Reference FretwellFretwell and others, 2013). Translating this uncertainty in ice-sheet geometry to uncertainty in subglacial water flux is more complicated than translating uncertainty to basal melt rates or lake volume change. Although stream etching makes this less likely, at least a few locations exist where small relative changes in the hydrologic potential may cause substantial redistribution in the water downstream. To address this possible instability, we developed a tuning method, described in the following subsection.
4.3. Tuning the bed topography
We initially run the model using the hydropotential output as described in the Appendix, to determine if reproducing the observed filling rates for each of the subglacial lakes (Fig. 6; supplementary material at http://www.igsoc.org/hyperlink/13j085) without modifications is possible. If a significant misfit exists, we compare our modeled water distribution with the lake location and hydropotential map to identify where realistic adjustments to the relative bed topography might redirect more or less water toward a given lake. All effort is taken to ensure that such adjustments are within the error limit of the ice thickness measurements, with larger adjustments allowed in areas where the bed topography is less well constrained. The model is rerun to calculate the impact of the adjustment on the target lake and lakes downstream, until modeled inflow for each lake is as close to observed inflation rates as could be reasonably accommodated. Around W6 and W7, we adjust the hydropotential by up to ∼5 m – equivalent to adjusting the ice thickness by ∼50 m, which lies within the stated combined uncertainty of measurements and gridding errors of ±59 m for cell-containing data in this region (Reference FretwellFretwell and others, 2013). Around W8, where input data are sparse and surrounding topography contains considerable relief, >20 m of adjustment to the elevation is necessary. The most critical adjustment is to the ‘flow-path trifurcation’ region (Figs 2 and 3) in the lower WIS trunk where data are sparse and previous reports (Reference Carter and FrickerCarter and Fricker, 2012) show evidence for diverging water flow (Fig. 3). Through this region in which the observed change in ice thickness (Section 3.1) could plausibly lead to switching of water flow between SLC (Southern Basin) and SLW (Central Basin) (Fig. 3), we construct a bed surface that does not deviate >100 m from the Bedmap2 interpolation.
4.4. Categorization of out-of-balance lakes
Reference Carter and FrickerCarter and Fricker (2012) report several lake-like features in their model for the Siple Coast for which no reasonable adjustment to the hydropotential surface upstream reproduces the observed filling rate, dividing them into four classes: ‘overfilling’, ‘underfilling’, ‘leaking’ and ‘through-flowing’. In this study, we retain this classification scheme for the purposes of evaluating the performance of the model and working toward a more realistic treatment of lake dynamics. ‘Overfilling’ and ‘underfilling’ lakes are lakes for which the modeled inflow is either significantly greater than or less than the observed filling rate for the lake. Although significant, these discrepancies do not comprise a significant fraction of the regional water budget, so the model treats the lakes as if they are in balance. In contrast, ‘leaking’ and ‘through-flowing’ lakes require significant adjustments to the way in which the final model handles them. Leaking lakes are defined as lakes for which the modeled inflow is in excess of the observed volume increase, especially as it approaches maximum volume, lying directly upstream of a lake for which the modeled inflow is substantially less than the observed volume increase. The leaking may result from one of two situations: (1) if the lake reaches maximum volume approximately halfway between two campaigns and begins draining, it might not drain below the level it was at during the previous campaign, and thus appear to be filling (i.e. aliasing); (2) if initial inflow exceeds outflow, the lake continues to increase in volume even as water escapes downstream. For either case, we include an additional parameterization: any modeled inflow in excess of the observed inflow is directed to the downstream side of the lake and added to the melt rate of the nearest gridcell downstream. In contrast, ‘through-flowing’ lakes are water bodies whose volumes appear to correlate with water flux from upstream, but do not appear to store a significant amount of that flux. Rather than lakes, these appear to be temporary storage during high flow. In the final model run (hereafter referred to as ‘main run’), such lakes are no longer treated as sinks or sources, but instead water is allowed to simply pass through them.
4.5. Model implementation
In addition to the main model run, we perform four variations: (1) a ‘no lakes’ run in which we do not include any subglacial lake dynamics, calculating water distribution in a manner similar to the steady-state model of Reference Le Brocq, Payne, Siegert and AlleyLe Brocq and others (2009); (2) a ‘static surface’ run, in which we run the main model but hold the hydropotential constant over the model run duration to test the significance of secular changes in ice thickness over the model time domain; (3) a ‘Whillans and Mercer only’ run, in which all flow from upper KIS is retained by a synthetic ‘lake’ just upstream of the confluence of waterways draining KIS and northern WIS, designed to test the impact of the hypothesized piracy of water from upper KIS to lower WIS (Reference Anandakrishnan and AlleyAnandakrishnan and Alley, 1997); and (4) the ‘K1 floods’ run, in which we allow K1 to drain after 2008. By comparing the inferred lake-filling rates for the lakes in lower WIS from each of these separate model runs with the observed volume change, we further constrain how activity in WIS trunk and further upstream impacts downstream hydrology.
4.6. Model results
4.6.1. Variations in the supply of meltwater from upper MIS, WIS and KIS
Between 2003 and 2009, water supplied to WIP from KIS, WIS and MIS varied between 12 ± 1.8 and 37 ± 4.2 m3 s–1 (1 m3 s–1 ≈ 0.0316 km3 a–1). The primary source of this variability came from the filling and draining of lakes K5 and K1 (with up to 5 ± 0.7 m3 s–1 supplied by the K5 branch of KIS and up to 15 ± 2.3 m3 s–1 supplied by the K1 branch (Figs 8 and 9)). For most of the study period, however, modeled inflow into WIP was 20–25 m3 s–1. This estimate for inflow was confirmed by the combined observed volume increase for SLE, SLW, SLC and SLM of 21 ± 2.7 m3 s– 1 for a total volume increase of 0.25 ± 0.033 km3 between October 2007 and March 2008 (Fig. 8c–f). The most significant uncertainty appears to be associated with the basin supplying K1. Previous models (e.g. Reference Johnson and FastookJohnson and Fastook, 2002; Reference Le Brocq, Payne, Siegert and AlleyLe Brocq and others, 2009) implied the presence of a flow bifurcation point upstream of K1 that directed water toward Bindschadler Ice Stream (BIS). Our model did not include this feature, so our estimates for flow toward WIP when K1 was draining may be high. The following subsections discuss in more detail how water was apportioned among the individual basins and how that may have varied over time.
Comparison between modeled lake volume change rates and lake volume change rates inferred from observations (see also Fig. 6). For each system of lakes we select results from the model runs most relevant to that system. Orange lines denote lake volume after leakage allowed.
Modeled apportioning of water flow from northern WIS/KIS among the three drainage basins of WIP. Solid line denotes ‘main run’, round dots denote ‘no Kamb’ run, square dots denote ‘static surface’ run and dashed line denotes ‘K1 floods’ run. Note that the ‘static surface’ run does not route any water to SLW.
4.6.2. SLW and the Central Basin
The variation of inflow into SLW over the 6 year observation period provides evidence for significant changes in flow routing and discharge upstream. Between November 2003 and November 2005, the volume of SLW increased at a mean rate of 0.11 ± 0.02 m3 s–1 (on the same order of magnitude as our modeled inflow of 0.25 ± 0.11 m3 s–1 (Fig. 8e)). This was consistent with the lake receiving water from only local sources downstream of the trifurcation region (Fig. 7a and b). During this period, our model predicted that water from upper WIS and KIS flows to the Northern and Southern Basins, bypassing the Central Basin. Beginning in November 2005, the observed rate of volume increase for SLW increased by over an order of magnitude to 5.6 ± 0.6 m3 s–1, reasonably close to the modeled inflow rate of 3.6 ± 0.7 m3 s–1. In the model, this change was primarily caused by a change in the hydropotential surface of the trifurcation region that results in a redirection of the water crossing WIS trunk away from SLC and the Southern Basin toward SLW and the Central Basin (Figs 7c and d and 10). Although the total magnitude of ice thickness change over the model’s time frame is <2.5 m, the relief of the ice and bed surface is sufficiently low that a small change results in the rerouting of a substantial portion of the discharge crossing from the northern margin of WIS trunk away from SLC and toward SLM. Without this change (i.e. the static surface run) SLW either received water for the entire observation period or did not fill at all, with water going instead toward SLC and SLM.
4.6.3. USLC, SLC, SLM and the Southern Basin
The observed surface changes in the trifurcation region with lake drainage in upper KIS were also significant in explaining the observed volume changes of USLC, SLC, SLM and L78 as being influenced by ongoing changes both in the trifurcation region and further upstream in the WIS and KIS basins. USLC, the uppermost lake in this basin, was supplied entirely by the southern WIS catchment (Fig. 3), and the modeled inflow varied between 2.6 ± 0.7 and 13.2 ± 2.1 m3 s–1 depending on whether its upstream neighbor W6 was filling or draining. During the first filling cycle, the mean modeled inflow rate (7.1 ± 1.1 m3 s–1) was close enough to the observed rate of volume increase (6.1 ± 1.2 m3 s–1) to confirm that water from southern WIS was not traveling to SLC during this time. The filling of SLC, immediately downstream, therefore requires an additional source of water, primarily the ∼8.4 m3 s–1 water crossing the WIS trunk from northern WIS at this time (Figs 7a and b, 8d and 10a and b). Consequently, the volume change time series for SLC and USLC are exceptional evidence for both water piracy between KIS and WIS and also flow of water across the lower WIS trunk between November 2003 and June 2005 when SLW was not undergoing significant volume change (Fig. 8b and d).
The volume change time series for SLC and SLM provides interesting clues about the nature of the subglacial drainage system. The modeled inflow for SLC between November 2004 and June 2005 exceeds the observed volume increase by 4.7 ± 2.4 m3 s–1, while the modeled inflow for SLM falls below the observed volume increase rate by 7.5 ± 1.3 m3 s–1 (Fig. 8g). This imbalance is consistent with water leaking out of SLC before it reaches maximum volume. Close connections of this nature were previously observed in MacAyeal Ice Stream (MacIS) where two lakes originally believed to be separate actually behaved as a single, larger combined lake (Reference CarterCarter and others, 2011). We explore mechanisms by which this occurred in Section 5.2.
We obtain further validation of our model for flow upstream of L78, the last lake in the system. This lake immediately downstream of SLM appears to swell and subside in response to the passage of water through the system (Fig. 8h), but does not retain water. With the new time series, modeling the water flow though this feature correlates better with observed surface uplift than it had in previous models (Reference Carter and FrickerCarter and Fricker, 2012)
4.6.4. SLE and the Northern Basin
Our observations of volume change for SLE complete our accounting for all water produced upstream of WIP and represent the best evidence to date for water piracy between the upper KIS and lower WIS. ICESat observations show that SLE began draining around March 2003 (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007) and continued draining until June 2006 (Figs 8c and 10a–c). The observed rate of volume change was 7.1 + 0.7 m3 s– 1 between June 2006 and March 2008, and then increased to 19.6 ± 2.0 m3 s– 1 until the end of the ICESat observation period (September 2009). The main model directs 5.7 ± 1.3 m3 s– 1 (0.18 ± 0.04 km3 a– 1) to the lake during the early part of the filling cycle, closely agreeing with observations, and stayed constant through to September 2009. Consequently, the main model is unable to reproduce the increase in filling rate observed at the end of the observation period. We were only able to reproduce the observed increase in the filling rate if K1 was set to begin draining in March 2008 (the K1 floods run). In this case, the modeled rate increased to 13.6 ± 2.9 m3 s– 1, providing the best evidence of K1 draining and additional support for the hypothesized piracy of water from upper KIS. For the ‘Whillans and Mercer only’ run, the filling rate was reduced by 0.5–4.5 m3 s–1 relative to the main model run.
5. Implications
5.1. Sensitivity of WIP lakes to changing conditions upstream
The ice-sheet geometry that currently directs water from KIS to WIS and to the three basins of WIP has undergone (and continues to undergo) dramatic changes. KIS, WIS, BIS and nearby MacIS have all shut down and restarted in the last millennium (Reference Catania, Hulbe, Conway, Scambos and RaymondCatania and others, 2012). Our results have provided a sense of the magnitude and (surprisingly) short time frame over which dramatic hydrological switches related to ice-flow switches occur. With only a few meters of ice thickness change over a few years, water that originally flowed from the northern margin of WIS across the trifurcation region to SLC and SLM began flowing to SLW. Ice flow and geometry in WIP continues to undergo spatially heterogeneous rates of slowdown and thickening (Reference Bindschadler, Vornberger and GrayBindschadler and others, 2005; Reference JoughinJoughin and others, 2005; Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012; Reference Scheuchl, Mouginot and RignotScheuchl and others, 2012). Our work suggests that the subglacial hydrology is simultaneously changing. In some locations, these changes were quite rapid, calling into question the long-term stability of many of the active lakes in the Siple Coast region (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007), as their presence requires a source of water and a hydropotential minimum.
5.2. Implications for more detailed hydrologic models of the ice sheet
Although our model is relatively simple, we use our results to begin to infer certain qualities about the mechanisms by which lakes drain and fill. Whereas previous modeling efforts for the region (e.g. Reference AlleyAlley, 1996; Reference Johnson and FastookJohnson and Fastook, 2002) focused exclusively on distributed systems, our work shows evidence for both distributed and channelized systems. The presence of several flow bifurcations or even trifurcations and the low threshold for flow diversion is consistent with a distributed system under a pressurized lid, such as that reported by Reference Catania and PaolaCatania and Paola (2001), or the linked cavities observed under WIS trunk and portions of lower KIS (Reference Kamb, Alley and BindschadlerKamb, 2001; Reference EngelhardtEngelhardt, 2004). Likewise, the presence of ‘through-flowing lakes’ (e.g. L78) suggests a distributed system in which water pressure and discharge are correlated.
Most of the lake drainage events in our study, however, appear to be more consistent with conduit R-channels (Reference RöthlisbergerRöthlisberger, 1972; Reference NyeNye, 1976) that form only when supplied by a reservoir large enough to maintain an ever-increasing discharge (Reference AlleyAlley, 1996). The initiation of drainage before the lake reaches the flotation height and the ongoing decline in water level during the lake drainage event requires a system where water pressures are somewhat below overburden pressure. The heat necessary to enlarge such a channel is generated by turbulent dissipation of water moving down the hydraulic gradient. Given low hydraulic gradients typical of the Siple Coast, we would expect these channels to take many months to reach capacities on the order of tens of m3 s–1,such that if a lake has sufficiently large inflow we will observe evidence for outflow in lakes downstream of it before it reaches maximum volume, as outflow will initially be less than inflow. Specifically, the observed dynamics of leaking lakes (e.g. SLC and SLM) provides the strongest evidence to date for a channelized system forming in a low hydraulic gradient. Modeling work by Reference Evatt, Fowler, Clark and HultonEvatt and others (2006) and Reference Carter, Blankenship, Young, Peters, Holt and SiegertCarter and others (2009) showed that for environments similar to the Siple Coast lake systems, the duration of floods via conduits under ice sheets was on the order of years to months. Understanding the respective roles of distributed and channelized systems’ water transport is as critical to understanding the subglacial hydrology of Antarctica as it is to Greenland (e.g. Reference Sundal, Shepherd, Nienow, Hanna, Palmer and HuybrechtsSundal and others, 2011). This is especially important given the degree to which a subtle shift in pressure associated with the transition between the two systems determined the direction of water flow.
We have shown that WIP has a low threshold for nonlinear changes in its subglacial water system, i.e. small changes in ice thickness can lead to large relative changes in hydropotential. Although incidences of water piracy have been suggested throughout Antarctica in response to relatively small relative changes in ice geometry (e.g. Reference Evatt, Fowler, Clark and HultonEvatt and others, 2006; Reference Wright, Siegert, Le Brocq and GoreWright and others, 2008; Reference Le Brocq, Payne, Siegert and AlleyLe Brocq and others, 2009; Reference Livingstone, Clark and WoodwardLivingstone and others, 2013; Reference SiegertSiegert and others, 2013), KIS and WIS are unique among modern drainage basins in Antarctica in that they have exceptionally low regional hydraulic gradients, very low basal relief and are undergoing rapid thickening of up to 1 m a−1 (Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009). In such an environment, water flow paths have very little bedrock constraint and the most efficient route down is constantly changing. Developing a predictive model of future ice discharge that accurately reproduces water piracy through the greater Siple Coast will close observations of secular changes in the ice surface and subglacial hydrology with observations of the hydrologic system for validation.
Although a number of studies have focused on how changing oceanographic conditions at the grounding line rapidly propagate inland through dynamic thinning and retreat, far fewer have explored how ice thickness change in one catchment affects neighboring catchments. Consequently, as changes to ice geometry propagate inland, they begin to affect adjacent basins through changes to the basal hydrology and redistribution of subglacial water. The interplay of the subglacial water system between adjacent basins as inferred for the Siple Coast (Reference Anandakrishnan and AlleyAnandakrishnan and Alley, 1997) is thought to have occurred elsewhere (e.g. Reference Wright, Siegert, Le Brocq and GoreWright and others 2008; Reference Livingstone, Clark and WoodwardLivingstone and others, 2013).
5.3. Implications for potential findings at the WISSARD drill site
Our study period (2003–09) is short compared to the ice flow history that experienced high variability on multiple timescales, from millennia (Reference Conway, Hall, Denton, Gades and WaddingtonConway and others, 1999) to centuries and decades (Reference Catania, Hulbe, Conway, Scambos and RaymondCatania and others, 2012). Rapid sliding of WIS started and stopped at least once during the past 1000 years. Our study shows that subglacial water-flow routing and distribution varies on shorter timescales. Any sediment present at SLW likely records some of that variability. Reference Christoffersen and TulaczykChristoffersen and Tulaczyk (2003) indicated that the presently occurring freezing at the base in lower WIS tends to induce erosion of the subglacial sediments, evidenced by packages of accreted sediment in boreholes drilled into KIS (Reference EngelhardtEngelhardt, 2004). The only areas not directly exposed to sediment erosion in this environment are the beds of subglacial lakes. Although the exact age of SLW is unclear, Reference Christianson, Jacobel, Horgan, Anandakrishnan and AlleyChristianson and others (2012) and Reference HorganHorgan and others (2012) have inferred a depth of 6 m for SLW (actually found to be ∼2 m by borehole access of the lake in January 2013 (personal communication from WISSARD field team, 2013)). Given the estimated pre-2005 inflow rate of 0.11 ± 0.1 m3 s–1, SLW is possibly less than a few decades old. The coring of the sediments and water chemistry may reveal many more clues regarding the relative recurrence of water-flow switching. Material within the lake may originate as far away as upper KIS and may provide information about the evolution of subglacial environments upstream. The subglacial hydrology feeding SLW is changing visibly on a timescale of years and decades, so we expect the WISSARD project to recover a record of biology, water chemistry and sediment history that is recent and constantly reworked.
6. Summary
We updated the ICESat-based lake volume change time series for WIS, MIS and KIS to cover the entire ICESat operational period (2003–09). We showed that the inferred hydrologic system is the ephemeral product of a region that has been undergoing progressive change for over a century. The volume of water entering lakes in lower WIS/MIS appears to require the entire meltwater budget estimated for these basins and much of the meltwater budget of upper KIS. The flow of this meltwater over time is regulated by a system of lakes distributed in the upper part of these catchments, resulting in a distribution that varies by a factor of four over the short time frame of our study (Fig. 10). Hydrologic variability over time also results from observed changes in surface elevation in the lower WIS trunk where water from northern WIS and KIS spreads out and is distributed among the three basins. In this area, localized surface inflation of ∼2 m in one key location causes a diversion of water ∼4 ± 0.9 m3 s–1 away from SLC and toward SLW resulting in a ∼20-fold increase in the filling rate of SLW. Given the widespread heterogeneous changes in ice thickness and surface velocity prevalent throughout the downstream sections of many of Antarctica’s ice streams and outlet glaciers, water piracy may be quite common in subglacial water systems, and may play a critical role in the coevolution of adjacent ice catchments. Ongoing efforts to sample SLW water, biota and the sedimentary biota beneath will likely sample material derived from a wide source region and capture a relatively short record of considerable degree of change in the system that supplied it.
Water flux distribution for select time-steps over the wider area (see box in Fig. 2 for location), focusing on areas affected by lake activity. Blue lakes are filling, yellow lakes are draining, and magenta lakes are overflowing.
Acknowledgements
ICESat data analysis was funded through NASA award NNX07AL18G. RES data collection was funded by NASA grant NNX08AN68G. The University of Texas Institute for Geophysics’ Support Office for Aerogeophysical Research, through which many of the RES data in this paper were collected, was supported by US National Science Foundation (NSF) grants OPP-9120464 and OPP-9319369. Funding during the analysis and writing was provided by the Scripps Institution of Oceanography Postdoctoral Program and by the Institute of Geophysics and Planetary Physics of Los Alamos National Laboratory LLC under subcontract No. 73593-001-09. We thank Ben Smith, Geir Moholdt and Fernando Paulo. Additional proofreading was performed by Judy Gaukel. We also thank Poul Christoffersen and an anonymous reviewer.
Appendix: Hydrologic Potential Grid
Subglacial water transport through a region is governed by the hydropotential (θh ); in Antarctica the hydropotential is typically calculated assuming water pressure is equal to the overburden pressure, giving
where ρ w and ρ i are the densities of water (1000 kg m–3) and ice (917 kg m–3), respectively, and g is acceleration due to gravity (9.81 m s–2), z srf is the measured surface elevation and h i is the ice thickness. This equation implies that surface elevation gradients are 11 times more significant than bed gradients for determining the direction of water flow (Reference ShreveShreve, 1972). Bedrock gradients, however, can locally exceed surface gradients by more than a factor of 11.
Ice thickness for the WIP region has been measured multiple times over the past five decades (Reference Bindschadler, Vornberger and GrayBindschadler and others, 2005), beginning with the International Geophysical Year (IGY) seismic traverses of 1956–57 (Reference Crary, Robinson, Bennett and BoydCrary and others, 1962), and continuing to the present day (Christianson and others, 2005). Our primary source of ice thickness came from the BEDMAP ice thickness database (Reference Lythe and VaughanLythe and others, 2001) and the US Support Office for Aerogeophysical Research (SOAR) archives (Reference BlankenshipBlankenship and others, 2001). Most of these RES ice thickness measurements have an accuracy of 20 m, with along-track sampling approaching every 30 m and track separation no greater than 5.5 km throughout much of the region. Our study area included one large region (40 km × 70 km) around USLC where no RES data were available at the time of writing (Fig. 2b). Unfortunately, this critical region is where the flow path from the northern WIS trifurcated to form the flow path entering the southern and central hydrologic basins. For this area and several smaller regions we rely upon the interpolation work done as part of the recently updated Bedmap2 dataset (Reference FretwellFretwell and others, 2013).
We calculated an initial hydrologic potential surface based on interpolated ice thickness data (Drewy, 1975; Reference Retzlaff, Lord and BentleyRetzlaff and Bentley, 1993; Reference BlankenshipBlankenship and others, 2001) and the ice surface data (Reference Haran and ScambosHaran and Scambos, 2007). Subglacial water flow may be controlled by features at length scales of a few hundred meters, which are not resolved by our grid (2.5 km). To incorporate the possible effects of narrow flow paths and optimize the high along-track resolution (30 m) of the RES ice thickness profile data, we traced several prominent drainages by following a contour map of the initial hydrologic potential using the method of Reference Dupain-TrielDupain-Triel (1791). We then located the intersection between these drainages and local minima in the hydrologic potential along every RES profile. We used a piecewise cubic hermitic polynomial to interpolate the geometry of the drainage path between these intersection points. The hydrologic potential of cells containing portions of this stream was set to the hydrologic potential of the points along the drainage path rather than averaging all measurements within the cell. This effectively etched a stream path into the hydrologic potential grid (Reference Saunders, Maidment and DjokicSaunders and others, 2000) (Figs 2 and 3).
