A new global ice sheet reconstruction for the past 80 000 years

The study addresses long-standing uncertainties in past global ice sheet evolution, particularly the Last Glacial Maximum (LGM) “missing ice” problem—an 8–28 m mismatch between far-field sea-level indicators and modelled global sea level from previous reconstructions. Prior reconstructions often tuned to far-field sea level and marine δ18O records, potentially introducing circularity and bias, especially before the LGM when direct sea-level data are sparse. The authors propose an ice sheet reconstruction built from glaciologically plausible physics and direct geological constraints (ice margins and flow indicators), avoiding δ18O and far-field sea-level tuning, to test whether: (1) the LGM lowstand can be matched without invoking extra ice volume; and (2) the assumed relationship between δ18O and global ice volume/sea level during MIS 3 is valid.

Previous GIA-based reconstructions fit constraints without strong ice-physics guidance or used scaled thicknesses and ensembles, often assuming far-field sea level as part of tuning. Marine δ18O has been used as a proxy for pre-LGM ice volume in the absence of direct sea-level data. Regional evidence from the Hudson Bay area suggests possible ice-free conditions during parts of MIS 3 and climatic conditions conducive to forest growth, implying a reduced Laurentide Ice Sheet; GIA studies of US east-coast sea levels also support reduced Laurentide extent. However, dating uncertainties near mid-MIS 3 challenge some interpretations, and proxy-based sea-level reconstructions often suggest much lower sea levels (-60 to -90 m) than some GIA and geological indicators allow. This background motivates a reconstruction independent of δ18O and far-field sea level to reassess global ice volume history.

Glacial isostatic adjustment (GIA) and sea-level calculations were performed with SELEN (solving the sea-level equation including shoreline migration, grounding-line adjustments, and rotational feedback) using a 1D, spherically symmetric Earth model: 120 km elastic lithosphere, upper mantle viscosity 4×10^20 Pa s, lower mantle viscosity 4×10^22 Pa s. Calculations used spherical harmonic degree 256 with three iterations to update the ocean function. Sea level and Earth deformation were output on a 1° grid and interpolated (Delaunay triangulation) to the 5 km ice-model grid. The reconstruction spans 80 ka BP to present with 2,500-year time steps (constant time step limitation in SELEN); a test with linearly interpolated 500-year steps improved near-field fits but yielded similar far-field responses. Ice sheets were reconstructed using ICESHEET (perfectly plastic, quasi-equilibrium assumption) requiring margin location, basal shear stress, and basal topography (RTopo-2). Workflow: (1) initial ice sheets on modern topography using shear stress informed by topography and surficial geology (tuned to reproduce modern Greenland/Antarctica, and prior studies for North America/Eurasia); (2) compute GIA and derive paleotopography; (3) recompute ice sheets on deformed paleotopography; (4) recompute GIA and sea level at indicator sites; (5) tune basal shear stress (including temporal reduction during deglaciation) to reduce misfit to sea-level indicators and fit broad-scale ice-flow directions; (6) iterate steps 4–5; (7) for pre-LGM periods with sparse sea-level data, match inferred ice-flow direction changes and geological constraints, adjusting margins and shear stress; (8) perform a final iteration without changing shear stresses to align ice load, deformation, and sea level. Margin reconstructions: LGM–present followed existing datasets; Eurasian pre-20 ka margins from published reconstructions indicating limited MIS 3 extent. North American and Antarctic pre-LGM margins were newly developed. The chronology of Laurentide advance/retreat events was aligned with Heinrich Events (e.g., MIS 4 maximum at 60 ka, HE6), acknowledging several-kyr uncertainties. Two MIS 3 scenarios were developed: minimal (temporary deglaciation of Hudson Bay around HE5, ~45 ka, with re-glaciation by 40 ka) and maximal (maintained broader coverage through MIS 3 with fluctuations timed to Heinrich events). Cordilleran ice had major MIS 4 and MIS 2 advances with limited MIS 3 extent; Greenland MIS 4 similar to MIS 2 with substantial MIS 3 retreats in some sectors; Innuitian pre-LGM extent limited to modern ice caps; Antarctic margins advanced gradually in West Antarctica to shelf-edge maximum by ~30 ka while East Antarctica stayed near present. Basal shear stress tuning principles: higher in mountainous terrain, lower over continuous sediment cover, and reduced in major ice-stream corridors (e.g., low in Hudson Strait). Deglacial shear stress reduced to reflect basal warming/meltwater and interactions with lakes/ocean. Lower mantle viscosity sensitivity tests showed higher viscosities (~10^23 Pa s) improved fits in formerly glaciated regions; 4×10^22 Pa s was adopted as a compromise enabling larger Laurentide core loads while fitting near-field Holocene sea level. Evaluation used global sea-level indicator databases, comparing site-by-site modeled relative sea level across multiple Earth models and time-step variants.

• PaleoMIST 1.0, a global ice sheet reconstruction for the past 80 ka built independently of far-field sea-level and δ18O proxies, matches far-field LGM sea-level constraints at multiple sites (Barbados, Sunda Shelf, Bonaparte Gulf, Great Barrier Reef) without requiring extra ice volume.
• Maximum LGM ice-volume equivalent sea-level fall is −116 m, with total ice volume 42.2×10^6 km^3; Antarctic contribution is ~10 m SLE. This is substantially less than some earlier estimates (−130 to −135 m; 53×10^6 km^3), removing the need for “missing ice.”
• The difference in global mean sea level between the minimal and maximal MIS 3 scenarios peaks at ~14 m SLE.
• MIS 4 global ice volume attains ~70% of the LGM value; Eurasian ice volume comparable to its LGM due to additional Kara Sea ice, while western Laurentide does not reach the Cordillera.
• Laurentide Ice Sheet geometry dominates global sea-level variability. Best-fitting Earth rheology for maximizing volume while fitting near-field sea level uses a high lower mantle viscosity (4×10^22 Pa s).
• During MIS 3, modeled sea levels are substantially higher (less negative) than δ18O-based reconstructions suggest; in the minimal scenario, temporary Hudson Bay deglaciation produces a sharp sea-level rise up to ~25 m. Even the maximal scenario generally remains above ~−50 m between 50–35 ka BP, inconsistent with proxy reconstructions of −60 to −90 m.
• Spatial variability at the LGM implies some far-field sites may record sea levels several meters below the global mean due to GIA and Earth rheology, emphasizing the need for site-specific interpretation rather than a single global target.

By constructing an ice sheet history constrained by ice physics, geological margins, flow indicators, and near-field sea-level data—while excluding far-field sea-level and δ18O tuning—the study demonstrates that the LGM lowstand can be reproduced with smaller total ice volume, eliminating the classic “missing ice” problem. The discrepancy with earlier larger-volume budgets likely stems from differing initial Antarctic loads and the compensating need for larger Northern Hemisphere ice in those studies, compounded by gravitational effects. The results further show that MIS 3 δ18O-based sea-level reconstructions are incompatible with geological and GIA constraints, indicating that the relationship between marine δ18O and global sea level/ice volume is more complex (e.g., influenced by temperature effects in tropical records or highstand durations too brief for proxy detection). Modeled spatial patterns highlight that local relative sea level at far-field sites can deviate from the global mean due to Earth rheology and loading geometry, suggesting reconstruction efforts should treat each site as an independent observation within a global network rather than target a single global lowstand. Overall, the findings address the research questions: they resolve the LGM budget without extra ice and call into question δ18O–sea level scaling during MIS 3.

The paper introduces PaleoMIST 1.0, a global reconstruction of ice sheets over the past 80 kyr that, without using far-field sea-level or δ18O proxies, reproduces far-field LGM sea levels with a smaller total ice volume, thereby resolving the LGM “missing ice” problem. It identifies a new inconsistency for MIS 3: proxy-based δ18O sea-level reconstructions conflict with geological and GIA-constrained ice configurations. The authors recommend future work to refine Earth rheology assumptions, ice margin chronologies, and to develop improved δ18O–sea level relationships (e.g., via water isotope modeling), along with higher temporal resolution GIA–ice coupling and expanded, site-specific sea-level indicator datasets to reduce non-uniqueness.

• Temporal resolution: 2,500-year time steps (due to SELEN’s constant time-step requirement) limit fidelity in deglacial near-field regions; although a 500-year interpolated test improves near-field fits, the main reconstruction remains coarse in time.
• Earth model simplification: 1D, spherically symmetric viscoelastic Earth model; lateral viscosity variations are not included, potentially affecting regional GIA responses.
• Non-uniqueness and sparse constraints: Pre-LGM margins (especially Antarctica and parts of North America) are sparsely constrained; MIS 3 chronology relies on alignment with Heinrich Events with multi-kyr uncertainties; flow-direction matching does not uniquely constrain thickness.
• Tuning choices: High lower mantle viscosity (4×10^22 Pa s) selected to balance maximizing LGM ice volume and near-field sea-level fits; alternative rheologies could yield different configurations.
• Exclusions in sea-level budget: Thermal expansion, groundwater storage, and smaller glaciers/ice caps (~3–4 m SLE) are not explicitly included, implying the reconstructed ice volume may be slightly overestimated.
• Antarctic shelf and shear stress parameterizations simplified to avoid excessive thickening; may under-represent some dynamics.