Significance Several areas of earth science require knowledge of the fluctuations in sea level and ice volume through glacial cycles. These include understanding past ice sheets and providing boundary conditions for paleoclimate models, calibrating marine-sediment isotopic records, and providing the background signal for evaluating anthropogenic contributions to sea level. From ∼1,000 observations of sea level, allowing for isostatic and tectonic contributions, we have quantified the rise and fall in global ocean and ice volumes for the past 35,000 years. Of particular note is that during the ∼6,000 y up to the start of the recent rise ∼100−150 y ago, there is no evidence for global oscillations in sea level on time scales exceeding ∼200 y duration or 15−20 cm amplitude.

Abstract The major cause of sea-level change during ice ages is the exchange of water between ice and ocean and the planet’s dynamic response to the changing surface load. Inversion of ∼1,000 observations for the past 35,000 y from localities far from former ice margins has provided new constraints on the fluctuation of ice volume in this interval. Key results are: (i) a rapid final fall in global sea level of ∼40 m in <2,000 y at the onset of the glacial maximum ∼30,000 y before present (30 ka BP); (ii) a slow fall to −134 m from 29 to 21 ka BP with a maximum grounded ice volume of ∼52 × 106 km3 greater than today; (iii) after an initial short duration rapid rise and a short interval of near-constant sea level, the main phase of deglaciation occurred from ∼16.5 ka BP to ∼8.2 ka BP at an average rate of rise of 12 m⋅ka−1 punctuated by periods of greater, particularly at 14.5–14.0 ka BP at ≥40 mm⋅y−1 (MWP-1A), and lesser, from 12.5 to 11.5 ka BP (Younger Dryas), rates; (iv) no evidence for a global MWP-1B event at ∼11.3 ka BP; and (v) a progressive decrease in the rate of rise from 8.2 ka to ∼2.5 ka BP, after which ocean volumes remained nearly constant until the renewed sea-level rise at 100–150 y ago, with no evidence of oscillations exceeding ∼15–20 cm in time intervals ≥200 y from 6 to 0.15 ka BP.

The understanding of the change in ocean volume during glacial cycles is pertinent to several areas of earth science: for estimating the volume of ice and its geographic distribution through time (1); for calibrating isotopic proxy indicators of ocean volume change (2, 3); for estimating vertical rates of land movement from geological data (4); for examining the response of reef development to changing sea level (5); and for reconstructing paleo topographies to test models of human and other migrations (6). Estimates of variations in global sea level come from direct observational evidence of past sea levels relative to present and less directly from temporal variations in the oxygen isotopic signal of ocean sediments (7). Both yield model-dependent estimates. The first requires assumptions about processes that govern how past sea levels are recorded in the coastal geology or geomorphology as well as about the tectonic, isostatic, and oceanographic contributions to sea level change. The second requires assumptions about the source of the isotopic or chemical signatures of marine sediments and about the relative importance of growth or decay of the ice sheets, of changes in ocean and atmospheric temperatures, or from local or regional factors that control the extent and time scales of mixing within ocean basins.

Both approaches are important and complementary. The direct observational evidence is restricted to time intervals or climatic and tectonic settings that favor preservation of the records through otherwise successive overprinting events. As a result, the records become increasingly fragmentary backward in time. The isotopic evidence, in contrast, being recorded in deep-water carbonate marine sediments, extends further back in time and often yields near-continuous records of high but imprecise temporal resolution (8). However, they are also subject to greater uncertainty because of the isotope signal’s dependence on other factors. Comparisons for the Holocene for which the direct measures of past sea level are relatively abundant, for example, exhibit differences both in phase and in noise characteristics between the two data [compare, for example, the Holocene parts of oxygen isotope records from the Pacific (9) and from two Red Sea cores (10)].

Past sea level is measured with respect to its present position and contains information on both land movement and changes in ocean volume. During glacial cycles of ∼105 y, the most important contribution with a global signature is the exchange of mass between the ice sheets and oceans, with tectonic vertical land movements being important mainly on local and regional scales. Global changes associated with mantle convection and surface processes are comparatively small on these time scales but become important in longer, e.g., Pliocene-scale, periods (11). Changes in ocean volume associated with changing ocean temperatures during a glacial period are also small (12).

The sea-level signal from the glacial cycle exhibits significant spatial variability from its globally averaged value because of the combined deformation and gravitational response of the Earth and ocean to the changing ice-water load. During ice-sheet decay, the crust rebounds beneath the ice sheets and subsides beneath the melt-water loaded ocean basins; the gravitational potential and ocean surface are modified by the deformation and changing surface load; and the planet’s inertia tensor and rotation changes, further modifying equipotential surfaces. Together, this response of the earth-ocean system to glacial cycles is referred to as the glacial isostatic adjustment (GIA) (13⇓⇓⇓–17). The pattern of the spatial variability is a function of the Earth's rheology and of the glacial history, both of which are only partly known. In particular, past ice thickness is rarely observed and questions remain about the timing and extent of the former ice sheets on the continental shelves. The sea-level response within, or close to, the former ice margins (near-field) is primarily a function of the underlying rheology and ice thickness while, far from the former ice margins (far-field), it is mainly a function of earth rheology and the change in total ice volume through time. By an iterative analysis of observational evidence of the past sea levels, it becomes possible to improve the understanding of the past ice history as well as the Earth's mantle response to forces on a 104 y to 105 y time scale.

In this paper, we address one part of the Earth’s response to the glacial cycle: the analysis of far-field evidence of sea-level change to estimate the variation in ice and ocean volumes from the lead into the Last Glacial Maximum (LGM) at ∼35,000 y ago (35 ka BP) to the start of the instrumental records. Such analyses can either be of high-resolution records of a single data type from a single location or of different sea-level indicators from many different locations. We have adopted the latter approach. Most sea-level indicators provide only lower (e.g., fossil coral) or upper (e.g., fossil terrestrial plants) limiting values, and multiple data-type analyses of both upper and lower limiting measurements are less likely to be biased toward one or the other limit. Tectonic displacement of the crust is always a potential contaminator of the sought signal, and a multiplicity of data from tectonically “stable” regions is more likely to average out any undetected tectonic effects as well as any uncertainties in the above-mentioned GIA contributions.

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 ∫ t 1 A o ( t ) d Δ M ice d t d t [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 ) + δ ζ ice iso ( ϕ , t ) + δ ζ water iso ( ϕ , t ) [2]where δ ζ ice iso ( ϕ , t ) and δ ζ water iso ( ϕ , 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, δ ζ water iso ( ϕ , 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 δ ζ water iso ( ϕ , 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, δ ζ water iso ( ϕ , t ) is about 15–25% of Δ ζ e s l ( t ) and leads to near-coastal observations of LGM sea level being systematically less than the corresponding Δ ζ e s l ( 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 δ ζ i c e i s o ( ϕ , 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, δ ζ i c e i s o ( ϕ , t ) shows less spatial variability than δ ζ w a t e r i s o ( ϕ , 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 Δ ζ r s l = Δ ζ e s l and using [1] estimate a first approximation of ice mass or volume function Δ V i c e ( 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 δ ζ e s l to the nominal Δ ζ e s l 0 corresponding to step ii, yielding an improved approximation of ice volume function Δ V i c e ( 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 Δ V i c e ( 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 14C 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 Δ ζ r s l ( ϕ , t ) = ζ ( ϕ , t ) − δ ζ d / 2 − δ ζ t i d e where δ ζ t i d e 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). 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−1 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 ) = Δ ζ esl 0 ( t ) + δ ζ esl ( t ) + δ ζ ice is o ( ϕ , t ) + δ ζ water iso ( ϕ , t ) [3]where ε obs are corrections to the observations of relative sea level and δ ζ e s l ( t ) is a corrective term to Δ ζ e s l 0 ( 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 δ ζ e s l ( t ) is by forward modeling through E space in which for any E k (k = 1…K), the corresponding δ ζ k , e s l ( t ) is estimated by weighted least squares. The minimum value of the variance function Ψ k 2 in E space, Ψ k 2 = 1 M ∑ m = 1 M { [ Δ ζ obs m − ( Δ ζ esl 0 + δ ζ k , esl + δ ζ k , iso ) ] / σ m } 2 , [4]is then sought, where δ ζ k , i s o 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 confines Effective lithospheric thickness : 30 ≤ H 1 ≤ 140 km . Effective upper-mantle viscosity : 10 19 ≤ η um ≤ 10 21 Pa s . Effective lower-mantle viscosity : 5 × 10 20 ≤ η lm ≤ 10 24 Pa s . [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 − Δ ζ predicte d ) / σ o b s | > 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 δ ζ i c e i s o ( ϕ , 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 Ψ k 2 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 Ψ k 2 > 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. Fig. 2. Minimum variance function Ψ k 2 [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 × 1021 Pa s (red dot) and a high lower-mantle viscosity solution at η lm ∼1023 Pa s (blue dot). The corresponding 95% confidence limits Φ k 2 [6] are defined by the red and blue bands. Confidence limits for the E across the space [5] are estimated using the statistic Φ k 2 = 1 M ∑ m = 1 M { ( Δ ζ k , predicted m − Δ ζ k ∗ , predicted m ) / σ m } 2 [6]where Δ ζ k ∗ , predicted m are the predicted relative sea levels for observation m and earth model E k ∗ that correspond to the least variance solution [4]. If observational variances are appropriate and the model is correct and complete, the contour Φ k 2 = 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 × 1020 Pa s [with 95% confidence limits of ∼(1, 2) × 1020 Pa s]. Resolution for lithospheric thickness is less satisfactory and, while the solutions point to a minimum Ψ k 2 at between 40 km and 70 km, the Φ k 2 ≤ 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 Δ ζ r s l ( 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 × 1021 (7 × 1020 – 4 × 1021) Pa s and a “high” viscosity solution of ∼7 × 1022 (1 × 1022 –2 × 1023) 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 Δ ζ e s l ( 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. Fig. 3. Low-definition solutions (1,000-y time bins) for the corrective term δ ζ e s l ( 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 Δ ζ o b s m ( ϕ , t ) provides an estimate of esl Δ ζ e s l m ( t ) = Δ ζ o b s m ( ϕ , t ) − [ Δ ζ k ∗ , p r e d i c t e d m ( ϕ , t ) − Δ ζ e s l 0 ( t ) ] [7a]with ( σ e s l m ) 2 = ( σ o b s m ( ϕ , t ) ) 2 + ( σ k ∗ , p r e d i c t e d m ( t ) ) 2 . [7b]The variance ( σ k ∗ , p r e d i c t e d m ( t ) ) 2 follows from the variance of the predicted values across the E space defined by the Φ k 2 ≤ 1 according to ( σ k ∗ , p r e d i c t e d m ) 2 = [ ∑ k ( Δ ζ k ∗ , p r e d i c t e d m − Δ ζ k , p r e d i c t e d m ) 2 / σ k 2 ] / ∑ k 1 / σ k 2 [8]with σ k 2 = [ ( 1 − Φ k 2 ) ] 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. 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 105 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 Δ ζ r s l ( ϕ , t ) = Δ ζ e s l ( ϕ , t ) . However, such contours do exhibit time and E dependence, and there are few observations that actually lie close to the Δ ζ i s o = 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 × 106 km3 (∼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 × 106 km3 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 ∼107 km3 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) × 1022 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 6 km 3 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 14 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 −1 , corresponding to a reduction of grounded ice volume of ∼45 × 10 6 km 3 ). 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 −1 . 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 −1 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 −1 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 −1 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 × 106 km3 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.

Acknowledgments In addition to the support from our home institutions, this research has been partially funded by the Chaires Internationales de Recherche Blaise Pascal de l'État Français et la Région Île de France, the International Balzan Foundation, and the Australian Research Council. The development of the Baysian partition modeling software was partly supported by the Auscope Inversion Laboratory, funded by the Australian Federal Government.

Footnotes This contribution is part of the special series of Inaugural Articles by members of the National Academy of Sciences elected in 2009.

Author contributions: K.L. designed research; K.L., H.R., A.P., Y.S., and M.S. performed research; Y.S. contributed field data; K.L., H.R., A.P., Y.S., and M.S. analyzed data; K.L. and H.R. wrote the paper; and A.P. and M.S. developed software.

Reviewers: E.B., CEREGE; J.X.M., Harvard University; and P.U.C., Oregon State University.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1411762111/-/DCSupplemental.