Sea-level trends across The Bahamas constrain peak last interglacial ice melt

Contributed by Maureen E. Raymo, June 7, 2021 (sent for review January 11, 2021; reviewed by Fiona D. Hibbert and Glenn A. Milne)
August 9, 2021
118 (33) e2026839118


The seas are rising as the planet warms, and reconstructions of past sea level provide critical insight into the sensitivity of ice sheets to warmer temperatures. Past sea level is reconstructed from the geologic record by measuring the elevations of fossilized marine sediments and coral reefs. However, the elevations of these features also record local uplift or subsidence due to the growth and decay of ice sheets since the time of deposition. We compare observations of paleo sea level across the Bahamian archipelago with a range of Earth deformation models to revise estimates of the last interglacial global mean sea level. Our results suggest that polar ice sheets may be less sensitive to high-latitude warming than previously thought.


During the last interglacial (LIG) period, global mean sea level (GMSL) was higher than at present, likely driven by greater high-latitude insolation. Past sea-level estimates require elevation measurements and age determination of marine sediments that formed at or near sea level, and those elevations must be corrected for glacial isostatic adjustment (GIA). However, this GIA correction is subject to uncertainties in the GIA model inputs, namely, Earth’s rheology and past ice history, which reduces precision and accuracy in estimates of past GMSL. To better constrain the GIA process, we compare our data and existing LIG sea-level data across the Bahamian archipelago with a suite of 576 GIA model predictions. We calculated weights for each GIA model based on how well the model fits spatial trends in the regional sea-level data and then used the weighted GIA corrections to revise estimates of GMSL during the LIG. During the LIG, we find a 95% probability that global sea level peaked at least 1.2 m higher than today, and it is very unlikely (5% probability) to have exceeded 5.3 m. Estimates increase by up to 30% (decrease by up to 20%) for portions of melt that originate from the Greenland ice sheet (West Antarctic ice sheet). Altogether, this work suggests that LIG GMSL may be lower than previously assumed.
The geologic record suggests that peak global mean sea level (GMSL) was +6 to 9 m higher than today during the Eemian or last interglacial [LIG; 130 to 116 ka (13)]. However, global mean temperature was only slightly warmer than preindustrial values, likely due to greater high-latitude summer insolation at that time [+1°C to 2°C (4)]. If such modest warming can cause a +6- to 9-m rise in sea level, then the sea-level response to ongoing and future warming will likely be drastic (5, 6). While +6 m of GMSL rise could be explained by thermal expansion of the oceans and the near-complete melting of the Greenland ice sheet (4), significant melting of ice on Antarctica must occur to reach +9 m (7). This significant uncertainty demonstrates that improving LIG GMSL estimates is needed to better understand past and hopefully, future ice sheet responses to warming.
Past GMSL estimates are derived by measuring the elevation and age of marine sediments that formed at or near sea level on passive margins and correcting those elevations for glacial isostatic adjustment [GIA; the elastic and viscoelastic response of the solid Earth, its gravity field, and rotation axis to the redistribution of ice and water (811)]. On glacial–interglacial timescales, topographic change due to GIA depends on the history of ice distribution on land (12) and the magnitude and rate of Earth deformation under that load [i.e., its rheology (9, 13)]. Changes in topography and hence, local sea level also can be caused by cooling of the lithosphere and mantle convection [dynamic topography (1418)], tectonics, or deformation associated with sediment redistribution.
For the LIG, uncertain knowledge of Earth’s internal viscosity, past ice sheet history, and dynamic topography hinder the determination of a precise and accurate estimate of GMSL (1820). Specifically, while the magnitude of postglacial rebound in the vicinity of former ice sheets has been used to constrain the viscosity of Earth’s interior (13, 21, 22), large uncertainties persist concerning both the most appropriate one-dimensional (1D) and the global three-dimensional (3D) viscosity structure of the mantle. When reconstructing ice history, it is often assumed that the penultimate glacial maximum (marine isotope stage or MIS 6) ice sheet history (e.g., location and mass of ice sheet loading) was very similar or identical to that of the last glacial maximum (LGM; MIS 2). However, even small differences in the ice sheet history between MIS 2 and MIS 6 could change GMSL inferred from tropical regions by several meters (19, 20). Lastly, models and observations of dynamic topography suggest that paleosea-level indicators in a single location could uplift or subside by several meters since the LIG, even on passive margins (18, 23).
The Bahamian archipelago is a tectonically stable region that has three important features that make it ideal for reducing GIA uncertainty and improving estimates of past GMSL. First, GIA models predict that this region, on the peripheral forebulge of the Laurentide ice sheet (LIS), should have large lateral variations in relative sea level (RSL) (2) (Fig. 1), and the magnitude and shape of these variations depend on mantle viscosity and past ice history. Observations of contemporaneous paleo sea level can therefore directly constrain these GIA model choices. Second, subsurface core data (24) can be used to constrain long-term subsidence and regional deformation that may have occurred due to dynamic topography or other slow processes such as tectonics, karst isostasy, sediment loading, or the cooling of the lithosphere. Third, the sediments of The Bahamas are excellent archives of past sea level (2529). In this paper, we compare our observations and existing observations of paleo sea level across the archipelago to a range of GIA model predictions to identify those GIA model inputs that are consistent with lateral variations in the sea-level observations; 576 unique GIA models are combined based on how well the model fits lateral variations in regional sea-level data, and this weighted GIA correction is used to revise estimates of GMSL during the LIG.
Fig. 1.
(A) Modeled sea level at the start of the LIG relative to today assuming no change in GMSL during the LIG. This number is equal to the elevation relative to sea level today at which one would find a sediment that records sea level from the start of the LIG if GMSL did not exceed present values. (B) Modeled sea level at the end of the LIG relative to today assuming no change in GMSL during the LIG. (C) Change in modeled local sea level for Great Inagua and Great Abaco. The GIA model used here has a lower-mantle viscosity (LMV) of 5 ×1021 Pa s, an upper-mantle viscosity of 0.5 ×1021 Pa s, and a lithosphere thickness of 96 km; the MIS 6 ice sheet size and geometry are equal to the LGM. Materials and Methods has more details.

Results and Discussion

Geology of the Bahamian Archipelago.

The sedimentary rocks that make up the islands of the Bahamian archipelago are primarily narrow aeolian or storm ridges that accumulated above sea level and broad limestone flats that accumulated underwater during times in the past when the sea level was higher than today. The elevations of these two types of sediments provide upper and lower limits, respectively, for sea level in the past (25, 3034).
The narrow ridges run along the windward (northeast) shorelines of most islands and are the highest landforms in the archipelago, commonly reaching elevations of +20 to 30 m. Occasional gaps in these shore-parallel ridges contain slightly lower–elevation U- or V-shaped ridges that point and extend in the landward direction (32, 33, 35). The limestone flats are typically on the western side of the islands at elevations of less than 7 m above modern sea level. Vertical exposures of these limestone flats along coastal cliffs or in sinkholes exhibit unambiguous marine features, such as fossil corals and traces of other marine organisms.
These exposures of LIG marine sediments can typically be described as shallowing-upward parasequences (25) where, over several meters, a decrease is observed in the apparent water depth in which the sediments were accumulating. The lowermost sediments in a parasequence often contain corals or burrow-filled carbonate sands and muds with marine fossils. Moving upward, mud and fine sand are rarer, replaced by coarser carbonate sands that are typically cross-stratified, which signifies transport by currents and waves. These cross-stratified sand units transition upward into well-sorted sands with nearly parallel millimeter-thin beds that dip toward the paleocoastline (27, 36, 37). This transition sometimes is marked by tabular intraclasts with worn or abraded surfaces. Such intraclasts are found along beaches in The Bahamas today, which also consist of nearly parallel-bedded and well-sorted sands. In most locations, the top of the parasequence contains cross-stratified sand with abundant root casts and fossilized land snails, an indication of aeolian or terrestrial deposition. The contact between these uppermost aeolian units and the lower-lying beach sediments is a snapshot of past sea level. More specifically, this contact marks the landward swash limit of constructive waves or the elevation of the ancient ordinary berm (38).
Modern observations suggest that the shallowing component of the shallowing-upward parasequence is likely a record of rapid accumulation of sediment, rather than a change in sea level. These deposits typically form along the margins of islands where the rate of sedimentation outpaces the generation of accommodation space (by subsidence, for instance). As sediment fills the available space, there is a coinciding decrease in water depth. For example, many Holocene beaches in The Bahamas are currently prograding seaward even though there has been sea-level rise in the region throughout the Holocene (39) (presumably due to GIA). A stratigraphic section through a modern beach would reveal sediments that record an upward decrease in water depth. Based on the above observations, we investigated LIG progradational deposits across The Bahamas and use the highest elevation of the contact between aeolian facies and underlying beach facies as an indication of local peak LIG sea level. Because the region can be expected to subside by up to 10 m due to GIA throughout an interglacial, we expect that the highest beach facies on each island will have formed toward the end of the LIG and that the age of these peak local sea-level deposits is similar across the archipelago.

Peak LIG Sea Level across the Archipelago.

An obvious multimeter decrease in the maximum elevation of LIG beach facies is observed across the archipelago from north to south (Fig. 2 and SI Appendix, Fig. S1). For example, on the island of Eleuthera (northern Bahamas), the highest beach facies are found at Whale Point at an elevation of 6.7 m. In contrast, on Great Inagua (southern Bahamas), the highest-elevation marine sediments are part of an aeolian capped shallowing-upward parasequence where the beach to aeolian transition is only 2.3 m above mean tide level (at Matthew Town). In other locations surveyed, the LIG peak local sea level changes gradually between islands, generally decreasing toward the south and southwest (Fig. 2).
Fig. 2.
The highest surveyed elevations of the contact between LIG beach and aeolian sediments from seven surveyed islands are annotated on a map of The Bahamas. Orthographic renders of six outcrops illustrate the typical shallowing-upward parasequences of The Bahamas.
On a single island, such as San Salvador, beach sediments within LIG parasequences record sea level at various elevations. At the Cockburn Town reef outcrop, the maximum elevation of beach sands is 3.7 ±0.1 m above mean tide level. However, 10 km to the south at Grotto Beach, there are beach sediments as high as 5.8 ±0.22 m. Both outcrops have corals in the basal units that grew during the LIG (27, 36). The difference in the elevation of these beach facies indicates that the sea level of San Salvador Island changed by several meters during the LIG. Local sea level is determined primarily by GIA and changes in GMSL. Each location records sea level at the moment when the accommodation space filled with sediment and transitioned to a terrestrial environment.
To determine peak LIG sea level on islands that were not directly surveyed in the field, peak LIG sea level was estimated from satellite-derived maps of elevation. The islands comprise terrestrial and marine landforms that were classified into three categories: Holocene beaches and marshes, aeolian dune ridges, and marine limestone flats (Materials and Methods). The elevation of marine limestone flats decreases across the archipelago to the south and west, generally matching the pattern observed in the outcrop (SI Appendix, Fig. S1). These elevations provide a lower bound for peak LIG local sea level (that is, these satellite-derived estimates are treated as marine limiting).

GIA in The Bahamas.

The islands of The Bahamas experience significant sea-level rise during interglacial periods due to GIA (Fig. 1) (4043). During glacial periods, the mass of the LIS depresses the underlying North American crust and forces mantle material to flow away from that depression. This flow to the ice periphery causes uplift in many regions, including The Bahamas. During an interglacial, the formerly glaciated regions rebound, and the displaced mantle material migrates northward, causing subsidence of the peripheral high. Thus, at the start of an interglacial, The Bahamas region is elevated relative to today, and sea-level indicators on Great Inagua from this time (at 128 ka for instance) (Fig. 1) would now be about 6 m below present day sea level (assuming that GIA is the only factor causing RSL change).
Furthermore, the predicted magnitude of GIA is greater in the northern islands than in the southern islands (Fig. 1A). As the archipelago subsides throughout the interglacial (Fig. 1), the magnitude of this subsidence decreases with distance from the former ice sheet (southward). Thus, sediments that recorded sea level at the start of the LIG will have subsided more in the north than in the south, and any observed elevation of paleo sea level from that time will also be lower in the north (Fig. 1A). Conversely, sediments (parasequences) that record sea level toward the end of the LIG will be higher in the north relative to coincident sediments in the south (Fig. 1B). Thus, the spatial pattern of peak local sea level, as inferred from parasequences in progradational LIG sedimentary systems, is a simple test on the relative timing of peak local sea level.
The GIA predicted late LIG RSL patterns (higher RSL in the north than the south) match the observed peak LIG sea-level patterns. Both the field and satellite observations show a clear decrease to the south and west in the maximum elevation of LIG marine sediments (Fig. 2 and SI Appendix, Fig. S1). The details of the GIA-induced paleosea-level gradient are sensitive to both the viscosity structure of the Earth as well as the size and distribution of the former (MIS 6) North American ice sheet. Uncertainty in these model parameters translates into uncertainty in the inferred paleo-GMSL predictions derived from regional data.

Determining the Most Likely GIA History.

We determine the accuracy of individual GIA models by identifying those input parameters (Earth viscosity and ice history) that result in RSL predictions that best match the elevation differences of observed sea level across the islands. For example, of a suite of GIA predictions we only retain those that yield RSL predictions late in the LIG that are high in the north and low in the southwest. We quantify how well a specific GIA prediction fits the sea-level trends and assign a relative weight to each model. Uncertainties in observed sea level and the formation age of the observations are fully incorporated into these model weights (Materials and Methods). We highlight that this approach does not make any assumptions about GMSL since sea-level trends across the archipelago at a given time are independent of the GMSL. Likewise if we were to use only two observations with very different ages (e.g., the beginning and the end of the LIG), the data would put no constraint on the GIA models as the sea-level difference could be due to changes in GMSL. The model weights are used to produce a single weighted mean GIA model, which we use to correct our RSL observations.
LIG RSL observations for the region include previously published coral data as well as our island-specific estimates of peak RSL (from both satellite-derived and field surveys). Where absolute ages from corals exist (West Caicos, San Salvador, New Providence, and Great Inagua), the ages of corals span the early interglacial through 119 ky (27, 29, 36, 37, 44). The sands that prograde over these coral outcrops are younger than the corals and often are the seawardmost strands of prograding beach systems. Invariably, these overlying beach deposits have been interpreted by previous studies as the final phase of the LIG, and the elevations of these features generally support interpretations of a late LIG local sea level high in the region inferred from stranded tidal notches (45, 46). Our observations and interpretations of the geology agree with these previous studies (Dataset S1); however, we will argue that most of the increase in local sea level during the LIG results from GIA subsidence.

Bayesian Inversion and a Synthetic Test.

Bayesian Gaussian process regression was used to estimate GMSL during the LIG from a set of our observations and existing Bahamian sea-level observations. Each sea-level observation was corrected for the indicative water depth range of the observation, the long-term subsidence of the archipelago, and GIA using a weighted combination of 576 distinct GIA models. Our fully Bayesian Gaussian process framework provides the statistical tools to incorporate data with very different age uncertainties and indicative water depth ranges to estimate how GMSL (and the uncertainty) changes over time (Materials and Methods). Long-term subsidence is assumed to be the same across the archipelago [2.5 ±1 m per 125 ky (24, 47)]. Large lateral variations in local sea level that are not caused by GIA are not considered in our current analysis, but data constraints on long-term subsidence across the archipelago could be incorporated into future inversions.
To determine the efficacy of our inversion approach, we designed a test where we prescribed a true GMSL curve (“Simulation GMSL”) to generate synthetic observations. We then apply our inversion approach to the synthetic data to infer GMSL, which ideally should reproduce the prescribed true GMSL curve. We emphasize that, like with real observations, the statistical model has no inherent knowledge about what the true GMSL is. First, we randomly selected 1 of the 576 GIA models to represent the true GIA model in this test. Next, we generated a synthetic dataset using the real locations, elevation uncertainties, ages, age uncertainties, and the indicative meaning of each type of geologic feature described in this paper. The only difference between the real and the synthetic data is that the elevation of each synthetic sea-level observation is obtained by combining the synthetic true GMSL and the randomly selected true GIA model.
This synthetic dataset was then analyzed in an identical manner to the real dataset. The model weights show that a range of GIA models, including the true one, perform well in fitting the predicted sea-level gradients as recorded by our synthetic data (SI Appendix, Fig. S2). If we were to assume that all GIA models are equally likely to have generated the synthetic data, the inversion results in a large uncertainty in the inferred GMSL (Fig. 3A). However, when we use the inferred weighting for each GIA model, which favors those models that match the spatial RSL trend, we significantly lower the uncertainty of our inferred GMSL and find that the true GMSL is within 1σ of the most likely GMSL (Fig. 3B). Even if the true GIA correction was known, our inferred GMSL will still not perfectly fit the true GMSL given the spatial and temporal distribution of our observations and their associated uncertainties (Fig. 3C). This test demonstrates that the approach used here can reproduce the underlying GMSL signal and that it indeed leads to a more accurate and precise inference compared with simply treating all GIA models equal.
Fig. 3.
GMSL inference from the synthetic dataset that differs from the real sea-level data analyzed in this paper only in that the elevations were changed to match a prescribed GMSL curve (shown as the dashed black line) and a single GIA model. (A) The estimated GMSL when all 576 GIA models are considered equally probable. (B) The solution where models are weighted by their ability to reproduce the spatial trends across the data (the approach used with the real data in this study). (C) The best solution possible if the true GIA in the region was perfectly known.

Comparison of Bahamas LIG Data and 576 Distinct GIA Models.

We considered a suite of 576 GIA predictions spanning 48 solid Earth configurations, six ice histories, and two penultimate deglacial onset timings (Materials and Methods). The magnitude of predicted gradients in RSL across the archipelago is most sensitive to the prescribed ice loading history and deglacial rate (SI Appendix, Fig. S3). When the modeled penultimate glacial LIS is large (i.e., comparable with its MIS 2 size), the expected RSL gradient in the late LIG data (as observed today) is mostly flat. That is, the solid Earth in The Bahamas would be expected to be in a similar state of isostatic disequilibrium during the late LIG as it is today. When the LIS is smaller or when the deglacial onset is earlier, the region will have relaxed to a state quite different from present day, resulting in mostly north–south gradients in RSL. Generally, models with slower penultimate deglacial rates have higher relative weights (i.e., produce a better fit to the observed RSL gradients) (the right two columns of SI Appendix, Fig. S3). Moreover, with a penultimate glacial ice configuration resembling the LGM, no combination of the modeled lithosphere thicknesses, upper-mantle viscosity, or LMV results in RSL gradients that match the observations from The Bahamas (row 1 of SI Appendix, Fig. S3).
The data–model agreement significantly improves when some of the Laurentide ice is redistributed to the Scandinavian ice sheet (smaller LIS going downward in SI Appendix, Fig. S3). These results are relatively insensitive to changes in lithosphere thickness (even vs. odd columns) (SI Appendix, Fig. S3) or upper-mantle viscosity (y axis) (SI Appendix, Fig. S3), within the range of values tested here. LMV is an important input (x axis) (SI Appendix, Fig. S3). Models with the smallest LIS (row 6 of SI Appendix, Fig. S3) are more consistent with the observations when LMV is in the lower end of the tested range (LMV of 5 to 10 ×1021 Pa s). However, models with the intermediate-sized LIS favor the higher end of the tested range (LMV of 10 to 30 ×1021 Pa s).
The results above suggest that the LIS was significantly smaller during the penultimate glacial maximum than during the LGM. Of the 576 GIA simulations, the sum of model weights with an LGM-like (specifically, ICE-6G from ref. 48) configuration for the penultimate glacial is significantly less than 0.1% (SI Appendix, Fig. S3). We find (at a 98% probability) that the LIS was between 16 and 23 m of sea-level equivalent (SLE) ice volume smaller during the penultimate glacial maximum than during the LGM. We find almost no probability of a larger (<16-m SLE reduction) or smaller (>23-m SLE reduction) LIS during the penultimate glacial maximum. Thus, paleosea-level data from The Bahamas strongly demonstrate that the LIS was smaller during the penultimate glacial maximum, consistent with evidence for either less overall ice volume at the time (as suggested in ref. 20) or larger ice sheets in Europe and Russia (4951). While there is evidence that the LIS covered a similar extent to the LGM during the penultimate glacial cycle (51), the specific timing and volume of that maximum ice load remain uncertain. Improvements to the specific timing of this LIS ice load history will be critically important for future GIA modeling and sea-level reconstruction efforts.
In addition to the large suite of GIA models that assume radial symmetry in viscosity (1D GIA models), we also included three model simulations with lateral variations in viscosity (3D GIA models) (Materials and Methods). Our inversion showed that these models did not yield a good fit to the observations, which could be driven by an incorrect viscosity structure or an incorrect ice history. A larger suite of 3D GIA models would be desirable, yet exploring the parameter space systematically is very computationally expensive and beyond the scope of this study. Yet, we did use the available 3D GIA simulations to test two first-order questions. First, we investigated whether the gradient that is predicted by the 3D GIA model could be matched by a 1D GIA prediction. We found that, in all cases, there were 1D models that fit the 3D GIA gradient across the region (the difference between the 3D GIA model and the most similar 1D GIA model was, on average, 10 cm). Second, we investigated whether the 3D GIA model could bias our result by repeating the synthetic test using one of the 3D GIA models as true input. We found that the true prediction is always within the 95% GMSL envelope and mostly within the 66% envelope of the inferred GMSL. We therefore conclude that 3D viscosity variations are unlikely to largely affect our results but also acknowledge that this result is based on a small set of 3D GIA simulations.


Gaussian process regression was used to interpolate in time across the set of GIA corrected sea-level observations to obtain our estimate of GMSL from this region (Materials and Methods and Fig. 4). Our posterior GMSL curve is lower than previous estimates throughout the LIG, with mean sea-level values generally between 0 and 4 m. We find a 95% probability that GMSL during the LIG period peaked at least 1.2 m higher than today, and it is very unlikely (5% probability) to have exceeded 5.3 m (Fig. 4B). It is very likely (95% probability) that this peak GMSL occurs before 121 ka, with the most likely peak occurring at 127 ka (Fig. 4C).
Fig. 4.
GMSL estimates for the LIG derived from analysis that includes a probabilistic assessment of 576 GIA models, previously published Bahamian archipelago corals, and the peak local sea-level observations reported in this paper. (A) Inferred GMSL change throughout the LIG; 68% of the solutions fall within the inner, darker envelope, and 95% of solutions fall within the outer, lighter envelope. (B) The distribution of maximum GMSL (at any point during the LIG). (C) Distribution in time for peak GMSL.
Our analysis supports previous geological arguments for a small oscillation in local Bahamian sea level within the LIG. We find it very unlikely that there was an oscillation of GMSL within the LIG that exceeded 1.8 m, and the average inferred oscillation is 0.5 m. These results are consistent with recently published work arguing for no significant ice sheet regrowth within the LIG (52). On Great Inagua, West Caicos, and San Salvador Island, LIG coral reefs are vertically separated into two units by a discontinuous marine-erosion surface that has been interpreted as a rapid GMSL oscillation (27, 37, 53, 54). Ref. 29 presented coral ages from the upper and lower units on San Salvador and Great Inagua and argued that the erosion occurred at 125 ky and may have lasted up to 1.5 ky (29, 55). RSL estimates for these sites are consistent with this past work (Fig. 5). All sites experience minor erosion followed by relatively stable sea level during the middle of the interglacial, followed by a later rise in local sea level. These two distinct phases of coral growth during the LIG in The Bahamas may reflect an early LIG GMSL peak (Fig. 4) and a regionally amplified (from GIA) late LIG peak (Fig. 5).
Fig. 5.
(AM) Inferred site-specific relative sea level; 68% of the solutions fall within the inner, darker envelope, and 95% of the solutions fall within the outer, lighter envelope. Coral data are plotted at the elevation of the coral sample today, and the RSL uncertainty (the range of possible water depths that the corals lived in) comes from previous work (27, 29, 44, 55, 56). Satellite-derived high-stand data are plotted at the observed elevations, which are interpreted as the minimum possible RSL (marine limiting). Outcrop observations of the boundary between aeolian sand and underlying beach facies are plotted at the observed elevation offset downward by estimates of the elevation of the landward swash limit of constructive waves (the ordinary berm) (38). Ordinary berm paleoelevations were estimated using models of site-specific modern day wave run-up and mean high higher water (57).
Before comparing our GMSL estimates with previously published values, we need to consider one additional effect. Our GIA correction does not assume any ice melt during the LIG beyond the present day ice configuration. However, as ice sheets melt, they cause an additional distinct spatial pattern of sea-level rise due to gravitational and viscoelastic deformation effects (58). For example, if all of the inferred positive GMSL during the LIG is derived from melting the Greenland ice sheet, sea level in The Bahamas will be around 25% lower than the global mean. In other words, an inferred GMSL value of 3 m in The Bahamas that is entirely sourced from the Greenland ice sheet would translate into 4 m of GMSL. Alternatively, if all of the positive GMSL is derived from the West Antarctic ice sheet, sea level in The Bahamas will be around 25% higher than the global mean.
High GMSL early in the LIG followed by a more stable climate is generally consistent with reconstructions of LIG atmospheric temperatures in Antarctica and Greenland (59, 60). While the relationship between temperature and ice sheet volume is not straightforward, coupled regional climate and ice sheet models of the Greenland ice sheet suggest that the maximum contribution to sea level from Greenland occurred at 121 ky and ranges from +2.1 (61) to +5.1 m (62). If the lower end-member is correct, we would expect to see around +1.7 m in The Bahamas, which is in line with our estimates at 121 ky, leaving little to no contribution from the Antarctic ice volume beyond its present value. However, if the higher Greenland contribution end-member is correct, then larger than present ice sheets must have existed in Antarctica (or elsewhere) at that time. The reconstructions of Greenland contributions to GMSL cited above are not consistent with the peak LIG GMSL at 127 ky (Fig. 4C), which suggests that the Antarctic ice sheet was smaller than today at this time. If susceptible sectors of the West Antarctic ice sheets collapsed, they would lead to 3.3 m of ice equivalent sea-level rise (63), which would be observed as 4.1 m in The Bahamas, an estimate that is in line with the upper end of our inference assuming little to no contribution from the Greenland ice sheet at this time. This conclusion leads us to hypothesize that the West Antarctic ice sheet collapsed early in the LIG (127 to 126 ka) but then started to expand, while the Greenland ice sheet was shrinking between 126 and 121 ka. Such asynchronicity is supported by the observation that peak temperature anomalies on Antarctica [128.7 ka (59)] occur earlier than on Greenland [126 ka (60, 62)]. Lastly, we note that ocean thermal expansion and mountain glaciers contributed to sea-level rise during the LIG but unlikely exceeded 0.4 ± 0.3 (4) and 0.3 ± 0.1 m (64), respectively.
Our estimates are at odds with the inference by ref. 1, which found that GMSL very likely (95% probability) peaked at least 6.6 m above present day values. Kopp et al. (1) included 105 sea-level index points, and there is no one subgroup of indicators (corals, erosional, facies, isotopic) that, if removed, lowers their inferred GMSL significantly. We note that their compilation includes several data points from The Bahamas, and they use a similar long-term subsidence rate for those data to what is used here. The inferred GMSL can nonetheless vary significantly depending on the GIA correction, which we constrain with gradients in LIG RSL across The Bahamas. To determine these constraints, we have assumed that no other process can explain spatial variations in RSL. While there is good agreement between RSL observations and GIA predictions, which may support a sole GIA origin for the observed lateral gradients, later variability in the long-term subsidence or GIA deformation beyond the range of Earth and ice models tested here can affect our inference. Our estimates are also lower than inferences from coral data from the Seychelles, which indicate that sea level was 7.6 ± 1.7 m higher than present at 125 ka (3). A higher sea level in the Seychelles could be explained by local uplift; however, there is no evidence to support or repudiate this. Lower GMSL is inferred based on a record of phreatic overgrowths on speleothems in a coastal cave on Mallorca (Spain), but uncertainties associated with the GIA process are on the order of multiple meters at this location, making it difficult to draw clear conclusions. Refs. 65 and 66 also argued for lower GMSL that peaked around 4 m during the LIG based on ice–ocean–atmosphere modeling, which is consistent with our findings.
We note that the GIA model weights constrained by Bahamian data are only partly useful for other locations. For example, given real 3D variations in Earth’s viscosity, other locations will require GIA models with a different Earth structure. Our inference about the size of the LIS during MIS 6 will affect the GIA corrections around the globe. However, this inference does not provide great constraints on the extent and ice history of the Fennoscandian ice sheet, which would be needed to improve GIA corrections closer to the Fennoscandian ice sheet.


This work demonstrates that laterally extensive paleosea-level data from a region sensitive to GIA model parameters can be used to better constrain those parameters and thus, provide more accurate estimates of past GMSL. Our documented southwestward decrease in late LIG sea-level indicators across The Bahamas archipelago, apparent in both outcrop and satellite elevation data, can only be reproduced with GIA models that assume a smaller LIS during the penultimate glacial maximum when compared with the LGM (a volume reduction of 16 to 23 m of sea-level equivalence). From an analysis that includes a probabilistic assessment of 576 GIA models, we find that sea level peaked 1.2 to 5.3 m above present levels (95 and 5% probability) and that this peak likely (95% probability) occurred before 121 ka. Our estimate is considerably lower than the currently accepted range of +6.6 to 9.4 (1, 3); however, it is in line with recent climate and sea-level modeling (66). We do not find evidence for a significant late rise of GMSL associated with a putative rapid collapse of ice sheets after prolonged warming. If true, our results suggest that polar ice sheets may be less sensitive to high-latitude warming than currently believed. Importantly, warming during the LIG was likely insolation driven, which may impact each hemisphere at different times. If the melting of Northern and Southern Hemisphere ice sheets was out of sync during the LIG, the individual ice sheets may still be very sensitive to local temperature change. Ongoing and future warming is driven primarily by greenhouse gases, which should affect both hemispheres more equally.

Materials and Methods

Outcrop Surveys.

Six islands that span the geographic range of The Bahamas were surveyed to investigate the maximum elevation of local sea level during the LIG. Each outcrop was extensively photographed using a full-frame digital camera. The high-resolution two-dimensional (2D) photographic data of each outcrop was then transformed into a 3D reconstruction using photogrammetry (from Agisoft Metashape Professional Edition). We measured the latitude, longitude, and elevation of at least five control points at each outcrop with differential global positioning system (GPS). The positions of these control points were identified manually on each photograph, and these positional data were used to georeference the 3D outcrop reconstructions. The relative uncertainty within a single outcrop reconstruction (the measured distance between two points or features) is less than a centimeter. The vertical uncertainty of the reconstruction when compared with measurements of local mean sea level or another outcrop is limited by the precision of the differential GPS and includes uncertainty from both measurements. The combined uncertainty is the square root of the sum of the squares of the measured uncertainty from each observation. The uncertainties reported in Fig. 2 reflect this combined uncertainty for each outcrop.
We provide a general description of the geologic sites used in this study (Fig. 2), which in part, have already been described elsewhere. At each site, we interpret the contact between the beach facies and the overlying aeolianite as the elevation of the landward swash limit of constructive waves [the ordinary berm (38)].

Clifton Pier, New Providence Island.

Exposed along the cliffs of the southwest corner of New Providence Island, there are excellent exposures of LIG shallowing-upward parasequences. The lower 2.5 m contain large-scale (>50-cm) trough cross-stratification and intraclast breccia blocks (roughly 50-cm diameter). The intraclasts disappear upward, and the cross-bed scale decreases slightly up to about meter 4.5, where the cross-stratification transitions into 10- to 20-cm tabular cross-beds. Farther upward, sparse tabular clasts mark the transition from this unit into roughly 2 m of well-sorted sands with nearly parallel millimeter-thin beds that dip toward the paleocoastline. The outcrop is capped by very fine sand that has abundant fossilized root casts. This unit is described in detail by ref. 67 and many others (25, 31, 39).

Grotto Beach, San Salvador Island.

At Grotto Beach, the LIG outcrop extends from the islandward limit of the modern beach and stretches across the grotto. Immediately adjacent to the modern beach is a paleopatch reef containing Porites astreoides, Psuedodiploria sp., and Montastrea sp. Around 2.5 m above mean sea level (MSL), these framework corals are capped by a thin layer (10 to 20 cm) of coralline red algae (Neogoniolithon strictum). Above, there are 1.5 m of cross-bedded and very poorly sorted (medium- to pebble-sized) angular carbonate sand. This sand unit is capped by intermittent breccia intraclasts that transition upward into well-sorted sands with nearly parallel millimeter-thin beds that dip toward the paleocoastline. These nearly parallel-bedded sands persist up to at least 5.8 m where bedding becomes difficult to discern. The uppermost sediments here have abundant root casts and are mostly very fine sand (likely aeolian). Ref. 36 contains a detailed description and interpretation of this outcrop.

Long Cay, Crooked Island.

Roughly 30 km south of Landrail Point along the western margin of Long Cay, there is outcrop of an LIG progradational strandplain that extends for 4.5 km. The lower 3 m of this sequence contain large-scale (50-cm) cross-stratification and abundant Ophiomorpha trace fossils (likely shrimp burrows). This lower unit contains occasional thumb-sized coral fragments (Siderastrea radians). There is an irregular and gradual contact between this lower facies and an overlying facies that consists of approximately 2 m of well-sorted sands with nearly parallel millimeter-thin beds that dip toward the paleocoastline. The uppermost part of this unit contains extremely well-preserved spherical beach fenestrae (the remnants of air escaping when waves crashed on the ancient beach). Above these fenestrae-bearing beach facies, there are large meter-scale fossilized aeolianite dunes. The dune height varies along the 4.5-km coastal outcrop from half a meter to more than 3 m.

Boat Cove, West Caicos Island.

This outcrop is part of an extensive exposure of an LIG reef that spans most of West Caicos Island. At this location, the lowermost 2 m are a heavily recrystallized framework of in situ corals (Acrapora palmata and Psuedodiploria sp.). Those corals are capped by 2 m of large-scale (50-cm) cross-stratification and abundant Ophiomorpha trace fossils (likely shrimp burrows). This cross-stratified unit sharply transitions upward into 1 m of well-sorted sands with nearly parallel millimeter-thin beds that dip toward the paleocoastline. The oceanward edge of these nearly parallel beds often contains tabular breccia blocks of the same unit (possibly ancient beach rock breccia). This unit is capped with intermittently spaced (10- to 20-m) large breccia-filled “flower pots.” These features are 1- to 2-m-diameter round bowl-shaped depressions carved downward into the subplanar beds. Each “pot” is filled with bright red limestone that contains clasts of recrystallized gray or black limestone. These features may indicate the past location of large plants or palm trees. Ref. 37 contains more detailed descriptions of the outcrop on West Caicos.

Old Gray’s, Long Island.

Northeast of Anderson’s settlement on Long Island, there is an outcrop of an LIG fossil reef that extends for 400 m. On the northern end of this reef, the modern coastline bends westward, and around that bend, there is a 3-m-high shallowing-upward parasequence. The lower meter contains large-scale (50-cm) cross-stratification. The contact between this lower unit and the overlying unit is difficult to pinpoint. The overlying unit contains well-sorted sands with nearly parallel millimeter-thin beds that dip toward the paleocoastline, and the uppermost beds in this unit contain well-preserved spherical beach fenestrae.

Matthew Town, Great Inagua Island.

In central Matthew Town, the coastal outcrop consists of an LIG shallowing-upward parasequence. The lowermost meter contains abundant large-scale (50-cm) cross-stratification and Ophiomorpha trace fossils. There is an irregular contact between this cross-bedded facies and the overlying facies that consist of well-sorted sands that are arranged in nearly parallel millimeter-thin beds that dip toward the paleocoastline. The top of this well-sorted unit contains a 5- to 8-cm-thick deposit with pervasive and irregular centimeter-scale macroporosity. This porous unit has been interpreted as a wrack line or plant accumulation layer at the contact between beach and dunes (68). The uppermost rocks in this location contain very fine sand and are full of fossilized root casts.

Quantifying mean sea level.

To compare one outcrop with another or to compare an outcrop with digital elevation maps, we developed a mean sea-level reference frame for the surveyed islands (SI Appendix, Fig. S4). There is a long-term tidal and elevation dataset available from Settlement Point on Grand Bahama Island (69). However, sea level relative to the reference ellipsoid measured by GPS varies by tens of meters across the archipelago. Therefore, we deployed tide gauges in each field area to measure water levels over several days to a week. The water-level records were clipped to include roughly the same number of high- and low-tide events, and mean water level was calculated from this subset of the data. The elevation of the water-level logger was measured with differential GPS to determine the elevation of the measured mean water level on the ellipsoid reference frame (specifically, world geodetic system revision 84 or WGS84).
Despite this approach, several days to a week is not enough time to precisely capture seasonal or multiyear signals in sea level and tides. The temperature of the water is also a major seasonal component to sea level, and in nontropical areas, this effect may be as large as 30 cm (70). In the tropical Bahamas where temperatures vary little throughout the year, this effect is likely much smaller (<10 cm).
To assess the higher-frequency mean water-level variability, the mean water-level calculation was bootstrapped with a range of tidal cycle lengths and starting points. This bootstrap results in a 9-cm range in the calculated mean water level for our Nassau water-level data. Without a multiyear tidal record in The Bahamas, measuring the uncertainty in estimating mean water levels from short weeklong datasets is not possible. However, it seems that seasonal signals should be relatively small in the region, and higher-frequency variations from our recorded intervals are similar to the differential GPS uncertainty. Therefore, we assume that the uncertainty associated with measuring mean sea level with data from a small set of tidal cycles is ±10 cm.
To estimate the elevation and uncertainty of mean tide level around the archipelago, the Settlement Point long-term record and our tidal data were combined with the EGM2008 (Earth Gravitational Model 2008) geoid model using a 2D Gaussian process regression. In addition to the tide gauge data, the input data to the Gaussian process regression included EGM2008 geoid sampled in a 0.04 latitude and longitude grid. The uncertainty ascribed to the geoid data was set to the SD of the residuals between the EGM2008 geoid and MSL at Settlement Point, Grand Bahama and our tide gauge stations (σ= 0.76 m). The kernel used for this regression is the sum of two Matèrn covariance functions with ν=32, and the variance and length-scale parameters were fit to the tidal and geoid data prior to calculating the regression. Using this estimated MSL surface, each GPS measurement was adjusted to local mean sea level, and the uncertainty in that MSL estimation was added to the GPS uncertainty in quadrature. Note that the uncertainty in MSL increases away from tidal measurements.

Mapping Landforms across the Archipelago with Neural Networks.

A landform map of aeolian ridges, LIG limestone flats, and Holocene marshes and strandplains was created for the Bahamian Islands and Turks and Caicos. This map is based on a combination of satellite-derived digital elevation data (TanDEM-X from the German Aerospace Center) and our own field surveys of West Caicos, San Salvador, New Providence, Great Inagua, Long Island, and Crooked Island. The TanDEM-X data were down sampled from a pixel resolution of 11 to 27 m for faster computational times. These elevation data were then converted from the provided geoid reference frame (WGS84) to elevation above mean sea level by subtracting our tidal surface map (in the previous section).
This elevation map was converted into a map of landforms using a convolutional neural network. Specific landforms (ridges, flats, marshes) have specific characteristic elevations and geographic patterns. A set of input data was created to teach the network how to distinguish these landforms. To create these input data, we first manually mapped landforms across the region. Ridge features were assumed to be aeolian dunes. Modern depositional systems, such as active strand plains and tidal marshes, were distinguished by eye on satellite images. The remaining land area, which was classified as limestone flats, is a mostly forested surface that tends to have abundant dissolution collapse features (sinkholes). Then, 15,000 random latitude and longitude pairs were selected from each landform class. The input data consist of a 51- × 51-pixel grid of elevation data centered on the latitude and longitude of each pair (pixel dimension 30 m). Ten percent of these input data were held back from the training process. These validation data offer some insight into the ability of the network to perform on data it never specifically learned from. The remaining 90% of the input data are the training data.
The training data were fed into a convolutional neural net with 10 processing layers between the input and the output. First, the observed training data pass through two convolutional layers that each have 32 convolutional kernels (or learned filters). The filter dimension is 3 × 3 pixels. Next, the partially processed information flows into a 2- × 2-pixel down-sampling layer (maximum pooling), after which a 25% dropout is applied. This down-sampled information passes through four more layers that are similar to the first four layers with one modification; the two convolutional layers have 64 kernels each. At this point, the information is flattened into a single array and fed through a fully connected or dense layer that contains 512 nodes (or neurons). The activation function for each node in this layer is a rectified linear unit. After the fully connected layer, another dropout layer is applied with a 50% dropout probability. The final output layer is another fully connected layer, this time with a softmax activation function. The output layer has three nodes, one for each landform classification.
This process transforms a 51- × 51-pixel map of elevation data surrounding a training point into three numbers that are the probability of the input data matching each of the three landform classifications. Before training the network, all of the layer parameters and filters are initialized randomly [Glorot uniform (71)]. In this state, the ability of the network to classify the training data will be as good as randomly drawing classifications. However, by comparing predictions with the correct answer, the backward propagation of errors is used to train the network parameters to improve performance. The network was optimized using an adaptive learning rate algorithm known as Adam (72). After one pass through the training data (an epoch), the network performance jumps from an initial accuracy of 33 to 81% on training data. Training was terminated after 73 epochs, where training accuracy was 99.7% and validation accuracy was 98.22%. This trained network was then used to convert every above sea-level pixel in the TanDEM-X digital elevation dataset into a landform classification. In the end, we produce SI Appendix, Figs. S1 and S2, which show the elevation and distribution of the LIG flats across the entire archipelago. In Results and Discussion, we relate these observations, and the variation in their elevation across the archipelago, to past sea level.
For the seven islands with outcrop surveys, the satellite-derived elevations of LIG marine flats (SI Appendix, Fig. S5, gray distributions) are compared with peak sea level determined from outcrop (SI Appendix, Fig. S5, white distributions). On New Providence Island (Nassau), the elevation of the highest limestone flats is very close to the outcrop estimate of peak LIG sea level, as expected. The elevation of the highest limestone flats on each of the remaining six surveyed islands is generally higher than the corresponding outcrop estimates of peak sea level by a few meters (SI Appendix, Fig. S1). We suggest that this difference exists because the interiors of these six islands are covered in secondary growth blackland coppice forests. On the heavily populated New Providence Island, most forests have been cleared, and the satellite-derived elevations more closely approximate the rocky surface. Because the other islands are forested, satellite-derived elevations represent the canopy top, and the relationship between satellite-derived elevations and LIG peak sea level for these islands should be similar to that observed on the six surveyed forested islands. The outcrop estimate of peak sea level on these six islands is equal to the 50 ±11 percentile of the distribution of limestone flat elevations (SI Appendix, Fig. S1). When this percentile is used to estimate peak LIG sea level around the archipelago, the estimates are consistent with the field observations discussed above. In the analysis, an uncertainty of ±2 m is applied to the satellite-derived estimates to account for the differences between satellite and field observations. Even with this large uncertainty, it is clear that both datasets exhibit a significant regional pattern across the archipelago where peak LIG sea level decreases to the south and west.

GIA Modeling.

To calculate the gravitationally self-consistent sea-level response to changes in ice and ocean load, we solve the sea-level equation following the algorithm in ref. 73. This formalism allows for the migration of shorelines and takes the motion of Earth’s rotation axis in response to surface load changes into account. This model requires as input viscoelastic properties of Earth’s interior that are assumed to be spherically symmetric. We use the density and elastic structure from PREM (74) and vary three parameters for the viscosity structure: upper-mantle viscosity, LMV, and elastic thickness of the lithosphere. The upper mantle is defined as the region below the lithosphere and extending to 660 km. The lower mantle extends from this depth to the core–mantle boundary. We vary the elastic thickness of the lithosphere between 71 and 96 km, the upper-mantle viscosity between 0.3 and 0.5×1021 Pa s, and the LMV between 3 and 40×1021 Pa s.
In addition to the Earth structure, GIA models also require an assumption about the temporal evolution of past ice sheets that drive the crustal/mantle deformation. We use the ICE-6G reconstruction (48) for the last deglaciation and follow the foraminifera isotope-based eustatic curve by ref. 75 prior to the LGM. Based on this curve, the assumed eustatic value at 128 ka is −75 m, which is at odds with coral evidence from the many locations that indicate sea level must have been close to present at that time (3, 29, 76). To account for this discrepancy, we shift the eustatic curve prior to the LIG back by 3.5 ka. This shift allows for a longer interglacial time period without changing the deglaciation pattern of the original curve and places the MIS 6 sea-level low stand at 135.5 ka. However, constraints from Tahitian corals (77) indicate that local sea level was already −85 m or higher by 137 ka and that the deglaciation had started by 142 ka.
To explore the uncertainty associated with this deglacial timing, we run simulations with two eustatic sea level curves, one that has its MIS 6 maximum at 135.5 ka and one that has it at 142 ka (SI Appendix, Fig. S7). We do not include a sea-level oscillation during the deglaciation as was proposed by refs. 77 and 78. In order to investigate deviations relative to present day sea level, we set the excess eustatic sea level during the LIG (assumed here to last from 128 to 117 ka) to 0 m (Fig. 1C).
We start the GIA model runs at MIS 11 (400 ky) following the shifted sea-level curve by ref. 75. To obtain ice sheet geometries from the eustatic curve prior to the LGM, we assume correspondence between ice sheet geometry and eustatic value during the last deglaciation to extrapolate back in time. This (typical) approach results in an ice sheet configuration during MIS 6 that obviously resembles the LGM (MIS 2). However, evidence exists that the Fennoscandian ice sheet was bigger (and consequently, the LIS likely smaller) during MIS 6 (4951). To explore this uncertainty, our GIA models were run with three ice sheet configurations: 1) a reconstruction where the MIS 6 ice sheet configuration is the same as MIS 2 (LIS = 89 m SLE ice volume), 2) the reconstruction of a larger Fennoscandian MIS 6 ice sheet by ref. 49 that is constrained by deformed shorelines and ice margins (LIS = 59 m SLE), and 3) a reconstruction of an even larger Fennoscandian ice sheet from ref. 50 that is derived from coupling climate models to paleoenvironmental proxy data from MIS 6 (LIS = 43 m). These are the same ice sheet configurations that were explored in ref. 19. We explore three additional ice distributions between scenario 1 and scenario 2 with the LIS equal to 81, 73, or 66 m SLE. For all of the ice sheet reconstructions, we assume that the increase in ice mass in the Fennoscandian ice sheet was perfectly compensated by a decrease in ice mass over North America. The contributions from the Antarctic and Greenland ice sheets are not varied (other than for different GMSL scenarios). Their contribution is small and will not lead to significant changes investigated here. However, their melt contribution will lead to distinct patterns over the course of the LIG, which is discussed in the text. The six ice-loading histories are paired with the range of viscoelastic structures and eustatic sea-level curves described above to result in a suite of GIA predictions for The Bahamas. Varying both Earth and ice history leads to 576 different GIA predictions for the Bahamian archipelago.
In addition to the suite of 1D GIA models, we investigated three 3D GIA models. These models use an ice history based on a Fennoscandian MIS 6 ice sheet (49) and a slower deglaciation based on ref. 77. The 3D viscosity structure is derived from the seismic tomography model SL2013sv by refs. 79 and 80 above the transition zone and SEMUCB-WM1 by ref. 81 below the transition zone. Seismic wave speed is converted into viscosity using a series of rheological relationships (82) in which conversion parameters are obtained by fitting first-order geodynamic observations (83). The boundary between the lithosphere and asthenosphere is based on the 1,175°C isothermal surface, and the 3D viscosity structure is imposed on a radial profile 40×1021 Pa s. The three different simulations differ by their radial profile. The first simulation assumes an average lithospheric thickness of 100 km and viscosities of 5×1020 and 5×1021 Pa s in the upper and lower mantles, respectively. The second simulation uses the same viscosity but assumes a thinner average lithospheric thickness of 80 km. The third simulation differs from the first by using a lower viscosity in the lower mantle (1.6×1021 Pa s down to a depth of 1,175 km and 3×1021 Pa s below that). The 3D GIA calculation is performed using the code described in ref. 84. More information on the viscosity structure and ice history can be found in ref. 85.

Estimating GMSL Change during the LIG.

Bayesian Gaussian process regression was used to estimate GMSL during the LIG from published elevations and ages of corals as well as our observations of late LIG peak sea level. Each sea-level observation is corrected for the indicative water depth range of the observation, the long-term subsidence of the archipelago, and for GIA using a weighted combination of 576 distinct GIA models. The Gaussian process framework allows us to incorporate data with very different age uncertainties and indicative water depth ranges to estimate how GMSL (and the uncertainty) changes over time.
GMSL is assumed to be a Gaussian process with a covariance function that is the sum of a radial basis function (RBF) and a white noise kernel. We set a prior distribution for the length scale for the RBF kernel to an inverse Gaussian distribution with μ=2 and λ=5. This length scale effectively ignores changes on short timescales (hundreds of years) and timescales longer than the LIG (>10 ka). A normal Gaussian distribution with μ=0 and σ=5 is used as a prior for the square root of the variance of the RBF function. A half (positive-only) Student’s t distribution with ν=1 and σ=0.1 is used as a prior for the variance of the white noise kernel.
The peak sea-level estimates from the neural network analysis are treated as sea-level marine-limiting points with a uniform flat prior that forces the point underwater. The prior for the age of these sea-level estimates is set to 116 ka plus an inverse Gaussian distribution with μ=2 and λ=5 (in other words, the ages are modeled as late LIG with the peak around 117 ka and an asymmetric higher tail that extends up to about 122 ka). This prior was designed to capture our geologic interpretations of the observations as described in the beginning of Results and Discussion while also agreeing with comparisons with the predicted gradients in GIA throughout the LIG. The elevations of the contact between aeolian and underlying beach facies were interpreted as the elevation of paleo-ordinary berms. This offset from mean sea level was modeled as a normal distribution where the mean ordinary berm elevation was calculated using IMCalc application (57) at each location today (μ is between 0.86 and 1.28 m), and the uncertainty prescribed is equal to the SD of berm elevations (σ=0.17m) across the entire dataset.
The age and uncertainty of the corals are set to the reported values bounded by the start (128 ka) and end (117 ka) of the LIG in our GIA model runs (27, 29, 44, 56). The elevations of the corals are normally distributed about the reported elevation and uncertainties from refs. 55 and 56. The water depth of each coral is modeled separately with an inverse Gaussian distribution, where μ and λ are scaled so that the 95th percentile matches reported maximum water depths. This prior for the coral water depth is consistent with the water depth interpretations of these sites presented in refs. 44 and 55. The long-term subsidence of the islands is considered to be the same across the archipelago and is modeled as a normal distribution centered at 2.5 m with an SD of 0.5 m [the observed rate from cores on Andros Island is 2.5 m per 125 ky (24, 47)].
Posterior distributions were sampled with the No U-Turn Sampler using the python probabilistic programming package PyMC3 (86, 87). The posterior GMSL distribution was calculated for each GIA model separately. The next section will describe how each of these statistical models is combined into a metamodel that contains our best estimate of past GMSL. The GMSL prior is a function of the priors on the covariance and variance parameters of the Gaussian process regression. The prior GMSL will be flat throughout the LIG with normally distributed uncertainty (2-σ)of ±10 m centered about zero (although with data, this prior is centered about the mean of the data).

Combining GMSL Estimates with Bayesian Model Averaging.

Bayesian model averaging was used to determine the relative probability that a given GIA model can explain the observations. For example, a perfect or true GIA model would completely remove all spatial variations in the paleosea-level data (assuming no other spatially variable processes). The statistical framework described above allows us to explicitly include uncertainties in age and sea level while determining the probability that each GIA model could generate the observations (among the 576 compared models). More specifically, we use Pareto smoothed importance sampling to approximate leave-one-out (LOO) cross-validation (88). The LOO cross-validation is used to generate relative weights that describe the probability of each model to have generated the data (89). These weights essentially describe the predictive power of each model (in other words, how well the model performs at predicting observations that are withheld from model training). These calculations were implemented using ArviZ (90). These weights are then used to create a metamodel that combines information learned from each model–data comparison. This metamodel can generally be thought of as a weighted mean. One realization of the posterior in this metamodel is generated by calculating a weighted mean that combines individual posterior realizations from each of the 576 GIA models. Generally speaking, we continue to generate realizations of this posterior until the statistics of interest stop changing significantly. In other words, when the posterior mean and variance start to converge, enough realizations have been generated.

Data Availability

All data, models, and code used in this work are available in GitHub ( and Zenodo (


This material is based upon work supported by NSF Grants OCE-1202632 “PLIOMAX - Pliocene maximum sea level” and OCE-1841888. J.A. and M.E.R. acknowledge support from the Vetlesen Foundation. The research of A.R. on the LIG was funded by Institutional Strategy of the University of Bremen, German Excellence Initiative Grant ABPZuK-03/2014 and the European Research Council under the European Union’s Horizon 2020 Research and Innovation Program Grant 802414. We thank The Bahamas Environment, Science & Technology Commission; the Gerace Research Center on San Salvador Island; and many Bahamian people for support of our field work. This work benefited from many discussions with the PLIOMAX community. We acknowledge PALSEA (Paleo constraints on sea level rise), a working group of the International Union for Quaternary Sciences and Past Global Changes, which in turn, received support from the Swiss Academy of Sciences and the Chinese Academy of Sciences. TanDEM-X digital elevation data are under copyright by the German Aerospace Center. All rights reserved; used with permission within Project DEM_GEOL1210 (A.R./M.E.R.).

Supporting Information

Appendix (PDF)
Dataset_S01 (CSV)


R. E. Kopp, F. J. Simons, J. X. Mitrovica, A. C. Maloof, M. Oppenheimer, Probabilistic assessment of sea level during the last interglacial stage. Nature 462, 863–869 (2009).
A. Dutton, K. Lambeck, Ice volume and sea level during the last interglacial. Science 337, 216–219 (2012).
A. Dutton et al., Sea-level rise due to polar ice-sheet mass loss during past warm periods. Science 349, aaa4019 (2015).
N. P. McKay, J. T. Overpeck, B. L. Otto-Bliesner, The role of ocean thermal expansion in Last Interglacial sea level rise. Geophys. Res. Lett. 38, L14605 (2011).
A. Levermann et al., The multimillennial sea-level commitment of global warming. Proc. Natl. Acad. Sci. U.S.A. 110, 13745–13750 (2013).
R. M. DeConto, D. Pollard, Contribution of Antarctica to past and future sea-level rise. Nature 531, 591–597 (2016).
K. M. Cuffey, S. J. Marshall, Substantial contribution to sea-level rise during the last interglacial from the Greenland ice sheet. Nature 404, 591–594 (2000).
W. Farrell, J. A. Clark, On postglacial sea level. Geophys. J. Int. 46, 647–667 (1976).
W. Peltier, J. Andrews, Glacial-isostatic adjustment. I. The forward problem. Geophys. J. Int. 46, 605–646 (1976).
K. Lambeck, M. Nakada, Constraints on the age and duration of the last interglacial period and on sea-level variations. Nature 357, 125–128 (1992).
C. P. Conrad, The solid Earth’s influence on sea level. Geol. Soc. Am. Bull. 125, 1027–1052 (2013).
W. R. Peltier, Ice age paleotopography. Science 265, 195–201 (1994).
N. Haskell, The motion of a viscous fluid under a surface load. Physics 6, 265–269 (1935).
C. L. Pekeris, Thermal convection in the interior of the Earth. Geophys. J. Int. 3, 343–367 (1935).
A. M. Forte, W. R. Peltier, The kinematics and dynamics of poloidal–toroidal coupling in mantle flow: The importance of surface plates and lateral viscosity variations. Adv. Geophys. 36, 1–119 (1994).
R. Moucha et al., Dynamic topography and long-term sea-level variations: There is no such thing as a stable continental platform. Earth Planet Sci. Lett. 271, 101–108 (2008).
M. E. Raymo, J. X. Mitrovica, M. J. O’Leary, R. M. DeConto, P. J. Hearty, Departures from eustasy in Pliocene sea-level records. Nat. Geosci. 4, 328–332 (2011).
J. Austermann, J. X. Mitrovica, P. Huybers, A. Rovere, Detection of a dynamic topography signal in last interglacial sea-level records. Sci. Adv. 3, e1700457 (2017).
S. Dendy, J. Austermann, J. Creveling, J. Mitrovica, Sensitivity of last interglacial sea-level high stands to ice sheet configuration during Marine Isotope Stage 6. Quat. Sci. Rev. 171, 234–244 (2017).
E. J. Rohling et al., Differences between the last two glacial maxima and implications for ice-sheet, δ 18O, and sea-level reconstructions. Quat. Sci. Rev. 176, 1–28 (2017).
J. X. Mitrovica, Haskell [1935] revisited. J. Geophys. Res. Solid Earth 101, 555–569 (1996).
H. C. Lau et al., Inferences of mantle viscosity based on ice age data sets: Radial structure. J. Geophys. Res. Solid Earth 121, 6991–7012 (2016).
A. Rovere et al., The Mid-Pliocene sea-level conundrum: Glacial isostasy, eustasy and dynamic topography. Earth Planet Sci. Lett. 387, 27–33 (2014).
D. F. McNeill, Accumulation rates from well-dated late Neogene carbonate platforms and margins. Sediment. Geol. 175, 73–87 (2005).
P. Garrett, S. J. Gould, Geology of New Providence Island, Bahamas. Geol. Soc. Am. Bull. 95, 209–220 (1984).
B. White, K. A. Kurkjy, H. A. Curran, K. A. Besom, Shallowing-upward sequence in a pleistocene coral reef and associated facies, San Salvador, Bahamas. Am. Assoc. Pet. Geol. Bull. 68, 539–539 (1984).
J. Chen, H. Curran, B. White, G. Wasserburg, Precise chronology of the last interglacial period: 234U-230Th data from fossil coral reefs in The Bahamas. Geol. Soc. Am. Bull. 103, 82–97 (1991).
J. L. Carew, J. E. Mylroie, Quaternary tectonic stability of the B/ahamian Archipelago: Evidence from fossil coral reefs and flank margin caves. Quat. Sci. Rev. 14, 145–153 (1995).
W. G. Thompson, H. A. Curran, M. A. Wilson, B. White, Sea-level oscillations during the last interglacial highstand recorded by Bahamas corals. Nat. Geosci. 4, 684–687 (2011).
L. V. Illing, Bahaman calcareous sands. AAPG Bull. 38, 1–95 (1954).
M. Ball, Carbonate sand bodies of Florida and The Bahamas. J. Sediment. Res. 37, 556–591 (1967).
P. J. Hearty, A. C. Neumann, D. S. Kaufman, Chevron ridges and runup deposits in The Bahamas from storms late in oxygen-isotope substage 5e. Quaternary Research 50, 309–322 (1998).
P. Kindler, A. Strasser, Palaeoclimatic significance of co-occurring wind-and water-induced sedimentary structures in the last-interglacial coastal deposits from Bermuda and The Bahamas. Sediment. Geol. 131, 1–7 (2000).
J. Bourgeois, R. Weiss, “Chevrons” are not mega-tsunami deposits: A sedimentologic assessment. Geology 37, 403–406 (2009).
L. Vimpere, P. Kindler, S. Castelltort, Chevrons: Origin and relevance for the reconstruction of past wind regimes. Earth Sci. Rev. 193, 317–332 (2019).
D. E. Hattin, V. L. Warren, Stratigraphic analysis of a fossil Neogoniolithon-capped patch reef and associated facies, San Salvador, Bahamas. Coral Reefs 8, 19–30 (1989).
C. Kerans, C. Zahm, S. L. Bachtel, P. Hearty, H. Cheng, Anatomy of a late Quaternary carbonate island: Constraints on timing and magnitude of sea-level fluctuations, West Caicos, Turks and Caicos Islands, BWI. Quat. Sci. Rev. 205, 193–223 (2019).
T. Tamura, Beach ridges and prograded beach deposits as palaeoenvironment records. Earth Sci. Rev. 114, 279–297 (2012).
A. Strasser, E. Davaud, Formation of Holocene limestone sequences by progradation, cementation, and erosion; two examples from The Bahamas. J. Sediment. Res. 56, 422–428 (1986).
E. K. Potter, K. Lambeck, Reconciliation of sea-level observations in the Western North Atlantic during the last glacial cycle. Earth Planet Sci. Lett. 217, 171–181 (2004).
G. A. Milne, M. Peros, Data–model comparison of Holocene sea-level change in the circum-Caribbean region. Global Planet. Change 107, 119–131 (2013).
R. Love et al., The contribution of glacial isostatic adjustment to projections of sea-level change along the Atlantic and Gulf coasts of North America. Earth’s Future 4, 440–464 (2016).
N. S. Khan et al., Drivers of Holocene sea-level change in the Caribbean. Quat. Sci. Rev. 155, 13–36 (2017).
D. R. Muhs, K. R. Simmons, R. R. Schumann, E. S. Schweig, M. P. Rowe, Testing glacial isostatic adjustment models of last-interglacial sea level history in The Bahamas and Bermuda. Quat. Sci. Rev. 233, 106212 (2020).
A. C. Neumann, P. J. Hearty, Rapid sea-level changes at the close of the last interglacial (substage 5e) recorded in Bahamian island geology. Geology 24, 775–778 (1996).
P. J. Hearty, A. C. Neumann, Rapid sea level and climate change at the close of the Last Interglaciation (MIS 5e): Evidence from the Bahama Islands. Quat. Sci. Rev. 20, 1881–1895 (2001).
G. W. Lynts, Conceptual model of the Bahamian Platform for the last 135 million years. Nature 225, 1226–1228 (1970).
W. Peltier, D. Argus, R. Drummond, Space geodesy constrains ice age terminal deglaciation: The global ICE-6G_C (VM5a) model. J. Geophys. Res. Solid Earth 120, 450–487 (2015).
K. Lambeck et al., Constraints on the Late Saalian to early Middle Weichselian ice sheet of Eurasia from field data and rebound modelling. Boreas 35, 539–575 (2006).
F. Colleoni, C. Wekerle, J. O. Näslund, J. Brandefelt, S. Masina, Constraint on the penultimate glacial maximum Northern Hemisphere ice topography (140 kyrs BP). Quat. Sci. Rev. 137, 97–112 (2016).
C. L. Batchelor et al., The configuration of Northern Hemisphere ice sheets through the Quaternary. Nat. Commun. 10, 3713 (2019).
N. L. M. Barlow et al., Lack of evidence for a substantial sea-level fluctuation within the Last Interglacial. Nat. Geosci. 11, 627–634 (2018).
B. White, H. A. Curran, M. A. Wilson, Bahamian coral reefs yield evidence of a brief sea-level lowstand during the last interglacial. Carbonates Evaporites 13, 10–22 (1998).
M. A. Wilson, H. A. Curran, B. White, Paleontological evidence of a brief global sea-level event during the last interglacial. Lethaia 31, 241–250 (1998).
A. Skrivanek, J. Li, A. Dutton, Relative sea-level change during the Last Interglacial as recorded in Bahamian fossil reefs. Quat. Sci. Rev. 200, 160–177 (2018).
F. D. Hibbert et al., Coral indicators of past sea-level change: A global repository of U-series dated benchmarks. Quat. Sci. Rev. 145, 1–56 (2016).
T. Lorscheid, A. Rovere, The indicative meaning calculator – quantification of paleo sea-level relationships by using global wave and tide datasets. Open Geospatial Data, Software Stand. 4, 10 (2019).
C. Hay et al., The sea-level fingerprints of ice-sheet collapse during interglacial periods. Quat. Sci. Rev. 87, 60–69 (2014).
J. Jouzel et al., Orbital and millennial Antarctic climate variability over the past 800,000 years. Science 317, 793–796 (2007).
D. Dahl-Jensen et al., Eemian interglacial reconstructed from a Greenland folded ice core. Nature 493, 489–494 (2013).
M. Helsen, W. Van De Berg, R. Van De Wal, M. Van Den Broeke, J. Oerlemans, Coupled regional climate–ice-sheet simulation shows limited Greenland ice loss during the Eemian. Clim. Past 9, 1773–1788 (2013).
A. M. Yau, M. L. Bender, A. Robinson, E. J. Brook, Reconstructing the last interglacial at Summit, Greenland: Insights from GISP2. Proc. Natl. Acad. Sci. U.S.A. 113, 9710–9715 (2016).
J. L. Bamber, R. E. Riva, B. L. Vermeersen, A. M. LeBrocq, Reassessment of the potential sea-level rise from a collapse of the West Antarctic Ice Sheet. Science 324, 901–903 (2009).
D. Farinotti et al., A consensus estimate for the ice thickness distribution of all glaciers on Earth. Nat. Geosci. 12, 168–173 (2019).
V. J. Polyak et al., A highly resolved record of relative sea level in the western Mediterranean Sea during the last interglacial period. Nat. Geosci. 11, 860–864 (2018).
P. U. Clark et al., Oceanic forcing of penultimate deglacial and last interglacial sea-level rise. Nature 577, 660–664 (2020).
M. Aurell, D. F. McNeill, T. Guyomard, P. Kindler, Pleistocene shallowing-upward sequences in New Providence, Bahamas; signature of high-frequency sea-level fluctuations in shallow carbonate platforms. J. Sediment. Res. 65, 170–182 (1995).
B. White, H. A. Curran, “Coral reef to Eolianite transition in the Pleistocene rocks of Great Inagua Island, Bahamas” in Proceedings of the Third Symposium on the Geology of The Bahamas, H. A. Curran, Ed. (CCFL Bahamian Field Station, Ft. Lauderdale, FL, 1987), pp. 165–171.
P. Caldwell, M. Merrifield, P. Thompson, Sea level measured by tide gauges from global oceans—the Joint Archive for Sea Level holdings (NCEI Accession 0019568) Dataset 10, V5V40S7W (Version 5.5, NOAA National Centers for Environmental Information, Silver Spring, MD, 2015).
B. C. Douglas, “Sea level change in the era of the recording tide gauge” in Sea Level Rise, B. Douglas, M. Kearney, S. Leatherman, Eds. (International Geophysics, v75, Academic Press, 2001), pp. 37–64.
X. Glorot, Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks” in Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics (AISTATS, 2010), vol. 9, pp. 249–256.
D. P. Kingma, J. Ba, Adam: A method for stochastic optimization. arXiv [Preprint] (2014). (Accessed 20 January 2018).
R. A. Kendall, J. X. Mitrovica, G. A. Milne, On post-glacial sea level. II. Numerical formulation and comparative results on spherically symmetric models. Geophys. J. Int. 161, 679–706 (2005).
A. M. Dziewonski, D. L. Anderson, Preliminary reference Earth model. Phys. Earth Planet. In. 25, 297–356 (1981).
C. Waelbroeck et al., Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records. Quat. Sci. Rev. 21, 295–305 (2002).
P. Blanchon, A. Eisenhauer, J. Fietzke, V. Liebetrau, Rapid sea-level rise and reef back-stepping at the close of the last interglacial highstand. Nature 458, 881–884 (2009).
A. L. Thomas, et al., Penultimate deglacial sea-level timing from uranium/thorium dating of Tahitian corals. Science 324, 1186–1189 (2009).
M. Siddall, E. Bard, E. J. Rohling, C. Hemleben, Sea-level reversal during termination II. Geology 34, 817–820 (2006).
A. J. Schaeffer, S. Lebedev, Global shear speed structure of the upper mantle and transition zone. Geophys. J. Int. 194, 417–449 (2013).
A. Schaeffer, S. Lebedev, Imaging the North American continent using waveform inversion of global and USArray data. Earth Planet Sci. Lett. 402, 26–41 (2014).
S. W. French, B. A. Romanowicz, Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography. Geophys. J. Int. 199, 1303–1327 (2014).
H. Yamauchi, Y. Takei, Polycrystal anelasticity at near-solidus temperatures. J. Geophys. Res. Solid Earth 121, 7790–7820 (2016).
F. D. Richards, M. J. Hoggard, N. White, S. Ghelichkhan, Quantifying the relationship between short-wavelength dynamic topography and thermomechanical structure of the upper mantle using calibrated parameterization of anelasticity. J. Geophys. Res. Solid Earth 125, e2019JB019062 (2020).
K. Latychev et al., Glacial isostatic adjustment on 3-D Earth models: A finite-volume formulation. Geophys. J. Int. 161, 421–444 (2005).
J. Austermann, M. J. Hoggard, K. Latychev, J. X. Mitrovica, F. D. Richards, The effect of lateral variations in Earth structure on last interglacial sea level. Earth ArXiv [Preprint] (2021). (Accessed 2 July 2021).
M. D. Hoffman, A. Gelman, The No-U-Turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res. 15, 1593–1623 (2014).
J. Salvatier, T. V. Wiecki, C. Fonnesbeck, Probabilistic programming in Python using PyMC3. PeerJ Comput. Sci. 2, e55 (2016).
A. Vehtari, A. Gelman, J. Gabry, Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat. Comput. 27, 1413–1432 (2017).
Y. Yao, A. Vehtari, D. Simpson, A. Gelman, Using stacking to average Bayesian predictive distributions (with discussion). Bayesian Anal. 13, 917–1007 (2018).
R. Kumar, C. Carroll, A. Hartikainen, O. A. Martín, ArviZ a unified library for exploratory analysis of Bayesian models in Python. J. Open Source Software 4, 1143 (2019).

Information & Authors


Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 118 | No. 33
August 17, 2021
PubMed: 34373328


Data Availability

All data, models, and code used in this work are available in GitHub ( and Zenodo (

Submission history

Published online: August 9, 2021
Published in issue: August 17, 2021


  1. last interglacial sea level
  2. glacial isostatic adjustment
  3. The Bahamas


This material is based upon work supported by NSF Grants OCE-1202632 “PLIOMAX - Pliocene maximum sea level” and OCE-1841888. J.A. and M.E.R. acknowledge support from the Vetlesen Foundation. The research of A.R. on the LIG was funded by Institutional Strategy of the University of Bremen, German Excellence Initiative Grant ABPZuK-03/2014 and the European Research Council under the European Union’s Horizon 2020 Research and Innovation Program Grant 802414. We thank The Bahamas Environment, Science & Technology Commission; the Gerace Research Center on San Salvador Island; and many Bahamian people for support of our field work. This work benefited from many discussions with the PLIOMAX community. We acknowledge PALSEA (Paleo constraints on sea level rise), a working group of the International Union for Quaternary Sciences and Past Global Changes, which in turn, received support from the Swiss Academy of Sciences and the Chinese Academy of Sciences. TanDEM-X digital elevation data are under copyright by the German Aerospace Center. All rights reserved; used with permission within Project DEM_GEOL1210 (A.R./M.E.R.).



Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964;
School of Earth and Ocean Sciences, University of Victoria, Victoria, BC V8W 2Y2, Canada;
Jacqueline Austermann
Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964;
William J. D’Andrea
Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964;
Roger C. Creel
Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964;
Michael R. Sandstrom
Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964;
Miranda Cashman
Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964;
Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964;
Center for Marine Environmental Sciences (MARUM), University of Bremen, D-28359 Bremen, Germany
Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964;


To whom correspondence may be addressed. Email: [email protected] or [email protected].
Author contributions: B.D., J.A., W.J.D., A.R., and M.E.R. designed research; B.D., J.A., W.J.D., R.C.C., M.R.S., M.C., A.R., and M.E.R. performed research; B.D. contributed new reagents/analytic tools; B.D. and J.A. analyzed data; and B.D., J.A., W.J.D., R.C.C., M.R.S., M.C., A.R., and M.E.R. wrote the paper.
Reviewers: F.D.H., University of York; and G.A.M., University of Ottawa.

Competing Interests

The authors declare no competing interest.

Metrics & Citations


Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Citation statements



If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by


    View Options

    View options

    PDF format

    Download this article as a PDF file


    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to get full access to it.

    Single Article Purchase

    Sea-level trends across The Bahamas constrain peak last interglacial ice melt
    Proceedings of the National Academy of Sciences
    • Vol. 118
    • No. 33







    Share article link

    Share on social media