The U.S. Mid-Atlantic sea-level record

The Albemarle Embayment geological record includes interfluvial, estuarine, intertidal and shallow marine lithofacies arranged in depositional sequences that record repeated sea-level highstands dated primarily by optically stimulated luminescence to MIS 5e, 5c, 5a and 3 (ref. 11) (Fig. 1; Supplementary Table 1; Supplementary Note 1). We adopt the minimum elevation of terrestrial facies and the maximum elevation of marine facies as upper and lower bounds, respectively, of MIS 5a (∼80 ka) and mid-MIS 3 (50–35 ka) sea level. For the MIS 5a data (Fig. 1; Supplementary Table 1), we bound a cluster of sea-level data from 2.5 to 7 m in agreement with previous assessments of sea-level records in the region20. We assume that rare terrestrial markers found at elevations below this range do not represent a constraint on the MIS 5a highstand, but rather a lower sea level reached during late MIS 5a or MIS 4. Furthermore, calculations described below (and detailed in Supplementary Note 3) demonstrate that RSL predictions for MIS 3 are relatively insensitive to the height of sea level during MIS 5a.

For the MIS 3 interval spanning 50–35 ka, three marine indicators constrain RSL to be above −0.9, −3 and −2 m (ref. 11). We thus adopt the elevation of the highest of these marine indicators, −0.9 m, as the lower bound. Regarding the upper bound, three terrestrial indicators, with ages between 50 and 35 ka, show a consistent constraint on the sea-level highstand of 1 m. Two terrestrial indicators dated to earlier in this time window13 are found at lower elevations, however these may represent deposition during a time of rising sea level rather than during the peak sea-level highstand. We adopt the terrestrial indicators at 1 m as the upper bound on sea level, yielding a range of −1 to 1 m. We apply an elevation error of ±3 m that reflects reconstructed paleo-tidal range for the region that may have been up to three times greater than the present amplitude of ∼1 m (ref. 21). The geological sea-level constraints we adopt below (for example, Fig. 2) incorporate these broader uncertainties.

Figure 2: Global mean sea-level curves and relative sea-level prediction. (a) Global mean sea-level curve for Version 1.2 of ICE-5G (dotted black line) and the ICE PC (blue line) ice histories. The GMSL curve for the ICE PC2 history is identical to the curve for ICE PC . (b) Relative sea-level predictions for the reference site in the Albemarle Embayment based on the ICE-5G (dotted black), ICE PC (solid blue) and ICE PC2 (dashed blue; ice-free eastern sector of the Laurentide Ice Sheet from 80 to 44 ka) ice histories. Orange rectangles span the observational constraints on peak MIS 5a and 3 sea level including the ±3 m paleotidal uncertainty (see Fig. 1). For MIS 5a, this range is −0.5 to 10.5 m, and at 44 ka during MIS 3 the range is −4 to 4 m. The grey-shaded region spans the MIS 3 time interval examined within the present analysis. Full size image

Models of glacio-isostatic adjustment

Ice sheet growth and melt produces a complex spatio-temporal pattern of sea-level change22. To predict the present elevation of sea-level markers, we perform calculations based on the sea-level theory and pseudo-spectral algorithm described by Kendall et al.23 with a spherical harmonic truncation at degree and order 256. The calculations include the impact of rotation changes on sea level24, evolving shorelines and the migration of grounded, marine-based ice23,25,26,27. We report RSL predictions at a representative site within the Albemarle Embayment (white star on inset of Fig. 1) for a representative time (44 ka) within the middle of MIS 3 (50–35 ka). This representative site lies within the latitudinal range of the reported geological sea-level markers used to define the bound on local peak MIS 3 sea level (Fig. 1). We have found that RSL highstand predictions for this reference site differ from field locations by less than 0.5 m. We deem simulations acceptable if they satisfy the aforementioned bounds for both MIS 5a and MIS 3 (Fig. 2b, orange rectangles).

Our numerical predictions require models for Earth’s viscoelastic structure and the history of global ice cover. We begin by adopting an Earth model with upper and lower mantle viscosities of 0.5 × 1021 Pa s and 1.5 × 1022 Pa s, respectively; this radial profile is consistent with inferences based on globally distributed ice age data sets28 and geological data along the U.S. mid-Atlantic29,30. Our initial GIA calculation adopts Version 1.2 of the ICE-5G ice history, characterized by a GMSL fall from −87 to −100 m throughout MIS 3 (ref. 8) (Fig. 2a; dotted black line); in this calculation, we make the standard assumption that, for any pre-LGM time step, the geometry of global ice cover was identical to the post-LGM ice distribution with the same GMSL value31. We explore alternatives to the GMSL history, ice geometry and viscosity profile in the discussion below. Using the combination of the ICE-5G model and Earth structure described above, we predict mid-MIS 3 sea level (at 44 ka) at the Albemarle Embayment reference site to be –67 m (Fig. 2b; dotted black line), grossly misfitting (by ∼70 m) the observational constraints (Fig. 1). The misfit is ∼25 m for the MIS 5a record (Fig. 2b). The level of misfit to the MIS 3 record highlights the enigmatic nature of the sea-level record in Fig. 1 and motivates the present study.

Many previous inferences of GMSL during the last glacial phase, particularly MIS 3, reconstruct higher peak sea level (smaller global ice volume) than adopted in the ICE-5G history9,32,33. To proceed in our analysis, we revise the ICE-5G ice history on the basis of results from two recent GIA analyses. First, following the Pico et al.10 analysis of sediment core records from the Bohai Sea, peak GMSL during MIS 3 is placed at −37.5 m at 44 ka. Second, we adopt GMSL values of −15 and −10 m for MIS 5a and 5c, respectively; these values are within bounds (5a: −18 to 0 m, 5c: −20 to 1 m) derived by Creveling et al.30 on the basis of globally distributed sea-level markers from both periods. The GMSL curve for the revised ice model, ICE PC , is shown in Fig. 2a (blue line). The RSL prediction based on this model (Fig. 2b, solid blue) maintains the assumption that the pre-LGM global ice geometry is equivalent to the deglacial phase whenever GMSL values are equal. This prediction is consistent with observational constraints for the MIS 5a highstand, but misfits MIS 3 data. Notably, this Earth–ice model pairing predicts a peak RSL of −12 m at 44 ka, well below the observational bounds of −4 to 4 m.

We performed a suite of simulations to explore the sensitivity of our predictions to the adopted ice history. Specifically, we generated 100 synthetic ice histories in which we varied GMSL randomly across the glacial phase but confined GMSL to −37.5±7 m at 44 ka (ref. 10), and to −15 m and −10 m, for MIS 5a and MIS 5c, respectively (Fig. 3a, blue lines; see Methods for detailed ice history construction). Using these ice histories, we predicted RSL at the reference site within the Albemarle Embayment using the standard treatment for the pre-LGM ice distribution; that is, this distribution matches the post-LGM geometry when the GMSL values are the same (similar to ICE PC ). In this case, the predicted RSL ranges from −26.5 to −7.5 m at 44 ka (blue lines, Fig. 3b), and thus all 100 simulations predict a peak RSL that falls outside the observational constraints.

Figure 3: The effect of varying ice history on relative sea-level prediction. (a) Global mean sea level curves for 100 randomly generated ice histories that pass through −37.5±7 m at 44 ka (blue lines), −15 m at MIS 5a and −10 m at MIS 5c. (b) Relative sea level predictions for the Albemarle Embayment reference site based on the ice histories shown in frame (a). The results of calculations that assume identical pre-LGM and post-LGM ice geometries when global mean sea level values are the same are plotted as blue lines; the calculations that assume an ice-free eastern Laurentide from 80 to 44 ka are shown by black lines. Orange rectangles span the adopted observational constraints on peak MIS 5a and MIS 3 sea level including the ±3 m paleotidal uncertainty. The grey-shaded region spans the MIS 3 time interval examined within the present analysis. Full size image

Revising the geometry of the Laurentide Ice Sheet during MIS 3

We next explored the impact of changing the geometry of the LIS on the local RSL predictions at the Albemarle Embayment. While few field data constrain the evolution of the LIS before the LGM2,4,5, a recent field-based study by Dalton et al.34 suggests that large portions of eastern Laurentia were ice-free during MIS 3 (Fig. 4). This conclusion implies limited or no ice growth from MIS 5a to MIS 3 within large areas of the sector of the LIS closest to the Albemarle Embayment. To investigate the effect of this revised ice geometry on RSL predictions, we constructed an ice model, ICE PC2 , with a GMSL history identical to that shown by the blue line in Fig. 2a, but that was distinguished from the ICE PC history in the following ways: (1) the eastern sector of the LIS is ice-free from 80 to 44 ka, consistent with the conclusions of Dalton et al.34; and (2) the ice removed in this exercise, equivalent to 6.8 m of GMSL, is distributed uniformly over the western sector of the LIS, and the Cordilleran and Fennoscandian Ice Sheets. The latter resulted in an increase in ice thickness of ∼170 m in each region (Fig. 4; See Methods for details on ice model construction). The post-LGM ice geometry remains identical to the ICE-5G and ICE PC models, and thus we no longer assume that global ice geometry prior and subsequent to the LGM were identical whenever the GMSL values match. The simulation, based on this revised ICE PC2 ice model and the viscoelastic Earth model discussed above, predicts a RSL of −3 m at 44 ka, consistent with the sea-level record at the Albemarle Embayment (Fig. 2b).