Research ArticleGEOLOGY

Deglacial grounding-line retreat in the Ross Embayment, Antarctica, controlled by ocean and atmosphere forcing

See allHide authors and affiliations

Science Advances  14 Aug 2019:
Vol. 5, no. 8, eaav8754
DOI: 10.1126/sciadv.aav8754


Modern observations appear to link warming oceanic conditions with Antarctic ice sheet grounding-line retreat. Yet, interpretations of past ice sheet retreat over the last deglaciation in the Ross Embayment, Antarctica’s largest catchment, differ considerably and imply either extremely high or very low sensitivity to environmental forcing. To investigate this, we perform regional ice sheet simulations using a wide range of atmosphere and ocean forcings. Constrained by marine and terrestrial geological data, these models predict earliest retreat in the central embayment and rapid terrestrial ice sheet thinning during the Early Holocene. We find that atmospheric conditions early in the deglacial period can enhance or diminish ice sheet sensitivity to rising ocean temperatures, thereby controlling the initial timing and spatial pattern of grounding-line retreat. Through the Holocene, however, grounding-line position is much more sensitive to subshelf melt rates, implicating ocean thermal forcing as the key driver of past ice sheet retreat.


Concerns about the stability of the West Antarctic Ice Sheet (WAIS) have existed for decades and continue to be a topic of discussion today as oceanic and atmospheric warming pose a serious threat to the floating ice shelves that currently buttress the grounded ice (13). As the largest ice shelf in Antarctica, the Ross Ice Shelf (RIS) plays an integral role in buttressing both the West and East Antarctic ice sheets and thus has strong potential to control future sea-level rise. Future projections from Antarctic ice sheet models suggest that reduced buttressing resulting from RIS thinning promotes grounding-line retreat into the deep basin of West Antarctica, thereby increasing the Antarctic contribution to global sea levels (4, 5). In the past, the grounding-line in the Ross Embayment has shown substantial retreat of more than 1000 km from its Last Glacial Maximum (LGM) position, extending nearly to the continental shelf edge, to its present location (6, 7). However, the sensitivity of the ice sheet system to environmental forcings remains poorly constrained in this region. As a result, conflicting hypotheses have been proposed for both the timing and the pathway of grounding-line migration during the deglaciation, indicating either weak or strong sensitivity to early deglacial increases in air and ocean temperatures (8, 9).

The traditional view of Ross Embayment grounding-line retreat is the “swinging gate” model (8), in which the LGM grounding line only started to retreat during the Holocene, thousands of years after atmospheric and oceanic temperatures, as well as global mean sea level, began to rise. This interpretation envisaged the most pronounced retreat occurring on the western side of the basin, synchronous with changes in the central embayment, and thus chronologies obtained from the Transantarctic Mountain front are reliable indicators of broader scale retreat in the Ross Sea. More recently, it has been posited that the grounding line retreated substantially earlier and more extensively than previously estimated (9), which would suggest considerably higher sensitivity to early deglacial environmental forcing. In this scenario, the grounding line retreated beyond its modern position along Siple Coast and readvanced during the Holocene. Others have also argued for more moderate and dynamic models of retreat, with a combination of oceanic and bathymetric controls driving a landward grounding-line migration from the central embayment (1014). This apparent conflict remains difficult to reconcile as few age constraints exist, and most are biased toward the western margin of the Ross Sea.

In terms of the drivers of retreat, differences in the seafloor bathymetry and substrate conditions likely resulted in asynchronous retreat behavior in the eastern, central, and western troughs of the Ross Sea (12). Although the specific role of climate remains relatively less understood, it has been noted that the ice shelf may be highly sensitive to warm atmospheric conditions and the incursion of warm ocean currents onto the continental shelf (15). Between the LGM and the Early Holocene, Antarctic climate and Southern Ocean conditions experienced substantial millennial-scale changes that likely affected Antarctic ice sheet mass balance, including the last glacial termination, Heinrich event 1, the Antarctic Cold Reversal (ACR), and the potential influences of the Meltwater Pulse 1a and 1b (MWP-1A and MWP-1B) events (10). However, it remains challenging to assess the effects of such climate changes on grounding-line retreat in the Ross Sea given the lack of age control from the central embayment (7) and because smaller-scale subcatchments, from where age control data may exist, can behave independently of regional grounding-line retreat (16).

Additionally, the precise magnitudes and timing of climate changes are difficult to determine on a regional scale given the sparse spatial and temporal coverage and large uncertainties of available paleoclimate proxy records. In particular, changes in subsurface ocean temperatures in the Ross Sea remain poorly constrained because of the poor preservation of carbonate in Southern Ocean sediments (17). Ice sheet modelers are tasked with selecting between climate forcings that are imperfect owing to these limitations of available proxy records, which may not be representative of the area of interest, and transient climate model simulations, of which there are relatively few and which may also have regional biases (18). To resolve this issue, here, we perform regional ice sheet model simulations using a wide range of climate forcings derived from various Antarctic ice core and marine sediment records and from two transient climate model simulations to investigate the sensitivity of the grounding-line position to changes in atmosphere and ocean conditions and to better constrain both the timing and spatial pattern of grounding-line retreat and ice shelf formation in the Ross Embayment. Comparisons of the simulations to the geologic record are, in turn, used to inform which deglacial climate forcings are the most reasonable for the Ross Embayment.

The regional ice sheet simulations are performed at 10-km resolution using the Parallel Ice Sheet Model (PISM), a sophisticated ice sheet model that allows for realistic ice streams that exhibit the full range of observed ice stream velocities (19, 20). As a result, it is one of the few ice sheet models capable of multimillennia simulations of ice sheets containing freely evolving ice streams. We consider a domain of the Ross Embayment in which the ice shelf and ice drainage basins of the grounded ice are allowed to dynamically evolve, and a larger area around the drainage basins is maintained at a constant thickness as a time-independent boundary condition of the model. This approach allows for higher computational efficiency due to the smaller domain size, while simultaneously allowing for the grounding-line position to change through time. Details of the regional model setup and initialization procedure are described in Materials and Methods.

For the deglacial simulations, which are run from 35 to 0 thousand years (ka), we apply combinations of the output of two transient climate model simulations, namely, LOVECLIM (LOch–Vecode-Ecbilt-CLio-agIsm Model) (21) and TraCE-21ka (22), the temperature reconstructions of the EPICA (European Project for Ice Coring in Antarctica) Dome C (EDC) ice core from East Antarctica [75°S, 123°E; (23)] and the WAIS Divide (WDC) ice core from West Antarctica [~79°S, 112°W; (24)], global and Southern Ocean benthic ocean temperature proxy reconstructions [GOT (global δ18O stack) and SOT (Southern Ocean δ18O)] (25, 26), and model-proxy averages (MPAs) as climate forcings, for a total of 36 deglacial climate forcing experiments (Table 1). The climate model simulations were performed with evolving boundary conditions to reproduce global climate of the last deglaciation (21, 22). The ice core temperature reconstructions are determined from isotopic measurements and also cover the full deglacial period. The GOT reconstruction is based on a compilation of 57 deep ocean records (25), whereas the SOT reconstruction is more specific to the Southern Ocean, derived from a marine sediment core from offshore New Zealand (26). The MPA atmosphere and ocean forcings are simple averages of each temperature proxy reconstruction and temperature output of both climate models. Additional details of the climate forcings are provided in Materials and Methods.

Table 1 Atmosphere-ocean forcing combinations applied to the ice sheet model.

The top row and the leftmost column, respectively, refer to the six atmosphere forcings and the six ocean forcings applied to the ice sheet model. In total, 36 simulations of the ensemble are performed, with each table entry referring to an individual simulation. The atmosphere forcing is applied in the form of surface air temperature (SAT), precipitation (7.0%/°C), and sea ice resistance anomalies to a regional climatology from the TraCE-21ka and LOVECLIM transient climate simulations, the EDC and WDC ice core records, and an MPA of each atmosphere forcing. The ocean forcing is applied in the form of subshelf mass flux (melt rate) offsets based on the climate model simulations, two benthic ocean temperature (OT) reconstructions, and an MPA (LOVECLIM 400m T and two benthic OT reconstructions). We define three scenarios, namely, warm, moderate, and cool, based on LGM SAT anomalies and timing and amount of early deglacial ocean warming (Figs. 1 and 2). SST, sea surface temperature.

View this table:

The climate forcings used in the ice sheet model simulations include both atmospheric and ocean components. For the atmosphere forcings, we uniformly apply surface air temperature (SAT) and paleoprecipitation based on the Clausius-Clapeyron relationship to temperature (7%/°C) as anomalies to a modern climatology (27). We also apply a heuristic time-varying back pressure (calving resistance) offset, intended to mimic the resistive effects of iceberg mélange or dense sea ice that are not specifically modeled (see Materials and Methods). Air temperature affects calving by influencing ice hardness, the presence of surface meltwater on ice shelves, the growth of sea ice, and the binding of iceberg clasts; hence, this back pressure forcing is an attempt to modulate calving based on atmospheric controls in a similar, though more simplistic, fashion to previous ice sheet modeling studies (5, 28, 29). The back pressure offset is applied as a scalar fraction that is included in the stress boundary condition of the ice shelf. The ocean forcing is in the form of scalar subshelf mass flux (melt rate) offset. The offset is a positive or negative basal mass flux (ice loss or gain, respectively) within the mass continuity equation and is applied uniformly to the base of the ice shelf, as in Golledge et al. (10). The subshelf melt rate and back pressure offsets are determined using scaling relationships between reconstructed LGM and present-day (PD) ocean temperature and SAT and applied uniformly as anomalies to PD conditions (see Materials and Methods). Last, we apply the same sea level forcing in each model run, which is based on statistical analysis of the global sea-level proxy record (3033). In addition to the main model ensemble runs, we perform an experiment using a higher mantle viscosity and a more aggressive ocean thermal forcing based on the WDC record, as in Kingslake et al. (9) (see Materials and Methods).


Sensitivity of grounding-line retreat to climate forcing

The differences between the climate forcings allow us to assess the sensitivity of grounding-line retreat to changes in atmosphere and ocean conditions. The surface temperature and precipitation forcings have a wide range of LGM and early deglacial anomalies (e.g., LGM anomalies of −7° to −13°C and −49 to −91% relative to PD, respectively), and differ substantially in terms of the timing and magnitude of the ACR, but are all relatively similar throughout the Holocene (i.e., differences <2°C and 14% from 11.7 to 0 ka, respectively; Fig. 1A). The subshelf melt rate and back pressure forcings vary most during the MWP-1A and ACR events, respectively (Fig. 1, B and C), but are all scaled to the same LGM anomaly (see Materials and Methods) and are likewise similar throughout the Holocene. The subshelf melt rate forcings are based on the modeled and reconstructed subsurface ocean temperatures. The models and proxy reconstructions suggest that MWP-1A generally has a subsurface warming effect; however, the degree and duration of this warming vary considerably among the forcings. We also consider the case of ocean cooling through the ACR (dark blue line in Fig. 1B), using the near-surface temperature of TraCE-21ka. The back pressure forcings are similarly applied as a scalar offset based on the SAT anomalies; hence, the forcings decrease from the LGM, a period of higher sea ice concentration, to PD (no offset). All simulations are performed with the same sea-level forcing (Fig. 1D).

Fig. 1 Climate forcings and modeled grounded ice volume.

Time series of the atmosphere forcings of (A) SAT and precipitation anomalies (°C relative to PD; % relative to PD) and (B) back pressure fraction offsets (λ; see Materials and Methods), (C) ocean thermal forcings of subshelf melt rate offsets (m year−1 relative to PD), (D) sea-level forcing (m), and (E) modeled grounded ice volume (m3). In (A), all time series correspond to both y axes because we use a temperature-precipitation scaling relationship for the precipitation forcing (7% °C−1). The blue bars in the panels correspond to the approximate timing of MWP-1A and MWP-1B, and the pink bar corresponds to the ACR. In (E), the six simulations forced with the TraCE-21ka SST (LOVECLIM) subshelf melt rate forcing in (C) are shown in dark blue (dark red). All other 24 simulations in the main model ensemble are shown in gray. BP, before the present.

The decline of grounded ice volume initiates between 15.5 and 13.2 ka in the model simulations (Fig. 1E). The model spread of grounded ice volume is largest through this interval of initial retreat and into the early Holocene as the ice shelf begins to form (~10 to 8 ka). Simulations forced with the relatively warmest ocean thermal forcing (i.e., LOVECLIM) show enhanced grounded ice volume decline ~14.5 to 13.8 ka, roughly synchronous with the timing of MWP-1A (red lines in Fig. 1E). In contrast, simulations with the relatively coolest ocean thermal forcing [i.e., TraCE-21ka sea surface temperature (SST)] display a rapid decrease in grounded ice volume ~11.2 to 10.5 ka, roughly synchronous with the timing of MWP-1B (blue lines in Fig. 1E). Simulations forced with ocean thermal forcing derived from the proxy record and the MPA generally exhibit more gradual decreases in grounded ice volume through this interval of increasing sea level (gray lines in Fig. 1E).

The atmosphere forcings (i.e., surface temperature, precipitation, and back pressure) exert a strong modulating effect on the sensitivity to the subshelf melt rate offset. In particular, simulations forced with cooler LGM SAT anomalies (e.g., TraCE-21ka SAT forcing) have higher ice hardness, lower temperatures at the bed, and decelerated ice flow, which makes the ice sheet more resistant to the increase in ocean thermal forcing from MWP-1A and in the Early Holocene, thereby causing a later retreat. The opposite is true for simulations forced with warmer LGM SAT anomalies (e.g., LOVECLIM-1), which are more responsive to ocean thermal forcing, and thus exhibit an earlier retreat. Similarly, steeper gradients of decreasing back pressure in the early deglacial period (e.g., WDC back pressure forcing) allows for higher surface ice flux, also inducing an earlier retreat. Simulations forced with relatively cooler (warmer) atmosphere forcings paired with stronger (weaker) ocean thermal forcing behave more similarly to those forced by moderate ocean thermal forcings (Fig. 1D). The sea-level forcing also contributes to the decline in grounded ice; however, in the absence of climate forcing, it is unable to drive ice sheet retreat (fig. S1).

In terms of the spatial pattern of grounding-line retreat, we consider individual scenarios: warm (LOVECLIM-1 atmosphere/LOVECLIM ocean), moderate (MPA atmosphere/MPA ocean), and cool (TraCE-21ka atmosphere/TraCE-21ka SST ocean) scenarios (Fig. 2, A to C). These scenarios are defined on the basis of LGM SAT anomalies and timing and amount of early deglacial ocean warming. The model bathymetry is based on the Bedmap2 dataset (34), although the bed is smoothed to improve computational efficiency using the bed roughness parameterization of Schoof (35) (Fig. 2D). In all three scenarios, the earliest retreat occurs in the central Ross Sea, and retreat pathways arising from the different forcings are generally consistent in this area. Some spatial variations exist in the eastern and western basins, with the cooler scenario displaying a delayed retreat in these basins relative to the central Ross Sea. The retreat behavior resembles the “saloon-door” and marine-based models of retreat proposed by McKay et al. (11), Halberstadt et al. (12), and Lee et al. (13), in which the grounding-line retreat initiated in the central embayment, with landward migration toward the modern East Antarctic Ice Sheet outlet glaciers. However, the warm scenario initially does exhibit earlier retreat in the eastern and western embayment. Each scenario shows ice remaining grounded on Pennell Bank and Crary Bank into the Holocene (~11.7 to 0 ka). The moderate scenario displays the most gradual retreat, whereas the warm (cool) scenario experiences periods of rapid retreat ~15 to 13 ka (~12 to 10 ka). Once the grounding line has retreated beyond Ross Island, each scenario shows a similar retreat pattern, with a steady progression along the Transantarctic Mountains toward the Siple Coast.

Fig. 2 Grounding-line retreat scenarios.

(A) Warm (LOVECLIM SAT/LOVECLIM 400m OT), (B) moderate (MPA SAT/MPA OT), and (C) cool (TraCE-21ka SAT/TraCE-21ka SST) grounding-line retreat scenarios. The black squares indicate marine sediment core locations with minimum age constraints of grounding-line migration from radiocarbon dating. The thin black line marks the grounding-line position in the eastern Ross Sea estimated by Bart et al. (14). (D) Smoothed bed elevation (m a.s.l.) used in the model simulations based on Bedmap2 (34). The gray lines show elevation contours of 200 m. The black line indicates the modern grounding-line position.

The timing of ice shelf formation varies considerably among the simulations and is related to the climate forcing, with an overall model spread of ~2 ka. However, the model simulations are consistent in terms of the process of ice shelf formation. As the grounding line retreats toward the Drygalski Trough in the western Ross Sea, the ice becomes ungrounded, and ice flow from South Victoria Land accelerates (Fig. 3A). The formation of the main ice shelf initiates as the grounding line retreats toward the Byrd ice stream (Fig. 3B). When the faster-flowing ice stream is no longer inhibited by thick grounded ice, its velocity increases, and the floating ice shelf begins to expand over an area of deeper topography. Grounding-line migration continues to propagate southwestward along the Transantarctic Mountains and eastward toward the Siple Coast.

Fig. 3 Ice shelf formation.

(A and B) Ice surface velocity (m year−1) before and after ice shelf formation of the moderate scenario (Fig. 2B). The modeled grounding-line position is shown in black. For context, the modern grounding (calving) line position is shown in pink (purple). IS, ice stream. (C) Floating ice shelf area (m2) of simulations forced with the MPA atmosphere forcing and different ocean forcings (orange lines) and simulations forced with the model-proxy ocean forcing and different atmosphere forcings (purple lines). The timing of (A) and (B) are indicated by the blue and pink bars in (C), respectively.

As the modeled ice shelf area reaches near-modern levels (11 to 7 ka), the ocean thermal forcing dominates in terms of the model spread, with small differences in the subshelf melt rate resulting in large variations in ice shelf area (Fig. 3C). With the development of the ocean cavity underneath the ice shelf, the subshelf melt rate forcing acts on a larger area and becomes the main driver of grounding-line retreat. This effect increases throughout the Holocene, even as the differences between the ocean thermal forcings decrease (Fig. 1C). The differences in ice shelf area between model simulations are primarily observed along the grounding line of the Siple Coast and the calving line rather than the grounding line along the Transantarctic Mountains.

On the basis of the ensemble of climate-forcing simulations, we can explore how the sensitivity to climate-forcing selection evolves both temporally and spatially by assessing model agreement in parameters such as ice thickness (Fig. 4). For example, if a given model grid cell is covered by thick grounded ice in some climate-forcing simulations at a specific time step but is covered by thinner floating ice shelf or open ocean in other simulations at the same time step, the SD of ice thickness in the ensemble as a whole will be higher at this point in time and space, indicating higher sensitivity to the climate-forcing selection. During the ACR (~14 ka; yellow line in Fig. 4), a large area of model disagreement (high SD of ice thickness), covering the western Ross Sea and Ross Island region, is observed, driven by the variation in the timing and spatial pattern of the initial retreat.

Fig. 4 Uncertainty due to climate forcing.

SD of modeled ice thickness (m) in the Ross Embayment of the 36 climate-forced deglacial simulations (Table 1). The grid color indicates the time-averaged uncertainty of a given grid cell through the interval of 18 to 4 ka. The colored lines outline the areas that include the top 30% SD in ice thickness at a given time slice in 2-ka intervals from 14 to 6 ka. EAIS, East Antarctic Ice Sheet.

The area of highest SD migrates southward through the Early Holocene on the western side of the embayment, toward the Transantarctic Mountains. The SDs reach peak values ~12 ka (red line in Fig. 4) and then decrease through time. As the climate forcings converge in the Holocene, the simulations show higher agreement (i.e., lower SDs of ice thickness), and the overall area of high uncertainty (top 30% SD) decreases and becomes concentrated along the grounding and calving lines.

Comparison of model simulations to the geologic record

Few geological constraints exist for the timing of grounding-line retreat in the Ross Sea given the poor preservation of carbonate for radiocarbon dating and large uncertainties in dating of glacimarine sediments (17). However, in the eastern Ross Sea, marine sedimentary facies succession and foraminifera-based radiocarbon chronology show initial grounding-line retreat by as early as 14.7 ka, but the grounding line stabilized on the outer continental shelf until ~11.5 ka (14). In Fig. 2, the warm scenario is the closest to capturing this early initial retreat in this area, but the cool scenario better matches the grounding-line stagnation and rapid Early Holocene retreat.

In the western Ross Sea, radiocarbon chronologies suggest open ocean conditions at Coulman High by 8.6 ka (11). This region in the western Ross Sea exhibits the highest sensitivity in the domain to climate forcing, as demonstrated in Fig. 4. However, none of the model ensemble simulations conflict with these minimum age constraints in the western Ross Sea (Fig. 2), with the grounding line reaching Coulman High by as early as 12 ka. The simulations forced with relatively cooler atmosphere/ocean forcings display the most similar timing to these estimated minimum ages (e.g., ~11 ka).

Ice thinning histories derived from cosmogenic nuclide surface-exposure records of marginal sites in the Ross Embayment generally show enhanced thinning rates in the early to Mid-Holocene, although the precise timing differs among individual glaciers over a large time range of ~14 to 3 ka (3641). Comparisons between ice sheet model simulations and surface-exposure records are difficult because of limitations in model resolution, which may not be sufficient for resolving smaller-scale features. With this caveat in mind, we compare ice thickness changes of the model ensemble in four regions from which cosmogenic nuclide surface-exposure records exist: the Northern Victoria Land (NVL) region, the McMurdo region, the Southern Transantarctic Mountain (STAM) region, and the northern Marie Byrd Land (MBL) region (Fig. 5). We use regional averages since individual model grid cells, which cannot resolve individual glaciers, may not be fully representative of the regional-scale ice sheet evolution.

Fig. 5 Model-data comparison.

Regional ice thickness anomalies relative to PD of the model simulations for the (A) NVL, (B) McMurdo, (C) STAM, and (D) MBL regions versus ice thickness anomalies of cosmogenic nuclide surface-exposure records of glaciers within each region. Individual model simulations are shown in gray, and the ensemble average is shown in black. Data for Tucker and Aviator glaciers were obtained from Goehring et al. (41). Data for Mt. Discovery and Mackay Glacier were obtained from Anderson et al. (40) and Jones et al. (38), respectively. The STAM glacier data were obtained from Spector et al. (39) and Todd et al. (37), and those of Mt. Rae were obtained from Stone et al. (36). All surface-exposure ages were recalculated following recent improvements to global production rate (53). Age uncertainties of the records are indicated by the horizontal bars. The inset in (A) shows the region locations (black box outline) and the glacier locations overlaid on the ensemble average of the final modeled RIS configuration (yellow shows floating ice, blue shows grounded ice, the thick black line shows modeled grounding-line position, white contours show surface elevation in 500-m intervals, the pink line shows observed modern grounding-line position, and the purple line shows observed modern calving line position). In the panels, the green line indicates the early retreat/Holocene readvance scenario (see Materials and Methods). CCR, Crater Cirque; SKF, Skua Basin; TAM, Transantarctic Mountain.

Even within these regions, surface-exposure ages of individual glaciers record differences in the timing of ice thinning of >1 ka. This is especially pronounced in the STAM region, where ice thinning of Beardmore Glacier precedes that of Scott Glacier by ~5 ka, and Reedy Glacier, based primarily on dated drift limits, shows relatively less overall thinning (Fig. 5C). The timing of ice thinning in the model ensemble runs generally occurs between that estimated by the Beardmore and Scott glaciers (i.e., ~13 to 10 ka), with simulations forced by relatively warmer (cooler) ocean/atmosphere forcings showing stronger agreement with Beardmore (Scott) Glacier.

In the NVL and McMurdo regions, most of the simulations exhibit precipitation-driven ice thickening in the early deglacial period (20 to 14.5 ka). Regional ice thinning occurs ~14 to 10 ka, with a substantial model spread, indicating that these regions are highly sensitive to the climate forcing selection (Fig. 5, A and B). A similar timing of ice thinning is observed for the MBL region (i.e., ~14 to 10 ka), although the rate of thinning is higher (Fig. 5D). Assessment of model performance is difficult in these regions given the lack of ice thickness estimations through the key time interval of 14 to 10 ka. The model ensemble shows similar ice thickness anomalies in the early deglacial period to those estimated at Tucker Glacier, Mackay Glacier, Mt. Rae, and Mt. Discovery. The simulations forced with relatively cooler atmosphere/ocean forcings are also generally consistent with Aviator Glacier and Mt. Discovery in the Early Holocene.

However, model-proxy mismatches of ice thickness do occur, as no simulation captures the increase and subsequent decrease in ice thinning estimated in the McMurdo region in the Mid-Holocene (~8 to 6 ka) or the decrease at Mt. Rae in the Late Holocene (~4 to 2 ka). The NVL glaciers also record a more gradual ice thinning in the early to Mid-Holocene than indicated by the model simulations. As explained above, some of this mismatch may be related to glaciers recording local signals instead of the broader regional signals. For example, the recorded rapid thinning of Mackay Glacier is likely the result of local retreat into an overdeepened basin rather than reflective of a broader regional deglaciation at this time (38). A model resolution of 10 km is unlikely to resolve such changes from local topography and may partly explain why the exact timing of recorded ice thinning at these sites lags the simulated regional-scale changes. In addition, the application of uniform precipitation and ocean forcing does not account for changes in atmospheric and oceanic circulation, which may be significant at the subcatchment scale, thereby contributing to the model-data mismatch.

Least consistent with the surface-exposure ages and marine radiocarbon dates of the Ross Sea is the simulation using higher mantle viscosity and the WDC record as a subshelf melt rate forcing (green lines in Fig. 5, A to D), paired with an average atmosphere forcing. In this simulation, which reproduces the early retreat and Holocene isostatic rebound-driven readvance proposed by Kingslake et al. (9), the ice sheet is highly sensitive to the subshelf ocean forcing, reaching a near-modern subshelf melt rate by ~15 ka (see WDC forcing in Fig. 1A). At this point, precipitous ice thinning occurs in all three regions, considerably earlier than estimated by any of the surface-exposure chronologies.


The deglacial simulations demonstrate that the ocean forcing is the primary control on the timing of grounding-line retreat in the Ross Embayment, an effect that is strongly modulated by the atmosphere forcing before ice shelf formation and the development of the ocean cavity. The greatest variation in the timing and spatial pattern of retreat, that is, the greatest sensitivity to the climate forcing selection, occurs during the period between the initial decline in grounded ice volume and the early formation of the ice shelf (i.e., ~14 to 8 ka). The retreat pathways of the climate forcing experiments are relatively consistent, suggesting that the seafloor bathymetry acts as a primary spatial control on the modeled grounding-line migration. In each simulation, the earliest retreat occurs in the central basin relative to the eastern and western basins, in agreement with the saloon-door model of retreat, most recently suggested in Lee et al. (13). The spatial variations that do exist are driven by differences in the relative timing of retreat in each basin, which is influenced by the ocean and atmosphere forcings, implicating these external forcings as secondary controls.

The retreat pathway observed in the simulations is in agreement with the asynchronous and dynamic behavior suggested by Halberstadt et al. (12), Lee et al. (13), and Greenwood et al. (16) based on multibeam bathymetry data of the Ross Sea; however, the model is limited in capturing some smaller-scale features. The simplified ocean forcing applied in these simulations likely contributes to this data-model mismatch, as it neglects the effects of ocean circulation, eddies, turbulence and focusing of Circumpolar Deep Water, and changes in continental shelf water mass formation processes, all of which is consequential considering the importance of these processes to basal melt rates (42). Although the climate models account for such circulation changes, we use an average temperature anomaly of the entire Ross Sea to apply the ocean thermal forcing to the domain (see Materials and Methods); thus, differences in subsurface warming and cooling between the eastern and western Ross Sea are not represented. This could explain why most of the simulations display a later grounding-line retreat than indicated in the Whales Deep basin in the eastern Ross Sea (14). In terms of the atmosphere forcings, the temperature-precipitation scaling relationship may not be appropriate over the full domain given the influence of atmosphere dynamics and synoptic-scale processes in millennial-scale changes in snow accumulation over West Antarctica (43). This inhibits comparisons to specific glaciers, as changes in accumulation that are not captured by the precipitation forcings may be relevant for glacier mass balance changes.

Additionally, given the abundance of ice rises in the western Ross Sea, which can act as pinning points that stabilize the grounding line (12, 44), it is conceivable that local grounding-line retreat occurred later than modeled in this region because of model resolution limitations, which result in a smoother than realistic seafloor bathymetry. It has recently been suggested that ice flow switching that occurred as the ice sheet retreated may have resulted in a Holocene readvance in the western Ross Sea (16). Although the simulations show enhanced flow and direction changes from South Victoria Land, this does not reproduce the Holocene ice sheet advance. These discrepancies between the simulations and observations are likely driven by resolution limitations and differences in seafloor substrates that are not adequately represented in the model but are important for till deformation and ice flow velocities. Despite these potential local biases, these simulations effectively demonstrate the climate sensitivity of grounding-line retreat on the broader regional scale. Grounding-line retreat displays the highest sensitivity to the climate forcing selection in the western Ross Sea, both spatially and temporally; therefore, establishing reliable age constraints in this region remains a high priority. From a modeling perspective, the roles of model resolution and spatial heterogeneity in substrate type and ocean and precipitation forcing require further exploration to understand and improve upon model-data discrepancies.

While most records from the Ross Embayment suggest the most substantial terrestrial ice thickness changes occurring in the early to Late Holocene, Kingslake et al. (9) recently proposed an early deglacial retreat scenario based on the presence of finite radiocarbon ages measured in the organic carbon from subglacial sediments. In this scenario, the grounding line retreated beyond the modern-day position along Siple Coast in the Early Holocene and readvanced to its current position because of isostatic rebound. In model simulations also using PISM, they simulate the furthest retreat in the Ross sector at 9.7 ka. Using a higher mantle viscosity term and higher early deglacial subshelf melt rates, we are able to reproduce a similar timing of ice sheet retreat and readvance as well as minimum grounding-line extent (9.8 ka) within our regional domain (fig. S2). However, we demonstrate here that this simulation is difficult to reconcile with the existing paleorecords because, in this scenario, the modeled TAM/MBL ice thinning and open ocean conditions in the Ross Sea occur substantially earlier than the empirical data suggest. While this does not eliminate the possibility of ice sheet retreat and readvance driven by isostatic uplift, the timing of retreat likely occurred later than suggested by Kingslake et al. (9), suggesting a less sensitive response to climate forcing. We also demonstrate here that the deglacial simulations exhibit high sensitivity to small changes in subshelf melt rates along the Siple Coast during the Holocene. Therefore, an additional mechanism of Siple Coast retreat/readvance that should also be considered is that incursions of warm (cool) ocean water into the RIS cavity in the Mid- to Late Holocene could contribute to more extensive retreat (advance). The effects of both mechanisms should be further explored with regard to the grounding-line position along the Siple Coast.

Geologic evidence and Antarctic ice sheet modeling indicate that the Ross Embayment was not a substantial contributor to MWP-1A (10), despite the continent as a whole most likely having made a multimeter contribution to this event (45). However, recent radiocarbon dates suggest that initial grounding-line retreat in the eastern Ross Sea may have occurred by as early as 14.7 ka (10). Differences in subsurface ocean warming associated with a large meltwater flux from the Antarctic Peninsula and the mid-to-outer Weddell Sea may explain how the grounding line retreated in Whales Deep Basin, with no coincident retreat observed in the western Ross Sea. In fact, Golledge et al. (10) demonstrate that Southern Ocean freshwater forcing in LOVECLIM results in subsurface ocean warming in the eastern and central Ross Sea and cooling in the western Ross Sea. TraCE-21ka likewise displays enhanced subsurface ocean warming in the eastern and central Ross Sea relative to the western Ross Sea resulting from MWP-1A freshwater forcing (fig. S3). It is therefore possible that in terms of the initial grounding-line retreat, the eastern Ross Sea followed a “warm” ocean scenario trajectory, potentially resulting from subsurface ocean warming associated with Heinrich event 1 and/or MWP-1A, as indicated in the climate simulations (10, 21, 22). In contrast, the western Ross Sea followed a “cool” ocean scenario trajectory, with a timing of retreat more consistent with MWP-1B. In the case of the eastern Ross Sea, the effect of the MWP-1A–related ocean warming was likely short-lived, however, as the grounding-line position remained relatively unchanged following the initial retreat until the Early Holocene (14).

We note that the atmosphere forcings that yield the best fit to the geologic record are derived from climate model simulations rather than from paleoclimate proxy reconstructions. In particular, the relatively cooler SAT and higher back pressure forcings of TraCE-21ka and LOVECLIM-2 are important for delaying the timing of retreat by reducing both the surface flux and the sensitivity to elevated subshelf melt rates. In the TraCE-21ka and LOVECLIM simulations, the continental margins and coastal seas exhibit some of the greatest temperature anomalies due to changes in sea ice concentration. The SAT anomalies of these climate model simulations range from 2.3° to 3.3°C cooler at 18 ka than the more continentally located EDC and WDC ice cores. Few marine ice cores exist in the region to determine if the climatic effects of sea ice changes are more strongly experienced along the continental margin of the Ross Embayment; however, the δ18O records of EDC and Talos Dome display strong consistency, despite the more coastal location of the latter (46).

Recent evidence from diatom oxygen isotopes, a proxy for ice discharge, and climate model simulations have implicated ocean thermal forcing as the main driver of Antarctic ice sheet retreat in the Holocene (47). We similarly show that changes in subshelf melt rates become increasingly important relative to the atmosphere forcings through the Holocene in terms of the grounding-line position in the Ross Embayment, particularly along the Siple Coast. This is notable given concerns for marine ice sheet instability and modern observations of dynamic thinning and grounding-line retreat of WAIS glaciers attributed to increased basal melt of ice shelves (48). The Siple Coast has experienced such changes in recent years, as the grounding-line position of the Kamb Ice Stream has shown rapid retreat related to post-stagnation thinning (49). The simulations here likewise suggest that the grounding line of the largest Antarctic ice shelf is potentially vulnerable to elevated subshelf melt rates from warm water incursions into the RIS cavity.


Past grounding-line retreat in the Ross Embayment has been substantial; however, the sensitivity of grounding-line dynamics to atmosphere and ocean forcing has remained poorly constrained in this area. To compensate for limitations in available paleoclimate proxy records and climate simulations, our regional ice sheet simulations use a wide range of climate forcings to highlight their effect on grounding-line migration. Differences between the ocean and atmosphere forcings are consequential for both the timing and spatial pattern of grounding-line retreat, particularly in the western and eastern Ross Sea, and therefore have implications for the sea-level contribution from the Ross sector over the last deglaciation as well as the overall deglacial climate evolution of this region. The model results indicate that early deglacial atmospheric conditions can enhance or diminish ice sheet sensitivity to rising ocean temperatures, thereby controlling the initial timing and spatial pattern of grounding-line retreat. Simulations using warm ocean forcing better match the early onset of retreat in the eastern Ross Sea indicated by marine radiocarbon dates, but those using cool atmosphere and ocean forcings generally reproduce the best overall agreement with the geological constraints of the grounding-line retreat in the western Ross Sea and terrestrial ice thinning. These cool climate scenarios show retreat initiating from the central embayment and enhanced Early Holocene ice thinning, synchronous with the timing of MWP-1B. High uncertainty with regard to the timing of grounding-line retreat exists, particularly in the western Ross Sea, where the simulations exhibit considerable differences related to climate forcing selection. Establishment of age controls in this region should be prioritized. The ocean forcing is the main driver of grounding-line retreat following the formation of the RIS and the development of an ocean cavity. In the Holocene, the grounding line along Siple Coast shows strong sensitivity to small changes in subshelf melt rates, which has implications for future sea-level rise.


Regional ice sheet model setup and initialization

Before the climate forcing experiments, we performed a full-continent spin-up run following the protocol of Martin et al. (20). This first phase is a 100-year smoothing run, in which the model evolves by internal deformation only to remove steep gradients and any inconsistencies in the input ice thickness data. This was followed by a 200-ka thermal evolution run, in which the geometry remains fixed and the three-dimensional temperature field evolves to equilibrium according to the conservation of energy equation. Last, we performed a 100-ka evolutionary run with constant climate forcing. This length of spin-up is necessary because the thermal and dynamic evolution require at least a full glacial cycle due to the ice sheet thickness and slow velocities in the ice interior.

We then constructed a regional drainage basin model of the Ross Embayment based on the ice sheet thickness and topography. The bathymetry used in the model was based on the Bedmap2 dataset (34). We used the grounding-line scheme developed by Feldmann et al. (50), which uses a subgrid interpolation of basal driving stress to improve the accuracy of grounding-line movements. The regional model was tuned according to PD conditions (i.e., grounding-line position, ice shelf extent, surface ice velocity, and ice thickness of the entire domain) by adjusting model parameters, namely, the calving thickness threshold and rate in the calving parameterization scheme, enhancement factors of the shallow ice and shallow shelf approximation equations, the effective pressure on the till, and the till friction angle.

The regional model was run at 10-km resolution from 115 to 0 ka (total spin-up of 215 ka: full continent, 100 ka; regional, 115 ka). For the regional component of the model spin-up, we applied a paleoclimate forcing using the EDC temperature record (23) with paleoprecipitation determined on the basis of the Clausius-Clapeyron relationship to temperature as anomalies to a modern climatology (precipitation-temperature scaling of 7.0%/°C). For the first 100 ka of the regional simulation, the ice shelf extent was fixed, and for the remainder of the run (15 to 0 ka), the calving scheme was incorporated. We used the eigen calving parameterization developed by Levermann et al. (51), in which the calving rates are proportional to the product of the principle components of the strain rates, derived from the shallow shelf approximation velocities. The resulting 10-km regional PD configuration served as the initial condition for all experimental runs.

Subshelf melt rate/back pressure forcing sensitivity experiments

In the deglacial climate forcing experiments, we applied ocean forcing in the form of subshelf mass flux (melt rate) offsets. In PISM, the subshelf mass flux was used as a source in the mass continuity equation (20)tH=MSQ

In the above equation, t is time, H is ice thickness, M is surface mass balance, S is basal mass balance (S > 0 is melting), and Q is the horizontal ice flux. The offset is a positive or negative flux applied to the basal mass balance of the ice shelf (ice loss or gain, respectively) and is measured in meters per year. The anomaly is applied uniformly, as in Golledge et al. (10).

In addition, we applied back pressure (calving resistance) offsets, intended to mimic the resistive effects of iceberg mélange or dense sea ice that are not specifically modeled. Sea ice growth stiffens the mélange matrix, binding iceberg clasts together, thereby exerting a resistance against calving (28). Therefore, we assume here that thicker and more extensive sea ice associated with cooler climate states, such as the LGM, will result in reduced calving. Southern Ocean proxy records indicate more extensive and higher concentration at the LGM (52), and deglacial climate simulations exhibit strong negative correlations between sea ice concentration and thickness and SAT, with similar temporal evolution observed in the Ross Sea (18).

The back pressure offset is applied as a scalar fraction (λ) that is included in the stress boundary condition of the ice shelf (PISM source code). It is assumed that the back pressure (psea ice) does not exceed the pressure of the ice column at the calving front (picepocean); thus, the offset is in relation to the default value of 0 and is less than 1psea ice=λ(picepocean)bh(pice(pocean+psea ice))dz=(1λ)bh(picepocean)dz

In the above equations, λ represents the back pressure fraction (from 0 to 1), and b and h are the ice base and top surface elevations, respectively. With this forcing applied, the stress condition written in terms of the velocity components becomes modified as follows2νH(2ux+uy)nx+2νH(uy+vx)ny=(1λ)bh(picepocean)dznx,2νH(uy+vx)nx+2νH(2vy+ux)ny=(1λ)bh(picepocean)dznywhere ν is the vertically averaged ice viscosity, ux, y and vx, y are velocity components, and n represents the horizontal normal vector pointing oceanward from the ice boundary.

The basal melt rate forcing was applied based on the Southern Ocean δ18O stack (SOT) (25), and the back pressure forcing was applied based on the EDC temperature reconstruction (23). To respectively convert the δ18O and temperature reconstructions into basal melt rate and back pressure forcings, we used scaling relationships between two end-member states, i.e., LGM and PD, following the protocol of Golledge et al. (10). The scaling relationships for each forcing were determined through a suite of sensitivity experiments (table S1and figs. S4 and S5). These simulations initiated from the onset of the last interglaciation (131 ka) and were run to 0 ka and were performed at 20-km resolution. In addition to the above forcings, we applied a paleoclimate forcing using the EDC temperature record with paleoprecipitation determined on the basis of the Clausius-Clapeyron relationship to temperature as anomalies to a modern climatology (7%/°C) and sea-level forcing on the basis of statistical analysis of the sea-level proxy record (3033). Although LGM constraints were limited, these experiments were assessed in comparison to LGM grounding-line position and ice thickness estimations (6, 7) and PD conditions (i.e., grounding and calving line position, ice thickness, and ice surface velocity). The subshelf melt rate-back pressure forcing combination applied in simulations that best fit LGM and PD grounding-line position and ice thickness, and PD grounding and calving line position, ice thickness, and ice surface velocity, was selected for the deglacial experiments (figs. S4 and S5).

Climate forcing details

In our deglacial climate-forcing experiments, we used the temperature reconstructions of two Antarctic ice cores: the EDC and WDC ice cores (23, 24). Temperature reconstructions of ice cores are commonly calculated as anomalies to modern temperatures based on deuterium (δD), deuterium excess (d), and oxygen isotope measurements (δ18O) using site-specific temperature-dependence relationships that assume that the modern-day calibration holds over the entire record (23). However, surface temperatures can also be more accurately determined with additional borehole temperature and nitrogen isotope data, as was the case for WDC (24). In addition to the Antarctic ice core records, we applied subshelf melt rates that were based on a global δ18O stack (GOT) and a Southern Ocean δ18O (SOT) record using the process described above (25, 26). The GOT reconstruction was based on a compilation of 57 deep ocean records (25). The SOT was derived from a marine sediment core from offshore New Zealand (26).

The remaining climate forcings were derived from two transient climate simulations of the last deglaciation. TraCE-21ka is a transient climate simulation of the last 21,000 years using the Community Climate System Model version 3, a coupled general circulation model with atmosphere, ocean, land surface, and sea ice components (22). The experiment included evolving orbital forcing according to Milankovitch theory, greenhouse gas concentrations (CO2, CH4, and N2O), ice sheets and paleogeography, and imposed meltwater fluxes into the North Atlantic and Southern Ocean. The LOVECLIM deglacial experiment was performed with the intermediate complexity LOVECLIM model version 1.1, which includes atmosphere, ocean, sea ice–ocean, ocean carbon cycle, and terrestrial vegetation components, for the period of 18 to 3 ka (21). Similar to TraCE 21ka, the transient experiment involved time-varying solar insolation, Northern Hemisphere ice sheet topography, and atmospheric CO2 determined from the EDC ice core, with CH4 and N2O fixed at LGM levels. Freshwater pulses were applied to the North Atlantic and Southern Ocean based on 231Pa/230Th data of the North Atlantic. Since the climate model runs do not cover the entire time period of the deglacial ice sheet simulations (i.e., 35 to 0 ka), we used averages of the proxy records for the preceding years of the model runs.

As the climate forcings were applied uniformly as scalar anomalies to a modern climatology based on RACMO2.3 output (27), we calculated regional average anomalies over the domain to obtain the atmosphere forcings of the climate model simulations. Since Antarctic ice sheet topography was not updated in the LOVECLIM deglacial experiment, the continental interior showed relatively lower LGM temperature anomalies than TraCE-21ka. Therefore, we considered two versions of the LOVECLIM atmosphere forcing: one covering the entire domain, which includes the likely biased continental interior, and one constrained to the continental margin. The subshelf melt rates were based on regional averages of ocean temperatures at 450-m depth in the Ross Sea (70°S–62°S, 168°E–160°W). We also considered a “cool” ocean case using the TraCE-21ka SST as the subshelf melt rate forcing. All subshelf melt rates were scaled to the same LGM anomaly of 5.0 m year−1, which is based on the sensitivity experiments described above.

Deglacial climate forcing experiments

Starting from the regional PD configuration described above, we performed an initial climate forcing simulation from 131 ka (last interglacial) to 35 ka (early glacial) at 10-km resolution. This simulation was forced by EDC-derived temperature, precipitation, and back pressure forcing; subshelf melt rate forcing derived from the SOT reconstruction; and sea-level forcing based on statistical analysis of multiple sea-level records (3033). From this 35-ka configuration, we then performed 36 climate forcing experiments using the atmosphere and ocean forcing combinations listed in Table 1. These experiments were run from 35 to 0 ka at 10-km resolution. All model parameter values applied in these experiments are listed in table S2, and the ensemble average deglacial ice sheet evolution is shown in fig. S6.

In addition to this main ensemble set, we performed a sensitivity experiment that is based on the scenario proposed by Kingslake et al. (9). In this experiment, which was also run from 35 to 0 ka at 10-km resolution, we used all of the same model parameters in the main ensemble set, with one exception. The mantle viscosity term in this simulation was 5 × 1020 Pa s, whereas in the main ensemble set, the term was 1 × 1019 Pa s. Additionally, we used a subshelf melt rate forcing based on the WDC temperature record, which has a steeper increase in the early deglacial period. The higher mantle viscosity and more aggressive ocean forcing were applied to reproduce the early retreat and isostatic readvance in the Holocene simulated by Kingslake et al. (9), and we showed similar timing of retreat and minimum grounding-line extent (fig. S2), i.e., 9.8 ka in our simulation versus 9.7 ka by Kingslake et al. (9).


Supplementary material for this article is available at

Fig. S1. Sea-level forcing effect on ice sheet–ice shelf evolution.

Fig. S2. High mantle viscosity/ocean forcing simulation.

Fig. S3. Southern Ocean freshwater forcing effect.

Fig. S4. LGM (20 ka) results of the subshelf melt rate (ssmr)–back pressure (bp) experiments.

Fig. S5. PD (0 ka) results of the ssmr-bp experiments.

Fig. S6. Moderate ocean/atmosphere simulation.

Table S1. Subshelf mass flux–sea ice back pressure sensitivity experimental setup.

Table S2. Model parameters used in the deglacial climate-forcing experiments.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge A. Aschwanden and C. Khroulev for constructive advice regarding the PISM, and the teams behind the TraCE-21ka and LOVECLIM DGns experiments for producing and sharing model output, publicly available via the NCAR Climate Data Gateway and the Asia-Pacific Data Research Center, respectively. We thank the Antarctic ice core and marine proxy communities for the use of their data. We also thank M. Kelly and four anonymous reviewers for their insightful comments on the manuscript. Funding: Funding for this project was provided by the Royal Society Te Aparangi Marsden Fund through Victoria University of Wellington (15-VUW-131); the New Zealand Ministry of Business, Innovation, and Employment Grant through GNS Science (540GCT32); and the New Zealand Antarctic Research Institute (NZARI2014-11). The development of PISM was supported by NASA grant NNX17AG65G and NSF grants PLR-1603799 and PLR-1644277. D.P.L. acknowledges support from the Antarctica New Zealand Doctoral Scholarship program. N.R.G. acknowledges support from the Royal Society Te Aparangi under contract VUW1501. R.S.J. was supported by a Junior Research Fellowship COFUNDed between Durham University and the European Union under grant agreement number 609412. R.M. acknowledges support from the Royal Society Te Aparangi Rutherford Discovery Fellowship (RDF-13-VUW-003). N.A.N.B. acknowledges support from the Royal Society Te Aparangi Rutherford Discovery Fellowship (RDF-VUW-1103). Author contributions: D.P.L., N.R.G., and N.A.N.B. devised the regional ice sheet model experiments. D.P.L. performed the model experiments and analysis. R.S.J. provided the cosmogenic nuclide surface-exposure data. All authors contributed to the development of ideas and writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article