The supply of subglacial meltwater to the grounding line of the Siple Coast, West Antarctica

Abstract

Recent satellite studies have shown that active subglacial lakes exist under the Antarctic ice streams and persist almost to their grounding lines. When the lowest-lying lakes flood, the water crosses the grounding line and enters the sub-ice-shelf cavity. Modeling results suggest that this additional freshwater influx may significantly enhance melting at the ice-shelf base. We examine the spatial and temporal variability in subglacial water supply to the grounding lines of the Siple Coast ice streams, by combining estimates for lake volume change derived from Ice, Cloud and land Elevation Satellite (ICESat) data with a model for subglacial water transport. Our results suggest that subglacial outflow tends to concentrate towards six embayments in the Siple Coast grounding line. Although mean grounding line outflow is ~60m³ s^–1 for the entire Siple Coast, maximum local grounding line outflow may temporarily exceed 300 m³ s^–1 during the synchronized flooding of multiple lakes in a hydrologic basin. Variability in subglacial outflow due to subglacial lake drainage may account for a substantial portion of the observed variability in freshwater flux out of the Ross Ice Shelf cavity. The temporal variability in grounding line outflow results in a net reduction in long-term average melt rate, but temporary peak melting rates may exceed the long-term average by a factor of three.

1. Introduction

Most of Antarctica’s ice mass loss takes place from its floating ice shelves; mass loss occurs primarily through basal melting near the grounding line and ice fronts and iceberg calving at the ice front (Jacobs and others, 1992; Rignot and Jacobs, 2002). For an ice shelf to be stable (i.e. mass balance remains close to zero), these mass losses must be approximately equal to the mass gained through influx of grounded ice, precipitation and basal refreezing downstream of the grounding line melt zone. Ice shelves can provide a substantial amount of back-stress to their tributary glaciers and ice streams, slowing the flow of grounded ice (Thomas and Bentley, 1978; Schoof, 2007; Goldberg and others, 2009), and there is evidence that the loss of ice shelves leads to accelerated ice flow and sea-level rise (e.g. De Angelis and Skvarca, 2003; Rignot and others, 2004; Dupont and Alley, 2005; Jenkins and others, 2010; Shuman and others, 2011). In the case of a marine-based ice sheet such as the West Antarctic ice sheet (WAIS), melting and thinning of the ice shelf near the grounding line may lead to runaway thinning and retreat (Weertman, 1974). Since ice shelves play such a key role in the stability of the Antarctic ice sheet, understanding the processes that shape their evolution is critical to assessing the stability of the ice sheet over time.

While basal melting can occur anywhere under an ice shelf given appropriate oceanic conditions, the highest rates are usually found near the deep grounding lines. The principal mechanism causing basal melting near the grounding line is generally well understood. Ocean water with temperature T _O that is higher than the in situ freezing point T _f flows to the grounding line where it melts ice, creating a buoyant, cold and relatively fresh meltwater plume that carries mass away from the grounding zone (Lewis and Perkin, 1986; Jacobs and others, 1992; Jenkins and Bombosch, 1995). The basal melt rate w _b will be a function of thermal forcing (T _O – T _f) and a measure of the oceanic turbulent kinetic energy (see, e.g., Holland and Jenkins, 1999; Mueller and others, 2012). For regions where the principal source of oceanic turbulence is the plume flow itself, w _b is approximately proportional to (T _O – T _f)² (Lewis and Perkin, 1986).

The delivery of ocean water (and its associated heat content) to the grounding zone can be accomplished through the thermohaline circulation itself and by advection forced by processes external to the sub-ice-shelf cavity. The cavities under Antarctic ice shelves experience a wide range of ocean inflows, ranging from very warm Circumpolar Deep Water (CDW; θ≈1˚C) (Payne and others, 2004, 2007; Jenkins and others, 2010; Pritchard and others, 2012) to High Salinity Shelf Water (HSSW) formed by the cooling and rejection of brine during winter sea-ice formation offshore of the ice front. The temperature of HSSW is near T _f at the surface (~–1.9˚C for typical ocean salinities); however, it contains sufficient heat to melt ice near the grounding line because of the pressure depression of T _f (Lewis and Perkin, 1986). This ~3˚C range of thermal forcing around Antarctica is reflected in basal melt rates near grounding lines varying from <1 to >40ma^–1 (Rignot and Jacobs, 2002).

Observations in Alaskan and Greenland fjords have identified an additional cause of basal melting in the grounding zone: discharge of subglacial fresh water across the grounding line (Motyka and others, 2003, 2011). The flux of subglacial fresh water can substantially increase basal melt rates in the grounding zone as the vigorous circulation associated with the buoyant plume increases the rate at which ocean heat is made available to the ice base through turbulent entrainment. It has generally been assumed that this process is less important in Antarctica than in Greenland because of the much lower rates of meltwater production on the surface and at the base of the Antarctic ice sheet; however, there may be locations where subglacial discharge is significant, at least on short timescales (Jenkins, 2011).

While much of the subglacial freshwater discharge in Antarctica is likely to occur as a quasi-continuous flux from the water-saturated till at the base of glaciers and ice streams, we now know that episodic drainage of subglacial lake systems that have been found near the grounding line (Fricker and others, 2007) can lead to large rates of freshwater flux to the adjacent sub-ice-shelf cavity (Fricker and Scambos, 2009). In this paper, we consider the potential impact of specific lake drainage events on the circulation of the sub-ice-shelf cavity near the grounding line. We focus specifically on drainage from the hydrologically active subglacial lakes described by Fricker and others (2007, 2010), which discharge into the cavity under the Ross Ice Shelf (RIS) along the Siple Coast. We use a model of subglacial water flux combined with a time series of lake volume change for the period 2003–08 for the Siple Coast ice streams to determine both the spatial and temporal distribution of meltwater crossing the RIS grounding line, and examine the role of lake drainage events on the freshwater budget of the sub-ice-shelf cavity.

2. Study Region and Scope

We focus on the RIS, which is the largest Antarctic ice shelf by area (Bentley and others, 1979) and is relatively stable (Zwally and others, 2004; Pritchard and others, 2012). The typical draft at the RIS grounding line is 400–1000m (Neal, 1979; Griggs and Bamber, 2011) where T _f is ~–2.2 to ~–2.7˚C. Typical basal melt rates near the Siple Coast grounding line are ~0.05–1ma^–1 according to modeling (MacAyeal, 1984; Holland and others, 2003; Dinniman and others, 2007) and observational studies (Catania and others, 2006, 2010). This is one to two orders of magnitude less than melt rates that have been reported for grounding lines where the CDW intrusion process is occurring, such as the floating part of Pine Island Glacier (tens of meters per year; e.g. Rignot and Jacobs, 2002). Since basal melt rates near the RIS grounding line are relatively low, subglacial meltwater potentially could constitute a significant component of the RIS cavity’s freshwater budget and influence local ‘background’ circulation and melt patterns where the water emerges at the grounding line.

The RIS is supplied with ice from both the East Antarctic ice sheet (EAIS) via outlet glaciers crossing the Transantarctic Mountains and from the WAIS via the Siple Coast ice streams: Mercer (MIS), Whillans (WIS), Kamb (KIS), Bindschadler (BIS) and MacAyeal (MacIS) (Fig. 1). The Siple Coast ice streams move at speeds of up to 800 ma^–1 (Joughin and others, 2002; Rignot and others, 2011a), despite low driving stresses in part due to the presence of saturated deformable sediments (Blankenship and others, 1987; Engelhardt and others, 1990; Tulaczyk and others, 2000). While a combination of relatively thick ice and shear heating results in basal melt rates on the order of a few cma^–1 in the tributaries, freezing conditions dominate the basal environment of the downstream regions of the Siple Coast ice streams (Joughin and others, 2004a; Fig. 1). Therefore the downstream transport of this water is essential to maintaining lubrication in areas where the basal thermal regime would otherwise result in freezing and hardening of the subglacial sediments (Christoffersen and Tulaczyk, 2003; Parizek and others, 2003).

Fig. 1. Map of the Siple Coast and adjacent Ross Ice Shelf; inset shows the location in Antarctica. The grounding line (Bindschadler and others, 2011) is shown as a yellow line. Color scale indicates melt rates as modeled by Joughin and others (2004a) landward of the grounding line and ice-shelf-cavity thickness seaward (Le Brocq and others, 2010). Contours indicate lithostatic hydropotential (kPa). The six grounding line embayments described in the text are shown in dark grey. ICESat tracks are displayed as thin black lines. Lake locations and outlines obtained from Fricker and Scambos (2009), Smith and others (2009) and Fricker and others (2010). Background MODIS Mosaic of Antarctica imagery from Scambos and others (2007). Upper insert shows location of the RIS with respect to the Siple Coast ice streams (ice velocity from Rignot and others, 2011a). Lower insert shows location of larger map with respect to the continent. All maps are projected using polar stereographic projection with minimum distortion at 71˚ S.

The distribution of subglacial water under the Siple Coast ice streams varies significantly on timescales of weeks to centuries. On centennial timescales, piracy of water from upper KIS toward WIS is associated with the shutdown of lower KIS that occurred ~160 years ago (Anandakrishnan and Alley, 1997). On timescales of weeks to decades the periodic filling and draining of subglacial lakes introduces temporal variability to subglacial water flux crossing the grounding line, even in the absence of changes in the ice-sheet geometry (Gray and others, 2005; Fricker and others, 2007; Smith and others, 2009). Recent modeling work by Jenkins (2011) has quantified the relationship between subglacial outflow and near-grounding-line melt rates for the ice-shelf base and shown it to be significant for lake drainage events along the Siple Coast. Consequently the contribution of subglacial meltwater to the RIS sub-ice-shelf cavity may vary spatially and temporally.

Our study explores how the variability in subglacial water distribution over the period 2003–08 relates to subglacial outflow, with the following objectives: (1) to identify the spatial distribution of subglacial outflow; (2) to quantify the temporal variability in subglacial outflow in response to the inferred lake volume change; (3) to assess the impact of the calculated outflow on the melt rates of the RIS near the Siple Coast grounding line; (4) to compare the predicted freshwater production associated with subglacial outflow from the Siple Coast to the greater freshwater budget of the RIS.

3. Methods

3.1. Model overview and inputs

Our model is a simple steady-state water transport model, described by Carter and others (2011). It has four inputs: (1) basal melt rate; (2) a hydropotential surface; (3) lake locations and time series of their estimated volume changes; and (4) an accurate grounding line. The source datasets for these inputs are described below. The model directs subglacial meltwater formed at the ice/bed interface via basal melting down the gradient of the hydropotential. The model estimates the change in subglacial water distribution associated with the drainage of known subglacial lakes and predicts their refilling rates from the background melt rate (and from lakes upstream of them). To estimate fluxes across the grounding line required some subtle modifications of the original model, as described in Section 3.2.

Data used to generate the four model inputs originated from several independent geophysical techniques, as described below.

3.1.1. Basal melt rate

Our input term for subglacial water is basal melt rate from Joughin and others (2004a) (Fig. 1); basal drag and associated shear heating in the ice column is calculated from ice-sheet geometry and observed surface ice velocity (Joughin and others, 2002, 2004b). The inferred value for shear heating, along with an assumed spatially homogeneous value of 70 mWm^–2 for geothermal flux, surface accumulation and calculated velocity distribution with depth, were then fed into a temperature model from which excess heat was converted into basal melt. The estimated basin-wide error was ±10% (Joughin and others, 2004a).

Since Joughin and others’ (2004a) study, there have been improved estimates for geothermal heat flow (e.g. Shapiro and Ritzwoller, 2004; Fox-Maule and others, 2005), ice surface geometry (Haran and others, 2006), ice thickness (Le Brocq and others, 2010; Carter and others, 2011) and, most recently, ice velocity (Rignot and others, 2011a). Furthermore, more recent work has shown that heat produced along the shear margins (neglected in Joughin and others, 2004a) may enhance basal melt in some areas (Beem and others, 2010); however, it is unclear what the significance of this is for the entire water budget. At present, we consider that a published dataset that is spatially continuous across the model domain and has been shown valid in part of our domain (Carter and others, 2011) is preferable to developing a new and untested geothermal heat flow model for this study.

3.1.2. Hydropotential

We estimated hydropotential, θ _h, from the sum of the overburden pressure and ice-base elevation. In this calculation we assumed that water pressure is equal to the overburden pressure, ρ _i gh _i, and used

(1)

where z _srf is the ice surface elevation, h _i is the ice thickness (equal to z _srf – z _bed, where z _bed is the elevation of the ice base), and ρ _i and ρ _w are the densities of ice and fresh water, respectively (Paterson, 1994). From Eqn (1) we see that surface gradients are about ten times more important than basal gradients for determining the regional direction of water flow (Shreve, 1972). Estimating the hydropotential requires a gridded dataset for ice surface elevation and ice thickness.

For z _srf we used a digital elevation model (DEM) derived from Ice, Cloud and land Elevation Satellite (ICESat) laser altimetry, acquired between 2003 and 2004 and enhanced with Moderate Resolution Imaging Spectroradiometer (MODIS) photoclinometry (Haran and others, 2006); this DEM has a resolution of 250m and each gridcell is subject to an error of 1–2 m. For h _i we used radio-echo sounding (RES) data acquired in 1971–2001, and seismic measurements from 1957–74 (Bentley and Ostenso, 1961; Bentley and Chang, 1971; Retzlaff and others, 1993; Blankenship and others, 2001; Lythe and others, 2001). Where there were no ice thickness estimates, we interpolated values using an algorithm described by Carter and others (2011). The ice thickness estimates have an error of up to 50 m due to factors such as: uncertainty in the radio/seismic wave propagation velocity in ice; navigation error; and off-nadir reflections (Lythe and others, 2001). Additional errors may have resulted from change in the ice thickness since the time of measurement of up to ±7m in areas where the data are either relatively old (i.e. seismic traverses around the tributaries of MacIS) or where rates of surface change are high (i.e. KIS) (Pritchard and others, 2009). From Eqn (1), an error of 50m in the thickness measurement, 7m in ice thickness change and 2 m in ice surface elevation equates to a total error of 46 kPa in the hydropotential. Translating this uncertainty in ice-sheet geometry to uncertainty in subglacial water flux at a given location is somewhat more complicated, as Wright and others (2008) have shown that small changes in the hydropotential surface in key locations may lead to large changes in water distribution downstream.

3.1.3. Lake locations and volume time series

We obtained lake locations and volume time-series information from three published datasets based on ICESat repeat track surface altimetry (2003–08); two of the datasets also used MODIS imagery. ICESat operated in campaign mode so that the time interval between each point in the volume time series was 3–4 months. For WIS, MIS and MacIS, we used lake outlines and volumes from Fricker and Scambos (2009) and Fricker and others (2010) because these incorporated either MODIS image differencing or MODIS imagery to improve the estimates of lake outline and area, and therefore of volume changes. Both of these studies calculated an error of ±20% for volume change undergone by each lake. Elsewhere, we used the data from Smith and others’ (2009) inventory, whose volume change estimates were based solely on ICESat data, and therefore limited spatially by the sampling of the ICESat altimeter. Consequently this study reported larger error estimates of ±50% for each lake.

The analyses for all three of the volume change datasets assumed a one-to-one ratio between observed surface displacement and lake volume change. Modeling work by Sergienko and others (2007) suggested this assumption may not be valid. They showed how surface uplift and subsidence could continue to occur long after filling and draining had ceased. They further suggested that localized surface uplift and drawdown could result even from processes unrelated to lake filling and draining. However, they considered primarily small (<10km diameter) subglacial lakes located in steep-sided bedrock depressions, which were surrounded by areas of high basal traction. In contrast, most of the lakes in our study are located in areas of low basal shear stress and occupy relatively shallow depressions in the basal topography. Consequently, Fricker and Scambos (2009) found secondary flow effects of lake drainage to be minimal. Furthermore surface uplift and drawdown due to ice dynamics would produce features of a different shape from those that resulted from lake drainage, with uplift in close proximity to drawdown. Smith and others (2009) attempted to eliminate any feature in which the shape of the surface change appeared to result from changes in basal traction rather than subglacial water volume change.

3.1.4. Grounding line position

For the location of the grounding line, we used the landward limit of tidal flexure (point ‘F’) compiled by the Antarctic Surface Accumulation and Ice Discharge (ASAID) project (Bindschadler and others, 2011). This grounding line was based on discrete estimates of F derived from repeat-track analysis of ICESat data by Brunt and others (2010), which were interpolated using Landsat imagery to form a continuous line. Although the recently published grounding line based on interferometric synthetic aperture radar (InSAR) data described by Rignot and others (2011b) is more precise and accurate, it is discontinuous in some of our study area; where it is present it does not depart significantly from the ASAID grounding line.

3.2. Model description

Our steady-state water model is initialized with the input basal melt rate (Section 3.1.1) and hydropotential (Section 3.1.2) fields. The model directs basal meltwater down the hydropotential surface using the steady-state 8 Directional (D8) formula (Quinn and others, 1998)

(2)

where Q _in is water flux into a cell, m _ is the basal melt rate (negative if water is freezing to the base) and Δx and Δy are the cell’s horizontal dimensions. The water flux out of a cell, Q _out, is apportioned among all ‘downstream’ cells:

where Q_i is the flux into one of up to eight adjacent cells that are down-gradient, k is the number of adjacent cells with lower hydropotential, (dθh/ds)i is the hydropotential gradient to the cell receiving Q_i , and is the sum of the hydropotential gradients to all adjacent downstream cells. This formulation distributes proportionally more water flow toward adjacent cells with steeper downward gradients.

The input lake location and volume time series (Section 3.1.3) are used to inform the model where the known subglacial lakes exist. For cells lying within known lakes, the model applies a time-dependent parameterization, depending on whether the lake is interpreted to be filling or draining. For time intervals during which a lake was filling, the lake cells are treated as sinks, setting Q _out to zero. For intervals during which a lake is draining, the observed volume loss is added to the melt rate for the cells within that lake, and water from upstream is allowed to pass through them, using Eqns (2) and (3). Effectively, the observations of lake drainage force the model, while observations of lake filling are used to validate the model’s treatment of flow from upstream.

We made two subtle modifications to the model as described by Carter and others (2011) for this study, to enable us to estimate subglacial water flux across the grounding line: (1) we allowed selected lakes to ‘overflow’ during their filling cycles (see Section 3.3.3); and (2) we divided the RIS sub-ice-shelf cavity, which had been previously treated as a single ‘lake’, into several separate ‘lakes’ to keep track of flux across the grounding line (see Section 3.2.4).

3.3. Model implementation

The implementation of the model was an iterative process. We started with the hydropotential surface (Section 3.1.2) and the time series of lake volume changes (Section 3.1.3) and ran the model to see if the water budget balanced. In the case of MacIS, it did balance (Carter and others, 2011); this was not the case, however, for every drainage system. If the water budget did not balance, we returned to the input data to establish whether there were some localized regions where a small adjustment in model input would improve the water budget; for example, if there was another lake nearby we examined its behavior to see if the water budget balanced more closely when both lakes were treated as a larger single lake; this was the case for lakes L7 and L8 in WIS/MIS (Fricker and Scambos, 2009), which we combined into L78. If there was a significant amount of water flowing in a gridcell adjacent to the lake, we examined the hydropotential grid to see whether a small adjustment in the hydropotential could direct sufficient water towards the lake.

3.3.1. Stream etching

To improve water routing in areas where water discharge was likely to be high, and minimize errors in modeled water distribution associated with uncertainties in the final hydropotential surface, we applied a ‘stream-etching’ technique to determine the hydropotential surface for continuous local minima, or ‘stream valleys’. From the initial map of hydropotential, we began by identifying several major drainage paths (Fig. 1).We then used an algorithm (described by Carter and others, 2011) which sets the hydropotential of any gridcell through which the stream valley passed to that of the local minimum in the hydropotential, as identified in the original RES ice thickness profiles. Additionally, where two separate drainage paths came close to one another we raised the hydropotential of cells between them to keep flows separated; this was necessary for MacIS (Fig. 2).

Fig. 2. Steady water flux from a control run in which no lakes were incorporated. Lake colors indicate comparison between modeled and observed filling rates. Lakes that exhibited close agreement between modeled and observed filling rates are colored cyan. For other lakes, there are substantial differences between modeled and observed lake-volume change. Lakes whose filling rates, derived from observations, were different from the model are colored yellow (‘underfilling’) or magenta (‘overfilling’). Lakes colored white (‘through-flowing’) were eliminated from the study. Lakes colored green (‘leaking’) were allowed to release some water downstream while still filling. Triangles denote locations where the hydropotential was either raised (upward-pointing) or lowered (downward-pointing) as part of the tuning process described in Section 3.3.2. Lake name abbreviations are: SLM = Lake Mercer, SLW = Lake Whillans, SLE = Lake Engelhardt, SLC = Lake Conway, USLC = Upper Subglacial Lake Conway, L78 = lakes 7 and 8 from Fricker and Scambos (2009); M1–M6 = Mac 1–Mac 6 from Carter and others (2011); B1–B6 = Bindschadler 1–Bindschadler 6, KT = Kamb Trunk, R1 = Raymond, K1–K12 = Kamb 1–Kamb 12, and W6–W8 = Whillans 6–Whillans 8 from Smith and others (2009).

3.3.2. Model tuning and validation

We validated and tuned the model by comparing the modeled flux of water through the subglacial lake system with observed estimates of subglacial lake volume change derived from satellite observations of the surface. We followed the approach described by Carter and others (2011), in which lake volume changes were shown to account for all the meltwater produced and transported in the MacIS catchment. By validating our modeled with observed filling rates from ICESat, we confirmed that our hydropotential model, our melt rate distribution, lake volume changes and water transport model were all self-consistent for this study; this gave us confidence in the accuracy of our estimates for water flux across the grounding line. More details of the model performance are provided in Section 4.1.

When initially comparing the model with the observations, we found that the model did not supply enough water to many of the lakes to match the observed filling rates. In many cases these lakes were found to be immediately adjacent to, but not intersected by, a stream valley. For each of these cases we experimented with raising or lowering the hydropotential in the immediate vicinity such that more water would be directed toward the lake (Fig. 2). For any given basin we began this process with the most upstream lake and then worked downstream.

A noteworthy case where we had to adjust the hydropotential was at Lake Whillans (SLW; orange triangle in Fig. 2). While initially we assumed that the hydropotential surface remained static over the study period, we found that in order to reproduce the observed filling rate for SLW it was necessary to increase the hydropotential surface by 14 kPa between late 2005 and late 2006.

We considered a lake to be in balance if, after adjustment to the hydropotential surface, the modeled and observed volume changes agreed to within a factor of 1.5, which is the error stated in Section 3.1.3.

3.3.3. Categorization of out-of-balance lakes

When no reasonable adjustment to the hydropotential surface upstream could reproduce the observed filling rate, the lake was considered to be out of balance. These lakes fell into one of four categories: ‘overfilling’, ‘underfilling’, ‘through-flowing’ and ‘leaking’. Overfilling and underfilling referred to lakes for which the model respectively either overestimated or underestimated the observed filling rate, by a factor of two or greater (Fig. 2). The misfit represented by these lakes became another source of uncertainty for our model of flux across the grounding line. This was usually because both the modeled and observed volume changes experienced by such lakes constituted a small portion of the total water budget (the budgets were all <0.2km³ and most were <0.05km³).

For lakes classified as ‘through-flowing’ or ‘leaking’, additional parameterization was necessary. Through-flowing lakes were typically observed to be filling at a rate less than one-tenth of the modeled rate. For such lakes it was often easier to balance the budget of lakes downstream by eliminating the lake completely. Similarly for leaking lakes it was difficult to balance the water budget of lakes downstream when the lake was treated as a perfect sink while filling. Unlike through-flowing lakes the observed volume change of leaking lakes comprised a significant fraction of the modeled influx; therefore, it was important that the volume change undergone by these lakes was accounted for in the total water budget. For these lakes, any modeled inflow in excess of the observed inflow was directed to the downstream side of the lake, and added to the melt rate of the nearest gridcell downstream of the lake (Fig. 2).

3.3.4. Estimation of flux across grounding line

After tuning, we ran the model to quantify the total flux of subglacial meltwater to the RIS cavity and the variability in that flux over time due to subglacial lake filling and draining. The RIS grounding line location (Section 3.2.4) was used to inform the model where water exits the subglacial system and enters the cavity. We identified six embayments in the grounding line, defined by distinctive concave sections, each supplied by a distinct water system (dark grey sections in Fig. 1 and 2). Each of these embayments was treated as a separate lake and was considered to be filling throughout the model run. In order to further assess the origin of water entering the sub-ice-shelf cavity, we partitioned flux across the grounding line into two categories: water derived from subglacial lake drainage, and water that was part of the background subglacial melt. While the latter may have passed through a lake, we assumed it was never stored in one.

3.3.5. Control run

We also performed one run of the model, hereafter referred to as the ‘control run’, in which there were no subglacial lakes, except for the embayments described above to simulate water flux in the absence of subglacial lake activity. The purpose of this run was (1) to provide a background estimate for subglacial outflow against which we could compare our other results and (2) to provide a diagnostic subglacial water distribution that would highlight major flow paths (see Section 3.3.2).

4. Model Performance

The control run revealed that the steady-state water distribution agreed well with the locations of known subglacial lakes. The water budgets for the downstream portions of MacIS, WIS and MIS all balanced. Most subglacial water flow was concentrated into pathways of elevated hydrologic flux that were 10–15km wide. All known subglacial lakes were located less than one 5 km gridcell from these pathways (Fig. 2).

When we forced the water source term with lake drainage observations and treated filling lakes as sinks (Section 3.2), the subglacial water distribution varied considerably over time throughout the model domain, in particular in the trunks of WIS and MacIS (Fig. 3). In this section, we evaluate the performance of our model with respect to reproducing observed lake volume changes (Section 4.1), and discuss the causes of differences between modeled and observed filling rates (Section 4.2). We then discuss the spatial and temporal variability in water distribution upstream of the grounding line (Section 4.3).

Fig. 3. Modeled water flux for time intervals corresponding to the temporal sampling of the ICESat campaigns relative to steady-state water flux predicted by the ‘control’ run shown in Figure 2.

4.1. Comparison between modeled and observed lake filling rates

For the 17 largest subglacial lakes in the model domain (including Mac 1, 2, 3, 4 and 5 (M1–M5), Bindschadler 1, 2 and 3 (B1–B3), Raymond (R), Kamb 4 and 5 (K4, K5), Whillans 6 and 7 (W6, W7), Upper Subglacial Lake Conway (USLC), SLW and Lake Engelhardt (SLE)) we found that the modeled and observed filling rates agreed to within a factor of 1.5; we considered these lakes to be in balance. In particular, for MacIS the combined modeled filling rates of all the lakes matched the sum of rate of meltwater production and net volume loss of lakes upstream in this basin. This confirmed that the input melt rates, hydropotential surface, water-routing method and observed volume changes were consistent and that the flux of water downstream of these lakes was realistic.

For 22 smaller lakes in the model domain, the model could not reproduce the observed filling rates; we classified these lakes as out of balance (Section 3.3.3). All of the overfilling and underfilling lakes were found in the upper parts of the WIS, KIS, BIS and MacIS catchments. While we only identified two through-flowing lakes (Kamb 2 (K2) and lake 7/lake 8 (L78)), we identified several overfilling and underfilling lakes in the upper parts of each hydrologic catchment (Fig. 2). There was some clustering within these two classifications. A number of the overfilling lakes were located along the same flow path that connects Kamb 10 (K10) with Kamb 1 (K1). Although the number of lakes that were out of balance was greater than the number of lakes in balance, the volumetric contribution of out-of-balance lakes accounted for only 5–10% of the total water budget and was never more than 30% of the observed volume change during any time interval. Most drainage basins contained at least one overfilling and one underfilling lake such that for the basin-wide water budget the errors were less than the sum of the individual errors, and the out-of-balance pair did not make a significant contribution to the errors in the water budgets of lakes further downstream. Consequently it was possible to balance the water budget of lakes downstream despite uncertainties in the filling and drainage rates of lakes that were out of balance.

The two primary leaking lakes (Conway (SLC) and Mercer (SLM)) were both located on the same flow path along the WIS/MIS trunk. We obtained better agreement between the modeled and observed filling rate of SLM immediately downstream when we allowed SLC to leak while filling, rather than treating it as a perfect sink. Likewise, when we allowed SLM to leak while filling, modeled water flux through L78 became more consistent with observations (Section 5.3) than when it was treated as a perfect sink. Overall the water budget of lower WIS balanced. There was one time interval (March–June 2006; Fig. 3h) during which SLC, SLW and SLE were all filling simultaneously; this made it possible to account for all the meltwater in the basin through their combined volume changes during that interval, without requiring any overflow from SLC. The sum of their filling rates was 33 ±3m³ s^–1 (modeled) and 38±4m³ s^–1 (observed) and balanced the sum of the modeled meltwater production of the greater basin and observed net volume loss of all lakes upstream (32±3m³ s^–1).

The agreement between our model and observations was greatly improved through several key adjustments to the hydropotential, in particular in the most-downstream sections of WIS and MacIS (Fig. 2). In both these locations, however, the hydropotential adjustment was primarily to keep adjacent flow paths separated from one another. Surface slopes, and therefore hydrologic gradients, are quite low in both these locations and there are many places where a small localized change in the hydropotential on the order of a few meters might reroute substantial amounts of water. One region where tuning of the hydropotential was required was around Whillans 8, and this was probably due to the lack of ice thickness information in this region, resulting in a higher uncertainty in the input hydropotential field. Two other places where tuning was necessary were near lakes Raymond (R1) and Bindschadler 3 (B3). Unlike other lakes in this study, these small lakes were both located one gridcell from flow paths along which steady discharge was high rather than directly on it, and thus tuning helped direct more water to them.

4.2. Causes of discrepancy between modeled and observed filling rates

Although we were able to account for most of the meltwater budget with lake volume change, many of the lakes in this study were out of balance, i.e. the modeled and observed filling rates were different. There are several possible explanations for this misfit, that vary depending on the location of the lake and nature of the misfit, which appears to lie both in the nature of the data and assumptions that went into our model:

1. 1. For most of the overfilling and underfilling lakes, the observed volume change over the time frame of the study was <0.2 km³ (Smith and others, 2009), and their shorelines were determined entirely from ICESat data without any MODIS imagery. Therefore, volume change estimates for such lakes had large uncertainties.

2. 2. We assumed a one-to-one ratio between vertical surface motion over a subglacial lake and volume change. Although this assumption appears to have been valid for the larger lakes in the study area, where rates of filling and draining were relatively high (tens of m³ s^–1) and basal surrounding topography is gradually sloped, it may not have been valid for some of the smaller lakes such as Whillans 8 (W8) and Kamb 12 (K12), which occupied relatively steep-sided bedrock depressions and experienced relatively low rates of water influx (<0.1m³ s^–1). For underfilling lakes in particular, some of the observed surface uplift may have resulted from ice flowing into the topographic depression containing the lake as predicted by previous modeling work (Sergienko and others, 2007; Pattyn, 2008). In particular, the basin containing Kamb Trunk (KT) may have been steadily filling with ice while periodically draining water throughout the study interval.

3. 3. Both of the through-flowing lakes in this study (K2 and L78) may not actually be storage features. These two lakes lie directly downstream of some of the largest lakes in their respective clusters (K2 is directly downstream of Kamb 5 (K5) and L78 is directly downstream of SLM). The modeled quantity of water flowing into either of these lakes showed better agreement with the total observed ‘volume’ than it did with the observed rate of volume change. Although they do not appear to have stored significant quantities, the fact that surface change over these features correlated so well with discharge makes them an interesting target for ongoing modeling and observations of the subglacial water system in this area.

4. 4. We assumed that a lake acted as a perfect sink while its volume was increasing. If, however, the rate of outflow was less than the rate of inflow, the lake could actually still increase in volume while simultaneously discharging a significant quantity of water downstream; this is the case for leaking lakes. The two largest leaking lakes in our study area (SLC and SLM) are both in the same interconnected system on WIS. They occupy a flow path along which the hydropotential gradient is low (Fricker and Scambos, 2009) and the influx of water from points upstream was high during this period. If either SLC or SLM were to drain via a semicircular channel as shown by Carter and others (2009), then the low hydropotential gradient would have limited the amount of turbulent viscous dissipation available to melt the channel walls (Nye, 1976; Evatt and others, 2006). During the early part of the drainage event, the discharge would have been less than the inflow; consequently water levels would have continued to rise.

4.3. Variability in water distribution upstream of grounding line

Variability in water distribution was most pronounced in locations where lakes are clustered (e.g. lower MacIS, upper KIS and lower WIS/MIS; Fig. 3). Outside these lake clusters, the greatest variability in water distribution with time occurred near lakes located in regions where basal freezing is believed to be widespread (e.g. KT in the trunk of KIS; Fig. 3). In contrast, on BIS where basal melting is widespread and subglacial lakes are absent in the downstream part of the trunk, our model predicted relatively constant water distribution over time despite the filling and drainage of several lakes in the upper part of the catchment.

Our modeled time-varying water distribution also predicted several ‘cascading floods’ where lakes upstream transferred water down the hydrologic gradient to lower lakes, which subsequently drained into lakes farther downstream. This behavior was difficult to verify with the observations, because the ~3–4 month ICESat sampling period precluded us from differentiating between two lakes draining simultaneously and two draining sequentially. We interpreted all events as sequential, with drainage initiating first at the most upstream lake in the chain. We inferred simultaneous drainage of four lakes in the MacIS catchment when the drainage of Mac 4 (M4) and Mac 5 (M5) (Fig. 3b) began during the same time interval as the drainage of Mac 3 (M3) (Fig. 3b) and Mac 1 (M1) (Fig. 3b–e) (Fricker and others, 2010; Carter and others, 2011). While the largest of these drainage events was only 100±20m³ s^–1, the combined discharge of several lakes draining simultaneously led to a flood exceeding 300 m³ s^–1. Meanwhile in the KIS–WIS system, we inferred that the water from two separate cascading floods traveled >300 km, with drainage of the uppermost lake preceding drainage of the most downstream lake by many months. In one case, drainage from Whillans 6 (W6) (Fig. 3c) was followed by the drainage of USLC (Fig. 3f), which was followed by the drainage of SLC and SLM (Fricker and Scambos, 2009) (Fig. 3i). Drainage of K5 and Whillans 7 (W7) (Fig. 3f and g) preceded drainage of SLW (Fig. 3i).

5. Subglacial Flux Across the Grounding Line

We now present the derived spatial and temporal distribution of the subglacial flux across the grounding line. Outflow occurred at only a few distinct locations along the grounding line. Where these locations are downstream of large subglacial lakes, the variability could be substantial (Section 5.1). We discuss how the spatial concentration and temporal variability of subglacial outflow can lead to temporary localized elevation of the ice-shelf basal melt rates above the background (Section 5.2). We explore how the enhanced melting may affect the evolution of the grounding line over time (Sections 5.2.3 and 5.2.4). In Section 5.3, we characterize the contribution of subglacial outflow to the freshwater budget of the RIS cavity and the embayments into which it first emerges. These results and inferences have implications for both the limitations of this study and ideas to further advance our understanding of the interactions between the subglacial water system and the ice-shelf cavity.

5.1. Emergence of subglacial water across the grounding line

Our model predicted that subglacial water crossed the grounding line at a few distinct locations along the Siple Coast. At these locations, subglacial water was concentrated into pathways that were a few gridcells (5–15 km) wide. The six locations where each of these pathways crossed the grounding line each coincided with prominent embayments (Fig. 1).

Our model showed that the flux of subglacial water over time among the six embayments varied considerably. The variability over time was strongly correlated with the presence of large subglacial lakes near the grounding line, as was the case for the MacIS (Fig. 4f), SLE (Fig. 4c) and SLM embayments (Fig. 4a). At times, flow rates exceeded 300 m³ s^–1, a factor of four higher than the highest rate of volume loss for any individual lake. In between lake drainage events, subglacial flux was often <5m³ s^–1. Similarly, the subglacial freshwater contribution to the SLW and KIS embayments (Fig. 4b and d) was primarily controlled by the drainage of lakes within these basins. The maximum outflow into these embayments was over an order of magnitude less than the MacIS, SLE or SLM embayments, and not significant when compared to the fresh water produced by background melting of the ice-shelf base. The BIS embayment (Fig. 4e) was unique among the six embayments in that, despite high mean subglacial outflow (19±4m³ s^–1), subglacial lake drainage did not contribute significantly to this total. Although some lakes in the BIS catchment were filling throughout the study period, their filling rates were on the order of 0.05 km³ a^–1 and made only a small impact on the total water budget. The absence of lakes in the trunk of BIS resulted in a relatively constant supply of water across the grounding line.

Fig. 4. Comparison of subglacial water flux across the grounding line versus order-of-magnitude estimate for ice-shelf meltwater production for each of the six embayments shown in Figure 1: (a) SLM drainage; (b) SLW drainage; (c) SLE drainage; (d) KIS; (e) BIS; and (f) MacIS. (g) Comparison of subglacial water transport across the Siple Coast grounding line against total ice-shelf meltwater production for the RIS (Loose and others, 2009).

Our control run estimated the total subglacial outflow for the Siple Coast to be 64 ±6 m³ s^–1 (Fig. 4g). The time-varying model showed that outflow varied by an order of magnitude from 26 ±4 m³ s^–1 (March–June 2006) to 519 ±53 m³ s^–1 (March–June 2004). The mean outflow over the study period was greater than that estimated by the control run, indicating that the total subglacial lake volume lost exceeded the total volume gained. As previously stated, however, this outflow was focused in a few distinct locations, and, during much of our study period, outflow at many of these locations was <10m³ s^–1.

5.2. Subglacial outflow and ice-shelf melting near the grounding line

Subglacial water not only influences the ice-shelf cavity through direct contribution to the freshwater budget in this environment, but also indirectly, through enhancing basal melting at the grounding line. The additional melting of the ice-shelf base associated with subglacial outflow may also constitute a significant part of the freshwater budget for the RIS cavity. We focus on the relationship between subglacial outflow and ice-shelf basal melting, comparing the contribution of both processes combined to the total freshwater budget of the ice-shelf cavity.

A recent modeling study by Jenkins (2011) predicted that for a range of environments under the RIS a tenfold increase in freshwater outflow would increase basal melt rates by a factor of two. We used this work to quantify the additional sub-ice-shelf melting resulting from subglacial water injection at the grounding line. We then compared our estimates for freshwater flux against previously published estimates for the background ice-shelf basal melt rate near the grounding line. Although we did not have precise estimates for the total melting of the ice-shelf base in each embayment, modeling results by Holland and others (2003) and Dinniman and others (2007) have suggested that a reasonable order-of-magnitude estimate for sub-ice-shelf melt rate in the ice shelf cavity is 0.33 ma^–1; this estimate has been validated in part by layer analysis from radar sounding data for this area by Catania and others (2010).

5.2.1. Spatial distribution of subglacial outflow

Our model showed that subglacial water flux across the grounding line is concentrated into channels ~10–15 km wide in a few limited locations. This result adds some validity to work that has previously neglected its influence on ocean circulation on a large scale (e.g. Holland, 2008). Additionally, it supports the outflow zone width estimate of 10km assumed by Jenkins (2011). The concentration of subglacial water flow into channels results in a locally enhanced basal melting where the channel intersects the grounding line, but a net reduction in melting elsewhere. For BIS the effect of concentrating the 19±2m³ s^–1 steady-state outflow along a 10 km section of the grounding line, rather than distributing it evenly over the entire 275 km, leads to a 90% decrease in the volume of ice shelf melted by subglacial outflow. In the vicinity of the outflow location, however, there is a threefold increase in the basal melt rate.

5.2.2. Effect of temporal variation in outflow rates on basal melt rate near the grounding line

Using the same parameters for ice-shelf thickness, basal slope, ocean salinity and temperature as Jenkins (2011) used for WIS, while varying subglacial outflow, we estimated the melt rates one would get from outflow via our control run and those associated with variable flow over time. Whereas Jenkins (2011) had values only for background flow and a rate of lake volume loss, we were able to incorporate estimates for the total outflow. For our maximum SLE outflow rate of 168±24m³ s^–1, melt rates were ~3ma^–1. Minimum melt rates were closer to 1ma^–1. When we used the outflow estimate from our control run in which lakes were excluded, we obtained a melt rate of 1.5 ma^–1. Although the mean melt rate (1.6ma^–1) for the study period (2003–08) was close to the mean melt rate predicted by the control run (1.5ma^–1), the mean outflow (62±10m³ s^–1) during this time period was over ten times the outflow from the control run (5±0.5m³ s^–1). This insignificant increase in average melt rate for such a substantial increase in average outflow relative to the control suggests that varying the outflow over time actually reduces the net melt in the long term. In other basins, there were similar reductions in net melting relative to the control run. This suggests that the net effect of lake filling and draining is a short-term increase in melt rate while lake drainage is occurring, but over the long term a net reduction in melting of the ice-shelf base.

5.2.3. Potential feedbacks between subglacial hydrology and grounding line evolution

Understanding the influence of our modeled outflow distribution on the evolution of the grounding line depends on knowing the relative importance of peak versus mean ice-shelf basal melt rates for grounding line stability. If maximum melt rates are more important than mean melt rates to grounding line retreat, then subglacial lake drainage may allow an unpinning of the ice sheet leading to an accelerated retreat within the embayment (Goldberg and others, 2009). Otherwise we would expect lower rates of long-term retreat for the grounding line downstream of an active subglacial lake system relative to outflow points where the discharge is stable.

Given the localized increases in basal melt rates at subglacial outflow points, the presence of pronounced embayments along the coastline around all the major outflow locations may result from a positive feedback. If enhanced grounding line retreat were focused on small embayments where subglacial outflow is high, then we would expect the drawdown at that location to propagate into surface lowering further into the interior. Even in low-sloping areas where drawdown associated with grounding line retreat would be small, this would still have the effect of drawing water in from a larger area, leading to a slight positive feedback between grounding line retreat and routing of subglacial water. At present, however, it is not clear how this feedback is initiated or what limits it.

Over time intervals of decades to centuries it has been shown that the ice-flow (Catania and others, 2012) patterns and subglacial hydrology (Anandakrishnan and Alley, 1997; Wright and others, 2008) have changed substantially for the Siple Coast. Changes to the subglacial hydrologic patterns in the interior may have led to substantial changes in the spatial distribution of subglacial outflow over time. When KIS was active and its grounding line was supplied presumably by a much larger water system, melt rates of the adjacent ice shelf were likely much higher. These higher melt rates in addition to changes in ice flow may explain some of the inferred changes of the grounding line position over the last few centuries (e.g. Catania and others 2006, 2012).

5.2.4. Co-evolution of lakes and ice shelves

Many of the lakes identified beneath the Siple Coast ice streams appear to have resulted from the ice moving rapidly over a bed of varying basal traction so as to create self-sustaining depressions in the ice surface that resulted in hydropotential depressions at the ice base, which then filled with water (Fricker and others, 2010; Sergienko and Hulbe, 2011). Although some of these depressions in the ice surface may also have resulted from corresponding depressions in the subglacial topography, given the low basal relief of the Siple Coast, the magnitude of surface depressions associated with along-flow variations in basal traction are sufficient to create localized depressions in the hydropotential in which water collected. This contrasts with lake clusters near the ice divide that formed by simple filling of hollows in the bedrock topography, as observed near ice divides (Tabacco and others, 2006). Given the dependence of lake locations on longitudinal variations in the stress field, we suspect that the additional back-stress provided by the ice shelf may also play a role in the formation of subglacial lakes near the grounding line. Once water begins collecting, the reduction in melting at the grounding line could then lead to further thickening and impoundment of greater quantities of water between floods. Meanwhile loss of an ice shelf may lead to surface steepening and subglacial lake drainage, as was observed by Scambos and others (2011) when a lake beneath Crane Glacier drained following the collapse of the Larsen B ice shelf. Continued monitoring of grounding zones downstream of large subglacial lakes (e.g. SLE), along with detailed modeling of ice flow in this environment, is necessary before this hypothesis can be confirmed.

5.3. Impact of subglacial lake drainage on the freshwater budget of the RIS

The total outflow predicted by the control run (no lakes) of 64±6m³ s^–1 is an order of magnitude less than the total freshwater budget for the RIS (33–50 km³ a^–1, or 1045– 1584m³ s^–1) as estimated by Loose and others (2009) using chemical tracers. This total value includes melting near the ice front, which is significant for the RIS (Horgan and others, 2011). When we introduced subglacial lake discharge we found that subglacial outflow to the RIS varied between 26 and 500 m³ s^–1, with the maximum outflow corresponding to a time when lakes in the MacIS water system and the SLE water system were flooding simultaneously between March and June 2004 (Figs 3b and 4g). If we assume the enhanced melting took place over a zone 10km wide extending a little less than 7 km from the grounding line (these values based on Jenkins, 2011) for each of the six outflow points, we could expect an additional 10–20m³ s^–1 of additional ice-shelf melt. These estimates suggest that variability in subglacial outflow could account for a significant portion of the variability of the freshwater budget of RIS over time; however, variability in additional ice-shelf melting, at least in the ~7 km zone in which Jenkins’ (2011) reduced model is applicable, does not contribute significantly to this.

Although subglacial outflow and enhanced melting do not comprise a significant portion of the total RIS freshwater budget, the subglacial contribution may be more significant within the individual embayments. Water discharges into all six embayments (Fig. 2), where the water column thickness is often <50m (Greischar and Bentley, 1980). A comparison between subglacial outflow and estimated rates of ice-shelf melt suggests three different environments in relation to the influence of subglacial water along the grounding zone.

In the first environment, typified by MIS (Fig. 4a), SLE (Fig. 4c) and MacIS (Fig. 4f), floods from subglacial lakes dominate the freshwater budget of the system and are comparable in magnitude to the basal melt rates within the cavity, at times dominating the freshwater budget of the region. Between flooding events, the amount of subglacial water supplied to the cavity is relatively low or even absent. All three of these embayments are fed by catchments with large subglacial lakes near the grounding line. In these embayments, the behavior of the subglacial lake may be an important control on the supply of buoyant fresh water and, likely, the circulation within the cavity.

In a second environment represented by the outflow of water from SLW (Fig. 4b) and KIS (Fig. 4d), very little subglacial water would emerge were it not for the periodic draining of subglacial lakes. While subglacial water may at times be an important source of fresh water, melting of the ice-shelf base still appears to be the dominant process for freshwater generation. Both of these embayments adjoin a portion of the ice sheet where the thermal balance causes widespread basal freezing (Fig. 1), so it is only during flood events that subglacial water reaches the grounding line.

The third environment, at the grounding line of BIS (Fig. 4e), is unique among the Siple Coast ice streams in that, despite high mean subglacial outflow (19±2m³ s^–1), subglacial lake drainage did not contribute greatly to this value. Although some lakes in the catchment were filling throughout the study period, their filling rates were on the order of 0.05 km³ a^–1 and made only a small impact on the total water budget. The absence of lakes in the trunk of BIS results in a relatively constant supply of water across the grounding line.

6. Conclusions

We have used a simple model for subglacial water transport combined with information on subglacial lake volume changes inferred from ICESat to quantify how the distribution of subglacial outflow across the grounding line of the Siple Coast varies over the ICESat period (2003–08) in response to the filling and draining of subglacial lakes in this area.

The majority of subglacial water appears to be channeled into waterways ~10–15 km wide, and emerges across the grounding line at a few limited locations corresponding to pronounced embayments along the coastline. Within the embayments, subglacial outflow varied by up to two orders of magnitude over the time frame of our study. This variability was most pronounced at locations immediately downstream of large subglacial lakes. Outflow variability was further increased when the water from cascading floods involving multiple lakes crossed the grounding line, which led to peak subglacial outflow rates higher than 300 m³ s^–1. Ice-shelf melt rates near the grounding line may have increased by a factor of three over background values during these events. Between lake drainage events, outflow and consequent basal melting was greatly reduced. Although focusing of outflow and episodic lake discharge raised the peak values for ice-shelf basal melt, the long-term average melt rate was substantially reduced.

The contribution of the subglacial water outflux to the freshwater budget of each embayment can be significant. In the case of hydrologic catchments containing large subglacial lakes near the grounding line, the lake drainage cycle can dominate the freshwater budget of the ice-shelf cavity, particularly for deeply recessed embayments along the grounding line. Although the average contribution of subglacial outflow and consequent enhanced melting of the ice-shelf base comprises a small fraction of the freshwater budget for the RIS cavity, the contribution of outflow variability to total variability in the freshwater budget may be substantial.

Given that water beneath the Siple Coast ice streams is generated entirely by melting at the base, subglacial outflow represents an input to the ice-shelf cavity system that responds more to changes in the dynamics of the interior ice than to changes in the ocean and atmosphere. This differs from marine-terminating outlet glaciers in Greenland or Alaska where the subglacial water system near the grounding line exhibits a seasonal cycle and responds rapidly to climate change (Motyka and others, 2003, 2011). Although subglacial outflow along the Siple Coast appears to be concentrated into small embayments, it may be quite significant to the evolution of the grounding line environment within these locales. The contribution of outflow variability to the local melt rate and freshwater budget appears to be both significant and detectable with current field techniques. Future efforts aimed at understanding the grounding line stability in fast-flowing regions should consider the effect of subglacial outflow more fully.