On the Non-Linear Nature of Long-Term Sea Surface Temperature Global Trends

Abstract

The year-to-year rate and acceleration of the changes in global sea surface temperatures (SSTs) for 1870–2022 are determined using non-linear techniques. Our methodology identifies the non-linear, long-term tendencies of the warming, revealing subtle but essential features of the SST changes that are impossible to identify by linear techniques. Our analysis identifies inhomogeneous patterns of SST evolution. For instance, in most equatorial regions, including the western Pacific, the acceleration of the warming is positive, and sustained warming replaces the cooling at the beginning of the records. On the contrary, the warming acceleration in the central and eastern equatorial Pacific is negative, indicating a sustained warming at the beginning of the record that weakens over the years, resulting in a stagnated warming or even a cooling in recent years. We identify a strong warming with a negative acceleration for the North Atlantic. We hypothesize that this warming is part of the multidecadal variability of SST rather than a long-term trend.

Introduction

Human influence on the climate system is a fact (IPCC, 2021). According to the Annual 2022 Global Climate Report of the National Aeronautics and Space Administration and the National Oceanic and Atmospheric Administration (NOAA), the globally averaged temperature has increased by about 1.15°C since the late nineteenth century (NOAA, 2023), warming related to increasing greenhouse gas emissions resulting from human activities (IPCC, 2021). Most of the warming occurred in the past 40 years, with 2015–2022 being the warmest period on record. One critical variable that permits the long-term characterization of climate evolution is the sea surface temperature (SST), which results from radiative and turbulent processes at or near the atmosphere-ocean interface (Deser, Alexander, et al., 2010).

Global or hemispheric analysis of SST trends supports a large body of available literature regarding climate change (IPCC, 2021). In the last four decades, available satellite data have provided detailed and accurate measurements of ocean temperatures, showing that there has been an increase in global SST close to 0.18°C per decade since 1981 (about 0.7°C until 2022), which has not been uniform in space or time (NOAA, 2023). Consequently, due to the non-uniform warming in ocean or land temperatures (Ji et al., 2014), critical questions remain regarding the physical interpretation of the long-term behavior of the climate system determined only in global or hemispheric analysis. In this regard, estimates of long-term SST trend patterns represent a valuable asset in climate research. Despite limitations inherent to poor spatial and temporal sampling, SST data is currently our best source of information available to test hypotheses and gain insight into the physical processes responsible for the local warming and cooling of distinct ocean basins.

Most long-term SST trend patterns rely on linear techniques, which help extract information from measurements satisfying two key assumptions: linearity and stationarity, implying that the SST can change only at a constant rate (Ji et al., 2014). The global warming resulting from economic activities, however, impacts the state of the World Ocean and the atmosphere, introducing changes in the climate that may result in oscillatory modes of variability of different frequencies, which, in turn, may induce non-stationary and non-linear pattern evolutions (Qin et al., 2012). Therefore, a simple straight line fitted to the data may be appropriate to analyze a linear and stationary world, but not in the natural and complex climate analysis, which focuses, for example, on the comprehension of global warming (Wu et al., 2007).

Because the spatially inhomogeneous warming is not time uniform, we need to diagnose the warming evolution on different spatial-temporal scales (Ji et al., 2014). Constructing a trend has many forms, but we can generalize it using linear and non-linear categories. Commonly, there is a priori determination of the shape of the trend (e.g., a time-constant linear trend or a time-varying exponential trend) (Ji et al., 2014). However, as noted by Wu et al. (2007), “a rigorous and satisfactory definition of either the trend of nonlinear nonstationary data or the corresponding detrending operation still is lacking, which leads to the awkward reality that determination of trend and detrending often are ad hoc operations.” These authors introduced a simple trend definition for non-linear, non-stationary time series, establishing their definition on an intrinsically determined monotonic function or a function that can have at most one extremum in the period considered. This trend definition was used by Ji et al. (2014), who applied the Ensemble Empirical Mode Decomposition (EEMD, developed by Wu & Huang, 2009) to global land surface air temperature. They found that the evolution of the zonally averaged trend of surface air temperature started to increase (>0.5°C) in 1920 in the northern sub-polar and sub-tropical regions, with the amplitude of the warming growing slowly and extending toward mid-latitude during the next three decades and, from 1955 onwards, an accelerated warming was present in all latitudes. Thus, Ji et al. (2014) showed that the EEMD technique allows the extraction of the evolution of coherent spatial structures at different timescales.

The central methodology used in this study is the Complementary Ensemble Empirical Mode Decomposition (CEEMD), introduced by Yeh et al. (2010), which improves EEMD results by eliminating the residue of added Gaussian white noise. Here, we use the CEEMD following Ji et al. (2014), who, as mentioned, applied the EEMD to spatial-temporal data such as gridded climate data, showing the advantages of this methodology to diagnosing climate system evolution. We use the CEEMD to extract non-linear, long-term trends from two available global SST data sets at each grid point. The technique adopted here is aptly fit to treat non-linear and non-stationary SST long-term records (e.g., Martinez-Lopez et al., 2018). Specifically, the nonlinearity of the extracted trends enables us to estimate their first and second derivatives, bringing additional information about the long-term evolution of the SST fields, particularly their accelerated rates of change.

This work examines the spatial and temporal structure of the non-linear, long-term evolution of SSTs over the World Ocean during the historical period, determining their year-to-year rates of change and acceleration. The last piece of information, from the non-linear technique, puts the evolution of dynamic processes in a new global perspective, linking different and remote ocean areas. This information is essential to contrast with those obtained from numerical models used to construct climate change scenarios for the near future.

Data and Methods SST Data

We use yearly averaged SST values in the period 1870–2022. Yearly values were constructed by averaging the monthly values provided by two historical reconstructions available from the NOAA Extended Reconstructed SST (ERSST), Version 5 (Huang et al., 2017) and from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST; Rayner et al., 2003).

In order to test some aspects of data uncertainty in the ERSST and HadISST data sets, we use the last four decades (1982–2022) data available from NOAA OI SST V2 High-Resolution Data set provided by the NOAA PSL, Boulder, Colorado, USA (https://psl.noaa.gov/data/gridded/) (Huang et al., 2021). Over these decades, the uncertainty of SST data has decreased significantly due to newly available measurements and technology from buoys, ships, and satellites (Huang et al., 2021; Kennedy, 2014).

Data Analysis

We separate, at each grid-data point, the yearly averaged SST time series into their different frequency and amplitude components using the CEEMD method, in which the residual, once we extracted all oscillatory components from the analyzed series, is either a monotonic function of time or has only an extremum. We repeat the CEEMD method 200 times in each grid-data point (changing in each grid-data point and iteration the settings of the used random number generator), thus obtaining an ensemble of CEEMDs and, from it, an ensemble average, which we referred to simply as CEEMD. The residual of this CEEMD is what we identified here as the non-linear, long-term SST trend in units of degrees Celsius. Figure 1 illustrates the methodology used in two locations over the equatorial Pacific to extract oscillatory components and the trend from the yearly SST data. The two upper rows correspond to a location in the west equatorial Pacific. Figures 1a–1c shows the six different frequency and amplitude components (CEEM₁ to CEEM₆) and one non-oscillatory component, which is either a monotonic function of time or has only an extremum, and it represents the non-linear, long-term SST trend (R). Figure 1d shows the reconstructed SST, obtained by adding the six oscillatory components plus R (black line), an estimate of the multidecadal variability (MDV) (resulting from adding CEEM_4–6 plus R), and R. By differentiating R and the MDV, we estimate the warming rate in two different time scales (Figure 1e). By differentiating again, we obtain the acceleration of the warming rate (Figure 1f). The two lower rows of Figure 1 correspond to a location in the eastern equatorial Pacific (Figures 1g–l).

Figure 1. Example of the Complementary Ensemble Empirical Mode Decomposition (CEEMD). The two upper rows show the CEEMD applied in the western equatorial Pacific. The two lower rows show the corresponding information in the eastern equatorial Pacific. See the main text for details.

Here, we mainly focus on the long-term trend analysis; therefore, it is essential to note that the long-term SST trend in the western equatorial Pacific is upward concave, whereas in the eastern equatorial Pacific, it is downward concave. Consequently, the second derivatives are positive and negative, respectively. Note that we can easily estimate the averaged acceleration in the longest time scale due to our long-term trend definition. In our case, we generally use the whole period (1870–2022), although the analysis of NOAA OI SST spans from 1982 to 2022. In this shorter-period analysis, however, the characteristics of the trend definition remain valid. If we consider the MDV, its acceleration shows complex oscillatory behavior (Figures 1f–1l), making the task difficult and thus requiring additional analysis. In the discussion section, we show an example of such analysis.

We estimate the statistical significance of the extracted non-linear trends at each grid-data point of the SST data sets using the methodology proposed by Ji et al. (2014), which uses a Monte Carlo method to determine the statistical significance of the secular trends; we give a detailed description of this methodology in Text S1 and Figures S1–S3 in Supporting Information S1. We also estimate linear SST trends using ordinary least squares to compare non-linear and linear results. We tested their statistical significance using a two-tailed Student's test with N-2 degrees of freedom, considering serial correlation following Santer et al. (2008). This methodology reduces the degrees of freedom when the regression residuals present serial correlation (Qian, 2016). Finally, we obtain estimates of the long-term rate of change of SST and its acceleration by applying a simple differentiation scheme to the SST trends once and twice, respectively.

Notice that we obtain respectively approximately straight sloped lines (Figures 1e–1k) and slightly changing values resembling a flat line (Figures 1f–1l) for the first and second derivatives, which is due to the smooth character of the residuals (long-term trends) resulting primarily from the trend definition and secondary from the averaging procedure applied to the CEEMDs. Therefore, the distribution of the second derivative values at each grid-data point is approximately normal (200 times 151 values), implying that the averaged-time second derivative is a reasonable estimate of the expected second derivative value. We confirm that using the statistical procedure of Bootstrapping, resampling the second derivative series of SST at each grid-data point to construct confidence intervals for the warming acceleration. The time-averaged second derivative was within the confidence limits in all grid data points, which means that the slightly changing acceleration values could also be a constant for the sake of simplicity. In the case of MDV, the oscillatory characteristics at this time scale (see Figure 1) imply that the corresponding first and second derivatives are also oscillatory. Hence, the time-averaged second derivative is not a reasonable estimate for MDV, and additional analyses are required, which, however, is beyond the scope of the present work.

Results Analysis of Yearly SST Data at Some Locations Around the World

Before showing the SST global trend patterns, it is instructive to look at the trend behavior of selected points close to the equator where HadISST and ERSST trends are comparable in magnitude, although with contrasting pattern evolution.

Figure 2 shows the time series of yearly SST obtained from HadISST, their corresponding non-linear trends, and the first and second derivatives at three locations over the equator: the Indian Ocean, IO (59.5°E, 0.5°S), the western Pacific Ocean, WP (146.5°E, 0.5°S), and the eastern Pacific Ocean, EP (100.5°W, 0.5°S). Close inspection of panels 1a, 1b reveals that in the first decades, the SST trend shows a slight cooling followed by sustained, long-term warming in both IO and WP locations (thick magenta and red curves). In contrast, the warming is less pronounced in the EP (thick blue curve), exhibiting a different behavior compared to WP and IO (Figure 2a). Notably, after the moderate warming, a cooling trend develops at the EP location during the last three decades of the observing period.

Figure 2. Time-series of yearly SST and their corresponding secular trends at three locations over the equator (a), as well as the yearly SST, secular trends with their initial value subtracted (b), their first derivative (c), and their second derivative (d). Units for yearly SST series and their secular trends are in Celsius degrees, Celsius degree per year for the first derivative, and Celsius degree per year2 for the second derivative.

Figure 2b shows all three non-linear trend curves with their initial values subtracted. At the IO location, there is a slight cooling trend from 1870 to 1900, initiating a sustained long-term SST increase, with a continually increasing warming rate, since 1900 (magenta curve). A long-term SST increase of about 1.23°C at this location from 1900 to 2022 represents the maximum increase over the equator, corresponding to a SST value of about 29°C in 2022. At the WP location, the SST shows a value of about 30°C in 2022, representing the maximum value across the equator. This location also shows a slight initial cooling from 1870 to 1916, followed by a sustained long-term SST increase but with a lower increasing warming rate than the IO location. The long-term SST increase since 1916 is about 0.86°C in 2022 (red curve). In contrast, at the EP, the SST increases faster in the first decades, then the warming slows down, reaching a maximum value of about 0.80°C by 1980, after which it finally starts to decrease (blue curve). Close inspection of the long-term SST in both WP and EP reveals an approximately symmetric behavior compared to a temperature average of the two locations (we will report details and further analysis on this issue elsewhere).

We obtain an estimate of the warming rate by using the first derivative of the non-linear SST-trends (Figure 2c). IO exhibits the fastest warming rate, with values progressively increasing from about −0.005°C/year in 1871 to 0.02°C/year in 2022. The negative values reflect the initial cooling observed at the IO location during the first years. In the WP, a similar behavior to that of the IO is evident, although with a slower warming rate. A faster-decreasing warming rate is present at the EP location, which has turned the warming from 1870 to 1980 into a cooling trend relative to the maximum SST increase of 1980 during the last 42 years (1981–2022). We note that in 1946, the WP and the EP warming rates were almost identical (circa 0.0045°C/year) and very close to their nearly similar time-averaged value (about 0.0046°C/year). Previous to 1946, the warming rate in the EP was greater than that of the WP. From 1947 onwards, the warming rates at the EP location were lower than WP's.

We can also detect subtle changes in SST trends by taking the second derivative of the time series at these three locations (Figure 2d). First, the IO and WP sites exhibit positive, almost time-invariant values, suggesting an upward concave behavior in the SST trend. In the IO, the second derivative values slightly decrease in time. At the same time, they increase in the WP, suggesting a decelerating upward trend in the IO and an accelerating upward trend in the WP. On the other hand, its eastern counterpart always remains negative, corresponding to a downward concave behavior. These contrasting trends point to different dynamical and thermodynamical processes occurring near or at the ocean's surface. Thus, regarding trend detection, second derivatives obtained using this non-linear methodology can guide the interpretation of the long-term behavior of the ocean and atmospheric interaction.

Next, we consider a second pair of points quite distant from each other but with a similar long-term SST trend behavior. We give such an example in Figure 3, which shows two locations: the northern Gulf of California (NGC, magenta lines) and the East China Sea (ECS, red lines). Both locations share a slight cooling trend during the first decades of the observing period (1870–1895 ECS; 1870–1916 NGC). However, after that, they also share a marked warming trend (Figures 3a and 3b) A 20-year delay is also clearly present in the onset of the ensuing warming period for these two locations. During the warming period, we note that ECS warms at a faster rate compared to NGC until about the early 1970s, when both sites warm at the rate of about 0.02°C/year (Figure 3c), which is a rather significant warming rate value compared to the rest of ocean. Then NGC warms faster than ECS, reaching warming rates of about 0.041°C/year and 0.036°C/year by 2022, respectively. As far as concavity is concerned, both sites exhibit almost constant positive second derivatives. However, the larger values (by about 37%) occur at the Gulf of California site (Figure 3d). We summarize this observed SST trend behavior as having two periods, the first one being a decelerated cooling in the first decades while the second period being an accelerated warming; however, these two periods share almost the same trend acceleration.

Figure 3. As Figure 2, but for two locations exhibiting similar, sustained warming. One is located in the Gulf of California (112.5°W, 29.5°N), and another one is located in the East China Sea (121.5°E, 28.5°N).

Finally, we consider some points in regions that have sustained cooling for the last 153 years. We give two examples in Figure 4. First, in the northern Atlantic (Figures 4a and 4b, blue curves), at the southern tip of Greenland, continuous cooling is evident, which is a well-known feature emerging from different linear analyses of available temperature data and climate simulations (Deser, Phillips, & Alexander, 2010; Rahmstorf et al., 2015). Less known, however, is the slight cooling present in the coastal waters of the northern Gulf of Mexico (Figures 4a and 4b, red lines). Notice that in both locations, the cooling rate has been increasing (Figure 4c). However, its rate of change (acceleration), which is more prominent in the North Atlantic than in the northern Gulf of Mexico, has slightly diminished over the analyzed period (Figure 4d).

Figure 4. As Figure 2, but for two locations exhibiting sustained cooling. One is located in the coastal waters of the northern Gulf of Mexico (89.5°W, 30.5°N), and another one is located in the North Atlantic (40.5°W, 52.5°N).

SST-Warming Patterns

Results from the last section prompt us to introduce the Warming Rate Patterns (WRP) concept built upon the linear and non-linear analysis of yearly SST series, showing how the surface waters have been heating (or cooling) in time average since 1870. In addition, we introduce the concept of Acceleration Warming Rate Patterns (AWRP), which describes how the WRP is evolving in a time-averaged fashion.

First, we obtain a linear estimate of the WRP for both data sets (Figures 5a and 5b). Some differences are apparent between HadISST and ERSST. Along the equator, surface waters warm faster in ERSST than in HadISST, particularly in the central Pacific; along the continental margins of the north Pacific Ocean, SST increases exceeding 1°C are present in HadISST, while in ERSST, they exhibit lower values. Then, we obtain the corresponding non-linear estimate of the WRP using the time average of the first derivative of the non-linear long-term SST trends (Figures 5c and 5d). By taking the difference, at each grid-data point, between the final and the initial values of the extracted non-linear SST trends, the results are similar to the non-linear WRP (not shown).

Figure 5. Warming Rate Patterns (WRP) derived from linear and non-linear trends of yearly SST data from HadISST (a, c) and Extended Reconstructed SST (ERSST) (b, d). The linear patterns show simply the warming rate multiplied by the analysis period (153 years). In contrast, we obtain the non-linear ones by multiplying the time average of the first derivative of the non-linear trends of yearly SST data by the analysis period. Positive WRP values indicate warming (long-term increase of SST over the analysis period), while negative values indicate cooling (long-term decrease of SST). Units in Celsius degree per year, multiplied by 153 years. Areas with overlaid hatching show statistically significant trends. For Acceleration Warming Rate Patterns, we use the time average of the second derivative of the non-linear trends of yearly SST data from HadISST (e) and ERSST (f). Positive values indicate that the SST-warming trend is concave up, implying sustained warming of the surface ocean waters (the time-averaged warming rate is increasing with time). Negative values indicate that the SST-warming trend is concave down, implying a stagnated SST-warming of the surface ocean waters or even an SST-cooling relative to its reached maximum (the time-averaged warming rate decreases with time).

The non-linear WRP from HadISST and ERSST share similar features and apparent regional differences (Figures 5c and 5d). For example, the vast cooling region in the northern North Atlantic is present in both data sets. A significant SST-warming is also present, which covers the northeastern North American continental shelf (from Cape Hatteras to Nova Scotia) and surrounding ocean waters. Notice that the most considerable cooling is observed in ERSST, while the most extensive warming occurs in HadISST. However, the spatial pattern of this warming/cooling dipole is consistent in both data sets. In contrast, in ERSST, a considerable SST-warming is found over the Arctic, with SST increases as high as 3°C in 153 years, while in HadISST, the Arctic increases do not exceed 2°C in the same period. Near the equator, in ERSST, the maximum values are located over the southward flowing warm Angola Current, with SST increases exceeding 4°C. In contrast, in HadISST, much lower values, around 1.7°C, are found in this region.

Some features are only present in HadISST. For example, a warming/cooling small dipole structure in the ECS and a cooling region cover part of the northern Pacific Ocean. Finally, notice that a large cooling area around the Antarctic is present in both data sets. The poor spatial and temporal sampling at this region increases uncertainty, hindering at the moment, its physical interpretation.

Figures 1–4, 5c and 5d illustrate the dynamical character of the long-term, non-linear SST-warming: the surface ocean waters are warming at different rates, generating spatial gradients that evolve over time. The time average of the first derivative at each grid-data point generates a simplified, two-dimensional WRP, showing the time-averaged SST-warming rates (Figures 5c and 5d). Visual inspection of both linear and non-linear WRP in Figures 5a–5d reveals similarity in both amplitude and patterns. There are, however, subtle and significant differences that are impossible to identify by limiting ourselves to linear analysis alone. For example, a visual comparison of non-linear WRPs (Figures 5c and 5d) indicates that the central and eastern equatorial Pacific is warming, with lightly higher warming rates in ERSST, resulting, in general, in SST increases not exceeding 1.2°C in the analyzed period. We can use the information of the second derivative, instead of the first one, to construct another pattern, the AWRP, which shows, in a simplified way, part of the long-term non-linear dynamical behavior of the surface ocean waters (Figures 5e and 5f). In contrast to the other analyzed patterns, in the AWRP, there are evident differences between ERSST and HadISST. In HadISST, a consistent negative pattern emerges in the central and eastern equatorial Pacific, indicating that a large part of the surface waters of the south and central Pacific do not exhibit sustained warming but rather a warming that has been declining over the years, leading to a stagnated SST-warming, or even to a recent cooling in some regions. In ERSST, a smaller negative pattern covers some central and tropical South Pacific regions; however, positive values in the eastern equatorial Pacific indicate sustained SST warming. Finally, regarding the cooling and warming trends at the North Atlantic and the eastern US coast, both ERSST and HadISST show an increasing accelerated warming at the eastern US coast and an almost constant accelerated cooling at the northern North Atlantic.

Satellite-Based Warming Patterns

This section shows WRP and AWRP extracted from NOAA OI SST V2 High-Resolution Data set for 1982–2022. To directly compare, we also show the corresponding patterns extracted from HadISST and ERSST using the same analysis period. We expect NOAA OI SST patterns (OISST) to share many common characteristics with both HadISST and ERSST since they all use satellite data in the analysis period.

Both linear and non-linear WRP from OISST (Figures 6a and 6b) show a dominant cooling structure in the Pacific Ocean, initiating close to the dateline over the equator and increasing toward the southern coast of South America, where maximum cooling values exceeds 1°C. The OISST WRP resembles the spatial structure of the HadISST negative AWRP observed in both the eastern equatorial and the southern, eastern Pacific Ocean (Figure 5e). The WRP of both ERSST and HadISST give us comparable information (Figures 6c and 6d), although the cooling is more pronounced in ERSST. In addition to the cooling pattern of OISST in the eastern Pacific, there are two significant warming areas in the tropical and extra-tropical South and North Pacific Oceans, as well as one in the North Atlantic. Comparable warming patterns are also present in HadISST and ERSST, although with some minor regional discrepancies. Finally, a significant cooling area around the Antarctic is present in all data sets, suggesting that this cooling is a robust characteristic during the period 1982–2022.

Figure 6. Warming Rate Patterns (WRP) derived using linear (a) and nonlinear trends (b) of yearly SST data from OISST, and nonlinear-trends of yearly SST data from Extended Reconstructed SST (c), and HadISST (d). The linear patterns show simply the warming rate multiplied by the analysis period (41 years), while the non-linear ones are obtained multiplying the time-average of the first derivative of the nonlinear-trends of yearly SST data by the analysis period. Areas with overlaid hatching show statistically significant changes. Positive WRP values indicate warming (multi-decadal increase of SST over the analysis period), while negative values indicate cooling (multi-decadal decrease of SST). Units in Celsius degree per year, multiplied by 41 years.

The AWRP from OISST (Figure 7a) shows a negative pattern in the south-eastern tropical Pacific, which is also present in both HadISST and ERSST (Figures 7b and 7c), although in these data sets, it is not so clearly defined as in OISST. However, the most dominant characteristic of the OISST AWRP is the negative pattern in the North Atlantic, covering part of the southern Caribbean Sea, the northern coast of South America and from there until Africa, and from there to the entire sub-polar gyre in the North Atlantic. This pattern is highly consistent with those in HadISST and ERSST. In the northeastern Pacific and western coast of the USA, OISST and ERSST agree, showing a positive AWRP, which is not so evident in HadISST.

Figure 7. Acceleration Warming Rate Patterns from OISST (a), HadISST (b), and Extended Reconstructed SST (c). Positive values indicate that the SST-warming trend is concave up, implying sustained warming of the surface ocean waters (the time-averaged warming rate is increasing with time). Negative values indicate that the SST-warming trend is concave down, implying a stagnated SST-warming of the surface ocean waters or even an SST-cooling relative to its reached maximum (the time-averaged warming rate decreases with time). (e) Evolution of yearly SST data (thick black line) at a location of the subpolar gyre where the cooling is significant (41.5°W, 53.5°N). Its long-term trend (thick dashed black line) and multi-decadal variability (thick black line). The yearly SST series from OISST (blue line) is consistent with that from HadISST during the last four decades, which also applies to their corresponding trends (thick blue and red lines, respectively).

Notice that time-averaged positive AWRP is associated with regions that could be experiencing first a cooling and then a warming or an intermediate state (see Text 2 and Figure S4 in Supporting Information S1). For example, a significant extension of the Gulf of Alaska shows a positive, well-defined AWRP pattern and an intense, positive WRP pattern, which together indicate that, on average, the Gulf of Alaska has warmed during the last four decades. A plot of the SST and its long-term trend reveals a slight cooling from 1982 to 1990 and then a sustained warming, leading to an increase of the SST of about 2°C from 1982 to 2022 (Figure 7d, black line). On the contrary, the Labrador Sea shows a negative, well-defined AWRP pattern, but it also shows, as in the Gulf of Alaska, an intense, positive WRP pattern. In this case, the time-averaged negative AWRP indicates that this region is experiencing warming and then cooling or an intermediate state (see Text 2 and Figure S4 in Supporting Information S1). Figure 7d also reveals the SST evolution in the Labrador Sea. First, intense, sustained warming occurs, and then stagnation dominates over the last 5 years.

This behavior in the Gulf of Alaska and the Labrador Sea is similar to that in the western and the eastern equatorial Pacific, although on a shorter time scale. In both cases, it is suggestive of an oscillation. What kind of forcing mechanism could link the multidecadal oscillation in the northeastern Pacific and the Labrador Sea must be elucidated, and it is beyond the scope of the present work. In the second case, the oscillation might result from the dynamic coupling between the surface waters of the western and eastern equatorial and the evolution of the diffusive thermocline. In case of occurrence, such long-term oscillation can affect the long-term behavior of the zonal contrast of equatorial SST, which plays a crucial role in climate evolution. This topic is beyond the scope of this work and will be part of future investigations.

Discussion

Our methodology allows us to identify long-term, non-linear SST-warming patterns: the surface ocean warms at different rates in different parts of the world, slowly evolving in time. The CEEMD permits isolating SST's non-linear, long-term behavior, thereby constructing the second derivative. This additional information allows us to identify subtle and essential differences not apparent by linear analysis alone. A significant result is that many central and eastern equatorial Pacific regions are warming, with higher warming rates in ERSST than in HadISST. Solomon and Newman (2012) calculated long-term SST trends in the tropical Indo-Pacific using least squares linear regressions and reported significant discrepancies among the available SST reconstructions. In particular, in HadISST their linear analysis reveals a cooling area over almost all of the eastern equatorial Pacific, which is not present in ERSST, in agreement with our results. They obtained, however, a more consistent and robust SST-trend among all data sets when they removed, by filtering, El Niño/Southern Oscillation signal from the data. We notice that in our results, the long-term warming over the equatorial Pacific is present in both HadISST and ERSST. Also, the ENSO-related variability is highly consistent in both data sets, which we get using components three to five of the six oscillatory modes obtained using the CEEMD. Our results agree with those of Lee et al. (2012); therefore, filtering ENSO-related variability should have had only minor effects on the long-term SST trends extracted from these data sets.

We construct the AWRP using the second derivative's information, which helps visualize part of the surface ocean waters' long-term, non-linear dynamical behavior. In contrast to other analyzed patterns, the AWRP allows us to identify subtle differences between ERSST and HadISST. For example, in HadISST, a negative AWRP emerges in the central and eastern equatorial Pacific, indicating that many regions of the south and central Pacific do not exhibit sustained warming, like the western equatorial Pacific, but rather a continually decreasing warming, leading to a stagnated SST-warming, or even to a recent cooling in some regions. In ERSST, a smaller negative pattern appears in some regions of the central and tropical south Pacific, in agreement with HadISST. However, positive values dominate the eastern equatorial Pacific, indicating SST warming. We explain these differences by the lower SST values present in ERSST before 1950 than those of HadISST (not shown).

Our results show that the observed North Atlantic warming/cooling dipole is a robust characteristic of the long-term evolution of the SST field, consistent in both data sets. This pattern is known and broadly discussed in the literature. For example, Caesar et al. (2018) used the results of a high-resolution climate model forced with increased atmospheric carbon dioxide concentrations. They indicated that the warm signature of this pattern is associated with the northward shift of the Gulf Stream, relating the cold signature with a reduction in the northward heat transport resulting from a slowdown in the Atlantic Meridional Overturning Circulation. In the last four decades, however, the entire North Atlantic sub-polar gyre has exhibited pronounced warming (Desbruyères et al., 2021; Marsh et al., 2008; Robson et al., 2012). Note that the time-averaged acceleration allows us to identify oceanic areas that vary coherently in the longest time scale inherent to the whole data set or the considered analysis period. In the North Atlantic, our long-term analysis pointed to a cooling trend. However, in the shorter satellite period, it is evident that almost the whole sub-polar gyre warms, although with an average negative acceleration, suggesting the possibility of an initial pronounced warming followed by an eventual cooling after the SST reached its maximum value, which could be part of the MDV. To show this, we chose a grid-data point with a cooling long-term trend to illustrate how the technique helps understand the regional observed warming in a shorter time than the whole period.

Figure 7e shows the evolution of yearly SST at a location of the sub-polar gyre where the cooling is significant (41.5°W, 53.5°N). The HadISST yearly data (black line) exhibits a negative long-term trend (thick dashed black line) superposed with well-defined multi-decadal variability (thick black line). The yearly OISST SST (blue line) is highly consistent with HadISST during the last four decades, which also applies to their corresponding trends (thick blue and red lines, respectively). Note that the HadISST multi-decadal variability (thick black line) shares many common characteristics with the OISST trend (thick blue line) and the HadISST trend for the shorter 1982–2022 period (thick red line), with the most considerable differences at the beginning of the analysis period. Thus, the OISST data set analysis shows that recent warming in the sub-polar gyre prevailed for a couple of decades, reaching its maximum in the first decade of the 21st century and starting to cool since then. Hence, in the long term, there has been cooling in that region, but in the last four decades, there has also been strong warming, which is already diminishing in intensity. This example nicely illustrates the utility of the concept of the second derivative.

Conclusions

The linear and non-linear analyses give us similar information about the rate of change of the long-term SST evolution in a time-averaged fashion. However, only non-linear techniques allow us to estimate second derivatives and identify coherent regions across the world ocean where possible low-frequency oscillations occur. Once an interesting spatial structure is identified, all of its points become amenable to a more detailed analysis that includes its long-term evolution and the time structure of a possible oscillation.

More sophisticated analyses are possible by including algorithms to find maximum, minimum, and inflection points in each long-term SST evolution series by automatically determining the shape of the oscillations in the analysis period. The algorithm would inform when a particular region experienced cooling (warming), followed by warming (cooling); thus, we can determine when the SST reaches its minimum or maximum values. Then, knowing that these extreme values are already reached, it is easy to obtain the per cent reduction (increase) of the actual value concerning that maximum (minimum). Such analysis should be appropriate for some other variables, like precipitation. This technique, therefore, helps analyze complex time series, isolating dynamically essential features that may facilitate their interpretation.

In the present work, we use time-averaged accelerated warming patterns to get compact information about the long-term non-linear evolution of the spatially inhomogeneous warming of the surface ocean waters. We consider these patterns to be suitable tools for analyzing the characteristics of the identified long-term trends, which are determined either by a monotonic function or a function that has at most one extremum in the period considered. Thus, the time-averaging effectively removes the slight fluctuations of the second derivative of the long-term trend, leading to the best estimate for the acceleration. In the case of MDV, its oscillatory characteristics represent a limitation for the time-average procedure. At this time scale, the time-average acceleration for a determined period is a preliminary estimate that needs to be complemented by individual analysis in selected grid-data points in the area of interest using the whole record length. In this regard, we are analyzing a zonal cross-section along the equator to understand better the multi-decadal SST variability observed in such crucial regions for the global climate.

Acknowledgments

We would like to thank the anonymous reviewers for their constructive comments.

Data Availability Statement

The Extended Reconstructed SST (ERSST), Version 5 (Boyin et al., 2017) is available from https://www.ncei.noaa.gov/pub/data/cmb/ersst/v5/netcdf/.

The Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST; Rayner et al., 2003) is available from: https://www.metoffice.gov.uk/hadobs/hadisst/. The NOAA OI SST V2 High-Resolution Dataset (Huang et al., 2021) is provided by the NOAA PSL, Boulder, Colorado, USA (https://psl.noaa.gov/data/gridded/).

We perform the non-linear decomposition using the fast EMD/EEMD/CEEMD Matlab code (Wang et al., 2014) provided by Prof. Yung-Hung Wang. https://sites.google.com/view/yhw-personal-homepage.