The past, present, and future evolution of Aurora Subglacial Basin’s subglacial drainage system

Domain: earth_science
Published: 2025-03-20
Editorial outcome: major_revision

Abstract
Aurora Subglacial Basin (ASB) in East Antarctica is among the most rapidly changing regions in Antarctica. Beneath the ASB, subglacial drainage networks allow for fast ice flow and facilitate mass losses by eroding the coastal ice shelves. The ASB underwent advance-retreat cycles since its inception 34 million years ago, and projections rely on these findings to predict future ice behaviour. Yet, these hindcasts and forecasts seldom consider the effects of subglacial drainage. Here, we model subglacial drainage networks and ice shelf basal melting of the ASB from 34 million years ago to 2100 CE. Our results indicate continual reorganisation of ASB drainage systems and ice shelf melt patterns, suggesting that the past does not fully serve as an analogue for the future. Since subglacial drainage accounts for up to 70% of ice shelf melt, our results imply that subglacial hydrology could trigger higher rates of mass loss than previously suggested.

Introduction
Ongoing climate warming has signalled concern about the fate of the Antarctic Ice Sheet and future sea levels, especially since the East Antarctic Ice Sheet contains 52 m of potential sea level rise1, the largest of all modern-day ice reservoirs. Within East Antarctica, Totten Glacier, in the Aurora Subglacial Basin (ASB, Fig. 1), is rapidly changing, displaying both the highest grounded ice thinning rates and highest rates of ice discharge of all East Antarctic outlet glaciers2. This area contains an ice volume equivalent of 3.9 m of sea level rise1 and is currently in a state of negative mass balance3. In recent decades, Totten Glacier has experienced periods of accelerated ice flow velocities, ice thinning, and grounding line retreat2–8. These changes are associated with intrusions of warm, modified Circumpolar Deep Water5,7–10 (mCDW), which erodes Totten Glacier’s Ice Shelf, the floating portion of Totten Glacier that serves as a girdle on ice flow.

Totten Glacier and the ASB may be inherently unstable due to a large proportion of the ice in the region grounded below sea level on a reverse bed slope1,11, making it susceptible to the Marine Ice Sheet Instability12 (MISI). Furthermore, collapse of buttressing ice shelves and unstoppable grounding line retreat may occur if the Marine Ice Cliff Instability (MICI) is set in motion due to ice shelf hydrofracture and ice cliff failure13. Models that include MICI suggest some of the highest projected ice loss and sea level rise estimates, but since there is limited observational evidence of past MICI behaviour14 and large uncertainties in MICI processes15, the plausibility of MICI continues to be debated. In the ASB, including MICI in modelled projections can cause thousands of square kilometres of grounded ice to be lost (Fig. S7). Regardless of the potential impacts of MICI, it is likely that mCDW will continue to affect ASB’s glaciers over the next century and push grounding lines farther inland, with the magnitude and extent of this retreat dependent on greenhouse gas emissions16. This situation could be exacerbated if instability mechanisms were to be initiated. Still, since most ice sheet models neglect evolving subglacial hydrology, it has seldom been considered that changes to subglacial drainage systems may be sufficient to initiate large-scale ice dynamic and geometric changes without the need for MICI.

The future trajectory of ice loss in ASB could potentially resemble that from the geologic past when climate conditions were similar to those predicted for the coming decades and centuries13,17. It has been suggested that the dynamic behaviour of the paleo-ASB can partially be explained by the subglacial drainage network configuration18–20. Present-day subglacial drainage in the ASB consists of a network of long, high-pressure channels that reach far into the interior of the ASB21. These channels play a critical role in regulating the flow of the overlying ice, while also pumping large volumes of freshwater beneath ice shelves, facilitating high and localised ice shelf basal melt rates22. If the volume of discharge exiting these channels, or their locations, were to change, there could be far-reaching impacts on the stability of ice flow in the ASB, such as accelerations of ice velocities and grounding line retreat rates23. However, there has not been a focused study on evolving subglacial drainage in the ASB during previous warm periods, how hydrology will be configured in future warm conditions, and the associated consequences of these changes for ice dynamics and ice shelf basal melt rates.

Here, we use the 2D, finite-element Glacier Drainage System24 (GlaDS, described in Methods) subglacial hydrology model to assess how ice sheet geometry at different time snapshots in the past and future affect the subglacial hydrology of the ASB’s grounded ice and infer the effects on ice shelf basal melt rates. We consider different concentrations of greenhouse gases and the impact of the MICI mechanism on the channel network and drainage system pressurisation. We show that there has been substantial reorganisation of the ASB’s basal hydrological networks since the inception of the Antarctic Ice Sheet 34 million years ago, and that drainage will continue to evolve over the next century in response to climate warming.

Results
Present-day subglacial hydrology
We first simulate modern-day subglacial hydrology using bed topography from BedMachine Antarctica1 and basal velocities and meltwater production rates from the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP625). With these datasets, we have revised and expanded the domain used in Dow et al.21 to capture up glacier activity and facilitate the comparison between our paleo and future simulations. Since Dow et al.’s21 model outputs were validated by radar-derived specularity content data, we consider their model outputs to be a benchmark for our model runs. We therefore adopt the same parameter values as used in Dow et al.’s21 base model run for bedrock obstacle sizes and distributed and channelised system conductivities. Like in Dow et al.21, we find that channel discharges are greatest through Totten Glacier, with fluxes reaching 33 m3 s-1 from the entire section of Totten’s grounding line (Fig. 2c). High water pressure occurs in the upland regions, another similarity between our simulations and those in Dow et al.21. Due to the high level of agreement between our modern-day results and those of Dow et al.’s21, the remainder of this paper focuses on how paleo and future subglacial hydrology compares to the contemporary drainage system. We direct readers to Dow et al.21 for an in-depth analysis of the present-day drainage system of the ASB.

Paleo subglacial hydrology
Our paleo subglacial hydrology simulations use grounding lines, ice sheet geometries, and basal inputs from Paxman et al.26. We focus on the Eocene-Oligocene Boundary (34 Mya) and the Mid-Miocene Climate Transition (14 Mya; hereafter Mid-Miocene) as these intervals correspond to ice volumes in ASB that are larger and smaller, respectively, than those in the present day26. The purpose of these experiments is to provide a large-scale overview of how changes in ice configurations and basal (melt, velocity, and topography) conditions influence ASB’s subglacial drainage system and are not intended to be interpreted as the exact subglacial conditions during each of these paleo time snapshots. In the Discussion, we compare our model results to geomorphological evidence.

Here, we focus on the subglacial hydrology and ice shelf basal melting results produced by adopting the ice geometries from Paxman et al.’s26 “No MICI, 500 ppm CO2” experiments. This level of CO2 represents (a) the likely lower bound during the Eocene-Oligocene Boundary, (b) the likely upper bound for the Mid-Miocene, and (c) the possible CO2 concentration for the future27–34. Due to limited observational evidence and divided findings from modelling efforts, is not yet clear if MICI has occurred in the geologic past and if it will occur in the next centuries14,15. We discuss how the ice geometries from Paxman et al.’s26 MICI and low CO2 experiments affect our results in the section titled “Carbon emissions and instability processes influence drainage morphology”.

Our Eocene-Oligocene Boundary results indicate that channels grew up to 370 km long (Fig. 2a), more than twice the length of present-day channels in the ASB (Fig. 2c and Dow et al.21). Nearly 160 m3 s-1 of channelised discharge exits the subglacial drainage system at the grounding line, shared between 9 outlet channels. 95% of this discharge is sourced from channels draining across Totten Glacier’s grounding line rather than the nearby Vanderford and Bond Glaciers. One of Totten’s channels pumps 117 m3 s-1 of freshwater into the ice shelf cavity, accounting for 75% of ASB’s cumulative grounding line discharge. Discharge from this channel alone exceeds the cumulative grounding line discharge from ASB’s current drainage system by 40 m3 s-1. All channels are fed by regions with water around ice overburden pressure. However, channels do not substantially reduce water pressure, and are within 5% of ice overburden pressure, consistent with modelling of present-day channels21. 95.47% of the domain is at or above 95% of ice overburden pressure; the highest proportion of all simulations.

We find a close correspondence with high pressure regions (Fig. 3) and thick (≥0.08 m) water sheets. These thicker sheets tend to be in basal troughs (Fig. S1) and show agreement with the up glacier initiation points of small channels. We do not observe an overlap between large water sheet thicknesses and high basal velocities or basal meltwater production rates. Since we do not couple our subglacial hydrology model to an ice sheet model, we use the basal meltwater production rates and basal velocities determined by Paxman et al.26 as inputs into our model. The highest velocities occur beneath Totten Glacier, which we capped at 800 m a-1 to prevent model instabilities (Table S1). Basal meltwater production rates beneath Totten Glacier reach 2 m a-1 at the grounding line, which is nearly six times larger than the present-day meltwater production rate at the grounding line21. Meltwater production remains elevated (≥10 cm a-1) up to 160 km from the grounding line, beyond which it approaches the catchment-wide mean of 6 mm a-1, exceeding the present-day catchment-wide average by 2.3 mm a-1. Overall, these results suggest spatially-variable controls on drainage beneath ASB. Topographic controls and basal melting have a greater influence on subglacial drainage up glacier, whereas high basal velocities and meltwater production rates play a dominant role in drainage morphology beneath Totten Glacier and closer to the grounding line. We performed sensitivity tests to determine the relative influence of bed and surface topography, basal melt, and basal velocity (Text S1). These tests indicated that the ice surface geometry is the primary control on where dominant drainage pathways form, but that changes in bed topography can influence channel lengths. Varying basal meltwater production rates and velocities caused more than a 10-fold change in grounding line channel fluxes and changes in the distribution and extent of high-pressure regions, both of which can affect ice dynamics through their influence ice shelf basal melt rates and on ice velocities, respectively.

In the Mid-Miocene, similar to the Eocene-Oligocene Boundary, high water pressures and thick water sheets are coincident and form in basal troughs. The largest channels also occur where basal velocities and rates of basal meltwater production are highest. Beyond this, the subglacial drainage morphology during the Mid-Miocene differs substantially from that at the Eocene-Oligocene Boundary. For instance, the maximum channel length for the Mid-Miocene is 270 km long (Fig. 2b), 100 km shorter than the longest channel in our Eocene-Oligocene runs (although ~100 km longer than present-day channels; Fig. 2c and Dow et al.21). The Mid-Miocene channels pump 35 m3 s-1 less water across the grounding line than the channels in the Eocene-Oligocene drainage system. However, there are more channels crossing the grounding line in the Mid-Miocene run since the ice in this simulation has retreated past Law Dome. This pattern of ice retreat results in a grounding line channel density (number of channels per km of grounding line) that is ~3.4 times higher than that of the Eocene-Oligocene Boundary. This suggests that retreat past Law Dome can trigger, or at least contribute to, reorganisation of water flow, as well as determining where the majority of the water exits the domain. We find that during the Mid-Miocene, discharge over the section of grounding line near Totten Glacier is relatively small compared to our Eocene-Oligocene simulations, and accounts for only a third of the total grounding line discharge. During the Mid-Miocene, our results indicate that the largest discharge crossing the grounding line occurs near the former location of Law Dome. Across this section of the grounding line, the maximum channel discharge is ~29 m3 s-1, with a cumulative channel discharge of 58.29 m3 s-1 in this region, which accounts for nearly half of the total channelised grounding line discharge.

Future subglacial hydrology
We next investigate the possible subglacial drainage morphologies of the ASB at 2100 CE, using ice geometries from Li et al.35, who considered how the spectrum of CMIP6 climate forcing impacts Antarctic Ice Sheet evolution. We selected the Li et al.35 ice geometries from CMIP6 forcings that resulted in the greatest ice loss (NorESM2-MM) and the largest ice gain (CMCC-CM2-SR5), representing end-member scenarios for ASB. Similar to our paleo experiments, we run two sets of simulations for each CMIP6 forcing, one in which the ice geometry was generated with MICI processes turned on, and one where these processes were excluded. The results from the MICI runs are presented in the supplementary information.

The degree of similarity of future subglacial drainage to present or past drainage networks is highly dependent on ice surface slope changes. For the maximum ice gain scenario, grounding line channel discharge at 2100 CE (Fig. 2e) is of a similar magnitude to the present (Fig. 2c, Table S3 and Dow et al.21) and has similar grounding line channel density. However, channels are up to 100 km longer in this future scenario. Under high rates of ice loss (Fig. 2d), channel length (303 km), discharges (163 m3 s-1), and grounding line channel density (0.0108 channels per km of groundling line) resemble that of the Eocene-Oligocene Boundary, albeit the organisation of the drainage system is more branching in the future simulation. These channel lengths, discharges, and densities are more than twice that of the present-day. For both future scenarios, most of the discharge exits at the grounding lines of Moscow University Ice Shelf and Totten Glaciers. However, while Totten Glacier accounted for three-quarters of total discharge during the Eocene-Oligocene Boundary, at 2100 CE, it only accounts for 35% of the grounding line discharge. This suggests that the changes in ice and bed geometry over time favour a more arborescent channelised drainage network. Furthermore, relative to the present-day simulation, we find that the partitioning of drainage in both of our future simulations tends to favour greater discharges beneath Bond and Vanderford Glaciers.

We find substantial differences between the channel discharges for the two future scenarios. The maximum ice loss scenario results in more than twice the channel density and discharge at the grounding line compared to the maximum ice gain scenario (0.0108 channels per km of groundling line for the maximum ice loss scenario versus 0.0053 channels per km of grounding line for the maximum ice gain scenario). The discharge velocity (channel flux divided by channel cross-sectional area) of the largest channel from the maximum ice loss scenario is 0.9741 m s-1, >30% faster than that observed in the maximum ice gain scenario (0.6347 m s-1) and >40 times faster than that from the Eocene-Oligocene (0.0233 m s-1). The steep surface slopes of the maximum ice loss run, especially near the grounding line, lead to higher hydraulic potential gradients (Fig. S13) that cause high basal velocities and meltwater production rates. Consequently, discharge volumes and associated water velocities are much larger than the others simulations, which drive enhanced rates of ice shelf basal melting at deep grounding line points by locally increasing rates of turbulent sub-ice shelf ocean mixing near the discharge outflow locations. We discuss these implications in the section titled “Ice shelf basal melt impacts”.

Similar to our paleo and present-day runs, our water pressure results for the future runs indicate that high pressures continue to persist in basal depressions. Compared to the paleo simulations, the region where water pressures reach or exceed ice overburden pressure (yellow contours in Fig. 3) has migrated further inland. These upland regions of over-pressured troughs do not occur in the present-day run. It is likely that erosion of the bed since the Eocene-Oligocene Boundary has heightened ASB’s sensitivity to climate-driven retreat26. Our simulations indicate that these bed changes over time have also influenced the distribution of subglacial water pressures and primary subglacial drainage pathways and support the findings of Wright et al.36 who suggested that a basal trough that meanders inland and near Totten Glacier connects up glacier subglacial drainage to the downstream drainage system. Our results show that pressures remain above overburden in the upland region, which, if using water pressure as a proxy for ice velocity, suggests that fast velocities may migrate further inland in the future.

Carbon emissions and instability processes influence drainage morphology
We also tested how paleo ice configurations that resulted from lower CO2 levels affected our outputs (Figs. S14–S21). Channels in our low (280 ppm) CO2 simulations were up to 37 km shorter than those from the 500 ppm simulations, with nearly 60 m3 s-1 less total grounding line discharge. Since the ice has not retreated as far back in the low CO2 simulations (Figs. S14–S15), the grounding line channel density is lower. This difference is most prominent in the Mid-Miocene 280 ppm No MICI run, where the grounding line channel density is half that of the 500 ppm simulation. Mean channel pressure in our low CO2 No MICI runs is also slightly lower (by up to 1.5%) than in our 500 ppm No MICI runs. However, the channel pressures are still within 95% of ice overburden, as in present-day Antarctic channels. This suggests that the effect of channels on ice-bed coupling in our paleo simulations is independent of the CO2 concentration. However, two-way coupled simulations of ice dynamics and hydrology are necessary to determine if this is the case. Similar to our 500 ppm simulations, thick and overpressured water sheets form in basal troughs. In the Eocene 280 ppm run, there is a lower proportion of the domain (54.7%) at and above 98% of ice overburden, compared to the 500 ppm run (58.48%). However, in the Mid-Miocene run, the opposite is true – raising the CO2 concentration lowers the percent of the domain at high pressure, from 49% in the 500 ppm run to 43.7% in the 280 ppm run.

The runs in which the ice geometries were influenced by MICI show considerable departures from the No MICI runs. In the MICI runs, grounding line channel densities were higher, channels tended to be shorter, and channel overburden pressures were lower than those in the No MICI runs. For the paleo runs, these results suggest that, in the absence of MICI processes, variability in CO2 concentrations can substantially affect subglacial drainage networks. This is likely because the steeper surface slopes in the higher-CO2 run cause steeper hydraulic potential gradients, resulting in longer channels. While MICI also tends to steepen ice surface, and therefore, hydraulic gradients, we do not observe longer channels in the paleo MICI runs compared to the No MICI runs because the grounding line has retreated much farther back in the runs where MICI is included. This pattern is not generalisable for our future simulations. While we do see shorter channels in the maximum ice loss MICI run, the maximum ice gain MICI simulation resulted in a maximum channel length of 308.73 km, more than 60 km longer than that of the corresponding No MICI run, and the longest channel of all the future runs. To investigate this further, we analysed the ratio of hydraulic potential gradients. In the MICI run, the ratio of hydraulic gradients for maximum ice loss to maximum ice gain is much more similar than in the No MICI runs. This is because the much steeper ice surface slopes in the maximum ice gain MICI run result in steeper hydraulic potential gradients. Relative to its respective No MICI simulation, the hydraulic gradients of the maximum ice gain MICI run are up to two times steeper. Therefore, we attribute the longer channels in the maximum ice gain MICI run to its steeper hydraulic gradients. Text S2 further discusses the interplay between ice geometry and hydraulic gradient.

We also found that the differences between total grounding line discharges in the No MICI versus MICI scenario were dependent on the time snapshot. In the paleo simulations, MICI lowered the total volume of grounding line discharge, with the largest decrease observed in the Mid-Miocene run, where discharges were 2.6 times lower in the MICI run compared to the No MICI run. In the future runs, however, incorporating MICI led to increases in channel discharge. While this increase was relatively small for the maximum ice loss run (3 m3 s-1), channel discharges in the maximum ice gain scenario increased nearly three-fold. The total channelised discharge volume in the maximum ice gain run with MICI run was the largest of all simulations considered, amounting to 198 m3 s-1. Taken together, our results highlight how MICI can substantially influence the drainage system. We motivate future work to build on what we have presented here by using coupled models that allow for evolving subglacial drainage and ice dynamics, as well as the option for MICI to be switched on and off. In doing so, a complete picture of how variations in subglacial drainage influence ice dynamics, and vice versa, will be obtained.

Ice shelf basal melt impacts
Subglacial outflow across Antarctic grounding lines and into sub-ice shelf ocean cavities has been posited to locally enhance ice shelf basal melting by both enhancing turbulence of the near-grounding line water mass and enhancing buoyancy of sub-ice shelf plumes that drive melting22,37,38. Using an ice shelf basal melt parameterisation that resolves melt driven by ocean circulation and subglacial freshwater outflow (Pelle et al.39; see Methods), we assess the sensitivity of ice shelf melt rates to the modelled paleo, modern, and future subglacial conditions discussed above. Parameterised ice shelf melt enhancements from subglacial freshwater outflow are strong but highly localised (Fig. 4). In the 500 ppm No MICI paleo simulations, the ice shelf melt enhancement at the Eocene-Oligocene Boundary (maximum enhancement of 38.43 m a-1) has a higher magnitude and is more spatially concentrated than the melt enhancement at the Mid-Miocene (maximum enhancement of 16.23 m a-1) due to the varying configurations of the subglacial channels draining across the grounding lines. In both paleo experiments, these melt enhancements comprise up to 70% of the total grounding line melt signal adjacent to the discharge channel (Fig. S32). In modern and future simulations, along Totten Glacier’s southern grounding zone, small changes in the subglacial channel distribution lead to considerable changes in the downstream limit of the ice shelf melt enhancement (5L’, see Methods), which reaches a maximum distance of 5.83 km, 10.04 km and 17.72 km in the modern, maximum ice gain, and maximum ice loss experiments, respectively. As ice shelf geometry and ocean conditions are consistent across these experiments, differences in the downstream limit of these melt enhancements are solely caused by differences in the configuration of the subglacial channels at the grounding line (i.e. fluxes and channel widths). In the modern, future maximum ice gain, and future maximum ice loss runs, subglacial outflow leads to enhanced ice shelf basal mass loss of 0.65 Gt, 2.23 Gt, and 3.62 Gt over the Totten Ice Shelf, respectively, and 0.21 Gt, 1.07 Gt, and 0.23 Gt over the Vanderford and Adams ice shelves, respectively (total basal mass loss over the Totten and Vanderford/Adams ice shelves in lieu of discharge are 178.44 Gt and 24.80 Gt, respectively). Although the magnitude of these shelf-integrated discharge-driven ice shelf melt enhancements seems inconsequential, these enhancements comprise up to 30% of the total Totten Ice Shelf grounding line melt signal in modern and maximum ice gain and up to 50% in maximum ice loss scenarios. Melt enhancements of Vincennes Bay outlet glaciers exhibit strong sensitivity to the simulated discharge conditions. Melt enhancements in the modern run are minimal due to the small amounts of simulated discharge (Fig. 4k); however, in 2100 (simulations with a more developed subglacial hydrology system) melt enhancements peak at 58.40 m a-1 and 25.64 m a-1 near the Vanderford and Adams grounding lines in maximum ice gain and maximum ice loss scenarios, respectively (Fig. 4l, m). We further explore the sensitivity of 2100 CE ice shelf melting of Vincennes Bay outlet glaciers (pink box in Fig. 4d) to subglacial hydrology model inputs of basal water production, ice velocity, and bed topography for the No MICI maximum ice gain runs (Fig. S33). Melt enhancements reach peak values near the southern limit of the Vanderford Ice Shelf and the Adams Ice Shelf in the modern bed topography experiment (Fig. S33a) and modern inputs experiment (Fig. S33c), respectively, owing to the differing distribution of grounding line discharge between these simulations. We note that in a complementary study, Pelle et al.23 developed a coupled ISSM-GlaDS configuration of ASB and, in projections to 2100, found that subglacial discharge-driven melt enhancements led to earlier onset of retreat of Totten Glacier’s southern grounding zone by ~25 years and continually accelerated retreat of Vanderford Glacier (up to 10 km of additional retreat compared to uncoupled “ISSM-only” simulations). That is, while the magnitude of these discharge melt enhancements is small relative to the background ocean-driven melt signal, their location at deep grounding line points has cascading impacts on the upstream stress regime of entire glacial basins. As such, the highly sensitive nature of the melt enhancements presented here highlights the importance of better constraining subglacial hydrology model inputs and boundary conditions, which will, in turn, improve simulations of past, present, and future ice sheet dynamics.

Methods
Subglacial hydrology modelling
For our subglacial hydrology simulations, we use the Glacier Drainage System (GlaDS) model24. GlaDS is a 2D finite-element model that allows for the co-evolution of distributed and channelised drainage over grounded ice. Distributed drainage is approximated as linked cavities that form on elements of a triangular mesh, which exchange water along element edges where channels form. Rather than being specified at the beginning of the run, these channels develop over the simulation time as fluxes and pressures evolve in the distributed system. The equation describing channel evolution is given by:

∂S/∂t = (Ξ/ρL) - (ASjNjn-1N) (1)

where S is the channel cross-sectional area, Ξ is the viscous dissipation of potential energy, Π describes the sensible heat changes due to changes in pressure melting point, ρi is the ice density, L is the latent heat of fusion, N is the effective pressure within the channel. A and n are empirically-derived constants, with the former relating to ice temperature, and the latter to and strain rate.

Although channels can form on any element edge within the domain, we restrict our analysis to channels with discharges >0.5 m3 s-1. To prevent model instabilities, we limit the maximum channel cross-sectional area to 500 m2. Close correspondence between GlaDS outputs (e.g., channel outlet locations and ponded water) and remotely sensed present-day basal conditions in East and West Antarctica suggests that the model is ideal for our investigation of paleo and future subglacial hydrology21,50–52.

We constructed 13 domains, one for each of the ice sheet configurations (8 for the paleo simulations, 1 for the present-day simulation, and 4 for the future simulations) considered in our study. These domains were determined by assuming that the subglacial water pressure is equal to ice overburden pressure. Since hydraulic potential analysis is based on the geometry of the ice and bed, the source data in this calculation influence the flow routing and, therefore, the domain size. From here, we construct our meshes, which are restricted to ~20,000 nodes due to computational limitations. We do not have refined regions in the meshes. In Table S1, we summarise key characteristics of the domains and the boundary conditions (e.g., the number of nodes, mesh edge lengths, maximum ice surface and bed elevation, basal water input rates, and basal velocities). We run the model up to steady-state over 30,000 days. Steady state is reached when the distributed water sheet thickness differs by <10-7 m day-1 between two consecutive timesteps21,50,53. After reaching steady state, we run the model for another 10,000 days. Results in the main text show the model outputs at the end of these 10,000 days. At this point, the sheet thickness differences between the two final timesteps range from 2.2 × 10-8 m day-1 to 8.52 × 10-7 m day-1 and are therefore within our criterion for reaching steady state.

GlaDS requires ice surface and basal topography, basal meltwater production rates (due to geothermal and frictional heating) and basal velocities as inputs. These data come from Paxman et al.26 for the paleo simulations, BedMachine Antarctica1 and initMIP-Antarctica25 for present day, and Li et al.35 for the future runs. We have summarised the range of values for these inputs in Table S1.

As mentioned in the main text, we selected paleo time snapshots when ice configurations in the ASB were larger (Eocene-Oligocene Boundary) and smaller (Mid-Miocene) than present day. These data are from Paxman et al.26, who explored the sensitivity of Antarctic Ice Sheet evolution to climate changes. Their experiments forced an ice sheet model with CO2 concentrations of either 280 ppm or 500 ppm, allowing MICI to be turned on or off. Using the Paxman et al. ice sheet model outputs as inputs to the GlaDS model, we consider all possible combinations for each time snapshot: 280 ppm without MICI, 280 ppm with MICI, 500 ppm without MICI, and 500 ppm with MICI. This resulted in a total of 8 simulations (4 for the Eocene-Oligocene Boundary, and 4 for the Mid-Miocene). Li et al.35 forced their ice sheet model with the 36 CMIP6 models available at the time of their publication. Since guidelines for selecting CMIP6 models for ice sheet modelling and analysis were unavailable at the time of writing, we select a subset of the model outputs from Li et al.35 based on a regional analysis of their performance over the ASB. Grounding line retreat distance was not an informative metric since the models showed the same grounding line locations for 2100, likely due to the coarse resolution of the models. Despite the similarities in grounding line locations, the patterns of ice loss and gain varied greatly between models. We, therefore, use the metric of total ice thickness difference at 2100 relative to 2020 over the ASB for the No MICI simulations. NorESM2-MM resulted in the largest ice loss, while CMCC-CM2-SR5 resulted in the greatest ice thickness gain over the ASB. These models therefore represent the endmembers for ice geometry change over the ASB based on Li et al.’s35 outputs. We also use the MICI complement of this subset of Li et al.35 model outputs, for a total of 4 future model runs.

The results presented in the main manuscript adopt the GlaDS parameter values that resulted in the closest agreement between present day hydrological modelling and remotely-sensed observations in the ASB21. These values are summarised in Table S2 in the supplement. In the Supplementary Material, we include the results of our sensitivity tests (Figs. S22–23 and Text S1). We focused on basal meltwater production rates, basal velocity, and basal topography. We ran simulations where we set the basal topography to that of BedMachine Antarctica1, but use the ice surface geometry, basal velocity, and basal melt data from the source publication. Since, in this test, we are only changing the bed topography but keeping everything else the same as in the source publication, it allows us to determine the level of control that basal topography has on the subglacial drainage system. Although the ice surface slope is known to exert a dominant influence on large-scale subglacial drainage40, water flow paths are also sensitive to basal topography54–56, with rerouting and water piracy possible under slight variations in bed topography54–56. Furthermore, large changes in bed slopes over relatively short distances can cause freeze-on as basal water moves from a region of high to low pressure; this effect has been incorporated in GlaDS57. Thus, our bed topography test allows us to examine the persistence and sensitivity of subglacial drainage pathways to bed geometry. Our second sensitivity test imposes modern-day values for basal topography, basal velocity and basal meltwater production to all model runs while keeping the ice surface geometry the same as in the source publication. This allows us to assess the role of ice surface topography on the drainage system evolution. Our final test examines the role of meltwater production rates by increasing the modern-day initMIP25 basal water production rates by a factor of 1.5, while keeping the basal velocities and ice and basal geometries the same as in the source publication. This provides an approximation of how melt could vary for ice sheet geometries governed by different ice sheet model physics.

Ice shelf basal melt parameterisation
For each subglacial hydrology simulation, we compute melting driven by sub-ice shelf ocean circulation and grounding line subglacial discharge with the buoyant plume parameterisation developed by Pelle et al.39. This parameterisation, based on the theory of buoyant sub-ice shelf plumes, applies a meltwater plume model58 that is modified near grounding line subglacial outflow locations to account for discharge-enhanced ice shelf melt using a discharge plume model37. The meltwater plume model simulates plumes that originate at the grounding line and travel upward along the base of the ice shelf, inducing melt up to the ice front or to the point at which the current loses buoyancy. The plume’s trajectory is controlled primarily by the ambient ocean temperature and salinity (Ta and Sa, respectively), the local slope angle of the ice shelf draft (α), and the depth the plume originated that reaches a particular shelf point (grounding line depth, zgl). Within a prescribed distance away from grounding line subglacial discharge locations, melt rates additionally become proportional to the grounding line subglacial freshwater flux (qsg) to the one-third power37,59. This distance (5L’, metres) that sets the cutoff between the meltwater and discharge plume models is set as a radius around grounding line discharge locations and is governed by the following equation:

5L’ = (500 * qsg^(2/3)) / (g * Δρi * C1 * Γ * Ts0 * G1 * G2 * λ3 * ΔTf,gl) (2)

where qsg^(2/3) is the grounding line subglacial discharge flux (m2 s-1), g is Earth’s gravitational constant (9.80 m s-2), Δρi is the density contrast between the plume and ambient ocean waters, C1=2Γ is a heat exchange parameter (6 * 10-4), G1 and G2 are geometric terms that provide melt rates the proper scaling to the slope angle of the ice shelf draft and grounding line depth field, and λ3 * ΔTf,gl is a measure of the ocean thermal forcing (°C). Given inputs of Ta, Sa, α, zgl, and qsg, the melt rate parameterisation outputs an ice shelf basal melt rate field that resolves melting driven by both large-scale ocean circulation and localised subglacial discharge. Note that melt enhancements from subglacial discharge were computed by differencing melt fields that were computed with qsg taken from GlaDS and with qsg set to zero.

Input ice shelf geometries (providing α and zgl) for both paleo, modern, and future subglacial hydrology simulations are consistent with the descriptions in the above subglacial hydrology modelling methods section. Input ocean conditions (Ta = -0.73 °C, Sa = 34.73 units on the practical salinity scale) for the 280 ppm paleo, modern, and future simulations were obtained by averaging observed60 bottom ocean temperature and salinity along the Sabrina Coast. For the 500 ppm paleo simulations, we apply 5 °C ocean warming to remain consistent with Paxman et al.26 and do not change the salinity (i.e., Ta = 4.27 °C and Sa = 34.73 units on the practical salinity scale).

Data availability
The subglacial hydrology model outputs and ice shelf basal melt rates generated in this study are available from the Zenodo repository: https://doi.org/10.5281/zenodo.14014169. Source data are provided with this paper.

Code availability
The Glacier Drainage System (GlaDS) model code is available from Mauro Werder (werder@vaw.baug.ethz.ch). GlaDS is also now implemented in the open-access Ice-sheet and Sea-level System Model (https://issm.jpl.nasa.gov/). The melt rate parameterization used in this study is published in Pelle et al.39 and uses the model outputs provided in the data submission.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of: The past, present, and future evolution of Aurora Subglacial Basin’s subglacial drainage system by Hayden et al\n– Nature Communications June 2024.\nOverall:\nThis paper uses a subglacial hydrology model to understand how the subglacial drainage system of the Aurora Subglacial\nBasin in East Antarctica would respond either in the past or in the future to changes in ice sheet configuration. A series of\nwarm time periods from the past are tested using input data from topographic reconstructions or modern topography and\nfrom a set of ice sheet models that have been used to understand warm-period response in this region. The model predicts\nwhere and how much subglacial discharge occurs, how it is routed to the edge of the continent and thereafter, using a plume\nmodel, how that water would translate into melt rates beneath the ice shelf. The outcomes of each model show quite\nsignificant differences in subglacial water flow pattern and rate, as well as differences in how the water reaches the margin\n(e.g. through a trunk vs. multiple outlets) which therefore has a knock-on impact on the pattern of melting under the ice\nshelf. A key point is that under retreated ice scenarios, the drainage system tends to become longer and produces more water,\ndelivering that to the grounding line where it could enhance ice shelf melting significantly. It is also useful that the authors\nhighlight that its not as simple as using a paleo warm-period scenario as an analogue for future warming impacts. The work\ndoes, as the authors say, represent the first time the hydraulic system has been investigated for this range of ice sheet\nchange scenarios and in that respect the key results are impactful. Indeed, it is an interesting paper and the quality of the\nwork is excellent. The work is largely valid in terms of predicting subglacial hydraulic changes, but as you can see with my\ntwo comments below, there are places where some additional consideration would help the manuscript. The approach and\ndata are appropriate, although there are also questions relating to some sensitivity tests (second comment below).\nSuggested Improvements:\nHowever, my first significant concern around the work is that it is not as impactful as it could be because the points about the\nimpact on ice velocity and buttressing which are used at the end of the paper to outline why this work is impactful are not\nactually tested. In other words, outputs of basal water conditions or ice shelf basal melting are not subsequently fed back into\nice sheet models in any way to demonstrate/quantify the importance. Thus although we can hypothesise that the impacts on\nice behaviour and mass change would be significant, this isn’t directly shown. I realise you signal that next steps would be to\ncouple things together in an ice sheet model, but there is a half-way house (via perturbation experiments) which would\nhammer home the significance of this work. Thus the impact of the changing hydrologic regime on the ice sheet feels\nsomewhat speculative at present because the impact on the ice sheet is not fully tested. If it could be firmed up in some way\nit would really demonstrate the significance of the results more clearly as well as becoming of interest to a broader\nreadership.\nSecondly, I felt that the work is not as clear as it could be about the relative influence of ice geometry vs. bed geometry than\nit could be (in the text its not so clearly justified why the bed geometry is less of a control than the ice surface geometry and\nthis may be a function of the fact that in sensitivity tests, there were varying basal topography whist holding ice sheet\ngeometries static for each scenario so I was less convinced that the bed topography influence had been clearly explored.\nThis is particularly important for the palaeo because ethe bed topography is quite different and quite smooth in these as\ncompared to the modern so some control for that needs to be tested.\nMany of the figures presented quite small maps and therefore there was significant challenge in seeing the variation\nbetween model outcomes – this would need addressed somehow so that patterns can be more clearly understood and\nlinked back to the text. I make specific comments on figures below.\nOther than the above points, I found a few minor issues in an otherwise interesting and well-written piece of work – very\ngood clarity and context and it uses literature appropriately.\nMinor points:\n41: mass imbalance – state in which direction (I assumed negative!).\nFig 1: text on figure is too small to read easily on colour bars. Red 50 m contour also hard to see in detail (can panels be\nmade a bit bigger?). Bed topograph also hard to make out in detail here (and we need that because its important for\nunderstanding the flow of basal water. Width of colourbar is too narrow to easily make out colour - expand width slightly.\nThere are red and green on same plot (e) - is this colourblind friendly? Guess the contour is dashed so it may be OK. I also\nwant to know what the ice surface actually looks like under these scenarios - its really not easy to tell here.\n60-62: On MICI, perhaps relating what projections with or without MICI would do in this area - e.g. is there much difference -\nexpress as uncertainty in ice extent and thus it becomes important to understand processes under both scenarios.\n76: Impacts on the stability of ice flow. Such as? Mention some impacts briefly.\nFig 2: Delineate grounding line to show difference between ice shelf and grounded ice? Does labelling need to be on all\npanels? Its not really clear where the ice shelf actually is - hence show grounding line as mentioned above. Can these\npanels be made bigger so we can see detail?\n119-123: I got a bit confused as to whether you were or were not modelling using MICI conditions. Your methods do indicate\nit, but in the text the water is muddied between what you say in line 119 and what you way on line 123.\n144: capped at 800 ma-1. What were the maximum values otherwise? I am just trying to see the extent to which you curtail\nthings here.\n149-150: It would be good to have a bit more info on the interplay between topo melt and velocity here - its a demonstration\nof what is causing the change and impacts of chance so seems odd to relegate it all to supplementary.\n195: ice configuration changes. by configuration, what do you mean? Extent? Thickness/surface? Flow? Which is more\nimportant as a control? Clarify please.\n221: This sentence sounds a bit speculative. Can you make a more concrete statement on this based on evidence?\n229-232: Were any sensitivity tests done whereby you use modern bed with palaeo ice surface/extent? In such a way it\nwould be possible to see what is more important for controlling the water conditions - the bed vs. the ice geometry. This\nfeeds into my second major comment above – it was less clear to me that you had controlled for the palaeo bed (which in\nitself is a model) and its impact on water flow.\n273-280: How much difference was there between palaeo and future MICI runs in terms of the ice geometry/thickness?\nCould be good to have some difference maps in supplement to highlight this and understand what it is about geometry that\nis controlling the hydraulic gradient.\nFig 4 (title): is it enhancements, or is it simply melt rates (from which we can then interpret enhancements)? I.e. the colour\nscale is about melt, so the figure is about melt, not melt enhancements. Alter the title accordingly.\n314 (and this paragraph): for melt and how it changes in each simulation, would it be worth stating how much the total melt\nwithin each system would be? Or the average beneath the ice shelf? That way we can understand more clearly not only the\nfocussing (or not) but also whether more melt happens overall or not.\n334: can the manner in which these melt enhancements impact glacier retreat and mass loss be expanded on a bit? This is\na key implication so slightly more detail is needed. It links to my first significant comment above.\n356: ice flow is likely sensitive. feels like speculation - e.g. not specifically shown in this paper.\n386: water systems facilitate ice motion – I agree but its not directly shown here. For a more impactful paper, this would need\nto be shown somehow, or at least expanded to really demonstrate what the impacts are in this case.\n391-392: weaker buttressing also not shown here - for more impact, this would need to be demonstrated.\n440: basal meltwater production rates. these can be quite different depending on ice sheet model used. Does the potential\nvariation introduced by ice sheet model differences need to be incorporated?\n472-472: Why choose to use ice thickness as opposed to ice surface geometry? If you use ice thickness, you would surely\nintroduce some quite different ice geometries which may not be realistic representations of those times. Perhaps a altering\nice surface would have been more appropriate? In addition, were tests done where bed was varied but other things were\nnot? That would help get at whether the bed is relatively more or less important than the geometry. Perhaps expanding the\nlogic here would help with response to my second major comment above.\n479-480: I am struggling to follow the logic here (sorry!) - I haven't seen you vary topography whilst keeping ice surface the\nsame - wouldn't that also be important to do?", "reviewer_2": "I read this well-prepared and clearly written manuscript with great interest because there is compelling geomorphologic\nevidence to suggest that subglacial meltwater, either from the drainage of subglacial lakes or surface melting, may be\ninfluencing/may have influenced ice dynamics in the ASB throughout the Cenozoic. As such, modeling studies are required\nto assess the potential impact of subglacial meltwater presence on ice behavior during past warm periods, particularly in the\nTotten Glacier catchment, where thinning and ice mass loss and an active subglacial hydrologic system are presently\nobserved.\nTwo past time intervals were chosen as part of this study to reflect times when 1) regional ice volume was larger than at\npresent with higher pCO2 (E-O boundary) and 2) possibly smaller than present with who knows what pCO2 (Middle\nMiocene). At least I think this was the goal because the introduction talks about needing to understand how the subglacial\nhydrologic system worked during past warm periods. However, the intervals listed are those of glacial expansion. Climate at\nthe EO boundary was relatively cold, and if the goal of the middle Miocene study was to hit the Miocene Climate Optimum,\nthe warmest time interval of the last 30 Ma, that date and ice sheet configuration should reflect ~16 Ma and not ~14 Ma,\nwhen ice expanded and stabilized. The authors conclude that there have been substantial reorganizations of the subglacial\ndrainage system over ~40 million years. Well yes, this seems like a pretty safe takeaway for a manuscript.\nAs an observational scientist, I understand that models do not often have the resolution that observers want and that slight\nchanges in boundary conditions dictated by observations may not make large differences in model results. With this in mind,\nI have a few questions that probably need clarification/addressing in the text: Why were the Paxman et al boundary\nconditions chosen for this study, rather than the more regionally specific Aitken et al 2016 geomorphologic conditions?\nFurther, how do the model results presented in this study compare to those of Aitken, which comment on basal velocities\nand channel systems visible in subglacial images, but do not prescribe a specific hydrologic system, rather wet and cold-\nbased regions based on geomorphology?\nPast grounding line configurations: The choice of grounding line configurations used in the paleo model runs is not clear.\nThis may be a place where model resolution is important, but because there is geological data to suggest that the chosen\npaleo grounding zone choices are incorrect, the choice of paleo grounding line location warrants further discussion for each\ntime interval. For instance, is the resolution of the model in this manuscript such that an adjustment of the grounding line out\nonto the continental shelf could make a difference to the model output? It might be that the model cannot resolve changes in\nthe grounding line of 100s of Km, and that is fair and could be mentioned. I think this is what the authors say in the methods\nand supplement, but that should be made clear in the main text. Further, if this is the case, how realistic are these model\nresults, given the complex and relatively fine-scale subglacial topography in the region?\nFor example, in Aitken’s and other studies, it is clear that the topographic high termed Ridge A might be essential for\nstabilizing regional grounding lines—past and present. Given the position of Ridge A at the mouth of the ASB, which means\nthat there is no reverse sloping bed until ice retreats landward of Ridge A, it might be useful for the authors to both address\nthis feature and the evidence for ice advance onto the shelf during early ice sheet evolution. Additionally, consideration of\nRidge A is important because (based on the channels in the bedrock on the seaward side of Ridge A; Fernandez et al.,\n2018), subglacial meltwater may be really important for getting ice over that high and out onto the shelf.\nIn the manuscript, it appears that Ridge A is an important constraint on the location of the EO grounding line. Is this a correct\nstatement? If it is, it should be noted that data from the continental shelf presented by Gulick, Shevenell, et al. (2017; note\nthat this is the correct citation for this paper) indicate grounded ice expansion onto the shelf, possibly at the EO boundary,\nbased on existing chronological constraints (not great, but data that exists until the shelf sediments are drilled). This\nobservation suggests that grounded ice had to make it over Ridge A and onto the shelf as ice in the ASB expanded. The\nauthors should mention the shelf evidence when discussing their choice of grounding line configuration for the EO\nboundary.\nThe authors should also address how Middle Miocene grounding zone locations were determined especially in light of data\nfrom the continental shelf that indicates at least 10 advances and retreats onto the shelf between what is likely the E-O\nboundary and the Late Miocene (Gulick, Shevenell et al., 2017). While the ages of these unconformities are unconstrained,\nthey are significant in that they also indicate large volumes of meltwater on the shelf based on the presence of tunnel valley\nfeatures. It seems that there is an unprecedented opportunity both address apparent data/model disagreement and compare\nmodel results with what we know about subglacial hydrology in the region from the geologic record.\nPast sea level: there is substantial data to suggest significant changes in sea level in the past. Was past sea level\nconsidered in the modeling exercise? How might past sea level impact the choice of grounding line locations in the paleo\nmodel runs? We do know that sea level dropped to -30 to -40 m relative to modern day across the EO boundary (Miller et al.,\n2020), which likely enabled grounded ice expansion from the ASB across the continental shelf, as observed in the geologic\nrecord. Further, during the Miocene Climate Optimum (~16 Ma), sea levels may have been +50m higher than present (likely\noverestimate), while at ~14 Ma, sea levels may have been -40 m lower than present. This is another reason why the authors\nreally need to clarify if they are talking about the MCO (~16 Ma) or ice expansion and stabilization associated with the\nMiddle Miocene Climate Transition at ~14 Ma.\nPaleo channel locations: there are paleo channel systems identified both in the subglacial landscapes (see Young et al.,\n2011, Aitken et al., 2016) inland of Ridge A and on the continental shelf seaward of the current study’s boundaries (see\nGulick, Shevenell et al., 2017; Fernandez et al., 2018; Montelli et al., 2020). How do the model results compare with the\nlocations of these channels? Some discussion of the modelled channel locations vs observed channel locations would\nstrengthen the discussion of model results and the impact of the manuscript.\nMICI as a viable mass loss mechanism: the authors admit the tenuous nature of the MICI mechanism in ice sheet models.\nThis reviewer is of the opinion that inclusion of MICI in the runs and in any discussion diminishes the results, as presented\n(mostly because the reader is left being skeptical of MICI, rather than thinking about the present modelling study). Would it\nbe better to think about and address how substantial subglacial meltwater presence might eliminate the need for a MICI-type\nmechanism to generate a responsive ice sheet? Such an argument, if possible, would really highlight the need to\nunderstand subglacial hydrologic systems, past and present.\nGenerally, after reading this manuscript, it seems that the authors have missed an opportunity to compare model runs with\nexisting observations of subglacial hydrologic systems in this catchment. This is one of the only locations in Antarctica\nwhere such a comparison is possible. Such an exercise would significantly strengthen the believability and impact of the\nmodeling studies. I recommend such a data-model comparison be undertaken before publication is considered.\nA few more minor questions/comments on the manuscript, in no particular order:\nThis study currently uses outdated CO2 reconstructions. I recommend using the latest CO2 compilation (CENCO2PIP\nConsortium, 2023).\nFigure 1 and its caption are quite confusing (there might be a typo for the E-O simulation- should it be C and H panels?). It\nwould be helpful to list the time intervals above the model domains and basal input as in Figure 4, and maybe place vertical\nlines between the time intervals.\nWhy does the Law Dome region disappear in the Mid-Miocene simulation? Is this related to sea level? And, if it is, the same\nwould be true for the Eocene, would it not?\nI don't understand the discussion of ice shelves in past ice configurations. There is no evidence of ice shelf presence in the\nEocene or the middle Miocene, but rather grounded ice on the continental shelves. How does this fact impact the results of\nthis study?\nFernandez et al (2018) present pretty interesting data for meltwater input to the shelf region during the last deglaciation (see\nbedrock channels and stacked grounding zone wedges termed feather moraines in the publication). Would this be a useful\ntime interval over which to test the model, considering that we have knowledge about the location of the ice sheet and of\nsubglacial meltwater output/presence during the last deglaciation just based on continental shelf geomorphology.\nWhat happens if surface-derived meltwater is included for the MCO run?", "reviewer_3": "Review of Hayden et al.\nThis is a very well written paper on changing Aurora subglacial drainage at different climate states that is based on a\nfinite-element glacier drainage system model which responds to climate forcing through a link between ice thickness and\nwater pressure at the bed. Using this methodology forced today with BedMachine for the topography and ISMIP6 for the\nclimate they simulate potential drainage pattern under Totten which seems reasonable.\nThey then attempt similar situations at the Eocene-Oligocene boundary and at the mid-Miocene. For the latter I think they are\nseeking to simulate the Mid-Miocene Climatic Optimum. There are issues with specifically matching the O-18 curve with\nCO2 as the authors are undoubtedly aware, but I don’t think this matters for the paper since you are forcing it at two CO2\nstates roughly equivalent to E-O and MMCO. The specific date does not matter for your study since as noted below, I don’t\nthink you are yet able to simulate what specifically occurred with the ice sheet and drainage system at the E-O and MMCO.\nThe larger issues however are three:\n1) bed topography has certainly evolved from the Eocene to today and thus what is the sensitivity of the model to this\nprocess since we do not have a E-O or Mid-Mio bedmachine? There are some sentences discussing how it might have\nchanged but I was left unclear as the expected magnitude and feedbacks if say the deeper basins under ASB we not yet\nfilled in and we had a subglacier lake scenario, or if the erosion through time has substantially changed the overall gradient\nto flatter today than in the E-O run. Can we be confident of these glacial drainage change results in absence of the\ntopography? Some sensitivity study/discussion is needed to place these results in context on this issue.\n2) The ice margin for the ASB portion of East Antarctica is not substantially changing between the E-O, Mid-Mio, Today, and\nFuture scenarios. There is certainly some retreat and some advance aspects but these are way too small to match the\noffshore shelf and slope geophysical data implying presence or absence of grounded ice (Gulick et al., 2017, Nature,\nMontelli et al. 2019 GSA Bulletin, McKay et al., 2022, Ant Clim Evol). For instance, at the E-O boundary there were tunnel\nvalleys > 125 km seaward of the current ice margin which require grounded ice and we have a glacial shelf crossing trough\n(Donda et al., 2020, Ant Sci; Hirano et al., 2023, Nature Comm) reaching the shelf edge during the advances. There is strong\ndata for grounded ice therefore over much of the shelf and within the trough to the shelf edge since there is a large trough\nmouth fan for the Totten system active in the Oligocene at least. Therefore your E-O simulation ice margin should be\nsubstantially more advanced to match this geologic evidence. In counterpoint, there appears to have been no evidence of\nproximal glacial sedimentary products in the MMCO on the shelf and thus one would infer a much more retreated ice margin\nto explain this. Your Mid-Mio simulation is somewhat retreated but probably far too little to explain the absence of glacial\nsedimentary facies at the MMCO on the shelf.\n3) The lack of greater changes in ice margin at warmer and colder periods imply the current models are missing some key\nfeedbacks. It could be the effect of ice cliff instability as raised by the authors. However another key change that should\noccur is going from colder to warmer based glaciers such as is seen in north versus south Greenland outlet glaciers where\nwarmer wetter ones moving larger volumes of sediment potentially building glacial moraines or morainal banks, and are\ncorrelated with tunnel valleys existing when advancing onto shelves (e.g., Alaska, North Sea). Colder states act more like\nmodern Antarctica or the colder based parts of Greenland forming grounding zone wedges instead and substantially less\nfeedbacks are likely. As Totten system transitioned to this more polar state only the Late Miocene-Pliocene, it is probably\nincorrect to assume the subglacial system worked the same way from a more modern polar state as it did in a past more\ndynamic state. In other words E-O and Mid-Mio likely dynamically behaved more like the North Sea or Alaska systems than\nPleistocene Antarctica in order to explain the greatly different shelf margin effects and glacial sedimentary products. What\nexactly is the change in feedbacks is that allows for much more dynamic advance and retreat cycles is not yet known of\ncourse but some combination of melt water, warm bases, and glacial sediment dynamics are likely involved.\nMoving forward with this paper I think the authors need to specifically note how the paleo model runs do not produce\ndynamic enough advances to reach the shelf edge during the E-O or leave the shelf entirely ice free during the Mid-Miocene\nand then discuss the possible additional feedbacks needed to result in greater ice margin changes. Then discuss what effect\nthis would have on the subglacial networks that are the focus of this paper. In summary, this paper does a good job with\nthe paleo models if the bed topography did not differ substantially and if the glacial dynamics did not change going into past\nclimate states. Issue is we know the ice was a lot further advanced and have evidence of a wetter more dynamic system that\ncreated tunnel valleys, trough mouth fans, and probably intermittent morainal banks … how is this kind of system, one more\nrepresentative of the LGM Alaska, North Sea, or West Greenland for instance, differ from what is modeled in this paper?\nGiven there is clearly a mismatch in ice edge position with the geologic record at least for maximum advance and retreat\nclimate states, then how can we use this paper’s results to inform us of process without yet being able to include all the\nprocesses required? I think the path forward is to use the paleo model results to discuss how different inputs such as CO2\naffect ice thickness and water pressure and drainage without making the argument that this study is faithfully\nreproducing the conditions at E-O and Mid-Mio. Consider what you models are not including and how you might then frame your model\nresults differently since its interesting the tradeoffs observed in water pressure and drainage length in your work, and\nconsidering future scenarios (see below), but these models clearly do not represent the state of the ASB glacial system at\nthe E-O or mid-Mio.\nFor the future model runs, I think these are fine and in fact really interesting suggestions for future\nstates, if the system stays in its current polar mode of dynamics. However, if it switches modes to the way it operated in the\nOligocene to Miocene then there would be expected a much greater dynamic change possibility that could be acknowledged.\nOn the writing, this paper is remarkably free of typos, has excellently written sentences, and was very clear to read. Only\nminor item on that front is that the word further (a concept) is used in two places where the word should be farther (a\ndistance).\nMy recommendation is therefore Major Revisions to consider the issues with mismatch to the geologic record and\nsuggestions of glacial dynamics having a warmer wetter regime with processes not yet included in the finite glacial drainage\nsystem modeling.\nHopefully this review will be helpful and I would be willing to do a re-review as needed,\nSean Gulick"}, "author_response": "Reply to reviewers\nThe past, present, and future evolution of Aurora Subglacial Basin’s subglacial drainage system\nHayden, Pelle, and Dow\nWe would like to thank the reviewers for their comments. We have addressed the comments point-\nfor-point below. The comments are laid out in normal font, and our replies are in italic green font.\nReviewer 1:\nThis paper uses a subglacial hydrology model to understand how the subglacial drainage system\nof the Aurora Subglacial Basin in East Antarctica would respond either in the past or in the future\nto changes in ice sheet configuration. A series of warm time periods from the past are tested using\ninput data from topographic reconstructions or modern topography and from a set of ice sheet\nmodels that have been used to understand warm-period response in this region. The model predicts\nwhere and how much subglacial discharge occurs, how it is routed to the edge of the continent and\nthereafter, using a plume model, how that water would translate into melt rates beneath the ice\nshelf. The outcomes of each model show quite significant differences in subglacial water flow\npattern and rate, as well as differences in how the water reaches the margin (e.g. through a trunk\nvs. multiple outlets) which therefore has a knock-on impact on the pattern of melting under the ice\nshelf. A key point is that under retreated ice scenarios, the drainage system tends to become longer\nand produces more water, delivering that to the grounding line where it could enhance ice shelf\nmelting significantly. It is also useful that the authors highlight that its not as simple as using a\npaleo warm-period scenario as an analogue for future warming impacts. The work does, as the\nauthors say, represent the first time the hydraulic system has been investigated for this range of ice\nsheet change scenarios and in that respect the key results are impactful. Indeed, it is an interesting\npaper and the quality of the work is excellent. The work is largely valid in terms of predicting\nsubglacial hydraulic changes, but as you can see with my two comments below, there are places\nwhere some additional consideration would help the manuscript. The approach and data are\nappropriate, although there are also questions relating to some sensitivity tests (second comment\nbelow).\nSuggested Improvements:\nHowever, my first significant concern around the work is that it is not as impactful as it could be\nbecause the points about the impact on ice velocity and buttressing which are used at the end of\nthe paper to outline why this work is impactful are not actually tested. In other words, outputs of\nbasal water conditions or ice shelf basal melting are not subsequently fed back into ice sheet\nmodels in any way to demonstrate/quantify the importance. Thus although we can hypothesise that\nthe impacts on ice behaviour and mass change would be significant, this isn’t directly shown. I\nrealise you signal that next steps would be to couple things together in an ice sheet model, but\nthere is a half-way house (via perturbation experiments) which would hammer home the\nsignificance of this work. Thus the impact of the changing hydrologic regime on the ice sheet feels\nsomewhat speculative at present because the impact on the ice sheet is not fully tested. If it could\nbe firmed up in some way it would really demonstrate the significance of the results more clearly\nas well as becoming of interest to a broader readership.\nThank you for raising this important point and for your help in maximising the impact of our\nresults. During the time this manuscript was in review, coauthor Pelle and colleagues published\ncomplementary work that presents future 2100 projections of the Aurora Subglacial Basin in a\ncoupled ice-subglacial hydrology model (ISSM-GlaDS) to test the sensitivity of outlet glacier\nbehaviour to future subglacial and ocean forcing (Pelle et al., 2024). In their projections, they\nfound that ice-subglacial water feedbacks (ice shelf melt enhancements and increased grounded\nice basal sliding) increased the ASB sea level contribution by 30%, drove earlier onset of retreat\nof Totten’s southern grounding zone by ~25 years, and led to continual accelerated retreat of\nVanderford and Adams Glaciers over the full 83-year projection. These new results confirm our\nhypothesis of the strong impact that subglacial processes can have on outlet glacier dynamics.\nThus, in the revised manuscript, we have framed our results in the context of these new findings\n(L360-370):\n“We note that in a complementary study, Pelle et al.24 developed a coupled ISSM-GlaDS\nconfiguration of ASB and, in projections to 2100, found that subglacial discharge-driven melt\nenhancements led to earlier onset of retreat of Totten Glacier’s southern grounding zone by ~25\nyears and continually accelerated retreat of Vanderford Glacier (up to 10 km of additional retreat\ncompared to uncoupled “ISSM-only” simulations). That is, while the magnitude of these discharge\nmelt enhancements is small relative to the background ocean-driven melt signal, their location at\ndeep grounding line points has cascading impacts on the upstream stress regime of entire glacial\nbasins. As such, the highly sensitive nature of the melt enhancements presented here highlights the\nimportance of better constraining subglacial hydrology model inputs and boundary conditions,\nwhich will, in turn, improve simulations of past, present, and future ice sheet dynamics.”\nWe note that Pelle et al. only ran two subglacial scenarios (a low and high emission climate\nforcing) due to the computational cost of running the coupled models through 2100 and thus did\nnot test the sensitivity of the system to variations in subglacial water input and changes in ice/bed\ngeometry as we have explored in our manuscript. They also did not perform paleo-simulations.\nNonetheless, we argue that the significance of our work is that we simulate paleo, modern, and\nfuture ASB subglacial hydrology in a consistent modelling framework and can highlight model\nsensitivity to the various inputs.\nPelle, T., Greenbaum, J. S., Ehrenfeucht, S., Dow, C. F., & McCormack, F. S. (2024). Subglacial\ndischarge accelerates dynamic retreat of aurora subglacial basin outlet glaciers, East Antarctica,\nover the 21st century. Journal of Geophysical Research: Earth Surface, 129, e2023JF007513.\nhttps://doi.org/10.1029/2023JF007513\nSecondly, I felt that the work is not as clear as it could be about the relative influence of ice\ngeometry vs. bed geometry than it could be (in the text its not so clearly justified why the bed\ngeometry is less of a control than the ice surface geometry and this may be a function of the fact\nthat in sensitivity tests, there were varying basal topography whist holding ice sheet geometries\nstatic for each scenario so I was less convinced that the bed topography influence had been clearly\nexplored. This is particularly important for the palaeo because ethe bed topography is quite\ndifferent and quite smooth in these as compared to the modern so some control for that needs to\nbe tested.\nWe thank the reviewer for raising this point. We have added clarification about the sensitivity tests\nwe perform in the main text (L169-175), methods (L541-561), and in our responses to the\ncomments below from reviewers 1 and 3 about our tests of ice and bed geometry and basal\nconditions.\nMAIN TEXT (L169-175):\n“We performed sensitivity tests to determine the relative influence of bed and surface topography,\nbasal melt, and basal velocity (Text S1). These tests indicated that the ice surface geometry is the\nprimary control on where dominant drainage pathways form, but that changes in bed topography\ncan influence channel lengths. Varying basal meltwater production rates and velocities caused\nmore than a 10-fold change in grounding line channel fluxes and changes in the distribution and\nextent of high-pressure regions, both of which can affect ice dynamics through their influence ice\nshelf basal melt rates and on ice velocities, respectively.”\nMETHODS (L544-564)\n“In the Supplementary Material, we include the results of our sensitivity tests (Figs. S22-23 and\nText S1). We focused on basal meltwater production rates, basal velocity, and basal topography.\nWe run simulations where we set the basal topography to that of BedMachine Antarctica1, but use\nthe ice surface geometry, basal velocity, and basal melt data from the source publication. Since,\nin this test, we are only changing the bed topography but keeping everything else the same as in\nthe source publication, it allows us to determine the level of control that basal topography has on\nthe subglacial drainage system. Although the ice surface slope is known to exert a dominant\ninfluence on large-scale subglacial drainage40, water flow paths are also sensitive to basal\ntopography54–56, with rerouting and water piracy possible under slight variations in bed\ntopography54–56. Furthermore, large changes in bed slopes over relatively short distances can\ncause freeze-on as basal water moves from a region of high to low pressure; this effect has been\nincorporated in GlaDS57. Thus, our bed topography test allows us to examine the persistence and\nsensitivity of subglacial drainage pathways to bed geometry. Our second sensitivity test imposes\nmodern-day values for basal topography, basal velocity and basal meltwater production to all\nmodel runs while keeping the ice surface geometry the same as in the source publication. This\nallows us to assess the role of ice surface topography on the drainage system evolution. Our final\ntest examines the role of meltwater production rates by increasing the modern-day initMIP11 basal\nwater production rates by a factor of 1.5, while keeping the basal velocities and ice and basal\ngeometries the same as in the source publication. This provides an approximation of how melt\ncould vary for ice sheet geometries governed by different ice sheet model physics.”\nMany of the figures presented quite small maps and therefore there was significant challenge in\nseeing the variation between model outcomes – this would need addressed somehow so that\npatterns can be more clearly understood and linked back to the text. I make specific comments on\nfigures below.\nWe have revised our figures to make the panels larger and text more legible. We hope the figures\nare now clearer to the reader.\nOther than the above points, I found a few minor issues in an otherwise interesting and well-written\npiece of work – very good clarity and context and it uses literature appropriately.\nMinor points:\n41: mass imbalance – state in which direction (I assumed negative!).\nWe have made this correction.\nFig 1: text on figure is too small to read easily on colour bars. Red 50 m contour also hard to see\nin detail (can panels be made a bit bigger?). Bed topograph also hard to make out in detail here\n(and we need that because its important for understanding the flow of basal water. Width of\ncolourbar is too narrow to easily make out colour - expand width slightly. There are red and green\non same plot (e) - is this colourblind friendly? Guess the contour is dashed so it may be OK. I also\nwant to know what the ice surface actually looks like under these scenarios - its really not easy to\ntell here.\nWe have updated Figure 1 with larger panels and text and wider colour bars. We have provided\na close up of the near-grounding line region for the bed topography and meltwater production\npanels; this region is outlined by the pink rectangle in panel a). In the supplement, we have now\nadded figures with the surface and basal topography for each time snapshot plotted side-by-side\n(see Figs S2-S6) and we have referenced these supplementary figures in the Figure 1 caption. We\nfeel our figure is much clearer now and we thank the reviewer for these suggestions.\n60-62: On MICI, perhaps relating what projections with or without MICI would do in this area -\ne.g. is there much difference - express as uncertainty in ice extent and thus it becomes important\nto understand processes under both scenarios.\nAlthough no specific study of the influence of MICI in the Aurora Subglacial Basin presently exists,\nwe have now included a figure in the supplement (Fig. S7) showing the geometry of the grounded\nice and the extent of future retreat for the CMIP6 forcings considered in our study with and without\nMICI. We have referred to this new figure in the manuscript on L65-66:\n“In the ASB, MICI can cause thousands of square kilometres of grounded ice to be lost (Fig. S7).”\n76: Impacts on the stability of ice flow. Such as? Mention some impacts briefly.\nWe have clarified this point in the manuscript and have added a reference to the appropriate paper\n(L80-82):\n“If the volume of discharge exiting these channels, or their locations, were to change, there could\nbe far-reaching impacts on the stability of ice flow in the ASB, such as accelerations of ice\nvelocities and grounding line retreat rates24.”\nFig 2: Delineate grounding line to show difference between ice shelf and grounded ice? Does\nlabelling need to be on all panels? Its not really clear where the ice shelf actually is - hence show\ngrounding line as mentioned above. Can these panels be bigger so we can see detail?\nOur GlaDS simulations are only performed over grounded ice; we have endeavoured to clarify\nthis in the main text and in the figure caption. We chose to label locations in all panels to\ndemonstrate how, at the same spatial scale, the drainage pathways have changed from one time\nsnapshot to another. We have revised the labelling so that the results are not blocked by the glacier\nnames. In the revised figure, we have used dotted lines to point to the glacier outlet as well as\nabbreviations of the glacier names. We have expanded on this in the figure caption. We have also\nincreased the size of the figure panel and added zoomed-in plots for individual simulations in the\nsupplement, which we have referenced in the figure 2 caption.\nMAIN TEXT (L86-89):\n“Here, we use the 2D, finite-element Glacier Drainage System24 (GlaDS, described in Methods)\nsubglacial hydrology model to assess how ice sheet geometry at different time snapshots in the\npast and future affect the subglacial hydrology of the ASB’s grounded ice and infer the effects on\nice shelf basal melt rates.”\nCAPTION:\nFigure 2. Modelled subglacial channel discharge from paleo, modern, and future time\nsnapshots. Panels show the channel discharge under grounded ice at a) the Eocene-Oligocene\nBoundary; b) the Mid-Miocene; c) present day; d) 2100 CE, under maximum volumes of ice loss;\nand e) 2100 CE under maximum volumes of ice gain. Black lines show the grounding lines,\nwhere water exits the domain, for each simulation. Dotted lines and abbreviations point to\nglacier locations. B: Bond Glacier; MU: Moscow University Ice Shelf; T: Totten Glacier; V:\nVanderford Glacier. Channel discharge for each glacier labelled in the figure panels at each\ntime snapshot is included in the supplementary material (Figs. S8-S12).\n119-123: I got a bit confused as to whether you were or were not modelling using MICI conditions.\nYour methods do indicate it, but in the text the water is muddied between what you say in line 119\nand what you way on line 123.\nWe thank the reviewer for pointing this out. Our study does not run an ice sheet model to generate\nthe ice sheet geometries. Rather, we used ice geometries from published ice sheet modelling that\nallowed MICI to be switched on or off: Paxman et al. (2020) for the paleo runs and Li et al. (2023)\nfor the future simulations. We have now clarified this point in the main text and in the methods:\nMAIN TEXT (L138-141):\n“We discuss how the ice geometries from Paxman et al.’s26 MICI and low CO experiments affect\n2\nour results in the section titled “Carbon emissions and instability processes influence drainage\nmorphology”.\nMETHODS (L523-528):\n“These data are from Paxman et al.26, who explored the sensitivity of Antarctic Ice Sheet evolution\nto climate changes. Their experiments forced an ice sheet model with CO concentrations of either\n2\n280 ppm or 500 ppm, allowing MICI to be turned on or off. Using the Paxman et al. ice sheet\nmodel outputs as inputs to the GlaDS model, we consider all possible combinations for each time\nsnapshot: 280 ppm without MICI, 280 ppm with MICI, 500 ppm without MICI, and 500 ppm with\nMICI.”\n144: capped at 800 ma-1. What were the maximum values otherwise? I am just trying to see the\nextent to which you curtail things here.\nThe maximum basal velocity for the Eocene-Oligocene Boundary 500 ppm No MICI run was\n6965.79 ma-1, 2685.67 ma-1 for the Mid-Miocene, 1717.18 ma-1 for the modern run, 892.18 ma-1\nfor the 2100 CE maximum ice loss scenario and 939.15 ma-1 for the 2100 CE maximum ice gain\nscenario. The exceptionally high ice velocities in the paleo runs are isolated points that occur at\nthe grounding line near Totten Glacier. We have expanded Table S1 to include a column of the\nuncapped maximum velocities.\n149-150: It would be good to have a bit more info on the interplay between topo melt and velocity\nhere - its a demonstration of what is causing the change and impacts of chance so seems odd to\nrelegate it all to supplementary.\nWe agree that expanding on these interactions in the main text would strengthen the manuscript,\nand we thank you for the suggestion. We have now included a summary of the sensitivity results in\nthe main text (L169-175):\n“We performed sensitivity tests to determine the relative influence of bed and surface topography,\nbasal melt, and basal velocity (Text S1). These tests indicated that the ice surface geometry is the\nprimary control on where dominant drainage pathways form, but that changes in bed topography\ncan influence channel lengths. Varying basal meltwater production rates and velocities caused\nmore than a 10-fold change in grounding line channel fluxes and changes in the distribution and\nextent of high-pressure regions, both of which can affect ice dynamics through their influence on\nice velocities and ice shelf basal melt rates.”\n195: ice configuration changes. by configuration, what do you mean? Extent? Thickness/surface?\nFlow? Which is more important as a control? Clarify please.\nWe have adjusted this for clarity, stating the importance of surface slope changes (L214-215).\n“The degree of similarity of future subglacial drainage to present or past drainage networks is\nhighly dependent on ice surface slope changes.”\n221: This sentence sounds a bit speculative. Can you make a more concrete statement on this based\non evidence?\nWe have updated this sentence to point directly to the melt results presented in this manuscript\nrather than the previous speculative statement (L240-243):\n“Consequently, discharge volumes and associated velocities are much larger than the other\nsimulations, which drive enhanced rates of ice shelf basal melting at deep grounding line points\nby locally increasing rates of turbulent sub-ice shelf ocean mixing near the discharge outflow\nlocations.”\n229-232: Were any sensitivity tests done whereby you use modern bed with palaeo ice\nsurface/extent? In such a way it would be possible to see what is more important for controlling\nthe water conditions - the bed vs. the ice geometry. This feeds into my second major comment\nabove – it was less clear to me that you had controlled for the palaeo bed (which in itself is a\nmodel) and its impact on water flow.\nYes, we performed sensitivity tests in which we used the BedMachine bed and the paleo ice surface.\nFor our model runs, including our sensitivity tests, ice thickness was always calculated as ice\nsurface elevation from the source publication minus bed topography (source publication bed\ntopography or BedMachine bed). In this way, the surface profile was always the same for a given\ntime snapshot, but the ice thickness varied depending on the bed geometry used. We have clarified\nthis in the supplementary material.\n273-280: How much difference was there between palaeo and future MICI runs in terms of the ice\ngeometry/thickness? Could be good to have some difference maps in supplement to highlight this\nand understand what it is about geometry that is controlling the hydraulic gradient.\nAs requested, in the supplement, we have now included ice thickness difference plots (No MICI\nminus MICI) for each time snapshot, as well as difference plots between paleo and future MICI\nsimulations (see Figs S24-S31). We have also added a table containing the maximum, minimum,\nand mean ice thickness differences for the pairs of differenced simulations (see Table S4).\nIn the supplement, we have added the following, which we reference in the main text as Text S2\n(L170-178):\n“We next investigated the influence of ice thickness changes in the MICI simulations on the\norganisation of the drainage system. The ice thickness changes between paleo and future MICI\nsimulations were greatest near the grounding line (Figs S24-S31). For example, at the Eocene-\nOligocene Boundary, ice thicknesses at the grounding line of Moscow University Ice Shelf were\nup to 1000 m greater than those suggested for the future. Despite these greater thicknesses,\nchannel discharges at Moscow University Ice Shelf were relatively low compared to what is\nanticipated for the future. This highlights the role of ice surface geometry, rather than absolute\nice thickness, in steering water flow, as greater ice thicknesses do not necessarily indicate\npreferential pathways for subglacial water flow.”\nFig 4 (title): is it enhancements, or is it simply melt rates (from which we can then interpret\nenhancements)? I.e. the colour scale is about melt, so the figure is about melt, not melt\nenhancements. Alter the title accordingly.\nWe have made this correction.\n314 (and this paragraph): for melt and how it changes in each simulation, would it be worth stating\nhow much the total melt within each system would be? Or the average beneath the ice shelf? That\nway we can understand more clearly not only the focussing (or not) but also whether more melt\nhappens overall or not.\nWe have added text in this paragraph to detail the magnitude of ice shelf-integrated discharge\nmelt enhancements over the Totten Ice Shelf as well as the Vanderford and Adams Ice Shelves\n(L342-346):\n“In the modern, future maximum ice gain, and future maximum ice loss runs, subglacial outflow\nleads to enhanced ice shelf basal mass loss of 0.65 Gt, 2.23 Gt, and 3.62 Gt over the Totten Ice\nShelf, respectively, and 0.21 Gt, 1.07 Gt, and 0.23 Gt over the Vanderford and Adams ice shelves,\nrespectively (total basal mass loss over the Totten and Vanderford/Adams ice shelves in lieu of\ndischarge are 178.44 Gt and 24.80 Gt, respectively).”\n334: can the manner in which these melt enhancements impact glacier retreat and mass loss be\nexpanded on a bit? This is a key implication so slightly more detail is needed. It links to my first\nsignificant comment above.\nWe have expanded on the discussion of how these discharge-driven melt enhancements can impact\nglacier dynamics, highlighting results from the new coupled ISSM-GlaDS paper (Pelle et al.,\n2024) that presents coupled projections of the Aurora Subglacial Basin through 2100 (see L359-\n367).\n“We note that in a complementary study, Pelle et al.24 developed a coupled ISSM-GlaDS\nconfiguration of ASB and, in projections to 2100, found that subglacial discharge-driven melt\nenhancements led to earlier onset of retreat of Totten Glacier’s southern grounding zone by ~25\nyears and continually accelerated retreat of Vanderford Glacier (up to 10 km of additional retreat\ncompared to uncoupled “ISSM-only” simulations). That is, while the magnitude of these discharge\nmelt enhancements is small relative to the background ocean-driven melt signal, their location at\ndeep grounding line points has cascading impacts on the upstream stress regime of entire glacial\nbasins”\n356: ice flow is likely sensitive. feels like speculation - e.g. not specifically shown in this paper.\nWe have revised this sentence to remove any speculation and have now referenced Pelle et al.\n(2024) (L388-390):\n“Our work complements the work by Pelle et al.24, who showed that ice flow is sensitive to the\npattern of melt enhancement, especially in regions such as the Aurora Subglacial Basin, which are\ngrounded below sea level and vulnerable to unstable and irreversible retreat.”\n386: water systems facilitate ice motion – I agree but its not directly shown here. For a more\nimpactful paper, this would need to be shown somehow, or at least expanded to really demonstrate\nwhat the impacts are in this case.\nWe have rephrased this sentence to better support our conclusion (L454-456):\n“Since subglacial meltwater affects basal lubrication and sliding speeds49, the higher meltwater\nproduction rates suggested in the past and future simulations imply knock-on effects on ice\ndynamics.”\n391-392: weaker buttressing also not shown here - for more impact, this would need to be\ndemonstrated.\nSimilarly, we have rephrased this section to better support this conclusion (L461-466):\n“Although our hydrology model was not coupled to ice and ocean models, our results shed light\non the effects of large-scale drainage reorganisation and the consequences of neglecting evolving\nsubglacial hydrology in ice sheet models. If high water pressure were to move upglacier and\nhigher grounding line channel discharges were to weaken ice shelf buttressing, subglacial water\nalone could drive sea level contributions from Antarctica that are larger than previously\nsuggested, without the need for instability mechanisms such as MICI.”\n440: basal meltwater production rates. these can be quite different depending on ice sheet model\nused. Does the potential variation introduced by ice sheet model differences need to be\nincorporated?\nWe thank the reviewer for raising this point.\nDifferences in the underlying physics in ice sheet models would likely influence the magnitude and\nspatial distribution of basal meltwater production rates. Coupled ice sheet and hydrology\nmodelling would be required to determine in what direction (increase/decrease) meltwater\nproduction rates would tend in long-term simulations (past and future) and where they would\noccur. Since our study aims to capture the overall characteristics of subglacial drainage and ice\nshelf basal melt in response to changes in basal conditions and ice geometries, we do not attempt\nan intercomparison of modelled subglacial hydrology from the suite of available ice sheet models\n(i.e., ISMIP6). However, in the updated submission, we have included an additional sensitivity test\nto approximate basal melt variability by taking the modern-day init-MIP basal melt rates and\nincreasing them everywhere by a factor of 1.5. We believe this gives an approximation of how melt\ncould vary for ice geometries with different ice sheet physics. The hydrology model outputs under\nthe 1.5x melt test and those under the standard configuration differ considerably – channel\ndischarges are up to 3.8x lower in the melt sensitivity test compared to the standard setup. This\nresult demonstrates the importance of accurately calculating and predicting basal melt rates. We\nhave added these new results to the supplementary material in Text S1 (L130-167) and Figs S22-\n23.\nSUPPLEMENTARY MATERIAL\n“To further explore the role of basal meltwater production rates, we performed a basal\nmelt rate sensitivity test in which we increased the modern-day initMIP basal meltwater\nproduction rates everywhere by a factor of 1.5, but used ice geometries, basal topographies, and\nbasal velocities from the source publication (Fig. S22p-t and Fig. S23p-t). The goal of this test is\nto examine how basal melt from models with ice sheet physics that differ from those that informed\nour subglacial hydrology modelling experiments (i.e., Li et al.4 and Paxman et al.3) could influence\nthe subglacial drainage system. Our results from this test indicated that total grounding line\ndischarge was up to 2.5 and 3.8 times lower than that from the standard simulation for the paleo\nand future runs, respectively. These much lower grounding line discharges are likely because,\neven after raising the modern melt rate by a factor of 1.5, the meltwater production rates in the\nsensitivity test were more than 5.5 times lower than that from the standard Eocene-Oligocene\nBoundary 500 ppm No MICI setup and 1.9 times lower in the Mid-Miocene 500 ppm No MICI run.\nCompared to the 1.5x melt test (Fig. S22t), the future maximum ice gain run with the standard\nsetup has higher rates of basal meltwater production across the entire domain, whereas the\nsensitivity test shows only isolated points of high basal melt production. There were smaller\ndifferences in channel discharges between the standard setup and the melt sensitivity test for the\nmaximum ice loss run (Fig. S22s). This may be due to computing limitations that restricted the\nsimulation time of the maximum ice loss melt sensitivity test to 35,000 days, rather than the 40,000\ndays used in all other simulations. We also found that the channel lengths were shorter in the 1.5x\nmelt run and that the location of the dominant channels was also affected. For the Eocene-\nOligocene boundary run, the maximum channel length for this sensitivity test was 122.3 km, ~250\nkm shorter than that of the standard setup. The maximum channel length for the Mid-Miocene run\nwas 107.2 km, ~163 km shorter than that of the standard run. Channel lengths were less affected\nfor the future maximum ice gain run; the longest channel was less than 60 km shorter than in the\nstandard setup. In the maximum ice loss scenario, the longest channel was ~291 km shorter than\nthat of the standard run. The spatial distribution of channels was also affected in this sensitivity\ntest. For both paleo runs, fewer channels crossed the grounding line and the channels that did\ndevelop were less anastomosing. This sensitivity test also resulted in substantial changes to the\nextent and distribution of high-pressure regions. For instance, for the Eocene-Oligocene\nBoundary, the fraction of the domain at or above 98% of ice overburden pressure was 2.67 times\nlower in the sensitivity test, compared to the standard setup. For the Mid-Miocene, only 9.3% of\nthe domain was at or above 98% of ice overburden pressure, compared to 43.7% in the standard\nrun. In the sensitivity tests for future runs, the fraction of the domain at high pressure decreased\nby a factor of 2.42 and 2.04 for the maximum ice gain and maximum ice loss runs, respectively.\nThe location of regions at overburden pressure also shifted closer to the grounding line in the\nsensitivity runs and did not extend as far inland as in the standard setup. These results underscore\nour finding from the modern inputs run, that basal meltwater production rates play a central role\nin the development and organisation of the subglacial drainage networks.\n472-472: Why choose to use ice thickness as opposed to ice surface geometry? If you use ice\nthickness, you would surely introduce some quite different ice geometries which may not be\nrealistic representations of those times. Perhaps a altering ice surface would have been more\nappropriate? In addition, were tests done where bed was varied but other things were not? That\nwould help get at whether the bed is relatively more or less important than the geometry. Perhaps\nexpanding the logic here would help with response to my second major comment above.\nWe thank the reviewer for this point. We did use ice surface geometry and ice thickness was\nalways calculated as ice surface elevation from the source publication minus bed topography\n(source publication bed topography or BedMachine bed). We have corrected this in the paper.\nYes, we included a test in which the bed was varied, but all other inputs were from the source\npublication. This was the BedMachine bed topography sensitivity experiment. We have also\nclarified this in the methods of the main text (L546-551):\n“We ran simulations where we set the basal topography to that of BedMachine Antarctica, but use\nthe ice surface geometry, basal velocity, and basal melt data from the source publication. Since,\nin this test, we are only changing the bed topography but keeping everything else the same as in\nthe source publication, this test allows us to determine the level of control that basal topography\nhas on the subglacial drainage system.”\n479-480: I am struggling to follow the logic here (sorry!) - I haven't seen you vary topography\nwhilst keeping ice surface the same - wouldn't that also be important to do?\nWe have now clarified this, see response to above comment.\nReviewer 2:\nI read this well-prepared and clearly written manuscript with great interest because there is\ncompelling geomorphologic evidence to suggest that subglacial meltwater, either from the\ndrainage of subglacial lakes or surface melting, may be influencing/may have influenced ice\ndynamics in the ASB throughout the Cenozoic. As such, modeling studies are required to assess\nthe potential impact of subglacial meltwater presence on ice behavior during past warm periods,\nparticularly in the Totten Glacier catchment, where thinning and ice mass loss and an active\nsubglacial hydrologic system are presently observed.\nTwo past time intervals were chosen as part of this study to reflect times when 1) regional ice\nvolume was larger than at present with higher pCO2 (E-O boundary) and 2) possibly smaller than\npresent with who knows what pCO2 (Middle Miocene). At least I think this was the goal because\nthe introduction talks about needing to understand how the subglacial hydrologic system worked\nduring past warm periods. However, the intervals listed are those of glacial expansion. Climate at\nthe EO boundary was relatively cold, and if the goal of the middle Miocene study was to hit the\nMiocene Climate Optimum, the warmest time interval of the last 30 Ma, that date and ice sheet\nconfiguration should reflect ~16 Ma and not ~14 Ma, when ice expanded and stabilized. The\nauthors conclude that there have been substantial reorganizations of the subglacial drainage system\nover ~40 million years. Well yes, this seems like a pretty safe takeaway for a manuscript.\nWe thank the reviewer for raising the point about the ice and climate conditions during the Mid-\nMiocene. In the manuscript, we have now clarified that, for the Miocene, we are looking at the\nMid-Miocene Climate Transition, as this is what was investigated in the Paxman et al. (2020)\npaper from which we acquired the modelled paleo ice geometries. Although, as the reviewer has\nmentioned, it would be ideal to model the hydrology during the Miocene Climate Optimum, the\nfull set of input data (i.e., basal velocity, basal meltwater production rates, ice and bed geometries)\nrequired for our subglacial hydrology model are as of yet unavailable for the Miocene Climate\nOptimum. These outputs would need to be generated through paleo ice sheet modelling\nexperiments, which we do not undertake here. We believe that examining the subglacial hydrology\nat the Mid-Miocene Climate Transition is a valuable output from our study since this time period\ncorresponds to highly variable climate conditions that had considerable impacts on the evolution\nof the Antarctic (Gasson et al., 2016; Levy et al., 2016; Shevenell et al., 2008).\nThere are often mismatches between geomorphological evidence and ice sheet modelling\nexperiments. Large-scale, long-term ice sheet modelling studies do not include evolving subglacial\nhydrology due to the high computational expense, with most ignoring hydrology entirely. Those\nfew studies that include hydrology only coarsely represent the system. It is conceivable that the\ndiscrepancies between geomorphological features and ice sheet model outputs are due to the lack\nof complex hydrological processes in ice sheet models. Our conclusion of the importance of\ntemporally changing hydrological systems is, therefore, key as it is an area frequently discounted\nin ice sheet modelling. In response to this specific comment and to emphasise our goals, in our\nrevised submission, we have indicated that our simulations are not intended to reproduce the exact\nconditions during each time snapshot, but rather to provide a large-scale overview of how the\nAurora Subglacial Basin’s subglacial drainage system is sensitive to changes in ice configurations\nand basal (melt, velocity, and topography) conditions, which is a key process for understanding\nglacial dynamic drivers in this region (L128-132):\nMAIN TEXT:\n“The purpose of these experiments is to provide a large-scale overview of how changes in ice\nconfigurations and basal (melt, velocity, and topography) conditions influence Aurora Subglacial\nBasin’s subglacial drainage system and are not intended to be interpreted as the exact subglacial\nconditions during each of these paleo time snapshots. In the Discussion, we compare our model\nresults to geomorphological evidence.”\nWe then go on to compare our subglacial hydrology outputs to the geologic evidence and highlight\nthat a closer match between these observations and ice sheet modelling is likely to be obtained\nwhen large-scale simulations of ice sheet evolution incorporate dynamic subglacial hydrology by\ncoupling the models together.\nAs an observational scientist, I understand that models do not often have the resolution that\nobservers want and that slight changes in boundary conditions dictated by observations may not\nmake large differences in model results. With this in mind, I have a few questions that probably\nneed clarification/addressing in the text:\nWhy were the Paxman et al boundary conditions chosen for this study, rather than the more\nregionally specific Aitken et al 2016 geomorphologic conditions?\nOur study aimed to explore the sensitivity of the subglacial drainage system to changes in ice\ngeometries that occur due to differences in climate forcing (i.e., different CO2 concentrations and\nclimate models) as well as highly debated mass loss mechanisms (MICI). To do this, we required\na full suite of inputs for each climate forcing to run the GlaDS hydrology model: basal topography,\nbasal meltwater production rates, basal velocities, and ice geometries. These data were available\nfrom Paxman et al. (2020) and Li et al. (2023) but not from the Aitken et al. (2016) paper.\nFurther, how do the model results presented in this study compare to those of Aitken, which\ncomment on basal velocities and channel systems visible in subglacial images, but do not prescribe\na specific hydrologic system, rather wet and cold-based regions based on geomorphology?\nThe reviewer raises an excellent point about comparing modelling results to observations, and\nthus, we have now added such a comparison in the discussion (L414-451) and in the\nsupplementary material (L181-212):\nMAIN TEXT (L414-451):\n“Geomorphological evidence indicates that the ASB experienced several advance and\nretreat cycles over the last 34 million years12,19–21,42,43. It has been suggested that surface\nmeltwater-rich regimes, such as those found in West Greenland44–48, may be responsible for the\ndynamic behaviour of the ice sheet margin19,21,42,43 (see also Text S3). However, constraining the\ntiming, volume, and location of meltwater inputs to the bed remains challenging due to the inherent\nuncertainties in interpreting the time-integrated paleo record, and these uncertainties translate to\ndifficulties in determining the most appropriate surface water routing scheme in hydrology models.\nBecause of these ambiguities, we performed a simple test in which the surface meltwater from\nPaxman et al.26 was added directly to the bed without seasonal variability. Channel discharge\nvolumes in this test increased, but the channel locations did not change since surface melt was\nconcentrated near the grounding line, where channels had already formed. If, in the coming\ncenturies, surface meltwater was to access the subglacial drainage network, a more dynamic\nsystem reminiscent of the paleo ASB and West Greenland may form, triggering highly variable\ngrounding line positions and ice flow patterns, even in the absence of MICI. We motivate future\nwork to build on what we have shown here by investigating variability in surface meltwater inputs\nto the bed, as this would provide insights into the stability of the ASB drainage network in the face\nof perturbations to melt supply.\nLarge-scale, long-term ice sheet modelling exercises are often run at coarse resolutions\nand with simplified ice dynamics to make the simulations computationally feasible. As a result,\nprocesses such as subglacial hydrology that are known to influence ice dynamics are generally\nneglected. The grounding line positions inferred by geomorphological records show larger\nvariations than those indicated by the ice sheet geometries that were produced by Paxman et al.26\nand, therefore, input to our GlaDS simulations. For instance, at the Eocene-Oligocene Boundary,\nice in the region may have expanded onto the continental shelf20,21,42, with grounded ice beyond\nour calculated domains. The incongruity between geomorphological data and Paxman et al.’s26\nice sheet geometries may be because the latter did not include a hydrology component and the\nsimulations were run at a relatively coarse resolution (10 km), two factors that could influence the\nrelative role of topographic pinning points on ice sheet retreat and expansion (e.g., Ridge A42; see\nalso Text S3). We do, however, find overlap between the inferred warm-based and fast-flowing\nregions indicated by Aitken et al.42 and areas of high velocities and melt production in the Paxman\net al.26 data that contribute to the high channel discharges in our subglacial hydrology\nsimulations. Isostatic adjustment may also influence long-term hydrology routing, but since we do\nnot couple our subglacial hydrology model to a solid earth model and perform long-term,\ncontinuous simulations, the direction or magnitude of potential feedbacks between solid earth\neffects and hydraulic gradients are unknown. A first step towards testing the suite of interactions\ncould be with fully-coupled, high-resolution studies of areas such as the ASB so that computational\nexpenses are minimised by focussing on one area while resolving the high degree of detail in basal\ntopography and subglacial drainage that is known to exist in the region12,19–21,24,22,42,43.”\nSUPPLEMENTARY MATERIAL (L181-212):\n“In glaciated and formerly glaciated regions, subglacial meltwater has left signatures in\nthe geomorphological record6. Tunnel valleys are features that are eroded into the bed from\nsubglacial meltwater and have been mapped across landscapes in Antarctica, Europe, and North\nAmerica, with each system of tunnel valleys displaying unique characteristics in terms of width,\ndepth, and geological setting6 . The ASB geomorphological record indicates the presence of tunnel\nvalleys on the continental shelf that likely formed during the Oligocene to late Miocene (~34 to ~8\nMya7). Gulick et al.7 suggest that the geometry of the ASB’s tunnel valleys indicates they formed\nbecause injections of surface meltwater to the bed exceeded the capacity of the subglacial drainage\nnetwork, which they argue consisted of a porous groundwater layer. The tunnel valleys mapped\nby Gulick et al.7 lie outside the limits of our GlaDS hydrology domains, which were informed by\nice sheet geometry data provided by Paxman et al.3. As a result, we cannot compare the tunnel\nvalleys with our basal channel routing. It is possible that the tunnel valleys are beyond the extent\nof our domains because Paxman et al.3 did not include surface or basal hydrology in their\nsimulations. Omitting hydrology may explain the mismatch between Paxman et al.’s3 modelled\ngrounding line positions and those inferred by the geomorphological record8 since Gulick et al.7\nindicated that surface meltwater inputs to the basal hydrology system likely played an important\nrole in determining the advance and retreat of ice in the ASB.\nWhile our subglacial model does not currently account for subglacial erosion processes\nresponsible for tunnel valley formation, there are bedrock features in our domain that could play\na role in determining the erosive potential of subglacial water and the larger-scale ice dynamics.\nFor instance, near the grounding line of the Moscow University Ice Shelf, there is a topographic\nridge, referred to as Ridge A/B in Aitken et al.8. If the subglacial hydrologic potential gradient is\ninsufficient to steer the water over the ridge, basal freeze-on can occur, which would prevent water\noutflow and erosion. Regarding ice dynamics, basal freeze-on at bedrock ridges suggests that ice\nflow would decelerate due to insufficient lubrication to allow for fast flow. If there was abundant\nsurface meltwater in the past, as suggested by the paleo record, this additional meltwater flux\ncould potentially alter the hydraulic gradient enough to drain water over the ridge, causing higher\nchannel discharge volumes than what we observe in our model results that could then carve the\ntunnel valleys. Changes in grounding line discharge because of altered hydraulic gradient at these\ntopographic ridges due to either freeze-on or surface meltwater inputs could, in turn, influence the\nmagnitude and distribution of ice shelf basal melt rates, thus feeding back onto grounding line\npositions and ice sheet stability.”\nPast grounding line configurations: The choice of grounding line configurations used in the paleo\nmodel runs is not clear. This may be a place where model resolution is important, but because\nthere is geological data to suggest that the chosen paleo grounding zone choices are incorrect, the\nchoice of paleo grounding line location warrants further discussion for each time interval. For\ninstance, is the resolution of the model in this manuscript such that an adjustment of the grounding line out\nonto the continental shelf could make a difference to the model output? It might be that the model cannot resolve changes in\nthe grounding line of 100s of Km, and that is fair and could be mentioned. I think this is what the authors say in the methods\nand supplement, but that should be made clear in the main text. Further, if this is the case, how realistic are these model\nresults, given the complex and relatively fine-scale subglacial topography in the region?\nWe agree that model resolution can be a limitation in reproducing features in the\ngeomorphological record. Our paleo grounding line locations are from Paxman et al. (2020),\nwhose results were posted at 10 km resolution. In the revised text, we have clarified that the paleo\ngrounding lines are from Paxman et al. (2020) (L124-125).\n“Our paleo subglacial hydrology simulations use grounding lines, ice sheet geometries, and basal\ninputs from Paxman et al.26”\nSince we are using pre-determined grounding lines, we are unable to vary their location or test\nthe influence of varying the grid resolution. Additionally, we have added a paragraph discussing\nthe implications of grid resolution, our results, and observations in the discussion (L414-451) in\nresponse to a comment from Reviewer 1 (see response above).\nWe note that the minimum edge length at the grounding line of our GlaDS hydrology model meshes\nis 0.25 km (Table S1).\nFor example, in Aitken’s and other studies, it is clear that the topographic high termed Ridge A\nmight be essential for stabilizing regional grounding lines—past and present. Given the position\nof Ridge A at the mouth of the ASB, which means that there is no reverse sloping bed until ice\nretreats landward of Ridge A, it might be useful for the authors to both address this feature and the\nevidence for ice advance onto the shelf during early ice sheet evolution. Additionally,\nconsideration of Ridge A is important because (based on the channels in the bedrock on the\nseaward side of Ridge A; Fernandez et al., 2018), subglacial meltwater may be really important\nfor getting ice over that high and out onto the shelf.\nIn our study, we are not modelling ice sheet evolution, and therefore cannot comment on the\ninfluence of Ridge A in determining where the grounding lines are located or how the ice advanced\nand retreated through time. However, importantly, we have now noted in the supplement how\nfeatures such as Ridge A may influence our hydrology results and the consequences for ice\ndynamics (see Text S3 and response to comment below).\nIn the manuscript, it appears that Ridge A is an important constraint on the location of the EO\ngrounding line. Is this a correct statement? If it is, it should be noted that data from the continental\nshelf presented by Gulick, Shevenell, et al. (2017; note that this is the correct citation for this\npaper) indicate grounded ice expansion onto the shelf, possibly at the EO boundary, based on\nexisting chronological constraints (not great, but data that exists until the shelf sediments are\ndrilled). This observation suggests that grounded ice had to make it over Ridge A and onto the\nshelf as ice in the ASB expanded. The authors should mention the shelf evidence when discussing\ntheir choice of grounding line configuration for the EO boundary.\nSince our grounding line configurations are from the Paxman et al. (2020) data, which we have\nnow made clear in the text in response to a comment from Reviewer 1, we cannot ascertain the\ninfluence of Ridge A on the location of the grounding line for the Eocene-Oligocene Boundary.\nHowever, in the discussion of our revised manuscript and the supplementary material, we have\ncommented on the expansion of ice during the Eocene-Oligocene and the potential role of Ridge\nA.\nMAIN TEXT (L434-442):\n“The grounding line positions inferred by geomorphological records show larger\nvariations than those indicated by the ice sheet geometries that were produced by Paxman et al.26\nand, therefore, input to our GlaDS simulations. For instance, at the Eocene-Oligocene Boundary,\nice in the region may have expanded onto the continental shelf20,21,42, with grounded ice beyond\nour calculated domains. The incongruity between geomorphological data and Paxman et al.’s26\nice sheet geometries may be because the latter did not include a hydrology component and the\nsimulations were run at a relatively coarse resolution (10 km), two factors that could influence the\nrelative role of topographic pinning points on ice sheet retreat and expansion (e.g., Ridge A42; see\nalso Text S3).”\nSUPPLEMENTARY MATERIAL (Text S3; L198-212):\n“While our subglacial model does not currently account for subglacial erosion processes\nresponsible for tunnel valley formation, there are bedrock features in our domain that could play\na role in determining the erosive potential of subglacial water and the larger-scale ice dynamics.\nFor instance, near the grounding line of the Moscow University Ice Shelf, there is a topographic\nridge, referred to as Ridge A/B in Aitken et al.8. If the subglacial hydrologic potential gradient is\ninsufficient to steer the water over the ridge, basal freeze-on can occur, which would prevent water\noutflow and erosion. Regarding ice dynamics, basal freeze-on at bedrock ridges suggests that ice\nflow would decelerate due to insufficient lubrication to allow for fast flow. If there was abundant\nsurface meltwater in the past, as suggested by the paleo record, this additional meltwater flux\ncould potentially alter the hydraulic gradient enough to drain water over the ridge, causing higher\nchannel discharge volumes than what we observe in our model results that could then carve the\ntunnel valleys. Changes in grounding line discharge because of altered hydraulic gradient at these\ntopographic ridges due to either freeze-on or surface meltwater inputs could, in turn, influence the\nmagnitude and distribution of ice shelf basal melt rates, thus feeding back onto grounding line\npositions and ice sheet stability.”\nThe authors should also address how Middle Miocene grounding zone locations were determined\nespecially in light of data from the continental shelf that indicates at least 10 advances and retreats\nonto the shelf between what is likely the E-O boundary and the Late Miocene (Gulick, Shevenell\net al., 2017). While the ages of these unconformities are unconstrained, they are significant in that\nthey also indicate large volumes of meltwater on the shelf based on the presence of tunnel valley\nfeatures. It seems that there is an unprecedented opportunity both address apparent data/model\ndisagreement and compare model results with what we know about subglacial hydrology in the\nregion from the geologic record.\nConcerning the choice of the Mid-Miocene grounding line location, we are unable to comment as\nthe grounding line configurations are from the Paxman et al. (2020) data. However, the reviewer\nraises an excellent point, that it is important to compare the model outputs to the geologic record.\nTherefore, we have now expanded on this in the discussion (L414-451, see response to above\ncomment) and supplementary material (Text S3; L181-212, see response to above comment).\nPast sea level: there is substantial data to suggest significant changes in sea level in the past. Was\npast sea level considered in the modeling exercise? How might past sea level impact the choice of\ngrounding line locations in the paleo model runs? We do know that sea level dropped to -30 to -\n40 m relative to modern day across the EO boundary (Miller et al., 2020), which likely enabled\ngrounded ice expansion from the ASB across the continental shelf, as observed in the geologic\nrecord. Further, during the Miocene Climate Optimum (~16 Ma), sea levels may have been +50m\nhigher than present (likely overestimate), while at ~14 Ma, sea levels may have been -40 m lower\nthan present. This is another reason why the authors really need to clarify if they are talking about\nthe MCO (~16 Ma) or ice expansion and stabilization associated with the Middle Miocene Climate\nTransition at ~14 Ma.\nOur subglacial hydrology model is run over grounded ice only, with the grounding line locations\nand ice geometries extracted from Paxman et al. (2020); we have now clarified this in the main\ntext and methods. As the grounding lines are determined by the ice thickness relative to sea level\nthe Paxman et al. (2020) data are our only control for sea level. Accounting for sea level changes\nwould require specialised ice sheet-sea level modelling experiments, which are beyond the scope\nof our study. We have also endeavoured to clarify in the main text that we are modelling the\nsubglacial hydrology at the Middle Miocene Climate Transition, as this is the time snapshot from\nPaxman et al.’s (2020) data (see L124-127 and response to above comment):\n“Our paleo subglacial hydrology simulations use grounding lines, ice sheet geometries, and basal\ninputs from Paxman et al.26. We focus on the Eocene-Oligocene Boundary (34 Mya) and the Mid-\nMiocene Climate Transition (14 Mya; hereafter Mid-Miocene) as these intervals correspond to\nice volumes in ASB that are larger and smaller, respectively, than those in the present day26.”\nPaleo channel locations: there are paleo channel systems identified both in the subglacial\nlandscapes (see Young et al., 2011, Aitken et al., 2016) inland of Ridge A and on the continental\nshelf seaward of the current study’s boundaries (see Gulick, Shevenell et al., 2017; Fernandez et\nal., 2018; Montelli et al., 2020). How do the model results compare with the locations of these\nchannels? Some discussion of the modelled channel locations vs observed channel locations would\nstrengthen the discussion of model results and the impact of the manuscript.\nThis is an excellent suggestion, and we agree that comparing our model results to observations\nwould strengthen the paper. Unfortunately, the data from the publications indicated above by the\nreviewer are either not publicly available, or in a format that we are unable to use. For instance,\nthe plotted data in Gulick et al.’s paper are from thousands of survey lines, the precise locations\nof the inferred tunnel valleys are not provided, and the depths are provided as two-way travel\ntimes. Converting the data into a usable form requires specialist knowledge of MB-System seafloor\nmapping software, which is outside the scope of our paper. Also, we approximated the tunnel\nvalley locations from Gulick et al.’s Figure 1b, line 21, and found that the inferred locations were\noutside of our model domain. However, in the supplementary material of our revised submission,\nwe have commented on the location of tunnel valleys and our modelled channel locations and the\npossible reasons for the mismatch (see Text S3).\nMICI as a viable mass loss mechanism: the authors admit the tenuous nature of the MICI mechanism in ice sheet models.\nThis reviewer is of the opinion that inclusion of MICI in the runs and in any discussion diminishes the results, as presented\n(mostly because the reader is left being skeptical of MICI, rather than thinking about the present modelling study). Would it\nbe better to think about and address how substantial subglacial meltwater presence might eliminate the need for a MICI-type\nmechanism to generate a responsive ice sheet? Such an argument, if possible, would really highlight the need to\nunderstand subglacial hydrologic systems, past and present.\nThe reviewer has raised concerns about our discussion of MICI, indicating that it detracts from\nthe impact of our results. We agree that, in our paper, we could have been clearer about the\nrationale for discussing MICI. In the revised paper’s introduction and final paragraph, while we\nhave kept the comparison of MICI and non-MICI results, we have now more explicitly mentioned\nthe implications of subglacial hydrology for initiating large-scale ice sheet changes without\nrequiring MICI. This is a key argument for incorporating the role of hydrology in ice sheet\ndynamic analysis and, therefore, we believe retaining the MICI experiments is a valuable output\nfor the community:\nL70-72 “Still, since most ice sheet models neglect evolving subglacial hydrology, it has seldom\nbeen considered that changes to subglacial drainage systems may be sufficient to initiate large-\nscale ice dynamic and geometric changes without the need for MICI.”\nL459-464 “Although our hydrology model was not coupled to ice and ocean models, our results\nallow us to shed light on the effects of large-scale drainage reorganisation and the consequences\nof neglecting evolving subglacial hydrology in ice sheet models. If high water pressure were to\nmove upglacier and higher grounding line channel discharges were to weaken ice shelf\nbuttressing, subglacial water alone could drive sea level contributions from Antarctica that are\nlarger than previously suggested, without the need for instability mechanisms such as MICI.”\nGenerally, after reading this manuscript, it seems that the authors have missed an opportunity to\ncompare model runs with existing observations of subglacial hydrologic systems in this catchment.\nThis is one of the only locations in Antarctica where such a comparison is possible. Such an\nexercise would significantly strengthen the believability and impact of the modeling studies. I\nrecommend such a data-model comparison be undertaken before publication is considered.\nWe agree that the manuscript is significantly strengthened with a comparison of model outputs to\nobservations and thus we have now added this to the revised manuscript (L414-451; see above\nresponse) and the supplementary material (Text S3; see above response).\nA few more minor questions/comments on the manuscript, in no particular order:\nThis study currently uses outdated CO2 reconstructions. I recommend using the latest CO2\ncompilation (CENCO2PIP Consortium, 2023).\nUnfortunately, since we are using reconstructions from Paxman et al. (2020), we are unable to\nuse alternative CO2 compilations. We hope this is now clearer in the text.\nFigure 1 and its caption are quite confusing (there might be a typo for the E-O simulation- should\nit be C and H panels?). It would be helpful to list the time intervals above the model domains and\nbasal input as in Figure 4, and maybe place vertical lines between the time intervals.\nReviewer 1 raised a similar point about Figure 1, and we have thus updated this figure.\nWhy does the Law Dome region disappear in the Mid-Miocene simulation? Is this related to sea\nlevel? And, if it is, the same would be true for the Eocene, would it not?\nWe thank the reviewer for raising this point. The domains for all simulations were calculated based\non hydraulic potential analysis using the ice and bed geometry from the source publication (in the\npaleo runs, these data were from Paxman et al. (2020). Because of this flow routing, Law Dome\nlies outside our calculated domain in the Mid-Miocene simulation, but within the domain for the\nEocene-Oligocene Boundary simulation. We have clarified this in the methods (L504-505):\n“Since hydraulic potential analysis is based on the geometry of the ice and bed, the source data\nin this calculation influence the flow routing and, therefore, the domain size.”\nI don't understand the discussion of ice shelves in past ice configurations. There is no evidence of\nice shelf presence in the Eocene or the middle Miocene, but rather grounded ice on the continental\nshelves. How does this fact impact the results of this study?\nWe apologise for causing confusion about the scope of our work and the influence of ice shelves.\nWe need to be clear that the goal of our study was not to model hydrology and ice shelf interactions\nexactly as they were and will be, nor to validate the exact configuration of ice during these time\nperiods, although Paxman et al. 2020 paleo topographies were constrained by independent\ngeological data. Rather, our aim was to show how the feedbacks between subglacial hydrology\nand ice shelves could influence the entire glacier system, and the importance of including these\ninteractions within a coupled ice sheet model, hydrology, and ice shelf framework. We were not\nsufficiently clear on this aim in the original manuscript, and we have endeavoured to make this\nclear in the revised version by amending the text in the introduction (L128-131):\n“The purpose of these experiments is to provide a large-scale overview of how changes in ice\nconfigurations and basal (melt, velocity, and topography) conditions influence ASB’s subglacial\ndrainage system and are not intended to be interpreted as the exact subglacial conditions during\neach of these paleo time snapshots.”\nFernandez et al (2018) present pretty interesting data for meltwater input to the shelf region during\nthe last deglaciation (see bedrock channels and stacked grounding zone wedges termed feather\nmoraines in the publication). Would this be a useful time interval over which to test the model,\nconsidering that we have knowledge about the location of the ice sheet and of subglacial meltwater\noutput/presence during the last deglaciation just based on continental shelf geomorphology.\nWe agree that this would be an interesting time snapshot to investigate. The Paxman et al. (2020)\nstudy did not perform simulations at the last deglaciation, and, therefore, we do not have ice sheet\nmodel data that will allow us to investigate the hydrology during this time period.\nWhat happens if surface-derived meltwater is included for the MCO run?\nWe were intrigued by this suggestion, and thus, we have investigated the distribution of surface\nmeltwater during the Mid-Miocene run. We ran a Mid-Miocene simulation where the surface\nmeltwater produced by Paxman et al.’s (2020) ice sheet modelling experiments was added directly\nto the bed. These additional meltwater inputs are spatially-variable, but temporally constant (i.e.,\nnot seasonally varying) and applied over the entire transient simulation of 10 000 days. Since this\nadditional melt is temporally fixed, there is no seasonal expansion and contraction of channels\nthat would accommodate time-varying volumes of meltwater, resulting in temporally varying\nchannel discharges. Therefore, this test likely represents the upper bound of the influence of\nsurface melt on basal hydrology. We adopt this simplified approach due to the uncertainties in the\ntiming, locations, and volumes of surface meltwater in the paleo record. We found that adding\nsurface meltwater increases the total discharge in channels but does not change the locations of\nchannels since the surface melt is concentrated near the grounding line, where channels have\nalready formed:\n[figure redacted]\nBased on this investigation, in the discussion, we have added a comment on the influence that\nadding surface melt has on the drainage system and noted it as a future research priority (L417-\n430):\n“However, constraining the timing, volume, and location of meltwater inputs to the bed remains\nchallenging due to the inherent uncertainties in interpreting the time-integrated paleo record, and\nthese uncertainties translate to difficulties in determining the most appropriate surface water\nrouting scheme in hydrology models. Because of these ambiguities, we performed a simple test in\nwhich the surface meltwater from Paxman et al.26 was added directly to the bed without seasonal\nvariability. Channel discharge volumes in this test increased, but the channel locations did not\nchange since surface melt was concentrated near the grounding line, where channels had already\nformed. If, in the coming centuries, surface meltwater was to access the subglacial drainage\nnetwork, a more dynamic system reminiscent of the paleo ASB and West Greenland may form,\ntriggering highly variable grounding line positions and ice flow patterns, even in the absence of\nMICI. We motivate future work to build on what we have shown here by investigating variability\nin surface meltwater inputs to the bed, as this would provide insights into the stability of the ASB\ndrainage network in the face of perturbations to melt supply.”\nReviewer 3:\nThis is a very well written paper on changing Aurora subglacial basin drainage at different climate\nstates that is based on a finite-element glacier drainage system model which responds to climate\nforcing through a link between ice thickness and water pressure at the bed. Using this methodology\nforced today with BedMachine for the topography and ISMIP6 for the climate they simulate\npotential drainage pattern under Totten which seems reasonable.\nThey then attempt similar situations at the Eocene-Oligocene boundary and at the mid-Miocene.\nFor the latter I think they are seeking to simulate the Mid-Miocene Climatic Optimum. There are\nissues with specifically matching the O-18 curve with CO2 as the authors are undoubtedly aware,\nbut I don’t think this matters for the paper since you are forcing it at two CO2 states roughly\nequivalent to E-O and MMCO. The specific date does not matter for your study since as noted\nbelow, I don’t think you are yet able to simulate what specifically occurred with the ice sheet and\ndrainage system at the E-O and MMCO. The larger issues however are three:\n1) bed topography has certainly evolved from the Eocene to today and thus what is the sensitivity\nof the model to this process since we do not have a E-O or Mid-Mio bedmachine? There are some\nsentences discussing how it might have changed but I was left unclear as the expected magnitude\nand feedbacks if say the deeper basins under ASB we not yet filled in and we had a subglacier\nlake scenario, or if the erosion through time has substantially changed the overall gradient to flatter\ntoday than in the E-O run. Can we be confident of these glacial drainage change results in absence\nof the topography? Some sensitivity study/discussion is needed to place these results in context on\nthis issue.\nWe thank the reviewer for raising these points about the bed topography. We realise we could have\nbeen clearer about the sensitivity tests that we performed; Reviewer 1 also mentioned this (please\nsee our response to their comment above), and in the revised version of the manuscript, we have\nnow summarised the key findings from Text S1, which details the results of our sensitivity tests.\nBelow, we provide a review of these sensitivity tests.\nOur first test investigated the sensitivity of drainage pathways to bed topography by keeping the\nice surface the same as in the source publication, but varying the bed topography, using either the\ntime-snapshot-specific bed topography or the modern BedMachine topography. For the paleo\nsimulations with the BedMachine basal topography, we found that channel lengths were affected,\nbut the main channel outlet locations remained the same.\nIn our second test, we kept the ice surface geometry the same as in the source publication but\nvaried the bed topography (using BedMachine), and modern-day basal velocities and meltwater\nproduction rates (the “modern inputs” sensitivity test). The outputs from this test allowed us to\nisolate the role of ice surface geometry in steering water flow. For the Eocene-Oligocene\nBoundary, this test indicated that channelised drainage beneath Totten remains dominant,\nhighlighting the steep surface slopes during this time period. In contrast, drainage through Totten\nis lower during the Mid-Miocene, which we attribute to weaker surface slopes. Despite the similar\ndrainage pathways to the runs with the time-snapshot specific data, the discharge volumes in the\nmodern inputs sensitivity test were lower than in the runs using time-snapshot-specific data.\nOur final sensitivity test, which is a new addition to our manuscript based on comments from\nReviewer 1, investigated the role of variability in basal meltwater production rates. In this test, we\nincreased the modern basal meltwater production rate everywhere by a factor of 1.5, and applied\nit to our paleo and future simulations while keeping the bed and ice geometries, as well as the\nbasal velocities, the same as those in the source publication. This test approximates the variability\nin meltwater production that may result from ice sheet models whose physics differ from those in\nthe ice sheet models of Paxman et al. (2020) and Li et al. (2023). This test resulted in shorter and\nless extensive channels with lower discharges than that of the standard setup, as well as a\nreduction in the extent of high pressure regions. Overall, this test highlights the role that basal\nmeltwater production can have in subglacial drainage development and that tighter constraints\non this basal boundary condition are required.\nRegarding the magnitude and feedbacks of bed topography changes, our subglacial hydrology\nmodel is not coupled to erosion and solid earth models, as this is outside the scope of our work.\nWe are, therefore, unable to provide an estimate of the direction and magnitude of feedbacks\nbetween solid earth uplift, erosional processes, and subglacial hydrology. However, we have\nmentioned this limitation in the revised manuscript (L431-451) and in Text S3 of the supplementary\nmaterial (see responses to the above comments from Reviewer 2). We note, however, that the input\ndata from Paxman et al. (2020) did account for isostatic adjustment, erosion, sedimentation,\nvolcanism, and plate motions in their simulations.\n2) The ice margin for the ASB portion of East Antarctica is not substantially changing between\nthe E-O, Mid-Mio, Today, and Future scenarios. There is certainly some retreat and some advance\naspects but these are way too small to match the offshore shelf and slope geophysical data implying\npresence or absence of grounded ice (Gulick et al., 2017, Nature, Montelli et al. 2019 GSA\nBulletin, McKay et al., 2022, Ant Clim Evol). For instance, at the E-O boundary there were tunnel\nvalleys > 125 km seaward of the current ice margin which require grounded ice and we have a\nglacial shelf crossing trough (Donda et al., 2020, Ant Sci; Hirano et al., 2023, Nature Comm)\nreaching the shelf edge during the advances. There is strong data for grounded ice therefore over\nmuch of the shelf and within the trough to the shelf edge since there is a large trough mouth fan\nfor the Totten system active in the Oligocene at least. Therefore your E-O simulation ice margin\nshould be substantially more advanced to match this geologic evidence. In counterpoint, there\nappears to have been no evidence of proximal glacial sedimentary products in the MMCO on the\nshelf and thus one would infer a much more retreated ice margin to explain this. Your Mid-Mio\nsimulation is somewhat retreated but probably far too little to explain the absence of glacial\nsedimentary facies at the MMCO on the shelf.\nWe thank the reviewer for this comment. Reviewer 2 raised similar concerns, which we have\naddressed above. Since we are not performing ice sheet modelling experiments, and our subglacial\nhydrology model is supplied with data from Paxman et al. (2020), we are unable to change the\nlocation of the grounding line. However, in the revised manuscript we have endeavoured to clarify\nthat the subglacial drainage systems we model are not intended to be “true” representations of\nwhat the hydrology system looked like in the past but rather to provide a large-scale overview of\nhow changes in ice configurations and basal conditions influence Aurora Subglacial Basin’s\nsubglacial drainage system (L128-131):\n“The purpose of these experiments is to provide a large-scale overview of how changes in ice\nconfigurations and basal (melt, velocity, and topography) conditions influence Aurora Subglacial\nBasin’s subglacial drainage system and are not intended to be interpreted as the exact subglacial\nconditions during each of these paleo time snapshots.”\nWe have also commented on the role that coupled subglacial hydrology and ice dynamics\nmodelling may have in reproducing paleo grounding line positions (see L434-451 and the\nresponses to Reviewer 2’s comments above). We hope our findings will inspire such coupled\nexperiments to better represent ice extent during these time periods.\n3) The lack of greater changes in ice margin at warmer and colder periods imply the current models\nare missing some key feedbacks. It could be the effect of ice cliff instability as raised by the\nauthors. However another key change that should occur is going from colder to warmer based\nglaciers such as is seen in north versus south Greenland outlet glaciers where warmer wetter ones\nmoving larger volumes of sediment potentially building glacial moraines or morainal banks, and are\ncorrelated with tunnel valleys existing when advancing onto shelves (e.g., Alaska, North Sea).\nColder states act more like modern Antarctica or the colder based parts of Greenland forming\ngrounding zone wedges instead and substantially less feedbacks are likely. As Totten system\ntransitioned to this more polar state only the Late Miocene-Pliocene, it is probably incorrect to\nassume the subglacial system worked the same way from a more modern polar state as it did in a\npast more dynamic state. In other words E-O and Mid-Mio likely dynamically behaved more like\nthe North Sea or Alaska systems than Pleistocene Antarctica in order to explain the greatly\ndifferent shelf margin effects and glacial sedimentary products. What exactly is the change in\nfeedbacks is that allows for much more dynamic advance and retreat cycles is not yet known of\ncourse but some combination of melt water, warm bases, and glacial sediment dynamics are likely\ninvolved.\nMoving forward with this paper I think the authors need to specifically note how the paleo model\nruns do not produce dynamic enough advances to reach the shelf edge during the E-O or leave the\nshelf entirely ice free during the Mid-Miocene and then discuss the possible additional feedbacks\nneeded to result in greater ice margin changes. Then discuss what effect this would have on the\nsubglacial networks that are the focus of this paper. In summary, this paper does a good job with\nthe paleo models if the bed topography did not differ substantially and if the glacial dynamics did\nnot change going into past climate states. Issue is we know the ice was a lot further advanced and\nhave evidence of a wetter more dynamic system that created tunnel valleys, trough mouth fans,\nand probably intermittent morainal banks … how is this kind of system, one more representative\nof the LGM Alaska, North Sea, or West Greenland for instance, differ from what is modeled in\nthis paper? Given there is clearly a mismatch in ice edge position with the geologic record at least\nfor maximum advance and retreat climate states, then how can we use this paper’s results to inform\nus of process without yet being able to include all the processes required? I think the path forward\nis to use the paleo model results to discuss how different inputs such as CO2 affect ice thickness\nand water pressure and drainage without making the argument that this study is faithfully\nreproducing the conditions at E-O and Mid-Mio. Consider what you models are not including and\nhow you might then frame your model results differently since its interesting the tradeoffs\nobserved in water pressure and drainage length in your work, and considering future scenarios (see\nbelow), but these models clearly do not represent the state of the ASB glacial system at the E-O or\nmid-Mio.\nWe thank the reviewer for raising these points. We’d like to emphasise that our work is exploratory\nin nature, being the first to compare past, present, and future hydrology and ice shelf melt using\nevolving 2-D models. Therefore, we are not attempting to recreate the exact conditions at each\ntime snapshot. We do, however, agree that acknowledging the limitations of our modelling\napproach is important for framing the results and their implications, and reviewer 2 raised similar\nconcerns about our modelled system and that suggested by the paleo record. In the discussion of\nour revised manuscript and in Text S3, we have noted the discrepancies between our modelling\nresults and observations. We have now provided suggestions for future work such that a closer\nmatch between model outputs and geomorphological evidence can be obtained. We have\naddressed the role of solid earth uplift, surface meltwater and erosion processes in the region and\nmotivated the use of coupled ice sheet-hydrology models as this may lead to grounding line\nlocations more representative of those documented in the paleo record. We hope these\namendments, along with our clarification of the aim and intention of this manuscript, have satisfied\nthese concerns.\nFor the future model runs, I think these are fine and in fact really interesting suggestions for future\nstates, if the system stays in its current polar mode of dynamics. However, if it switches modes to\nthe way it operated in the Oligocene to Miocene then there would be expected a much greater\ndynamic change possibility that could be acknowledged.\nWe thank the reviewer for the positive feedback, and we have now made this comparison in the\nrevised version (L424-427):\n“If, in the coming centuries, surface meltwater were to access the subglacial drainage network,\na more dynamic system reminiscent of the paleo ASB and West Greenland may form, triggering\nhighly variable grounding line positions and ice flow patterns, even in the absence of MICI.”\nOn the writing, this paper is remarkably free of typos, has excellently written sentences, and was\nvery clear to read. Only minor item on that front is that the word further (a concept) is used in two\nplaces where the word should be farther (a distance).\nWe have made this correction.\nMy recommendation is therefore Major Revisions to consider the issues with mismatch to the\ngeologic record and suggestions of glacial dynamics having a warmer wetter regime with processes\nnot yet included in the finite glacial drainage system modeling.\nHopefully this review will be helpful and I would be willing to do a re-review as needed,\nSean Gulick\nReply to reviewers\nThe past, present, and future evolution of Aurora Subglacial Basin’s subglacial drainage system\nHayden, Pelle, and Dow\nWe would like to thank the reviewers for their comments on our revised manuscript. We have\naddressed the comments point-for-point below. The comments are laid out in normal font, and our\nreplies are in italic green font.\nReviewer #1 (Remarks to the Author):\nThis is a re-review of Hayden et als paper on the subglacial drainage system of the Aurora\nsubglacial basin. Thank you to the authors for taking the time and care to complete a series of\nadjustments to the paper. I think the responses are fairly comprehensive and I am struggling to find\nissue with any of them - I think you have been balanced and fair in your responses and its good to\nsee some additional detail of other model runs come in to back up some key points. There is good\ndetail in the additional text, both in the manuscript and in the more comprehensive supplementary\ninfo.\nThus overall I think this paper is a very useful study that demonstrates the potential drainage\npatterns in this region and how and why they may evolve. I won't re-iterate comments relating to\nthe strength and importance of the work that I made in my previous review (reviewer 1) - those all\nstill stand and I think this will make an excellent contribution to this journal. I really don't have\nany further comments to make in terms of requesting further changes - it looks great to me.\nStewart Jamieson\nWe thank the reviewer for re-appraising and supporting our study.\nReviewer #3 (Remarks to the Author):\nThis is a revised submission of a very interesting and well written modeling study by Hayden et al\non the effects of subglacial hydrology on ice sheet dynamics in the Aurora subglacial basin, East\nAntarctica. In general, I accept a lot of the response to reviewers aspects that the authors raised\nand appreciate the edited discussion in the main text and supplemental.\nThank you.\nHowever, I recommend an additional minor set of edits which are a suggestion to be more specific\nin the caveats starting with line 430. Right now the paper lumps the geological evidence under the\nterminology \"geomorphologic evidence\" (or some variant) which is partly correct but misses also\nthe glacial sequence stratigraphy on the shelf (i.e. the geologic record) published in several of the\npapers you cite in this section.\nWe agree that expanding the discussion to include a more detailed comparison between the\nstratigraphic record and our results would strengthen the manuscript. In the revised manuscript\nand summarised below, we address the reviewer’s suggestions.\nSo I strongly recommend not just mentioning that the modeling of ice extents is known to be\ninconsistent with shelf morphology (e.g. presence of a shelf crossing trough and lateral moraines)\nduring ice advances, but also note that it is inconsistent with presence of stratigraphic indicators\nof grounded ice (e.g. tunnel valleys and sedimentary deposits) both for the advanced and retreated\npositions of the ASB marine terminating ice.\nWe thank the reviewer for pointing this out. In the revised manuscript, we have explicitly stated\nthe geologic markers in the sedimentary record and how our domains and model results compare\nto this evidence.\nL434-440: “In the sediment stratigraphy of the ASB’s continental margin, there are repeat\nsignatures of sediment wedges and bedrock channels shaped by subglacial meltwater flows20,21,43,\ndating from the Eocene to the Miocene. These features are located beyond the grounding lines of\nthe Paxman et al.26 data (and, hence, our GlaDS domains) and indicate that fluctuations in ASB\nterminus positions were much larger during both cold and warm climate states than what Paxman\net al.’s26 results suggest. These mismatches in grounding line locations, in turn, affect our ability\nto compare our modelled hydrologic pathways to those preserved in the geologic record.”\nI understand ice extents are baked into the boundary conditions of your model so you cannot\nchange this aspect, but you can make it clear that your modeling therefore greatly under-predicts\nthe terminus response due warm states (e.g. the MMCO) and cold states (e.g. E/O boundary)\nWe have endeavoured to clarify this in the response to the preceding comment.\nand a more complete list of the key factors/feedbacks could matter and need to be integrated into\nfuture studies. For instance I recommend adding in that one key feedback not included within this\nmodeling or in Paxman et al is the role of subglacial sediments in addition to subglacial hydrology.\nFor instance your long hydrologic conduits very likely should become tunnel valleys in the model\nwith a feedback built in that these are erosive sub glacially but depositional near the terminus and\nthus can provide a direct feedback to the ground line stability, which likely then changes under\ndifferent flow regimes/with connection to surface melt/with different hydrologic pressures (as you\nnote).\nWe agree that expanding our discussion of the feedback mechanisms at play but not currently\nincluded in our or Paxman et al.’s modelling approaches is important for highlighting future\nresearch priorities. The revised manuscript describes the interactions between erosional and\ndepositional processes, subglacial hydrology, and ice dynamics.\nL444-458: “The absence of subglacial hydrology in Paxman et al.’s26 model could have also\naffected their computed rates and patterns of subglacial erosion and sediment deposition, with\ncascading effects on the locations and advance/retreat rates of terminus positions. For example,\nsubglacial meltwater that is continuously routed along the same path could create the above-\nmentioned bedrock channels under grounded ice. The sediment mobilised by subglacial meltwater\ncould then be deposited further offshore, where it could lead to the formation of topographic\npinning points that stabilise grounding line positions and affect ice dynamics. Whereas Paxman et\nal.26 did not incorporate hydrology, our hydrology model does not resolve erosional and\ndepositional processes, and therefore, we are unable to capture the redistribution of sediment in\nthe system by subglacial meltwater and the resulting geomorphological features. The limitations\nof both Paxman et al.’s26 and our model highlight an opportunity for future work to incorporate\nerosion, sedimentation, and hydrology in ice sheet models since the bed is not static but instead\nevolves in concert with, and feeds back on, the ice sheet and hydrology systems. By accounting for\nthese processes, modelling stands to match paleo records better and produce more tightly\nconstrained projections of ice sheet evolution.”\nI don't think you can ignore the role of erosion/redistribution of sediment within the subglacial\nhydrologic system and so should at least include this within the discussion as to potential reasons\nwhy the ice terminus is not as responsive as the geologic record has evidence for and help bridge\nthe gap between the current modeling of a stable terminus (despite lots of mass balance changes)\nand what the geologic record indicates.\nWe have addressed this in our response to the previous reviewer comment.\nFor instance in 470-474, could add \"with changes to the bed and sediment redistribution being a\npotentially important feedback\".\nWe have edited this sentence to include the suggested modification.\nL487-490: “It is possible that a limited number of channels with high discharges and focused over\na small length of the grounding line could have the same impacts on the (de) stabilisation of ice\ndynamics as a scenario where many smaller channels occupy the entire grounding line, with\nchanges to the bed and sediment redistribution being a potentially important feedback.”\nFor 474 the list \"ice sheet, ocean, subglacial hydrology\" is great but also should include \" and \"the\nbed/subglacial geology/sediments\" pick your term but the feedbacks don't just include water and\nice.\nWe have made this change. The sentence now reads:\nL490-493: “Determining the nature of these relationships requires coupling between the ice sheet,\nthe ocean, subglacial sediments, and the subglacial drainage system and will provide crucial\ninformation for assessing past, contemporary, and future sea level contributions.”\nWhile this may seem like a minor suggest for editing, I personally think its critical that the reader\nof this important paper understand the real differences in ice sheet size/terminus position present\nin the geologic record but not yet accounted for in this work and that it be communicated feedbacks\nthat ultimately affect sea level do not end at the bed since the bed itself is also dynamic and part\nof the system\nWe would like to thank the reviewer for their in-depth comments and suggestions. We have revised\nour manuscript to incorporate the reviewer’s recommendations and believe that these edits have\nstrengthened the manuscript.\nReply to reviewers\nThe past, present, and future evolution of Aurora Subglacial Basin’s subglacial drainage system\nHayden, Pelle, and Dow\nReviewer comments are laid out in normal font, and our replies are in italic green font.\nReviewer #3 (Remarks to the Author):\nI wish to congratulate the authors on a very thorough job of revisions throughout this submission\nand review process. This latest version of their manuscript has satisfied all my earlier concerns\n(Reviewer #3) and has done so eloquently. I suggest that this mature version is ready for\npublication.\nSean Gulick\nWe would like to thank the reviewer for supporting our study and for their constructive feedback\nthat strengthened the manuscript."}, {"round": 2, "reviewer_comments": {"reviewer_1": "This is a re-review of Hayden et als paper on the subglacial drainage system of the Aurora subglacial basin. Thank you to the authors for taking the time and care to complete a series of adjustments to the paper. I think the responses are fairly comprehensive and I am struggling to find issue with any of them - I think you have been balanced and fair in your responses and its good to see some additional detail of other model runs come in to back up some key points. There is good detail in the additional text, both in the manuscript and in the more comprehensive supplementary info.\nThus overall I think this paper is a very useful study that demonstrates the potential drainage patterns in this region and how and why they may evolve. I won't re-iterate comments relating to the strength and importance of the work that I made in my previous review (reviewer 1) - those all still stand and I think this will make an excellent contribution to this journal. I really don't have any further comments to make in terms of requesting further changes - it looks great to me.\nStewart Jamieson", "reviewer_3": "This is a revised submission of a very interesting and well written modeling study by Hayden et al on the effects of subglacial hydrology on ice sheet dynamics in the Aurora subglacial basin, East Antarctica. In general, I accept a lot of the response to reviewers aspects that the authors raised and appreciate the edited discussion in the main text and supplemental.\nHowever, I recommend an additional minor set of edits which are a suggestion to be more specific in the caveats starting with line 430. Right now the paper lumps the geological evidence under the terminology \"geomorphologic evidence\" (or some variant) which is partly correct but misses also the glacial sequence stratigraphy on the shelf (i.e. the geologic record) published in several of the papers you cite in this section. So I strongly recommend not just mentioning that the modeling of ice extents is known to be inconsistent with shelf morphology (e.g. presence of a shelf crossing trough and lateral moraines) during ice advances, but also note that it is inconsistent with presence of stratigraphic indicators of grounded ice (e.g. tunnel valleys and sedimentary deposits) both for the advanced and retreated positions of the ASB marine terminating ice. I understand ice extents are baked into the boundary conditions of your model so you cannot change this aspect, but you can make it clear that your modeling therefore greatly under-predicts the terminus response due warm states (e.g. the MMCO) and cold states (e.g. E/O boundary) and a more complete list of the key factors/feedbacks could matter and need to be integrated into future studies. For instance I recommend adding in that one key feedback not included within this modeling or in Paxman et al is the role of subglacial sediments in addition to subglacial hydrology. For instance your long hydrologic conduits very likely should become tunnel valleys in the model with a feedback built in that these are erosive sub glacially but depositional near the terminus and thus can provide a direct feedback to the ground line stability, which likely then changes under different flow regimes/with connection to surface melt/with different hydrologic pressures (as you note). I don't think you can ignore the role of erosion/redistribution of sediment within the subglacial hydrologic system and so should at least include this within the discussion as to potential reasons why the ice terminus is not as responsive as the geologic record has evidence for and help bridge the gap between the current modeling of a stable terminus (despite lots of mass balance changes) and what the geologic record indicates. For instance in 470-474, could add \"with changes to the bed and sediment redistribution being a potentially important feedback\". For 474 the list \"ice sheet, ocean, subglacial hydrology\" is great but also should include \" and \"the bed/subglacial geology/sediments\" pick your term but the feedbacks don't just include water and ice.\nWhile this may seem like a minor suggest for editing, I personally think its critical that the reader of this important paper understand the real differences in ice sheet size/terminus position present in the geologic record but not yet accounted for in this work and that it be communicated feedbacks that ultimately affect sea level do not end at the bed since the bed itself is also dynamic and part of the system."}, "author_response": "Reply to reviewers\nThe past, present, and future evolution of Aurora Subglacial Basin’s subglacial drainage system\nHayden, Pelle, and Dow\nWe would like to thank the reviewers for their comments on our revised manuscript. We have\naddressed the comments point-for-point below. The comments are laid out in normal font, and our\nreplies are in italic green font.\nReviewer #1 (Remarks to the Author):\nThis is a re-review of Hayden et als paper on the subglacial drainage system of the Aurora\nsubglacial basin. Thank you to the authors for taking the time and care to complete a series of\nadjustments to the paper. I think the responses are fairly comprehensive and I am struggling to find\nissue with any of them - I think you have been balanced and fair in your responses and its good to\nsee some additional detail of other model runs come in to back up some key points. There is good\ndetail in the additional text, both in the manuscript and in the more comprehensive supplementary\ninfo.\nThus overall I think this paper is a very useful study that demonstrates the potential drainage\npatterns in this region and how and why they may evolve. I won't re-iterate comments relating to\nthe strength and importance of the work that I made in my previous review (reviewer 1) - those all\nstill stand and I think this will make an excellent contribution to this journal. I really don't have\nany further comments to make in terms of requesting further changes - it looks great to me.\nStewart Jamieson\nWe thank the reviewer for re-appraising and supporting our study.\nReviewer #3 (Remarks to the Author):\nThis is a revised submission of a very interesting and well written modeling study by Hayden et al\non the effects of subglacial hydrology on ice sheet dynamics in the Aurora subglacial basin, East\nAntarctica. In general, I accept a lot of the response to reviewers aspects that the authors raised\nand appreciate the edited discussion in the main text and supplemental.\nThank you.\nHowever, I recommend an additional minor set of edits which are a suggestion to be more specific\nin the caveats starting with line 430. Right now the paper lumps the geological evidence under the\nterminology \"geomorphologic evidence\" (or some variant) which is partly correct but misses also\nthe glacial sequence stratigraphy on the shelf (i.e. the geologic record) published in several of the\npapers you cite in this section.\nWe agree that expanding the discussion to include a more detailed comparison between the\nstratigraphic record and our results would strengthen the manuscript. In the revised manuscript\nand summarised below, we address the reviewer’s suggestions.\nSo I strongly recommend not just mentioning that the modeling of ice extents is known to be\ninconsistent with shelf morphology (e.g. presence of a shelf crossing trough and lateral moraines)\nduring ice advances, but also note that it is inconsistent with presence of stratigraphic indicators\nof grounded ice (e.g. tunnel valleys and sedimentary deposits) both for the advanced and retreated\npositions of the ASB marine terminating ice.\nWe thank the reviewer for pointing this out. In the revised manuscript, we have explicitly stated\nthe geologic markers in the sedimentary record and how our domains and model results compare\nto this evidence.\nL434-440: “In the sediment stratigraphy of the ASB’s continental margin, there are repeat\nsignatures of sediment wedges and bedrock channels shaped by subglacial meltwater flows20,21,43,\ndating from the Eocene to the Miocene. These features are located beyond the grounding lines of\nthe Paxman et al.26 data (and, hence, our GlaDS domains) and indicate that fluctuations in ASB\nterminus positions were much larger during both cold and warm climate states than what Paxman\net al.’s26 results suggest. These mismatches in grounding line locations, in turn, affect our ability\nto compare our modelled hydrologic pathways to those preserved in the geologic record.”\nI understand ice extents are baked into the boundary conditions of your model so you cannot\nchange this aspect, but you can make it clear that your modeling therefore greatly under-predicts\nthe terminus response due warm states (e.g. the MMCO) and cold states (e.g. E/O boundary)\nWe have endeavoured to clarify this in the response to the preceding comment.\nand a more complete list of the key factors/feedbacks could matter and need to be integrated into\nfuture studies. For instance I recommend adding in that one key feedback not included within this\nmodeling or in Paxman et al is the role of subglacial sediments in addition to subglacial hydrology.\nFor instance your long hydrologic conduits very likely should become tunnel valleys in the model\nwith a feedback built in that these are erosive sub glacially but depositional near the terminus and\nthus can provide a direct feedback to the ground line stability, which likely then changes under\ndifferent flow regimes/with connection to surface melt/with different hydrologic pressures (as you\nnote).\nWe agree that expanding our discussion of the feedback mechanisms at play but not currently\nincluded in our or Paxman et al.’s modelling approaches is important for highlighting future\nresearch priorities. The revised manuscript describes the interactions between erosional and\ndepositional processes, subglacial hydrology, and ice dynamics.\nL444-458: “The absence of subglacial hydrology in Paxman et al.’s26 model could have also\naffected their computed rates and patterns of subglacial erosion and sediment deposition, with\ncascading effects on the locations and advance/retreat rates of terminus positions. For example,\nsubglacial meltwater that is continuously routed along the same path could create the above-\nmentioned bedrock channels under grounded ice. The sediment mobilised by subglacial meltwater\ncould then be deposited further offshore, where it could lead to the formation of topographic\npinning points that stabilise grounding line positions and affect ice dynamics. Whereas Paxman et\nal.26 did not incorporate hydrology, our hydrology model does not resolve erosional and\ndepositional processes, and therefore, we are unable to capture the redistribution of sediment in\nthe system by subglacial meltwater and the resulting geomorphological features. The limitations\nof both Paxman et al.’s26 and our model highlight an opportunity for future work to incorporate\nerosion, sedimentation, and hydrology in ice sheet models since the bed is not static but instead\nevolves in concert with, and feeds back on, the ice sheet and hydrology systems. By accounting for\nthese processes, modelling stands to match paleo records better and produce more tightly\nconstrained projections of ice sheet evolution.”\nI don't think you can ignore the role of erosion/redistribution of sediment within the subglacial\nhydrologic system and so should at least include this within the discussion as to potential reasons\nwhy the ice terminus is not as responsive as the geologic record has evidence for and help bridge\nthe gap between the current modeling of a stable terminus (despite lots of mass balance changes)\nand what the geologic record indicates.\nWe have addressed this in our response to the previous reviewer comment.\nFor instance in 470-474, could add \"with changes to the bed and sediment redistribution being a\npotentially important feedback\".\nWe have edited this sentence to include the suggested modification.\nL487-490: “It is possible that a limited number of channels with high discharges and focused over\na small length of the grounding line could have the same impacts on the (de) stabilisation of ice\ndynamics as a scenario where many smaller channels occupy the entire grounding line, with\nchanges to the bed and sediment redistribution being a potentially important feedback.”\nFor 474 the list \"ice sheet, ocean, subglacial hydrology\" is great but also should include \" and \"the\nbed/subglacial geology/sediments\" pick your term but the feedbacks don't just include water and\nice.\nWe have made this change. The sentence now reads:\nL490-493: “Determining the nature of these relationships requires coupling between the ice sheet,\nthe ocean, subglacial sediments, and the subglacial drainage system and will provide crucial\ninformation for assessing past, contemporary, and future sea level contributions.”\nWhile this may seem like a minor suggest for editing, I personally think its critical that the reader\nof this important paper understand the real differences in ice sheet size/terminus position present\nin the geologic record but not yet accounted for in this work and that it be communicated feedbacks\nthat ultimately affect sea level do not end at the bed since the bed itself is also dynamic and part\nof the system\nWe would like to thank the reviewer for their in-depth comments and suggestions. We have revised\nour manuscript to incorporate the reviewer’s recommendations and believe that these edits have\nstrengthened the manuscript.\nReply to reviewers\nThe past, present, and future evolution of Aurora Subglacial Basin’s subglacial drainage system\nHayden, Pelle, and Dow\nReviewer comments are laid out in normal font, and our replies are in italic green font.\nReviewer #3 (Remarks to the Author):\nI wish to congratulate the authors on a very thorough job of revisions throughout this submission\nand review process. This latest version of their manuscript has satisfied all my earlier concerns\n(Reviewer #3) and has done so eloquently. I suggest that this mature version is ready for\npublication.\nSean Gulick\nWe would like to thank the reviewer for supporting our study and for their constructive feedback\nthat strengthened the manuscript."}, {"round": 3, "reviewer_comments": {"reviewer_3": "I wish to congratulate the authors on a very thorough job of revisions throughout this submission and review process. This latest version of their manuscript has satisfied all my earlier concerns (Reviewer #3) and has done so eloquently. I suggest that this mature version is ready for publication.\nSean Gulick"}, "author_response": "Reply to reviewers\nThe past, present, and future evolution of Aurora Subglacial Basin’s subglacial drainage system\nHayden, Pelle, and Dow\nReviewer comments are laid out in normal font, and our replies are in italic green font.\nReviewer #3 (Remarks to the Author):\nI wish to congratulate the authors on a very thorough job of revisions throughout this submission\nand review process. This latest version of their manuscript has satisfied all my earlier concerns\n(Reviewer #3) and has done so eloquently. I suggest that this mature version is ready for\npublication.\nSean Gulick\nWe would like to thank the reviewer for supporting our study and for their constructive feedback\nthat strengthened the manuscript."}]

Decision letter