To analyse post-net zero temperature changes, we used several different methods, including Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)18, Locally Weighted Scatter Plot Smoothing (LOWESS), linear splines and a Hodrick-Prescott filter. These methods are applied to a single realisation of the nine ESMs.
CEEMDAN is a signal decomposition technique designed to handle non-stationary and non-linear data. CEEMDAN builds upon Empirical Mode Decomposition (EMD)19 by addressing the limitations of mode mixing and ensuring more accurate extraction of oscillatory modes.
CEEMDAN decomposes temperature time series into intrinsic mode functions (IMFs) and a residual. The residual reflects the longer-term trend in the data, whilst the IMFs represent oscillations across various time scales—from inter-annual to multi-decadal variability (such as ENSO and IPO). The decomposition process is adaptive and data-driven, making no prior assumptions about the data provided. IMFs are extracted iteratively by identifying and removing local extrema, producing oscillatory modes that vary over time.
The steps of EMD involve first identifying the local extrema (minima and maxima) of the signal x(t), and then using cubic spline interpolation to construct the upper envelope \({L}_{1}\left(t\right)\) and lower envelope \({L}_{2}\left(t\right)\). The mean envelope, \(M\left(t\right)\), is then calculated as:
$$M\left(t\right)=\frac{{L}_{1}\left(t\right)+{L}_{2}\left(t\right)}{2}$$
(1)
The difference between the signal and the mean envelope produces the first component:
$$H\left(t\right)=x\left(t\right)-M\left(t\right)$$
(2)
If \(H\left(t\right)\) satisfies the conditions for an IMF, it becomes the first IMF C1(t). Otherwise, the process is repeated until an IMF is obtained. After extracting the first IMF, the residual is calculated as:
$${r}_{1}\left(t\right)=x\left(t\right)-{C}_{1}\left(t\right)$$
(3)
This process is repeated for the residuals until no local extrema can be found, with the leftover series being called the residual (\({r}_{n}\left(t\right)\)), which represents the long-term behaviour or trend. The final result yields additional IMFs and a residual \({r}_{n}\left(t\right)\), so that the signal can be reconstructed as:
$$x\left(t\right)=\mathop{\varSigma }\limits_{i=1}^{N}{C}_{i}\left(t\right)+{r}_{n}\left(t\right)$$
(4)
EMD can suffer from mode mixing, where different frequency components overlap within a single IMF. CEEMDAN enhances this by adding adaptive white noise to each decomposition step. At each iteration, Gaussian white noise \(\epsilon {d}_{i}\left(t\right)\) is added to the signal \(x\left(t\right)\) where \(\epsilon\) is the noise amplitude and \({d}_{i}\left(t\right)\) represents the noise realization. The noisy signal becomes:
$${x}_{i}\left(t\right)=x\left(t\right)+\epsilon {d}_{i}\left(t\right)$$
(5)
IMFs are then calculated for each noisy realization and the mean IMF across all realizations is taken as the final IMF for that mode. The residual \({r}_{1}\left(t\right)\)) is updated as:
$${r}_{1}\left(t\right)=x\left(t\right)-\frac{1}{k}\mathop{\varSigma }\limits_{i=1}^{k}{C}_{i}\left(t\right)$$
(6)
This process is iterated over each mode until all IMFs are obtained. CEEMDAN reduces mode aliasing by ensuring that each IMF represents a distinct oscillatory mode, making it better suited for handling complex climate signals. Endpoint issues are addressed using a reflection method25,26.
As the post-net zero temperature changes in ZECMIP may not be monotonically increasing or decreasing, we estimate the post-net zero changes are estimated by not only using the residual, but also adding IMFs with wavelengths longer than 50 years. IMFs with wavelengths greater than 50 years were isolated using a Fast Fourier Transform (FFT). The FFT is defined as:
$$X\left(f\right)=\mathop{\varSigma }\limits_{n=0\,}^{N-1}x\left(n\right){e}^{-2i\pi \frac{{fn}}{N}}$$
(7)
where \(X\left(f\right)\) represents the frequency-domain signal, \(\left(n\right)\) is the original time-domain signal, M is the total number of points, f is the frequency. By combining FFT with CEEMDAN, the trend was estimated without natural variability.
We then compare this approach to several alternative methods for estimating GMST changes: Locally Weighted Scatterplot Smoothing (LOWESS), which captures non-linear trends by fitting multiple locally weighted regressions; the Hodrick–Prescott (HP) Filter, a time-series decomposition tool that separates a time-series into trend and cyclical components; and linear splines, where the data is modelled as piecewise linear segments over fixed 25-year intervals, allowing for abrupt changes in slope at these specified breakpoints. As the other methods tend to have strong agreement, for the local level, only CEEDMAN and LOWESS were used.