The long-term warming of the climate system is indisputable, widely covering the atmosphere, ocean, and cryosphere. In the last decade, the global-mean surface temperature (GMST) was 1.09°C higher than that in 1850-1900, approaching the targets of the Paris Agreement (IPCC, 2021). However, over the past century, the observed GMST has shown a stepwise upward evolution, which should be attributed mainly to the combined effects of the continuously increasing concentration of human-emitted greenhouse gases (GHGs) and the modulations of multi-scale internal variabilities. (Kosaka and Xie, 2013; Kosaka and Xie, 2016; Yao et al., 2016a; Huang et al., 2017b; Yao et al., 2017; Chen and Tung, 2018a).

As an essential role in the climate system, the oceans occupy more than 70% of the Earth’s surface, highly regulating global and regional climate by exchanging momentum, heat, and water vapor with the atmosphere. Energy fluxes at the sea-air interface greatly depend on the sea surface temperature (SST) and several vital atmospheric parameters, including wind speed, air temperature, and humidity (Deser et al., 2010), making SST well suited for monitoring climate change and revealing critical variability patterns. By coupling with the atmosphere, the ocean spontaneously triggers multi-scale variabilities of the GSST (Brown et al., 2015), such as El Niño-Southern Oscillation (ENSO) (Bjerknes, 1972; Rasmusson and Carpenter, 1982; Trenberth and Hoar, 1997; Rayner et al., 2003; McPhaden et al., 2006; Sarachik and Cane, 2010), which can be further classified as eastern Pacific (EP) and central Pacific (CP) El Niño (Trenberth and Stepaniak, 2001; Ashok et al., 2007; Ashok and Yamagata, 2009; Kug et al., 2009; Yeh et al., 2009; Takahashi et al., 2011; Sullivan et al., 2016; Timmermann et al., 2018), and Indian Ocean Dipole (IOD) (Saji et al., 1999; Webster et al., 1999) on the interannual time scale; Pacific Decadal Oscillation (PDO) (Mantua et al., 1997; Minobe, 1999; Newman et al., 2016) and Atlantic Multidecadal Oscillation (AMO) (Folland et al., 1986; Schlesinger and Ramankutty, 1994; Delworth and Mann, 2000; Schlesinger et al., 2000; Knight et al., 2005) on the decadal to multi-decadal time scales (Yao et al., 2016a; Chen and Tung, 2018a). These internal variabilities at different time scales, superimposed on the long-term trend, vividly compose the diverse climate change processes.

Great efforts have been made to reveal the multi-scale features of GSST. There is, however, no indication that long-term trends should be linear, despite the fact that some earlier research on GSST trends were predominantly based on linear assumptions (Wu et al., 2007; IPCC, 2013). And the results related to the linearity approach are also highly dependent on the subjective decision of the time interval (Ji et al., 2014; Xu et al., 2021; Xu et al., 2022). Likewise, there is no prior function of any oscillatory components on the other time scales. The period and frequency should be time-varying with the non-stationary nature of time series (Huang et al., 1998). Yet the studies based on the Fourier Transform or low-pass-filter are also heavily limited by arbitrary assumptions, entirely departing from the core principle of ‘data-driven’. Moreover, some studies focusing on multi-scale variabilities of surface temperature are discussed from a global or regional mean perspective (e.g., Yao et al., 2016a; Zhang et al., 2019; Coleman, 2022). However, analysis based on the GMST or global-mean sea surface temperature (GMSST) may eliminate a number of local characteristics and then arrive at a general conclusion that the secular trend has contributed to the majority of the deviations. At any grid points globally, does the long-term trend predominate the variance contribution? Although Mann and Park (1994) developed a method named Multi-taper Method Singular Value Decomposition (MTM-SVD), which could carry out a local frequency domain spatio-temporal decomposition of data variance, evaluating if a particular large-scale pattern exists within a narrow frequency band and reconstructing. They argued that there is no consistent evidence to prove the existence of decadal or longer internal variabilities (Mann et al., 2020), and volcanic forcing may be a potential driver (Mann et al., 2021). However, it is hard to fully demonstrate the local features according to the specific frequency determined from the entire variable field, and motivated by previous studies, we aim to uncover multi-scale variabilities’ spatio-temporal features using observational data instead of arguing whether they are part of the climate system or a response to external forcings.

Additionally, analyses using the Empirical Orthogonal Function (EOF) method are constantly plagued by “mode mixing,” which is partly caused by the blending of the signatures of global warming and other dynamical patterns. Alternatively, several dynamical patterns linked to their dominant frequency bands appear in the same spatial mode, which greatly confuses the real physical facts (Zhang et al., 2010; Chen et al., 2017; Chen and Tung, 2018a; Feng et al., 2020a; Feng and Tung, 2020b). Besides, though Huang et al. (2017b) proposed the decadal modulated oscillation (DMO), covering the ENSO, PDO, AMO, and Arctic Oscillation (AO), they do not separate these climate variabilities well in their Figure 7 and Table 2 from a frequency perspective. Thus, perhaps it is more appropriate to separate the frequency bands of the GSST series grid-by-grid before using EOF analysis. Then apply the rotation algorithm to eliminate the confounded signatures, which may be more consistent with the actual physical facts.

To better investigate the multi-scale variability features of the GSST over the past century, here we take advantage of an adaptive analysis method named Ensemble Empirical Mode Decomposition (EEMD) (Huang et al., 1998; Huang and Wu, 2008; Wu and Huang, 2009). By now, this advanced method has been widely taken advantage of in climate and environment research, driving several novel conclusions (e.g., Franzke, 2009; Qian et al., 2010; Wu et al., 2011; Chen et al., 2014; Ji et al., 2014; Niu et al., 2017; Xu et al., 2021; Wu et al., 2022; Xu et al., 2022). Benefiting from this method, the high-frequency components (HFC), low-frequency components (LFC), and a secular trend (ST) of the SST series can be extracted (See Section 2 for more details), the variance contributions of which will be further investigated. Subsequently, EOF analysis is applied to separate the principal spatio-temporal structures in each frequency band. Furthermore, the Pairwise-Rotated EOF (REOF) analysis (Chen and Wallace, 2016; Chen et al., 2017; Chen and Tung, 2018a) is then employed to eliminate ‘mode mixing,’ making it approachable to some well-known climate patterns. Overall, a broader understanding and interpretation of these multi-scale variability features of the GSST could enrich and deepen our comprehension of the global surface climate change even the climate system.

The other sections are arranged as follows: Section 2 describes the SST dataset and research methodologies. Section 3 discusses the main features of the GMSST. Section 4 displays the local variance contributions of the multi-scale variabilities. Section 5 demonstrates the main space-time modes of the GSST tied to different frequency bands. Finally, Section 6 summarizes and discusses this paper.

Since the industrial revolution, a large amount of carbon dioxide emitted by human activities has been driving up the Earth’s surface temperature due to its greenhouse effect. As a key indicator for monitoring climate change, the GMSST shows a long-term warming trend over the past century, while overlays several oscillations on relatively higher frequency bands. Although GMSST is of smaller amplitude than GMST owing to the smaller heap capacity of the land, both of them have similar characteristics. As presented in Figure 1A , the evolution of GMSST resembles a rising staircase (Kosaka and Xie, 2013; Kosaka and Xie, 2016; Yao et al., 2017; Chen and Tung, 2018b; Huang et al., 2022): the secular upward trend was interrupted occasionally by several hiatus periods (the blue shading), defined by where the instantaneous warming rate of multi-decadal variabilities (MDV) is zero, occurring in the early twentieth century (1900-1908), the mid-twentieth-century (1942-1974), and the junction of the twentieth and twenty-first centuries (2003-2010), respectively. Overall, the global warming slowdown period we obtain highly matches the previous studies (e.g. Kosaka and Xie, 2016; Yao et al., 2016a; Yao et al., 2017). A common feature is that global warming hiatus stages are often accompanied by the opposite orientations between LFC and ST, partially canceling each other out. While during rapid warming, LFC remains relatively stable or slightly increases, thus highlighting the upward trend of the ST. Figure 2A shows the evolution of the multi-scale variabilities of mean SST for several ocean basins, with their correlation coefficients (Figures 2B–D). Their long-term trends exhibit extremely similar consistency (see Figure 2D ). Moreover, it can be observed that global warming hiatus tends to correspond with declining stages of LFC in the major four oceans. However, no similar features are present during rapid warming. Our conclusion does not contradict the results of Yao et al. (2017), because the long-term trend is not separated in their work. The spatial patterns obtained by regressing the GSST onto the HFC, LFC, and ST are shown in Figures 1B–D , respectively. Regarding the HFC ( Figure 1B ), the northern tropical Atlantic mode (NTAM) and the southern subtropical Atlantic mode (SSAM) in the Atlantic, as well as the Indian Ocean Basin Mode (IOBM) of the Indian Ocean basin, are all mentioned as being associated with ENSO variability at low latitudes (Rasmusson and Carpenter, 1982; Lau et al., 1997; Xie and Carton, 2004; Huang and Shukla, 2005). At mid-latitudes, the PDO-like and Southern Hemisphere PDO (SPDO)-like patterns are also linked to the tropical ENSO variability through the atmospheric bridge. While for the LFC ( Figure 1C ), two dominant spatial patterns emerge, characterized by the horseshoe-shaped pattern in North Pacific and an interhemispheric dipole in the Atlantic, that is, uniform warming in the North Hemisphere while cooling in the South Hemisphere, and vice versa, which bear resemblance to the PDO and AMO, respectively. The results shown above are roughly similar to Wu et al. (2011) and Chen et al. (2017). Moreover, since we extend the frequency range from interdecadal to multi-decadal, the oscillations longer than 10 years in the South Pacific and southwest Indian Ocean are also captured (Mantua et al., 1997; Yao et al., 2016b). When turning to the ST ( Figure 1D ), the spatial pattern is almost universally positive, except for the subpolar North Atlantic, which is associated with the so-called “warm hole” (Drijfhout et al., 2012; Woollings et al., 2012; Marshall et al., 2015; Xu et al., 2021; Xu et al., 2022). Meanwhile, the larger values are mainly distributed across the tropical Indian Ocean in the northern hemisphere, the Southern Ocean around the southern side of Africa, the southeastern side of Australia, east of the continents in both hemispheres linked to the western boundary currents with their mid-latitude extensions, and the eastern tropical Atlantic.

Previous studies show that the evolution of the GMSST is typically dominated by the ST (e.g., Chen et al., 2017; Zhang et al., 2019; Xu et al., 2022). However, no more attention was paid to the local features. Here, benefiting from the spatial locality of the MEEMD method (Wu et al., 2009; Ji et al., 2014; Xu et al., 2021), the explained variance of the HFC, LFC, and ST will be analyzed grid-by-grid, respectively, based on the annual-mean SST data. The annual-mean data were used for the EEMD decomposition rather than monthly data since the explained variance of the extratropical regions’ annual cycle component might reach more than 90% (Stine, 2010; Vecchio et al., 2010), which would make it difficult to comprehend the other components.

Figure 3 exhibits the variance contributions (or explained variance) of these three components mentioned above at each grid point. The explained variance of the HFC in the equatorial central-east and south mid-latitude Pacific could reach more than 60%, and that in regions such as the North Pacific, the northern tropical Atlantic, the mid-latitude Atlantic in both hemispheres, and the Indo-western Pacific also exceed 40%. These grid points, which are primarily located in the subpolar North Atlantic and Southern Ocean in both the Pacific and Atlantic sectors, are those where the LFC variance contribution is greater than 40%. Meanwhile, the variance contribution of those grid points over the North Pacific and South Indian Ocean are also around 20-40%. The most significant discrepancies between the HFC and LFC patterns are centralized in the Pacific, with slight deviations in the tropical and subtropical Atlantic (see Figure S1 ). Turning to the ST, it is noted that the key characteristics are somewhat similar to those in Figure 1D , showing that the values greater than 60% are displayed in the tropical Indian Ocean in the northern hemisphere, the Southern Ocean, which extends from the tip of southwest Africa to the southern side of Australia, the Indo-western Pacific, east of the continents in both hemispheres, and the tropical Atlantic, whose areas are also the most noticeable regions where the ST differs markedly from the other two components (shown in Figure S1 ).

The explained variances tied to different frequency bands among different ocean basins have remarkable discrepancies, so here we further investigate the multi-scale variability features of the mean SST in several regions, which is shown in Figures 4 , 5 . From a regional perspective, the secular trends of SST contribute to the largest part of the variance, suggesting that global warming still remains the most conspicuous feature ( Figure 4A ). It is also notable that the variance contributions of the HFC in the Pacific and the Indian Ocean are as many as two and three times higher than the GMSST, respectively ( Figures 4A, B ), which may be related to the interannual variabilities such as the ENSO and IOD. What’s more, it seems that the amplitude of the HFC in the Southern Ocean is smaller than that in the remaining three ocean basins, according to Figure 2A . Another feature is that the HFC of the Southern Ocean has little correlation with other regions (see Figure 2B ), possibly as a result of its isolation, whereas the HFC of the Pacific, Atlantic, and Indian Oceans can interact or connect via some tropical physical processes (Li et al., 2016; Cai et al., 2019). In addition, the explained variances of the LFC in four basins are all approximately two times higher than that of the GMSST, with the largest one in the Atlantic ( Figure 4A ). The reason why the LFC of the GMSST has a smaller explained variance ( Figure 4B ) may be attributed to the superposition of positive and negative phases of the decadal to multi-decadal signals among different ocean basins.

More specifically, Figure 5A presents the variance contribution of all the IMFs. Since the secular trend accounts for most of the explained variance for region-averaged SST, it is not involved in the discussion again here. As mentioned in the methods section, HFC is composed of IMF1 and IMF2. It shows that the HFC variance contribution in the Pacific and Indian Ocean, which are greater than GMSST, could be primarily traced back to the IMF1, with an estimated period of about 3-4 years (see Figure 5B and Table S1 ), roughly three times longer than the sampling frequency of the raw data. There is a well-recognized decadal oscillation signal, PDO, but the explained variances of IMF3 and IMF4 in the Pacific basin are not dramatically greater than the others (see Figure 5A ). The reason may be that the explained variance of the LFC in the Pacific basin is highly diminished due to the offset of compensating patterns of PDO in the northwest and northeast Pacific (Chen and Tung, 2018a). However, the variance contribution of IMF5 is notably higher than the GMSST as well as the other ocean basins, which could be related to the AMO.

In this section, EOF analysis will be used to reveal the main spatio-temporal characteristics of monthly GSST in different frequency bands. Regarding monthly SST series, nine IMFs and a long-term trend can be obtained by the EEMD decomposition, according to Eq. (2). Based on the binary filter property of the EEMD method, the periods of IMF1 are approximately 3 months and IMF2 to IMF9 are twice as long as the former, respectively. In general, IMF1 represents the ultra-high frequency information of the intra-annual scale, while IMF2 and IMF3 contain annual cycle components, all of which are not involved in EOF analysis. Therefore, IMF3 to IMF6 are summed as the HFC, while IMF7 to IMF9 are summed as the LFC.

The three leading EOFs and PCs of the HFC of the monthly GSST are presented in Figure 6 . The explained variance of the first mode (29.4%) is larger than the others. However, due to the orthogonality of the EOF approach, these spatial patterns have overlapping centers of action and perturbation polarities juxtaposed from visual observation, making it difficult to distinguish the signals at tropical and extratropical Pacific. To estimate the optimal angle for rotating EOF1-EOF3 patterns, an evaluation of the SASC metric proposed by Chen et al. (2017) is regarded as a critical reference. As mentioned in Section 2, the other two metrics named FASC and ‘trend transfer’ do not properly suit the discussion here, as the MEEMD decomposition has already separated the components of the GSST series at different time scales for each grid point. The rotation of the EOFs and PCs was performed in two steps. It is clarified that after each rotation step, the (R)EOFs and (R)PCs will be reordered according to the newest explained variance. Figures 7A, B ) exhibit the values of the SASC as a function of the rotation angle, which are associated with the rotation between PC1 and PC3 (step (i)) as well as RPC2 and PC3 (step (ii)), respectively. All the curves in Figure 7 act as sinusoidal-like functions, but at different phase. The minimum evaluated values corresponding to step (i) and (ii) are at about 43° and 78°, with selected regions in (18°N-28°N, 172°W-152°W) and (10°S-10°N, 138°E-178°E), respectively. The modes after rotation are shown in Figure S2 (after the first step) and Figure 8 . After pairwise rotation, it is found that the spatial patterns of REOF1 and REOF2 are predominantly concentrated in the central-eastern and central Pacific along the equator, respectively, representing ENSO diversity, i.e., EP and CP El Niño. The explained variances of these two rotated modes are very close (18.4% and 16.3%) and maintain orthogonality, indicating that besides being the dominant modes at interannual timescale, EP and CP El Niño also exhibit independent properties of each other, which aligns with the previous studies (e.g., Ashok et al., 2007; Ren and Jin, 2011; Timmermann et al., 2018).

Parallelly, Figure 9 exhibits the two leading EOFs and PCs of the LFC, which explaines around 66% variance in total. It is reported that due to the orthogonality constraint of spatial modes, although EOF1 pattern has some AMO features, the signals in Pacific basin are also irrationally assigned, resulting in the mode mixing of the Atlantic and Pacific. Similarly, the SASC criterion is adopted to determine the optimal rotation angle under the target region (0°-66°N, 32°W-2°W), which is plotted in Figure 7C , for detaching the PDO and AMO. It is shown that the optimal rotation angle is 56°. The REOFs shown in Figure 10 are associated with the AMO-like and PDO-like patterns. However, as presented in REOF1, the positive AMO mode is accompanied by a cooling pattern of the equatorial Pacific, and vice versa. By analyzing the low-pass-filtered GSST data, Chen et al. (2017) argued that there exists an Atlantic-Pacific coupling, not an illusion caused by mode mixing, so that REOF cannot eliminate it through any angle in the two-dimensional phase space, but conversely, the PDO-like mode does not demonstrate an AMO-like pattern.

Finally, the main features of the secular trend are displayed in Figure 11 , which substantially matches the findings revealed by Xu et al. (2021, 2022) by portraying the EEMD trends. The two leading EOFs have explained almost all of the variance here, especially the EOF1 even reaching more than 95%. We argue that the EOF1 pattern, as well as the PC1, are mainly related to global warming, with almost everywhere showing the warming trend, except for the subpolar North Atlantic (cooling pattern), the tropical central Pacific, and the Southern Ocean in the Pacific sector (close to zero, implying no signal, according to Equation (3)). Unlike the monotonic increase of PC1, PC2 resembles a parabola with an upward opening. By combining the EOF2, the second mode depicts a secular trend of decrease followed by increasing in the equatorial central Pacific and the Southern Ocean in the Pacific sector (even very slightly), with the turning point around the middle of the last century. Just the opposite, several regions such as the North Pacific show a first upward followed by a downward trend. Although this feature is also represented in Figure 2 given by Xu et al. (2021), the discrepancies still exist and probably can be attributed to the fact that the secular trends are analyzed directly here, whereas they defined the EEMD trend first.

In this study, the EEMD method is applied grid-by-grid to extract the multi-scale variabilities of the GSST. By superimposing the IMFs in nearly frequency bands, the HFC, LFC, and ST at each grid point are obtained, respectively. It is found that the global warming slowdown is highly modulated by the LFC signals stuck in a downward stage, as the decadal to multi-decadal oscillations in the four main ocean basins are mostly in the downward or negative phases. By regressing the GSST onto the HFC, LFC, and ST, we also found that the pattern of the HFC is associated with the ENSO, NTAM, SSAM, and IOBM at low latitudes, as well as the PDO-like and SPDO-like patterns at mid-latitudes, which are induced by tropical Pacific variabilities through the atmospheric bridge. The LFC shows a horseshoe-shaped pattern in North Pacific and an interhemispheric dipole in the Atlantic, related to the PDO and AMO, respectively. With the exception of the subpolar North Atlantic, the ST’s regressing pattern is almost always positive. Larger values are primarily found in the tropical Indian Ocean in the northern hemisphere, the Southern Ocean surrounding southern Africa, the southeast coast of Australia, east of the continents in both hemispheres connected to the western boundary currents with their mid-latitude extensions, and the eastern tropical Atlantic.

By taking advantage of the spatial locality of the MEEMD method, the variance contributions of multi-scale variabilities for each grid point are explored. It has been also revealed that the explained variance of the HFC in the equatorial central-eastern Pacific, and the mid-latitude Pacific in the southern hemisphere could reach more than 60%, whereas the pattern of the LFC is mainly concentrated in subpolar North Atlantic and the Southern Ocean in both Pacific and Atlantic sectors, larger than 40%. The explained variance of the ST hitting beyond 60% is widely throughout the Indian Ocean in the northern hemisphere, the Southern Ocean from the tip of southwest Africa expanded to the southern side of Australia, the Indo-western Pacific, east of the continents in both hemispheres, and the tropical Atlantic. Moreover, the relative variance contributions of the multi-scale variabilities of specific IMFs are further analyzed in several ocean basins, exhibiting characteristics from a regional mean perspective.

The main space-time modes of the variabilities in different frequency bands obtained by EOF analysis show unjustified mixing. Therefore, pairwise rotation is then applied for eliminating the overlapping centers with opposite polarities, which is regarded as a compensator for the regions with the same polarity, enabling the orthogonality of the spatial modes. It is identified that on the high frequency band, ENSO diversity is the dominant, characterized by the presence of CP and EP El Niño in the first two REOF modes, while regarding the low frequency band, the AMO and PDO emerge. It can thus be concluded that the diversity of ENSO, manifesting as EP and CP El Niño, constitutes two primary orthogonal axes on the interannual scale. Meanwhile, when considering the decadal to multi-decadal scales, the AMO and PDO stand out as the key axes, which could be viewed as an extension to the findings of Chen et al. (2017).

There are still some shortages that need to discuss. On the one hand, it has to be admitted that in the early record period as well as during World War II, the scarcity of raw data and the backwardness of observation tools resulted in the poor quality of reconstructed SST data jointly, which perhaps caused a large uncertainty. On the other hand, for those signals with periods longer than decadal, it is so hard to meet statistically significant when we calculate their correlation coefficients with the RPCs due to an extremely few statistical degrees of freedom caused by limited historical record length. Despite these, we still insist that the results gained from the pairwise rotated EOF analysis are enlightening. Besides, deeper mechanisms of these features still need to be further explored, and we also face more blank about how these nonlinear modulated among these ocean basins.

All in all, the purpose of this work is trying to find a novel solution or perspective for the analysis of the multi-scale variabilities, which we are strongly convinced is beneficial to have a better comprehension of the global climate change for both the past and future, even in the slightest hint. It is also worth noting that the multi-scale variability of GSST holds significant implications for atmospheric variables. For instance, overall, under global warming, the long-term increasing trend of the SST may contribute to the intensification of tropical cyclones (TCs), which has been substantiated through observations (e.g., Sobel et al., 2016; Wang et al., 2022), theoretical analyses (Emanuel, 1987), and numerical experiments (Knutson and Tuleya, 2004; Knutson et al., 2020; Emanuel, 2021). However, on a regional scale, internal variability, such as the Interdecadal Pacific Oscillation (IPO), can also exert a significant regulatory impact on TCs’ activity (Zhao et al., 2020). Consequently, how the multi-scale variabilities of GSST influence atmospheric variables and its future evolution, still deserves further research.

Publicly available datasets were analyzed in this study. This data can be found here: NOAA_ERSST_V5 data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their website at https://www.esrl.noaa.gov/psd/.

ZX, GH, and FJ conceived this study. ZX led to writing the initial manuscript in discussion with FJ. ZX did the data analysis and produced figures. FJ improved this manuscript. All authors contributed to the article and approved the submitted version.