## Abstract

Estimating minimum ice volume during the last interglacial based on local sea-level indicators requires that these indicators are corrected for processes that alter local sea level relative to the global average. Although glacial isostatic adjustment is generally accounted for, global scale dynamic changes in topography driven by convective mantle flow are generally not considered. We use numerical models of mantle flow to quantify vertical deflections caused by dynamic topography and compare predictions at passive margins to a globally distributed set of last interglacial sea-level markers. The deflections predicted as a result of dynamic topography are significantly correlated with marker elevations (>95% probability) and are consistent with construction and preservation attributes across marker types. We conclude that a dynamic topography signal is present in the elevation of last interglacial sea-level records and that the signal must be accounted for in any effort to determine peak global mean sea level during the last interglacial to within an accuracy of several meters.

## INTRODUCTION

To gain insight into ice sheet stability in the face of global warming, scientists have turned to studies of the geologic past and, in particular, periods during which temperatures may have been elevated relative to the present day (*1*). The last interglacial (LIG), or marine isotope stage (MIS) 5e, ~125 thousand years ago (ka), is of particular interest in this regard. Although CO_{2} concentrations during the LIG were comparable to preindustrial values (*2*), global average temperatures appear to have been 1° to 2°C above modern values (*3*).

On the basis of a probabilistic assessment of widely distributed sea-level and isotopic records from the LIG, Kopp *et al*. (*4*) reported that global mean sea level (GMSL) during the LIG very likely peaked at least 6.6 m (95% probability) higher than that today but was unlikely (33% probability) to have exceeded 9.4 m. Analyses of paleo sea level indicators in a small number of tectonically stable sites (most notably the Seychelles and Western Australia) also indicate that peak LIG sea level lies between 6 and 9 m above present levels (*5*–*7*).

A complication in all these studies is that various geodynamic processes contribute to the present elevation of paleo sea level records (*8*). A notable example of these processes is tectonic uplift or subsidence at active plate boundaries [for example, see Zazo *et al*. (*9*)], which often leads to the exclusion of these sites in reconstructions of past GMSL. Another important deformational process is the response of the Earth system to changes in ice and ocean loading during ice age cycles (*10*, *11*), or glacial isostatic adjustment (GIA), a process first studied in the context of the LIG by Lambeck and Nakada (*12*). The accuracy of model-derived corrections for this global process is subject to uncertainties in ice history and mantle viscoelastic structure [for example, see Lambeck *et al*. (*13*)]. Additionally, the redistribution of sediment can lead to sea level changes through the buildup of topography and loading-induced deformation of the solid Earth and gravity field (*14*, *15*).

Earth’s surface is further deflected by viscous stresses associated with buoyancy variations and flow within Earth’s convective mantle that can alter the elevation of sea-level markers subsequent to their formation (*16*–*22*). Effects of this so-called dynamic topography (DT) on the current elevation of sea-level markers dating to the mid-Pliocene (~3 million years ago) have been documented (*23*–*25*) and imply meter-scale displacements for LIG sea-level markers (*24*). Kopp *et al*. (*4*) incorporated uncertainties due to vertical land movement and applied nonzero rates in several passive margin regions. Although this correction may implicitly include the DT process, effects of DT are generally not addressed in sea level studies of Pleistocene interglacials and have not previously been shown to be detectable during the LIG. Here, we quantify and analyze the effects of DT on globally distributed markers of local peak sea level during the LIG.

## RESULTS AND DISCUSSION

We use a recent compilation of LIG sea-level indicators (*26*) drawn from published databases [(*4*, *27*–*29*); see Materials and Methods for details] and distinguish four groups based on their tectonic setting: passive margins, active margins, oceanic islands characterized by recent volcanic activity, and oceanic islands of nonvolcanic origin or characterized by extinct volcanism (Fig. 1A). Note that many age determinations for these markers, especially for sites in extratropical regions, are based on either relative (for example, amino acid racemization or biostratigraphy) or analytic (for example, optically stimulated luminescence and electron spin resonance dating) methods that are less precise than uranium-series ages. We consider data to be from the same location if they are within ~20 km of one another, and we retain only the highest elevation when multiple observations are available in one location. Data from active margin sites have likely been perturbed by processes that are not modeled in this study, such as plate boundary deformation and hot spot dynamics, and these sites are therefore not considered further. The final database consists of 298 passive margin sites that can be divided into four different landform categories (*30*): marine terraces, fossil coral reefs, beach deposits, and beach ridges (see Fig. 1B for sea-level marker locations and fig. S1 for field examples).

The sample elevation distribution associated with sea-level indicators differs markedly according to category (Fig. 2, A to E). Sample distributions are estimated accounting for the reported measurement uncertainty at each site by randomly sampling the data’s respective uncertainty range. For the data taken from Hibbert *et al*. (*27*), we assign uncertainty ranges based on parameterized modern depth ranges of the respective coral genera (see Materials and Methods and fig. S2), and, for the remaining data, we assume that uncertainties follow a normal distribution. An offset between the elevation of a sea-level marker and the actual sea level including the uncertainty thereof is referred to as the indicative meaning and must be included for each data point. We use the paleo sea level elevations including the assumed indicative meaning, as reported in previous work (*4*, *26*–*29*). To account for spatial clustering in the data set when estimating global sea-level marker distributions, we sample each site with a frequency that is inversely proportional to its covariance with surrounding data (see Materials and Methods and fig. S3).

Although 281 of the 298 sea-level markers are represented as having normal distributions, the combined distribution (Fig. 2A) is positively skewed because of differences in the mean and standard deviation (SD) associated with individual markers. This skew is reflected in a mean value of 7 m, which is higher than the maximum likelihood value of 5.4 m, and is also present if, instead of retaining the highest elevation when multiple observations are available in one location, the mean elevation for this location is used (see fig. S4).

The positive skew in global sea-level markers may reflect greater rates of preservation and identification of LIG markers that are emplaced, or subsequently elevated, further above sea level. Figure 2 (B to E) shows the distributions associated with the four main types of sea-level indicators that together account for 96% of the observations (marine terraces, fossil coral reefs, beach deposits, and beach ridges). The thin negative tail in the distribution of the full data set primarily represents drowned coral reef records (Fig. 2C), whereas occurrences in the stronger upper tail between 10 and 20 m are due to beach ridge indicators (Fig. 2E) and marine terraces (Fig. 2B). This offset could be associated with systematic differences in deformation (uplift or subsidence) between these indicator types or due to an inaccurate assessment of the indicative meaning in the data interpretation. For example, coral reefs generally constitute a lower bound on sea level, whereas beach ridges most often form 1 to 6 m above local sea level, depending on local wave conditions. A reevaluation of differences in the indicative meaning of various landforms appears warranted and will have implications for estimates of peak GMSL during the LIG if these data sets are used.

Although we later use numerical models to correct for DT, one component in our definition of DT—thermal subsidence of oceanic lithosphere as it moves away from mid-ocean ridges and cools—can be estimated based on observed bathymetry and oceanic lithosphere age alone [fig. S5; (*31*, *32*)]. Although not usually accounted for explicitly in studies of LIG sea level (*4*–*7*, *33*), this thermal subsidence has been shown to lower the elevation of certain sea-level markers in the Pacific [for example, see Dickinson (*34*)]. In our database, we detect a significant correlation between the thermal subsidence of the lithosphere and coral reef sea-level indicators that are located within ocean basins (fig. S6M). However, the correlation between the former and the remaining indicator types is low (fig. S6, L, N, and O).

To account for the full DT signal, including sublithospheric buoyancy variations, we perform simulations using the mantle convection code ASPECT that are based on four different Earth models paired with three different top boundary conditions [(*35*, *36*); see Materials and Methods and fig. S7 for modeling details]. The simulations are constrained to fit present-day geophysical observables, including long-wavelength gravity, heat flow, long-wavelength present-day DT, plate motions, and core-mantle boundary topography (*37*, *38*). Simulations use buoyancy variations within the entire mantle and therefore include the signal associated with ocean lithosphere subsidence, which is particularly apparent around intermediate- to fast-spreading ridges (for example, the Southeast Indian Ridge and the East Pacific Ridge; Fig. 3A). Thus, we do not apply a separate ocean-cooling correction.

Our DT simulations indicate that changes of several meters have occurred since the LIG (Fig. 3A), a magnitude that is consistent with inferences from observations [for example, see Rovere *et al*. (*24*), Richards *et al*. (*39*), and Czarnota *et al*. (*40*)]. The SD at LIG data sites across the 12 simulations averages to 2.6 m (Fig. 3B) and reflects uncertainties in Earth’s viscoelastic and buoyancy structure, as well as the coupling between mantle flow and surface plates.

Despite the uncertainties inherent to our DT simulations, the average of our 12 DT simulations captures notable geographic variability in the observations (Fig. 3C). Markers from the Maldives and other islands south of India show large negative LIG highstand elevations of approximately −20 m that are in line with (albeit larger than) predictions of −10 m by the DT models, which are driven partly by thermal subsidence. A 10-m tilt across southeastern Australia and Tasmania is also captured by the DT predictions, as is some of the ~20-m uplift in East Antarctica and 3-m uplift in Siberia. Observed highstand elevations in Angola are also directionally consistent with DT predictions, although underestimated. None of our DT simulations capture the uplift inferred for the South Atlantic islands. We also note that predicted changes due to DT at the Seychelles are −0.5 ± 2 m (1σ), and, in Western Australia, they range from −2 ± 2 m in the north to 4 ± 7 m (1σ) in the south. DT may thus play a role in the discrepant elevations of LIG markers noted by Dutton and Lambeck (*5*) between these locations. There are large variations in highstand elevations at sites within previously glaciated regions, such as Scandinavia and the Gulf of Saint Lawrence (Fig. 3C). This variability is expected on account of deformation associated with GIA, which is sufficiently large to obscure contributions from DT in this region (see Materials and Methods for a discussion of the role of GIA in deforming the data locations used in this study).

To quantify global goodness of fit, we examine the overall correlation between our DT predictions and the observations. The 12 DT simulations produce spatial fields that correlate with observed marker elevations with an average correlation coefficient of 0.24. If only the four highest correlating DT simulations are retained, the correlation is 0.32 (Fig. 2F). Accounting for regional coherence, we find that the magnitudes of these correlations are significant at the 95% level (see Materials and Methods). The correlation between predicted DT in these four simulations and individual categories of indicator types is also significant at the 90% level for marine terraces, coral reefs, and beach deposits (Fig. 2, G to I). Beach ridges are not correlated significantly given the much higher required *R* value, which is due to the small sample size (21 data points) and limited spatial extent of this indicator type (see Fig. 1B). These results constitute the first detection of a DT signal in LIG sea-level highstand records.

As noted for individual marker sites, the magnitude of predicted DT changes tends to be lower than that of observed changes, and a linear trend fit between observations and predictions has a slope of only 0.19 m/m. The fact that the slope is not closer to one may reflect an incorrect interpretation of the indicative meaning of some sea-level indicators or missing corrections, such as those required for GIA or sediment redistribution and loading. The magnitude of the DT prediction may also be underestimated by the numerical simulations given the long-wavelength structure inherent to the Earth models used in this study (*41*), our incomplete knowledge of asthenospheric viscosity, and lateral density variations that are potentially underestimated due to the low resolving power of amplitudes in seismic tomography models.

Correcting the sea-level indicator elevations for the effect of DT (using the four highest correlating simulations) results in a global marker distribution that has lower variance and is less skewed than the uncorrected distribution, implying greater accuracy and a smaller difference between the mean and maximum likelihood value (Fig. 2F). This tightening of the overall sea-level marker distribution helps to resolve inconsistencies in marker elevations, where the range of the distribution means according to marker types reduces from 4.7 to 13.7 m (Fig. 2, B to E) to 4.5 to 11.1 m (Fig. 2, G to J).

The appearance of systematic corrections to the mean associated with marker type may reflect a preservation bias in the observation. Uplifted beach ridges are more likely to be preserved because they are otherwise easily eroded by wave activity during storms. Albeit somewhat more resilient to storms, similar preservational concerns hold for consolidated beach deposits. Consistent with uplift improving the odds of preserving these deposits in the geologic record, our DT estimates indicate 2.6 m of uplift, on average, for beach ridges and 1.3 m of uplift for beach deposits. In contrast, fossil coral reefs are more robust and are, based on our predictions, almost equally likely to be found in uplifting and subsiding locations, having a mean and SD DT across marker locations of −0.6 ± 5.0 m (1σ). Marine terraces are similarly robust and are estimated to experience 0.1 ± 3.4 m (1σ) of deformation due to DT. These results emphasize the importance of considering DT not only for applying corrections but also for determining patterns of preservation.

## CONCLUSION

Surface motions driven by mantle convective flow have the potential to produce important vertical deflections that alter the elevations of Pleistocene sea-level markers. Numerical models of mantle convection predict changes in DT since 125 ka that are detected at a high level of statistical significance in a database of LIG highstands. These changes affect the elevation of local sea-level markers that are used to estimate peak GMSL during the LIG and can influence these estimates by several meters, depending on the data used. To illustrate this uncertainty, note that 2 m of GMSL rise is equivalent to more than half of the ice volume expected in a collapse of the West Antarctic Ice Sheet. Predicted rates of DT change are relatively constant across the million-year numerical mantle convection simulations, indicating that these corrections increase by a factor of ~3 to 4 for the earlier MIS 11 interglacial. Improvement in the accuracy of Pleistocene sea-level reconstruction thus substantially depends on progress in modeling DT, including reducing current uncertainties in mantle buoyancy and viscosity, and in the coupling between mantle flow and surface plates.

## MATERIALS AND METHODS

### Database of relative sea-level indicators

Here, we used an updated version of the LIG sea-level database compiled by Rovere *et al*. (*26*) (available online at http://pliomax.weebly.com/data--products.html). Our database combines observations provided by Kopp *et al*. (*4*), Hibbert *et al*. (*27*), Pedoja *et al*. (*28*), and Ferranti *et al*. (*29*). The sea-level indicators described by each author were divided into one of the landform categories identified by Rovere *et al*. (*30*). We maintained the original paleo relative sea-level elevation reported by each study, as discussed below. Note that we excluded any tectonic or vertical land movement correction that might have been applied in the original publication.

We extracted data at 71 sites worldwide from Kopp *et al*. (*4*). In their database, there is a clear distinction between marker elevation (and associated accuracy), reference water level, and indicative range (*30*). The authors reconstructed these values using data contained in the original publications, as well as generic descriptions of similar landforms. For the paleo relative sea-level values and related uncertainty, we used the original columns “PaleoSL” and “PaleoSLSD.”

Hibbert *et al*. (*27*) compiled a database of uranium-series dated corals in which each single coral measurement was treated as a sea-level index point, and paleo sea level was calculated following the modern depth distribution of coral genera. The database contains 2496 index points. We removed index points that (i) fall outside MIS 5e (130 to 115 ka), (ii) have an elevation marked as “unknown” in the original database, (iii) have no age assignment, and (iv) do not specify the species or genus. This culling yielded 494 data points, 319 of which are located at 26 different locations on passive margins. For genera that have a present living depth range that is well approximated by a normal distribution, we used the median depth listed in the database as an approximate value for paleo water depth and assumed a normal distribution as the indicative range, with an SD calculated as the mean of the lower and upper 68% error (fig. S2A). This includes genera such as *Montastraea* and *Siderastrea*. For genera that are located close to the sea surface and have a living range with an asymmetric distribution, we used a gamma function to parameterize the indicative meaning (fig. S2B). The gamma function is fit to the observed modern depth distribution provided by Hibbert *et al*. (*27*). This includes genera such as *Acropora*, *Porites*, *Goniastrea*, and *Faviidae*.

Pedoja *et al*. (*28*) reviewed records from a large number (942) of sites worldwide. The database contains a significantly larger number of data than Kopp *et al*. (*4*), and its structure is more complicated. First, in regard to the geographic location of samples, the authors derived mean elevations within primary tectonic zones with more abundant data. As an example, Ferranti *et al*. (*29*) reported 246 sites in Italy that were combined in the database of Pedoja *et al*. (*28*) into 71 sites. For this reason, the database of Pedoja *et al*. (*28*) includes four columns for latitude/longitude coordinates to demarcate zones within which they averaged observations. Here, we averaged these coordinates and associated them with a single point. The paleo relative sea-level estimates were reported alongside their uncertainties, but the indicative meaning of each landform was not explicitly considered. In regard to the uncertainty in the sea-level elevation, the authors state, “in this compilation, we systematically attributed a minimum error range of one meter to the measurements on elevation published without any margin of error. Where authors provided altitude ranges, we took the mean value of the range; for example, with a shoreline angle between 20 and 30 m, Error = 25 +/− 5 m.” From this description, it is apparent that the error indicated for each elevation corresponds to a measurement error and does not include a consideration of the indicative range. Therefore, the uncertainty of paleo sea level could be underestimated in this study. For the paleo relative sea-level values and related uncertainty, we used the columns “max elevation (m)” and “MoE (m).”

Ferranti *et al*. (*29*) reviewed MIS 5e sea-level markers in Italy; the 14 contributing authors are all geologists active in different regions of the Italian peninsula. We reported this database as, in contrary to the other ones we compiled, each marker has been reviewed or remeasured and inserted in the database by a regional expert, coordinated by one author that visited most of the outcrops. For the paleo relative sea-level values, we used the columns termed “shoreline elevation.” The data column named “Uncertainty” represents an estimate of measurement error plus indicative range uncertainties, as discussed by Ferranti *et al*. (*29*). Regarding the measurement errors, the authors state, “At most sites we performed direct measurement of the shoreline height using a tape referenced to the mean sea- level, and in few cases geodetical stations were used.” The authors did not quantify this error, which could be significant (1 to 2 m, depending on survey conditions). We note that all of the data by Ferranti *et al*. (*29*) are located at active margins and therefore not included in the final data set used for the analysis in this study.

A compilation of all data and model results used in this study can be found in the Supplementary Materials. For detailed information on each site, we refer the reader to the original publications and to the source databases. The studies on which each source database is based are listed alphabetically in an additional spreadsheet attached to the database.

### DT modeling

To calculate DT, we solved the continuity, momentum, and energy equations for mantle convection, assuming a simple equation of state that relates changes in density to changes in temperature through thermal expansion [for example, see Schubert *et al*. (*42*)]. The convection code ASPECT solves the relevant equations for compressible flow within Earth’s mantle. We used ASPECT 1.4.0-pre (*35*, *36*, *43*) published under the GPL2 license. As input, we adopted four different Earth models: The first two models are based on the density model TX2008, derived to match (in conjunction with a viscosity model) present-day geophysical observables, such as plate velocities, long-wavelength present-day DT, long-wavelength geoid, and the excess ellipticity of the core mantle boundary (*37*). These two Earth models differ by the choice of the radial viscosity profiles, V1 and V2 (fig. S7). These have been inferred from simultaneous inversions of these convection-related data sets plus a suite of observations associated with glacial isostatic adjustment (*44*, *45*). The two additional Earth models that were considered are based on the shear-wave tomography models S40RTS (*46*) and Savani (*47*). Density perturbations in each case were calculated from shear-wave perturbations using the depth-dependent conversion factor given by Steinberger (*38*). The surface density was assumed to be 3300 kg/m^{3}. In addition, we determined the depth and buoyancy of the continental lithosphere using the parameterization by Steinberger (*38*) and a depth-dependent viscosity profile specific to each of these two tomography models (fig. S7). The conversion factors and viscosity profiles were chosen to minimize the misfit of predictions to the heat flux, geoid, and present-day DT (*38*). The viscosity profiles were additionally constrained to satisfy the Haskell constraint (*38*, *48*, *49*). Note that Steinberger and Calderwood (*49*) and Steinberger (*38*) only sometimes include the oceanic lithosphere in their reconstructed and modeled DT. In our analysis, we retained the density perturbation within the oceanic lithosphere because it is an important component of the convective system.

Our simulations did not include phase changes, thermal boundary layers, internal radiogenic heat production, or the deflection of internal boundaries. The one-dimensional temperature and density profiles follow an adiabat. Thermal conductivity, thermal expansivity, and heat capacity vary with depth alone and are adopted from Glišović and Forte (*50*). ASPECT did not calculate gravity self-consistently, and we therefore imposed the radially varying gravity profile from Glišović and Forte (*50*). We applied a free-slip boundary condition at the core-mantle boundary. For the top boundary condition, we adopted either free-slip, no-slip, or prescribed present-day plate velocities from Seton *et al*. (*51*). The depth resolution of our numerical grid is ~30 km in the upper 1000 km and varies from 65 to 115 km below this depth. The lateral resolution is on the order of 80 km in the upper 1000 km and 175 to 250 km below that depth. To calculate DT from the radial stress field computed by the mantle convection model, we used the theory described by Austermann and Mitrovica (*52*). The topography calculations did not include self-gravity. The convection simulation was run forward in time for 1.5 million years. To avoid numerical artifacts evident in the first few time steps of the simulation, we computed the rate of change in DT from the mean change in DT between 0.5 and 1.5 million years of the simulation, and this value was used to calculate the change in DT over the past 125 thousand years. To calculate the change in DT at a specific site, we combined perturbations in topography associated with evolving viscous stresses in the mantle with a signal associated with the translation of plates across Earth’s surface. We used the no-net-rotation angular plate velocity tabulation of DeMets *et al*. (*53*) to translate the DT field calculated for the LIG into its present-day coordinates and calculated the difference between the predicted present-day DT field and the rotated predicted LIG DT field. This approach does not account for changes in the ocean basin volume, which could contribute to the peak GMSL estimate for the LIG.

### Details on the statistical analysis: Outliers, clustering, and statistical significance

We removed outliers in the data set to avoid bias in the statistical characteristics, such as the mean and the variance. We identified outliers as sites at which for all the 12 DT corrections, the corrected highstand elevation is further away from the mean than two times the SD (which is above 31 m or below −17 m). This results in 3% (8 of 298) of the sites being regarded as outliers.

The frequency at which we sampled each site was chosen to be a function of how isolated the site is to remove any potential bias introduced by spatial clustering in the data set. We adopted a sampling frequency based on kriging weights, which corrects for spatial clustering in a manner that reflects the structural properties of the DT predictions. We first calculated experimental variograms for each of the 12 DT predictions, where the distance between two points was measured as the length of the great circle that connects them. We then took the mean of the 12 normalized variograms and fit a theoretical variogram, including a nugget effect (see fig. S3A). Next, we used ordinary kriging and the theoretical variogram to calculate the kriging weights. Finally, we obtained the sampling frequency by offsetting and scaling the weights so that the SD of the weights is approximately 20% of the mean. The sample frequency ranged from 7,600 to 18,500, with an average of 10,000 samples per LIG site (see fig. S3B). Correlation coefficients computed in this paper include this weighting for clustering as well as uncertainties in the observed elevation of sea-level indicators.

For a sample size of 298 independent data points, the correlation is significant at a level of 90% (95%) if *R* > 0.075 (0.096). However, this significance level is increased if there is a covariance between model predictions for different sites. The latter is the case, for example, if two sites are located close to one another relative to the spatial wavelength of the model prediction. To account for this covariance, we used rotational sampling of the modeled fields. We generated 1000 random rotations of the 12 original DT fields. We then calculated the synthetic predictions of the rotated models (at the observed locations) and determined outliers as described above. We continued to calculate the correlation between the synthetic prediction and the observed data. The significance levels quoted in this paper are percentiles of the resulting distribution of correlation coefficients. This procedure preserves the spectral and statistical characteristics of the original fields and provides a sensible null hypothesis in testing the statistical significance of correlation. We used an analogous approach to determine statistical significance for the indicator subsets and GIA and ocean subsidence fields.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/7/e1700457/DC1

Supplementary Text

fig. S1. Field examples of four distinct types of sea-level indicators.

fig. S2. Observed and parameterized indicative meaning for corals from the database of Hibbert *et al*. (*27*).

fig. S3. Semivariogram and kriging weights for LIG sites.

fig. S4. Distributions of observed and corrected sea-level estimates when choosing the mean elevation across indicators at the same location (instead of the maximum).

fig. S5. Analysis of ocean basin subsidence due to cooling.

fig. S6. Distributions of observed and corrected sea-level highstands.

fig. S7. Radial viscosity profiles.

fig. S8. Relative performance of different DT models in matching variability of highstand observations.

fig. S9. Variability among different DT predictions.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

**Acknowledgments:**We thank A. Simms and five anonymous reviewers for their constructive comments on this manuscript, as well as R. Moucha for discussions on this topic. We acknowledge PALSEA2 [a PAGES International Geosphere-Biosphere Programme/International Union for Quaternary Science (INQUA) working group] and MOPP (INQUA Project 1603P) for useful discussions at yearly meetings.

**Funding:**Support for this research was provided by Harvard University (to J.A. and J.X.M.), including a Merit Fellowship (to J.A.) from the Graduate School of Arts and Sciences, NSF grant OCE-0825293 “PLIOMAX” (to J.X.M., J.A., and A.R.), and NSF grant AGS-1338832 “VOICE” (to P.H.). We thank the Computational Infrastructure for Geodynamics (http://geodynamics.org), which is funded by the NSF under awards EAR-0949446 and EAR-1550901. J.A. acknowledges funding from the Royal Society. A.R.’s research is supported by the Institutional Strategy of the University of Bremen and funded by the German Excellence Initiative (ABPZuK-03/2014) and ZMT, Leibniz Centre for Tropical Marine Research.

**Author contributions:**J.A. and J.X.M. developed the initial concept. A.R. and J.A. compiled the sea-level database. J.A. conducted the mantle convection modeling. J.X.M. conducted the GIA modeling. J.A. and P.H. performed the statistical analysis. All authors contributed to the analysis and writing of the paper.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional code and input files related to the mantle convection modeling may be requested from the authors.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).