Hay et al. (2015)

Title:

Probabilistic reanalysis of twentieth-century sea-level rise

Key Points:

• Reconstruction of the global mean sea-level changes since 1900, from combining the probability distributions of a set of tide gauge records, refines estimates and is more consitent with estimates based on the sum of contributions

• The global mean sea level rose by about 1.2 ± 0.2 mm/year between 1901 and 1990 and by 3.0 ± 0.7 mm/year between 1993 and 2010.

• The rate of sea-level rise has been accelerating

Keywords:

sea level reconstruction, tide gauge, global mean sea level rise, probabilistic framework

Corresponding author:

Carling C. Hay.

Citation:

Hay, C. C., Morrow, E., Kopp, R. E., & Mitrovica, J. X. (2015). Probabilistic reanalysis of twentieth-century sea-level rise. Nature, 517(7535), 481–484. doi:10.1038/nature14093

URL:

https://www.nature.com/articles/nature14093

Abstract

Estimating and accounting for twentieth-century global mean sea-level (GMSL) rise is critical to characterizing current and future human-induced sea-level change. Several previous analyses of tide gauge records [1,2,3,4,5,6] — employing different methods to accommodate the spatial sparsity and temporal incompleteness of the data and to constrain the geometry of long-term sea-level change — have concluded that GMSL rose over the twentieth century at a mean rate of 1.6 to 1.9 millimetres per year. Efforts to account for this rate by summing estimates of individual contributions from glacier and ice-sheet mass loss, ocean thermal expansion, and changes in land water storage fall significantly short in the period before 1990 [7]. The failure to close the budget of GMSL during this period has led to suggestions that several contributions may have been systematically underestimated [8]. However, the extent to which the limitations of tide gauge analyses have affected estimates of the GMSL rate of change is unclear. Here we revisit estimates of twentieth-century GMSL rise using probabilistic techniques [9,10] and find a rate of GMSL rise from 1901 to 1990 of 1.2 ± 0.2 millimetres per year (90% confidence interval). Based on individual contributions tabulated in the Fifth Assessment Report [7] of the Intergovernmental Panel on Climate Change, this estimate closes the twentieth-century sea-level budget. Our analysis, which combines tide gauge records with physics-based and model-derived geometries of the various contributing signals, also indicates that GMSL rose at a rate of 3.0 ± 0.7 millimetres per year between 1993 and 2010, consistent with prior estimates from tide gauge records4. The increase in rate relative to the 1901–90 trend is accordingly larger than previously thought; this revision may affect some projections11 of future sea-level rise.

Editorial Summary: Twentieth century sea levels revisited

Rates of sea-level rise calculated from tide gauge data tend to exceed bottom-up estimates derived from summing loss of ice mass, thermal expansion and changes in land storage. Carling Hay et al. provide a statistical reassessment of the tide gauge record — which is subject to bias due to sparse and non-uniform geographic coverage and other uncertainties — and conclude that sea-level rose by about 1.2 millimetres per year from 1901 to 1990. This is slightly lower than prior estimates and is consistent with the bottom-up estimates. The same analysis applied to the period 1993–2010, however, indicates a sea-level rise of about three millimetres per year, consistent with other work and suggesting that the recent acceleration in sea-level rise has been greater than previously thought.

Main

Tide gauges provide records of local sea-level changes that, in the case of some sites, extend back to the eighteenth century [12,13,14]. However, using the database of tide gauge records [15] to estimate historical GMSL rise (defined as the increase in ocean volume normalized by ocean area) is challenging. Tide gauges sample the ocean sparsely and non-uniformly, with a bias towards coastal sites and the Northern Hemisphere, and with few sites at latitudes greater than 60° (see, for example, refs 4, 9). In addition, tide gauge time series show significant inter-annual to decadal variability, and they are characterized by missing data (that is, intervals without observations at the start, middle or end of a time series). From the perspective of estimating GMSL changes, the data are contaminated by local and regional signals due to ongoing glacial isostatic adjustment (GIA) associated with past ice ages [16,17], the spatially non-uniform pattern of sea-level rise associated with changes in contemporary land ice sources[18,19,20,21], ocean/atmosphere dynamics [22], and other local factors including tectonics, sediment compaction, groundwater pumping and harbour development.

Different approaches have been used to address these complexities in efforts to estimate twentieth-century GMSL rise [23]. These include averaging rates at sites with the longest records [1,2], averaging rates determined from regional binning of records [3], incorporating shorter records into the analysis to distinguish between secular trends and decadal-scale variability [3], and using altimetry records to determine dominant sea-level geometries and then using tide gauge records to estimate the time-varying amplitudes of these geometries [4,5]. In most cases, other criteria were applied to cull the tide gauge sites adopted in the analysis (for example, excluding sites near tectonic activity or major urban centres).

Estimates of twentieth-century GMSL rise from these previous analyses range from 1.6 to 1.9 mm yr−1 (refs 1, 2, 3, 4, 5, 6) and define an important enigma. Independent model- and data-based estimates of the individual sources of GMSL, including mass flux from glaciers and ice sheets, thermal expansion of oceans, and changes in land water storage, are insufficient to account for the GMSL rise estimated from tide gauge records8, particularly before 19907. For example, a tabulation of contributions to GMSL rise from 1901 to 1990 in the Fifth Assessment Report (AR5; ref. 7) of the Intergovernmental Panel of Climate Change (IPCC) total 0.5 ± 0.4 mm yr−1 (90% confidence interval, CI) less than a recent tide gauge derived rate of 1.5 ± 0.2 mm yr−1 (90% CI) estimated by Church and White4 for the same period (the confidence range for this estimate is taken from AR5; refs 7 and 23). Using IPCC terminology, the latter suggests that it is ‘extremely likely’ (probability P = 95%) that GMSL rise from 1901 to 1990 was greater than 1.3 mm yr−1, although the bottom-up sum of contributions is ‘likely’ (P > 67%) below this level. The above discrepancy has been attributed to underestimation of almost all possible sources: thermal expansion, glacier mass balance, and Greenland or Antarctic ice sheet mass balance7,8.

In this Letter, we revisit the analysis of GMSL since the start of the twentieth century using Kalman smoothing9 (KS; see Methods). This statistical technique naturally accommodates spatially sparse and temporally incomplete sampling of a global sea-level field, provides a rigorous, probabilistic framework for uncertainty propagation, and can correct for a distribution of GIA and ocean models. We applied the approach to analyse annual records from 622 tide gauges included in the Permanent Service for Mean Sea Level (PSMSL) Revised Local Reference database15,24 and reconstruct the global field of sea-level change for each year from 1900 to 2010.

To examine the skill with which the KS reconstruction reproduces the tide gauge observations, we compute the time series of residuals at each tide gauge site and examine the distribution of the mean residual (that is, bias) for each site (Fig. 1a). The mean of the mean residuals across all 622 observations is 0.3 mm, with a standard deviation of 5.1 mm, indicating minimal systemic bias.

Figure 1: Fit of the KS-based reconstruction of sea level to the tide gauge record. a, Histogram of mean residuals (mm) between the sea-level reconstruction and the tide gauge observations at all 622 sites. The mean of all mean residuals is 0.3 ± 5.1 mm (±1 s.d.). b–f, Time series of reconstructed annual sea level (black lines, KS mean estimate; grey shading, 1σ uncertainty) at New York, USA (b), Fremantle, Australia (c), Zemlia Bunge, Russia (d), Vaasa, Finland (e), and Champlain, Canada (f), together with the associated annual mean tide gauge observations (red lines).

Comparing reconstructions and tide gauge observations at a selection of individual sites (Fig. 1b–f) shows generally excellent agreement, although there are a small number of outliers. An example outlier is the Champlain tide gauge (Fig. 1f), which has a mean residual of 52 mm. This particular misfit (also evident at other sites in the vicinity) can be attributed to the St Lawrence being a regulated water system where flow is dominated by anthropogenic control rather than global-scale climate dynamics25. The eight sites that have mean residuals greater than ±3σ (15 mm) from the mean exhibit an average interannual sea-level variability (estimated as the standard deviation after detrending the tide gauge observations) of ±130 mm, more than triple the mean inter-annual variability of ±40 mm across all sites. Although these outliers have large inter-annual variability, the site-specific variability is incorporated into the covariances computed in the probabilistic reconstruction, and the uncertainties in the estimated sea-level trends at these sites reflect this.

The sum of the KS-estimated GMSL changes associated with the mass balance of the Greenland and Antarctic ice sheets, the mass balance of 18 mountain glacier regions, and thermal expansion (Fig. 2, blue line and shading; see Methods) is characterized by an average GMSL rate of 1.2 ± 0.2 mm yr−1 (90% CI) for 1901–90. As shown in Fig. 3, this is significantly lower than the estimates of 1.5 ± 0.2 mm yr−1 from Church and White4 (magenta line in Fig. 2) and 1.9 mm yr−1 from Jevrejeva et al.3 (red line in Fig. 2). The KS-estimated acceleration is 0.017 ± 0.003 mm yr−2, larger than our estimates based on the Church and White4 (0.009 ± 0.002 mm yr−2) and Jevrejeva et al.3 (0.011 ± 0.006 mm yr−2) time series (see Methods).

Figure 2: Time series of GMSL for the period 1900–2010. Shown are estimates of GMSL based on KS (blue line), GPR (black line), Church and White4 (magenta line) and Jevrejeva et al.3 (red line). Shaded regions show ±1σ pointwise uncertainty. Inset, trends for 1901–90 and 1993–2010, and accelerations, all with 90% CI. Confidence intervals for Church and White4 are from refs 7 and 23. Confidence intervals were not available for Jevrejeva et al.3; data in this reference ends in 2002, so the rate quoted here for 1993–2010 is actually for 1993–2002. Since the GPR methodology outputs decadal sea level, no trend is estimated for 1993–2010. Accelerations are consistently estimated from the KS, GPR, and GMSL time series in refs 3 and 4 (see Methods) from 1901 to the end of each reconstruction.

Figure 3: Comparison of mean GMSL rates for 1901–90. Shown are estimates of GMSL rise for the period 1901–90 obtained from six different sampling methods along with previously published rates (see main text for description of each). The box covers the 1σ uncertainty range, while the bars represent the 90% CI. In the case of the Jevrejeva et al.3 estimate, uncertainties and confidence intervals were not available.

Church and White [4] combined stationary empirical orthogonal functions (EOFs), computed from ∼20 years of satellite altimetry data spanning latitudes up to about ±60°, with amplitudes estimated from sparse tide gauge observations. Given the relatively short duration of the altimeter record, the EOFs may be dominated by patterns due to interannual variability rather than the geometry associated with long-term sea-level change26,27. Jevrejeva et al.3 used tide gauge records to compute regional sea-level means and from these computed a global average. Both methodologies involve spatially sparse, temporally incomplete sampling of the global sea-level field, which introduces a potentially significant bias into estimates of GMSL. The KS technique differs from these approaches by using the spatial information inherent in the observations to infer the weights associated with the individual, underlying contributions to the sea-level change. The method extracts global information from the sparse field by taking advantage of the physics-based and model-derived geometry of the contributing processes, thereby reducing the potential for sampling bias.

To understand the origin of the differences between the KS estimate and the higher values of refs 3 and 4, and in particular to quantify the impact of regional binning, spatial sparsity and missing data, we performed several tests.

First, we applied to the KS global sea-level reconstruction a regional binning algorithm similar to that of Jevrejeva et al.3. In particular, we sampled the reconstruction at the locations of the 622 tide gauge sites, imposed sections of missing data consistent with the PSMSL data availability15, binned the tide gauges into 12 ocean regions, and averaged across these regions to compute a GMSL curve. The resulting estimate of the mean GMSL rate from 1901 to 1990 (Fig. 3; ‘KS PSMSL sampling’), 1.6 ± 0.4 mm yr−1 (90% CI), is significantly closer to the estimate of Jevrejeva et al.3, indicating that combined spatial sparsity and missing data generate an upward bias in estimates of GMSL rates (Fig. 3). Second, we performed a bootstrapping test that repeated the above algorithm for tide gauge subsets ranging from 25 to 600 sites that confirmed this result (see Methods and Extended Data Fig. 3). We also implemented a test to estimate the possible bias in the estimate of GMSL rate introduced in the EOF analysis of Church and White4 (see Methods; Fig. 3; ‘KS EOF’); the result was consistent with the difference between the KS and Church and White4 results in Fig. 2.

We performed several other tests to explore the impact of sparsity and missing data on the estimates. Specifically, we applied the binning algorithm as described above but without imposing sections of missing data. The resulting mean GMSL rate estimate for 1901–90 was 1.0 ± 0.4 mm yr−1, close to the KS result (Fig. 3; ‘KS 622 sites, no missing data’). Third, we sampled the full reconstruction at a large number of globally distributed sites—that is, the sampling was not confined to the tide gauge sites and no sections of missing data were imposed on the time series—and performed the same regional binning and averaging (‘KS global reconstruction’). The resulting rate estimate, 1.2 ± 0.1 mm yr−1, was identical to the KS result (Fig. 3). This indicates that regional binning of estimates, in the absence of sparsity and missing data, does not introduce a significant bias.

To assess the robustness of our probabilistic reanalysis, we also performed a second, independent statistical analysis based on Gaussian process regression28 (GPR), a technique that also naturally accommodates data sparsity and gaps, and incorporates a suite of GIA and ocean models (see Methods; black line in Fig. 2). The mean GMSL rate for 1901–90 estimated from the GPR analysis, 1.1 ± 0.4 mm yr−1, is consistent with the results of the KS analysis (Fig. 3).

Previous analyses appear to have overestimated the mean GMSL rate over the twentieth century. The KS estimate for the period 1901–90 indicates that it is ‘very likely’ (probability P = 90%) that the rate of GMSL rise during this period was between 1.0 and 1.4 mm yr−1. This estimate closes the sea-level budget for 1901–90 estimated in AR5 (ref. 7) without appealing to an underestimation of individual contributions from ocean thermal expansion, glacier melting, or ice sheet mass balance. Moreover, it may contribute to the ultimate resolution of Munk’s sea-level enigma28 (defined by the argument that Earth rotation measurements and bounds on ocean warming are inconsistent with a rate of sea-level rise beginning in the late nineteenth century of 1.5–2.0 mm yr−1), since it may lower the signal of twentieth century ice melting in Earth rotation measurements.

In contrast, for the period 1993–2010—which coincides with the era of satellite altimetry measurements of sea surface height changes29—the KS estimate is consistent with previous results (Fig. 2). The KS estimate, 3.0 ± 0.7 mm yr−1 (90% CI), is essentially identical to the tide gauge analysis of Church and White4 (2.8 ± 0.5 mm yr−1; ref. 23). It is also consistent with the estimate based on TOPEX and Jason altimeter measurements (3.2 ± 0.4 mm yr−1; ref. 29 as cited by ref. 23 for the period 1993–2010, see also ref. 7).

To assess the anomalous nature of recent sea-level change, we compute 15-year rates through the KS-derived GMSL time series in Fig. 2 from 1901 to 2010. Figure 4 shows both the time series and distribution of these 96 rates, where the 5 most recent time windows are shown in red. The former is in qualitative agreement with a previous inference of multi-decadal trends in acceleration during the twentieth century30. While the rates show significant variability, the rate for the 1996–2010 time window, 3.1 mm yr−1, is the largest of all computed 15-year rates.

Figure 4: Moving 15-year averages of GMSL rate estimated using the KS reconstruction of sea level across the entire interval 1901–2010. The x-axis represents the mid-point of each 15-year averaging window, and the shading gives the 1σ uncertainty range. Inset, histogram of 15-year mean GMSL rate estimates (mm yr−1) for all time windows. The five most recent windows are shown in red.

We have revisited twentieth century GMSL rise using probabilistic techniques that combine sea-level records with physics-based and model-derived geometries of the contributing processes. Our estimated GMSL trend for the period 1901–90 (1.2 ± 0.2 mm yr−1) is lower than previous estimates, indicating that the rate of GMSL rise during the last two decades represents a more significant increase than previously recognized. Projections of future sea-level rise based on the time series of historical GMSL, notably semi-empirical approaches11, should accordingly be revisited.

Methods

Probabilistic estimation methods

Kalman smoothing (KS) and Gaussian process regression (GPR), both discussed in detail below, share three advantages over the approaches taken in traditional tide gauge analyses. First, the Bayesian nature of both approaches naturally accommodates the spatiotemporal changes in the availability of the sea-level records (that is, sparsity and missing data). Second, the probabilistic approaches correct for a distribution of GIA and ocean models rather than adopting a specific model for each process, and they thus reduce a potentially important bias in previous estimates of the GMSL change17,31. Last, as both methods are fully probabilistic, they allow for the propagation of measurement and inferential uncertainties and correlations throughout the complete analysis time period. Despite these commonalities, the implementations of KS and GPR differ significantly.

Kalman smoother

The KS methodology is divided into four steps9, the first three of which are repeated by employing the spatial fields of GIA and ocean dynamic models from all possible combinations of 161 different Earth rheological models and 6 global climate model (GCM) simulations from CMIP5 (ref. 32) (see below for details of the rheological and climate models). First, a priori model estimates of both local sea level and the individual mass contributions from the Greenland, West Antarctic and East Antarctic ice sheets, as well as 18 major mountain glacier regions, are recursively corrected by tide gauge observations as the estimates are propagated forward through time. The local sea level is linked to the individual mass contributions through the unique spatial patterns, or ‘fingerprints,’ of sea-level change associated with rapid mass loss from land-based ice18,19,20,21. The forward step yields an estimate of local sea level and land ice contributions at each time slice, conditional on all earlier observations and a particular combination of GIA and GCM models. Second, the procedure is run backward in time, with the initial state estimate being the last estimate from the first step. The third, smoothing step optimally combines the results of the first two passes based upon the uncertainties of the respective estimates. The result is an estimate of local sea levels and land ice contributions conditional upon the entire set of observations and specific pairings of GIA and GCM models. Finally, the results from different GIA/GCM combinations are linearly combined, weighted by their likelihood, to yield an a posteriori probability distribution for local sea levels and land ice contributions, conditional upon the tide gauge observations.

A comprehensive discussion of our application of the KS technique to the analysis of tide gauge measurements is given in ref. 9, which also includes synthetic tests to assess the performance of the procedure. Several subsequent refinements of this approach are summarized below.

Reference 9 defined the state vector to include estimates of sea level at every tide gauge site, the mass loss rates of three ice sheets, and the temporally correlated noise in the sea-level observations. Using only tide gauge observations limits our ability to separate estimates of sea level from estimates of the temporally correlated noise. This led us to modify the KS approach in two ways.

First, the state vector includes only an estimate of total sea level at every tide gauge site in addition to the desired mass loss rates. This yields the following state vector, xk, at every time step, k:

xk = [hk Bk]^T

where hk is a vector of sea level at the 622 tide gauge sites, and Bk is a vector containing the scalar weightings of 3 ice sheets and 18 mountain glacier regions (see below), as well as a uniform component that accounts for global mean thermal expansion and any additional mass contributions from smaller mountain glaciers.

Second, while in ref. 9 the observation model consisted of the sum of the estimated sea level, correlated noise, and white noise, here, the observation model consists only of the estimated sea level and white noise at each tide gauge site. Temporal correlations due to ocean dynamics are now modelled by the annual, spatial, CMIP5 ocean model fields (see below for a more detailed description of the CMIP5 model fields).

Sea level is modelled as the Euler integration of the contributions from melt sources, Bk-1y, (with y being the matrix of sea-level fingerprints associated with rapid land-ice mass loss), the ongoing rate of sea-level change due to GIA, G, and the rate of change of sea level due to ocean dynamics, Sdot_{k-1} , from the spatial fields in the CMIP5 model outputs:

hk = h_{k-1} + ∆t(B_{k-1} y + G + Sdot_{k-1}) + wh

where wh represents a zero-mean, white noise term associated with sea level.

The scalar weightings of the fingerprints are modelled as a random walk:

Bk = B_{k-1} + wB

where wB represents a zero-mean, white noise term associated with the melt contributions. The forward filtering pass of the Kalman smoother follows the steps outlined in ref. 9. A final departure from the methodology presented in ref. 9 is that we implemented a three-pass fixed-interval smoother33 in place of a Rauch-Tung-Stiebel two-pass smoother [34].

Gaussian process regression

The GPR approach, in contrast, models sea level as a multivariate Gaussian field defined by spatiotemporal mean and covariance functions that describe the underlying processes responsible for sea-level variability. Specifically, Gaussian process priors describing the contributions from land ice, GIA, and ocean models are conditioned simultaneously upon the available observations to produce the conditional, posterior distribution of sea level at decadal intervals throughout the twentieth century. In contrast to the KS, the GPR approach directly estimates the intertemporal covariance of the posterior; the associated computational demands require the use of decadal rather than annual means. Rather than being based upon discrete GIA and GCM models as in the KS approach, the GPR approach employs Gaussian process priors for the GIA and ocean dynamics contributions that are estimated, respectively, from the 161 GIA model predictions and 6 GCM outputs (see below). The distribution describing each land ice mass contribution is modelled assuming a prior spatio-temporal covariance, with the temporal component estimated from previous, non-sea level based estimates of land ice melt and the spatial component from the sea-level fingerprints associated with the melt source.

We model decadal-average sea level as a spatiotemporal field:

f(x,t) = fGIA(x,t) + fM(x,t) + fLSL(x,t)

where fGIA, fM, and fLSL are respectively the components of sea level due to ongoing GIA, land ice mass loss, and local effects associated with ocean dynamics, tectonics and other non-climatic factors, each as a function of location, x, and time, t. Each sea-level component is modelled as a Gaussian process with a prior mean function, μi(x,t), and covariance function, Ki(x,t,x′,t′).

The total field can be partitioned into observed sites, f1, and unobserved sites, f2, and subsequently written as a joint, multivariate distribution, such that:

[f1, f2] ~ N([μ1, μ2],[K11 K12, K^T12 K22])

Observations, y, are modelled as the underlying sea-level field with additive white noise characterized by zero mean and a covariance Σp, such that the joint distribution becomes:

[y, f2] ~ N([μ1, f2],[K11+Σp K12, K^T12 K22])

Using standard statistical results (see, for example, ref. 35), the posterior mean and covariance, f2 and V2, of the unobserved field conditioned upon the observations are:

f2 = f2 + K^T12 [K11+Σp]^{-1} y

and

V2 = K22 - K^T12 [K11+Σp]^{-1} K12

To estimate the underlying constituents of the total sea-level field, the prior mean and covariance of the unobserved field (that is, μ2, K12, K22) are set to the distribution of the desired quantity alone. For example, setting μ2, K12, and K22 equal to μ2M, K12M, K22M, returns the posterior mean and covariance of sea-level change due to the melt contributions. Once all the underlying constituent sea-level fields are determined, the global mean of those components can be computed and added to estimate GMSL.

The elements of the prior covariance matrix of the melt contribution, KM, are defined as:

K^Mi,j = Σ^n_{a=1}(A^{M,L}_{i,j,a} + A^{M,RQ}_{i,j,a})(B^M_{i,j,α})

where the subscripts indicate the ith row and jth column element of the ath ice sheet or mountain glacier. The time dependence of the covariance matrix is taken to be the sum of a linear component, AM,L, which accounts for secular changes in the melt contributions, and a rational quadratic term, AM,RQ, that represents a smoothly-varying function of variability:

A^{M,L}(tq,tp) = k1 tq tp

and

A^{M,RQ}(tq,tp) = k2 (1 + ∆t^2q,p / 2 α τs^2)^{-α}

Here, tq and tp represent the time at the qth and pth time step, Δtq,p represents the time difference between these steps, and k1, k2, α, and τs are hyperparameters that define the linear amplitude, rational quadratic amplitude, roughness, and characteristic timescale of the covariance functions35. To estimate the hyperparameters we adopt an empirical Bayesian approach where we compute the parameters that maximize the likelihood of reconstructed time series of previous mountain glacier estimates36 and ice sheet estimates37.

The spatial weighting of the prior covariance, BM, is computed as the outer product of the unique fingerprint associated with melt from the corresponding land-based ice source.

The prior spatiotemporal mean and covariance for the GIA contribution to sea-level change, μGIA(x,t) and KGIA(x,t), respectively, are taken as the sample mean and covariance of the 161 predictions of sea-level change described below.

The distribution of the contribution to sea-level changes from thermosteric and ocean dynamic effects is partially modelled as the sample mean and covariance of the CMIP5 model outputs32. However, since a small number of models are used to compute the distribution statistics, the estimated distribution may not be representative of the parent distribution. Consequently, we augment the sample covariance with a space-time separable covariance structure consisting of the product of two Matérn functions [35], C: one representing the temporal distribution and the other representing the spatial, such that the total prior covariance describing local sea-level change is given by:

KLSL(x,t) = KCMIP5(x,t) + C(t,ν1,τ) C(x,ν2,L)

where KCMIP5 is the sample covariance of the CMIP5 model outputs, ν1 and τ are the smoothness parameter and characteristic timescale of the temporal Matérn function, respectively, and ν2 and L are the smoothness parameter and characteristic length scale of the spatial Matérn function, respectively. For the exponents within the Matérn functions we follow ref. 10 and set the exponent on the spatial component to ν2 = 5/2 (reflecting a relatively smooth, twice-differentiable field) and the exponent on the temporal component to ν1 = 3/2 (reflecting a once-differentiable time series, in which rate is always defined but can change abruptly). As with the melt covariance hyperparameters, we use an empirical Bayesian approach to estimate the maximum-likelihood time and length scales of the Matérn functions to be 46 years and 90 km, respectively. Note that there is some trade-off between the Matérn exponent values and the hyperparameter characteristic scales: the selection of, say, a lower exponent (giving rise to a less smooth functional form) would result in a longer length scale.

In addition to capturing the inaccuracies of the ocean dynamics distribution, the Matérn functions also model local tectonic, geomorphological and other non-climatic contributions to local sea-level change. These hyperparameters, and a white-noise variance, are computed by finding the parameters that maximize the likelihood of the available tide gauge observations given the complete sea-level model.

Sea-level fingerprints

Extended Data Fig. 1a and b shows global maps of sea-level change, known as sea-level fingerprints, associated with rapid, uniform mass loss across the Greenland Ice Sheet (GIS) and the West Antarctic Ice Sheet (WAIS), respectively. The sea-level changes are normalized by the equivalent GMSL change. Both fingerprints are characterized by a large amplitude sea-level fall in the region adjacent to the melting ice sheet with a gradual rise in sea level moving away from the ice sheet. The computation of the fingerprints is based upon a gravitationally self-consistent sea-level theory that takes into account shoreline migration and changes in grounded, marine-based ice cover as well as the impact on sea level of perturbations in the Earth’s rotation axis38,39,40.

In addition to the GIS and WAIS, fingerprints were computed for the East Antarctic Ice Sheet (EAIS) and glaciers of Alaska, the Alps, Baffin Island, the Caucasus, Ellesmere Island, Franz Josef Land, High Mountain Asia, Altai, Iceland, Kamchatka, the low-latitude Andes, New Zealand, Novaya Zemlya, Patagonia, Scandinavia, Severnaya Zemlya, Svalbard, and Western Canada/US.

We also include a spatially uniform pattern to account for changes in GMSL due to land ice sources not included in the above set of glaciers. In the Kalman smoother, this uniform ‘fingerprint’ also captures changes in GMSL due to globally uniform thermal expansion and terrestrial water storage variations9.

GIA models

The first step when analysing tide gauge records is to correct for sea-level contributions due to the ongoing GIA of the Earth in response to the ice age cycles. Predictions of GIA are dependent on the geometry and deglaciation history of the Late Pleistocene ice sheets and the Earth’s viscoelastic structure. In this study, we computed 160 different GIA predictions distinguished on the basis of the adopted lower-mantle viscosity, upper-mantle viscosity, and thickness of a high-viscosity (effectively elastic) lithosphere. Additionally, we computed a GIA prediction using the VM2 viscosity profile41. These were combined with the ICE-5G (Ref. 41) global ice sheet reconstruction for the last glacial cycle. A detailed description of physical processes that contribute to the total GIA signal can be found in ref. 42.

We adopted values for the three rheological model parameters that encompass all recent estimates of the Earth’s structure. The lower-mantle viscosity was varied in the range (2–100) × 1021 Pa s, upper-mantle viscosity in the range (0.3–1) × 1021 Pa s, and lithospheric thickness in the range 72–150 km. Extended Data Fig. 2a and b shows the mean and standard deviation of the model predictions. The largest variance is seen in the region within the near field of the former ice sheets, including areas of ancient ice cover and the so-called peripheral bulges.

Ocean dynamics models

We treat the thermosteric and ocean dynamic contributions to sea level using the historical experiment output from 6 global climate models of the World Climate Research Programme’s (WCRP) Coupled Model Intercomparison Project phase 5 (CMIP5) data set32. Following ref. 9, the models we use are: bcc-csm1-1 from the Beijing Climate Center, CanESM2 from Environment Canada, the NOAA-GFDL model GFSL-ESM2M, the Institut Pierre Simone Laplace IPSL-CM5A-LR model, MRI-CGCM3 from the Japanese Meteorological Institute, and NorESM1-M from the Norwegian Climate Centre. For the KS methodology, we use the zero-mean spatial field ‘zos’ that is supplied by all the models. In the GPR, we add to ‘zos’ each model’s estimated globally averaged sea-level change due to thermal expansion: ‘zossga’.

While the CMIP5 model outputs are provided as global ocean grids, the field values at the specific locations of tide gauges are required, as input, to both the KS and GPR analyses. Where the tide gauges are coincident with model grid points, the associated value of the model output is used. Otherwise, an inverse distance weighting interpolation scheme is used to estimate the field at the desired location.

We examined three alternative interpolation schemes to assess the sensitivity of the KS GMSL estimate to this choice: (1) a nearest-neighbour approach, selecting the value on the CMIP5 grid that is closest to the tide gauge site; (2) a Delaunay interpolant, computing a linear interpolation between the irregularly spaced model cells along the coastlines; and (3) a Gaussian process (or simple kriging) methodology. For the Gaussian process interpolation, we employed a Gaussian process prior with a mean equal to the mean of the model grid values within a 200 km radius of the tide gauge location and a Matérn covariance function with smoothness parameter equal to 5/2. Since we are interested in the variability of the ocean models immediately surrounding each tide gauge site, the length scale of the Matérn covariance function was set to 1° (∼110 km). Neither the nearest-neighbour approach nor the Delaunay interpolated altered the estimate of the GMSL rate over the time period 1901–90. The Gaussian process interpolation scheme changed the GMSL estimate by less than 2%, significantly smaller than the estimated ±0.2 mm yr−1 90% CI on the estimate.

Computation of GMSL rates and accelerations

The mean and uncertainty of GMSL rates are estimated using a generalized least squares regression of a linear trend to the reconstructed GMSL time series. While the GPR methodology outputs a full temporal covariance matrix, the KS methodology does not. For this purpose, we adopt a temporal covariance matrix Σ with elements having the form:

Σ_{i,j} = σi σj exp(-(tj - ti)/τ

where σi and σj are the instantaneous uncertainties in GMSL at time i and j, respectively, derived in the multi-model KS analysis. To estimate the decorrelation timescale, τ, we examined the annual PSMSL tide gauge data and computed the mean temporal correlation coefficient across all tide gauges. This coefficient approaches zero after 2 years, and we set τ to 3 years. Estimates of acceleration in GMSL cited in the main text for the two probabilistic analyses are computed using a generalized least squares fit of a quadratic through the associated GMSL time series. Estimates of acceleration for the Church and White4 and Jevrejeva et al.3 time series listed in Fig. 2 are based on a weighted least squares regression through the published time series (see figure legend).

Analysis of bias introduced by using a subset of tide gauges

We used a bootstrapping technique to assess the potential biases introduced in estimates of GMSL rates when only a subset of tide gauge records is used. We randomly sampled our global sea-level reconstruction based on the Kalman smoother at a specific number of tide gauge sites in the database of 622 sites, computed the associated GMSL curve by binning the sites into 12 regions and averaging the result, and then used this curve to determine the rate of sea-level change over the time period 1901–90. The time series of the sea-level reconstruction at any given tide gauge site were sampled to match any missing data at that site in the PSMSL database. We repeated the analysis 100 times for subsets ranging in size from 25 to 600 sites. The mean sea-level rate we computed in this exercise and its associated uncertainty are shown in Extended Data Fig. 3 as a function of the number of sites. The horizontal blue line and shading is the mean rate of sea-level rise from 1901 to 1990, and its associated uncertainty, respectively, obtained from the KS-derived time series (1.2 ± 0.2 mm yr−1; Figs 2, 3).

The mean sea-level rate obtained from this analysis asymptotes towards its final value and the spread in rates decreases monotonically as the number of tide gauges used in the analysis increases. The asymptote lies ∼0.4 mm yr−1 above the KS estimate, which is consistent with the difference between the KS and ‘KS PSMSL sampling’ rate estimates for 1901–90 shown in Fig. 3. This result suggests that the combined effects of data sparsity and missing data introduce an upward bias into the estimate of GMSL. This bias is reduced in the KS (and GPR) methodologies because these techniques extract global information by using the observations, together with model-based geometries (or covariances) associated with the underlying contributions, to estimate (and sum) these contributions.

Analysis of bias introduced by an EOF analysis of altimetry records

To compare our results with EOF-based reconstructions of sea level4,43, we computed the GMSL time series following the approach adopted by Church and White4, but replacing altimetry and tide gauge observations with our KS reconstruction. The EOFs were computed using the KS sea-level reconstruction from 1993 to 2010, limited to the latitudinal observation range of satellite altimetry (65° N to 65° S). As in ref. 43, a spatially uniform EOF was added to the basis set to account for changes in mean sea level within the altimetry data (here the KS reconstruction), while the weights of the EOFs were computed using the first differences of the KS reconstruction at the tide gauge locations (sampled to reflect missing data in the PSMSL database) in order to eliminate dependence on a consistent datum. The GMSL time series was computed using an area-weighted mean of the EOF-reconstruction. To compute the uncertainty in our estimated GMSL, we sampled our distribution for each KS reconstructed tide gauge 1,000 times and computed the corresponding EOF-derived GMSL time series. We used this distribution of GMSL curves with a generalized least squares regression to compute a trend and uncertainty. This analysis yielded a linear trend of 1.4 ± 0.4 mm yr−1, demonstrating the existence of a bias since the ‘true’ underlying reconstruction has a trend of 1.2 ± 0.2 mm yr−1 (see Fig. 3, ‘KS EOF’).

Inverted barometer correction

The results in the manuscript were obtained using tide gauge observations that were not corrected for the inverse barometer (IB) effect. Previous studies (for example, refs 44 and 45) have shown that the sea-level response to atmospheric pressure changes can be non-negligible on regional scales.

In order to investigate the potential effect that atmospheric pressure changes have on our probabilistic estimate of GMSL, we repeated the KS analysis on the full tide gauge data set after we corrected these records for the IB effect. Specifically, we used the HadSLP2 global reconstructed atmospheric pressure data set46 to compute the IB correction. We next applied the correction to the observations at the 622 tide gauge sites and then re-ran the KS analysis. The 1901–90 GMSL rate of change associated with this analysis is 1.2 ± 0.2 mm yr−1, consistent with the value cited in the main text. We conclude that while the IB effect can impact regional sea-level histories, it has a negligible effect on our probabilistic estimates of GMSL.

Optimality of the Kalman smoother

Local sea levels observed by tide gauges reveal significant interannual and decadal variability. This variability can lead to temporal correlation in the sea-level time series that needs to be considered if one seeks to obtain optimal estimates of the underlying GMSL contributions. In order to test the optimality of the Kalman smoother, we investigated the properties of the innovation sequence by computing the residuals between the observations and the KS model estimate of sea level at every tide gauge site. Since every KS estimate of sea level is accompanied by its associated uncertainty, we randomly sampled from each sea-level distribution to obtain 100 time series of residuals for every site. Following the optimality test described in ref. 47, we computed the mean AR(1) coefficient across the 100 samples at each tide gauge site. An optimal Kalman smoother is characterized by a white noise innovation sequence. In practice, this means that, within uncertainty, the AR(1) coefficients of the innovation sequences will be close to zero. In the exercise above, we obtained a mean AR(1) coefficient of 0.2 ± 0.3 (90% confidence). This indicates that our innovation sequence is (within uncertainty) white noise and that the smoother is, or is close to, optimal.

Sensitivity of GMSL estimates to limitations of the CMIP5 climate simulations

The presence of unmodelled ocean dynamics can also affect the smoother performance. As described above, the limitations of the CMIP5 simulations as models for the true dynamic variability of the oceans is addressed in the GPR analysis by augmenting the covariance computed from the climate runs with two additional terms: a covariance modelled with two Matérn functions, and a white noise variance.

To assess the sensitivity of the KS analysis to unmodelled ocean dynamics, we examined its response to (1) a known synthetic ocean dynamic signal and (2) the inclusion of the dynamic response to freshwater hosing of the North Atlantic.

We used the mean KS estimates of the ice sheet melt rates and uniform sea-level contribution, as well as the multi-model estimate of the GIA contribution, to construct synthetic sea-level observations at the 622 tide gauge sites. We then added the dynamic sea-level change associated with one of the six CMIP5 climate models and ran the multi-model KS using the five remaining climate models to obtain an estimate of the GMSL rate. We repeated this analysis for each of the CMIP5 simulations. By not including the climate model used in constructing the synthetics in the multi-model component of the KS methodology, we tested the ability of the smoother to account for unmodelled dynamics. The 1901–90 GMSL rates determined from the complete set of 6 analyses ranged from 1.1 to 1.3 mm yr−1. Five of these analyses yielded a 90% CI of 0.2 mm yr−1, while the sixth yielded a 90% CI of 0.3 mm yr−1. These values are consistent with the results for the KS analysis cited in the main text (1.2 ± 0.2 mm yr−1).

To assess the sensitivity of our GMSL results to ocean dynamic effects due to freshwater input (‘hosing’) from GIS melt, we used the results of a previous study48 to investigate the dynamic sea-level signal arising from North Atlantic freshwater ‘hosing’ simulations. Specifically, we computed the difference between the results of the 0.1 Sv hosing run and the control (no-hosing) simulation described in ref. 48 and scaled this difference by 0.05 to approximate a synthetic dynamic signal for a GIS melt rate equivalent to 0.5 mm yr−1 GMSL rise over the twentieth century. After subtracting a uniform 0.5 mm yr−1 from the spatial pattern, we calculated time series of this signal at all 622 tide gauge sites, added these to the observed record, and repeated the KS analysis. The presence of these unmodelled dynamics has negligible effect on our estimate of GMSL. The 1901–90 rate estimated in the above test agrees with the value presented in the manuscript (1.2 ± 0.2 mm yr−1).

While the above sensitivity tests indicate that the probabilistic analyses have quantified, with reasonable accuracy, the impact of uncertainties in CMIP5 models of ocean dynamic variability, improving such models is an important requirement in any effort to further refine estimates of GMSL rates.

Kalman smoother reconstruction of sea level at sites with no observations

To investigate how well the KS is able to estimate sea level at sites without observations, we ran the Kalman smoother using data from 450 randomly chosen tide gauge sites and estimated the sea level at the remaining 172 sites. Extended Data Fig. 4 shows the GMSL time series estimated in this new analysis as well as a comparison of the estimated and observed sea level at a representative subset of 5 of these 172 sites (the remaining sites show similar fits). We calculate a 1901–90 GMSL rate of 1.2 ± 0.2 mm yr−1, consistent with the results presented in the manuscript when all 622 tide gauge sites are used in the analysis. The consistency between the estimated and observed values at the 172 tide gauge sites also indicates that limitations of the CMIP5 simulations in modelling ocean dynamics are not degrading the ability to predict sea-level trends at sites without observations. More generally, the analysis demonstrates the power of the KS method in reconstructing sea level when the method is applied with physics-based and model-derived geometries of the underlying physical processes.

Extended data figures and tables

Extended Data Figure 1: Illustrative sea-level fingerprints. a, b, Normalized sea-level changes due to rapid melting of the Greenland Ice Sheet (a) and the West Antarctic Ice Sheet (b). The variable ‘normalized sea-level change’ on the colour scale is formally dimensionless, but may be interpreted as having the unit of metres of sea-level change per metre of the equivalent GMSL change associated with the melt event.

Extended Data Figure 2: The present-day rate of change of sea level in mm yr−1 due to GIA for a suite of Earth models. a, b, Mean sea-level change (a) and standard deviation (b) computed from the output of 161 GIA model simulations (see text). In both frames, the colour scale saturates in the near field, which includes areas of post-glacial rebound and peripheral subsidence.

Extended Data Figure 3: Bootstrapping analysis of GMSL rate for 1901–90 obtained by sampling the global reconstruction of sea level. Data points show the mean computed from a bootstrapping analysis of the 1901–90 GMSL rate as a function of the number of geographic sites used in the analysis (ranging from 25 to 600). Error bars, ±1s.d. Sites are obtained by randomly sampling the global KS reconstruction at a subset of tide gauge sites and introducing data gaps that are consistent with those that exist in the PSMSL database15. The analysis was repeated 100 times for each choice of the number of sites. Also shown (horizontal blue line and shading) is the 1901–90 rate and its 90% CI computed from the KS GMSL curve in Fig. 2 (1.2 ± 0.2 mm yr−1; Figs 2 and 3).

Extended Data Figure 4: Results of the KS analysis performed using a random subset of 450 tide gauges. a, KS-estimated GMSL curve derived using a subset of 450 of the 622 tide gauge records discussed in the main text (blue line) and the reconstruction of Church and White4 (magenta line) and Jevrejeva et al.3 (red line). The shaded regions represent the 1σ certainty range. Panels b–f show the KS reconstructions (black lines) at a representative set of 5 of the 122 sites that were not used in the estimation procedure. The observations are shown in red.