Cyclostationary empirical orthogonal function sea‐level reconstruction

Since 1993, satellite altimetry has provided accurate measurements of sea surface height with near‐global coverage. These measurements led to the first definitive estimates of global mean sea‐level (GMSL) rise and have improved understanding of how sea levels are changing regionally at decadal time scales. These relatively short records, however, provide no information about the state of the ocean prior to 1993, and with the modern altimetry record spanning only 20 years, the lower frequency signals that are known to be present in the ocean are difficult or impossible to resolve. Tide gauges, on the other hand, have measured sea level over the last 200 years, with some records extending back to 1807. While providing longer records, the spatial resolution of tide gauge sampling is poor, making studies of the large‐scale patterns of ocean variability and estimates of GMSL difficult. Combining the satellite altimetry with the tide gauges using a technique known as sea‐level reconstruction results in a dataset with the record length of the tide gauges and the near‐global coverage of satellite altimetry. Cyclostationary empirical orthogonal functions (CSEOFs), derived from satellite altimetry, are combined with historical sea‐level measurements from tide gauges to create the Reconstructed Sea Level dataset spanning from 1950 to 2009. Previous sea‐level reconstructions have utilized empirical orthogonal functions (EOFs) as basis functions, but by using CSEOFs and by addressing other aspects of the reconstruction procedure, an alternative sea‐level reconstruction can be computed. The resulting reconstructed sea‐level dataset has weekly temporal resolution and half‐degree spatial resolution.


Introduction
Sea level is a measurement of considerable interest and importance for the study of climate because it reflects both the mass and heat storage changes in the global ocean. Since 1993, satellite altimetry has provided accurate measurements of sea level. The near-global coverage of the satellite altimeters has improved the understanding of how sea level is changing on regional and global scales. In addition, the first definitive estimates of global mean sea-level (GMSL) rise have been made using satellite altimetry data (e.g. Cazenave & Nerem, 2004;Leuliette et al., 2004;Beckley et al., 2007).
The modern satellite altimetry record, however, spans only 20 years, making it difficult to separate longer term secular trends and accelerations from natural climate variability (Church & White, 2006Woodworth et al., 2009;Hamlington et al., 2011a). The large positive trend in GMSL in the past 20 years, when compared to the estimated trend over the past 60 or even 110 years, suggests the altimetric rate of sea-level rise may be exceptional relative to time periods in the past (Church & White, 2011). Furthermore, sea-level variations on decadaland longertimescales are known to contribute to the sea-level trends on short timescales. Understanding and evaluating these trends in sea level is challenging using only the short satellite altimeter record. To make accurate comparisons between sea-level variations over different time periods, a long and consistent data record is necessary. In general, studying past sea-level variations is vital to improving the understanding of current and future sea-level change, as predicting future change is predicated on the ability to accurately quantify and represent the state of sea level in both the past and present.
Sea-level reconstructions provide a possible solution to this problem (Chambers et al., 2002;Church et al., 2004;Hamlington et al., 2011bHamlington et al., , 2012Meyssignac et al., 2012). Historical measurements of sea level from tide gauges extend back as far as the beginning of 19th century. While providing long records, the measurements of tide gauges are generally sparse, particularly before 1950. By combining the dense spatial coverage of satellite altimetry with the long record length of the tide gauges, however, it is possible to examine longer time scale climate signals and assess their contribution to sea-level trends both regionally and globally. Furthermore, it is possible to determine whether the current rate and spatial pattern of sealevel change are exceptional or instead are simply a recurrence of multi-decadal climate oscillations.
Here, we present a reconstructed sea-level dataset created by fitting cyclostationary empirical orthogonal functions (CSEOFs), derived from satellite altimetry, to historical sea-level measurements from tide gauges, to reconstruct sea-level fields from 1950 to 2010. Previous sea-level reconstructions have primarily utilized empirical orthogonal functions (EOFs) as basis functions, but by using CSEOFs and by addressing other aspects of the reconstruction procedure, an alternative andin some waysimproved sea-level reconstruction can be computed.

Cyclostationary EOFs
Previous reconstructionsboth sea surface temperature (SST) and sea levelhave generally relied on EOFs to form the basis for the reconstruction. When compared to CSEOFs, however, EOFs have characteristics that make them suboptimal for use as basis functions for sea-level reconstruction. EOFs enforce a stationarity on the spatial variability in the resulting reconstruction. A single spatial map defines the basis function, and the reconstruction procedure simply computes the amplitude modulation of this map through time. Given the evidence that many signals in geophysical data are cyclostationary, CSEOFs provide significant advantages over EOFs when dealing with signals such as modulated annual cycle (MAC) and El Nino-Southern Oscillation (ENSO) signals. Although EOFs are able to represent cyclostationary features through a superposition of multiple modes, CSEOFs are able to explain cyclostationary signals in a single mode, increasing the opportunity for interpretability.
The decomposition of data in terms of a set of basis functions is often very useful in understanding the complicated response of a physical system. By decomposing into less complicated patterns, it may be easier to understand and shed light into the nature of the variability in a dataset. While theoretical basis functions have been studied extensively, exact theoretical basis functions are very difficult to find and in general, computational basis functions are sought instead. Perhaps the simplest and most common computational basis functions are EOFs. Consider a simple system defined by the equation: where LV(r) is a physical process (termed to be the loading vector) modulated by a stochastic time series PC(t), which is called the principal component time series (PCTS). Each loading vector and PCTS pair represents a single EOF mode. As mentioned above, however, physical processes and the corresponding statistics are time-dependent. Representing the data with stationary EOFs can lead to erroneous and difficult interpretation of the data (Dommenget & Latif, 2002). Kim et al. (1996), Kim and North (1997), Kim and Wu (1999), Kim and Chung (2001) introduced the concept of CSEOF analysis to more compactly capture the time-varying spatial patterns and longer timescale fluctuations present in geophysical signals. The significant difference between CSEOF and EOF analysis is the LVs' time dependence, which allows the spatial pattern of each CSEOF mode to vary in time, with the temporal evolution of the spatial pattern of the CSEOF LVs constrained to be periodic with a selected "nested period". In other words, the system is defined as: where the loading vectors are now time-dependent, and are periodic with the nested period, d. As a result, each CSEOF mode is composed of 12 LVs and one PCTS when, for example, using monthly data and a 1year nested period. CSEOF analysis minimizes mode mixing, which is a common problem with EOF decomposition. When mode mixing occurs, the annual cycle frequently spreads across several modes, which is one reason the signal is usually removed from the data by some other means. Recent studies, however, have demonstrated the efficacy of using CSEOFs to extract robust modes representing the MAC and ENSO variability (Hamlington et al., 2011a,b). This leads to robust estimates of the MAC or ENSO variability from satellite altimetry data without affecting signals associated with other ocean variability. For further details on CSEOFs and the procedure for computing CSEOFs, the reader should refer to Kim and North (1997) in which a detailed description of the computation of CSEOFs is provided. By fitting CSEOFs in place of EOFs to tide gauges to reconstruct sea level for the 1950-2010 time period (Hamlington et al., 2011b), we have found that it is possible to create an alternate and, we believe, improved reconstruction to those based on EOFs. Our motivation for using CSEOFs, as described in Hamlington et al. (2011b), in place of EOFs is fourfold: (1) EOFs are not an optimal basis for non-stationary signals with nested oscillations that are undergoing lowfrequency oscillation; (2) CSEOFs account for both the high-and low-frequency components of the annual cycle and do not require the removal of the annual signal from the satellite altimetry nor the tide gauge records before reconstruction; (3) specific signals, such as those relating to the MAC and ENSO can be reconstructed individually using CSEOFs with little mixing of variability between modes; and (4) the reconstruction procedure using CSEOFs spans data gaps smaller than the nested period and fits a larger window of data, allowing for the use of fewer historical data to obtain meaningful results.

Satellite altimetry data
CSEOF sea-level basis functions for our reconstruction were estimated from the Archiving, Validation and Interpretation of Satellite Oceanographic data (AVISO) quarter-degree resolution multiple altimeter product based on satellite altimeter measurements spanning 1993-2011 collected by the Topex/Poseidon, ERS-1&2, Geosat Follow-On, Envisat, Jason-1 and OSTM satellites. This updated and reprocessed gridded data product, which was released in December 2011, was created using the delayed time Ssalto/DUACS multimission altimeter data processing system with improved homogenous corrections and inter-calibration applied to the entire data record. Global crossover minimizations and local inverse methods are used to derive inter-calibrated highly accurate along track data that are referenced to a consistent mean. The along track data were then merged through a global spacetime objective mapping technique that takes into account correlated noise (Le Traon et al., 1998).
We applied very little additional processing other than removing the mean and a linear least square fit from the time series at each spatial grid point. A CSEOF decomposition of the satellite altimetry data is not able to extract the change in mean sea level into a single mode due to the short record length associated with the satellite altimeters. It is therefore necessary to remove mean sea level from the satellite altimetry data before computing the basis functions to avoid putting low-frequency power into each CSEOF mode. Although our technique does result in the removal of the spatial pattern of sea-level trends, it is unlikely that the regional distribution of sea-level trends over the past two decades is the same as that since 1950 and it is, therefore, unwise to force this stationary pattern upon the reconstruction. Removing the trend from each grid point does not significantly affect the ability of the reconstruction to capture decadal time-scale signals. Using the study of Tai (1989), the percentage of signal reduction at various periods caused by removing the linear trend from a 17-year record can be computed. Signals with periods of~10 years undergo an RMS signal reduction of less than 5% and even signals at 20-year periods are reduced by only 30%. Linear detrending of the altimeter data before computing the CSEOFs should therefore have little effect on decadalscale variations and our ability to capture them in the reconstruction. The seasonal signal was not removed from the time series prior to computing the CSEOFs. We do not use the CSEOF mode representing the seasonal signal in the reconstruction; we find that there is little difference in the reconstruction regardless of whether or not the annual cycle is removed from both the historical data and satellite altimetry data. Before CSEOF decomposition, the data were weighted using the square root of the cosine of latitude to yield areaweighted variance decomposition.

Tide gauge data
The central sea-level dataset we use for the period 1950-2010 is monthly mean sea-level records gathered from the data archive at the Permanent Service for Mean Sea Level (PSMSL). We use only the Revised Local Reference (RLR) data, which are measured sea level at each site relative to a constant local datum over the complete record. At present, we have not included the metric data offered by PSMSL as they can have substantial unknown datum shifts and their use in time series analysis is generally not recommended.
The editing procedures of Church et al. (2004) were followed closely, and a similar set of tide gauges is used in our reconstruction. We upsampled the monthly tide gauge data to weekly sampling using linear interpolation to match the 1-week temporal resolution of the altimeter-derived CSEOFs. We found the nearest grid point for each tide gauge as the basis functions obtained from the satellite altimetry are on a half-degree resolution spatial grid. Finally, the available tide gauge records were averaged to produce a single time series if there were multiple tide gauges associated with a single spatial grid point. After the editing procedure, 407 tide gauges remained for use in our reconstruction.
The PSMSL sea-level data are relative sea level; therefore, the records must be corrected for the ongoing glacial isostatic adjustment (GIA). We use the ICE-5G VM2 model (Peltier, 2004) to perform the GIA correction. As an inverted barometer correction is applied to the satellite altimetry data, the sea-level measurements from the tide gauges were corrected using an inverted barometer response of sea level to atmospheric loading based on the pressure fields from the NCEP/NCAR reanalysis. It should be noted, however, that neither of these corrections were found to have a noticeable effect on the resulting sea-level reconstruction.

Methods
Using basis functions computed from a short, spatially dense dataset to interpolate a long time series of spatially sparse observations were first implemented in SST studies. Smith et al. (1996) computed EOFs from 12 years of satellite-derived SST data and used them as basis functions to estimate global SST temperature fields from 1950 to 1992. Kaplan et al. (2000), improved on this procedure by adding weighting dependent on known errors in the data to the reconstruction procedure. Sea-level reconstructions soon followed using the techniques developed for SSTs.

Reconstruction procedure
To perform the CSEOF reconstruction, a CSEOF decomposition of the AVISO merged satellite altimetry data from 1993 to 2011 was performed using a nested period of 1 year. As mentioned in section 1.2, before decomposition, the time series at each grid point is detrended by removing a least squares fit linear trend, and the resulting value is area weighted by multiplying by the square root of the cosine of the latitude. The first CSEOF mode represents the MAC and is not included in the reconstruction. The second mode represents the variability associated with ENSO (Figure 1). The top panels of Figure 1 show the CSEOF loading vector, with one spatial pattern for each month (the decomposition is conducted weekly, but for ease of display, the weekly maps are averaged to 12 monthly maps). The middle panel shows the PCTS for the second CSEOF mode, describing the amplitude modulation (strength) of ENSO since 1993. Finally, the bottom panel for Figure 1 shows the contribution of the second CSEOF modeand thus ENSOto GMSL.
Once obtained, the CSEOF modes are used as basis functions and fit to the tide gauge data using weighted least squares. Each CSEOF basis function spans the 1-year nested period, so it is necessary to fit each basis function to a full year of tide gauge data. This leads to a loss of a half a year at each end of the reconstruction, resulting from the 1-year windowing. The goal of the reconstruction procedure is to reproduce the PCTS of each CSEOF mode back through time. The mathematical details of this computation are found in the works by Kaplan et al. (2000) and Church et al. (2004) and thus are not repeated here. Once the PCTS have been reconstructed, they are recombined with the CSEOF basis functions to create the CSEOF sea-level reconstructed dataset. Figure 2 shows the comparison of the PCTS from both the reconstruction and the actual CSEOF decomposition of the satellite altimetry for the second CSEOF mode corresponding to ENSO (Figure 1). From 1993 onwards, the reconstructed PCTS and altimetry PCTS agree very well, highlighting the ability of the recon-struction to represent the satellite altimetry data. A comparison between the Multivariate ENSO Index (MEI) (Wolter, 2010) and the reconstructed PCTS is shown from 1950 to 2010. The agreement over the full record demonstrates the ability of the sea-level reconstruction to capture the ENSO variability even with the limited tide gauge sampling in 1950.

Global mean sea level
Estimating GMSL using reconstruction techniques is not a trivial task. Christiansen et al. (2010) discussed the difficulties of estimating GMSL using EOF reconstruction techniques, albeit using model data instead of satellite altimetry and tide gauge data. While CSEOF basis functions describe the cyclostationary variability in sea level, as a result of the short length of the satellite altimetry record, no single CSEOF mode captures the secular trend. EOF analysis of the satellite altimetry record has similar difficulties in extracting the secular trend. Even if the CSEOF analysis was capable of extracting the secular trend from the AVISO data, we cannot assume that the resulting spatial pattern is stationary over the entire 60-year time period. Church et al., 2004; approximate the trend in their reconstruction by introducing a constant basis function.
Rather than introducing an additional basis function in an attempt to account for the secular trend, we separate the computation of the secular trend from the actual reconstruction procedure. By secular trend, we refer not only to the linear portion of the trend but also to the nonperiodic variations that are unexplained by the reconstruction computed using the CSEOF basis functions. The altimetry-derived CSEOF basis functions are insensitive to the secular trend in the tide gauges largely because of the annual nested periodicity imposed on the spatial and temporal variability in the relatively short altimetric record. Separating the computation of the trend from the reconstruction computation is also necessary because of the complexity of implementing a CSEOF reconstruction with differenced tide gauge data. The time dependence of the LVs would introduce a term in the least squares procedure necessitating the computation of the time derivative of the LVs. For details of how GMSL is included in the reconstruction, the reader is referred to Hamlington et al. (2011b). In short, we have used the CSEOFs to estimate the MAC and ENSO (and other climate variability) and then removed this variability from the weighted average of the tide gauge data. The resulting time series is then added back into the reconstruction, and the GMSL time series (including signals such as ENSO) is computed by averaging the CSEOF reconstruction globally at each point in time.

Dataset location and format
The CSEOF reconstructed sea-level dataset contains sea-level data derived from satellite altimetry and tide gauges using CSEOFs. The data are hosted and main-  tained by NASA Jet Propulsion Laboratory Physical Oceanography Distributed Active Archive Center (JPL PO DAAC).The data span from January 1950 to June 2009 and is periodically updated to extend the time series. The data are in NetCDF format, and each file contains~10 years worth of data. As a result of the reconstruction procedure, the dataset has the temporal and spatial resolution of the satellite altimetry and the record length of the tide gauges. The resulting reconstructed sea-level dataset has weekly temporal resolution and half-degree spatial resolution. The data are provided in units of centimetres.

Dataset use and reuse
With the modern satellite altimetry record spanning only 20 years, separating longer term secular trends and accelerations from natural climate variability is a challenge. Sea-level variations on a range of timescales are known to contribute to sea-level trends on decadal and longer timescales (e.g. Feng et al., 2004;Woodworth et al., 2009;Sturges & Douglas, 2011;Chambers et al., 2012). The "red" nature of the sealevel spectrum makes it difficult to separate trends without a longer record of sea level than that provided by the satellite altimetry. The sea-level reconstruction discussed here provides the opportunity to extract and separate lower frequency climate variability and to study its effect on sea level.
For instance, a recent use of this dataset showed how the Pacific Decadal Oscillation (Mantua et al., 1997;Mantua & Hare, 2002) contributes significantly to 20-year trends in GMSL (Hamlington et al., 2013). With the proven ability to capture signals like ENSO and the PDO over a 60-year sea-level record, a wide range of studies are possible with this dataset. Specifically, one could study how these signals contribute to sea-level variability both regionally and globally. The study of shorter timescale signals is also possible with this dataset, although such research has not been conducted to date. In short, most studies that are possible with a gridded dataset from satellite altimetry can be conducted with the CSEOF sea-level reconstructed dataset with the added benefit of a longer record length. Consideration, however, must be given to (1) limitations imposed by the tide gauge sampling back through time, (2) limitations imposed by the decomposition of satellite altimetry into cyclostationary basis functions, and (3) to the limitations of the techniques required to account for GMSL in the reconstructed dataset. The CSEOF reconstruction is not a direct substitute for the satellite altimetry, but does provide the opportunity to compare sea level in the last two decades to sea level over the past 60 years.