Dynamics of the last glacial maximum Antarctic ice-sheet and its response to ocean forcing

Geological data provide essential constraints to our ice-sheet modeling experiments. Pan-Antarctic cosmogenic surface exposure ages from bedrock samples and glacial erratics indicate that ice thicknesses close to the present coast were several hundred meters greater during the LGM than at present, and tapered inland (SI Appendix, Table S1), in agreement with inferences from ice cores (18). This terrestrial record constrains the vertical dimensions of the former ice sheet, whereas marine geologic interpretations from swath bathymetric and seismic surveys allow the lateral extent of the ice sheet to be reconstructed (19, 20).

Through systematic iteration of model parameters, we achieve an optimum simulation in which both the surface elevation and lateral extent of the ice sheet accords with the majority of geological data (Fig. 1A and SI Appendix, Table S1). Across much of East Antarctica, surface elevations of the simulated ice sheet are comparable to present values, although some inland ice divides are modeled to have been lower (Fig. 1B). Considerable thickening occurs around the coast to the extent that the present-day Filchner–Ronne, Ross, and Amery ice shelves are replaced by grounded ice that extends across much of the continental shelf (Fig. 1A). Modeled ice thicknesses over Ross Island are almost identical to empirically inferred values (21), but our simulation does not reproduce the ice-free McMurdo Dry Valleys, and so overestimates ice-surface elevations in this area. Surface exposure ages in the southern Ross Sea (22) and eastern Weddell Sea (23) embayments suggest that our modeled ice sheet is also too thick in these areas. However, the limited thickening implied by these empirical data, coupled with the greatly advanced grounding-line position interpreted from marine geological data, can only be reconciled with a surface slope of the LGM grounded ice sheet that is similar to that of the present ice shelf, requiring extremely low basal shear stress (< 15 kPa). This disagreement with the observations thus requires further investigation.

Isostatic depression of the bed is greatest in the Weddell Sea sector, with lesser amounts of loading in the Ross Sea and Amery embayments (Fig. 1 C–E)—consistent with global positioning system (GPS) measurements of present-day bedrock uplift in Antarctica (SI Appendix, Fig. S7).

Throughout the domain, modeled flowlines agree closely, but not perfectly, with LGM flow directions inferred from orientations of mega-scale glacial lineations (20, 24, 25) and from marine sediment core mineral provenance studies (26) (Fig. 2 and SI Appendix, Fig. S1). Highest modeled sliding velocities and maximum discharge rates occur in the Thiel/Crary Trough (eastern Weddell Sea), and in the eastern and central Amundsen Sea where Pine Island and Thwaites glaciers coalesce. Modeled discharge through the Amery Trough is also high, and sliding velocities here exceed 500 ma^-1 in the outer trough. Although sliding velocity and discharge rates throughout the Ross Sea are lower and less well-partitioned than any of these three other outlets, fast flow influences a much larger sector of this embayment (Fig. 2). The geometry and pattern of flow of our steady-state LGM Antarctic ice sheet therefore fit well with both terrestrial and marine geological constraints, and accounts for a net increase in grounded ice volume of 2.702 × 10⁶ km³ (6.67 m sea level equivalent). The physics and fine resolution of our model of the LGM Antarctic ice sheet permit us to resolve far greater spatial variability in the balance between dominant mechanisms of flow than previously simulated (27, 28). Where ice flow is primarily driven by viscous deformation in response to gravitational driving stress, basal velocities are close to zero, but where longitudinal coupling dominates the force balance of the modeled ice sheet, basal velocities are much higher. According to our model, basal sliding is the sole contributor to glacier motion in large sectors of West Antarctica and in coastal areas of East Antarctica (Fig. 3A). In the eastern Weddell Sea sector, ice draining the East Antarctic ice sheet is organized into discrete catchments that each nourish narrow, fast-flowing, conduits that anastamose around islands of slower ice (Fig. 2). A transect perpendicular to flow spanning 84 to 83 °S in this quadrant of the domain reveals that these fast-flowing conduits exhibit basal velocities that increase abruptly by up to three orders of magnitude across the creep-sliding transition (Fig. 3B). In contrast, surface velocity changes between fast and slow zones are typically half.

To establish likely ice-sheet response to changing oceanic boundary conditions, we carried out a suite of sensitivity experiments involving sea-level and ocean temperature perturbations. The results show that highly partitioned flow, with abrupt lateral boundaries separating sliding and nonsliding ice, leads to spatially variable responses (Fig. 4). The fast-flowing ice-sheet outlets in our model respond instantly to perturbations at oceanic margins and propagate changes inland very rapidly—a dynamical sensitivity that is due largely to the effects of longitudinal coupling (14). Significantly, although our simulations show glacier acceleration confined to the conduits where enhanced flow occurs, changes in ice thickness are witnessed across far more extensive areas, reflecting a substantial drawing down of the surrounding ice-sheet surface. The conduits thus become the foci of greatest mass loss from the ice sheet. Assuming spatially uniform increases in sea level and oceanic heat flux, we identify that the outlets most susceptible to oceanic changes during deglaciation would likely have been those in the Thiel/Crary Trough, the eastern and central Amundsen Sea, the central Ross Sea, and to a lesser extent the Amery Trough (Fig. 4). In our model, these sectors respond more sensitively to oceanic forcings than slower-flowing or deformation-dominated sectors of the ice sheet, by accelerating and drawing down ice from their entire catchment areas. Geochronological data lend some support to this notion—relatively early deglaciation (> 16 ka) is recorded close to dynamic outlets in these areas (8, 9, 19, 29, 30), whereas retreat occurred later (< 15 ka) in neighboring, but significantly less mobile areas, including marginal areas of the Ross Sea (6, 7, 31–34) (Fig. 2). Despite this encouraging agreement, however, we cannot entirely exclude the possibility that some coastal areas respond early simply because of their proximity to open ocean, or that some inland sites respond later because of their greater distance from the ice-sheet margin. A greater density of reliable retreat ages from around the continent may help test these hypotheses in the future.

Further experimentation may also help to determine whether the deglacial behavior of an LGM ice sheet that was out of equilibrium, continually adjusting to time-transgressive changes in ocean temperature, sea-level, and atmospheric conditions, would be the same as that of our steady-state ice sheet. If the LGM ice sheet were still expanding when deglaciation began, for example, it may have taken longer to respond to negative mass balance at the margin than an ice sheet that was already retreating, perhaps leading to delayed retreat in some areas. Nonetheless, in our steady-state experiments we find that the rate of grounding-line migration is controlled by the magnitude of the imposed forcings, whereas the regional pattern of ice margin retreat depends on the rate and volume of ice discharge from inland areas, the presence of pinning points, and the bathymetry of the embayments where ice loss takes place. Thus we infer that an indented margin characterized by calving bays that were fed by fast-flowing conduits, as predicted by our model in embayments such as the Ross Sea, is likely representative of the morphology of the retreating LGM ice-sheet margin, irrespective of whether deglaciation initiated from an equilibrium configuration or not.

By using an empirically constrained high-resolution ice-sheet model to simulate the Antarctic ice sheet at the LGM, and by forcing this model with oceanic perturbations, we conclude that spatial contrasts in Antarctic ice-sheet dynamics played a more important role in modulating southern hemisphere ice-sheet sensitivity to ocean forcings than previously realized, and may explain patterns of retreat not otherwise accounted for by regional enhancement of sea levels or grounding-line stabilization from sediment accumulation. That widespread ice-sheet thinning arises from localized acceleration of ice streams is particularly significant when considering recent observations of velocity increases in both Antarctic and Greenland ice sheets (35).

We use a three-dimensional, thermomechanical, continental ice-sheet model constrained by geological data that define lateral and vertical extents of the expanded Antarctic ice sheets around the time of the LGM. As in previous studies, we employ boundary distributions from modified BEDMAP topography (36), temperature and precipitation fields from gridded datasets (37, 38), and a spatially varying geothermal heat flux interpolation (39). Our model computes ice thickness and temperature changes, isostatic depression of topography, migration of grounding lines, and the growth of ice shelves. Interaction between modeled ice shelves and their surrounding ocean is accounted for using a mass balance determination based on heat flux across the ice-water boundary. We employ a stress balance that includes longitudinal (membrane) stresses (SI Appendix, Materials and Methods), impose boundary conditions (sea-level lowering, precipitation reduction, and atmospheric temperature perturbations) representative of the LGM, based on ice and marine sediment core isotopic deviations (40, 41) (SI Appendix, Fig. S5), and adjust model parameters affecting bed traction, ice rheology, and ice-shelf mass balance. In contrast to other studies (42), however, the geometry and dynamics of our modeled ice sheet are able to evolve naturally, because we do not prescribe grounding line or ice stream locations. Furthermore, we make use of parallel processing to implement our model at a uniquely high (5 km) resolution, achieving a 16–64 times increase in detail compared to other Antarctic simulations (27, 42, 43).

We are grateful to Ed Bueler, Constantine Khroulev, and Andy Aschwanden for help with the Parallel Ice Sheet Model, and to Tony Dale and Vladimir Mencl (University of Canterbury) for access to and assistance with the Bluefern Supercomputer. Tim Naish, Peter Barrett, Rob McKay and two anonymous reviewers are gratefully acknowledged for comments on previous versions of this manuscript. N.R.G. and A.N.M. acknowledge financial support from Victoria University Foundation Grant, "Antarctic Research Centre Climate and Ice-Sheet Modelling”. C.J.F. is supported by Australian Research Council Fellowships FL100100195 and FT120100004.