Sea-level records from the U.S. mid-Atlantic constrain Laurentide Ice Sheet extent during Marine Isotope Stage 3

The U.S. mid-Atlantic sea-level record is sensitive to the history of the Laurentide Ice Sheet as the coastline lies along the ice sheet's peripheral bulge. However, paleo sea-level markers on the present-day shoreline of Virginia and North Carolina dated to Marine Isotope Stage (MIS) 3, from 50 to 35 ka, are surprisingly high for this glacial interval, and remain unexplained by previous models of ice age adjustment or other local (for example, tectonic) effects. Here, we reconcile this sea-level record using a revised model of glacial isostatic adjustment characterized by a peak global mean sea level during MIS 3 of approximately −40 m, and far less ice volume within the eastern sector of the Laurentide Ice Sheet than traditional reconstructions for this interval. We conclude that the Laurentide Ice Sheet experienced a phase of very rapid growth in the 15 kyr leading into the Last Glacial Maximum, thus highlighting the potential of mid-field sea-level records to constrain areal extent of ice cover during glacial intervals with sparse geological observables.

The data of MIS 5a age from the mid-Atlantic coast of the U.S. suggest that local relative sea level (RSL) was between +2.55 and +7.5 m. Ages determined by Parham et al. 1 include samples FNB-02 and GH10, which are dated to 84.4 ± 8.6 ka and 81.5 ± 6.9 ka, at elevations of 2.55 m and 7.5 m respectively. FNB-02 is a marine limiting indicator, corresponding to a shelly muddy sand, interpreted as shallow marine. GH10-01 is a terrestrial-limiting indicator in a sand facies, interpreted as "shoreline complex" 1 . While other sea level markers are lower than this range during this time interval, given the associated large age uncertainties, we assume the higher range represents the MIS 5a highstand. Moreover, this sea-level highstand is in agreement with previous studies of MIS 5a RSL in the region 3 .
We add an uncertainty of ± 3 m to the elevation of these sea-level indicators since the tidal range may have been different in the past relative to the Holocene. For instance, Hill et al. 4 (2011) showed that the tidal range of the mid-Atlantic U.S. coast might have been 2-3 times the magnitude at 9 ka. Present day tidal range in the region is ~1 m 5 , and we thus assign a maximum value of ± 3 m to account for the possible paleo tide range.

Supplementary Note 2: Decomposing Sea Level Change into Deformational and Direct Gravitational Effects
Relative sea level (RSL) can be decomposed into global mean sea level (GMSL) plus a term associated with glacial isostatic adjustment (GIA): and the GIA term can be further decomposed into sea level contributions from crustal deformation (including the gravitational perturbation associated with this deformation; R) and the direct gravitational effects of the surface mass load (GD): Figure S1 shows the decomposition of RSL predictions at our reference site in the Albemarle Embayment region (white star, Figure 1) into the three terms GMSL, R and GD. Results are shown for simulations based on ice histories ICEPC and ICEPC2. Delaying the growth of ice in the eastern sector of the Laurentide Ice Sheet (LIS) until MIS 3 (44 ka), suppresses deformational effects, yielding a higher RSL prediction for the ICEPC2 model at 44 ka relative to the ICEPC simulation.
The relative contribution of crustal deformation to predicted RSL will be a function of the Earth model, and in particular the values adopted for lithospheric thickness and lower and upper mantle viscosities. Using the Earth model described in the main text, we noted that crustal deformation dominates the RSL signal compared to direct gravitational effects. We explore the sensitivity of this decomposition to variations in Earth structure by running simulations for additional Earth models. In Supplementary Figure 2

Supplementary Note 3: Sensitivity to ice history and geographic distribution
To test the sensitivity of our simulations to the GMSL values at MIS 5a and MIS 5c, we constructed ice histories as described above, with the exception that we allowed GMSL to range from -16 to 0 m at 80 ka (MIS 5a), and from -20 to 0 m at 100 ka (MIS 5c). We ran simulations with 100 of these ice models and found that the predicted RSL at 44 ka is perturbed by no more than ~0.7 m.
Next we explored the sensitivity of sea-level predictions in the Albemarle Embayment region to variations in the ice cover over the eastern sector of the LIS ( Figure S3). Specifically, we: (1) shifted the northern latitudinal boundary of the no-ice zone in eastern Laurentia from 60 °N (ICEPC2) to 57°N (ICEPC3) and 55°N (ICEPC4); and (2) modified the location of ice in eastern Laurentia, such that Newfoundland and Northern Quebec are glaciated (Geometry 1,2, and 3). For this series of ice geometries, we include one ice history where Baffin Island is glaciated (Geometry 1), and two where the southern portion of the island is not (Geometries 2 and 3). The latter two ice geometries are distinguished by their latitudinal bounds: in Geometry 2, ice cover extends to 51°N in the eastern sector of the LIS, whereas in Geometry 3 it extends to 52°N. Supplementary Table 2 gives the change in ice volume compared with the standard assumption of ice distribution for each ice model listed, and the resulting relative sea level predicted. 6.53 -3.8

Supplementary
Supplementary Figure 3 | Geographic distribution of ice cover for various ice models. A. ICEPC (standard ice distribution, as in Figure 2 and 3B) B. ICEPC2 (alternate ice distribution, as in Figure 2 and 3B). ICEPC3 and ICEPC4 are distinguished by the latitudinal extent of ice cover over the eastern sector of the LIS (see text). C. Geometry 1. D. Geometry 2 & 3 (these two models are distinguished by the latitudinal extent of ice cover -see text). In frames B-D, the difference in ice volume between the given ice model and model ICEPC in eastern Laurentia is distributed uniformly over the LGM extent of western Laurentia, the Cordilleran Ice Sheet, and Fennoscandia.

Supplementary Note 4: Earth model sensitivity
We ran a large suite of simulations in which we adopted ice histories with the GMSL curves shown in Figure 3A. The calculations, shown in Figure 3B, were based on an Earth model with upper and lower mantle viscosities of 0.5 x 10 21 Pa s and 15 x 10 21 Pa s, respectively. These runs yielded predicted peak RSL values within the Albemarle Embayment region during MIS 3 that ranged from -28.6 to -5.7 m when pre-and post- LGM ice geometries were assumed to match when GMSL values were equal (as in ice model ICEPC) and -14.8 m to 2.4 m when it was assumed that the eastern sector of the LIS remained ice free in an extended period leading to MIS 3 (as in ice model ICEPC2). When we reran these calculations using an Earth model in which the lower mantle viscosity was increased to 2 x 10 22 Pa s, the above ranges were perturbed upwards by ~1.5 m. Decreasing the lower mantle viscosity to 10 22 Pa s, perturbed the above RSL range downwards by ~3 m (Supplementary Figure 3).
To explore the sensitivity of our results to other Earth model parameters, we ran simulations using a suite of Earth models where we varied values of the lithospheric thickness and upper and lower mantle viscosity. The resulting RSL predictions for MIS 3, based on the ice model ICEPC2 (Figure 2A), for the representative site in Albemarle Embayment (white star, Figure 1) are plotted in Figure S5.  Figure 3B). The orange bars at 80 ka and 44 ka show the observational constraints ( Figure 1). All calculations are based on ice models that assume that the eastern sector of the LIS remained ice free in the period 80-44 ka (similar to the ICEPC2 model).
Supplementary Figure 5 | Relative sea level predictions at location shown by white star (Figure 1) are based on the ice history ICEPC2 (Figure 2A) with seven different Earth models (including that adopted in main text). A. Lithospheric thickness is increased to 127 km (pink) and decreased to 72 km (dark blue). B. Upper mantle viscosity is increased to 1 x10 21 Pa s (pink) and decreased to 0.3x10 21 Pa s (dark blue). C. Lower mantle viscosity is increased to 30 x10 21 Pa s (pink) and decreased to 5 x10 21 Pa s (dark blue). The Earth model adopted in the main text is shown by the white circle, and the observational constraints on the MIS 3 sea-level highstand is shown by the orange box ( Figure 1).