Evidence of rapid subglacial water piracy under Whillans Ice Stream, West Antarctica

Abstract

The subglacial water system of lower Whillans Ice Stream on the Siple Coast, West Antarctica, contains numerous connected subglacial lakes in three hydrological basins (northern, central and southern). We use Ice, Cloud and land Elevation Satellite (ICESat) data to derive estimates of lake volume change and regional thickness changes. By combining these results with a water budget model, we show that a uniform, localized thickness increase perturbed the hydropotential, resulting in a change in course of a major flow path within the system in 2005. Water originating from upper Whillans and Kamb Ice Streams that previously supplied the southern basin became diverted toward Subglacial Lake Whillans (SLW). This diversion led to a tenfold filling rate increase of SLW. Our observation suggests that water piracy may be common in the Siple Coast region, where the gentle basal relief makes the basal hydropotential particularly sensitive to small changes in ice thickness. Given the previously inferred connections between water piracy and ice-stream slowdown elsewhere in the region, the subtle and complex nature of this system presents new challenges for numerical models.

1. Introduction

There are 379 known subglacial lakes under the Antarctic ice sheet, which have been inferred through a variety of remote-sensing techniques (Wright and Siegert, 2012). Before 2005, airborne radar sounding (e.g. Oswald and Robin, 1973; Siegert and others, 2005; Carter and others, 2007; Wright and Siegert, 2012) was the primary tool used to identify subglacial lakes. Recently, satellite measurements of surface elevation change (e.g. Gray and others, 2005; Wingham and others, 2006; Fricker and others, 2007; Smith and others, 2009) and satellite image differencing (e.g. Fricker and others, 2007; Fricker 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 km³ on monthly to decadal timescales (Fricker and others 2007; Fricker and Scambos, 2009; Smith and others, 2009). In one location (Byrd Glacier), a lake drainage event was concurrent with temporary acceleration of ice downstream (Stearns 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 (Carter and others, 2011; Carter 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) (Fricker and others, 2011; Priscu 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 (Bindschadler and others, 2005) and likely has undergone much larger-scale flow reorganization in the region over the last millennium (Catania 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.

Fig. 1. 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 Rignot 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 (Rignot and others, 2011) despite low surface and basal slopes (generally <0.001), low driving stresses (Bindschadler and others, 2005), and widespread freezing conditions at the ice base with accretion rates approaching 5–6 mm a^–1 (Joughin and others, 2004) (Fig. 2). Although small quantities of meltwater are generated by frictional heating over regions of enhanced basal traction or ‘sticky spots’ (Stokes and others, 2007; Sergienko and Hulbe, 2011), most of the basal lubrication responsible for sustaining such high flow rates comes from basal meltwater generated upstream (Bougamont and others, 2003; Parizek 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 (Anandakrishnan 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 (Joughin and others, 2005; Scheuchl 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 (Pritchard and others, 2009).

Fig. 2. 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 Haran 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 (Shreve, 1972). In this framework (Appendix), surface elevation gradients are 11 times more significant than bed gradients in determining the direction of water flow (Shreve, 1972; Iken 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 (Shabtaie and Bentley, 1987; Fricker and Scambos, 2009), and connected lake events inferred from ICESat data (Figs 2 and 3) (Fricker and others, 2007; Fricker 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).

Fig. 3. 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).

Carter 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 (Fretwell 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 Fricker and others (2009) and Smith and others (2009) reported additional lakes in our study area, work by Carter 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; Siegfried 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 (Fricker and others, 2007; Fricker and Scambos, 2009). Previous studies on lake volume change (e.g. Fricker and Scambos, 2009; Smith and others, 2009) and regional surface elevation change (e.g. Pritchard 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 Fricker and Scambos (2009), while for (2) we constrained our analysis to the trifurcation region where our understanding of the hydropotential surface from previous work (Carter 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. Fricker and Padman, 2006). We also applied a correction for inter-campaign biases (personal communication from T. Urban, 2013, updated from Urban and Schutz, 2005; Siegfried 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 (Borsa 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 (Moholdt and others, 2010).

Fig. 4. 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 Fricker and Scambos (2009) and Smith and others (2009). Collect all available repeats along those tracks and identify the limits of the along-track elevation change anomalies as described by Fricker and others (2007).

• 2. Within these limits, calculate the average elevation for each repeat (Fricker and others, 2007), applying the corrections identified in Section 3.1 and an additional correction for regional ice thickness change (Smith 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, Δh_k is the average elevation anomaly across the lake for track k, l_k 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 Fricker and others (2007), Fricker and Scambos (2009) or Smith 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).

Fig. 5. 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.

Fig. 6. Volume time series for all lakes in this study. Previous time series (Smith and others, 2009, for K1, K5, W6, W7 and W8; Fricker 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 (Fricker 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% (Smith and others, 2009) while the errors for the remaining lakes were closer to ±20% as their shorelines were better constrained (Fricker 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.

Fig. 7. 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 Ligtenberg 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 (Pritchard and others, 2012). From these assumptions the relation between surface elevation change and hydropotential change is

(2)

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 Fricker 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 Fricker 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 Smith 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; Pritchard 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 (Pritchard and others, 2012; Scheuchl and others, 2012), and has been reported in previous studies (e.g. Bindschadler and others, 2005; Pritchard 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 Carter and others (2011) and Carter 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. Evatt and others, 2006; Fowler, 2009), nor does it attempt to identify a particular flow mechanism (i.e. focused conduits versus distributed sheets or cavities; e.g. Flowers and Clarke 2002; Creyts 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’ (Quinn and others, 1998), in which flux out of a cell is equal to the sum of incoming flux and local melt:

(3)

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

(4)

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 (Joughin 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. Lythe and others, 2001 vs Fretwell and others, 2013), and its exclusion of longitudinal stress effects (e.g. Beem 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 (Haran and Scambos, 2007), ice thickness data obtained by airborne radio-echo sounding (RES) campaigns from 1971 to 1999 (Drewy, 1975; Retzlaff and Bentley, 1993; Blankenship and others, 2001) and the recently published Bedmap2 ice thickness interpolation for h _i (Fretwell 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 Carter 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 (Bindschadler and others, 2005), but may be up to ±10 m across large gaps in the sampling (Fretwell 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 (Fretwell 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 (Carter 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

Carter 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 Le 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 (Anandakrishnan 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 m³ s^–1 (1 m³ s^–1 ≈ 0.0316 km³ 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 m³ s^–1 supplied by the K5 branch of KIS and up to 15 ± 2.3 m³ 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 m³ s^–1. This estimate for inflow was confirmed by the combined observed volume increase for SLE, SLW, SLC and SLM of 21 ± 2.7 m³ s^{– 1} for a total volume increase of 0.25 ± 0.033 km³ 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. Johnson and Fastook, 2002; Le 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.

Fig. 8. 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.

Fig. 9. 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 m³ s^–1 (on the same order of magnitude as our modeled inflow of 0.25 ± 0.11 m³ 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 m³ s^–1, reasonably close to the modeled inflow rate of 3.6 ± 0.7 m³ 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 m³ 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 m³ s^–1) was close enough to the observed rate of volume increase (6.1 ± 1.2 m³ 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 m³ 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 m³ s^–1, while the modeled inflow for SLM falls below the observed volume increase rate by 7.5 ± 1.3 m³ 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 (Carter 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 (Carter 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 (Fricker 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 m³ s^{– 1} between June 2006 and March 2008, and then increased to 19.6 ± 2.0 m³ s^{– 1} until the end of the ICESat observation period (September 2009). The main model directs 5.7 ± 1.3 m³ s^{– 1} (0.18 ± 0.04 km³ 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 m³ 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 m³ 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 (Catania 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 (Bindschadler and others, 2005; Joughin and others, 2005; Pritchard and others, 2012; Scheuchl 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 (Fricker 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. Alley, 1996; Johnson 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 Catania and Paola (2001), or the linked cavities observed under WIS trunk and portions of lower KIS (Kamb, 2001; Engelhardt, 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 (Röthlisberger, 1972; Nye, 1976) that form only when supplied by a reservoir large enough to maintain an ever-increasing discharge (Alley, 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 m³ 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 Evatt and others (2006) and Carter 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. Sundal 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. Evatt and others, 2006; Wright and others, 2008; Le Brocq and others, 2009; Livingstone and others, 2013; Siegert 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⁻¹ (Pritchard 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 (Anandakrishnan and Alley, 1997) is thought to have occurred elsewhere (e.g. Wright and others 2008; Livingstone 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 (Conway and others, 1999) to centuries and decades (Catania 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. Christoffersen 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 (Engelhardt, 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, Christianson and others (2012) and Horgan 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 m³ 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 m³ 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.

Fig. 10. 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.