Articles | Volume 16, issue 9
The Cryosphere, 16, 3575–3599, 2022
The Cryosphere, 16, 3575–3599, 2022
Research article
02 Sep 2022
Research article | 02 Sep 2022

Persistent, extensive channelized drainage modeled beneath Thwaites Glacier, West Antarctica

Persistent, extensive channelized drainage modeled beneath Thwaites Glacier, West Antarctica
Alexander O. Hager1,2, Matthew J. Hoffman2, Stephen F. Price2, and Dustin M. Schroeder3,4 Alexander O. Hager et al.
  • 1Department of Earth Sciences, University of Oregon, Eugene, OR, USA
  • 2Fluid Dynamics and Solid Mechanics Group, Los Alamos National Laboratory, Los Alamos, NM, USA
  • 3Department of Electrical Engineering, Stanford University, Stanford, CA, USA
  • 4Department of Geophysics, Stanford University, Stanford, CA, USA

Correspondence: Alexander O. Hager (


Subglacial hydrology is a leading control on basal friction and the dynamics of glaciers and ice sheets. At low discharge, subglacial water flows through high-pressure, sheet-like systems that lead to low effective pressures. However, at high discharge, subglacial water melts the overlying ice into localized channels that efficiently remove water from the bed, thereby increasing effective pressure and basal friction. Recent observations suggest channelized subglacial flow exists beneath Thwaites Glacier, yet it remains unclear if stable channelization is feasible in West Antarctica, where surface melting is nonexistent and water at the bed is limited. Here, we use the MPAS-Albany Land Ice model to run a suite of over 130 subglacial hydrology simulations of Thwaites Glacier across a wide range of physical parameter choices to assess the likelihood of channelization. We then narrow our range of viable simulations by comparing modeled water thicknesses to previously observed radar specularity content, which indicates flat, spatially extensive water bodies at the bed. In all of our data-compatible simulations, stable channels reliably form within 100–200 km of the grounding line and reach individual discharge rates of 35–110 m3 s−1 at the ice–ocean boundary. While only one to two channels typically form across the 200 km width of the glacier in our simulations, their high efficiency drains water across the entire lateral extent of the glacier. We posit the large catchment size of Thwaites Glacier, its funnel-like geometry, and high basal melt rates together accumulate enough water to form stable channels. No simulations resembled observed specularity content when channelization was disabled. Our results suggest channelized subglacial hydrology has two consequences for Thwaites Glacier dynamics: (i) amplifying submarine melting of the terminus and ice shelf while (ii) simultaneously raising effective pressure within 100 km of the grounding line and increasing basal friction. The distribution of effective pressure implied from our modeling differs from parameterizations typically used in large-scale ice sheet models, suggesting the development of more process-based parameterizations may be necessary.

1 Introduction

Subglacial hydrology is a leading control on basal friction and frontal ablation rates of tidewater glacier termini, yet the morphology of subglacial drainage systems beneath the Antarctic Ice Sheet is poorly characterized. Subglacial water can either flow through a highly pressurized, distributed network of bedrock cavities (Walder1986; Kamb1987), sediment canals (Walder and Fowler1994), films (Weertman1972), and porous till (Clarke1987) or efficiently drain through arborescent channels melted upward into basal ice (Röthlisberger1972). Water flow through a distributed system creates low effective pressures contributing to fast basal sliding (Walder1986; Kamb1987), whereas channelized drainage increases effective pressures (Röthlisberger1972; Schoof2010; Hewitt2011) and local submarine melt rates at the ice–ocean boundary (Slater et al.2015). To date, most models of basin- or ice-sheet-scale Antarctic subglacial drainage have focused on hydropotential mapping (e.g., Stearns et al.2008; Carter and Fricker2012; Le Brocq et al.2013; Livingstone et al.2013; Smith et al.2017) and have only recently distinguished between conduit types under Antarctic glaciers (Dow et al.2020; Wei et al.2020). However, a growing body of work suggests a variety of drainage styles may be important in Antarctica, with obvious relevance to ice sheet dynamics.

In Antarctica, shallow hydropotential gradients and the lack of significant surface melt has led to the conventional paradigm that subglacial water fluxes are too small to permit stable channelized drainage beneath the ice sheets (e.g., Weertman1972; Alley1989; Walder and Fowler1994; Carter et al.2017). This assumption has led to the use of purely distributed subglacial hydrology models (e.g., Alley1996; Le Brocq et al.2009) or simplifying approximations of effective pressure in large-scale Antarctic ice sheet models (e.g., Leguy et al.2014; Asay-Davis et al.2016; Yu et al.2018; Nias et al.2018; Cornford et al.2020). However, channelized drainage under Antarctic ice sheets has recently been inferred through observations of ice shelf basal melt channels (Le Brocq et al.2013; Marsh et al.2016; Drews et al.2017), radar specularity content (Schroeder et al.2013), and subglacial hydrology models (Dow et al.2020; Wei et al.2020). In the absence of surface meltwater, subglacial channels must be sustained through basal melting, and the presence of basal melt channels under ice shelves suggests that their grounded counterparts must persist stably for decades or centuries (Le Brocq et al.2013; Marsh et al.2016).

Thwaites Glacier contains enough ice to raise sea level by 65 cm (Rignot et al.2019) and may currently be undergoing an unstable retreat, likely triggered by increased melting of its ice shelf and terminus (Joughin et al.2014; Rignot et al.2014; Seroussi et al.2017; Milillo et al.2019; Hoffman et al.2019). Ice flux from Thwaites Glacier increased 76 % between 1976–2013 (Mouginot et al.2014), coinciding with thinning rates of up to 10 m yr−1 and a surface acceleration of 100m yr−1 near the grounding line (Pritchard et al.2009; Helm et al.2014; Gardner et al.2018). While bed topography primarily regulates Thwaites Glacier retreat, uncertainty in basal friction laws, ice flow models, and ice shelf melt parameterizations could affect mass loss projections for this century by up to 300 % (Yu et al.2018). As a prominent control on both basal friction and submarine melting, subglacial hydrology has the potential to be a critical component of Thwaites Glacier dynamics, yet the configuration of its drainage network is poorly understood.

Using a recent survey of radar specularity content, Schroeder et al. (2013) hypothesized that channelized subglacial drainage is pervasive within 75–100 km of the Thwaites Glacier grounding line. However, subsequent satellite detection of subglacial lakes led to the interpretation that such channels may only be ephemeral, forming only during lake drainage events (Smith et al.2017). Here, we pair remote sensing with the two-dimensional subglacial hydrology model implemented within the MPAS-Albany Land Ice Model (MALI) (Hoffman et al.2018) to provide a more complete picture of the likely configuration of the Thwaites Glacier subglacial drainage system. We run a suite of 138 modeling simulations, then compare our results with the observed radar specularity content of Schroeder et al. (2013) to define a subset of scenarios as possible representations of reality. Results from this subset are then collated with ice shelf basal melt rates and common parameterizations of basal friction to explore the significance of channelization on submarine melt rates and ice dynamics.

2 Methods

2.1 Model framework

Here, we use only the subglacial hydrology component of MALI, which contains both distributed and channelized flow components and operates on an unstructured, two-dimensional Voronoi grid. Velocities and fluxes are calculated on the edge midpoints of each cell, and all other variables are located at cell centers. Channel segments connect the centers of neighboring cells. The distributed system is treated as a macroporous sheet that is designed to resemble the bulk flow of water through cavities on the lee-sides of bedrock bumps (Flowers and Clarke2002; Hewitt2011; Flowers2015) but may also reasonably describe flow through other porous media, such as till or till canals (Hewitt2011; Flowers2015; Hoffman et al.2016). The distributed system discharge is given by

(1) q = - k q W α 1 ϕ α 2 ϕ ,

where kq is the conductivity coefficient of the distributed system, W is the water thickness, and α1 and α2 are 54 and -12, respectively, to resemble a Darcy–Weisbach flow law. The hydropotential, ϕ, is defined as

(2) ϕ = ρ w g Z b + P w ,

where ρw is the water density, g is the gravitational acceleration, Zb is the bed topography (Fig. 1a), and Pw is the distributed water pressure. It is assumed all basal cavities remain filled, and thus water thickness is a function of cavity opening from basal sliding over bedrock bumps and creep closure:

(3) d W d t = c s u b ( W r - W ) - c cd A b N 3 W ,

where cs is a bed roughness parameter, ub is the ice basal sliding velocity (Fig. 1b), Wr is the maximum bed bump height, ccd is a creep scaling parameter for the distributed system, and Ab is the temperature-dependent ice flow rate parameter for basal ice. The effective pressure, N, is defined as the difference between the ice overburden and water pressures: N=ρigH-Pw, for ice thickness H and ice density ρi.

The channelized system formulation resembles that of Werder et al. (2013), where channel discharge is given by

(4) Q = - k Q S α 1 ϕ α 2 ϕ ,

where kQ is the channel conductivity coefficient. Channel cross-sectional area, S, is a function of creep closure and melting/freezing due to the dissipation of potential energy, Ξ, and pressure-dependent changes to the sensible heat of water, Π:

(5) d S d t = 1 ρ i L ( Ξ - Π ) - C cc A b N 3 S .

Here, L is the latent heat of melting, and Ccc is a creep scaling parameter for channels. Ξ includes dissipation terms for both the distributed and channelized systems, so that

(6) Ξ = Q d ϕ d s + l c q c d ϕ d s ,

where s is the along-channel spatial coordinate, and qc is the discharge in the distributed system within a distance, lc, from the channel. Using this formulation, channels may only develop if there exists sufficient discharge in the distributed system for melting to overcome creep closure. In our experiments, we disabled the pressure-dependent melting/freezing term, Π, to avoid nonphysical instabilities arising from intricate bed topography. The implications of neglecting this term are discussed in Sect. 4.3.

Closing the system of equations requires the conservation of water mass within the combined distributed and channelized subglacial drainage systems and a conservation of energy equation for the production of basal meltwater. Conservation of mass is written as

(7) d W d t = - q - S t + Q s δ ( x c ) + m b ρ w ,

where δ(xc) is the Dirac delta function applied along the locations of the linear channels, and mb is the production of basal meltwater (Fig. 1d). Conservation of energy is written as

(8) m b L = G + u b τ b

for basal shear τb and geothermal flux G.

Time derivatives are discretized using an explicit forward Eulerian method that fulfills advective and diffusive Courant–Friedrichs–Lewy (CFL) conditions for the distributed system and advective CFL conditions for the channelized system. Model outputs are written at 1-month intervals, and all reported results are averaged over 5 years of model time to smooth any minor oscillations remaining in the system.

Figure 1(a) Bed topography (Zb), (b) basal sliding speed (|ub|), (c) basal friction heat flux (Ff), and (d) the production of basal meltwater (mb) used as inputs for the subglacial hydrology model. Transects spaced every 50 km from the terminus (used for determination of flux steady-state and in Fig. 6) are shown as black lines, with the dotted lines spanning the transition zone of Schroeder et al. (2013). The locations of map corners are given in Standard Antarctic Polar Stereographic coordinates. The inset in (d) depicts the location of Thwaites Glacier (blue) within Antarctica.

2.2 Thwaites model domain

We ran the majority of our simulations on a variable resolution domain of Thwaites Glacier that has a 4 km cell spacing over the fast-flowing regions and coarsens to 14 km at the interior ice divide for a total of 4267 grid cells. An additional simulation was performed with a higher-resolution mesh that uses 1 km cell spacing in fast-flowing regions, coarsening to 8 km at the interior ice divide, for a total of 75 500 cells. The bedrock and ice geometry were interpolated onto the model mesh using conservative remapping from the BedMachine Antarctica v1 ice thickness and bed elevation dataset (Morlighem et al.2020). However, a maximum bed elevation of 1200 m and a ice thickness of 550 m were imposed over Mt. Takahe (>250km from the terminus) to avoid instabilities arising from steep bed topography. The resulting thickness gradients were then smoothed by running only the ice dynamics and geometry evolution portions of MALI for 15 years. The geothermal flux was interpolated from the 15 km resolution dataset of Martos et al. (2017). The ice sliding velocity (ub) and basal shear stress (τb) fields required by the subglacial hydrology model follow the methods used by Hoffman et al. (2018) to generate a present-day initial condition, where a basal friction parameter is optimized in order to minimize the misfit between modeled and observed ice surface velocity (Perego et al.2014).

Within the subglacial hydrology model, no flow lateral boundary conditions were applied at the ice-covered lateral boundaries of the model domain. At the glacier grounding line, a Dirichlet boundary condition on the hydropotential (ϕ) was applied equal to the hydropotential of seawater at each grid cell seaward of the grounding line,

(9) ϕ o = ρ w g Z b - ρ o g Z b ,

where ρo=1028kg m−3 is the density of ocean water. Note this boundary condition results in hydropotential values close to zero but spatially varying, as ocean pressure varies along the grounding line with the thickness of the ocean water column. Additionally, inflow from the ocean to the subglacial drainage system is disallowed if the hydropotential underneath the grounded ice falls below the ocean hydropotential. This condition can occur due to a spatially variable, ocean-lateral boundary condition and the assumption of constant density within the subglacial drainage system, which in combination with subglacial channelization can locally result in the modeled unstable inflow of ocean water. The model was spun-up with channelization disabled and a kq value of 1.5×10-3m7/4kg-1/2 to allow water pressures to equilibrate at >90 % overburden pressure. All other simulations were then initialized from the steady-state solution of this run.

2.3 Parameter sweep and sensitivity analysis

Four primary yet poorly constrained parameters exist in Eqs. (1), (3), and (4): kq, kQ, Wr, and cs. While some theoretical and observational basis exists for the values of these parameters, the appropriate values are uncertain and likely vary by glacier basin. A few recent studies have addressed this uncertainty by using inversion techniques to infer values of hydraulic parameters (e.g., Brinkerhoff et al.2016; Koziol and Arnold2017; Irarrazaval et al.2021; Brinkerhoff et al.2021). Here, we used an ensemble approach and compared results to multiple limiting criteria to identify the most realistic parameter combinations. Our ensemble consisted of 113 different channel-enabled simulations and 25 simulations disallowing channelization. All runs were within a plausible parameter space based on observations and theory, as described below.

Observations of jökulhlaups suggest the typical Manning roughness, n, of subglacial channels ranges from 0.023–0.12 sm-1/3 (Nye1976; Clarke1982; Bjornsson1992; Clarke2003). We can translate these Manning roughness values to the equivalent channel conductivity range of 0.03–0.17 m7/4kg-1/2 using (Werder et al.2013)

(10) k Q 2 = 1 ρ w g n 2 ( 2 π ) 2 / 3 ( π + 2 ) 4 / 3 .

However, jökulhlaups do not provide an exhaustive range of roughness characteristics for channel flow, and dye-trace breakthrough curves have indicated that n values for low-discharge, high-friction subglacial channels could be as low as n=0.68sm-1/3 (Gulley et al.2012) or kQ=0.006m7/4kg-1/2. On the other extreme, the Manning roughness of a smooth brass pipe is 0.009 sm-1/3 (Chow1959) or kQ=0.44m7/4kg-1/2, which we consider a generous upper end-member for kQ. We therefore ran our model with kQ ranging from 0.005–0.5 m7/4kg-1/2 to encompass the full set of plausible values.

Because kq may be chosen to portray porous flow through cavities in till or bedrock, we selected kq values to be within the appropriate range of till or greater. Estimates for the hydraulic conductivity, κ, of subglacial till ranges widely from 10−125×10-4m s−1 (Fountain and Walder1998), which can be converted to an equivalent distributed conductivity coefficient in our model via (Bueler and van Pelt2015)

(11) k q = κ ρ w g W 1 / 4 ϕ - 1 / 2 .

Using a characteristic W of 0.1 (see below) and |∇ϕ| of 100 Pa m−1 (approximated from our model domain), we estimate the conductivity coefficient of subglacial till in our model would be 10−1510−6m7/4kg-1/2, which should span our lower limit for kq. In practice, however, simulations with kq<1.5×10-5m7/4kg-1/2 were over-pressurized and did not reach steady-state. Although no proper upper bound exists for kq, we attempted to limit our kq parameter sweep to values that kept the average water pressure >90 % flotation, which typically occurred for kq5×10-3m7/4kg-1/2 across different bed roughness combinations. This choice was based off of near-flotation water pressures observed at Ice Stream B (Engelhardt and Kamb1997) and estimated for Pine Island Glacier (Gillet-Chaulet et al.2016), which we assume are similar to those beneath Thwaites Glacier.

In theory, Wr represents the characteristic bed bump height (decimeter-scale), while cs represents the characteristic meter-scale bed bump spacing (Fig. 2). Typical values used for Wr and cs are ∼0.1m (e.g., Schoof2010; Hewitt2011; Schoof et al.2012; Werder et al.2013; de Fleurian et al.2018; Dow et al.2020) and ∼0.5m−1 (e.g., Schoof et al.2012; Werder et al.2013; Hoffman and Price2014; de Fleurian et al.2018; Dow et al.2020), respectively. We tested the sensitivity of our results to these parameters by running the model with six different combinations of Wr=0.05, 0.1, 0.2, and 1.0 m and cs=0.25, 0.5, and 1.0 m−1, holding one at the default value of Wr=0.1m or cs=0.5m−1 and varying the other parameter. We spaced kq and kQ samples at consistent intervals and stopped sampling conductivity parameter space when runs failed to reach steady-state or were under-pressurized (<90 % flotation). As a result, we conducted a different number of runs for each bed parameter combination, ranging between 9–29 channel-enabled simulations with kq and kQ values within their plausible ranges (Appendix A).

Additionally, for each pair of bed roughness parameters, we ran four to five simulations with channelization disabled across a similar range of kq values (25 runs total). These were used as counter-examples to explore the impact of subglacial channel drainage under Thwaites Glacier.

By design, the parameter sweep forces our model to operate at the limit of its ability to remain stable, and thus some runs failed to reach a true steady-state. This occurred for two main reasons: either local numerical instabilities developed in the channel model or the domain became over-pressurized so that the adaptive time step became impractically small to meet the pressure CFL condition. We thus found it useful to define two separate steady-state criteria that allowed us to identify which information was usable from each run, and we categorized runs as either reaching a pressure steady-state or a flux steady-state. Pressure steady-state was defined as NijtNij-10.5 %, where 〈〉 denotes an average over all grounded grid cells j and time steps i over 5 years of model time. Flux steady-state was attained when the area-integrated version of Eq. (7) upstream of a specified cross-glacier transect was met within 0.5 % when averaged over 5 years. Transects were defined every 50 km within 200 km of the grounding line (Fig. 1). Runs that failed to reach flux steady-state did not represent steady systems where the subglacial discharge realistically balanced the production of meltwater, and so it was not possible to accurately assess the relative fraction of channel discharge to distributed system discharge. Therefore, we report results regarding water thickness and water pressure from pressure steady-state runs but only report discharge results from runs that also reached flux steady-state at each transect.

We use this approach because water pressure and thickness fields from pressure steady-state runs strongly resemble their flux steady-state neighbors in parameter space, yet the channel model fails to reach equilibrium in some runs due to local channel instabilities that do not affect area-averaged water pressure or water thickness. We thus have confidence that pressure steady-state runs still yield useful information about water pressure and thickness. In some cases, instabilities could be avoided by changing the englacial porosity, which acts as a buffer between meltwater production and the subglacial system but does not affect the steady-state configuration. As our goal was to explore as much of parameter space as possible, runs were continually restarted until they either reached flux steady-state, formed an unpreventable numerical instability, or became computationally untenable to keep running. Simulations that did not reach either of the steady-state criteria were discarded. The sensitivity of our results to our steady-state criteria is discussed in Appendix A.

2.4 Model comparison with observed specularity content

All simulations that reached a pressure steady-state were compared with observed radar specularity content from Thwaites Glacier (Schroeder et al.2013) to further narrow the range of viable parameter combinations. Specularity content determined from airborne ice-penetrating radar is commonly used for detecting subglacial water bodies beneath ice sheets (e.g., Schroeder et al.2013, 2015; Young et al.2016, 2017; Dow et al.2020), and it has recently been used to validate a subglacial hydrology model of Totten Glacier, East Antarctica (Dow et al.2020). Although our methods differ, we rely on the same concepts that make specularity content a useful tool for subglacial hydrology model validation.

Ponding within the subglacial drainage system creates flat, reflective surfaces that cause bright specular returns, as opposed to bedrock, which has a lower dielectric contrast to ice and whose rough texture scatters energy (Schroeder et al.2015). Similarly, the curved surface of less uniform conduits such as channels or rough linked cavities scatters energy uniformly in all directions, creating areas of low specularity content despite the presence of water (Schroeder et al.2013). High specularity content, therefore, unequivocally depicts flat-surfaced water bodies in an inefficient distributed system, while low specularity content can represent either a distributed system below its capacity (bedrock cavities are smaller than their maximum size allowed by bed roughness) or the existence of water in rougher, more variably shaped conduits, such as channels. However, by comparing specularity content with a numerical model, we are able to determine which of these two features is responsible for creating the weakly specular regions beneath Thwaites Glacier.

To compare specularity content with our model output, we first averaged the specularity content from the north–south and east–west radar transects from Schroeder et al. (2013) onto a 5 km grid. We then defined a water thickness to bed bump height ratio, Rwt, which indicates the degree to which modeled conditions would produce flat and extensive interfaces between water and ice at the glacier bed and therefore highly specular surfaces:

(12) R wt = W W r .

For Rwt1, distributed water thickness nears or exceeds bed bump height, thus creating a flat, highly specular surface of water. However, for Rwt≪1 bedrock geometry determines the roughness of the lower interface, and the location is considered rough-surfaced and non-specular (Fig. 2). Additionally, with a proper value of kq, Rwt can also parameterize till saturation, with low and high Rwt indicating undersaturated (non-specular) till and saturated (specular) till, respectively. For easy comparison, Rwt was calculated for each model grid cell and then interpolated onto the same 5 km grid as the specularity content data. Note that a spatially uniform Wr is likely unrealistic but is an assumption commonly used in subglacial hydrology models. As applied here, Eq. (12) is used as a relative metric of how close to maximum size a linked cavity system is, and this interpretation would apply to both uniform or spatially variable bump heights.

Figure 2Schematic of a specular and non-specular distributed system, as defined by the water thickness ratio, Rwt. Physical representations of bed roughness parameters are included.


Measured specularity content and modeled Rwt both represent broad, flat areas of pooled water at high values but should not be expected to covary when their values are low due to nonlinearities in the measurements and model formulations, as well as ambiguity in the physical representation of low specularity content. This makes comparing the two difficult, and a simple spatial correlation is unlikely to work as a comparison method. Instead, we rely on binary masks that map where specularity content and Rwt are high/low, as determined by their value being above/below a threshold value. Unfortunately, this method requires choosing thresholds that are considered high for each quantity, which we address by creating a population of masks for each variable, each using a different threshold within a reasonable range.

Specularity content depends on the geometry of ice thickness, survey geometry, radar processing, and subglacial water geometry (Schroeder et al.2013, 2015; Young et al.2016; Haynes et al.2018). As a result, specularity content can be interpreted as the relative amount of the bed that is covered by flat subglacial water bodies, which gradually transitions from non-specular to specular with the addition of water. Therefore, the classification of high or low specularity content is determined relative to a specific survey, and we base the threshold value used for creating specularity masks on the cumulative distribution of specularity within our dataset (Fig. B2). As masks are sensitive to the exact choice of threshold, we created 11 specularity masks with thresholds, Scrt, ranging from 0.15–0.25 at evenly spaced intervals of 0.01, which selects for the greatest ∼5 %–20 % of our specularity data. Similarly, we assume specular surfaces require cavities that are near maximum size (Fig. 2), so there must be a range of Rwt near 1 that could plausibly represent high specularity. Again, we account for this range by creating six masks of Rwt using thresholds, Rwtcrt, between 0.95–1.0 at intervals of 0.01. The resultant 66 combinations of specularity content and Rwt masks were then compared using two criteria:

  1. The masks were divided into four zones based on Schroeder et al. (2013): a near-terminus non-specular zone thought to have channelized flow, a lower specular zone approximately at the transition zone of Schroeder et al. (2013), an upper specular zone where ponding is thought to occur, and an upper non-specular zone likely containing little basal water (Fig. 3). The specularity content and Rwt masks had to agree for a majority of the cells within each zone.

  2. The two masks needed to have an overall correlation coefficient of ≥0.35, which was empirically tuned to select for similar patterns between masks when paired with the first criterion.

Model runs that had at least one Rwt mask that met these comparison criteria with at least one specularity mask were deemed data compatible and used for further analysis. By admitting runs that satisfy the comparison criteria for even a single set of masks out of the 66 compared, we make the selection highly inclusive so that conclusions about extent of channelization consider the widest range of parameters compatible with specularity observations. Hereafter, runs that additionally met flux steady-state criteria will be referred to as data-compatible flux steady-state (FSS) runs. See Appendix B for more information about these comparison criteria, as well as a flow chart illustrating the comparison process (Fig. B1).

Figure 3An example comparison of catchment-scale features identified with binary masks (black) of observed specularity content and modeled Rwt. (a) Radar specularity content (Schroeder et al.2013) and (c) Rwt for a data-compatible flux steady-state model run, together with their coinciding binary masks, (b) Scrt=0.19 and (d) Rwtcrt=0.98, respectively. The dashed pink line in (a) marks the transition between highly specular, distributed drainage, and channel-dominated drainage, as hypothesized in Schroeder et al. (2013). The four zones used for comparison between specularity content and Rwt are color-coded in (b) and (d)​​​​​​​. Light and dark gray lines in (c) are the 50 % and 90 % Rwt contours, respectively. The percent match between masks within each zone and the overall correlation are given in (d). The locations of map corners are given in Standard Antarctic Polar Stereographic coordinates.

3 Results

3.1 Channel-enabled parameter sweep

3.1.1 Model tuning and correspondence with specularity content

Of our 113 channel-enabled runs, 39 met our pressure steady-state criterion, while 23 of those also met our flux steady-state criterion across all transects. A total of 20 pressure steady-state runs, including 13 flux steady-state runs, had at least one Rwt and specularity mask combination that met our comparison criteria, and they were therefore considered possible representations of reality. Each of these runs showed a strong resemblance between Rwt and specularity content masks (Fig. 3). Average water pressures in data-compatible runs were between 91 %–96 % flotation, and in general, runs that did not correspond with specularity content had water pressures outside of this range.

All 66 combinations of Scrt and Rwtcrt masks yielded successful comparisons for some sets of parameters, although successful pairings varied with model parameters. Across all runs, comparison success rate exponentially increased with higher values of Rwtcrt, with Rwtcrt of 0.99 or 1.0 accounting for 60 % of all matches. Conversely, masks with Rwtcrt=0.95 only accounted for 4 % of the 713 successful mask combinations. The few runs that had successful matches with an Rwtcrt of 0.95 also had successful matches using higher Rwtcrt thresholds, indicating this choice of lower bound does not influence our results. Match success rate was not sensitive to Scrt, and each threshold value was responsible for 7 %–10 % of successful matches.

Data-compatible runs either had kq values of 1.5×10-4 or 5×10-4m7/4kg-1/2 (Fig. 4), with the only exceptions occurring when Wr=0.05m or Wr=1.0m, in which data-compatible kq values reached 1.5×10-3 and 5×10-5m7/4kg-1/2, respectively. The range of kq in data-compatible runs is above that of pure glacial till and is consistent with a bed composed of both till and bedrock, as is thought to be the case for Thwaites Glacier (Joughin et al.2009; Muto et al.2019b, a). For the channelized conductivity values, all data-compatible runs had kQ values of 0.005–0.1 m7/4kg-1/2, coinciding with the expected range given by dye-trace breakthrough curves and jökulhlaup observations (Nye1976; Clarke1982; Bjornsson1992; Clarke2003; Gulley et al.2012). No runs with kQ=0.5m7/4kg-1/2, outside of our brass pipe upper limit, reached either steady-state criterion. Typical channel velocities in our data-compatible runs do not exceed the typical observed jökulhlaup range of 0.6–2.7 m s−1 (Magnusson et al.2007; Werder and Funk2009, Fig. 4), which provides an additional loose constraint on the validity of our channel model, although currently no observations of subglacial flow velocities exist from Antarctica.

Figure 4The conductivity parameter sweep for bed roughness parameters Wr=0.1m and cs=0.25m−1. Stars represent runs that reached flux (and pressure) steady-state, triangles symbolize pressure steady-state simulations, and filled black circles depict runs that did not reach either steady-state criterion. Symbols for steady-state runs are color-coded by the average flotation percentage of grounded ice. Circles around stars or triangles indicate runs that matched observed specularity content and are considered data compatible. Gray lines are 95th percentile channel velocity contours for channels with Q>5m3 s−1. kQ limits determined from a brass pipe and dye-trace breakthrough curves are plotted as dashed brown and dark-blue lines, respectively, and the blue shaded area represents the typical observed jökulhlaup kQ range.


3.1.2 Extent of channelization in data-compatible simulations

Subglacial channels were ubiquitous in all data-compatible FSS runs. In most of these runs, channels with discharges over 5 m3 s−1 extended at least 150 km from the glacier terminus, with some channels reaching farther than 200 km (Fig. 5). The initiation of these channels generally coincided with the upper specular zone observed in Schroeder et al. (2013). However, channel discharge between 150–200 km was divided between two to four small channels, each with an individual discharge of less than 20 m3 s−1. At 150 km from the terminus, distributed discharge was still the dominant mode of drainage, with average channelized and distributed discharges of 27 ± 18 and 42 ± 19 m3 s−1 (± indicates standard deviations), respectively, across data-compatible runs.

Figure 5(a) Average effective pressure and channel discharge across all data-compatible FSS runs. (b–c) Effective pressure and channel discharge for (c) the high-resolution model and (b) its low-resolution counterpart. The insets are enlarged views of the black boxes, and the star in (a) indicates the location of the secondary channel seen in one data-compatible FSS run. Sub-ice-shelf melt rates from Adusumilli et al. (2020) are plotted in all frames. For clarity, only channels with Q>5m3 s−1 are pictured in each frame. Again, transects spaced every 50 km from the terminus (used for determination of flux steady-state and in Fig. 6) are shown as black lines, with the dotted lines spanning the transition zone of Schroeder et al. (2013). The locations of map corners are given in Standard Antarctic Polar Stereographic coordinates.

A transition occurs between 50–100 km from the terminus from a distributed-dominated to a channel-dominated system, coinciding with the region where Schroeder et al. (2013) hypothesized channelization begins under Thwaites Glacier. In our model, all data-compatible runs had formed at least one channel transporting >10m3 s−1 by 100 km from the terminus, and by 50 km, these channels had grown and converged into one to two primary channels, each draining up to 50 m3 s−1 of water. Our 50 km transect is the first at which channelized drainage slightly outweighs distributed drainage, with discharges of 55 ± 21 and 47 ± 20 m3 s−1, respectively (Fig. 6). Consistent with Joughin et al. (2009), basal friction melting is the primary contributor of melt in our model, and the 50–100 km transition to channelized flow coincides with a substantial increase in basal friction melt rate (Figs. 1, 5, and 6).

Figure 6(a) Total distributed (blue) and channel discharge (gray), as well as the discharge of the largest channel (red), across each transect (see Fig. 5) for all data-compatible FSS runs (circles). Boxplots indicate the maximum, minimum, mean, and standard deviations. The stars indicate the high-resolution model, and the white-edged circles designate its low-resolution counterpart. (b) The number of channels with Q>5m3 s−1 (gray) and Q>10m3 s−1 (blue) at each transect for all data-compatible FSS runs.


Channelized discharge grows rapidly within 50 km of the terminus. By the point at which water reaches the grounding line, channelized drainage accounts for 127 ± 24 m3 s−1 of runoff into the ocean, whereas only 25 ± 21 m3 s−1 is expelled through the distributed system (Fig. 6). In all data-compatible FSS runs, the majority of channel discharge at the grounding line occurred through one primary channel with a discharge of 80 ± 24 m3 s−1 near the center of the grounding line (-1.5369×106, -4.7298×105m; standard Antarctic Polar Stereographic). This location corresponds to the region of high basal melting observed at the Thwaites Ice Shelf in Adusumilli et al. (2020) (Fig. 5). In one simulation, a secondary channel intersects the grounding line with a discharge of 38 m3 s−1 at -1.5310×106, -4.8585×105m, where we lack basal melt data (Fig. 5a). Other channelized discharge across the grounding line occurs through very small channels (10m3 s−1) scattered along the marine boundary.

3.2 Grid resolution sensitivity analysis

One data-compatible FSS simulation (kQ=0.05m7/4kg-1/2, kq=4×10-4m7/4kg-1/2, cs=0.5m−1, Wr=0.1m) was rerun to flux steady-state with the high-resolution domain. The high-resolution run matched observed specularity content and produced effective pressures and water fluxes that closely resembled its low-resolution counterpart. High-resolution channels followed very similar pathways as those in the low-resolution model (Fig. 5b–c), and distributed and channelized discharges at each transect were approximately equal to those at low resolution (Fig. 6a). The main exception occurred at the grounding line, where the two main channels reached the ocean independently in the high-resolution model but merged just above the grounding line with lower resolution (Fig. 5b–c). This explains the almost twofold discrepancy of maximum channel discharge at the grounding line between the two resolutions (Fig. 6a). Additionally, the high-resolution run had lower effective pressures near the upper domain boundary, although effective pressures within 300 km of the terminus are in strong agreement with the low-resolution model (Figs. 5, 8a).

3.3 Distributed-only model configuration

Average water pressures in our 25 distributed-only simulations ranged from 74 %–98 % flotation, and all met our flux steady-state criteria. However, no distributed-only run had a Rwt field that matched observed specularity content. In particular, the greatest mismatch occurred between 0–50 and 100–150 km of the grounding line, where Rwt was consistently over Rwtcrt but where observed specularity content was low (Fig. 7). In other cases where the average flotation percentage was below 90 %, water thicknesses were too low to produce any regions of RwtRwtcrt. Furthermore, distributed-only simulations had unrealistically low effective pressures within 150 km of the terminus. Of the runs with an average water pressure over 90 % flotation, many were at or near flotation within 200 km of the terminus (Fig. 8a). Within 50 km of the terminus, the average effective pressure across these distributed-only runs was one-third that of data-compatible channel-enabled scenarios.

Figure 7Three typical (a–c) Rwt configurations and (d–f) coinciding binary masks (black) for distributed-only runs. Masks depict regions where Rwt is above its threshold value, and thus the distributed system is at or above its capacity. Rwt in distributed-only runs generally resembled one of these three patterns. Light and dark gray lines in (a)(c) are the 50 % and 90 % Rwt contours, respectively. Color-coding in (d)(f) corresponds to the same zones as in Fig. 3. Purple line in (a) is the center-line transect used in Fig. 8. kq values used in each run, along with the Rwtcrt used to create the coinciding binary mask, are provided in (d)(f). All three runs had bed roughness parameters Wr=0.1m and cs=0.5m−1.

Figure 8(a) The range and mean (solid line) of effective pressures along the center-line transect in Fig. 7a for all data-compatible FSS, channel-enabled runs (magenta), and all distributed-only runs above 90 % flotation (gray). The black line depicts transect effective pressures from the high-resolution run. Shown in blue is the calculated effective pressure if assuming a perfect hydrostatic connection with the ocean. Note the different y-axis scales in (a). (b) Basal shear stress used as input in our model (red) plotted with reconstructed basal shear stress using a Budd-style friction law (blue). Blue hues represent different exponents used in the friction law. All lines follow the same center-line transect as in (a).


4 Discussion

4.1 A reconciled framework for channelization beneath Thwaites Glacier

The key result of our study is the likely existence of stable subglacial channels beneath Thwaites Glacier. In our model, channels typically extended over 100–200 km inland and had grounding line discharges of 80±24m3 s−1, much larger than the maximum discharges of 1–5 and <25m3 s−1 modeled at Getz (Wei et al.2020) and Totten (Dow et al.2020) glaciers, respectively. No distributed-only experiments matched observed specularity content, and all had unrealistically high water pressures within 100 km of the terminus. This strongly argues that channelized drainage is necessary to explain observed radar specularity content.

Certain geometric and hydrologic conditions at Thwaites Glacier are unfavorable to the development of subglacial channels, and thus the extent of channelization in our model is somewhat surprising. In theory, subglacial channels should develop when the distributed system reaches a critical discharge that is inversely proportional to the hydropotential gradient (Schoof2010; Hewitt2011). In Greenland, it is believed that glaciers are unable to reach this critical threshold farther inland where gentle surface slopes weaken the hydropotential gradient and thick ice may expedite creep closure (Chandler et al.2013; Meierbachtol et al.2013; Dow et al.2014). Similar logic could also apply to the thicker and broader Antarctic ice sheets, especially given their insignificant surface melt input. Yet, our model consistently depicts subglacial channels extending 100–200 km inland in all parameter choices. These channels could be explained by the large catchment size (189 000 km2) of Thwaites Glacier (Joughin et al.2009), its funnel-like geometry, and high basal melt rates of 3.5 km3 yr−1 (Joughin et al.2009), which together accumulate enough water to exceed the critical discharge threshold within 100–200 km from the grounding line. At first, the critical discharge may only be met locally (e.g., Hewitt2011) through the accrual of water in topographic depressions, which the subglacial channels tend to follow. High basal friction melt rates of 100–1000 mm yr−1 in the terminal 100 km, as calculated for our model input and by Joughin et al. (2009), are then likely responsible for the increased channelization near the grounding line.

Previous work has offered contrasting hypotheses on the persistence of subglacial channels beneath Thwaites Glacier. Originally, Schroeder et al. (2013) argued radar scattering from widespread concave channels produced the near-terminus, non-specular region they observed. However, an extensive channelized system may not allow for the isolation of subglacial lakes, and the discovery of subglacial lakes beneath Thwaites Glacier suggested channels may be ephemeral, forming only during subglacial lake drainage events (Smith et al.2017). Based on our model, we here present a refinement of the hypothesis of Schroeder et al. (2013) that leaves room for the development of the subglacial lakes observed by Smith et al. (2017).

In agreement with Schroeder et al. (2013), we interpret the overlapping regions of observed high specularity content and high Rwt between 100–250 km from the terminus to unequivocally indicate the pooling of broad, flat water bodies in a distributed system near or at its capacity. This distributed-dominated system then transitions to a channel-dominated system between 50–100 km from the terminus. Schroeder et al. (2013) hypothesized this transition to channelized flow occurs through the development of many channels spread across the glacier width, which scatter radar energy and lower specularity; however, our modeling instead suggests that the near-terminus, non-specular zone of Schroeder et al. (2013) depicts a below-capacity distributed system, whose water has been partially evacuated by a small number of large, stable channels. Such a configuration would produce non-specular radar returns due to a rough surface of discontinuous water cavities at a variety of orientations (Fig. 2).

In such a sparsely channelized system, it is expected that isolated areas of the bed exist in which subglacial lakes may form. Disconnected portions of the drainage network are common beneath alpine and Greenland glaciers, particularly in the summer when channels draw water from the surrounding distributed system, leading to the isolation of poorly connected basal cavities (Murray and Clarke1995; Gordon et al.1998; Andrews et al.2014; Hoffman et al.2016; Chu et al.2016; Rada and Schoof2018). Disconnected areas may exist year-round or may reconnect following a reconfiguration of the channelized system or the collapse of channels in the winter (Hoffman et al.2016; Rada and Schoof2018). However, substantial subannual reshaping of the drainage system should not occur in the absence of a seasonal melt cycle, like at Thwaites Glacier, and thus parts of the bed may remain disconnected for extended periods of time. This would allow disconnected water to gradually pool into lakes that drain when they periodically exceed their hydropotential seal (Fowler1999). Such drainage events could act as similar catalysts for drainage network reconfigurations as the seasonal melt cycles of alpine and Greenland glaciers. Our model lacks the complete physics to properly simulate the filling and draining of subglacial lakes (e.g., Carter et al.2017); however, it is evident that persistent and extensive subglacial channels can exist concurrently with subglacial lakes beneath Thwaites Glacier, and further work is needed to understand the interaction between the two drainage features.

4.2 Implications of channelization on Thwaites Glacier dynamics

4.2.1 Channelization and submarine melting at the grounding line

The rapid and potentially unstable retreat of Thwaites Glacier is likely driven by enhanced sub-ice-shelf melting (Rignot et al.2014; Joughin et al.2014; Seroussi et al.2017; Yu et al.2018; Milillo et al.2019; Hoffman et al.2019), resulting in part from intruding warm Circumpolar Deep Water (CDW) flowing along bathymetric troughs to the grounding line (Nakayama et al.2019; Milillo et al.2019; Hogan et al.2020). The most rapid retreat (12–18 km between 1992–2011) was recorded at the glacier's central, fast-flowing core (Rignot et al.2014), where the retreat had continued at a rate of 0.6 km yr−1 until at least 2017 (Milillo et al.2019). Ice shelf submarine melt rates exceed 200m yr−1 at the fast-flowing core, coincident with the recent formation of a prominent sub-shelf cavity (Adusumilli et al.2020; Bevan et al.2021).

In all but one of our low-resolution data-compatible FSS runs, both main channels converge near the grounding line directly above the subshelf cavity described in Bevan et al. (2021) (Fig. 5a–b). In our high-resolution model, one channel intersects the grounding line at this location, while the second reaches the ocean 16 km to the east, also in the region of high subshelf melting (Fig. 5c). Subglacial discharge plumes, formed from channelized subglacial water entering the ocean, amplify local submarine melting through turbulent heating and the entrainment of deep and often warm water, such as CDW, along the terminus and ice shelf (Jenkins2011; Slater et al.2015; Asay-Davis et al.2017). While it would be an over-interpretation of our model to regard the exact locations of subglacial channels as reality, the ubiquitous conjunction of large channels (33–106 m3 s−1) with high subshelf melt rates at the grounding line in all data-compatible scenarios strongly suggests channelized subglacial discharge augments submarine melting in this region. Recent ocean modeling of the Pine Island Ice Shelf cavity supports this assertion and indicates that subglacial discharge localized at the grounding line and of similar magnitude to what occurs in our model can explain the local ice shelf melt rates of ∼200m yr−1 observed at Pine Island Glacier (Nakayama et al.2021). Similar results have also been reported for the nearby Getz Ice Shelf, where subglacial discharge accelerates subshelf submarine melting by entraining and displacing CDW along the base of the ice shelf (Wei et al.2020).

Additionally, CDW reaches the Thwaites Glacier grounding line through a series of bathymetric troughs and sills that moderate its flow (Nakayama et al.2019; Hogan et al.2020), and it is possible the entrainment of ambient water into subglacial discharge plumes may further enhance CDW flushing of the Thwaites subshelf cavity, similar to the subglacial plume-driven renewal of Greenland fjords (Gladish et al.2015; Carroll et al.2017; Zhao et al.2021). However, plume-driven buoyancy forcing may only have a minimal effect on cavity circulation beneath the Pine Island Ice Shelf (Nakayama et al.2021), and thus it could be assumed that the comparable grounding line fluxes given by our model are still too weak to significantly enhance CDW advection to Thwaites Glacier.

4.2.2 Implications of channelization for effective pressure and basal sliding

Despite contributing to high ice shelf basal melt rates and potential loss of ice shelf buttressing, our model suggests subglacial channels may have a stabilizing effect on basal drag near the grounding line. Effective pressures are 3 times higher within 50 km of the grounding line in channel-enabled runs than in distributed-only runs (Fig. 8d). This region of high effective pressure coincides with a distributed system that is operating below its capacity (Fig. 3), something not reproducible in distributed-only simulations (Fig. 8a–c). Only one to three principal channels exist within the terminal 100 km; nevertheless, comparison with distributed-only experiments indicates that a small number of channels are still able to efficiently evacuate water from the entire region due to their lower hydropotential compared to the surrounding area. Higher effective pressure in the terminal 100 km implies higher basal friction, which has been shown to be a leading control on the retreat and mass loss of Thwaites Glacier (Yu et al.2018) and surface velocities at the neighboring Pine Island Glacier (Gillet-Chaulet et al.2016; Joughin et al.2019). High basal shear stress associated with competent bedrock is already thought to exist within 80 km of the grounding line (Joughin et al.2009) and may work in tandem with channelized subglacial drainage to help buttress against further retreat.

Effective pressures decrease substantially further inland where channelization is minimal. In the upper highly specular area, average effective pressures in data-compatible runs range between 200–600 kPa, almost an order of magnitude less than the near-terminus region (Fig. 8). Effective pressures in highly specular areas are similar to the −30–150 kPa effective pressures observed at Ice Stream B (Engelhardt and Kamb1997), which to our knowledge remain the only direct observations of effective pressures in West Antarctica.

Smith et al. (2017) noted that the small (<10 %) increase in ice velocity observed after subglacial lake drainage events may indicate an insensitivity of Thwaites Glacier dynamics to its subglacial hydrology. However, the linked subglacial lake drainage event measured by Smith et al. (2017) beneath Thwaites Glacier in 2013–2014 had an average discharge of 160240m3 s−1 over 6 months – only 3–5 times greater than modeled channel discharge 50 km from the terminus and 1–2 times greater than the largest modeled channels at the grounding line. Any pre-existing channels of similar size to those in our model could, therefore, help accommodate the additional flux from lake drainage events, which may explain the relatively minor increase in ice velocity they observed. Thus, this lake drainage event could also be interpreted as evidence of channelized drainage stabilizing glacier dynamics, as is indicated by our model. As Thwaites Glacier continues to thin and retreat, we expect the subsequent changes in glacier geometry and meltwater input to continually reshape its subglacial drainage network. Our results suggest this will alter ice dynamics, and it should be taken into account when considering the uncertainty in model projections.

Ice dynamics models have recently started implementing effective pressure-dependent sliding laws supported by current theory. However, a challenging problem is how to best parameterize effective pressure in order to solve for basal shear stress. A common approach is to approximate effective pressure by assuming a perfect hydrostatic connection with the ocean (e.g., Leguy et al.2014; Asay-Davis et al.2016; Yu et al.2018; Nias et al.2018; Cornford et al.2020, and others), shown for our model domain in Fig. 8a. Effective pressure using an ocean connection assumption is in fair agreement with our channel-enabled runs within 5 km of the grounding line but is up to an order of magnitude too high further inland, indicating a parameterization based on an open ocean connection may only be realistic near the terminus. This suggests a regularized-Coulomb friction law (e.g., Joughin et al.2019) may be appropriate for Thwaites Glacier, as it only accounts for effective pressure where effective pressure is low and basal sliding speeds are high, such as near the grounding line (Schoof2005). However, our channel-enabled model indicates effective pressure actually decreases between 5–100 km from the grounding line and maintains its proportionality to basal shear stress throughout the entire domain (Fig. 8b). This implies basal shear stress stays in the Coulomb regime even within the glacier interior, and thus a yield stress or semi-plastic Budd-type law may work equally well for Thwaites Glacier, as has previously been successful at Pine Island Glacier in reproducing observed surface velocities (Gillet-Chaulet et al.2016).

To test this hypothesis we attempt to reconstruct our input basal shear stress using a Budd-style friction law of the form τb=CNub1/m, where ub is a model input, N is solved for by the hydrology model, and C is a tunable basal slipperiness coefficient. Here, m is the bed-dependent stress exponent that is likely between 5–10 for Pine Island Glacier (Gillet-Chaulet et al.2016; Nias et al.2018; Joughin et al.2019), which is assumed to have similar basal properties to Thwaites Glacier. Figure 8b illustrates the results using four plausible values of m and accompanying C values that minimize the root mean square error with the model input. All four versions effectively recover the input basal shear stress, with the best agreement using m=5 or m=8, which is consistent with previous work (Gillet-Chaulet et al.2016; Nias et al.2018; Joughin et al.2019). Therefore, we assert that a Budd-style friction law is appropriate for Thwaites Glacier, assuming accurate knowledge of the effective pressure field. Based on these results we caution against the continued usage of the hydrostatic ocean connection parameterization for effective pressures beyond the marginal 5 km for Thwaites Glacier, which may produce unrealistically slow sliding velocities.

4.3 Model considerations

Our results highlight the need for validation of subglacial hydrology models across the entirety of a glacier. We found a wide range of parameter values resulted in steady-state configurations, and most had some degree of channelization coincident with the location of observed anomalously high sub-ice-shelf melting. However, many simulations had water pressures and discharges that were either too low or too high to be realistic, and without comparison with radar specularity content, it would have been easy to arbitrarily choose the wrong parameters and base our conclusions on an unrealistic model. Borehole validation has been previously attempted for a small alpine glacier (Rada and Schoof2018), but the scale of Antarctic and Greenland glaciers makes this unattainable for ice sheets. We therefore suggest that ice-penetrating radar, such as that used in this paper and in Dow et al. (2020), or other broad-scale proxies for basal water, is the best approach for validation of ice sheet subglacial hydrology models. While our comparison between Rwt and specularity content is somewhat ad hoc, it selected for a coherent grouping of parameters, water pressures, and channel velocities within the expected realistic range, which gives us confidence in its effectiveness. Comparison criteria may need customization to be applicable to other glaciers, but the overall methodology presented in this paper should be beneficial in many settings. Bed conditions differ within and between glacier basins, and we stress our parameter choices should not be extrapolated to other glaciers without validation.

Many assumptions built into subglacial hydrology models remain unsupported, and it is uncertain how such assumptions may influence our results. We therefore deem it necessary to consider the primary underlying simplifications that may impact this paper. Our choice to ignore pressure-dependent melting/freezing in Eq. (5) neglects the effects of supercooling, which would lead to the abatement of R-channels and the expansion of the distributed system as water flows out of a prominent overdeepening. Supercooling has been shown to decrease channelization in other subglacial hydrology models (de Fleurian et al.2018). However, the overdeepening within 100 km from the grounding line (Fig. 1a), in which channelization becomes pronounced, is far from meeting the supercooling threshold of Werder (2016). Furthermore, the upward bed slope in the terminal 100 km is only 60 % of the downward surface slope and should therefore allow for sufficient dissipative heating to continually grow channels (Alley et al.1998). We therefore do not expect the neglect of Π in Eq. (5) to significantly affect our conclusions.

Uniform parameterizations of the distributed system do not account for realistic heterogeneity in bed geometry or lithology, both of which can locally influence distributed connectivity (Murray and Clarke1995; Gordon et al.1998; Andrews et al.2014; Hoffman et al.2016; Rada and Schoof2018; Downs et al.2018). The bed of Thwaites Glacier is thought to consist of alternating regions of bedrock and glacial till (Joughin et al.2009; Muto et al.2019a, b; Holschuh et al.2020) that could potentially affect the connectivity of the distributed system and thus conductivity and discharge. Currently, all subglacial hydrology models assume a consistent kq across their domains, although allowing kq to vary with bed lithology may account for spatial differences in connectivity and produce more realistic results (Hoffman et al.2016).

Modeling (Joughin et al.2009) and seismic data (Muto et al.2019b, a) suggest bed elevation could serve as a reasonable proxy for bed lithology under Thwaites Glacier, where subglacial till (low conductivity) accumulates in depressions, and exposed bedrock (high conductivity) primarily exists at topographic highs. Regions of high specularity content coincide with low-lying troughs, and it is therefore conceivable that imposing a high kq above these troughs, and low kq within them, could reproduce the observed specularity content without the need for channelization. However, our results suggest the minimum kq necessary to prevent channelization would still be high enough over a majority of the domain to drop water pressures below realistic levels. Lowering kq within troughs but maintaining the same kq at higher elevations as used in our data-compatible FSS runs could help pool water into subglacial lakes in till-laden depressions (see Sect. 4.1), but it seems unlikely this would divert enough water to preclude the overall growth of channels in the terminal 100 km. Furthermore, the location of modeled channelized flow at the grounding line presents a convincing explanation for the anomalously high sub-ice-shelf melt rates observed at the same position, something that would be lacking in a purely distributed system. We acknowledge the neglect of a spatially variable kq could create some uncertainty in our discharge results but is likely minimal, and our kq parameter sweep may already account for this variability.

As described in Downs et al. (2018), the value of kq used in subglacial hydrology models is a proxy for the connectivity of orifices linking cavities in the bed. Models assume the orifices scale with cavity size; however, in their original conception, orifices behave like small ​​​​​​​R-channels that may enlarge with turbulent melting (Kamb1987; Fowler1987). Downs et al. (2018) used this argument to scale kq with meltwater input, which better captured seasonal water pressures. Although Thwaites Glacier lacks a seasonal meltwater cycle, we could use the same argument to justify the use of a different distributed system flow law.

Darcy or Darcy–Weisbach flow laws are used almost ubiquitously in subglacial hydrology models (e.g., Schoof2010; Hewitt2011; Werder et al.2013; Hewitt2013; Hoffman and Price2014; Downs et al.2018; Hoffman et al.2018; de Fleurian et al.2018; Dow et al.2020, and others), yet these laws are largely unvalidated in the subglacial environment. Distributed discharge with a Darcy–Weisbach turbulent flow law, as used in this paper, has a 54 power dependency with water thickness. However, in other flow laws, such as Darcy porous media flow or Poiseuille laminar flow, the exponent may vary between 1 and 3 (e.g., Hewitt2011, 2013; Kyrke-Smith and Fowler2014; Kyrke-Smith et al.2014). In practice, the use of a higher exponent could produce similar behavior to a melt-dependent kq and could account for a larger connectivity with increased meltwater, driven by the dissipative melting and opening of orifices. Although such a flow law would increase efficiency of the distributed system and potentially minimize channelization, we do not believe its use would dramatically change our results. Water thicknesses using a Darcy–Weisbach law are fairly uniform within 200 km of the grounding line (Fig. 3c), which suggests an increased dependency of discharge on water thickness may make little difference in our model.

5 Conclusions

This paper leverages observations from a variety of sources to select for the subglacial hydrology model scenarios that are the most likely representations of reality. Our range of possible steady-state scenarios highlights the need for thorough parameter sweeps in subglacial hydrology models, which are then winnowed to the most realistic grouping of simulations based on extensive observations. We emphasize validation of subglacial hydrology models within the glacier interior, and not just at its terminus, is necessary to properly constrain realistic drainage behavior. Furthermore, our work demonstrates subglacial hydrology models still produce a range of results that are compatible with data, and thus model results should be reported as a suite of possible scenarios instead of one feasible configuration.

Our work presents an updated conceptual model for the subglacial drainage system beneath Thwaites Glacier. Our model indicates a few stable channels exist within 200 km of the grounding line and coalesce into one to two large stable channels within the terminal 50–100 km. These channels intersect the ice–ocean boundary directly at the location of highest sub-ice-shelf melt rates, suggesting they play an important role in frontal ablation and grounding line retreat. However, in the interior of the glacier, subglacial channels efficiently evacuate water from a broad portion of the bed, thereby increasing basal friction within 100 km of the grounding line and potentially buttressing against further retreat. At this point, it remains unclear how common such drainage systems are in Antarctica or what impact subglacial channels have on sub-ice-shelf cavity circulation and ice dynamics. We expect the subglacial drainage network to continually reconfigure with future changes in meltwater production and glacier geometry, which will subsequently lead to spatially and temporally evolving basal shear stress and frontal ablation rates. Further work with a fully coupled ice dynamics–subglacial hydrology model will be necessary to determine the exact influence of subglacial channels on future retreat and mass loss.

Appendix A: Parameter sweep, sensitivity analysis, and steady-state criteria

A full sweep of realistic conductivity parameter space was conducted for each set of bed roughness parameters; however, our method for determining bed roughness parameters was closer to that of a sensitivity analysis (varying one parameter at a time). This choice was made because real physical constraints exist for conductivity parameters, while bed roughness parameters are theoretical quantities approximating general bed characteristics that only have indirect physical corollaries. A sensitivity analysis is thus more suitable for bed roughness parameters and allowed us to ease the complexity of sampling a four-dimensional parameter space. Results for each conductivity parameter sweep (in addition to Fig. 4) are depicted in Figs. A1A5.

Figure A1Same as Fig. 4 but for bed roughness parameters Wr=0.1m and cs=0.5m−1. Stars represent runs that reached flux (and pressure) steady-state, triangles symbolize pressure steady-state simulations, and filled black circles depict runs that did not reach either steady-state criterion. Symbols for steady-state runs are color-coded by the average flotation percentage of grounded ice. Circles around stars or triangles indicate runs that matched observed specularity content and are considered data compatible. Gray lines are 95th percentile channel velocity contours for channels with Q>5m3 s−1. kQ limits determined from a brass pipe and dye-trace breakthrough curves are plotted as dashed brown and dark-blue lines, respectively, and the blue shaded area represents the typical observed jökulhlaup kQ range.


Establishing steady-state criteria inherently involves defining a cutoff threshold for acceptable noise remaining in the model. For our pressure steady-state runs, effective pressure at each cell is allowed to fluctuate 0.5 % of its value on average. This equates to an allowable fluctuation of roughly 1 kPa where effective pressure is lowest (∼200kPa) and 10 kPa where effective pressure is highest (∼2000kPa). For flux steady-state runs, meltwater production above each transect must equal the total discharge across the transect within 0.5 %. Total melt production above the grounding line is roughly 155 m3 s−1, so our flux steady-state criteria require that we know the total grounding line discharge within 0.8 m3 s−1, which is orders of magnitude less than the uncertainty between data-compatible FSS runs.

Figure A2Same as Fig. 4 but for bed roughness parameters Wr=0.1m and cs=1.0m−1.


Figure A3Same as Fig. 4 but for bed roughness parameters Wr=0.05m and cs=0.5m−1.


Figure A4Same as Fig. 4 but for bed roughness parameters Wr=0.2m and cs=0.5m−1.


Figure A5Same as Fig. 4 but for bed roughness parameters Wr=1.0m and cs=0.5m−1.


Appendix B: Comparison criteria between modeled and observed specularity content

Matching specularity and Rwt masks is a comparison between two spatial point patterns, which can be challenging as it requires a global statistic that can recognize local patterns of point clusters. Other comparisons of spatial point patterns have relied on segmenting the domain into areal units and determining an overall similarity statistic across all units (e.g., Andresen2009, 2016). The method developed in the current paper shares the concept of areal units by defining four physically based zones within which we assess similarity between the two masks. These zones are intentionally chosen to loosely encompass regions of specularity or non-specularity, which allows for some spatial variability between masks and decreases the sensitivity to the zonal boundaries. We then require the two specularity and Rwt masks to match at 50 % or more of grid points within each zone. As low specularity can occur for a variety of reasons, segmenting the domain into specularity-based zones does not predetermine a specific drainage style but preserves the specular pattern of interest and allows us to test hypotheses concerning its formation.

While the first criterion does well by itself in selecting positive matches, it also selects many false positives. This occurs when the Rwt mask is almost entirely non-specular and when over 50 % of the cells in each zone are non-specular in the observed specularity mask (Fig. B3h–i). It was therefore necessary to include a second criterion that can remove these false positives, which we do by requiring an overall correlation coefficient of r≥0.35. Correlations are calculated with

(B1) r = m n ( S m n - S ) ( R m n - R ) ( m n ( S m n - S ) 2 ) ( m n ( R m n - R ) 2 ) ,

where S and R are the specularity and Rwt masks, respectively. Again, correlation by itself does a fair job at identifying positive matches, but it also identifies false positives when the Rwt mask is overly specular (Fig. B3a–b). As the two criteria fail for opposing reasons, they can check and balance each other if the thresholds are tuned appropriately (Fig. B3d–e). We acknowledge this comparison method is sensitive to multiple choices of thresholds, so we attempt to make our criteria for selecting data-compatible runs as generous and inclusive as possible while still removing runs that clearly do a poor job at resembling observations. We empirically determined that requiring ≥50 % of cells in each zone to agree and r≥0.35 works well at identifying positive matches and is sufficiently general to allow a reasonable variety of Rwt masks to pass this filtering process.

Figure B1Flow chart illustrating the step-by-step process for determining which model runs were compatible with observed specularity content.


Figure B2Cumulative density function of observed specularity data from Schroeder et al. (2013). The green band highlights the range of specularity values used to create our 11 specularity masks, which are in the 81st to 94th percentile of our dataset.


Figure B3Select Scrt (observed) and Rwtcrt (modeled) mask combinations from the Wr=0.1m, cs=0.5m−1, kq=5×10-2m7/4kg-1/2, and kQ=5×10-4m7/4kg-1/2 model run, plotted over the four zones used for the first comparison criterion (also shown in Figs. 3, 7). Nzns indicates the number of zones that meet criterion 1, and r is the overall correlation between mask pairs. Background color indicates successful (green) and unsuccessful (red) matches. Values of Scrt and Rwtcrt used to make each mask are displayed above each plot.

Code and data availability

Model output, radar data, processing files, and the MALI model code used to perform the simulations described can be obtained online at (Hager et al.2021).

Author contributions

AOH and MJH conceived of the study and designed the simulation plan. AOH conducted the model simulations, developed and carried out the analysis, wrote the majority of the manuscript, and created the figures. MJH also assisted with analysis and writing of the paper, designed and implemented the MALI subglacial hydrology model, and created the Thwaites model domain. SFP established funding for the research, provided experience in ice sheet modeling, and gave extensive guidance throughout the research process. DMS contributed his expertise in radar specularity analysis and interpretation, which was critical for comparison with model results. All authors contributed to editing the manuscript and discussing methodology.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank Trevor Hillebrand for processing the observational datasets used to generate the model initial condition and providing computing time on a Los Alamos National Laboratory Institutional Computing Program allocation. We thank Mauro Perego for contributing an optimized ice velocity and stress field for the model domain and for reviewing a version of the manuscript. We also thank our editor, Kang Yang, as well as Doug Brinkerhoff, Ankit Pramanik, and the anonymous reviewer, for their insights for improving the manuscript. This work developed out of a student project at the 2018 International Summer School in Glaciology in McCarthy, Alaska, organized by the University of Alaska Fairbanks.

Financial support

Support for this work was provided through the Scientific Discovery through Advanced Computing (SciDAC) program and the Energy Exascale Earth System Model (E3SM) project funded by the U.S. Department of Energy (DOE), Office of Science, Biological and Environmental Research, and Advanced Scientific Computing Research programs. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science user facility supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-05CH11231, and resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under contract DE-AC52-06NA25396. This work was also partially supported by the US National Science Foundation Office of Polar Programs under grant 1543012.

Review statement

This paper was edited by Kang Yang and reviewed by Douglas Brinkerhoff, Ankit Pramanik, and one anonymous referee.


Adusumilli, S., Fricker, H. A., Medley, B., Padman, L., and Siegfried, M. R.: Interannual variations in meltwater input to the Southern Ocean from Antarctic ice shelves, Nat. Geosci., 13, 616–620,, 2020. a, b, c

Alley, R.: Water-pressure coupling of sliding and bed deformation: I. Water system, J. Glaciol., 35, 108–118,, 1989. a

Alley, R. B.: Towards a hydrological model for computerized ice-sheet simulations, Hydrol. Process., 10, 649–660,<649::AID-HYP397>3.0.CO;2-1, 1996. a

Alley, R. B., Lawson, D. E., Evenson, E. B., Strasser, J. C., and Larson, G. J.: Glaciohydraulic supercooling: a freeze-on mechanism to create stratified, debris-rich basal ice: II. Theory, J. Glaciol., 44, 563–569,, 1998. a

Andresen, M. A.: Testing for similarity in area-based spatial patterns: A nonparametric Monte Carlo approach, Appl. Geogr., 29, 333–345,, 2009. a

Andresen, M. A.: An area-based nonparametric spatial point pattern test: The test, its applications, and the future, Methodological Innovations, 9, 2059799116630659,, 2016. a

Andrews, L. C., Catania, G. A., Hoffman, M. J., Gulley, J. D., Lüthi, M. P., Ryser, C., Hawley, R. L., and Neumann, T. A.: Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet, Nature, 514, 80–83,, 2014. a, b

Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone, R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.: Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), Geosci. Model Dev., 9, 2471–2497,, 2016. a, b

Asay-Davis, X. S., Jourdain, N. C., and Nakayama, Y.: Developments in simulating and parameterizing interactions between the Southern Ocean and the Antarctic ice sheet, Current Climate Change Reports, 3, 316–329,, 2017. a

Bevan, S. L., Luckman, A. J., Benn, D. I., Adusumilli, S., and Crawford, A.: Brief communication: Thwaites Glacier cavity evolution, The Cryosphere, 15, 3317–3328,, 2021. a, b

Bjornsson, H.: Jokulhlaups in Iceland: prediction, characteristics and simulation, Ann. Glaciol., 16, 95–106,, 1992. a, b

Brinkerhoff, D., Aschwanden, A., and Fahnestock, M.: Constraining subglacial processes from surface velocity observations using surrogate-based Bayesian inference, J. Glaciol., 67, 385–403,, 2021. a

Brinkerhoff, D. J., Meyer, C. R., Bueler, E., Truffer, M., and Bartholomaus, T. C.: Inversion of a glacier hydrology model, Ann. Glaciol., 57, 84–95,, 2016. a

Bueler, E. and van Pelt, W.: Mass-conserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635,, 2015. a

Carroll, D., Sutherland, D. A., Shroyer, E. L., Nash, J. D., Catania, G. A., and Stearns, L. A.: Subglacial discharge-driven renewal of tidewater glacier fjords, J. Geophys. Res.-Oceans, 122, 6611–6629,, 2017. a

Carter, S. and Fricker, H.: The supply of subglacial meltwater to the grounding line of the Siple Coast, West Antarctica, Ann. Glaciol., 53, 267–280,, 2012. a

Carter, S. P., Fricker, H. A., and Siegfried, M. R.: Antarctic subglacial lakes drain through sediment-floored canals: theory and model testing on real and idealized domains, The Cryosphere, 11, 381–405,, 2017. a, b

Chandler, D. M., Wadham, J. L., Lis, G. P., Cowton, T., Sole, A., Bartholomew, I., Telling, J., Nienow, P., Bagshaw, E. B., Mair, D., Vinen, S., and Hubbard, A.: Evolution of the subglacial drainage system beneath the Greenland Ice Sheet revealed by tracers, Nat. Geosci., 6, 195–198,, 2013. a

Chow, V. T.: Open-channel hydraulics, McGraw-Hill, New York, ISBN 9780070859067, 1959. a

Chu, W., Schroeder, D. M., Seroussi, H., Creyts, T. T., Palmer, S. J., and Bell, R. E.: Extensive winter subglacial water storage beneath the Greenland Ice Sheet, Geophys. Res. Lett., 43, 12–484,, 2016. a

Clarke, G. K.: Glacier outburst floods from “Hazard Lake”, Yukon Territory, and the problem of flood magnitude prediction, J. Glaciol., 28, 3–21,, 1982. a, b

Clarke, G. K.: Subglacial till: a physical framework for its properties and processes, J. Geophys. Res.-Sol. Ea., 92, 9023–9036,, 1987. a

Clarke, G. K.: Hydraulics of subglacial outburst floods: new insights from the Spring–Hutter formulation, J. Glaciol., 49, 299–313,, 2003. a, b

Cornford, S. L., Seroussi, H., Asay-Davis, X. S., Gudmundsson, G. H., Arthern, R., Borstad, C., Christmann, J., Dias dos Santos, T., Feldmann, J., Goldberg, D., Hoffman, M. J., Humbert, A., Kleiner, T., Leguy, G., Lipscomb, W. H., Merino, N., Durand, G., Morlighem, M., Pollard, D., Rückamp, M., Williams, C. R., and Yu, H.: Results of the third Marine Ice Sheet Model Intercomparison Project (MISMIP+), The Cryosphere, 14, 2283–2301,, 2020. a, b

de Fluerian, B., Werder, M. A., Beyer, S., Brinkerhoff, G. J., Delaney, I., Dow, C. F., Downs, J., Gagliardini, O., Hoffman, M., LeB Hooke, R., Seguinot, J., and Sommers, A. N.​​​​​​​: SHMIP The subglacial hydrology model intercomparison Project, J. Glaciol., 64, 897–916,, 2018. a, b, c, d

Dow, C., McCormack, F., Young, D., Greenbaum, J., Roberts, J., and Blankenship, D.: Totten Glacier subglacial hydrology determined from geophysics and modeling, Earth Planet. Sc. Lett., 531, 115961,, 2020. a, b, c, d, e, f, g, h, i

Dow, C. F., Kulessa, B., Rutt, I., Doyle, S., and Hubbard, A.: Upper bounds on subglacial channel development for interior regions of the Greenland ice sheet, J. Glaciol., 60, 1044–1052,, 2014. a

Downs, J. Z., Johnson, J. V., Harper, J. T., Meierbachtol, T., and Werder, M. A.: Dynamic hydraulic conductivity reconciles mismatch between modeled and observed winter subglacial water pressure, J. Geophys. Res.-Earth, 123, 818–836,, 2018. a, b, c, d

Drews, R., Pattyn, F., Hewitt, I. J., Ng, F. S. L., Berger, S., Matsuoka, K., Helm, V., Bergeot, N., Favier, L., and Neckel, N.: Actively evolving subglacial conduits and eskers initiate ice shelf channels at an Antarctic grounding line, Nat. Commun., 8, 15228,, 2017. a

Engelhardt, H. and Kamb, B.: Basal hydraulic system of a West Antarctic ice stream: constraints from borehole observations, J. Glaciol., 43, 207–230,, 1997. a, b

Flowers, G. E.: Modelling water flow under glaciers and ice sheets, P. Roy. Soc. A-Math. Phy., 471, 20140907,, 2015. a, b

Flowers, G. E. and Clarke, G. K.: A multicomponent coupled model of glacier hydrology 1. Theory and synthetic examples, J. Geophys. Res.-Sol. Ea., 107, 2287​​​​​​​,, 2002. a

Fountain, A. G. and Walder, J. S.: Water flow through temperate glaciers, Rev. Geophys., 36, 299–328,, 1998. a

Fowler, A.: Sliding with cavity formation, J. Glaciol., 33, 255–267,, 1987. a

Fowler, A.: Breaking the seal at Grímsvötn, Iceland, J. Glaciol., 45, 506–516,, 1999. a

Gardner, A. S., Moholdt, G., Scambos, T., Fahnstock, M., Ligtenberg, S., van den Broeke, M., and Nilsson, J.: Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years, The Cryosphere, 12, 521–547,, 2018. a

Gillet-Chaulet, F., Durand, G., Gagliardini, O., Mosbeux, C., Mouginot, J., Rémy, F., and Ritz, C.: Assimilation of surface velocities acquired between 1996 and 2010 to constrain the form of the basal friction law under Pine Island Glacier, Geophys. Res. Lett., 43, 10–311,, 2016. a, b, c, d, e

Gladish, C. V., Holland, D. M., Rosing-Asvid, A., Behrens, J. W., and Boje, J.: Oceanic boundary conditions for Jakobshavn Glacier. Part I: Variability and renewal of Ilulissat Icefjord waters, 2001–14, J. Phys. Oceanogr., 45, 3–32,, 2015. a

Gordon, S., Sharp, M., Hubbard, B., Smart, C., Ketterling, B., and Willis, I.: Seasonal reorganization of subglacial drainage inferred from measurements in boreholes, Hydrol. Process., 12, 105–133,<105::AID-HYP566>3.3.CO;2-R, 1998. a, b

Gulley, J., Walthard, P., Martin, J., a.F. Banwell, Benn, D., and Catania, G.: Conduit roughness and dye-trace breakthrough curves: why slow velocity and high dispersivity may not reflect flow in distributed systems, J. Glaciol., 58, 915–925,, 2012. a, b

Hager, A. O., Hoffman, M. J., Price, S. F., and Schroeder, D. M.: Persistent, Extensive Channelized Drainage Modeled Beneath Thwaites Glacier, West Antarctica, Zenodo [code and data set],, 2021. a

Haynes, M. S., Chapin, E., and Schroeder, D. M.: Geometric power fall-off in radar sounding, IEEE T. Geosci. Remote, 56, 6571–6585,, 2018. a

Helm, V., Humbert, A., and Miller, H.: Elevation and elevation change of Greenland and Antarctica derived from CryoSat-2, The Cryosphere, 8, 1539–1559,, 2014. a

Hewitt, I.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sc. Lett., 371, 16–25,, 2013. a, b

Hewitt, I. J.: Modelling distributed and channelized subglacial drainage: the spacing of channels, J. Glaciol., 57, 302–314,, 2011. a, b, c, d, e, f, g, h

Hoffman, M. J. and Price, S.: Feedbacks between coupled subglacial hydrology and glacier dynamics, J. Geophys. Res.-Earth, 119, 414–436​​​​​​​,, 2014. a, b

Hoffman, M. J., Andrews, L. C., Price, S. F., Catania, G. A., Neumann, T. A., Lüthi, M. P., Gulley, J., Ryser, C., Hawley, R. L., and Morriss, B.: Greenland subglacial drainage evolution regulated by weakly connected regions of the bed, Nat. Commun., 7, 13903​​​​​​​,, 2016. a, b, c, d, e

Hoffman, M. J., Perego, M., Price, S. F., Lipscomb, W. H., Zhang, T., Jacobsen, D., Tezaur, I., Salinger, A. G., Tuminaro, R., and Bertagna, L.: MPAS-Albany Land Ice (MALI): a variable-resolution ice sheet model for Earth system modeling using Voronoi grids, Geosci. Model Dev., 11, 3747–3780,, 2018. a, b, c

Hoffman, M. J., Asay‐Davis, X., Price, S. F., Fyke, J., and Perego, M.: Effect of Subshelf Melt Variability on Sea Level Rise Contribution From Thwaites Glacier, Antarctica, J. Geophys. Res.-Earth, 124, 2798–2822,, 2019. a, b

Hogan, K. A., Larter, R. D., Graham, A. G. C., Arthern, R., Kirkham, J. D., Totten Minzoni, R., Jordan, T. A., Clark, R., Fitzgerald, V., Wåhlin, A. K., Anderson, J. B., Hillenbrand, C.-D., Nitsche, F. O., Simkins, L., Smith, J. A., Gohl, K., Arndt, J. E., Hong, J., and Wellner, J.: Revealing the former bed of Thwaites Glacier using sea-floor bathymetry: implications for warm-water routing and bed controls on ice flow and buttressing, The Cryosphere, 14, 2883–2908,, 2020. a, b

Holschuh, N., Christianson, K., Paden, J., Alley, R., and Anandakrishnan, S.: Linking postglacial landscapes to glacier dynamics using swath radar at Thwaites Glacier, Antarctica, Geology, 48, 268–272,, 2020. a

Irarrazaval, I., Werder, M. A., Huss, M., Herman, F., and Mariethoz, G.: Determining the evolution of an alpine glacier drainage system by solving inverse problems, J. Glaciol., 67, 421–434,, 2021. a

Jenkins, A.: Convection-driven melting near the grounding lines of ice shelves and tidewater glaciers, J. Phys. Oceanogr., 41, 2279–2294,, 2011. a

Joughin, I., Tulaczyk, S., Bamber, J. L., Blankenship, D., Holt, J. W., Scambos, T., and Vaughan, D. G.: Basal conditions for Pine Island and Thwaites Glaciers, West Antarctica, determined using satellite and airborne data, J. Glaciol., 55, 245–257,, 2009. a, b, c, d, e, f, g, h

Joughin, I., Smith, B. E., and Medley, B.: Marine ice sheet collapse potentially underway for the Thwaites Glacier basin, West Antarctica, Science, 344, 735–738,, 2014. a, b

Joughin, I., Smith, B. E., and Schoof, C. G.: Regularized Coulomb Friction Laws for Ice Sheet Sliding: Application to Pine Island Glacier, Antarctica, Geophys. Res. Lett., 46, 4764–4771,, 2019. a, b, c, d

Kamb, B.: Glacier surge mechanism based on linked cavity configuration of the basal water conduit system, J. Geophys. Res.-Sol. Ea., 92, 9083–9100,, 1987. a, b, c

Koziol, C. P. and Arnold, N.: Incorporating modelled subglacial hydrology into inversions for basal drag, The Cryosphere, 11, 2783–2797,, 2017. a

Kyrke-Smith, T. and Fowler, A. C.: Subglacial swamps, P. Roy. Soc. A-Math. Phy., 470, 20140340,, 2014. a

Kyrke-Smith, T., Katz, R., and Fowler, A.: Subglacial hydrology and the formation of ice streams, P. Roy. Soc. A-Math. Phy., 470, 20130494,, 2014. a

Le Brocq, A. M., Payne, A., Siegert, M., and Alley, R.: A subglacial water-flow model for West Antarctica, J. Glaciol., 55, 879–888,, 2009. a

Le Brocq, A. M., Ross, N., Griggs, J. A., Bingham, R. G., Corr, H. F., Ferraccioli, F., Jenkins, A., Jordan, T. A., Payne, A. J., Rippin, D. M., et al.: Evidence from ice shelves for channelized meltwater flow beneath the Antarctic Ice Sheet, Nat. Geosci., 6, 945–948,, 2013. a, b, c

Leguy, G. R., Asay-Davis, X. S., and Lipscomb, W. H.: Parameterization of basal friction near grounding lines in a one-dimensional ice sheet model, The Cryosphere, 8, 1239–1259,, 2014. a, b

Livingstone, S. J., Clark, C. D., Woodward, J., and Kingslake, J.: Potential subglacial lake locations and meltwater drainage pathways beneath the Antarctic and Greenland ice sheets, The Cryosphere, 7, 1721–1740,, 2013. a

Magnusson, E., Rott, H., Bjornsson, H., and Palsson, F.: The impact of jokulhlaups on basal sliding observed by SAR interferometry on Vatnajokull, Iceland, J. Glaciol., 53, 232–240,, 2007. a

Marsh, O. J., Fricker, H. A., Siegfried, M. R., Christianson, K., Nicholls, K. W., Corr, H. F., and Catania, G.: High basal melting forming a channel at the grounding line of Ross Ice Shelf, Antarctica, Geophys. Res. Lett., 43, 250–255,, 2016. a, b

Martos, Y. M., Catalán, M., Jordan, T. A., Golynsky, A., Golynsky, D., Eagles, G., and Vaughan, D. G.: Heat Flux Distribution of Antarctica Unveiled, Geophys. Res. Lett., 44, 11417–11426,, 2017. a

Meierbachtol, T., Harper, J., and Humphrey, N.: Basal drainage system response to increasing surface melt on the Greenland ice sheet, Science, 341, 777–779,, 2013. a

Milillo, P., Rignot, E., Rizzoli, P., Scheuchl, B., Mouginot, J., Bueso-Bello, J., and Prats-Iraola, P.: Heterogeneous retreat and ice melt of Thwaites Glacier, West Antarctica, Science Advances, 5, eaau3433,, 2019. a, b, c, d

Morlighem, M., Rignot, E., Binder, T., Blankenship, D., Drews, R., Eagles, G., Eisen, O., Ferraccioli, F., Forsberg, R., Fretwell, P., Goel, V., Greenbaum, J. S., Gudmundsson, H., Guo, J., Helm, V., Hofstede, C., Howat, I., Humbert, A., Jokat, W., Karlsson, N. B., Lee, W. S., Matsuoka, K., Millan, R., Mouginot, J., Paden, J., Pattyn, F., Roberts, J., Rosier, S., Ruppel, A., Seroussi, H., Smith, E. C., Steinhage, D., Sun, B., den Broeke, M. R., Ommen, T. D., van Wessem, M., and Young, D. A.: Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nat. Geosci., 13, 132–137,, 2020. a

Mouginot, J., Rignot, E., and Scheuchl, B.: Sustained increase in ice discharge from the Amundsen Sea Embayment, West Antarctica, from 1973 to 2013, Geophys. Res. Lett., 41, 1576–1584,​​​​​​​, 2014. a

Murray, T. and Clarke, G. K.: Black-box modeling of the subglacial water system, J. Geophys. Res.-Sol. Ea., 100, 10231–10245,, 1995. a, b

Muto, A., Alley, R. B., Parizek, B. R., and Anandakrishnan, S.: Bed-type variability and till (dis)continuity beneath Thwaites Glacier, West Antarctica, Ann. Glaciol., 60, 82–90​​​​​​​,, 2019a. a, b, c

Muto, A., Anandakrishnan, S., Alley, R. B., Horgan, H. J., Parizek, B. R., Koellner, S., Christianson, K., and Holschuh, N.: Relating bed character and subglacial morphology using seismic data from Thwaites Glacier, West Antarctica, Earth Planet. Sc. Lett., 507, 199–206,, 2019b. a, b, c

Nakayama, Y., Manucharyan, G., Zhang, H., Dutrieux, P., Torres, H. S., Klein, P., Seroussi, H., Schodlok, M., Rignot, E., and Menemenlis, D.: Pathways of ocean heat towards Pine Island and Thwaites grounding lines, Scientific Reports, 9, 16649​​​​​​​,, 2019. a, b

Nakayama, Y., Cai, C., and Seroussi, H.: Impact of subglacial freshwater discharge on Pine Island Ice Shelf, Geophys. Res. Lett., 48, e2021GL093923,, 2021. a, b

Nias, I., Cornford, S., and Payne, A.: New mass-conserving bedrock topography for Pine Island Glacier impacts simulated decadal rates of mass loss, Geophys. Res. Lett., 45, 3173–3181,, 2018. a, b, c, d

Nye, J. F.: Water Flow in Glaciers: Jökulhlaups, Tunnels and Veins, J. Glaciol., 17, 181–207,, 1976. a, b

Perego, M., Price, S., and Stadler, G.: Optimal initial conditions for coupling ice sheet models to Earth system models, J. Geophys. Res.-Earth, 119, 1894–1917​​​​​​​,, 2014. a

Pritchard, H. D., Arthern, R. J., Vaughan, D. G., and Edwards, L. A.: Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets, Nature, 461, 971–975,, 2009. a

Rada, C. and Schoof, C.: Channelized, distributed, and disconnected: subglacial drainage under a valley glacier in the Yukon, The Cryosphere, 12, 2609–2636,, 2018. a, b, c, d

Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H., and Scheuchl, B.: Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011, Geophys. Res. Lett., 41, 3502–3509,, 2014. a, b, c

Rignot, E., Mouginot, J., Scheuchl, B., Van Den Broeke, M., Van Wessem, M. J., and Morlighem, M.: Four decades of Antarctic Ice Sheet mass balance from 1979–2017, P. Natl. Acad. Sci. USA, 116, 1095–1103,, 2019. a

Röthlisberger, H.: Water pressure in intra-and subglacial channels, J. Glaciol., 11, 177–203,, 1972. a, b

Schoof, C.: The effect of cavitation on glacier sliding, P. Roy. Soc. A-Math. Phy., 461, 609–627,, 2005. a

Schoof, C.: Ice-sheet acceleration driven by melt supply variability, Nature, 468, 803–806,, 2010. a, b, c, d

Schoof, C., Hewitt, I. J., and Werder, M. A.: Flotation and free surface flow in a model for subglacial drainage. Part 1. Distributed drainage, J. Fluid Mech., 702, 126–156,, 2012. a, b

Schroeder, D. M., Blankenship, D. D., and Young, D. A.: Evidence for a water system transition beneath Thwaites Glacier, West Antarctica, P. Natl. Acad. Sci. USA, 110, 12225–12228​​​​​​​,, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v

Schroeder, D. M., Blankenship, D. D., Raney, R. K., and Grima, C.: Estimating subglacial water geometry using radar bed echo specularity: application to Thwaites Glacier, West Antarctica, IEEE Geosci. Remote S., 12, 443–447,, 2015. a, b, c

Seroussi, H., Nakayama, Y., Larour, E., Menemenlis, D., Morlighem, M., Rignot, E., and Khazendar, A.: Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation, Geophys. Res. Lett., 44, 6191–6199​​​​​​​,, 2017. a, b

Slater, D., Nienow, P., Cowton, T., Goldberg, D., and Sole, A.: Effect of near-terminus subglacial hydrology on tidewater glacier submarine melt rates, Geophys. Res. Lett., 42, 2861–2868,, 2015.  a, b

Smith, B. E., Gourmelen, N., Huth, A., and Joughin, I.: Connected subglacial lake drainage beneath Thwaites Glacier, West Antarctica, The Cryosphere, 11, 451–467,, 2017. a, b, c, d, e, f

Stearns, L. A., Smith, B. E., and Hamilton, G. S.: Increased flow speed on a large East Antarctic outlet glacier caused by subglacial floods, Nat. Geosci., 1, 827–831,, 2008. a

Walder, J. S.: Hydraulics of subglacial cavities, J. Glaciol., 32, 439–445,, 1986. a, b

Walder, J. S. and Fowler, A.: Channelized subglacial drainage over a deformable bed, J. Glaciol., 40, 3–15,, 1994. a, b

Weertman, J.: General theory of water flow at the base of a glacier or ice sheet, Rev. Geophys., 10, 287–333,, 1972. a, b

Wei, W., Blankenship, D. D., Greenbaum, J. S., Gourmelen, N., Dow, C. F., Richter, T. G., Greene, C. A., Young, D. A., Lee, S., Kim, T.-W., Lee, W. S., and Assmann, K. M.: Getz Ice Shelf melt enhanced by freshwater discharge from beneath the West Antarctic Ice Sheet, The Cryosphere, 14, 1399–1408,, 2020. a, b, c, d

Werder, M. A.: The hydrology of subglacial overdeepenings: A new supercooling threshold formula, Geophys. Res. Lett., 43, 2045–2052,, 2016. a

Werder, M. A. and Funk, M.: Dye tracing a jokulhlaup: II. Testing a jokulhlaup model against flow speeds inferred from measurements, J. Glaciol., 55, 899–908,, 2009. a

Werder, M. A., Hewitt, I. J., Schoof, C. G., and Flowers, G. E.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res.-Earth, 118, 2140–2158,, 2013. a, b, c, d, e

Young, D., Schroeder, D., Blankenship, D., Kempf, S. D., and Quartini, E.: The distribution of basal water between Antarctic subglacial lakes from radar sounding, Philos. T. R. Soc. A, 374, 20140297,, 2016. a, b

Young, D. A., Roberts, J. L., Ritz, C., Frezzotti, M., Quartini, E., Cavitte, M. G. P., Tozer, C. R., Steinhage, D., Urbini, S., Corr, H. F. J., van Ommen, T., and Blankenship, D. D.: High-resolution boundary conditions of an old ice target near Dome C, Antarctica, The Cryosphere, 11, 1897–1911,, 2017. a

Yu, H., Rignot, E., Seroussi, H., and Morlighem, M.: Retreat of Thwaites Glacier, West Antarctica, over the next 100 years using various ice flow models, ice shelf melt scenarios and basal friction laws, The Cryosphere, 12, 3861–3876,, 2018. a, b, c, d, e

Zhao, K. X., Stewart, A. L., and McWilliams, J. C.: Geometric Constraints on Glacial Fjord–Shelf Exchange, J. Phys. Oceanogr. 51, 1223–1246, 2021. a

Short summary
The presence of water beneath glaciers is a control on glacier speed and ocean-caused melting, yet it has been unclear whether sizable volumes of water can exist beneath Antarctic glaciers or how this water may flow along the glacier bed. We use computer simulations, supported by observations, to show that enough water exists at the base of Thwaites Glacier, Antarctica, to form "rivers" beneath the glacier. These rivers likely moderate glacier speed and may influence its rate of retreat.