Two-dimensional reconstruction of past sea level (1950–2003) from tide gauge data and an Ocean General Circulation Model

Abstract. A two-dimensional reconstruction of past sea level is proposed at yearly interval over the period 1950–2003 using tide gauge records from 99 selected sites and 44-year long (1960–2003) 2°×2° sea level grids from the OPA/NEMO ocean general circulation model with data assimilation. We focus on the regional variability and do not attempt to compute the global mean trend. An Empirical Orthogonal Function decomposition of the reconstructed sea level grids over 1950–2003 displays leading modes that reflect two main components: (1) a long-term (multi-decadal), regionally variable signal and (2) an interannual, regionally variable signal dominated by the signature of El Nino-Southern Oscillation. Tests show that spatial trend patterns of the 54-year long reconstructed sea level significantly depend on the temporal length of the two-dimensional sea level signal used for the reconstruction (i.e., the length of the gridded OPA/NEMO sea level time series). On the other hand, interannual variability is well reconstructed, even when only ~10-years of model grids are used. The robustness of the results is assessed, leaving out successively each of the 99 tide gauges used for the reconstruction and comparing observed and reconstructed time series at the non considered tide gauge site. The reconstruction performs well at most tide gauges, especially at interannual frequency.


Introduction
Sea level is an indicator of climate change because it integrates the response of many components of the Earth sys-Correspondence to: W. Llovel (william.llovel@legos.obs-mip.fr) tem: the ocean and its interaction with the atmosphere, land ice, terrestrial waters. Even the solid Earth has some impact on sea level. Since the beginning of the 1990s, sea level is precisely measured by satellite altimetry systems (i.e., Topex/Poseidon, Jason-1 and now Jason-2) with global coverage and short revisit time. The satellite observations have revealed that sea level does not rise uniformly: some regions rise faster than the global mean; in some other regions sea level rise slower (Bindoff et al., 2007). It has been shown that the main cause of regional variability in rates of sea level change is non uniform thermal expansion of the oceans (Cabanes et al., 2001), although other processes may also give rise to regional sea level trends (e.g., the solid Earth response to last deglaciation and gravitational effects of on-going land ice melt). Studies have established that trend patterns in thermal expansion fluctuate both in space and time in response to ENSO (El Nino-Southern Oscillation), NAO (North Atlantic Oscillation) and PDO (Pacific Decadal Oscillation) (Lombard et al., 2005). Thus, it is important to know past regional variability and see how it evolves with time. Climate model projections suggest significant regional variability with respect to the global sea level rise for the end of this century (Meehl et al., 2007). But account for the interannual/decadal variability associated with ENSO and other phenomena by coupled climate models is still imperfect insight into past regional variability over time spans longer than the altimetry record may be helpful to improve the climate models.
Unfortunately, for the last century, information about sea level is sparse and essentially based on tide gauge records along islands and continental coastlines. This data set cannot alone inform on open ocean regional variability. For that reason, a number of previous studies have attempted to reconstruct past decades sea level in two dimensions (2-D), combining sparse but long tide gauge records with global Published by Copernicus Publications on behalf of the European Geosciences Union. gridded (i.e., 2-D) sea level (or sea level proxies) time series of limited temporal coverage (Chambers et al., 2002;Church et al., 2004;Berge-Nguyen et al., 2008). The present study has a similar objective: it expands an earlier work by Bergé-Nguyen et al. (2008) (hereafter denoted as BN08) but makes use of different information for the 2-D fields used for the reconstruction. Previous studies used global sea level grids based on Topex/Poseidon satellite altimetry of limited (10 to 15 years) temporal coverage (e.g., Chambers et al., 2002;Church et al., 2004) or long but spatially inhomogeneous gridded time series of thermal expansion based on in situ hydrographic data (e.g., BN08). In this study, we use global dynamic heights grids from an Ocean General Circulation Model (OGCM), the OPA/NEMO model constrained by data assimilation (Madec et al., 1998). These model outputs, available over a 46-year time span , are combined with tide gauge records (that cover the period . We consider 44-year (1960-2003) time span for the OPA/NEMO outputs to be in line with the tide gauge records length (that ends in 2003). The advantage of using such a long 2-D data set is twofold: (1) the 44-year long coverage minimizes the probably non stationarity of altimetry-based spatial patterns (see BN08 for a discussion), (2) the combination of model-data resulting from an assimilation approach solves the problem of poor geographical and deep ocean coverage of in situ hydrographic data. The resulting sea level reconstruction is presented below for the 1950-2003 time span.

Method
Several studies have developed methods for reconstructing past time series of oceanographic (e.g., sea surface temperature, sea surface height) or atmospheric (e.g., surface wind speed, surface pressure) fields by combining 2-D grids of limited temporal coverage (in general available from remote sensing observations over the last 2 decades or less) with historical (several decade-long), sparse 1-D records (e.g., Smith et al., 1994Smith et al., , 1996Kaplan et al., 1997Kaplan et al., , 1998Kaplan et al., , 2000Church et al., 2004Church et al., , 2006Chambers et al., 2002;Beckers and Rixen, 2003;Alvera-Azcarate et al., 2005;Rayner et al., 2003). The general approach uses Empirical Orthogonal Empirical Functions (EOF) decomposition (e.g., Preisendorfer, 1988) of the 2-D time series to extract the dominant modes of spatial variability of the signal. These EOF spatial modes are then fitted to the 1-D records to provide reconstructed multidecade-long 2-D fields. Different computational variants of the method have been developed to estimate the reconstructed long-term 2-D fields depending on the use of a priori information and data errors (e.g., Kaplan et al., 2000;Rayner et al., 2003;Church et al., 2004) or not (Smith et al., 1996). When applied to long-term past reconstruction, an implicit assumption of the method is the temporal stationarity of the spatial patterns recovered from the EOF modes of the short-term 2-D fields. The spatial covariance obtained from the 2-D fields must indeed be able to describe the spatial covariance over the entire period of the reconstruction (e.g., Smith et al., 1996). If the dominant modes of spatial variability have characteristic time scales longer than the time interval covered by the 2-D fields, EOF spatial modes may incompletely capture the relevant long-term signal. In particular the multi-decadal components of the reconstructed signal may be in error.
We briefly summarize below the methodology. Let us call Fo(x,y,t) and Go(x,y,t) observed global gridded short-term fields and sparse, incomplete long-term sea level data respectively, with x, y and t being Cartesian coordinates and time. The time span T f covered by the Fo(x,y,t) fields is basically shorter than that -called Tg-of the reconstructed fields. Here, Fo(x,y,t) corresponds to gridded sea level from the OPA/NEMO models over 1960-2003, while the Go(x,y,t) -spatially incomplete-fields correspond to tide gauge records over 1950 function is expressed as a sum of combined X n (x, y) spatial modes and e n (t) principal components using an EOF decomposition (with zero global mean trend). The objective of the reconstruction is to compute 2-D Go(x,y,t) fields with global spatial coverage -hereafter denoted as G R (x,y,t) -over the Tg time span (here 1950-2003). The 2-D reconstructed sea level fields are written as where Y n (t) are new principal components computed at each time step t and mode n, through a least-squares fit that minimizes the quantity ε expressed by ε=[Go(x,y,t)− [X n (x, y)Y n (t)]] 2 . For more details, see BN08.

Tide gauge data
The tide gauge data used in this study are extracted from the Permanent Service for Mean Sea Level (PSMSL) database (Woodworth and Player, 2003). We use Revised Local Reference (RLR) tide gauge records (annual averages). Detailed descriptions of these time series are available at www.pol.ac. uk/psmsl. From the whole set of records available, we consider stations that have almost complete temporal coverage over 1950-2003. A very careful selection of sites has been realized. Compared to the 118 sites considered in BN08, here we use only 99 sites, deleting from the 118 set a number of tide gauge records with suspect behaviors (e.g., offsets). Search of information on internet about the deleted sites indicates that in almost all cases, tectonic (e.g., co seismic offset or post seismic relaxation of the crust), volcanic or ground subsidence of anthropogenic origin, could be identified as causes of the spurious bias or trends. The tide gauge records are corrected for glacial isostatic adjustment (GIA) using the ICE-5G, VM-2 model (Peltier, 2004). We also correct the tide gauge time series for the inverted barometer response of sea level to atmospheric loading using surface pressure fields from the National Centers for Environmental Project (NCEP) (Kalnay et al., 1996). One problem with tide gauge records is that measurements are made in local datum that varies from one site to another. By working with the derivatives, this problem can be overcome (e.g., Holgate and Woodworth, 2004;Holgate, 2007). Here we choose a different approach (e.g., Kuo et al., 2008) consisting of subtracting from each sea level record a mean value computed over the 1950-2003 time span (note that the 99 tide gauge records are almost complete; when small gaps,<3 years, are observed, we linearly interpolate missing data). Figure 1 shows the distribution of the tide gauge sites used in this study (superimposed on a map of OPA/NEMO dynamic height spatial trends).

OPA/NEMO Ocean General Circulation Model outputs
The ocean reanalysis used in this study has been produced with a 3-dimensional variational assimilation system (Daget and Weaver, 2008) applied to the OPA/NEMO OGCM (Madec et al., 1998). The model resolution, 2 • on average, with a latitudinal refinement in the tropics, is coarse but appropriate for the multidecadal period of interest here. The model is forced by the standard reanalyse ERA40 heat fluxes (Uppala et al., 2005) and corrected water fluxes (Troccoli and Kalberg, 2004). From September 2002 onwards, when ERA40 terminates, ECMWF (European Centre for Meteorological Forecast) operational surface fluxes are used as forcing. The assimilation system minimizes the discrepancy between model and observations by constraining the model to remain close to the a priori model state. In this iterative procedure, the distance to the observations (respectively to the model state) is taken as the norm defined by the observation error (respectively the model state error). Whilst the observational error is straightforward to characterise, a lot of work has been done to properly define the model state error (Ricci  Weaver et al., 2006). This ensures the propagation of information to the model variables not directly observed (e.g., sea level and velocity) and hence the realism of the analyses. Quality controlled temperature and salinity profiles from the EN3 oceanographic data base (Ingelby and Huddleston, 2007) are assimilated every 10 days from January 1960 to December 2006. Note that these profiles were not corrected for the recently discovered instrumentation problems. The model outputs are on a monthly basis but only annual averages are considered hereafter. As described in Roullet and Madec (2000), the model is formulated with a prognostic free surface, constant volume and salt preserving scheme, assuming that the mean sea level does not vary. Hence, both forcing and data assimilation have been designed to ensure that no drift in the mean sea level occurs. Water flux balance between precipitation, evaporation and runoff is set to zero in the free surface equation, and assimilation increments (i.e., the model corrections applied every 10 days to the model using temperature and salinity observations) are built under the constraint that the sea level increment is zero on global average. In addition, the model is relaxed towards a climatology (Levitus et al., 1994) poleward of 60 • and in the semi-enclosed seas, so that interannual variations in those regions are almost suppressed. On the other hand, no constraint is applied over open ocean areas. In Fig. 1 are displayed sea level trends from the OPA/NEMO model with data assimilation over 1960-2003. For the reasons explained above, the map has zero global mean (uniform) trend. Important regional variability is observed, especially in the western parts of the basins. The strongest spatial patterns are located close to the western boundary currents (e.g., Gulf Stream in the Northwestern Atlantic, Kuroshio in the Northwestern Pacific; Malvinas Current in the Southwestern Atlantic).
Here we focus on the regional variation and do not attempt to reconstruct the global mean trend. However, we must pay attention on a possible contamination of any global mean trend to the reconstructed sea level. We first checked that the model EOFs have zero global mean trend (this is expected as the model does not contain any uniform trend signal). We further checked that EOFs of the reconstruction do not contain either significant non-zero global average.

Reconstructed spatial trend patterns over 1950-2003
We In the following, we consider as nominal case (called case 1), the first 10 modes when reconstructing sea level (with 75% of the total signal variance). Corresponding reconstructed spatial trend map over 1950-2003 is presented in Fig. 2a.
We note that trend amplitudes are everywhere higher than in the model trend map (Fig. 1) but extrema are located in the same regions (e.g., Gulf Stream, Kuroshio; Malvinas Current, etc.). For comparison we also show the spatial trend map with the first 20 modes used for the reconstruction (Fig. 2b). The 20 modes case contains more energy but is also noisier (as discussed in BN08).
In order to assess the benefit of using multi-decadal gridded OPA/NEMO time series for reconstructing past sea level, we performed tests with shorter model time series. Four cases are considered: 31 years of OPA/NEMO grids  -case 2-, 21 years of OPA/NEMO grids  Clim. Past, 5, 217-227, 2009 www.clim-past.net/5/217/2009/  -case 3-, and 11 years of OPA/NEMO grids (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) -case 4. Case 4 can be compared with studies that use decade-long Topex/Poseidon altimetry grids for sea level reconstruction (e.g., Chambers et al., 2002;Church et al., 2004). Corresponding reconstructed spatial trend maps are presented in Fig. 3a, b, c (using for each case, the number of modes that correspond to 75% of the total variance). Significant regional differences are noticed between the first two cases (44 years and 31 years of OPA/NEMO grids, respectively) and cases 3 and 4 (21 years and 11 years of OPA/NEMO grids, respectively) , in particular in the North Atlantic, Indian and Austral oceans, and Northeast Pacific.   from ∼20 years to ∼30 years. As discussed below, the lowfrequency sea level signal may not be well captured in cases that use short spatial grids for the reconstruction (cases 3 to 5). In the latter cases, patterns representing interannual variability dominate the reconstructed spatial trends.
To compare with case 4, we also show in Fig.3d reconstructed sea level trends (over 1950-2003) using 11 years (1993-2003)   We performed EOF decompositions of the reconstructed sea level grids over 1950-2003 for cases 1 to 4. Figure 4 a, b, c, d shows the two leading EOF modes each case. Mode 1 temporal curve (principal component) of nominal case is dominated by a positive slope. Associated mode 1 spatial map closely resembles case 1 reconstructed trend map (Fig. 2), with strong signal in the Austral Ocean (especially southeast of Africa). This suggests that the reconstructed trend map for the nominal case mostly reflects a long-term (multi-decadal), regionally variable signal. Mode 2 of case 1 (Fig. 4a right panel) is dominated by the interannual variability, and displays clear signature of ENSO in the tropical Pacific: the Southern Oscillation Index (SOI) -a proxy of ENSO-is significantly correlated (61%) to the temporal curve on which is it superimposed. It is worth mentioning that the spatial map of mode 2 (case 1) is very similar to the satellite altimetry-based sea level trend map over 1993-2003 (see below). Hence, whilst mode 1 reflects long-term, multidecadal signal, mode 2 mostly reflects ENSO-type interannual variability.
The two leading modes of case 2 (Fig. 4b) closely resemble those of case 1. The first EOF mode of cases 3 and 4 reflects interannual variability (as mode 2 of cases 1 and 2). This is illustrated in Fig. 5 which shows energy spectra of the interannual signal for cases 1, 4 and 5 (i.e., temporal curves spectra of corresponding EOF modes). SOI spectrum is superimposed. The agreement between the four curves is striking. Peaks in the 3-5 yr and 10-15 yr wavebands dominate. They mainly reflect ENSO frequency and associated decadal modulation. On the other hand, looking at the long-term signal (e.g., comparing modes 1 of case 1 and 2 with mode 2 of cases 3 and 4), we note significant difference in the spatial maps, suggesting that multidecadal fluctuations are only partly recovered when using short-term gridded time series for the reconstruction (e.g., cases 3 and 4).

Sea level reconstruction over the altimetry (1993-2003) period using the OPA/NEMO EOFs
A way to check the validity of the reconstruction is to look at the reconstructed sea level trends over the altimetry period (here 1993-2003) for which we trust the spatial trend patterns. Figure 6a shows the reconstructed spatial trend map over 1993-2003 based on 44-years of OPA/NEMO grids. Comparing with Fig. 6b which shows observed satellite altimetry spatial trend map over 1993-2003 (uniform trend removed) indicate very good agreement, although the reconstructed map is smoother (as expected considering the low resolution of the OPA/NEMO model).

Cross-validation of reconstructed series and tide gauge records
Another way to check the robustness of the reconstruction consists of reconstructing sea level, leaving out successively each one of the 99 tide gauge records (thus each of these 99 reconstructions now uses a set of 98 tide gauges, with   different distribution from one case to another). For each deleted tide gauge, we compare the reconstructed sea level time series at the tide gauge site with the PSMSL data (reconstructed sea level is averaged within 2 • around the tide gauge site). For this test we consider two cases: case 1 (reconstruction with 44 years of gridded OPA/NEMO grids) and case 5 (reconstruction with 11 years of Topex/Poseidon gridded data). Figure 7a shows a subset of 18 comparisons (corresponding site locations are enhanced by stars in Fig. 1). Figure 7b is similar to Fig. 7a except for the trend which has been removed. We note that in general interannual to decadal variability is well reproduced (Fig. 7b). The average correlation at the 99 sites amounts to ∼60%. The trends also agree well in most cases, although not everywhere. Sites where trends disagree concern mostly northeast Atlantic areas (e.g., North Shields, Lowestoft, Santander, La Coruna, and Vigo). We suspect that this is due to local underestimated variability in the OPA/NEMO reanalysis, especially in the first 20 years of the period. Similar test performed for case 5 (not shown) shows very similar results for the detrended time series.   To see the above results in another way, we have computed the root mean squared (rms) differences between reconstructed and observed (detrended) sea level time series as well as between trends. Corresponding histograms for the cases 1 and 5 are shown in Fig. 8. Rms (detrended) sea level differences are included in the 15-60 mm but histogram for case 1 is more spread than for case 5. Figure 9 shows plots of reconstructed sea level trends at the 99 tide gauges as a function of observed trends for the two cases (cases 1 and 5).

Figure 7b
Here we see that case 1 gives better results than case 5, with higher correlation. These comparisons indicate that the reconstruction performs well at the interannual time scale (as previously found by Chambers et al., 2002 andChurch et al., 2004), with better results for case 5 than case 1. On the other hand, case 1 gives much better results for the trends (hence at multidecadal time scale) than case 5.

Conclusions
We have performed a new 2-D sea level reconstruction over 1950-2003. The main change compared to previous published results (e.g., Chambers et al., 2002;Church et al., 2004;Berge-Nguyen et al., 2008) is the use of 44-year long gridded sea level data set from the OPA/NEMO OGCM with assimilation. This allows us to compute spatial EOFs over a >40 yr time span, long enough to capture the multidecadal variability in regional sea level (in addition to the interannual variability). Hopefully, this may prevent from problems due to the use of short-term altimetry grids. Another advantage is the good spatial sampling of the OPA/NEMO outputs compared to in situ-based thermal expansion data (as in BN08). The main conclusion of this study is that using spatial EOFs computed over a short time span (e.g., ∼10-20 years) leads to 2-D reconstructed trend patterns significantly different than when 40 years of EOFs are used. The latter case, not only captures the interannual variability (related to ENSO and possibly NAO and PDO) but also the multidecadal variability. Even if they have local deficiencies, such 2-D reconstructed sea level time series based on low resolution ocean reanalysis are able to bring physically consistent large-scale, low-frequency patterns associated with major climatic modes of variability. As a final remark, we think that the use of OGCM outputs is a step towards better reconstruction of long-term sea level time series (at least waiting for global multidecadal altimetry records). Future improvements is expected by using new generation of eddy-permitting OGCM outputs with higher spatial resolution (e.g., 0.25 • ×0.25 • ) in which some local misbehaviors can be corrected. This would provide finer description of the spatial trend patterns. Progress in this direction is already underway.