Sedimentary basins reduce stability of Antarctic ice streams through groundwater feedbacks

Antarctica preserves Earth’s largest ice-sheet, which in response to climate warming, may lose ice mass and raise sea level by several metres. The ice-sheet bed exerts critical controls on dynamic mass loss through feedbacks between water and heat fluxes, topographic forcing, till deformation and basal sliding. Here we show that sedimentary basins may amplify critical feedbacks that are known to impact ice-sheet retreat dynamics. We create a high-resolution subglacial geology classification for Antarctica by applying a supervised machine-learning method to geophysical data, revealing the distribution of sedimentary basins. Hydro-mechanical numerical modelling demonstrates that during glacial retreat, where sedimentary basins exist, the groundwater discharge rate scales with the rate of ice unloading. Antarctica’s most dynamic ice streams, including Thwaites and Pine Island glaciers, possess sedimentary basins in their upper catchments. Enhanced groundwater discharge and its associated feedbacks are likely to amplify basal sliding and increase the vulnerability of these catchments to rapid ice retreat and enhanced dynamic mass loss. A machine-learning-based mapping of Antarctic subglacial geology suggests sedimentary basins lie beneath some of the most dynamic ice streams, increasing their vulnerability to rapid ice retreat.


1
School of Earth Sciences, The University of Western Australia, Perth, Western Australia, Australia. 2 The Australian Centre for Excellence in Antarctic Science, The University of Western Australia, Perth, Western Australia, Australia. 3 CSIRO Mineral Resources, Perth, Western Australia, Australia. 4 ARC Industrial Transformation and Training Centre in Data Analytics for Resources and the Environment (DARE), Sydney, New South Wales, Australia. 5 School of Biosciences, Geography and Physics, Swansea University, Swansea, Wales, UK. 6 School of Geography, Planning, and Spatial Sciences, University of Tasmania, Hobart, Tasmania, Australia. ✉ e-mail: lu.li@research.uwa.edu.au P alaeoclimate data indicate that mass loss of the Antarctic ice sheet is tightly coupled with warm climate forcing 1 that may cause a multi-metre contribution to sea-level change on centennial timescales 2 . Satellite observations show a continuous mass loss over the past four decades in both West and East Antarctica 3 . However, in predictions of future sea-level change, the spatial extent and magnitude of dynamic ice-mass loss in Antarctica remains highly uncertain [4][5][6] . This uncertainty reflects the identification of critical thresholds for dynamic mass loss, which are often associated with changing basal conditions 7,8 .
A sedimentary basin comprises the development of accommodation space into which accumulations of sediment are deposited that, over time, may be consolidated into sedimentary rock. A subglacial sedimentary basin may be defined by its fill and its characteristic geomorphology, which can be detected by geophysical techniques. Ice-sheet bed conditions in sedimentary basins differ from igneous or metamorphic basement because of their higher primary porosity, intrinsic permeability, lower relief surfaces, rheological stratification and increased susceptibility to glacial erosion 9 .
With 99% of Antarctica covered by ice, the understanding of subglacial geology relies on interpretations of geophysical data. The current understanding of sedimentary basins (Fig. 1) is derived from numerous individual studies (Supplementary Section 1.2). Direct constraints to define the extents of sedimentary basins are lacking in many areas. A systematic understanding of basins has also been limited by the diverse data and methods used and inconsistent mapping criteria between studies.
In this Article, we develop a sedimentary-basin likelihood map for Antarctica using the supervised machine-learning method random forest (RF) 10 . RF has proved to be a valid tool in lithology classification, demonstrating high predictive performance with limited training information 11 . We train the model using the current understanding of Antarctic sedimentary-basin and basement distributions ( Fig. 1). Due to variable spatial data density and interpretation uncertainty, training points are derived at random sample locations within 100 km cells. In each cell, a single training point may exist for the basin class, basement class or both (Supplementary Section 1.4.3). Input variables are sourced from continental-scale geophysical datasets including bed topography 12 , gravity field 13 , magnetic field 14 and their derivative products. Highly inaccurate data are relatively rare; however, errors in these input features may derive from interpolation errors due to variable data density and quality, and so areas with low confidence are masked (Supplementary Section 1.3.1, 1.4.1).
The likelihood map output from RF analysis (Fig. 2) shows the predominance of sedimentary-basin classification outcomes over basement classification outcomes. Likelihood >0.5 indicates a probable sedimentary basin. On this criterion, our map correctly predicts 90% of training information. RF results are stable using ten different training information sampling runs with 82% of area returning a consistent result across all runs. To assess the model robustness in areas lacking training information, tenfold block cross validation is applied by excluding training information from one fold. The classification returned for the untrained fold is consistent with the training information with 78% classification accuracy. Inconsistent classifications are associated with high variability in likelihood across runs and indicate locally complex geological settings either at the margins of basins or where small-scale variability may dominate on scales of tens of kilometres (Supplementary Section 1.5.3, 1.6).
Strikingly, the map defines sedimentary basins in some of Antarctica's most dynamic ice-sheet catchments, including the Amundsen Coast and Siple Coast regions in West Antarctica and Wilkes Land and the Recovery regions in East Antarctica.

East Antarctic basins
Most of East Antarctica is classified as basement but possesses some major basin regions. An extensive sedimentary basin is predicted in Antarctica preserves Earth's largest ice-sheet, which in response to climate warming, may lose ice mass and raise sea level by several metres. The ice-sheet bed exerts critical controls on dynamic mass loss through feedbacks between water and heat fluxes, topographic forcing, till deformation and basal sliding. Here we show that sedimentary basins may amplify critical feedbacks that are known to impact ice-sheet retreat dynamics. We create a high-resolution subglacial geology classification for Antarctica by applying a supervised machine-learning method to geophysical data, revealing the distribution of sedimentary basins. Hydro-mechanical numerical modelling demonstrates that during glacial retreat, where sedimentary basins exist, the groundwater discharge rate scales with the rate of ice unloading. Antarctica's most dynamic ice streams, including Thwaites and Pine Island glaciers, possess sedimentary basins in their upper catchments. Enhanced groundwater discharge and its associated feedbacks are likely to amplify basal sliding and increase the vulnerability of these catchments to rapid ice retreat and enhanced dynamic mass loss.
the southern Wilkes Subglacial Basin (WSB), transitioning to less extensive cover in the northern WSB (Fig. 3a). In western Wilkes Land, the inland region contains broadly distributed sedimentary basins (Fig. 2). The predicted sedimentary-basin distribution in these regions is consistent with differential subsidence during formation overprinted by spatially varying erosion from glacial retreat and advance cycles 15 . In the broader Recovery region, major basins are predicted in the upper catchments of the Bailey, Slessor and Recovery ice streams, linking to coastal basins (Fig. 3b). We also predict two distinct sedimentary basins near the South Pole: the Pensacola Pole Basin extending into the Academy and Support Force Glacier catchments and the South Pole Basin extending towards the WSB (Fig. 2). Elsewhere, many smaller-scale basins are predicted, several linked to the East Antarctic Rift System 16 and Lambert Rift.

West Antarctic basins
In the Ross Embayment of West Antarctica, our analysis resolves widespread basin coverage (Fig. 2). Predicted sedimentary basins continue, via the Siple Coast region, into the inland area of the Amundsen Sea Embayment. A transition from predicted broad basin coverage to more discontinuous coverage occurs in the Siple Coast region (Fig. 3c). In the Amundsen Sea, basement in Pine Island Bay transitions to sedimentary basin near the grounding line (Fig. 3d). The lower catchments of Pine Island and Thwaites glaciers are dominated by basement, but their upper catchments each have sedimentary-basin cover. Extensive sedimentary basins are predicted in the Weddell Sea sector beneath the Ronne-Filchner Ice Shelf and extending into the Robin Subglacial Basin (Fig. 2).

influence of sedimentary basins on ice-sheet dynamics
Boundary conditions of the ice-sheet bed, including till conditions, heat flux and hydrology, exert primary controls on ice-sheet dynamics 6 . The distinct properties of sedimentary basins allow key processes to occur that strongly influence these boundary conditions and increase the opportunity for instability-driven variations in ice-sheet dynamics. In sedimentary basins, the ice-sheet bed may be characterized by a hydro-mechanical stratigraphy incorporating a hard and impermeable basement beneath an overlying sedimentary sequence with weaker and more permeable strata. Unconsolidated sediments at the ice-sheet bed often comprise a rheologically stiff lower layer and a water-saturated flow-capable upper layer 17 that can facilitate basal sliding 18 . It is challenging to define the presence of such a feature at a regional scale 19   regions, thick and widespread coverage of glacial sediments is associated with the presence of sedimentary basins 20 . Continuous porous till can allow for fast ice flow accommodated by deformation of the ice-sheet bed, with sediment sources, hydrology and ice-loading processes all important 17,18 . Till continuity requires a consistent supply of sediment from either upstream erosion or in situ marine sediment deposition in previous ice-free periods 21 . In contrast to basement, sedimentary basins have structures that promote higher erodibility 9 , and so the presence of sedimentary basins in the upper catchment of an ice stream may favour a sustained supply of till to the lower catchment 22 .
Subglacial hydrology is a critical factor in ice-sheet dynamics, with wet-based glaciers showing faster flow and enhanced basal sliding compared with dry-based glaciers 8,23,24 . Subglacial hydrology also strongly impacts till conditions 25 . The presence of a permeable sedimentary basin at the ice-sheet bed introduces the potential for interaction between the basal water system and deep groundwater aquifer systems [26][27][28] . Excess water discharged from groundwater may weaken the till strength at the bed 25,29 , and upward-directed water flux may cause sediment liquefaction 30 . These processes reduce the coupling at the ice-sheet bed, potentially leading to increased ice-flow speeds. In addition, freshwater that flows across the grounding line may reduce ice-shelf stability with feedbacks that affect grounded ice dynamics 31 .
Ice-sheet retreat causes unloading (Extended Data Fig. 1) and is likely to increase groundwater discharge at the ice-sheet bed. We use a two-dimensional control-volume finite-element model 32 to test whether a permeable sedimentary basin in the upstream catchment causes a substantial increase in basal water-discharge rate relative to an entirely low-permeability basement. Furthermore, we investigate the scaling between the grounding line retreat rate and the basal water-discharge rate.
Our base scenario has a 3-km-thick sedimentary basin with basement, confining unit (clay and shale) and aquifer units (sandstone), with vertical permeabilities of 10 −19 , 10 −17 and 10 −15 m 2 , respectively 28 . Our simulations show that ice unloading during ice-sheet retreat causes excess groundwater discharge from sedimentary basins into the basal water system. We define the excess subglacial discharge rate (ESDR) as the difference in groundwater flux relative to groundwater flux in the control scenario with only the basement unit (Fig. 4). Increasing permeability and thickness of aquifer units facilitates higher ESDR during ice-sheet unloading, except for extremely high permeability, for example, in gravels or fractured rock aquifers (Supplementary Section 2.3.1).
With a grounding line retreat rate of 130 m yr −1 , ESDR reaches 1.96 mm yr −1 after 5.9 thousand years (kyr) (Fig. 4a). A scenario with a retreat rate of 325 m yr −1 , approximately that observed for Thwaites Glacier 33 , reaches a mean ESDR of 5 mm yr −1 after 2.3 kyr. Such discharge is an important contribution to the hydrologic budget. A model for Antarctica 34 indicates the mean basal melt rate for grounded ice is 5.3 mm yr −1 , and modelled basal melt rates in the upper Thwaites catchment are similar 35 . These rates are small compared with basal melt rates in fast-flowing glaciers 35 , which may reach 20-100 mm yr −1 . Crucially, very high retreat rates may have even more marked effects on groundwater discharge rate (Fig. 4c,d), with a fast retreat of the ice sheet (1,300 m yr −1 )     generating ESDR of almost 20 mm yr −1 after 0.6 kyr. These results show unloading-induced groundwater discharge may be an important source of subglacial water even in high basal melt-rate environments.
Water motion in the sedimentary basin advects geothermal heat and may transport heat from depth to the ice-sheet bed, with impacts on till and ice rheology and hydrology 26 . For example, with a 60 mW m −2 basal heat-flux boundary condition, we model additional heat flux due to groundwater discharge during the retreat phase increasing from 1 to 5 mW m −2 with faster retreat rate (Supplementary Section 2.3.3), leading to a potential enhancement of basal melt rate (0.1-0.5 mm yr −1 ).

Key vulnerabilities associated with sedimentary basins
Several ice streams in West Antarctica may be especially susceptible to rapid changes in flow 2,3,36,37 . Although East Antarctic ice-sheet behaviour has historically been relatively stable, observations from past warm periods and models of future change scenarios indicate a rapid dynamic response for several East Antarctic ice streams 38,39 . Many of these catchments possess extensive sedimentary basins, particularly in their upper catchments (Fig. 2).
The ice streams in the Wilkes and Recovery sectors are not currently associated with major thinning (Extended Data Fig. 2), but they are efficient retreat pathways and are key to the future vulnerability of the East Antarctic ice sheet 38 . In both cases, an extensive sedimentary basin is predicted in the upper catchment, with more spatially restricted basins in the lower catchment (Fig. 3a,b).
Recovery Ice Stream possesses an active lake system in its upper catchment 24 that may interact with a deep groundwater system. While identified lakes for the WSB are fewer, a high potential for large volumes of water to be released into the lower catchment is recognized 26 . The lower catchments of these regions have narrow sedimentary basins beneath major ice streams that are likely to be hydrologically active 40 . This basin arrangement suggests a high potential for inland sedimentary basins to influence regional ice dynamics through enhanced basal sliding associated with groundwater discharge.
The Siple Coast ice streams are notable for rapid transitions in ice streaming, including abrupt accelerations and slow downs associated with changes in subglacial conditions 37 . The ice-sheet state is controlled by till rheology variations coupled with subglacial hydrology 8 , with major contributions from groundwater systems 41 . Marine sediments were deposited in a broad sedimentary basin during previous deglaciations (Fig. 3c). In this case, basins are not extensive in the upper catchment but are a recognised factor for controlling the onset of fast flow 22 . In a similar setting, the Robin Subglacial Basin of the Weddell Sea defines subglacial conditions for the Institute Ice Stream, with fast ice flow and ice thinning occurring above a sedimentary basin 36 . The cold ocean water and buttressing effects of the Ross and Ronne-Filchner ice shelves may presently protect these regions from retreat. However, model scenarios of centennial glacial retreat show major thinning focused in these areas 39 , with the retreat organization being highly sensitive to subglacial hydrology variations 8,36 . Thinning-induced groundwater discharge is likely to be an important factor in these regions for organizing ice-sheet flow and potentially destabilizing the ice shelves.
In the Amundsen Sea sector, Pine Island and Thwaites glaciers are characterized by enhanced thinning and retreat driven by warm ocean water and marine ice instability processes 3,12,33,42 . The subglacial environment of these ice streams is complex and includes, in Thwaites Glacier, the transition from the upper catchment possessing a distributed subglacial hydrological system to the lower catchment dominated by a basal network of channels 43 . This transition is associated with increasing basal shear 43 and with the transition from sedimentary basin to basement 44 (Fig. 3d). In Pine Island Glacier, a sedimentary sequence is identified underlying the upper catchment and tributaries 45 associated with sedimentary basin (Fig. 3d). Downstream, the ice-stream trunk has widespread sediments at its bed, with basement beneath 45 . We suggest that upstream basins in both ice streams provide a consistent supply of sediment for a till-rich setting, as well as subglacial water and heat to influence ice-sheet dynamics through enhanced basal sliding.
The Pine Island and Thwaites Glacier grounding zones are associated with a transition from sedimentary basin inboard to basement outboard (Fig. 3d). Grounding line retreat may have stabilized at this transition due to a weak plastic bed extending ice thinning farther inland, rather than the focused thinning near the grounding line seen for a viscous bed 46,47 . Basal sliding allows for fast inland flow that stabilizes the grounding line 46 . With further forcing, retreat on the plastic bed into the already-thinned ice is faster than for a viscous bed with thicker ice, inducing a vulnerability to late-stage retreat exceeding 500 m yr -1 (ref. 7 ). With retreat-induced groundwater discharge, this is likely to be enhanced further (Fig. 4e).
Although basal sliding through deformation of soft sediment is possible in both basin and non-basin settings, the presence of thick and widespread sedimentary cover, coupled with the potential effects of groundwater discharge, mean that enhanced basal sliding is focused in catchments possessing sedimentary basins. Where other conditions also promote ice-sheet retreat and mass loss, especially if rapid, catchments with upstream sedimentary basins may exhibit excess subglacial groundwater discharge to the subglacial hydrology system. Excess groundwater discharge into the ice-sheet bed is coupled with ice-retreat rate. Slow to moderate retreat Ice grows over the initial 19 kyr and reaches glacial maximum extents of 1,300 km. The grey shaded area shows along-flow ice-sheet geometry during glacial maximum conditions. Vertical water flux at glacial maximum is mostly negative (blue areas), indicating groundwater recharge. From 20 kyr, the ice-sheet margin retreats at a linear rate indicated by the black dashed line, with positive vertical water flux (red areas) indicating groundwater discharge. We show sedimentary-basin geometry at the top. Aquifer unit (yellow) has horizontal permeability κ x = 10 −14 m 2 and vertical permeability κ z = 10 −15 m 2 ; the confining unit (light green) has κ x = 10 −16 m 2 and κ z = 10 −17 m 2 ; the basement (darker green) has κ x = κ z = 10 −19 m 2 . e, Excess subglacial water flux relative to basement for different retreat rates. The permeable bed promotes enhanced discharge from the upstream basin during ice-sheet retreat with discharge rate scaling with retreat rate.
scenarios generate groundwater discharge rates comparable to average basal melt rates and may be an important part of the overall budget 35 . In fast retreat scenarios, groundwater flux rapidly exceeds 10 mm yr −1 , potentially sustaining a high contribution to subglacial water even in high basal melt-rate environments. Excess water flux may cause substantial weakening of subglacial till, at the same time potentially increasing heat flux and basal melt rate, further enhancing basal sliding and inland dynamic thinning. This hydrogeological instability is a key feedback mechanism for several key catchments that focus the vulnerability of the Antarctic ice sheet, and its activation may contribute to accelerated mass loss in those catchments.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-022-00992-5.

Methods
Sedimentary-basin mapping with RF prediction. The RF prediction algorithm 10 is built on the basis of an ensemble of decision trees 50 . The randomness of RF is guaranteed by the uniqueness of each decision tree. Every tree in RF is constructed by a bootstrapped sample method (select nearly two-thirds of total training data with replacement) and grown by a random subset of variables at each split. For geological classification, each tree 'votes' for the geological class (sedimentary basin or basement), and the final result is assembled from all uncorrelated trees. In each location, the averaged vote represents the probability of each class being present. For the binary type classification problem in this study, a 0.5 likelihood is a natural boundary representing sedimentary basin or the basement. The R core 4.0.2 and package randomForest 51 was used to perform RF analysis.
We generate training information on the basis of the current understanding of the sedimentary-basin and basement distribution in Antarctica (Fig. 1). We choose several geophysical datasets as input features and then select the most important yet uncorrelated features to build the RF model. Sparse and spatial equally distributed training points are valid to construct RF model based on previous knowledge 11,52 . Here we subsample entire training information by a reference fishnet and then use the sampled training points to build a sub_RF model. The entire training information subsampling and RF training process are repeated ten times in ten different seeds to reduce bias during the random sampling process. We use the mean of all ten sub_RF model results to represent sedimentary-basin likelihood in Antarctica. Finally, we use tenfold block cross validation to validate the model, with a 78% classification accuracy.
We use the standard deviation of ten sub_RF models to assess the model uncertainty. We find major areas show low standard deviations, indicating a robust classification result. When comparing the cross-validation results with standard deviations, areas with low standard deviations show consistent predictions. Inconsistent predictions are associated with high standard deviations, which occur at geological boundaries or in complex geological settings.
For further details of the RF method, see Supplementary Section 1.
Hydro-mechanical model. Control-volume finite-element model CVFEM_Rift2D is a multi-physical modelling package including geomechanical modelling, hydrologic modelling, solute transportation and heat transformation 32 . Geomechanical modelling calculates the deformation and stress change during the ice-sheet advance and retreat. The mean normal stress-change rates act as a source-loading term in the hydrologic modelling. In our simulation, the geomechanical and hydrologic modelling are partially coupled. We use a coarse model (281 nodal columns and 11 nodal rows with 10 km intervals) to solve the lithosphere deformation and stress change associated with ice-sheet loading and unloading. The geomechanical modelling result is then interpolated into a refined hydrologic model domain to drive the groundwater flow. The fine-sized hydrologic model has a horizontal interval of 16 km and a vertical interval of 0.1 km. It gives the hydrologic model domain with 101 nodal columns and 71 nodal rows.
For the base case glacial retreat scenario, we simulate a 30 kyr glacial cycle with linear expansion and retreat over a 3-km-thick sedimentary basin. The 3-km-thick sedimentary basin includes five-layered hydro-mechanical stratigraphy with two aquifers confined by three confining units. The aquifer unit has horizontal permeability κ x = 10 −14 m 2 and vertical permeability κ z = 10 −15 m 2 . The confining unit has κ x = 10 −16 m 2 and κ z = 10 −17 m 2 . The basement has κ x = κ z = 10 −19 m 2 .
We simulate ice retreat in a specified period to investigate the relationship between ice-retreat rates and water-discharge rates. From 0 kyr to 19 kyr, the ice sheet grows from 0 to 3 km thickness, reaching 1,300 km lateral extent. The ice-sheet thickness remains stable during the glacial maximum for another 1 kyr. After glacial maximum, the ice sheet retreats to 0 km over a specified period, with the slowest occurring over 10 kyr and the fastest over 1 kyr.
We further investigate the impact of aquifer permeability and sedimentary-basin thickness on the water-discharge rates. In our simulation, the vertical permeability of aquifer units ranges from 10 −13 m 2 to 10 −17 m 2 . Sedimentary-basin thickness varies from 1 km to 4.5 km. We find increasing basal aquifer permeability and thickness facilitates higher groundwater discharge rates during ice-sheet unloading, although extremely high permeability may sink major amounts of basal water into the groundwater system.
More technical details and model parameters of hydro-mechanical modelling are in Supplementary Section 2.

Data availability
The results of RF classification and hydro-mechanical modelling can be accessed at https://doi.org/10.5281/zenodo.6611940.

code availability
The R notebook for RF classification can be accessed at https://doi.org/10.5281/ zenodo.6611940. The code of CVEFM_Rift2D can be accessed in the supplementary material from the original paper 32 .