The evolution of the Antarctic ice sheet during the mid-Pliocene warm period (MPWP) remains uncertain and has important implications for our understanding of ice sheet response to modern global warming. The extent to which marine-based sectors of the East Antarctic Ice Sheet (EAIS) retreated during the MPWP is particularly contentious, with geological observations and geochemical analyses being cited to argue for either a relatively minor or a significant ice sheet retreat in response to mid-Pliocene warming. The stability of marine-based ice sheets is intimately linked to bedrock elevation at their grounding lines, and previous ice sheet modeling assumed that Antarctic bedrock elevation during the MPWP was the same as today with the exception of a correction for the crustal response to ice loading. However, various processes may have perturbed bedrock elevation over the past 3 m.y., most notably vertical deflections of the crust driven by mantle convective flow, or dynamic topography. Here we present simulations of mantle convective flow that are consistent with a wide range of present-day observables and use them to predict changes in dynamic topography and reconstruct bedrock elevations during the MPWP. We incorporate these elevations into a simulation of the Antarctic ice sheet during the MPWP and find that the correction for dynamic topography change has a significant effect on the stability of the EAIS within the marine-based Wilkes Basin, with the ice margin in that sector retreating considerably further inland (200–560 km) relative to simulations that do not include this correction for bedrock elevation.
The mid-Pliocene warm period (MPWP), from 3.264 Ma to 3.025 Ma (Dowsett et al., 2010), was characterized by atmospheric CO2 levels similar to that of today (400 ppm; Pagani et al., 2010) and by global mean temperatures that were on average elevated by ∼2–3 °C relative to today (Dowsett et al., 2009). Accordingly, the period is considered an important case study for investigating polar ice sheet response in a modestly warmer world. There is general agreement that the Greenland Ice Sheet and the West Antarctic Ice Sheet experienced widespread deglaciation during the MPWP (Masson-Delmotte et al., 2013; Naish et al., 2009; Dutton et al., 2015); however, the extent to which the East Antarctic Ice Sheet (EAIS) retreated during the same period remains contentious, with arguments ranging from a cold-based and stable EAIS (e.g., Sugden et al., 1993) to a warmer and highly variable EAIS (e.g., Webb et al., 1984). Sediments recovered from cores in the Ross Sea (Naish et al., 2009) and across the Antarctic margin (Cook et al., 2013; Patterson et al., 2014) point to a variable, obliquity-paced EAIS in the vicinity of the Transantarctic Mountains and the Wilkes Basin (Fig. 1).
Simulations of Antarctic ice sheet evolution under MPWP conditions are generally characterized by a much smaller West Antarctic Ice Sheet with little or no marine-based ice, but a stable EAIS with minor melting limited to the outer margins of some marine-based sectors (de Boer et al., 2015). Calculations by Pollard and DeConto (2009) suggest a net ice mass difference relative to the present day equivalent to ∼7 m of global mean sea-level rise. Combining this scenario with modeling of the Greenland Ice Sheet during MPWP yields a total global mean sea-level rise of ∼14 m (Masson-Delmotte et al., 2013). Pollard et al. (2015) have recently augmented their modeling to consider two additional ice sheet instability mechanisms, namely hydrofracturing due to surface melt drainage into crevasses and gravitational failure of ice cliffs. The actual impact of these processes on Antarctic ice sheet stability during the MPWP is uncertain, but acting in tandem, they could have significantly enhanced EAIS retreat, yielding a total Antarctic contribution of as much as ∼17 m.
Ultimately, the stability of a grounded, marine-based ice sheet is largely governed by a balance between mass added to the ice sheet through precipitation and mass lost by outflow at the grounding line. The amount of outflow strongly depends on the ice thickness at the grounding line (Schoof, 2007), and therefore any changes in relative sea level, or equivalently bedrock elevation, at the grounding line will perturb the mass balance. Over the multi-million-year time span from the mid-Pliocene to present, evolving viscous stresses and buoyancy associated with mantle convection (e.g., Mitrovica et al., 1989; Gurnis, 1990) could have driven large changes in bedrock elevation. These long-term mantle-driven changes in surface topography, termed dynamic topography, could either stabilize or destabilize the ice sheet depending on the sign of the crustal deflection. There are reasons to believe that mantle convective flow may have played a role in the evolution of Antarctic topography, including beneath some parts of the EAIS (Faccenna et al., 2008), and influenced ice mass loss across local grounding lines, however these effects have never been incorporated in ice sheet reconstructions.
In this paper, we reconstruct Antarctic bedrock elevation during the mid-Pliocene using viscous flow modeling of dynamic topography constrained by a broad suite of seismic, mineral physics, and geodynamic data sets (Forte et al., 2015). We then explore to what extent this change in bedrock elevation could have impacted EAIS stability during the MPWP.
MANTLE FLOW BENEATH THE ROSS ICE SHELF AND WILKES BASIN
We use the finite element code ASPECT (aspect.dealii.org; Kronbichler et al., 2012) to solve the equations governing conservation of momentum, mass, and energy for compressible mantle flow driven by temperature-induced buoyancy (see the GSA Data Repository1 for additional modeling details). Our simulations require, as input, the present-day density and viscosity structure within the mantle. We use the density model TX2008, which was derived by a joint inversion of a large database of shear-wave travel times and geodynamic data, with a scaling between seismic wave speeds and density based on constraints from mineral physics (Simmons et al., 2009). The inversion was also based on the radially varying mantle viscosity profile V2 (Fig. DR1 in the Data Repository), which itself was derived from a joint inversion of the geodynamic data set and a suite of observables related to glacial isostatic adjustment (Forte et al., 2010). The mantle flow computed on the basis of the density and viscosity field yields normal stresses at the surface of the model that are used to compute dynamic topography.
We additionally incorporate a temperature-dependent viscosity field in the flow simulations; our choice of parameters yields lateral variations in viscosity of approximately ±2 orders of magnitude at each depth (see the Data Repository for details). Introducing this more realistic viscosity structure has only a minor effect on the instantaneous solution for present-day mantle flow (and thus we preserve the fits discussed above), but it has a significant impact on the predictions of the rate of change of dynamic topography (Gurnis et al., 2000; Moucha et al., 2007).
Figures 2A and 2B show lateral temperature variations (relative to the average at each depth) as well as mantle flow computed in two vertical cross-sections that are oriented to pass through the Wilkes Basin. The temperature fields are dominated by a hot (buoyant) anomaly that extends from the top of the lower mantle into the shallow upper mantle; the anomaly is located beneath the Ross Ice Shelf and the Transantarctic Mountains and extends to the northwest around the thick East Antarctic craton, ultimately leading to shallow corner flow and upwelling at the edge of the craton. This same anomaly is a robust feature in a suite of different seismic tomography models (see Fig. DR3). The computed flow beneath the Ross Ice Shelf is characterized by two counter-rotating convection cells and associated divergent horizontal flow in the shallow mantle, consistent with arguments that an active mantle upwelling is the driving force for Cenozoic extension in the West Antarctic rift system (Faccenna et al., 2008) as well as volcanism in several nearby locations (Faccenna et al., 2008; Gupta et al., 2009). Test calculations (not shown) indicate that density variations below the lithosphere drive uplift in the area of the Transantarctic Mountains, however the shallow corner flow north of the Wilkes Basin is significantly reduced if density variations in the upper 200 km are removed from the calculation.
DYNAMIC TOPOGRAPHY AND UPLIFT
Changes in local dynamic topography (i.e., changes measured at a site fixed to a tectonic plate) are computed by combining the change in dynamic topography in the flow model domain with the motion of a plate through this topography field. In our calculation of dynamic topography, we adopt the Euler pole derived by Sella et al. (2002) to describe the motion of the Antarctic plate relative to the underlying mantle.
The computed change in local dynamic topography over the past 3 m.y. is shown in Figure 2C. The large, thermally buoyant anomaly evident in Figures 2A and 2B drives uplift of the Transantarctic Mountains that extends into both the western and northern margins of the Wilkes Basin. The western margin has contributions from both the geometry of the upwelling as well as plate motion over this upwelling. The predicted change in dynamic topography in the Dominion Range is ∼50 m, which is consistent with the bound on uplift inferred from geomorphological analyses and surface exposure dating (Ackert and Kurz, 2004). The computed uplift in the Dry Valleys and Adare Peninsula is significantly larger, on the order of 250 m. This value is consistent with Pliocene uplift estimates based on submarine and subaerially erupted volcanic rocks (Wilch et al., 1993; Mortimer et al., 2007).
We have performed a series of analyses to test the sensitivity of our predictions of dynamic topography to a suite of input variables (see Fig. DR4). These analyses indicate that the computed uplift is a robust feature of all simulations, although the magnitude of the change in dynamic topography is sensitive to the adopted viscosity structure. We note that our reconstruction of the bedrock elevation does not include changes associated with sediment transport (Hill et al., 2007; Wilson et al., 2012), which are likely small on the time scale we are considering.
ICE SHEET STABILITY
We use an established Antarctic ice sheet model (Pollard and DeConto, 2009; DeConto et al., 2012) to investigate whether reconstructions of EAIS stability during the MPWP are sensitive to the incorporation of dynamic topography changes in bedrock elevation. This model tracks ice thickness and temperature distributions in Antarctica that result from gravitationally driven ice sheet deformation as well as mass addition and removal due to precipitation, basal melt and runoff, oceanic melt, and calving of floating ice (for details, see Pollard and DeConto, 2012). Furthermore, the model includes a crustal rebound term associated with ice unloading that assumes an elastic lithosphere and a viscous relaxation time of 3000 yr, and it applies the parameterization for ice flux across the grounding line derived by Schoof (2007). We force the model with Pliocene conditions (see Pollard et al., 2015, and references therein), including the adoption of a warm atmospheric climate from a slightly modified version of the RegCM3 regional climate model (users.ictp.it/RegCNET/model.html), atmospheric CO2 concentration of 400 ppmv, ocean warming of 2 °C, and fixed, orbitally driven changes in summer insolation (DeConto et al., 2012). The ice model is first spun up to modern conditions, which is followed by an instantaneous change in forcing to a warmer Pliocene climate. The model is then run for an additional 5000 yr to allow for an equilibrium state to establish.
We ran two simulations of Antarctic ice cover during the MPWP, one that does not include dynamic topography changes in reconstructing Pliocene bedrock elevation (Fig. 3A) and one that does (Fig. 3B). The western and northern margins of the Wilkes Basin have been uplifting due to changes in dynamic topography for at least the past 3 m.y. (Fig. 2C), in contrast to other marine-based sectors of the EAIS, which have been subject to much smaller dynamic topography changes (e.g., Aurora Basin). Thus, the elevation correction for changes in dynamic topography lowers the estimated mid-Pliocene bedrock elevation on the margins of the Wilkes Basin and makes the local, marine-based ice sheet more susceptible to retreat. In particular, grounding line retreat in the simulation that includes a dynamic topography correction extends ∼400 km further inland into the Wilkes Basin than the simulation based on present-day bedrock elevation. The integrated difference in ice melting in the two simulations is equivalent to ∼2 m of global mean sea-level change.
We ran ice sheet simulations for a suite of mantle flow calculations that were performed to test the sensitivity of our predictions of dynamic topography change to model inputs (see Fig. DR5). In this set of simulations, the correction of the mid-Pliocene bedrock elevation for dynamic topography leads to grounding line retreat of 200–560 km at the margin of the Wilkes Basin and an ice volume change equivalent to a global mean sea-level difference of 1.1–2.8 m relative to a simulation with no change in dynamic topography (or 5.7–7.4 m relative to present-day Antarctic ice cover).
Pollard et al. (2015) have argued that the combined effects of ice cliff failure and hydrofracturing may have enhanced the retreat of the EAIS, including within both the Wilkes Basin and adjacent Aurora Basin, during the MPWP. We performed an additional series of simulations that include these mechanisms (Fig. DR6). These calculations indicate that changes in dynamic topography will enhance the retreat of grounded ice within the Wilkes Basin even when the cliff failure and hydrofracturing processes are included. We note that if one were to model each of these three mechanisms in isolation, only the change in bedrock elevation due to dynamic topography would have been sufficient to drive significant ice sheet retreat in the Wilkes Basin during the MPWP.
In contrast to these destabilizing mechanisms, gravitationally self-consistent sea-level changes at the margin of an evolving EAIS will act to stabilize the ice sheet (Gomez et al., 2013). This effect will likely be small in the ice sheet modeling described here, in which the system is assumed to be close to isostatic equilibrium, but it may play a role in determining ice sheet extent during time-evolving Pliocene ice age cycles.
We have demonstrated that incorporating bedrock elevation changes due to dynamic topography is important when trying to reconstruct the EAIS in the vicinity of the Wilkes Basin during the MPWP. Mantle flow modeling of dynamic topography predicts that bedrock elevation was ∼100–200 m lower on the western and northern margins of the Wilkes Basin during the mid-Pliocene, and including this change in ice sheet model simulations introduces a destabilizing influence on the local ice margin causing additional 200–560 km of grounding line retreat into the Wilkes Basin. The amplitude of dynamic topography change in other marine-based sectors of the EAIS (e.g., Aurora Basin, Prydz Bay; Fig. 1) is relatively small in our simulations and has therefore little impact on the ice sheet stability in these sectors during the MPWP.
Our prediction is consistent with geochemical analysis of Integrated Ocean Drilling Program core U1361, which has been used to argue for a retreat of the ice margin of several hundred kilometers in the Wilkes Basin during the Pliocene (Cook et al., 2013). However, this argument assumes that the enhanced erosion inferred from the geochemical analysis was sourced from the ice margin. Direct geologic evidence of ice sheet retreat in the Wilkes Basin is lacking, making it an important target for future field research. Finally, the impact of dynamic topography changes on ice sheet stability is not limited to the Pliocene and has likely affected the evolution of the Antarctic ice sheet since its inception.
We thank Tim Naish, Daniel Hill, and one anonymous reviewer for their constructive comments on this manuscript, as well as Robert Ackert, Jr., for a number of discussions. We thank Wolfgang Bangerth and Timo Heister for making the ASPECT code available to the community, and the Computational Infrastructure for Geodynamics (geodynamics.org) for supporting Austermann’s participation at two ASPECT hackathons. Support for this research was provided by National Science Foundation Division of Ocean Sciences grant OCE-0825293 “PLIOMAX” and Harvard University. Forte acknowledges support provided by the Natural Sciences and Engineering Research Council of Canada.
↵1GSA Data Repository item 2015310, additional information on the solid Earth modeling (the file 2015310-Model-output.zip contains ASCII files of the dynamic topography predictions along with a readme file that explains the format and a MATLAB script that can be used to plot the model output), and figures exploring uncertainty, is available online at www.geosociety.org/pubs/ft2015.htm, or on request from or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.
- Received 18 May 2015.
- Revision received 12 August 2015.
- Accepted 13 August 2015.
- © 2015 Geological Society of America