# Estimating the sources of global sea level rise with data assimilation techniques

See allHide authors and affiliations

Edited by Jeffrey Shaman, Columbia University, New York, NY, and accepted by the Editorial Board March 22, 2012 (received for review October 26, 2011)

## Abstract

A rapidly melting ice sheet produces a distinctive geometry, or fingerprint, of sea level (SL) change. Thus, a network of SL observations may, in principle, be used to infer sources of meltwater flux. We outline a formalism, based on a modified Kalman smoother, for using tide gauge observations to estimate the individual sources of global SL change. We also report on a series of detection experiments based on synthetic SL data that explore the feasibility of extracting source information from SL records. The Kalman smoother technique iteratively calculates the maximum-likelihood estimate of Greenland ice sheet (GIS) and West Antarctic ice sheet (WAIS) melt at each time step, and it accommodates data gaps while also permitting the estimation of nonlinear trends. Our synthetic tests indicate that when all tide gauge records are used in the analysis, it should be possible to estimate GIS and WAIS melt rates greater than ∼0.3 and ∼0.4 mm of equivalent eustatic sea level rise per year, respectively. We have also implemented a multimodel Kalman filter that allows us to account rigorously for additional contributions to SL changes and their associated uncertainty. The multimodel filter uses 72 glacial isostatic adjustment models and 3 ocean dynamic models to estimate the most likely models for these processes given the synthetic observations. We conclude that our modified Kalman smoother procedure provides a powerful method for inferring melt rates in a warming world.

Sea level (SL) changes occur on timescales ranging from seconds to millions of years. At the longest of these scales, over millions of years, changes are driven by processes such as continental collision, sedimentation, variations in the mean spreading rate of the sea floor, and dynamic topography fluctuations due to mantle convection (1, 2). At intermediate timescales, from hundreds to tens of thousands of years, glacial isostatic adjustment (GIA) in response to the Late Pleistocene glaciation yields a global pattern of SL change that varies systematically from the near field to the far field of ancient ice cover (3). Finally, at shorter timescales, a suite of processes contributes to the complex spatiotemporal pattern of modern SL change (4). These latter processes can be separated into “dynamic effects,” which include ocean tides, ocean–atmosphere interactions, salinity and temperature variations, and ocean circulation, and “static equilibrium effects,” which are the gravitational, deformational, and rotational signatures of mass flux from polar ice sheets and mountain glaciers (5). In this paper we explore the connection between the processes and sources that give rise to SL change and the geographic variability inherent in the associated observational constraints to investigate whether—and to what accuracy—the signal from recent ice mass fluxes may be estimated from modern tide gauge records of SL.

Twentieth century rates of SL change have been estimated using data from tide gauges: devices that measure relative sea level (i.e., the difference in height between the sea surface and the sea floor) at a specific geographic location. The Permanent Service for Mean Sea Level (PSMSL) compiles tide gauge records for >1,000 sites distributed around the globe. Estimates of globally averaged 20th century SL change rates based on various subsets of these tide gauge records (6⇓⇓⇓⇓–11) fall within the range of 1–2 mm/y. In the last decade of the 20th century and in the current century, tide gauges have been complemented by satellite altimetry measurements of sea surface height and ocean float measurements of temperature and salinity (12). These records will be critical for understanding recent and future SL rise, but they are too short to inform our understanding of 20th century sea levels.

Estimating mean global SL alone ignores the information content inherent to geographic variability. As an example, it is now well understood that the rapid melting of any individual grounded ice reservoir gives rise to a distinct and highly nonuniform pattern—or fingerprint—of static equilibrium sea level change (13⇓⇓⇓–17). In particular, in the near field of a melting ice sheet (within ∼2,000 km of the margin), SL will fall due to both crustal uplift and the reduction of the gravitational pull on the ocean from the ice sheet. This fall can be an order of magnitude greater than the equivalent rise in mean global SL associated with the meltwater addition to the ocean. In contrast, in the far field of the ice sheet, SL will rise with (generally) greater amplitude as the distance from the ice sheet increases; this rise can exceed the global mean value by ∼30%. These fingerprints provide a framework for moving beyond inferences of globally averaged SL rise to estimate the contributions from individual meltwater sources. As an example, Mitrovica et al. (16) argued that the lower rates of SL change inferred from tide gauge records at European sites relative to the global average were consistent with 20th century melting from Greenland at a level equivalent to a eustatic SL rise of ∼0.5 mm/y.

Efforts to fingerprint modern records of SL change face potential contamination from dynamic processes, and thus robust estimates of ice mass flux must address the extent to which this contamination can be accounted for. Dynamic SL processes can be subdivided into three categories, due respectively to (*i*) internal, unforced dynamics (e.g., ocean tides and some air–sea interactions), (*ii*) atmospheric and ocean warming, and (*iii*) freshwater addition.

In this study, we outline and evaluate a statistical formalism for estimating global sea level change and the individual contributions to this change in the presence of the three dynamic SL effects listed above. Most previous work has estimated mean global SL change by computing the linear least-squares fit to an assumed constant, long-term trend in select tide gauge records (8, 10). In contrast, our formulation makes no assumption in regard to the stationarity of the trends. This generalization is accomplished by using a Kalman filter (18) and smoother (19) to recursively compute the maximum-likelihood estimate of the total SL change as well as an estimate of its associated uncertainty.

There are several advantages in using this approach. Due to the recursive form of the algorithm, the Kalman smoother naturally accommodates measurements with data gaps. This is an important feature, given that tide gauge records extend back into the 19th century but are characterized by pervasive and extensive data gaps. Accommodating data gaps permits the use of a larger selection of data than would be available when selecting and truncating records on the basis of, among other criteria, record continuity. The Kalman smoother assimilates observations when they are available, and in the absence of observations the smoother relies on the model dynamics. The recursive nature of the estimator approach is also useful because it reduces the storage requirements necessary when a large volume of observations are available. Additionally, although in essence a least-squares estimator, the Kalman smoother permits the determination of time-varying nonlinear trends in the measurement data that are overlooked using a simple linear least-squares approach. Finally, as our assumed model dynamics are linear, the Kalman smoother computes the optimal estimate of the state at every time step after convergence. We note that all of these advantages hold when the Kalman smoother is applied in applications that move beyond the determination of the total SL change to estimate time-varying contributions from individual sources of meltwater, such as the Greenland ice sheet (GIS) and the West Antarctic ice sheet (WAIS). Finally, we extend the Kalman smoother by including a multimodel filter (20), which provides a straightforward methodology for including uncertain contributions to the SL changes from sources other than meltwater, including ongoing GIA and thermosteric effects.

This study has two overarching goals. The first is to provide a rigorous and computationally efficient framework for analyzing the growing and disparate database of SL records. The second is to assess the accuracy with which individual sources of meltwater are detectable in synthetic SL records that have the same temporal and geographic data extent as presently available SL observations and are also characterized by realistic levels of dynamic variability. Our analysis is motivated by the following interrelated questions: For a suite of ice sheet melting scenarios, can we robustly identify individual ice sheets’ static equilibrium SL signals in the presence of dynamic variability over the last century given the current spatiotemporal sampling of the SL field? Are there specific areas in the world where observations are particularly effective in constraining these contributions? To what extent are these estimates degraded by uncertainties associated with other contributors to the SL change, including GIA and thermosteric signals? Is it possible to use 20th century SL measurements to estimate GIA and thermosteric signals?

This paper focuses on a description of the modified Kalman smoother procedure and a discussion of its performance when applied to synthetic data derived from tide gauge observations. *Statistical Formalism* describes in detail the formalism and the *Implementation* section describes the coverage and noise characteristics of the existing tide gauge dataset, the sea level fingerprint predictions that are adopted in the construction of the synthetics, the models used to account for dynamic SL and GIA, and the implementation of our statistical procedure as applied to tide gauge data. Finally, *Results* discusses the performance of the Kalman smoother in a large suite of synthetic tests.

## Statistical Formalism

### Kalman Smoother.

The smoothing algorithm implemented in this study can be divided into two main steps: a forward filtering pass and a backward smoothing pass. In the forward pass, the filter uses all past information to obtain an optimal estimate of the state vector, whereas the smoothing step runs backward in time, computing a smoothed state estimate using information from both the forward and the backward passes.

In the forward pass, the objective of the Kalman filter (21) is to use observations (or measurements), **z**_{k}, to estimate the state vector, **x**_{k}, at every time step (*k*). We model the observations with five different contributions to SL in the state vector: (*i*) SL changes due to melting of the GIS at a rate given by α_{k}; (*ii*) SL changes due WAIS melting at a rate of β_{k}; (*iii*) a geographically uniform SL rise, *U _{k}*, that accounts for any unmodeled sources of SL rise such as mountain glaciers; (

*iv*) SL trends due to GIA (

*G*); and (

_{i}*v*) the SL contribution,

*S*, due to ocean dynamics including thermal expansion.

_{i}In the filter, the observation model is defined as

where **H**_{k} is the observation matrix, **x**_{k} is the state vector, and ε_{k} is the white observation noise, arising primarily from ocean dynamics (*Sea Level Predictions* section), defined to have a mean of zero and a covariance **R**.

At every time step we define our state vector as

where *U _{k}*, α

_{k}, and β

_{k}are the scalar weightings of the first three contributors to sea level in our synthetic datasets, as listed above; the vector

**h**

_{k}is the SL at each site

*i*that includes contributions from GIS and WAIS melting, uniform SL rise, GIA, and ocean dynamics; and

**v**

_{k}accounts for unforced dynamic SL variability. We account for any temporal correlation in the dynamic variability by modeling

**v**

_{k}as an autoregressive (AR) process of order

*r*(

*SI Text*). The sum of

**h**

_{k}and

**v**

_{k}at every time step are the model’s estimate of the observations . Thus,

**h**

_{k}is a row vector whose length is equal to the number of tide gauge sites,

*n*

_{sites}, used in the analysis, and

**v**

_{k}is a row vector whose length is equal to

*n*·

_{sites}*r*.

At every time step we construct the observation matrix **H**_{k} that maps the model state vector into the observation space. **H**_{k}, defined as

generates the observation matrix

where **H _{v}** is the following <

*n*× (

_{k}*n*·

_{k}*r*) > matrix,

where each block has the dimension < *n _{k}* ×

*r*>. The sizes of both matrices depend on the number of observations,

*n*, available at time step

_{k}*k*.

Each time step, *k*, contains a forecast step and an analysis step. In the forecast step, predictions of a new estimate of the state vector, , and its covariance, , are based on the analysis estimates from the previous time step using the following equations:

The propagation from time *k* − 1 to time *k* depends on the state transition matrix **Φ**, the model noise **w**, defined to have zero mean and covariance **Q**, the input control parameter **u**, and the matrix **B** that maps **u** into the state space. The input control matrix is determined a priori and does not depend on the filter’s estimate of the state vector at any time step.

The state transition matrix, **Φ**, predicts how the state variables evolve between time steps,

which assumes that the globally uniform SL change and both ice sheet melt rates remain constant between time steps. The matrix **ρ** contains the AR coefficients that characterize the noise at every site (*SI Text*).

The input control parameter **u**_{k} governs how SL will vary due to modeled GIA and ocean dynamics:

Here **G** is a vector of GIA rates at every site *i*, whereas is the instantaneous rate of SL change due to ocean dynamics. **G** is set to be one GIA model (*Sea Level Predictions*) and is derived from one historical climate model (*Dynamic Sea Level*). Although one GIA model and one climate model must be chosen in the implementation of the filter, these models can be varied in repeated tests to examine the sensitivity of the filter on model choice. This procedure is discussed in detail in *Multimodel Kalman Filter*.

With the input control parameter introduced, we define the matrix **B** that relates **u**_{k} to the state space:

Substituting **Φ**, **u**_{k}, and **B** into Eq. **6** results in the following recursive equations for our five state variables:

The SL height, , thus depends on the SL at the previous time step, ; the estimated melt rates, , , and from the previous time step; and the control inputs **G** and . The predictions for the melt rates , , and are modeled as random processes with AR1 coefficients equal to 1. As previously noted, the temporally correlated measurement noise, , is also modeled as an AR process, of order *r*, with coefficients in the matrix **ρ**. We acknowledge that aspects of this state propagation model may be simplistic, particularly the time-invariant melt rates, but we rely on the available observations to correct the state propagation.

The model noise for all state variables is characterized by **w**, which has zero mean and covariance **Q**. Although **Q** can take on a time-varying structure, for simplicity we choose the following constant form:

The first upper diagonal is an *n*_{sites} × *n*_{sites} matrix corresponding to the noise in the modeled heights. The **V**_{h} structure is an exponential covariance matrix:

Here τ is the decorrelation length scale and **D** is the distance between tide gauge sites (21). The parameter σ_{h}^{2} gives the variance of each SL height, and σ^{2} is the variance of the three melt rates (which is assumed to be the same for all three). The final element of **Q** is the covariance of **w _{v}** that is associated with the correlated noise. This element of

**Q**is set to

**R**′, which contains

**R**(the covariance of the detrended tide gauge observations—see

*SI Text*for a description of

**R**). The elements of

**R**are assigned within

**R**′ to the positions corresponding to the first lagged correlated noise state for each site.

Following the forecast step is the analysis step when the observations are assimilated. The analysis state vector, , and the analysis error covariance, , can be updated using Eqs. **14** and **15**, respectively:

In Eq. **14**, the analysis state vector is a linear combination of the forecast estimate of the state vector, , and the weighted difference between the observations and the measurement prediction, . The weighting **K**_{k}, known as the Kalman gain matrix, minimizes the *posteriori* error covariance in Eq. **15** and is calculated following Eq. **16**:

If observations are unavailable, then , and .

Once the forward pass through the observations is complete, the backward smoothing pass is implemented. Whereas the Kalman filter pass uses all observations at and before the time at which the state vector is estimated, the smoothing pass runs backward in time, using intermediate results from the forward pass so that the final state estimate depends on all available observations over time. This type of smoothing algorithm is known as a fixed-interval Rauch–Tung–Striebel (RTS) two-pass smoother (19).

The smoothed estimate at the final time step *t*_{f} is initialized in the following manner:

is then recursively calculated using Eq. **18**,

where

The associated uncertainty on the smoothed estimate, , is also calculated on the backward pass:

It is important to emphasize that although the backward pass does not involve the processing of actual observations, it does use the full forward filtering solution. As we are estimating historical melt rates rather than predicting future ones, the smoothing pass ensures that at every point in time our estimated state vector and its associated uncertainty include all information available between the first and the final time steps.

### Multimodel Kalman Filter.

As evident in the description of the forward filter pass, the Kalman filter forecast step requires the construction of the input control parameter **u**_{k} (Eq. **6**), which depends on the choice of models both for GIA and the ocean dynamic signal. Although there is only one correct GIA model and one correct ocean dynamics model, we are, in practice, unsure which models those are. Filter convergence can occur in the case of a poorly modeled **u**_{k}; however, this situation will produce residuals between the state estimate and the observations that can grow over time (22). This problem can be addressed by either estimating the unknown model parameters by adding them to the state vector or extending our Kalman smoother to include and test multiple models (20). We implement a multimodel filter where we define each model, **μ**_{j}, to be a combination of one GIA model and one ocean dynamics model, resulting in a total of *N* different model pairings. The objective of the multimodel filter is to use the residuals between each filtered-model estimate and the observations to determine the most likely combination of GIA and ocean dynamics models.

For each model, **μ**_{j}, we define the likelihood of one model given the observations at time step *k* as the conditional probability density function,

i.e., is normally distributed with mean and covariance , where is defined as

All variables are consistent with their previous definitions and *n*_{states} is the size of the state vector. The calculation of *f* requires the forecasted estimate and covariance of the state vector because these estimates depend on the chosen GIA and ocean dynamics models. We compute the likelihood of each model at each time step, using the filtered state estimates and from Eqs. **6** and **7**. Calculating *f* also requires computing the inverse and determinant of the sometimes ill-conditioned matrix **C**. To estimate the inverse when the covariance matrix has a condition number exceeding a threshold value, we compute the eigenvalues of **C** and then regularize by retaining only the largest eigenvalues and eigenvectors that contain 99.5% of the variance of the covariance matrix. We then reconstruct **C** using the truncated eigenvectors and eigenvalues.

At each time step *k*, we calculate the probability, *p _{j}*, that the model,

**μ**

_{j}, is the correct model at that time:

The recursive computation of Eq. **23** requires initialization of the probabilities at the first time step. While prior knowledge of the likelihood of individual models can be used to initialize *p*, in the absence of this information we can consider all *N* models to be equally likely; i.e., *p _{j}*(0) = 1/

*N*. Although the likelihood and probabilities could be computed in one batch, we calculate them recursively to avoid large matrix inversions and storage. We then construct the model-conditioned estimate of the state vector, , and the associated error covariance using Eqs.

**24**and

**25**, respectively. These quantities are defined to be the weighted combination of each model’s smoothed state and covariance estimates with the weightings set to the model probability estimates at the final time step,

*p*(

*t*

_{f}):

## Implementation

### Datasets.

There are 1,166 tide gauges in the PSMSL revised local reference (RLR) database (www.psmsl.org). Although tide gauge records exist for sites in all continents (Fig. 1), there is a significant sampling bias, with the vast majority of tide gauges found along the coasts of North America and Eurasia. Examination of the time series at individual sites reveals that, in addition to spatial biases, the length of the datasets varies significantly throughout the database, with the longest records located in the Northern Hemisphere, particularly in the United Kingdom and Scandinavia (Fig. S1).

For the present study, we use mean annual tide gauge records to avoid seasonal variability. The largest interannual variability (Fig. 1) occurs along the northern coast of Russia, in closed basins (e.g., the Mediterranean and Black Seas), and in the western Pacific. We have chosen to use tide gauges that have at least 20 y of data in the last 50 y, leaving 596 records for analysis. We have binned the 596 tide gauges into the five regions illustrated in Fig. 1. Region 1 includes all tide gauges located poleward of 60°N, regions 2–4 move southward in 20° latitudinal bins, and region 5 includes all sites south of the equator.

### Dynamic Sea level.

To estimate long-term, warming-forced dynamic sea level trends, we use three climate models that participated in the World Climate Research Programme's (WCRP's) Coupled Model Intercomparison Project phase 3 (CMIP3) multi-model dataset (23): the Goddard Institute for Space Studies Model-E with the Russell Ocean Model (GISS-ER), the Model for Interdisciplinary Research on Climate Medium Resolution Model 3.2 (MIROC-3.2 medres), and the Meteorological Institute of the University of Bonn ECHAM4 HOPE-G Model (MIUB ECHO-G). The climate models are solely used to determine long-term trends in SL due to changing ocean and atmosphere dynamics; decadal-scale and interannual variability is included in the noise term.

Although the climate model outputs are provided on a global grid, we require coastal points for our tide gauge analysis. We estimate the model output at each site by using the model cell that is coincident with the tide gauge position if it exists. If the tide gauge location lies between multiple model grid points, we use an inverse distance interpolation on all available surrounding grid cells that are within a 1,000-km radius of the site. The 1,000-km radius is appropriate given that the correlation length scale of the tide gauge observations is ∼1,500 km (*SI Text*) and the southern hemisphere is data sparse, thereby requiring a large radius to ensure sufficient observations for interpolation. Despite the apparently large radius, the inverse distance relationship results in significantly higher weightings being given to sites close to the interpolated location.

### Sea Level Predictions.

The nonuniform patterns of SL change that arise due to rapidly melting ice sheets are commonly referred to as SL fingerprints. The fingerprints we adopt in this study are calculated by assuming geographically uniform ice melting across both GIS and WAIS (24). They are based on a gravitationally self-consistent SL theory (25) that incorporates shoreline migration due to local transgressions and regressions and changes in grounded, marine-based ice cover, as well as feedback into SL of contemporaneous perturbations in the Earth’s rotation vector (26). Our calculations adopt a 1D elastic Earth model with structure prescribed by the seismically inferred Preliminary Reference Earth Model (PREM) (27). This elastic (i.e., instantaneous) treatment of load-induced deformation and gravity changes is valid for melting events with timescales of several centuries or less (see Fig. S2 for GIS and WAIS SL fingerprints).

The effect of ongoing GIA due to the last ice age on relative SL observations is modeled using a more general version of the SL theory described above that incorporates viscoelastic adjustments. GIA predictions are sensitive to both the space–time geometry of the Late Pleistocene ice cover and the Earth’s viscoelastic structure. Most previous analyses of modern SL datasets were corrected using a GIA calculation based on a radial viscosity profile that increases only moderately with depth from the base of the lithosphere to values of ∼2–3 × 10^{21} Pa·s in the deep mantle (6⇓–8, 28, 29). However, Mitrovica and Davis (30) showed that relatively minor variations in the viscoelastic Earth structure map into uncertainties in estimates of globally averaged SL rise of order of 0.5 mm/y. Furthermore, an increase in the lower mantle viscosity to values of 5–10 × 10^{21} Pa·s, preferred on the basis of several analyses of GIA-related datasets (31, 32), significantly reduces systematic geographic trends in GIA-corrected tide gauge rates along the US East Coast (9).

In the formulation described above, uncertainties in the GIA signal are explored by performing the analysis using a suite of GIA models. We have run 72 GIA models, all of which are based on the global ice sheet reconstruction model ICE-5G (33). The models are distinguished on the basis of the choice for the lithospheric thickness and the viscosities in both the upper and the lower mantle. Specifically, we use lithosphere thicknesses in the range 72–125 km, upper mantle viscosities between 0.3 and 1 × 10^{21} Pa·s, and lower mantle viscosities between 2 and 20 × 10^{21} Pa·s. Illustrations of the SL prediction for one GIA model and the variance between the 72 models are presented in *SI Text* (Fig. S3).

### Synthetic Datasets.

To evaluate the performance of the Kalman smoother with known signals, we construct a suite of synthetic tide gauge data in the following manner. First, we remove long-term trends from the 596 RLR time series. For tide gauges containing <50 y of data we remove a linear trend, and for time series spanning >50 y we remove both a linear and a quadratic signal. Removing only a linear trend from the shorter time series was necessary to avoid the unrealistically large quadratic terms that occasionally resulted from estimating higher-order polynomial coefficients from insufficient data. The intent of the detrending is to remove all long-term trends associated with melting ice sheets, thermosteric and halosteric signals, GIA, and tectonics, leaving primarily decadal and interannual variability, which we treat as the observation noise. The detrended time series thus capture the shorter-term dynamic variability evident in the observations that is due to ocean circulation changes, long period tides, air–sea interactions, temperature and salinity variations, etc., as well as any regional variability. Second, we add SL fingerprints computed for various scenarios of GIS and WAIS melting. The SL fingerprints are normalized using the equivalent eustatic SL change and therefore the construction of the synthetic observations requires only that a weighting (in units of mm/y) be applied to each of the two fingerprints. A geographically uniform trend is also added to the synthetics to account for any additional, unmodeled, sources of SL rise such as mountain glacier and small ice sheet melting. Finally, we also add at each site SL contributions due to GIA and ocean dynamics, which include thermal expansion of the ocean.

The general expression for the synthetic tide gauge record at a given site, *i*, is given by

where *z _{i}* is the time series of relative SL at site

*i*;

*U*is a geographically uniform SL rise at time

_{k}*t*=

*k*;

*y*

_{i}_{,GIS}and

*y*

_{i}_{,WAIS}are the normalized GIS and WAIS SL fingerprints, sampled at site

*i*, with prescribed melt rates at time

*t*=

*k*given by α

_{k}and β

_{k}, respectively;

*G*is the sea level trend due to GIA at site

_{i}*i*; Δ

*t*is time step;

*S*is the SL contribution from ocean dynamics at each site

_{i}*i*; and

*v*is the site-specific dynamic variability (considered to be noise) in the observations, with covariance

_{i}**R**. Because the noise

*v*(

*t*) is derived from PSMSL data, the synthetic measurements contain the same irregular and pervasive data gaps as the PSMSL data. Although the longest tide gauge record begins in 1807 and all records end by the year 2008, resulting in the full dataset spanning 202 y, we begin our synthetic dataset in 1861 to coincide with the climate model runs. The GIA model added to the synthetics is arbitrarily chosen to have the Earth’s structure characterized by a lithospheric thickness of 125 km and upper and lower mantle viscosities of 0.5 × 10

^{21}Pa·s and 5 × 10

^{21}Pa·s, respectively. As discussed above, the dynamic SL,

*S*, is taken from 140 y of a randomly selected historical climate simulation. Specifically, we set the heights from 1861 to 2000 to be the heights from the MIROC3.2 medres historical climate model realization 1 and all dynamic heights after this time as constant at the year 2000 values. From each tide gauge site we use only the linear and quadratic trends in the model-derived dynamic SL time series, which ensures that all short term, interannual, and decadal variability in the synthetic time series arises from the detrended tide gauge observations.

_{i}With the synthetic datasets constructed, we can define the covariance matrices **R** and **Q** and initialize the state vector and its associated uncertainty (*SI Text*). We then recursively apply the algorithm to the full synthetic dataset, computed using one GIA–ocean dynamics model pairing, and run it forward and backward in time. The multimodel filter is implemented using the 72 GIA models (*Sea Level Predictions*) and 3 ocean dynamics models (Dynamic Sea Level), resulting in 216 GIA–ocean dynamic model pairings. For all tests, we calculate and plot the probabilities of each model as a function of time and present the final melt rate estimates and their associated uncertainties.

## Results

Results based on applying our Kalman smoother methodology to synthetic tide gauge data are presented in Figs. 2 and 3. The detection analysis is performed for each individual geographic region (Fig. 1) and for the combination of all regions. We begin by presenting the probabilities, calculated using Eq. **23**, of each of the 216 GIA–dynamic ocean pairing models as a function of time (Fig. 2). Each frame of Fig. 2 corresponds to results associated with a specific subset of tide gauges, beginning with region 1 in Fig. 2*A* and ending in Fig. 2*F*, where all 596 tide gauges are considered. In this scenario, the synthetics are computed using constant GIS and WAIS melt rates of 0.3 and 0.5 mm eustatic sea level (esl)/y, respectively, and a uniform rise in SL of 0.8 mm/y. The thick blue line in Fig. 2 *A–F* corresponds to the computed probability of the correct GIA and ocean dynamics pairing, and the gray lines correspond to the remaining 215 pairings. In addition to these probabilities, the melt estimates at the final time step along with their associated uncertainty (±1σ) are provided in Table 1.

To begin, consider the results for region 1 (Fig. 2*A* and Table 1). This region is in the near field of the GIS and is therefore a location where the GIS SL fingerprint has a large signal (and signal gradient) due to both the decrease in the gravitational attraction of the GIS on the surrounding ocean and the crustal uplift associated with the ice and water unloading. The probabilities indicate that although the correct GIA–dynamic ocean model pairing becomes more likely over time, this particular pairing has a peak probability of only 0.2 by the final time step. The spread in probabilities associated with the remaining 215 models suggests that this region would not be the optimal subset to adopt in a detection effort. Indeed, although the WAIS melt rate is inferred with high accuracy, it has very low precision. The opposite it true for the detection of the GIS melt rate. The GIS melt rate is inferred with greater precision than the WAIS melt rate (±0.33 compared with ±0.79) because the latter fingerprint is nearly uniform within region 1. This geometry makes it difficult to detect the signal in the presence of dynamic variability or to separate it from the geographically uniform component (estimated uncertainty of ±0.92) of the SL synthetics. Detection of the two fingerprints is clearly challenging in this region, even in the case (as in this study) where local dynamical effects associated with the influx of fresh water are not accounted for.

There is no improvement in the accuracy of the GIS detection when data from region 2 are considered (Table 1). In fact, whereas the precision remains almost constant between regions 1 and 2, the estimated GIS melt rate moves farther from the true value. The significant uncertainty in the WAIS result is, once again, due to the difficulty in separating the WAIS fingerprint from a globally uniform SL trend. There is large GIA model variability in region 2 (Fig. S3). This variability makes the analysis of SL data in this region particularly sensitive to the GIA signal, and therefore data from this region provide a potentially strong constraint in estimating the GIA signal. We can explore this in more detail by considering the results from the multimodel filter (Fig. 2*B*). Examination of Fig. 2*B* reveals that whereas the correct model combination (thick blue line) does not emerge as the most probable, the two most likely model combinations pair the correct ocean dynamics model with GIA models having the correct lithospheric thickness and upper mantle viscosity but incorrect lower mantle viscosity (10 × 10^{21} Pa·s and 8 × 10^{21} Pa·s for the most and second most probable models, respectively, instead of the correct value of 5 × 10^{21} Pa·s) (thin red lines). The filter’s bias against estimating the lower mantle viscosity of the GIA model may be due to the fact that region 2 includes a large number of tide gauges situated close to the location of the former Fennoscandian ice sheet (Fig. 1). The depth of the postglacial adjustments in the near field of an ice sheet depends on the size and weight of the ice sheet. Compared with the large Laurentide ice sheet, the solid earth deformation associated with the smaller Fennoscandian ice sheet was relatively shallow (limited to the top half of the mantle); therefore, GIA corrections in this region will be more sensitive to variations in the viscosity of the upper mantle than the lower. It is important to note that these sensitivities are highly dependent on the GIA model used in the synthetics and also on the binning of tide gauges.

As we consider tide gauge data in the vicinity of the equator (regions 3 and 4), it becomes increasingly difficult to distinguish between the SL fingerprints of GIS and WAIS melting and the uniform SL change, making detection of the fingerprints using data from this region alone highly problematic. This difficulty is reflected in the large errors in the inferred melt rates (Table 1). In addition, it also becomes difficult to distinguish between the 216 GIA–dynamic ocean model pairings using this data subset; indeed the probabilities of the model parings remain nearly constant over time (Fig. 2 *C* and *D*). As we move farther south we find that detection of WAIS melting has a smaller estimated error in region 5 than in any other region (Table 1). This difference is due to the large gradient in the static WAIS fingerprint that this region samples. However, the computed probabilities in this region (Fig. 2*E*) indicate that we cannot distinguish among the 216 GIA–dynamic ocean model pairings.

When all tide gauges are used in the analysis, the increase in precision of both GIS and WAIS melt rates is significant (Table 1). Because the Kalman smoother computes the state estimate by weighting the variance from the least-squares estimator, the best detection occurs when all available data are included in the analysis. Using the full, global dataset limits contamination associated with long-term regional periodic signals such as the El Niño-Southern Oscillation and also makes it possible to constrain the GIA and ocean dynamics signal (Fig. 2*F*). In Fig. 2*F*, the correct model emerges out of the ensemble with a final probability value of ∼0.75. The second most probable model pairs the correct ocean dynamics model with a GIA model that has the correct lithosphere thickness and upper mantle viscosity, but an incorrect lower mantle viscosity (8 × 10^{21} Pa·s instead of the correct value of 5 × 10^{21} Pa·s).

To illustrate the time-evolving performance of the multimodel Kalman smoother, we present in Fig. 3 the results for one specific melting scenario. Fig. 3 shows the melt rates (black lines) and error (±1σ, gray shading) estimates as a function of time.

To test the filter’s ability to track nonlinear trends, we carried out analyses using synthetics that include accelerations in both the GIS and the WAIS melt rates (Fig. 3 *A–C*). The initial GIS and WAIS melt rates are both set to 0. In 1910, the GIS melt rates increase to 0.1 mm esl/y. This increase is followed 50 y later by an acceleration of 0.009 mm esl/y^{2}. For the WAIS we impose two jumps in the melt rate: one in 1920 and one in 1935, to values 0.3 mm esl/y and 1 mm esl/y, respectively. We also linearly increase the uniform SL rise from an initial value of 0 mm/y to a final value of 1 mm/y. The synthetics computed using these accelerations are not meant to be a recreation of historical tide gauge records, but rather serve as a test of our methodology.

In Fig. 3 *A–C*, we note that the estimated errors are nearly constant over time. Note that the results in Fig. 3 are the multimodel smoothed estimates and uncertainties (Eqs. 25 and 26), which include the full forward and backward solutions. Therefore, the estimates at each time step include all assimilated observations from 1861 until 2008. In this case, we found that the Kalman smoother was able to track the accelerating rates of all three contributors to SL. Note that the estimate of the WAIS melt rate smooths gradients in the melt rates adopted in the synthetics; once again this result is because the Kalman smoother includes both a forward and a backward component (Fig. 3*B*). Our estimated uncertainties (±1σ) for the GIS and WAIS melt rates suggest that historical tide gauge records have information sufficient to robustly infer 20th century GIS and WAIS melt rates exceeding ∼0.3 and ∼0.4 mm esl/y, respectively.

## Conclusions

Most past studies of recent SL change have focused on determining mean SL trends, generally disregarding the distinct geographic patterns in SL that rapidly melting ice sheets produce. Our tide gauge analysis has provided a series of detection thresholds for inferring individual contributions of the GIS and the WAIS to SL change during the last century.

We have implemented a multimodel Kalman smoother data assimilation technique that iteratively performs a least-squares analysis whenever observations are available and in the absence of observations relies on the model dynamics to compute the best estimate of the state variables. In this study, we apply the technique to estimate three contributions to synthetic SL trends: melting from the Greenland and West Antarctic ice sheets and a globally uniform rise representing unmodeled components. The Kalman smoother is well suited for this problem because it accommodates data gaps and generates minimum mean square estimates of the SL changes. Also, in contrast to ordinary least-squares techniques, the use of the Kalman smoother permits the estimate of a nonstationary signal, which is clearly advantageous for efforts to detect changes in ice sheet melt rates. Because the Kalman smoother is a recursive technique, it also avoids storage constraints when dealing with large datasets. Additionally, the multimodel filter component of our formalism is not only able to account for uncertainties in SL signals due to GIA and ocean dynamics, but also provides the possibility of using SL observations to constrain models of GIA and ocean dynamics probabilistically.

In our tide gauge analysis we used 596 detrended tide gauge time series from the PSMSL database to construct a series of synthetic datasets for combined GIS and WAIS melt scenarios. Constructing our synthetic datasets using detrended tide gauge data ensures that our synthetics have the same data gaps as the historical tide gauge records. Having applied the Kalman smoother to subsets of tide gauges as well as the full set of synthetics, we have determined that when all available tide gauges are included in the analysis, it is possible, in principle, to rigorously estimate stationary and accelerating GIS and WAIS melt rates with 1σ uncertainty values of ±0.3 and ±0.4 mm esl/y, respectively. The size of these uncertainty values depends on the adopted model error covariance matrix **Q** and the observation noise covariance matrix **R**. Initializing these matrices establishes the trade-off between the relative weighting of a priori model constraints and the observations in the inference procedure. Therefore, additional tuning of these covariances (i.e., of the trade-off) may be necessary to optimize filter performance.

The Kalman smoother explicitly includes information on dynamic SL and GIA, thereby providing more robust estimates of melt rates than previous studies. The results from the multimodel filter indicate, perhaps not surprisingly, that data from some regions of the ocean will be particularly useful for constraining models of GIA. Although we somewhat arbitrarily divided the 596 tide gauges into latitudinal regions, more careful binning could be used to select tide gauges that have the greatest sensitivities to the Earth structure governing GIA. Constraining the dynamic SL model was also possible, particularly when the full set of tide gauges was used. Expanding this suite of dynamic ocean models to include a broader group of climate models will be important in the application of this methodology on 20th century tide gauge observations. In addition, we chose to implement a fixed multimodel filter that assumes that there is a single correct model for GIA and ocean dynamics over the full time span of the data. This assumption is certainly valid when considering GIA because the Earth’s structure will not change over centennial timescales. However, it may not hold for the ocean dynamics model when the methodology is applied to historical tide gauge records; that is, one of the suite of models for ocean dynamics may best capture the dynamics over one portion of the time history, whereas a second model may be most accurate in a different time window. One option is to include a dynamic multimodel filter, which would allow for the preferred ocean model to vary continuously over time.

In the methodology presented here we have included a uniform component of SL, which we attribute to unmodeled sources of SL change including mountain glaciers. In reality, each melting glacier produces a unique SL fingerprint. The Kalman smoother could easily be extended to estimate these additional source contributions; however, our ability to isolate these individual contributions from a uniform SL rise will be limited unless observations exist in the near field of the glaciers. Where such observations do exist, our methodology should be extended to incorporate the associated SL fingerprint.

This study has focused on the application of the multimodel Kalman smoother to tide gauge data; however, it is also possible to apply our statistical formalism to analyze satellite altimetry records of sea surface height variations. Although the current time span of altimetry data, ∼20 y, is a limiting factor, the near-global coverage of the dataset can provide potentially powerful constraints on SL fingerprints. The analysis of historical SL will ultimately combine the sparse tide gauge observations with the shorter time window, but nearly global, altimetry observations.

## Acknowledgments

Tide gauge data were provided by the Permanent Service for Mean Sea Level (www.psmsl.org). The code for the regularized expectation-maximization technique was provided by Tapio Schneider, and the code for computing the autoregressive coefficients was provided by Alois Schlögl. We thank the modeling groups, the Program for Climate Model Diagnosis and Intercomparison, and the World Climate Research Programme's (WCRP’s) Working Group on Coupled Modeling for their roles in making available the WCRP Coupled Model Intercomparison Project phase 3 (CMIP3) multimodel dataset. Support for this dataset is provided by the Office of Science, US Department of Energy. We also thank the anonymous reviewers for their helpful comments. R.E.K. was supported in part by an appointment to the US Department of Energy American Association for the Advancement of Science Fellowship Program administered by Oak Ridge Institute for Science and Education. J.X.M. acknowledges funding from Harvard University and the Canadian Institute for Advanced Research.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: jxm{at}eps.harvard.edu.

Author contributions: C.C.H., E.M., R.E.K., and J.X.M. designed research; C.C.H. and E.M. performed research; C.C.H., E.M., and R.E.K. contributed new reagents/analytic tools; C.C.H. and E.M. analyzed data; and C.C.H., E.M., R.E.K., and J.X.M. wrote the paper.

The authors declare no conflict of interest.

This paper results from the Arthur M. Sackler Colloquium of the National Academy of Sciences, “Fostering Advances in Interdisciplinary Climate Science” held March 31–April 2, 2011, at the AAAS Auditorium in Washington, DC. The complete program and audio files of most presentations are available on the NAS Web site at www.nasonline.org/climate-science.html.

This article is a PNAS Direct Submission. J.S. is a guest editor invited by the Editorial Board.

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

## References

- ↵
- ↵
- Müller RD,
- Sdrolias M,
- Gaina C,
- Steinberger B,
- Heine C

- ↵
- ↵
- Cazenave A,
- Nerem RS

- ↵
- ↵
- Peltier WR,
- Tushingham AM

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Woodward R

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Athans M,
- Chang CB

*Adaptive Estimation and Parameter Identification Using Multiple Model Estimation Algorithm*(MIT Lincoln Lab, Lexington, MA), Tech Rep 1976-28. - ↵
- Sherman M

- ↵
- Grewal MS,
- Andrews AP

- ↵
- ↵
- Mitrovica JX,
- et al.

- ↵
- Kendall RA,
- Mitrovica JX,
- Milne GA

- ↵
- Mitrovica JX,
- Wahr J,
- Matsuyama I,
- Paulson A

- ↵
- Dziewonski AM,
- Anderson DL

- ↵
- ↵
- Trupin A,
- Wahr J

- ↵
- ↵
- Lambeck K,
- Smither C,
- Johnston P

- ↵
- Mitrovica JX,
- Forte AM

- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Earth, Atmospheric, and Planetary Sciences