Radar altimeters permanently transmit signals to the Earth, and receive the reflected echo from the sea surface. The satellite orbit has to be accurately tracked and its position is determined relative to a reference surface (an ellipsoid). The sea surface height (SSH) is calculated by subtracting the measured distance between the satellite and sea surface from the precise orbit of the satellite. The sea level anomalies (SLA), defined as variations of the sea surface height (SSH) with respect to an a priori mean sea surface (MSS), are generally used as precious and main indicators for the development of scientific applications that aim to study the ocean variability (mesoscale circulation, seasonal variations, El Niño…).

In our research, we investigate a temporal series of averaged sea level anomalies (SLA) over all oceans and seas for the entire altimetry period, obtained by combining the time series from all three TOPEX, Jason-1 and Jason-2 missions (Nerem et al., 2010). These datasets were provided by the CU Sea Level Research Group of the University of Colorado¹. This series covers the time period from the beginning of 1993 to the end of 2012 (20 years), with a sampling rate of 10 days.

The used mean sea level anomalies series is represented on Fig. 1. The main input data for computing this series are the Geophysical Data Records produced by NASA and CNES (TOPEX, Jason-1, Jason-2), which are therefore of the highest quality, notably in orbit determination. The three missions are linked together during the “verification” phases of the Jason-1 and Jason-2 missions in order to calculate very precisely the bias in GMSL between these missions. TOPEX and Jason-1 were connected by applying a bias of 88.48 mm to Jason-1's GMSL. Similarly, Jason-2's GMSL was connected to Jason-1's GMSL by applying a bias of 58.52 mm to Jason-2's GMSL and also by adding the bias between Jason-1 and TOPEX/Poseidon. All of the standard corrections to the altimeter range were applied to the sea surface height (SSH), including removal of ocean tides and an inverted barometer correction. The used mean sea surface (MSS) is the CLS01 model (Hernandez and Schaeffer, 2001), obtained from 7 years of TOPEX/Poseidon data, 5 years of ERS-1/2 data, 2 years of Geosat data added to ERS-1 geodetic data, which were processed and homogenized in order to produce a (1/30° × 1/30°) gridded map. Table 1 describes the parameters used by the University of Colorado for computing global mean sea level reference series.

Monitoring of ENSO conditions primarily focuses on the sea surface temperature anomaly (SSTA) index (variations from an average or other statistical reference value) in the equatorial Pacific. SSTA equal to or greater than 0.5 °C in the Niño3 region (5°N to 5°S, 150°W to 90°W), or exceeding 0.4 °C in the Niño3.4 region (5°N to 5°S, 170°W to 120°W), could indicate ENSO warm phase (El Niño) conditions (Trenberth, 1997), while anomalies less than or equal to −0.5 °C (or −0.4 °C for the Niño3.4 region) are associated with cool phase (La Niña) conditions.

In this paper, the used ENSO index was the average of SSTA computed over the Niño3 region (5°S to 5°N, 90°W to 150°W), with respect to a 1971–2000 month climatology (Xue et al., 2003). These data sets were obtained from the most recent version (v3b) of the NCDC Extended Reconstruction Sea Surface Temperatures (ERSST) analysis², monthly from January 1880 to June 2010. ERSST.v3b was generated using in-situ Sea Surface Temperature (SST) data and improved statistical methods that allow stable reconstruction using sparse data. The monthly analysis extends from January 1854 to the present, but because of sparse data in the early years, the analyzed signal was damped before 1880. After 1880, the strength of the signal was more consistent over time. ERSST is suitable for long-term global and basin wide studies; local and short-term variations have been smoothed in ERSST. The ERSST.v3b analysis is exactly as described in the ERSST.v3 paper (Smith et al., 2008), with one exception: satellite SST data are not used in ERSST.v3b.

The UCM are models in which time series are decomposed as the sum or a product of a number of other simple time series with economic or physical meaning. One widely accepted univariate version of UCM is (Young et al., 1999):

In practice, not all the components in eq. (1) are necessary: indeed, the simultaneous presence of all these components can induce identification problems and so use of the complete model (1) is not advisable in practice unless adequate precautions are taken. In this paper, the simplified used model is one such decomposition that contains only the trend plus an AR component, seasonal and white noise components.

The trend component is formulated in a stochastic state space context and then solved by using an integrated random walk (IRW) plus noise modelling approach. The full state space description is given by the state equation (Pedregal and Trapero, 2007):

x_1t is the first state variable known as the ‘level’, x_2t is the second state variable and generally known as the ‘slope’, η_t is a zero mean, serially uncorrelated white noise process with a certain variance. Young et al. (1999, 2007) also defined two parameters for computing convenience: noise variance ratio (NVR) Q_r and pˆt. Q_r is a normalized matrix of the diagonal noise variance matrix Q, and pˆt is the error covariance matrix associated with the estimation state. Young et al. (1999) propose one final simplification using the estimate of the trend from an autoregressive (AR) model. The AR order is identified by the Akaike's Information Criterion (AIC), whereas, the seasonal term is modelled as a sum of the fundamental signal and of its associated harmonics:

The cross wavelet transform and wavelet coherence analysis used in this paper follows the method of Grinsted et al. (2004). Although the basic components of this analysis are reviewed here, readers are referred to Grinsted et al. (2004) for a more detailed explanation.

The continuous wavelet transform (CWT) allows analyzing the temporal evolution of the frequency content of a given signal or timing series. The CWT of a time series (x_n, n = 1, …, N), with uniform time steps δt, is defined as the convolution of x_n with the scaled and normalized wavelet:

The wavelet power is defined as Wnxs2. The complex argument of Wnxs can be interpreted as the local phase (Torrence and Compo, 1998).

From two CWTs we construct the cross wavelet transform (XWT), which will expose their common power and relative phase in time–frequency space. The XWT of two time series x_n and y_n is defined as:

Following Torrence and Webster (1999), the wavelet coherence (WTC) of two time series, which can find significant coherence even though the common power is low, is defined as:

All of the results in this section are obtained by means of the Captain Toolbox v7.5³ (Taylor et al., 2007). The Captain Toolbox is a set of MATLAB functions for non-stationary time series analysis and forecasting. It allows the identification of unobserved components models, time variable parameter models, state-dependent parameter models and multiple-input discrete and continuous time-transfer function models. The toolbox is useful for system identification, signal extraction, interpolation, backcasting, forecasting and databased mechanistic analysis of a wide range of linear and nonlinear stochastic systems.

The global mean sea level series was decomposed into trend and seasonal components using a model based on a trend plus an AR component. Apart from the decomposition of time series into trend and seasonal components, this methodology also provides a time series model of the data, which is used for the interpolation of missing data. Note that this model can also be used for backcasting and forecasting.

In the first place, we estimated a smooth trend using an IRW with an initial noise variance ratio (NVR) hyper-parameter of 10⁻⁴. The green line and the red line on Fig. 2 correspond to the smooth trend and the residuals series, respectively.

Then, the order of the autoregressive (AR) model for the perturbations is identified using the AIC. We searched for models ranging from 1st to the 10th order. The best-fit AR polynomial is obtained for an order of 9. The AR model for the perturbations is estimated, together with the associated optimal NVR for the IRW trend.

Table 2 represents the estimation results of the perturbational component about the trend, were P is the AR order for perturbations, ARp is the AR polynomials for the perturbations and ARpse is the standard error of AR parameters for perturbations.

The final trend NVR is estimated at 2.7554 × 10⁻⁴, with a trend integration order of 2. Finally, both components (trend and seasonal component) are estimated simultaneously.

Figs. 3, 4 and 5 represent the estimated trend, the total seasonal component and the fitted data, respectively. A zoom on the final 2 years, with the standard errors, is shown on Fig. 6.

The analysis of the estimated trend indicates that the global mean sea level is subject to significant rise, from −6 mm to 62 mm during the period 1993–2012. Applying a least squares linear regression analysis to the extracted trend gives a rate of 3.22 mm/year between 1993 and 2012, with a formal error of 0.02 mm/year (see Fig. 7).

The estimated rate of GMSL rise agrees with the CU GMSL rate⁴ over the same period of 3.2 ± 0.4 mm/year, where the realistic error assigned in this quantity (± 0.4 mm/year) is based on comparisons of altimeter heights with tide gauge-based sea level measurements (Mitchum, 2000). For further information on errors assessment affecting the calculation of the GMSL rate and comparison of altimetry-based sea level with in-situ measurements, the reader is referred to Ablain et al. (2009).

One can see on Fig. 7 evidence that there are a number of fluctuations over short periods in the estimated trend. This variability is at least partly related to El Niño/La Niña episodes (sea level rises during El Niño and falls during La Niña) and associated changes in the hydrological cycle.

To compare the global mean sea level to the sea surface temperature anomalies index over the Niño3 region, we removed from original global mean sea level time series (series of mean sea level anomalies) the extracted trend and the total seasonal components (both estimated in Section 4.1). Then, since ERSST.v3b SSTA data are monthly distributed, we created from the detrended GMSL datasets a lower resolution series with a time scale of one month to be used for relationship analysis.

Fig. 8 shows the used averaged SSTA time series (1993–2012) over the Niño3 region and the detrended and seasonally adjusted GMSL series used for correlation analysis. The plotted series show a strong correlation between the global mean sea level and the SSTA Niño3 index.

One can see from Fig. 8 evidence that there was a particularly strong warm event of large amplitude during the 1997–1998 period. The El Niño of 1997–98 was extremely severe. This event was a so-called “Type-1” El Niño (Fu et al., 1986), having the strongest averaged SSTA > 3.0 °C from October 1997 to January 1998 (see Fig. 8). Although it started in April to May 1997, its effects extended into the early summer of 1998.

In early 2006, sea surface temperature positive anomalies are measured in the central equatorial Pacific. These high-temperature anomalies characterize a new El Niño event. In early summer 2009, sea surface height positive anomalies are measured from August 2009, with a maximum during November 2009 (SSTA > 1.1 °C), characterizing a moderate El Niño event (see Fig. 8).

Severe La Niña was also formed in 1998 and 2012. From June 2007, data indicated a moderate La Niña event, which strengthened until early 2008. A new La Niña episode was developed in mid-2010, and lasted until early 2011. It intensified again in the mid-2011 (see Fig. 7). Between the beginning of 2010 and the middle of 2011, sea level fell until −3 mm. This occurred concurrently with the strong phase, La Niña 2011.

For examining strictly speaking relationships between the two time series SSTA and GMSL (detrended and total seasonal component removed) series in time–frequency space, first the CWT of both datasets was calculated. The XWT and WTC tests are then performed.

The CWT of the SSTA index and GMSL series are shown on Fig. 9. There are clearly common features in the wavelet power of the two time series in the 2–4-year band in the period from 1996 to 2000. However, the similarity between the two patterns for other periods is quite low and it is therefore hard to say whether it is significant. The cross wavelet transform helps in this regard; it will expose the common power and relative phase in the time-frequency space.

The XWT of SSTA index and GMSL and their relative phase relationship (with in-phase: arrows pointing right, anti-phase: arrows pointing left) are shown on Fig. 10. Here we note that there is a significant common power in the 1–4-year band from 1996 to 2000. Also, there is an appreciable power in the 1–2-year band from 2005–2010. All these areas show an in-phase relationship, suggesting that GMSL and SSTA vary synchronously.

The squared WTC of SSTA and GMSL is shown on Fig. 11. Compared with the XWT, a larger section stands out as being significant and all these areas show an in-phase relationship between SSTA and GMSL. Oscillations in SSTA are manifested in the GMSL on wavelengths varying from 1–4 years in the period from 1996 to 2000 (El Niño 1997–98), 1–2 years in the period from 2005 to 2008 (2006–2007 El Niño event and 2007–2008 La Niña episode). Also, there is an appreciable section in the 3–4-year band from 2007–2012 (2009–2010 El Niño, 2010–2011 La Niña events), though this area is outside of the cone of influence in which edge effects are considerable.