Sea level and global ice volumes from the Last Glacial Maximum to the Holocene

Analysis Strategy

In a zero-order approximation, the change in sea level during glacial cycles is the ice-volume equivalent sea level (esl) Δζ_esl defined asΔζesl(t)=−1ρo∫t1Ao(t)dΔMicedtdt[1]where ΔM_ice is the change in ice mass on the continents and grounded on the shelves at time t with respect to present, A_o is the ocean area defined by the coastline and ice grounding line at t, and ρ_o is the average density of the ocean. Superimposed on this is the response of the earth and ocean surface to the changes in loading. Hence, as a first approximation, the relative sea-level Δζ_rsl at a location ϕ and time t isΔζrsl(ϕ,t)=Δζesl(t)+δζiceiso(ϕ,t)+δζwateriso(ϕ,t)[2]where δζiceiso(ϕ,t) and δζwateriso(ϕ,t) are the isostatic contributions from the changing ice and water loads, respectively. Because of the direct gravitational attraction of the ice sheets and deformation of the ocean basin during the glacial cycle, δζwateriso(ϕ,t) is not independent of ice sheet geometry. Consequently, higher-order iterative solutions are necessary to solve [2].

In the absence of tectonics, the dominant departure of far-field sea level from its global mean [1] is due to δζwateriso(ϕ,t) and results in a strong spatial variability of Δζrsl(ϕ,t) across continental margins and within partially enclosed ocean basins that is dependent on mantle rheology, the amount of water added into the oceans, and the rate at which it is discharged, but is relatively insensitive to where this melt water originated. Typically, during the melting phase, δζwateriso(ϕ,t) is about 15–25% of Δζesl(t) and leads to near-coastal observations of LGM sea level being systematically less than the corresponding Δζesl(t). When deglaciation ceases, the seafloor continues to deform because of the viscous nature of the mantle, and sea levels particularly at continental edges fall slowly until all load stresses have relaxed. This mechanism often produces small highstands at ∼5–7 ka BP that mark the end of the dominant melting period (15, 16, 18).

The δζiceiso(ϕ,t) reflects mantle-scale flow induced by the changes in ice load. This signal includes the evolution of a broad trough and bulge system that extends several thousand kilometers beyond the ice margin, as a combined result of surface deformation and gravitational potential change (the “geoidal” bulge). Beyond these features, δζiceiso(ϕ,t) shows less spatial variability than δζwateriso(ϕ,t) and is relatively insensitive to the details of the ice load, but the amplitudes of these signals are not negligible and are earth-rheology dependent.

The challenge is, in the presence of imperfect and incomplete data, to invert the complete formulation of [2] for both the earth and ice unknowns. This challenge is approached here through the following iterative procedure: (i) Start with far-field sea-level data, assume Δζrsl=Δζesl and using [1] estimate a first approximation of ice mass or volume function ΔVice(t). (ii) Distribute this ice between the known ice sheets, guided by either published information or by glaciological hypotheses and with ice advance and retreat histories, depending on the state of knowledge of the particular ice sheet. (iii) Invert far-field sea-level data for earth rheology parameters E (mantle viscosity, elastic thickness of the lithosphere) and for a corrective term δζesl to the nominal Δζesl0 corresponding to step ii, yielding an improved approximation of ice volume function ΔVice(t). (iv) Analyze sea-level data separately for each of the major ice sheets (near-field analyses) and, starting with the ice models from step ii, invert for both E and corrections to ice-thickness functions I. (v) Impose the condition that the sum of the individual ice volume functions equals ΔVice(t) from step iii and distribute any discrepancies between the ice sheets. (vi) Repeat steps iii–v until convergence has been achieved.

Advantages of this approach include the ability to analyze individual ice sheets with different resolutions depending on the state of a priori knowledge of the ice sheet and quality of the observational data, as well as to make a first-order estimate for lateral variation in mantle viscosity.

Observational Evidence

The principal sources for quantitative sea-level information are from sediment and coral records whose depositional environments relative to mean sea level (MSL) are assumed known and whose ages have been determined either by radiocarbon or uranium-series dating. In the former case, all ¹⁴C ages have been calibrated using either the calibrations provided by the original authors or the IntCal09 calibration (19). Reservoir and isotopic fractionation corrections have been applied where appropriate.

Corals provide lower-limit estimates of mean-low-water-spring (MLWS) tide, but their growth depth range δζ_d is species and environment dependent (20). Assuming that the modern range is representative of the past growth range, the adopted mean sea level from in situ fossil corals of age t at elevation ζ (all vertical measures are positive upward) isΔζrsl(ϕ,t)=ζ(ϕ,t)−δζd/2−δζtidewhere δζtide is the MLWS level with respect to MSL (with the further assumption that, in the absence of paleo tide amplitude information, there has been no significant change in tides). The adopted precision estimate of the observation is Δζ_d/2, added in quadrature to other error sources. Greater uncertainties can arise at times of very rapid sea-level rise if reef formation cannot keep up. In some special morphological forms, the in situ coral can be related more directly to MLWS, as is the case of microatolls where the coral grows up to this limiting level before growing radially outward at the MLWS level (21). Observed records are restricted to the past ∼6–7 ka, when far-field sea-level change was slow enough for the microatoll development to be able to keep up with change (22). The most complete record is from Kiritimati Atoll (23). The principal coral records for the late glacial period are from Barbados (24⇓–26), Tahiti (27, 28), Huon Peninsula of Papua New Guinea (29, 30), the western Indian Ocean (31), and the Great Barrier Reef of Australia, supplemented for the past 7 ka with other observations from the Australian region and the Indian Ocean (Fig. 1 and SI Appendix, Table S1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Distribution of far-field sea-level data for the past 35 ka. (A) Depth−age relationship of all data with 2σ error estimates. (B) Time distribution of the data. (C) Geographic distribution of all far-field coral (red dot) and sediment (black triangles) data.

The sedimentological evidence consists of age−depth relationships of material formed within a known part of the tidal range, such as terrestrial plants representative of the highest tidal zone or microfaunal material representative of a lower tidal range (32) (observations that provide only limiting estimates have been excluded). This evidence has the ability to provide high-precision results if restricted to materials with well-understood depositional ranges, but they also present their own difficulties. Is the dated material within the sediment layer in situ? Has there been age contamination by transport of older carbon into the sediment sequence? Has compaction of the sediment column subsequent to deposition been important? Significant data sources included are from the Sunda Shelf (33, 34), the Indian Ocean (35, 36), the Bonaparte Gulf (37, 38), New Zealand (39), Singapore (40, 41), Malaysia (42), and from the South China Sea to the Bohai Sea. Many of the Asian observations are from large delta systems, and comparisons of such data with observations from adjacent sites often indicate lower levels of the former, corresponding to differential subsidence rates of the order of 1 mm⋅y⁻¹ for the past ∼8 ka. Hence data from the deltas of Chao Phraya of Thailand, Mekong and Red rivers of northern Vietnam, and Pearl and Yangtze rivers of China, have not been used for the purpose of estimating the esl function.

Where possible, the selected data are from tectonically stable regions: away from plate margins and marked by an absence of seismic activity and recent faulting. Where available, evidence for the Last Interglacial (marine isotope stage 5.5) shoreline, formed when, globally, sea levels were ∼4–8 m above present (43, 44), has been used to assess the long-term stability of the region or to correct for tectonic displacement. Both Huon Peninsula and Barbados have been subject to uplift, and Tahiti to subsidence (SI Appendix, Text S1). The eastern Asian margin in many locations is characterized by rifting of Cenozoic age and postrift thermal subsidence (45), but the expected rates are small on time scales of the past 20 ka compared with the GIA signals and observational uncertainties, and no corrections have been applied.

Inversion Results: Mantle Rheology

The observation equation isΔζobs(ϕ,t)+εobs=Δζpredicted(ϕ,t)=Δζesl0(t)+δζesl(t)+δζiceiso(ϕ,t)+δζwateriso(ϕ,t)[3]where ε_obs are corrections to the observations of relative sea level and δζesl(t) is a corrective term to Δζesl0(t), parameterized in the first instance as mean values in time bins of 1,000 y from present to 22 ka BP, with larger bins for the data-sparse intervals 22–26 ka BP, 26–31 ka BP, and 31–36 ka BP (Fig. 1). The inversion for the earth-model parameters E and δζesl(t) is by forward modeling through E space in which for any E_k (k = 1…K), the corresponding δζk,esl(t) is estimated by weighted least squares. The minimum value of the variance function Ψk2 in E space,Ψk2=1M∑m=1M{[Δζobsm−(Δζesl0+δζk,esl+δζk,iso)]/σm}2,[4]is then sought, where δζk,iso is the total isostatic correction for E_k and σm is the SD of the mth observation (m = 1…M). The E space is initially restricted to a three-layered model of the mantle: elastic lithosphere (including the crust) defined by an effective elastic thickness H, upper mantle from the base of the lithosphere to the 670 km seismic discontinuity with a depth-averaged effective viscosity η_um, and a lower mantle viscosity extending down to the core-mantle boundary of depth-averaged effective viscosity η_lm. Elastic moduli and density are based on realistic depth profiles determined from seismic data inversions. These approximate models nevertheless describe well the response of the Earth to surface loading at periods and mantle-lithosphere stress levels of glacial cycles (46, 47). The E search is conducted within the confinesEffective lithospheric thickness:   30≤H1≤140 km.Effective upper-mantle viscosity:   1019≤ηum≤1021 Pas.Effective lower-mantle viscosity:   5×1020≤ηlm≤1024 Pas.[5]that span parameters found in earlier far- and near-field analyses. The theory for the viscoelastic solution of earth deformation is based on the transformation of the elastic formulation into the Laplace domain using the correspondence principle and then inverting the Laplace-domain solution back into the time domain (13, 15, 48). This latter inversion is carried out using a pure collocation method (49) that has the potential to become unstable when the ratio of the relaxation times of the two mantle layers becomes very large, and, for numerical reasons, the above E space is restricted to models for which η_lm/η_um ≤ 500 (see Discussion and Conclusions).

In the first iteration solution, a conservative data-quality criterion |(Δζobs−Δζpredicted)/σobs|>6 has been adopted to examine gross discrepancies without rejecting observations that may point to rapid changes in ice volumes not captured by the parameterized δζ_esl(t) function. Only where such large discrepancies are incompatible with other observations close in location and time are they excluded, and 18 observations, of a total of 992, have been rejected. The Barbados record was not used in early iterations because it is from the outer edge of the geoidal bulge around the North American ice sheet and sea-level response there may be more sensitive to δζiceiso(ϕ,t) than sites further from the ice margin, as well as to potential differences in mantle viscosity beneath continents and oceans and local mantle structure associated with the descending lithospheric slab (50). However, tests with and without this data yielded near-identical results for both E and δζ_esl(t), and, in all subsequent iterations, observations from Barbados (but not from other Caribbean sites closer to the former ice margin) have been included.

The resulting variance function Ψk2 is illustrated in Fig. 2. The minimum value of ∼2.8 exceeds the expected value of unity. About 10% of the observations contribute nearly 50% to this variance, but excluding this information from the analysis does not lead to different E. That Ψk2>1 indicates either an underestimation of the observational uncertainties or unmodeled contributions to sea level. Of the former, conservative estimates of observational accuracies have already been made, but, in some instances, local tectonic or subsidence contributions, or corrections for the coral growth range, may be more important than assumed. Part of the latter may be a consequence of rapid changes in ice volume not well captured by the adopted bin sizes, but solutions with smaller bins give identical results for E. Part may be a consequence of the oversimplification of the depth dependence of the viscosity layering, but solutions with further viscosity layering in the upper mantle do not lead to a variance reduction, and two layered upper-mantle models with their boundary at the ∼400-km seismic discontinuity lead to a minimum variance solution for near-zero viscosity contrast across this boundary. A further possible contributing factor is the assumption of lateral uniformity of viscosity in the mantle. Separate solutions for continental-margin and midocean-island data do not require lateral variation, partly because the island data are limited in their distribution and partly because the resolution of island data for depth dependence of viscosity is intrinsically poor. Solutions excluding the eastern margin of Asia also lead to the same E. Hence the resulting E parameters are descriptive of the mantle response for continental margins and midocean environments across the Indian and eastern Pacific Oceans and the margins of Australia, East Africa, and southern and eastern Asia.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Minimum variance function Ψk2 [4] as function of (A) lithospheric thickness H, (B) upper-mantle viscosity η_um, and (C) lower-mantle viscosity η_lm across E space defined by [5] with η_lm/η_um ≤ 500. Unique solutions are found for H and η_um but two minima are identified for η_lm: a low lower-mantle viscosity solution at η_lm ∼2 × 10²¹ Pa s (red dot) and a high lower-mantle viscosity solution at η_lm ∼10²³ Pa s (blue dot). The corresponding 95% confidence limits Φk2 [6] are defined by the red and blue bands.

Confidence limits for the E across the space [5] are estimated using the statisticΦk2=1M∑m=1M{(Δζk,predictedm−Δζk∗,predictedm)/σm}2[6]where Δζk∗,predictedm are the predicted relative sea levels for observation m and earth model Ek∗ that correspond to the least variance solution [4]. If observational variances are appropriate and the model is correct and complete, the contour Φk2=1.0 defines the 67% confidence limits of the minimum variance solution. Of the E_k, η_um is well constrained (Fig. 2), with the preferred minimum variance occurring at 1.5 × 10²⁰ Pa s [with 95% confidence limits of ∼(1, 2) × 10²⁰ Pa s]. Resolution for lithospheric thickness is less satisfactory and, while the solutions point to a minimum Ψk2 at between 40 km and 70 km, the Φk2≤1 criterion indicates only that very thick lithospheres can be excluded. This lack of resolution is because (i) midocean small-island data have little resolution for H and (ii) the distribution of the observations from near the present continent-ocean boundary is insufficient to fully reflect the gradients in Δζrsl(t) predicted across this boundary. The solution for η_lm is least satisfactory of all (Fig. 2C) and points to two local minima, a “low” viscosity solution at 2 × 10²¹ (7 × 10²⁰ – 4 × 10²¹) Pa s and a “high” viscosity solution of ∼7 × 10²² (1 × 10²² –2 × 10²³) Pa s: An addition of a uniform layer of melt water, whose spatial characteristics are dominated by ocean-basin-scale wavelengths, predominantly stresses the lower mantle. Thus, the question remains whether the geographic variability of the data set is adequate for a complete separation of E and Δζesl(t) or whether a unique separation of ice and earth parameters is even possible. We keep both solutions for the present and examine their implications on the ice volume estimates below.

The δζ_esl(t) for the two solutions (Fig. 3) implies that the global estimate of ice volume needs a correction of ∼5–10% for some epochs, but the far-field analysis alone does not allow this change to be attributed to one ice sheet or another. To evaluate the possible dependence of the estimates of both E and Δζesl(t) on how δζ_esl(t) is distributed, we have considered three options in which: (i) the δζ_esl(t) for the LGM and late-glacial period is attributed to the North American ice sheet, (ii) the same as option i but with the Holocene δζ_esl(t) attributed to Antarctica, and (iii) with δζ_esl(t) attributed to Antarctica for the entire interval. The solution for E and a further corrective term to Δζ_esl is then repeated for each of these so-modified ice sheets. This results in equivalent results for E and to a rapid convergence for the Δζ_esl function (SI Appendix, Table S2), justifying thereby the basic assumption that the far-field analysis is not critically dependent on the distribution of the ice between the component ice sheets.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Low-definition solutions (1,000-y time bins) for the corrective term δζesl(t) for the two lower-mantle viscosity solutions (low-viscosity solution in red with yellow error bars and high-viscosity solution in blue with pale blue error bars).

Inversion Results: Ice-Volume Equivalent Sea Level

The above parameterization of δζ_esl(t) (Fig. 3) results in a low temporal resolution of δζ_esl(t) that could preclude detection of changes in sea level in intervals <1,000 y. Thus, once optimum effective earth parameters are established that are independent of the details of the ice model, we adopt a post E solution processing approach to solve for the incremental δζ_esl(t) in which each observation Δζobsm(ϕ,t) provides an estimate of eslΔζeslm(t)=Δζobsm(ϕ,t)−[Δζk∗,predictedm(ϕ,t)−Δζesl0(t)][7a]with(σeslm)2=(σobsm(ϕ,t))2+(σk∗,predictedm(t))2.[7b]The variance (σk∗,predictedm(t))2 follows from the variance of the predicted values across the E space defined by the Φk2≤1 according to(σk∗,predictedm)2=[∑k(Δζk∗,predictedm−Δζk,predictedm)2/σk2]/∑k1/σk2[8]withσk2=[(1−Φk2)]2.The underlying signal in the resulting noisy time series is then estimated using the transdimensional Markov chain Monte Carlo approach (51) that allows abrupt or rapid changes to be quantified against a background of long-term trends and that infers the probability distributions on the number and location of change points, as well as the trends between change points. This produces an ensemble of candidate time series curves, the collective density of which follows a Bayesian a posteriori probability density function (PDF) defined by the product of a likelihood and a priori PDF on control parameters. The algorithm avoids the need to impose artificial constraints on model complexity by using a flexible parameterization of the time series with a variable number of unknowns but at the same time remains parsimonious in that it eliminates unwarranted detail in the reconstructed signal. The mean of the ensemble is taken as an objective estimate of the underlying “denoised” time series while uncertainty estimates are obtained through appropriate projections of the ensemble. Uncertainty estimates, in terms of probability density functions, are obtained for the entire time signal as well as location of the change points representing abrupt changes in gradient (SI Appendix, Text S2 and Fig. S1).

Fig. 4 illustrates the result corresponding to the high-viscosity lower mantle model, and that for the low-viscosity lower mantle is given in SI Appendix, Fig. S2. Solution with and without the China data yield indistinguishable results, and the introduction of this generally lower-quality material has not distorted the inferences that can be drawn about rates of global sea-level change but, rather, reinforces them.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Solution for the ice-volume esl function and change in ice volume. (A) Individual esl estimates (blue) and the objective estimate of the denoised time series (red line). The Inset gives an expanded scale for the last 9,000 y. (B) The same esl estimate and its 95% probability limiting values. Also shown are the major climate events in the interval [the Last Glacial Maximum (LGM), Heinrich events H1 to H3, the Bølling-Allerød warm period (B-A), and the Younger Dryas cold period (Y-D)] as well as the timing of MWP-1A, 1B, and the 8.2 ka BP cooling event. (C) The 95% probability estimates of the esl estimates. (D) Estimates of sea-level rate of change.

Discussion and Conclusions

On time scales of 10⁵ years and less, sea-level change at tectonically stable regions is primarily a function of changing ice volumes and the Earth’s response to the changing ice-water load, but neither the ice history nor the response function is independently known with the requisite precision for developing predictive models. Observations of sea level through time do provide constraints on the ice and rheology functions, but a complete separation of the two groups of parameters has not yet been achieved. Separation of the analysis into far-field and near-field areas provides some resolution, but ambiguities remain: a consequence of inadequate a priori information on ice margin evolution and ice thickness, observational data that deteriorates in distribution and accuracy back in time, the likelihood of lateral variations in the planet’s rheological response, and the ever-present possibility of tectonic contributions.

We have focused here on the far-field analysis of nearly 1,000 observations of relative sea level for the past 35 ka, underpinned by independent near-field analyses for the individual major ice sheets of the Northern Hemisphere and inferences about the past Antarctic ice sheet. Within the confines of model assumptions and with one caveat, a separation can be achieved of effective rheology parameters for the mantle beneath the oceans and continental margins and the change in total ice volume or ice-volume esl. The caveat is that two solutions are possible, characterized as high-viscosity and low-viscosity lower-mantle models, each with its own esl function (SI Appendix, Fig. S2). One approach to separate earth and ice parameters could be to search for “dip-stick” sites (52) where the glacioisostatic and hydroisostatic components cancel and Δζrsl(ϕ,t)=Δζesl(ϕ,t). However, such contours do exhibit time and E dependence, and there are few observations that actually lie close to the Δζiso=0 contour at any epoch (SI Appendix, Fig. S3).

The esl obtained for the two solutions (Fig. 4 and SI Appendix, Fig. S2) differ in two important respects: (i) at the LGM, the low-viscosity model requires ∼2.7 × 10⁶ km³ (∼7 m esl) more ice at the maximum glaciation than the high-viscosity model, and (ii) during the mid to late Holocene (the past ∼7 ka), the low-viscosity model requires less ice (∼0.5 × 10⁶ km³ or ∼1.3 m esl at 6 ka BP) than the high-viscosity model (the relaxation in the low-viscosity model is more complete than for a high-viscosity model). Of these, the first difference may provide some guide to the choice of solution.

The changes in Antarctic ice volume since the LGM remain poorly known, and published estimates differ greatly: ∼25–35 m of esl from the difference in far-field and Northern Hemisphere near-field estimates of changes in ice volume (53), ∼20 m from the combination of such methods with the inversion of rebound data from the Antarctic margin (54), ∼10–18 m from glaciological models (55, 56), and ∼10 m from combined glaciological and geological modeling (57). In developing our component ice sheets, we have used an iterative approach in which, at any step, the Antarctic ice-volume function is the difference between the far-field derived global estimate and the sum of the Northern Hemisphere and mountain glacier ice volumes for that iteration. This “missing” ice is then distributed within the ice sheet to respect the LGM Antarctic margin (58) and with the assumption that ice elevation profiles across the shelf and into the interior approximate a quasi-parabolic form (59). Thus, the intent of this ice sheet is only to preserve the global ocean−ice mass balance and not to produce a realistic rebound model for far-southern latitudes. The starting iteration of the present far-field solution has an Antarctic contribution to esl of 28 m (equivalent to ∼10⁷ km³ of ice, with a large fraction of this ice grounded on the shelves; see Eq. 1). For the high-viscosity solution this is reduced to ∼23 m whereas it is increased to ∼30 m for the low-viscosity solution. The choice of model could then be made on the basis of how much ice can plausibly be stored in Antarctica during the LGM, consistent with other geological and glaciological constraints. On the basis of the comparison with independent estimates (55⇓–57), we favor the high-viscosity model.

Both solutions point to an ongoing slow increase in ocean volume after 7 ka BP, even though melting of the large Northern Hemisphere ice sheets had largely ceased by this time. The high northern latitude Holocene Climate Optimum peaked between 6 and 4 ka BP, and some Greenland ice-margin retreat occurred at least locally (60, 61). There is also evidence that some Antarctic melting occurred after ∼6 ka BP (62). Furthermore, high and midlatitude mountain glaciers, including remnants of Late Pleistocene arctic and subarctic ice caps, may also have continued to decrease in volume (63). However, in all cases, there is not enough observational evidence to independently constrain global estimates of ice-mass fluctuations during this post-7-ka period to permit discrimination between the two esl solutions.

One option is to impose constraints on the lower-mantle viscosity from sources independent of the far-field results. Our analyses of rebound data from formerly glaciated regions (64) and of deglaciation-induced changes in the dynamic flattening and rotation of the Earth (65⇓–67) are consistent with the high-viscosity results, although these solutions are also ice-model dependent. Inversions of geoid and seismic tomographic data are less sensitive to the choice of ice models (the observed geoid needs corrections for the GIA contribution) and point to an increase in depth-averaged mantle viscosity of one to three orders of magnitude from average upper to lower mantle (68⇓⇓⇓⇓–73). Likewise, inferences of viscosity from the sinking speed of subducted lithosphere also point to high (3–5) × 10²² Pa s values for the lower-mantle viscosity (74). On this basis, and because it yields the smaller volume for the LGM Antarctic ice volume, we adopt the high-viscosity lower-mantle solution.

The adopted esl function illustrated in Fig. 4 and SI Appendix, Table S3, shows many features previously identified in one or more individual records but which, because it is based on geographically well-spaced observations corrected for the isostatic contributions, will reflect more accurately the changes in globally averaged sea level and ice volume. These features are:

• i) A period of a relatively slow fall in sea level from 35 to 31 ka BP followed by a rapid fall during 31–29 ka. This is based on data from Barbados, Bonaparte Gulf, Huon Peninsula (Papua New Guinea), and a few isolated observations from the Malay Peninsula and the Bengal Fan. It points to a period of rapid ice growth of ∼25 m esl within ∼1,000 y to mark the onset of the peak glaciation, consistent with the transition out of the Scandinavian Ålesund Interstadial into the glacial maximum (75), although this ice sheet alone is inadequate to contribute 25 m to esl. Chronologically, the timing of the rapid fall corresponds to the nominal age for the Heinrich H3 event (76, 77).

• ii) Approximately constant or slowly increasing ice volumes from ∼29–21 ka BP. The data for this interval are sparse but are from geographically well-distributed sites (Barbados, Bonaparte, Bengal Fan, East China Sea, and Maldives). The slow increase in ice volume is consistent with eastward and southward expansion of the Scandinavian ice sheet during the LGM (78) as well as with the southward advance of the Laurentide ice sheet (79). The esl reaches its lowest value of ∼134 m at the end of this interval, corresponding to ∼52 × 10⁶ km³ more grounded ice—including on shelves—at the LGM than today. This is greater than the frequently cited −125 m (e.g., 24, 80) that is usually based on observations uncorrected for isostatic effects. Heinrich event H2 at ∼19.5–22 ¹⁴C ka (∼25 ka BP) (77, 81, 82) is not associated with a recognizable sea-level signal.

• iii) Onset of deglaciation at ∼21–20 ka BP with a short-lived global sea-level rise of ∼10–15 m before 18 ka. The evidence comes from the Bonaparte Gulf (37, 38), has been identified elsewhere (83, 84), and is supported by isolated observations from five other locations (Bengal Bay, Cape St Francis (South Africa), offshore Sydney, Barbados, Maldives) which, although less precise than the principal data set, spread the rise over a longer time interval than originally suggested. Chronologically, this rise occurs substantially later than the H2 event.

• iv) A short period of near-constant sea level from ∼18–16.5 ka BP. Support for this comes from observations from Barbados, the Sunda Shelf, Bonaparte Gulf, Mayotte, and Cape St Francis.

• v) A major phase of deglaciation from ∼16.5–7 ka BP. The total esl change in this interval is ∼120 m, at an average rate of ∼12 m⋅ka⁻¹, corresponding to a reduction of grounded ice volume of ∼45 × 10⁶ km³). Within this interval, significant departures from the average occur.

• vi) A rise of ∼25 m from ∼16.5–15 ka BP at the long-term average rate of ∼12 m⋅ka⁻¹. The data are from Sunda, Tahiti, the East China Sea, Mayotte, and Australia. Chronologically, the onset of this rise occurs at the time of the H1 event dated at 16.8 ka (85) or 16 ka (86). This period of rising sea level is followed by a short period (∼500–600 y) of near-constant sea level.

• vii) A high rate of sea-level rise starting at ∼14.5 ka BP of ∼500 y duration. The onset occurs at the start of the Bølling−Allerød warm period. Its duration could be <500 y because of uncertainties in chronology, and the globally averaged rise in sea level of ∼20 m occurs at a rate of ∼40 mm⋅y⁻¹ or greater. This pulse, MWP-1A, has been identified separately in the records of Barbados (24), Sunda (33), and Tahiti (28, 87). Spatial variation in its amplitude can be expected because of the planet’s elastic and gravitational response to rapid unloading of ice in either or both of the two hemispheres (88) with, based on the ice−earth models used here, model-predicted values ranging from ∼14 m for Barbados to ∼20 m for Tahiti (SI Appendix, Fig. S4). This compares with observational values of ∼15–20 m (24, 28) for Barbados and 12–22 m for Tahiti (28). Observational uncertainties remain large, including differences in the timing of this event as recorded at the different localities, and it is not possible from this evidence to ascertain the relative importance of the contribution of the two hemispheres to MWP-1A.

• viii) A period of sea-level rise from ∼14 to ∼12.5 ka BP of ∼20 m in 1,500 y. The rate of rise is near the long-term average. Data are relatively dense in this interval and come from well-distributed sites (Barbados, Tahiti, Sunda, Huon Peninsula, Australia and New Zealand, Indian Ocean, and the Yellow and East China seas).

• ix) A period of a much reduced rate of rise from ∼12.5–11.5 ka BP. This short duration pause in the sea-level rise has been tenuously noted before in both composite (89) and individual (27) records. The chronology corresponds to the timing of the Younger Dryas stadial of the Northern Hemisphere when retreat of the Northern Hemisphere ice sheets ceased momentarily.

• x) A period from ∼11.4–8.2 ka BP of near-uniform global rise. The average rate of rise during this 3.3 ka interval was ∼15 m⋅ka⁻¹ with little convincing evidence of variations in this rate. A rapid rise, MWP-1B, has been reported at ∼11.3 ka but remains elusive (27) and is not seen in the composite record other than as a slightly higher rate of increase to ∼16.5 mm⋅y⁻¹ for a 500-y period immediately after the Younger Dryas period.

• xi) A reduced rate of sea-level rise for 8.2–6.7 ka BP. This is consistent with the final phase of North American deglaciation at ∼7 ka BP. A marked cooling event has been recorded at 8.2 ka BP in Greenland and North Atlantic cores (90), but there is no suggestion in the sea-level record of a corresponding fall or slowdown in global sea-level rise. The detailed local record from Singapore from 8.5 to 6 ka BP (40, 41) is consistent with the global rates within this interval except that a period of near-zero rise from 7.8 to 7.4 ka is not seen globally, possibly lost in the noise of other observations at around this time, possibly because it reflects local phenomena (Fig. 4).

• xii) A progressive decrease in rate of rise from 6.7 ka to recent time. This interval comprises nearly 60% of the database (Fig. 1). The total global rise for the past 6.7 ka was ∼4 m (∼1.2 × 10⁶ km³ of grounded ice), of which ∼3 m occurred in the interval 6.7–4.2 ka BP with a further rise of ≤1 m up to the time of onset of recent sea-level rise ∼100–150 y ago (91, 92). In this interval of 4.2 ka to ∼0.15 ka, there is no evidence for oscillations in global-mean sea level of amplitudes exceeding 15–20 cm on time scales of ∼200 y (about equal to the accuracy of radiocarbon ages for this period, taking into consideration reservoir uncertainties; also, bins of 200 y contain an average of ∼15 observations/bin). This absence of oscillations in sea level for this period is consistent with the most complete record of microatoll data from Kiritimati (23). The record for the past 1,000 y is sparse compared with that from 1 to 6.7 ka BP, but there is no evidence in this data set to indicate that regional climate fluctuations, such as the Medieval warm period followed by the Little Ice Age, are associated with significant global sea-level oscillations.