Chen et al. (2017)

Title:

The increasing rate of global mean sea-level rise during 1993–2014

Key Points:

The global mean sea level (GMSL) rose from 2.2 ± 0.3 mm yr−1 in 1993 to 3.3 ± 0.3 mm yr−1 in 2014, indicating an acceleration in sea-level rise during this period.
The mass loss from the Greenland Ice Sheet significantly increased during this period, accounting for less than 5% of the GMSL rate in 1993 but more than 25% in 201
The mass contributions to GMSL, which include glacier mass loss, Greenland and Antarctic ice sheets, and anthropogenic terrestrial water storage, increased from about 50% in 1993 to 70% in 2014, indicating their growing role in sea-level rise.

Keywords:

Global mean sea-level rise, acceleration, Greenland Ice Sheet, Mass contributions, Satellite altimetry

Corresponding author:

Xianyao Chen

Citation:

Chen, X., Zhang, X., Church, J. A., Watson, C. S., King, M. A., Monselesan, D., et al. (2017). The increasing rate of global mean sea-level rise during 1993–2014. Nature Climate Change, 7(7), 492–495. doi:10.1038/nclimate3325

URL:

https://www.nature.com/articles/nclimate3325

Abstract

Global mean sea level (GMSL) has been rising at a faster rate during the satellite altimetry period (1993–2014) than previous decades, and is expected to accelerate further over the coming century1. However, the accelerations observed over century and longer periods2 have not been clearly detected in altimeter data spanning the past two decades3,4,5. Here we show that the rise, from the sum of all observed contributions to GMSL, increases from 2.2 ± 0.3 mm yr−1 in 1993 to 3.3 ± 0.3 mm yr−1 in 2014. This is in approximate agreement with observed increase in GMSL rise, 2.4 ± 0.2 mm yr−1 (1993) to 2.9 ± 0.3 mm yr−1 (2014), from satellite observations that have been adjusted for small systematic drift, particularly affecting the first decade of satellite observations6. The mass contributions to GMSL increase from about 50% in 1993 to 70% in 2014 with the largest, and statistically significant, increase coming from the contribution from the Greenland ice sheet, which is less than 5% of the GMSL rate during 1993 but more than 25% during 2014. The suggested acceleration and improved closure of the sea-level budget highlights the importance and urgency of mitigating climate change and formulating coastal adaption plans to mitigate the impacts of ongoing sea-level rise.

Main

Projections of future sea levels must be based on a sound understanding of historical changes in GMSL and its underlying processes, as well as recent changes in the rate of rise1. In a previous study, the apparent decrease in the rate of GMSL rise from 3.2 mm yr−1 in the first decade of satellite altimetry to 2.8 mm yr−1 in the second was suggested to be primarily a result of natural interannual variability, related to water exchange between ocean and land during El Niño/Southern Oscillation (ENSO) cycles3. After removing this variability, the underlying rate of GMSL rise was 3.3 ± 0.4 mm yr−1 for both decades, with neither deceleration nor acceleration of GMSL inferred over 1993 to 2014. This lack of observed acceleration of GMSL contrasts with a simultaneously increased contribution from the Greenland ice sheet (GIS) and a less certain increase from the Antarctica ice sheet (AIS) overall7, and is inconsistent with the positive acceleration presented in century-long tide gauge data8 and global mean sea-level reconstructions2.

By comparing tide gauge and satellite altimeter sea-level observations, a recent study6 identified a possible systematic drift within the altimeter record, particularly affecting the first six years (1993–1999). This systematic error erroneously elevates the GMSL trend during 1993–1998 by between 0.9 ± 0.5 and 1.5 ± 0.5 mm yr−1, depending on whether a glacial isostatic adjustment (GIA) or Global Positioning System (GPS) data set was used to correct for the effects of land motion at tide gauges used in the bias estimation process. After removing these biases, the estimated rate of GMSL rise from 1993 to mid-2014 was between 2.6 ± 0.4 and 2.9 ± 0.4 mm yr−1, with a positive but not statistically significant acceleration of 0.041 ± 0.058 mm yr−2, compared with a not statistically significant deceleration of −0.057 ± 0.058 mm yr−2 for unadjusted data.

GMSL rise results from the ocean thermal expansion, loss of mass from glaciers9, the GIS and the AIS7, and changes in land water storage from climate variability and anthropogenic effects10,11. To study how the rate of the GMSL rise varies during the satellite period, we investigate the time-varying intrinsic trend in GMSL and these contributing components by separating them from their interannual variability using an adaptive data analysis approach.

An intrinsic trend is defined as ‘an intrinsically fitted monotonic function or a function in which there can be at most one extremum within a given data span’12. Unlike the commonly used linear polynomial trend that requires a priori assumptions regarding stationarity and linearity of time series, the intrinsic trend is not defined by a predetermined functional form of the trend and, hence, is more adaptive to the underlying physical properties of observations. The method we apply to derive the intrinsic trend in GMSL and its components is ensemble empirical mode decomposition (EEMD)13, which is based on the empirical mode decomposition (EMD) method designed for adaptive analysis of nonlinear and non-stationary time series14 (see Methods and Supplementary Information) and has only recently been applied to sea-level trend estimation8. The main benefit of using EMD is that it can separate non-stationary oscillations (such as natural variations on different timescales) from the long-term trend, and the trend is found empirically without any assumptions about its shape.

Figure 1a presents the unadjusted GMSL from four different processing groups, and the adjusted CSIRO GMSL derived using vertical land motion (VLM) estimates at tide gauges based on GIA and GPS6. We note there is not yet a homogeneously reprocessed altimeter data set that addresses the likely systematic bias estimated by ref. 6. In the absence of such a data set, an assessment of all available records (including the adjusted record from ref. 6) is entirely appropriate. The intrinsic trend and interannual variability of each time series are shown in Fig. 1a, b, respectively. The significance of the intrinsic trend is tested against a null hypothesis of red background noise with the same lag-1 autocorrelation as the raw time series (Supplementary Fig. 1), and it is shown that the increase in the intrinsic trend of the GPS-based adjusted GMSL record is statistically significant during the recent decade (Supplementary Fig. 2).

Figure 1: Global mean sea level (GMSL). Curves show the GPS-based and GIA-based adjusted GMSL, and unadjusted GMSL from four different groups. a, GMSL and the time-varying secular trend from EEMD analysis. b, The interannual variability of GMSL. c, The instantaneous rate of GMSL rise. The uncertainties of the derived interannual variability, intrinsic secular trend, and its instantaneous rate are shown in coloured shades.

We derive the rate of the time-varying intrinsic trend by calculating its first-order temporal derivative (Fig. 1c). Consistent with previous studies3, the unadjusted GMSL exhibits a slightly decreasing rate of rise from about 3.5 mm yr−1 during the first decade to 3.0–3.3 mm yr−1 during the second. In contrast, the rate of the GPS-based adjusted GMSL rise increases by 0.5 mm yr−1 from about 2.4 ± 0.2 (1σ) mm yr−1 in 1993 to around 2.9 ± 0.3 mm yr−1 in 2014 (2.8 ± 0.2 to 3.2 ± 0.3 mm yr−1 for the GIA-based adjusted GMSL). That is, the time-varying trend of the adjusted altimeter data suggests an acceleration in GMSL in agreement with ref. 6, with the dominant increase in the rate of rise occurring in the recent decade.

Figure 1b shows the interannual variability of all GMSL records derived by EEMD. It includes a significant drop of water level related to the transfer of water from the ocean to the land during the strong La Niña event in 2011 and the subsequent rapid recovery in the following two years15. This interannual variability agrees well with the interannual variability of the land water storage based on the global hydrological model16 and the interannual variability of the thermosteric sea level17, and is significantly correlated (0.42 for the unadjusted CSIRO-based time series, and up to 0.56 for the others) with the ENSO index18. By definition of the EMD method, this interannual variability does not contribute to the intrinsic trend over the whole period.

We determined global mean steric sea-level (GMSSL) anomalies using a range of subsurface measurements of temperature and salinity data sets19. The GMSSLs from these data sets exhibit wide discrepancies owing to the inhomogeneous observations in the ocean, different data quality control procedures, XBT (expendable bathythermograph) bias corrections, mapping methods and model structures20. These discrepancies are especially pronounced until 2005 when sufficient spatial data coverage was obtained from Argo floats. We select ocean temperature–salinity data sets that do not have obvious discontinuities in the GMSSL time series and whose linear trend of GMSSL during the Argo period (2005–2014) remains within the 2σ range of that derived from three Argo gridded data sets (Supplementary Table 1). Figure 2 shows monthly GMSSL anomalies, interannual variability and the intrinsic trends from seven data sets based on the above selection criteria.

Figure 2: Global mean steric sea level (GMSSL). Coloured curves show the global mean steric sea level from seven data sets. a, The GMSSL and its time-varying secular trend. b, The interannual variability of GMSSL. c, The instantaneous rate of GMSSL rise. The uncertainties of the derived interannual variability, intrinsic secular trend, and its instantaneous rate are shown in coloured shades. In c, the dots denote the median of all GMSSL records at each year, with the uncertainty estimated using the median statistical method. Note that the different length of the GMSSL time series affects the median rates over the last few years, and consequently affects the budget over the last few years as shown in Fig. 4.

Even with these relatively strict criteria, the GMSSLs of the selected ocean temperature–salinity data sets still exhibit remarkable differences over the whole period. The instantaneous rate of the GMSSL of some models indicates acceleration, whereas others not. To reduce the impact of the skewness, we estimate the instantaneous rate of GMSSL rise as the median of that derived from seven data sets at each year. The derived mean thermal expansion contribution is about 0.94 ± 0.16 mm yr−1 during 1993–2014, which is equivalent to about 0.48 ± 0.08 W m−2 net surface heat flux into the ocean, and consistent with the observed top-of-atmosphere heat imbalance21. The ensemble estimate of the GMSSL rise suggests little acceleration during the satellite altimetry period.

The main contributions to the global ocean mass changes are from the GIS, the AIS and glaciers. The GIS and AIS mass changes are investigated using the estimates based on altimetry, gravimetry and mass flux data for 1993–20127, and the GRACE observations during 2003–2014 by adjusting its trend to match the published data over 2003–200922 (Supplementary Fig. 5). The glacier data are estimated from a glacier mass balance model driven by gridded climate observations9.

Figure 3 shows that all three sources of mass loss exhibit an increasing contribution to GMSL. The rate of glacier mass loss increased over 1993 to 2005, from 0.60 ± 0.15 to 0.87 ± 0.21 mm yr−1 GMSL equivalent, but is then nearly unchanged up to 2013 (Fig. 3c). The GIS mass loss increased from around 0.11 ± 0.03 mm yr−1 in 1993 to around 0.85 ± 0.03 mm yr−1 in 2014, approaching an average acceleration of 0.03 mm yr−2. The rate of the AIS mass loss is around 0.22 ± 0.02 mm yr−1 in 1993, and only slightly increases to 0.31 ± 0.02 mm yr−1 in 2014. These trends agree quantitatively with previous linear estimates7,9 over the whole satellite period, and contribute to the acceleration of GMSL.

Figure 3: Global mean ocean mass change. Curves show the land-glacier, AIS and GIS and anthropogenic TWS contributions to GMSL. a–c, Each mass component and its secular trend (a), the interannual variability (b), and the instantaneous rate (c) of the ocean mass change. In c, the black dots and their error bars show the rate of thermal expansion, GMSSL, from Fig. 2 for comparison with the mass change rate. The uncertainties of the derived interannual variability, intrinsic secular trend, and their instantaneous rates are shown in coloured shades.

Another contribution to changes in the global ocean mass is from terrestrial water storage (TWS), including that associated with anthropogenic activities (groundwater extraction, irrigation, impoundment in reservoirs, wetland drainage, and deforestation) and natural climate variability. Here, the anthropogenic TWS changes are based on the estimates of ref. 11, with their groundwater depletion being replaced with the estimates of ref. 10, which are 20% smaller. This smaller estimate is consistent with 80% of the extracted ground water making its way to the ocean23. The intrinsic trend and its instantaneous rate of the anthropogenic TWS show a slightly increased contribution to GMSL from around 0.11 ± 0.04 mm yr−1 during the first decade to about 0.24 ± 0.06 mm yr−1 during the second (Fig. 3c).

Regarding the natural variability of TWS, global values are not reliable before the GRACE mission in 2002. Interannual fluctuations of TWS based on the continental water balance model are estimated as about 0.25 mm yr−1 (GMSL equivalent) during 1993–199824, whereas the GRACE observations during 2002–2012 suggested a natural TWS contribution to GMSL of around −0.71 ± 0.20 mm yr−1 (ref. 25). This rate is approximately consistent with the 5.5 mm fall in GMSL over 2002–2012 in the interannual variability (Fig. 1b), which is highly correlated with the La Niña-like variability in the Pacific15, when precipitation decreases over the ocean and increases over the land. Because of the strong ENSO-related interannual variability, there can be significant trends in TWS over periods of a decade or shorter. Therefore, short-period linear trend estimates do not adequately represent the time series over the whole satellite altimeter period, and it is likely that the total trend is small (but poorly quantified).

Using time series of GMSL, GMSSL and all components of global ocean mass change, Fig. 4 shows the instantaneous budget of GMSL over the satellite period. The thermal expansion component is about 50% of the total contributions in 1993. Although the rate of this contribution did not change much throughout the record, by the end of the record it is reduced to about 30% of the sum of the contributions because of the acceleration in the global ocean mass component, consistent with a previous estimate of the changing relative roles of ocean thermal expansion and ocean mass26. The ocean mass change is initially dominated by the contribution of glacier mass loss, with smaller contributions from the GIS and AIS mass loss and anthropogenic TWS changes. But in the recent decade, the acceleration of the mass loss from the GIS was the largest, and its contribution to the GMSL became almost equal to that from thermal expansion and glaciers by 2014. The year-by-year contribution from the AIS mass loss is nearly constant while the glacier contribution increases slowly.

Figure 4: Instantaneous closure of the global mean sea-level budget. Yearly instantaneous rate of change of the GPS-based adjusted (black dots) and mean unadjusted GMSL (grey stars) and that of the GMSSL, and ocean mass contributions from the GIS, the AIS, the anthropogenic TWS and glaciers, and ocean thermal expansion (each shown in coloured shades, ordered from top to bottom). The blue dots denote the sum of the instantaneous rates of change of each component with its uncertainty estimated as the square root of the sum of the squares of the uncertainty in each instantaneous rate, as shown in previous figures. The time series of the loss of mass from the glaciers and the anthropogenic TWS stops in 2012, and 2009, respectively. Their rates in the years up to 2014 are assumed unchanged and shown in a lighter colour.

In all future projection scenarios of the Fifth Assessment Report of the Intergovernmental Panel on Climate Change27, the largest contribution to changes in GMSL is the ocean thermal expansion, accounting for 30–55% of the projection, whereas the glaciers are the second largest, accounting for 15–35%. Our analysis of recent observations shows that the acceleration of ocean thermal expansion during 1993–2014 is not significant. Climate model simulations indicate the fall in ocean heat content following the 1991 volcanic eruption of Mount Pinatubo and the subsequent recovery has probably resulted in a rate of thermal expansion about 0.5 mm yr−1 higher than would be expected from greenhouse gas forcing alone28. The recovery in ocean thermal expansion following major volcanic eruption takes more than 15 years28,29. Thus, the underlying acceleration of thermal expansion in response to the anthropogenic forcing may emerge over the next decade or so, resulting in a further acceleration in the rate from that reported here and recent estimates30. In contrast to the lack of observed acceleration in the ocean thermal expansion, there has been a significant acceleration in the mass contributions, dominated by the increased GIS mass loss. This results in an approximate closure of the sea level budget throughout the study period from 1993 to 2014 and, importantly, both the sum of contributions and the GPS (and GIA)-based adjusted altimeter rates indicate an acceleration in sea level over the satellite altimeter period.

This approximate but improved closure of the sea-level budget throughout 1993–2014 is progress with respect to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, and increases confidence in our observations and understanding of recent changes in sea level. The study period is still short and ongoing observations are required to understand the longer-term significance of this finding, and to identify the contributions of decadal and multi-decadal variations that are unresolved in the 20-year-long records. The estimated increase in the rate of rise has important implications for projections of sea-level rise and for society.

Methods

Data

Satellite altimetry

We use six different altimetry-based monthly sea-level data from four processing groups: Archiving Validation and Interpretation Satellite Oceanographic Center (AVISO; https://podaac.jpl.nasa.gov/dataset/MERGED_TP_J1_OSTM_OST_GMSL_ASCII_V4); Colorado University (CU, Release 4; http://sealevel.colorado.edu/files/2015_rel4/sl_ns_global.txt); Goddard Space Flight Center (GSFC; https://www.aviso.altimetry.fr/en/data/products/ocean-indicators-products/mean-sea-level.html); Commonwealth Scientific and Industrial Research Organization (CSIRO) and University of Tasmania (ftp://ftp.marine.csiro.au/pub/legresy/gmsl_files/CSIRO_Alt_refined.csv); adjusted CSIRO data set using a model of glacial isostatic adjustment to estimate vertical land motion at tide gauges; and adjusted CSIRO data set using GPS data to estimate vertical land motion at tide gauges. All six data sets are based on TOPEX, Jason-1 and OSTM/Jason-2 data. The global average is computed over 66° S–66° N for AVISO, CU and GSFC, but over 65° S–65° N for CSIRO. Detailed descriptions of each data set are available in the corresponding websites.

Steric sea-level data sets

We use seven products describing the global ocean monthly temperature–salinity to compute the global mean steric sea level. Two of these are objective analyses based on optimal interpolation without constraints from ocean model dynamics, and five are reanalyses based on data assimilation with models. From 20 global ocean temperature–salinity data sets [19], we select a subset of seven ocean temperature–salinity data sets that do not have obvious discontinuities in the GMSSL time series and whose linear trends of GMSSL during the Argo period (2005–2014) remain within the 2σ range of that derived from three Argo gridded data sets. Supplementary Table 1 provides the basic information of these data sets.

Land glaciers

The global yearly glacier mass data set used in this paper is produced with a glacier model driven by gridded climate observations [9].

Greenland and Antarctic ice sheets

Greenland and Antarctic ice sheet records during 1993–2012 available at http://imbie.org/data-downloads [7] are used in this study. To extend the records to the end of 2014, observations based on the GRACE satellite are used. Noting the potential GIA error in GRACE, especially for Antarctica31, we adjusted the trend of the GRACE records so that they agreed with the published trends over 2003–2009. Then two records are connected in 2003. After 2003, the GRACE record is used. The derived monthly time series are shown in Supplementary Fig. 5.

Anthropogenic terrestrial water storage

The yearly data of anthropogenic terrestrial water storage are extracted from the century-long time series [11], after their groundwater depletion is replaced with the most recent estimates [10].

Ensemble empirical mode decomposition method

The main method used in this study is ensemble empirical mode decomposition (EEMD) [13], which was developed on the basis of the empirical mode decomposition (EMD)14 method. The EMD and EEMD methods have been applied to oceanic and climatic time series analysis32, and are also used to study regional [8,33,34] and global sea-level variability35. Here we briefly introduce the general decomposition procedure and mainly introduce the tests used to assess statistical significance and the estimation of the uncertainty of the intrinsic secular trend derived by the EMD/EEMD method.

Decomposition and intrinsic secular trend

In EMD, a time series x(t) is decomposed into a set of amplitude–frequency-modulated oscillatory functions (so-called intrinsic mode functions, IMFs) Cj(t), j = 1,2, … n and a residual R(t): x(t) = ∑ j=1nCj(t) + R(t) through a sifting process. The following is undertaken: (1) locate all maxima and minima and connect all maxima (minima) with a cubic spline as an upper (lower) envelope of the time series; (2) compute the difference between the time series and the mean of the upper and lower envelopes to yield a new time series h(t); (3) for the time series h(t), repeat steps 1 and 2 until upper and lower envelopes are symmetric with respect to zero mean under the stopping criteria13,14, then an IMF, Cj(t), is derived as the time series h(t); and (4) subtract Cj(t) from original time series to yield a residual R(t) and treat R(t) as the original time series and repeat steps 1–3 until the residual R(t) becomes a monotonic function or a function with only one extremum; then, the whole sifting process is completed and all IMFs and the residual function, namely, the intrinsic trend of x(t), are obtained.

The EEMD approach is based on EMD13. In EEMD, multiple noise realizations are added to the time series x(t), from which an ensemble average of the corresponding IMFs is extracted to yield scale-consistent signals. The main steps in the EEMD are as follows: (1) add a white noise series to the targeted data; (2) decompose the data with the added white noise into IMFs; (3) repeat step 1 and 2 again and again, but with different white noise series each time; and (4) obtain (ensemble) means of the respective IMFs of the decompositions as the final result. The advantage of the EEMD is that by using an ensemble mean, non-physical oscillations due to random data errors are reduced and thus low-frequency modes are more accurate.

It is proved that the added white noise with the variance σ will have at most impact on the resulting IMFs13, where N is the number of ensemble members. When N increases, this impact is negligibly small. In this paper, we always use the white noise with variance σ = 0.2 relative to the variance of the original time series, and N = 1,000 ensemble members.

As demonstrated in previous applications of the EEMD method involving trend analysis of global mean surface temperature36, global land surface air temperature37, and the sea-level observations along the eastern US coast8,33, the residual function R(t) is derived by removing any, but not a predetermined, variability on shorter timescales than the length of time series, as represented by the IMFs Cj(t). Consequently, R(t) can take any unspecified shape and will preserve the potential variability on longer timescales than the length of time series. For the observations used here (22 year duration for the altimetry), the relatively high-frequency variability on the interannual timescales will be shown by summing all of the IMFs, that is, ∑ j=1nCj(t) and the residual R(t) will be regarded as the intrinsic trend.

It should be noted that the intrinsic trend R(t) is in the same unit as the raw time series (in this case for the global mean sea level used here, the unit is millimetres). Taking the first-order time derivative of the time-varying intrinsic trend yields the instantaneous rate of the trend, in units of millimetres per year for global mean sea-level rise, which provides more time-varying information on how the intrinsic trend has evolved within the given time region, compared with a typical fitted polynomial and the time-varying estimations based on the sliding window approach38,39.

The validity of the intrinsic trend is strongly based on the specified data duration. The properties of the trend beyond the data length require further investigation with more observations. In this study, the potential decadal variability of GMSL on the timescale longer than the length of satellite altimetry cannot be separated from the secular trend, which implies that the accelerated GMSL may partially reflect the internal decadal variability, as well as the effects of the anthropogenic forcing.

Significant test of the intrinsic secular trend

To test the statistical significance of the intrinsic secular trend, one needs to reject a null hypothesis that it arises by chance for stochastic processes with zero means at given confidence levels. In climate sciences, two widely used null hypotheses are that the underlying processes are noise characterized as white (that is, no temporal correlation) or red (with lag temporal correlation). There are many methods to test the statistical significance of a linear, curve-fitted, or time-varying trend against a white or red noise null hypothesis33,40,41,42. Here we applied one approach developed for testing the time-varying trend derived by the EMD/EEMD method37. Although the detail of the statistical significance test is given in their Supplementary Section 2, the general procedures of this approach are introduced as follows, in order for the integrity of this work.

For any time series x(t) with time-varying secular trend R(t) derived by the EEMD method, the statistical significance test includes: (1) computing the lag-1 autocorrelation α of the time series x(t). If α = 0 then the null hypothesis that the white background noise is applied; if α > 0, then the null hypothesis that the red background noise is applied; (2) generating 5,000 samples of red noise time series with the same temporal data length of x(t) and the lag-1 autocorrelation α; (3) deriving the intrinsic trend of each generated red noise time series by using the EEMD method. This yields an empirical probability density function of the intrinsic trends, which is approximately normally distributed, at any temporal locations; (4) comparing the intrinsic trend of the studied time series with the two-standard-deviation spread value of the trends of the red noise time series (around 95% of confidence) at any temporal locations. If the former is larger, the intrinsic trend of the studied time series is considered statistically significant and the null hypothesis that the intrinsic trend of the time series is from noise could be rejected; (5) taking the first-order time derivative of the intrinsic trend yields its instantaneous rate. If the intrinsic trend is statistically significant, we will consider its instantaneous rate is significant.

In this approach, the noise time series are generated on the basis of the AR(1) process x(t) = αx(t − Δt) + w(t), where the coefficient α is the autoregressive parameter (that is, lag-1 correlation coefficient), Δt is the sampling rate and w(t) is the white noise with the unit standard deviation. A more general case is to generate the red noise on the basis of the ARMA(p, q) process43,44. Empirically, the selection of the red noise model may change the probability density function shape of the EEMD trends of each red noise time series, but it will not change the general patterns of the statistical significance test if we go through all possibilities of the lag-1 autocorrelation from 0 to 0.99. In this study, we adopt the AR(1) model.

Taking the GPS-based adjusted GMSL time series as an example, Supplementary Figs 2–4 present the analysis of its intrinsic trend, and the test of its statistical significance. Supplementary Fig. 2 shows that the lag-1 autocorrelation coefficient of the GPS-based adjusted GMSL is α = 0.84. With this lag-1 autocorrelation coefficient, 5,000 AR(1) time series are generated and then decomposed by the EEMD method to derive their intrinsic trend (Supplementary Fig. 2a). Both thick black lines in Supplementary Fig. 2a are two-standard-deviation spread lines of the trend of AR(1) time series. Note that these intrinsic trends of the background noise are dimensionless. To test the statistical significance of the intrinsic trend of GPS-based adjusted GMSL time series, which is in the same unit of millimetres, we divide it by the standard deviation of the linearly detrended GMSL time series (red line in Supplementary Fig. 2a) and compare it with two standard deviations of the intrinsic trend of the red background noise at each temporal location. If the former is larger, the trend of GMSL is considered statistically significant. Supplementary Fig. 2b, c presents the comparison at two randomly selected years 1999 and 2009, respectively, in which the trend of GMSL time series is not statistically different from the trend of AR(1) red noise with 0.84 lag-1 autocorrelation in 1999, but significantly different in 2009. In Supplementary Fig. 2a, comparing the temporal variability of the intrinsic secular trend of the GPS-based adjusted GMSL time series with the two-standard-deviation line (around 95% confidence) of the intrinsic trend of the 5,000 AR(1) processes with 0.84 lag-1 autocorrelation shows that the intrinsic trend of GPS-based adjusted GMSL becomes statistically significant since 2005.

As for the definition of the intrinsic trend12, this statistical significance test is also valid only during the studied period, because the properties of the intrinsic trend beyond the data length may change a lot, and so their statistical significance.

To test the time series with different lag-1 autocorrelation, Supplementary Fig. 3 shows the two-standard-deviation line (around 95% confidence) of the intrinsic secular trend of the 5,000 AR(1) processes with the lag-1 autocorrelation ranging from 0 to 0.99. The trend spread depends on the value of lag-1 autocorrelation. When the noise is getting redder (larger lag-1 autocorrelation), the corresponding spreads become wider37.

Based on the 95% confidence levels of different AR(1) time series with lag-1 autocorrelation ranging from 0 to 0.99, the statistical significance of the intrinsic secular trend of the GMSL, the GMSSL, and the global ocean mass change can be tested, as shown in Supplementary Fig. 4a–c, respectively.

Estimation of the uncertainty of the intrinsic trend

To estimate the uncertainties in the intrinsic trend, the down sampling approach36 is applied. This method is also used to study the increasing flooding hazard in Miami Beach, Florida [34].

For the monthly time series of global mean sea level, global mean steric sea level, or the Greenland and Antarctic ice sheet mass loss, we randomly pick a value of the time series for each calendar year to represent the entire annual average. This step can theoretically yield 12 (ref. 22) different time series for 22 years of monthly data. We randomly choose 1,000 series and re-compute their intrinsic secular trend, and then obtain the mean of the trend, and the spread of the trends provides the confidence interval. For the yearly time series of glacier and anthropogenic terrestrial water storage time series, we randomly pick up a value of the time series within two standard deviations of the time series to represent the spread of the time series, and choose 1,000 different series and re-compute their intrinsic secular trend. The rate of the intrinsic trend is obtained by computing the mean of the time derivative of each intrinsic trend and its uncertainty. This approach can also be applied to estimate the uncertainty of the other IMFs of the time series if the timescale of the function is longer than a month.

Since the higher-frequency variability of the time series is gradually separated using EEMD, the uncertainty of the intrinsic trend is generally less than the uncertainty estimation of the linear trend of original whole time series. The uncertainty of the intrinsic trend at the start and end of the data range is relatively large given the edge effects. This is unavoidable for any temporally local analysis method, such as the Gibbs effect of the Fourier transform and the ‘cone of influence’ of wavelet analysis45. Compared with these methods, the temporal locality of EEMD is smaller, and so are the uncertainties at the two ends, as discussed in the study of uncertainty of the global sea surface temperature [36].

The above uncertainty estimation of the intrinsic trend does not take into account systematic (non-averaging) error terms present in the time series. For the case of the adjusted GMSL time series, uncertainty associated with the bias drift estimation for each mission needs to be considered. The above uncertainty estimation method can be applied to each realization of the time series, with each realization generated sampling the bias drifts and intra/inter-mission bias and associated uncertainties. Applying the EEMD analysis to each time series and estimating the uncertainty of the derived intrinsic trend gives a joint uncertainty estimation of the instantaneous rate of the GPS-based adjusted GMSL rise (not shown). The slightly higher uncertainty of the instantaneous rate of the GMSL change (from ±0.6 mm yr−1 during 1993 to ±0.4 mm yr−1 during 2014) reflects the different mission lengths and the split between TOPEX side A and side B in the early part of the record.