The global mean sea level (MSL) has been rising since the beginning of the twentieth century (IPCC 2013) and more rapidly during the last 20 years (Chen et al. 2017). Quantifying the rate and patterns of MSL change is of great importance (Jevrejeva et al. 2009, Marcos and Amores 2014, Luu et al. 2018) but is still challenging (Hay et al. 2015), mostly due to uncertainty in the measurement (Cazenave et al. 2014, Bos et al. 2014, Dieng et al. 2017) and interpretation (Visser et al. 2015, Dangendorf et al. 2017, Royston et al. 2018) of sea level data. Dating back to the eighteenth century, tide gauge records have long been a source of the spatio-temporal rates of MSL change. However, their limitations are sparseness, uneven geographic distribution, high sensitivity to local geodynamic and hydrological influences and, notably, unavailability in the open oceans. The launch of oceanographic satellites since 1993 provided an unprecedented opportunity to reveal the spatial patterns of sea level change with global coverage, capturing the non-uniformity in trend and variability (Ablain et al. 2015, 2017, Nerem et al. 2018).

In some regions of the western tropical Pacific Ocean, the MSL change rate has been shown to be 3 to 4 times higher than the global rate for the period 1993-2011, thanks to the exploitation of altimetric data (Zhang and Church 2012, Frankcombe et al. 2015). Regional change of MSL results from a complicated combination of several climatic factors at different timescales and geographical distributions (Hughes and Williams 2010, Stammer et al. 2013, Chen et al. 2018). The El Niño Southern Oscillation (ENSO) is a major global driver of interannual sea level variability, being prominent in the Pacific Ocean and the Indian Ocean (Landerer et al. 2008, Chen et al. 2010, Boening et al. 2012), which during some extreme events show a local rise in level of approximately 30 cm (Becker et al. 2012, McGregor et al. 2012, Widlansky et al. 2014). The Pacific Decadal Oscillations (PDO) also affect sea level variability, especially in the North Pacific region, where this climate fluctuation could modulate the sea level by 10 cm in some regions (Hamlington et al. 2013, Luu and Tkalich 2014, Frankcombe et al. 2015). Qualifying the regional MSL change may also have an alternative social implication, such as determining land areas subject to a sovereignty claim in the South China Sea (Lyons et al. 2018).

While many efforts have been made to derive the globally averaged MSL rise rate linked to global warming (e.g. Dangendorf et al. 2015, Slangen et al. 2016, Wu et al. 2017), its regional patterns have received less attention. The separation of regional trend and variability required a sufficiently long period for two reasons. Firstly, climate oscillations may increase the MSL by 50 mm or higher within a year, which was significantly greater than the annual increase (<5 mm) linked to global warming (Han et al. 2010, Tkalich et al. 2013, Luu et al. 2015). Secondly, each climate factor has a different spatial impact on the Earth’s oceans. Many works (Zhang and Church 2012, Palanisamy et al. 2015) attribute a large part (>12 mm year^–1) of the MSL rise observed in the western tropical Pacific to climate variability. In a recent analysis using climate model ensembles, Fasullo and Nerem (2018) showed that climate forcing associated with the ENSO and the PDO contribute significantly to the spatial patterns of global MSL rise.

In this study, we extend the work of Zhang and Church (2012) using satellite altimetry data in three aspects. First, our domain consists of not only the Pacific Ocean, but also the Indian and Atlantic oceans. Second, we took four dominant modes of sea level variability into consideration instead of two modes, accompanied by a better statistical model to correct for autocorrelation. Lastly, we further considered the lagged times to represent the delayed response of climate impact on sea level.

Monthly sea level data from the Ssalto/Duacs delay-time product provided by the Archiving, Validation and Interpretation of Satellite Oceanographic Data service (AVISO, http://www.aviso.altimetry.fr/en/data.html) for the period 1993-2015 were obtained. The product was reconstructed from ten satellite missions having fine spatial bins and delivering data at a grid resolution of 1/4°×1/4°. We removed the seasonal cycle, applied a 5-month moving average for time series and smoothed spatial patterns to a resolution of 1°×1° covering the domain 0-360°, 50°S-60°N. To correlate with climate fluctuations, we decomposed the time series of sea level at each grid point into decadal and interannual components, as suggested by previous studies (Vimont 2005, Zhang and Church 2012). The decadal dataset was obtained by applying 25-month moving averages followed by another 37-month smoothing on the original time series at each grid point, while the interannual dataset was achieved by subtracting the decadal components from the time series, as suggested by Zhang and Church (2012). For discussion, we further used the sea level data provided by the Commonwealth Scientific and Industrial Research Organization (CSIRO), which combined data from four satellite measurements TOPEX/Poseidon, Jason-1, Jason-2/OSTM and Jason-3, available from http://www.cmar.csiro.au/sealevel/sl_data_cmar.html).

We used different climate indices to associate with dominant modes of sea level variability. The highest correlated forcings are the ENSO, the Central Pacific ENSO (CP ENSO) and the PDO. The multivariate ENSO index (MEI; http://www.esrl.noaa.gov/psd/enso/mei/) was chosen as a proxy of the ENSO (Wolter and Timlin 1998) that has a strong correlation with interannual changes in sea level (Luu et al. 2015, Wu et al. 2017). The PDO index was adopted (from http://research.jisao.washington.edu/pdo/) alongside CP ENSO, as an alternative climate factor akin to the ENSO (Kao and Yu 2009, available at http://www.ess.uci.edu/~yu/2OSC/). The solar radiation data used by Luu et al. (2018) are not considered in this study.

To ease the autocorrelation problem arising from the original least square fitting, the first-order autoregressive and first-order moving average (ARMA(1, 1)) was applied instead of the commonly used first-order autoregressive (AR(1)) model. As pointed out by Foster and Brown (2015), its advantage over the AR(1) model is its capability to resolve the underestimation of standard errors. All statistical values in our study were computed for a two-tailed Student t-test at a 95% significance interval. Note that the corrected confidence intervals are associated with smoothed data and might consist of unavoidable dominant uncertainties, including observational errors in the satellite orbit corrections applied (Ablain et al. 2015).

The empirical orthogonal functions (EOF) technique has been applied successfully in various climate studies (e.g. Bjornsson and Venegas 1997, Church and White 2011). Its advantage is its ability to decompose the spatio-temporal signals into dominant modes associated with spatial patterns and time series by means of orthogonal basis functions. We used the PCAtool developed in Matlab software by Guillaume Maze (given at https://au.mathworks.com/matlabcentral/fileexchange/17915-pcatool) to implement a regular (non-rotating) EOF analysis on sea level variability. Four dominant modes were detected from this analysis: the leading mode from decadal time series and three modes from interannual signals. They were then associated with well-identified climatic drivers that have high correlations, as described below.

For the decadal dataset, the first three leading EOF modes explained 61.2% of total variance. Accounting for 39.8% of the signals, the most dominant mode (EOF1-D) exhibited seesaw patterns in the Pacific Ocean extending symmetrically to mid-latitudes, and in the Indian Ocean (Fig. 1A). Positive values were observed in the eastern Pacific region (180-100°W; 15°S-20°N) and the western Indian Ocean (50-90°E; 0-20°S), while negative values were observed near the western Pacific-eastern Indian areas (110-150°E; 20°S-20°N) and in the central North Pacific (160°E-140°W; 20-40°N). Its temporal evolution was highly correlated (Pearson’s correlation coefficient, r=0.92) with the fluctuations of the low-pass-filtered PDO index (Fig. 1B). We attributed this mode to the PDO, which was consistent with previous findings (Zhang and Church 2012, Frankcombe et al. 2015).

Three leading modes revealed from the interannual dataset explained 39.9% of total variance. The largest principal mode (EOF1-I) accounted for 26.8% of net signals and was related to the ENSO influence due to a significant correlation (r=0.94) between spatio-temporal patterns of this mode and the high-pass-filtered ENSO index (Fig. 2B). For example, in the Pacific Ocean, a narrow seesaw pattern was found, comprising a negative anomaly in the west of the tropical equator (120-170°E; 5°S-15°N) and a positive anomaly extending from the middle of the Pacific toward the western coast of the American continent (180-65°W; 10°S-10°N). The results were similar to the ENSO-induced sea level in the tropical Pacific Ocean revealed in other studies (e.g. Hamlington et al. 2011, Widlansky et al. 2014). In contrast, a high anomaly was depicted around the western part of the Indian Ocean (40-100°E; 10°S-10°N), while a negative one was detected in its eastern part. These patterns were linked to the ENSO in earlier studies (Hamlington et al. 2011, Zhang and Church 2012).

The second principal mode (EOF2-I) showed a smaller contribution to the interannual variance (8.8%) and had a complicated structure. In the Pacific, it comprised a narrow negative equatorial belt bounded between 5°S and 5°N, a positive bar mirrored over the latitudinal line of 10°N, and a significantly negative trace appearing in the South Pacific Convergence Zone (SPCZ). Adoption of the MEI yielded the best correlation of 0.90 for a 7-month lag (Fig. 2D). This half-year lag was consistent with the findings of Widlansky et al. (2014), who observed an interaction between sea level in the northwest and southwest Pacific tropical regions and the ENSO extreme events. Meanwhile, the third principal mode (EOF3-I) depicted a strong positive anomaly in the Niño-4 region (5°S-5°N, 160°E-150°W) shown in Fig. 2E. This anomaly was previously reported by Kug et al. (2009), who pointed to the CP ENSO influence. The highest correlation between the temporal evolution of the CP ENSO and the mode was 0.56 after a 2-month lag (Fig. 2F and Table 1), contributing 4.3% to total variance.

The time lags between climate indices and principal modes of sea level variability are summarized in Table 1. To further examine the sensitivity of these time lags, we repeated the EOF analysis on the sea level dataset provided by CSIRO and obtained similar results (Table 1).

Regional changes in sea level over a 20-year period were partitioned by taking advantage of high-resolution satellite altimeter data. To determine the components of sea level variability, its dominant modes were identified from interannual and decadal time series by means of an EOF method, which linked three leading interannual modes to the ENSO and the CP ENSO, and the principal decadal mode to the PDO. MVLR models were used to separate the trend and variability, which showed an improvement versus the SVLR in two statistical measures. Firstly, the use of the multivariate model improved the coefficients of determination for two datasets used, namely AVISO and CSIRO. Secondly, the MVLR model also significantly eased the statistical uncertainty in estimating the sea level change rate. The MVLR model was then employed to display the patterns and temporal variability of regional sea level change.

Significant adjustments were visible in the Pacific Ocean, where the positive change in sea level (>10 mm year^–1) shifted from its western basin to the tropical areas (150°E-125°W; 20°S-20°N), while the negative rate on its eastern side diminished. The rise became stronger in some areas of the Indian Ocean (around 15°S) and Atlantic tropics (0-60°W; 0-20°N). Along the meridional axis, the largest discrepancy was observed in the Pacific Ocean, which is mainly contributed by the PDO. The MSL rose faster in the Southern Hemisphere than in the Northern Hemisphere. In the temporal change, the ENSO exhibited prevailing impacts in interannual sea level variations in all oceanic basins. Despite being insignificant in the long term, the CP ENSO made a noticeable contribution at the interannual scale, for instance, over the period 1997-1998.

We derived the regional patterns of sea level change by means of a statistical approach using satellite data. Our results are consistent with the conclusion from a recent study by Fasullo and Nerem (2018), who used climate model ensembles and came to the similar conclusion that external climate forcing including the ENSO and the PDO is probably a main contributor to the observed non-uniform patterns of sea level rise in the global oceans. They further suggested that these patterns might remain for decades, and will probably be intensified under the changing climate.

For future studies, the fairly low significant level in some regions and the small variances explained by interannual modes suggested that improvements are desirable. These may include the introduction of other climatic and non-climatic factors and a more suitable regression model, as there are still a lot of interannual and decadal variabilities left after applying the MVLR model, which alias into the trend estimates. In fact, this is a great challenge, since a quantitative relationship between regional sea level change and climate fluctuations is undeterminable outside the context of the Earth’s sophisticated climate system, which is driven by both deterministic and stochastic processes. As a result, further continuing efforts are needed to give deeper insights into not only the statistical decomposition technique to attain the rate, but also the establishment of non-deterministic mathematical relationships between sea level variability and other physical components of the climate system, including ice sheet melts and volcano eruptions (Slangen et al. 2016, Marcos et al. 2017, Huang et al. 2018).