Research ArticleCLIMATOLOGY

Aerosol-forced multidecadal variations across all ocean basins in models and observations since 1920

See allHide authors and affiliations

Science Advances  17 Jul 2020:
Vol. 6, no. 29, eabb0425
DOI: 10.1126/sciadv.abb0425


Earth’s climate fluctuates considerably on decadal-multidecadal time scales, often causing large damages to our society and environment. These fluctuations usually result from internal dynamics, and many studies have linked them to internal climate modes in the North Atlantic and Pacific oceans. Here, we show that variations in volcanic and anthropogenic aerosols have caused in-phase, multidecadal SST variations since 1920 across all ocean basins. These forced variations resemble the Atlantic Multidecadal Oscillation (AMO) in time. Unlike the North Atlantic, where indirect and direct aerosol effects on surface solar radiation drive the multidecadal SST variations, over the tropical central and western Pacific atmospheric circulation response to aerosol forcing plays an important role, whereas aerosol-induced radiation change is small. Our new finding implies that AMO-like climate variations in Eurasia, North America, and other regions may be partly caused by the aerosol forcing, rather than being originated from the North Atlantic SST variations as previously thought.


Superimposed on a warming trend, sea surface temperatures (SSTs) in the Atlantic, Pacific, and Indian oceans have undergone distinct decadal to multidecadal variations during the last 100 years (1). The Interdecadal Pacific Oscillation (IPO) and Atlantic Multidecadal Oscillation (AMO) are the leading mode of these SST variations in their respective basins (2) that are thought to result from internal dynamics of our climate system (3), while the Indian Ocean features a basin-wide warming pattern with small interdecadal variations (4). Temperature and precipitation variations over many regions around the globe are correlated with the recent IPO or AMO cycles, and they are often attributed to the IPO- or AMO-associated SST variations in their respective basins through atmospheric teleconnections (58). However, recent studies (911) suggest that volcanic aerosols (VAs) and anthropogenic aerosols (AAs) may have contributed to the AMO cycles since the late 19th century in the North Atlantic. This raises two important questions: (i) Can the aerosols generate AMO-like variations [i.e., in-phase multidecadal SST fluctuations similar to those in North Atlantic SSTs (NASSTs)] outside the North Atlantic? (ii) If yes, then would it be more natural to consider many of the AMO-like climate variations outside the North Atlantic as a local response to the aerosol forcing rather than due to the remote influence from NASSTs, as is commonly thought in the literature?

Here, we attempt to answer the first question and discuss the implications for the second question by analyzing historical SST datasets and model simulations. We aim to quantify the role of external forcing and internal variability in recent decadal to multidecadal SST variations across all three ocean basins. Our results should improve current understanding of the recent decadal to multidecadal variations in the Atlantic, Pacific, and Indian oceans, and their interbasin connections, especially regarding the role of external forcing in generating any in-phase SST variations across all basins. They also have important implications for attributing AMO-like climate variations outside the North Atlantic. Previous studies (1014) have revealed an AMO-like aerosol-forced signal in the North Atlantic only, and other studies (4, 15, 16) also showed a same-sign response to aerosol forcing in the Pacific and other ocean basins. Here, we further show that multidecadal variations in both VAs and AAs may have caused AMO-like (in terms of temporal phasing) multidecadal climate variations across all the oceans and possibly over many land areas as well.


To focus on externally forced variations unrelated to the steady increase in greenhouse gases (GHGs), we low-pass-filtered all the data and removed the GHG-induced changes in the CMIP5 (phase 5 of the Coupled Model Intercomparison Project) all-forcing multimodel ensemble mean (MMM) SST fields through local linear regression against the global-mean SST from the GHG-only CMIP5 simulations (see Materials and Methods). We then applied an empirical orthogonal function (EOF) analysis to the residual SSTs. The GHG-detrended SST fields from 1920 to 2012 are dominated by a multidecadal mode, with negative anomalies from the 1960s to mid-1990s and positive anomalies before the early 1960s and after the late 1990s (Fig. 1A), which are roughly in phase with the observed AMO in the North Atlantic (with a correlation coefficient r = 0.81). It is also anti-correlated with the multidecadal variations in oceanic aerosol loading (r = −0.71 after linear detrending; Fig. 1A). This mode has the same sign of loading, and thus, it is in phase across the three oceans, with the largest amplitude over the North Pacific and North Atlantic (Fig. 1B), as would be expected from the relatively high aerosol loadings over these regions due to AA emissions from East Asia, North America, and Europe (15).

Fig. 1 In-phase, multidecadal SST variations across all ocean basins in observations and models.

(A) Temporal (black line) and (B) spatial patterns of the leading EOF (92% variance) in global (60°S to 60°N) low pass–filtered SST fields from CMIP5 MMM of all-forcing simulations from 1920 to 2012 after the GHG-induced warming (pink line) being removed through linear regression. The smoothed tropospheric aerosol optical depth (AOD) averaged over the oceans (red, gray = trend), linearly detrended annual-mean SST anomalies averaged over the North Atlantic (80°W to 0°W, 0° to 60°N) as in (18) (blue), and annual SST anomalies averaged over the North Atlantic after removing the global-mean SST as in (19) (green) are also shown in (A). (C) Correlation between the black line in (A) and local low pass–filtered SST anomalies from observations during 1920–2012 after removing the GHG-induced warming through regression. The black (gray) stippling indicates statistically significant correlations at the 0.05 (0.1) level based on Student’s t test considering autocorrelation. The data near the right end in (A) [marked by dashed lines and excluded in the calculation for (C)] were derived with mirrored data in the filtering. Results are similar using 13-year moving averages. The outlined boxes in (C) are analyzed further below.

To examine how the forced AMO-like variations are reflected in observed SSTs, which contain large internal variability, similar to the CMIP5 all-forcing case, we removed the component in the observed SSTs associated with the GHG-forced change via linear regression (see Materials and Methods) and then correlated the residual with the forced AMO-like variations (i.e., the black line in Fig. 1A). The GHG-detrended SSTs from observations are strongly correlated positively (r = 0.6 to 0.9) with the forced AMO-like variations in the North Atlantic, tropical western Pacific (TWP), subtropical North Pacific (SNP), and subtropical South Pacific (SSP), whereas they are negatively correlated (r = −0.4 to −0.7) over the tropical central Pacific (TCP) and the subtropical South Atlantic (Fig. 1C). Positive correlations are also seen in the Indian Ocean, although mostly insignificant (Fig. 1C; see Supplementary Text). Similar spatial and temporal patterns of the forced multidecadal variations and their correlations with the observed GHG-detrended SSTs are also found using the CESM1 (Community Earth System Model version 1) large ensemble mean (EM; fig. S1). The positive correlations imply that the externally forced AMO-like signal is evident in observations and in phase with the observed SST variations over those regions.

Regionally averaged SST time series of the GHG-detrended SST from observations and its externally forced and internally generated components (Fig. 2) confirm the correlations shown in Fig. 1C. The forced AMO-like variations are roughly in phase and comparable in magnitude (since the 1950s) with the internally generated (IV) component in the North Atlantic, TWP, SNP, and SSP, leading to strong positive correlations between the forced component and the observed SSTs over these regions. In contrast, SST variations in the TCP and the subtropical South Atlantic are dominated by the IV component [e.g., IPO-like variability; (7, 17)] that is roughly out of phase with the forced AMO-like signal, leading to negative correlations there, although the AMO-like signal does exist in these two regions. In the Indian Ocean (Fig. 2F), the AMO-like signal is fairly strong, but it is not quite in phase with the IV component, leading to relatively weak correlation with the observed variations. Such a phase aliasing between the externally forced and internally generated SST variations in the North Atlantic and other regions complicates the quantification and separation of the internally generated AMO using observations alone. Existing AMO indices (18, 19) do not exclude the externally forced component (red line in Fig. 2A).

Fig. 2 Contributions of internal variability and external forcing.

Regionally averaged time series of the low pass–filtered SST anomalies (°C) from observations after removing the GHG-forced changes via linear regression (at each point) with the global-mean SST time series derived from the CMIP5 GHG-only MMM (Obs*, black line). To estimate internally generated decadal variations due to internal variability (IV, blue line), we removed the changes and variations in observed SSTs via linear regression associated with all the forced signal, which is represented by the global-mean SST time series from the CMIP5 all-forcing historical MMM. The difference between the black and blue lines represents the estimated externally forced decadal variations in observed SSTs (EX*, red line). The last 9 years, plotted using dashed lines, used mirrored data in the filtering and thus are less reliable and were excluded in the correlation calculations for the r shown on this figure. The gray line shows the percentage of the grid boxes within each region with observations in HadSST3. The six regions, outlined in Fig. 1C, are the (A) North Atlantic (ATL; 0° to 55°N, 78° to 2°W), (B) TWP (0° to 25°N, 125° to 160°E), (C) SNP (20° to 32°N, 160° to 192°E), (D) TCP (5°S to 15°N, 180° to 245°E), (E) SSP (40° to 15°S, 185° to 255°E), and (F) tropical Indian Ocean (IO; 30°S to 5°N, 60° to 100°E).

Over the North Atlantic, about 42% of the decadal SST anomalies since 1950 come from the externally forced signal, with the remaining 58% resulting from internal variability (fig. S2), such as that associated with the Atlantic meridional overturning circulation (20, 21). In contrast, over the TWP, the multidecadal anomalies since 1950 result primarily from external forcing, with little (<5%) contribution from internal variability (Fig. 2B and fig. S2). Thus, the IV-induced NASST variations did not cause similar variations in TWP, and the AMO-like variations in TWP are likely a direct response to historical external forcing that caused the AMO-like variations in both TWP and the North Atlantic (and other oceans). To further verify this, we analyzed the TWP SST versus NASST correlations on decadal-multidecadal time scales in a long preindustrial control run by the CESM1 and compared them with those from its all-forcing historical runs. The strong correlations between the TWP SST and NASST seen in the historical runs disappear in the control run that contains only internal variability (Fig. 3). Thus, without external forcing, the AMO-like variations would disappear in TWP, and the IV-induced NASST variations cannot generate AMO-like variations in TWP. We conclude that the AMO-like SST variations since 1920 in TWP in both models and observations are caused by external forcing, rather than originated from NASST variations through atmospheric teleconnection, as suggested previously (8).

Fig. 3 The North Atlantic–TWP SST connection in models.

The histogram of the correlation coefficients between the linearly detrended annual-mean low pass–filtered SST anomalies averaged over the North Atlantic (0°to 55°N, 78° to 2°W) and TWP (0° to 25°N, 125° to 160°E) from a 1800-year CESM1 preindustrial control run (blue) and the 40-member CESM1 all-forcing historical simulations for 1920–2012 (red). For the preindustrial control simulation, 100 overlapped 93-year segments, similar to the length of the observations from 1920 to 2012, are used in estimating the histogram, which have a mean of 0.13. Results are similar if different overlapping is used in selecting the segments. Values above the vertical dashed line are statistically significant at the 0.05 level. The decadal SST variations between the North Atlantic and TWP are correlated significantly for most (37 of 40, >90%, P < 0.05) of the all-forcing simulations, but they are uncorrelated for most of the segments from the control run. This suggests that the correlation results from the external forcing. A cubic spline is fitted to the histograms to produce the smooth curves shown here.

CMIP5 single-forcing simulations show that the forced AMO-like SST variations are mainly caused by VAs and AAs (figs. S2 and S3). For example, VA (AA) has contributed ~50% (~43%) of the simulated NASST decadal variations since the 1950s, while VA has contributed a larger part (>60%) in the tropical Pacific and Indian oceans (fig. S2). A volcanic influence on Pacific and NASSTs has been shown previously (9, 11, 12, 17, 22, 23). Recent volcanic eruptions, such as the 1963 Agung, 1982 El Chichón, and 1991 Pinatubo, caused a cooling response in the pan-tropics during the few years following the eruption. The leading EOF mode of the detrended SST fields from the MMM of the CMIP5 VA forcing–only simulations (Fig. 4A) exhibits variations similar to the AMO-like signal and the observed AMO index (r = 0.65 to 0.87, P < 0.05; Fig. 1A). The results are similar when using only those CMIP5 simulations that extended to 2012 or using CMIP6 VA forcing–only simulations. Thus, we conclude that VA led to AMO-like in-phase multidecadal variations since 1920 across the three oceans, especially in the tropics.

Fig. 4 Impacts of aerosol forcing on global SST.

(A) Time series of the principal component (PC; black line) and (B) spatial pattern of the leading EOF mode, which explains 73% of the variance, in the near-global (60°S to 60°N) linearly detrended, low pass–filtered SST fields from the MMM of CMIP5 VA forcing–only simulations from 1920 to 2005. The blue line in (A) represents the PC1 from the CMIP6 model VA forcing–only simulations from 1920 to 2012. The red line in (A) represents the global mean of stratospheric AOD at λ = 550 nm as implemented in CMIP5 model simulations. The data near the right end were derived with mirrored data in the filtering and thus are less reliable; they are marked by the dashed lines. (C and D) Same as (A) and (B) but for the MMM of CMIP5 AA forcing–only simulations from 1920 to 2005. The blue line in (C) represents the PC1 from the CMIP6 AA forcing–only simulations. The EOF patterns for the CMIP6 VA and AA single simulations are comparable to those shown in (B) and (D), respectively. The red (pink) line in (C) represents the smoothed ambient AOD at λ = 550 nm averaged over the near-global oceans from the CMIP5 (CMIP6) AA forcing–only simulations.

The AA forcing and its model-simulated SST response also exhibit multidecadal variations since 1920, with negative SST anomalies from the late (for CMIP5) or early (for CMIP6) 1960s to the 1990s and positive anomalies outside this cold period (but after ~1930; Fig. 4, C and D). This forced signal is also roughly in phase with the observed AMO (r ≈ 0.78). The local correlations between the CMIP5 model-simulated (under VA or AA forcing) and observed SST decadal-multidecadal variations (fig. S3) reveal that the VA- and AA-induced SST variations are evident in the observations mostly over the North Atlantic and most of the Pacific except the tropical central and eastern Pacific. In particular, the Pacific SST response to the VA forcing reflected in the observations shows a La Niña–like pattern in the Pacific (fig. S3).

The AA forcing influences NASSTs mainly through its direct and indirect (i.e., through changes in clouds) effects on surface solar radiation (10), although CMIP5 models still have large uncertainties in simulating these effects from aerosols (24, 25). Our analyses (Fig. 5) also indicate that this is the case in the North Pacific, as surface solar radiation change patterns roughly match the SST change patterns. The opposite changes in solar radiation in the tropical western and subtropical Pacific and Indian oceans are due to the reversed changes of local aerosol loading (Fig. 5D). However, the linkage among the solar radiation changes, clouds, and SSTs (through indirect effects) in response to the AA forcing appears to be weak in the tropical central and western Pacific with little influence by the reversed aerosol loading (Fig. 5C).

Fig. 5 Differences in spatial response between warm and cold phases.

Annual-mean anomalies for the cold period (1965–1995) relative to the warm period (1935–1960) for (A) surface winds (vectors, m s−1) and SST (°C, shading), (B) rainfall (mm day−1, shading) and SLP (hPa, contours; black line represents the zero line, and purple and red lines indicate the negative and positive contours, respectively), (C) surface net downward solar radiation (W m−2), and (D) net clear-sky surface downward solar radiation (W m−2) from the CMIP5 MMM AA forcing–only simulations. The model simulations were linearly detrended before calculating the anomalies. The black (gray) stippling indicates that the anomalies are statistically significant at the 0.05 (0.1) level based on Student’s t test. Note the SST anomalies in (A) resemble those shown by Fig. 4D with the sign reversed.

Another possible cause for the decadal SST changes in the tropical Pacific is the response of the large-scale atmospheric circulation to the aerosol forcing (Fig. 5B and fig. S4). Despite geographically uneven emissions of AAs, the climate response is relatively insensitive to the spatial distribution of the aerosol concentrations (26). Along with the increasing AA loading over global oceans (Fig. 1A), the spatial pattern and long-term trend caused by the aerosol forcing in SSTs, rainfall, and atmospheric circulation exhibit similar features (but with an opposite sign) to those resulting from GHG changes (27). In the equatorial Pacific, increased AAs lead to lower SSTs and decreased rainfall (Fig. 5, A and B), following the opposite of the warmer-get-wetter mechanism (26). The North Pacific shows high sensitivity to the AA forcing, with a pronounced SST cooling and intensified Aleutian low (Fig. 5B). Previous studies also found a pronounced interhemispheric asymmetry in aerosol-forced simulations (27). In the South Pacific, the southward shift of the tropical rain belt is associated with southward cross-equator surface winds, increased local rainfall, and weakened southeasterly trade winds (Fig. 5, A and B). Surface westerly winds weaken over the Southern Ocean in response to aerosol forcing, although the local aerosol forcing is small (Fig. 5A); this suggests a remote response through atmospheric circulation, such as the effect of the interhemispheric anomalous Hadley circulation (26). In addition, the AA forcing causes rainfall to decrease over East Asia, the Sahel, and South America, which is consistent with previous studies (Fig. 5B) (11, 27). Overall, the spatial characteristics of the climate response shown in Fig. 5 are similar to those shown in previous studies (26, 27), and they show an in-phase, global response to the AA forcing. However, our study extends beyond the previous studies (26, 27) by showing the multidecadal variations in SSTs in global oceans, and in regional climates and large-scale atmospheric circulations resulting from AAs (Fig. 5).

Furthermore, atmospheric circulation patterns in the pan-tropics show significant decadal changes, with negative (positive) sea level pressure (SLP) anomalies over the tropical and subtropical Pacific (tropical and subtropical Atlantic) under increased aerosols (Fig. 5B). The zonally varying SLP anomalies are associated with an intensified Walker circulation in the Pacific (fig. S4) and increased surface easterlies over the Equator in the eastern Pacific during the high AA and cold period (Fig. 5A), leading to cooling there or a La Niña–like pattern under increased AA (Fig. 5A). Previous studies (28, 29) suggest that enhanced Atlantic tropical warming since the early 1990s is a key driver of the intensification of the Pacific Walker circulation, eastern Pacific cooling, and strong easterly winds. However, the cause of the recent Atlantic warming is under debate as the AMO may also affect recent Atlantic SSTs, besides GHG-induced warming. Here, we show that the historical AA forcing can produce zonally varying tropical SST and SLP anomalies that can alter the Walker circulation. Our results suggest that the local indirect and direct effects of aerosols on surface radiation dominate over the North Atlantic and North Pacific, whereas atmospheric circulation response to the AA forcing plays an important role in the SST decadal variations in the tropical Pacific.

Note that Asian aerosols show a monotonic increase with little decadal variation (fig. S5); their impact is largely removed by the detrending, leaving mainly the influence from the aerosols emitted from North America and Europe (fig. S6). Therefore, the AMO-like SST response is mainly forced by decadal variations in AAs emitted from North America and Europe, while Asian aerosols have little decadal variation and thus contribute little to the AMO-like SST response. As a result, the AA-forced global SST change pattern remains largely unchanged when the effect of Asian aerosols is excluded (fig. S6).

Because the effect of the AA forcing remains a live controversy, and excessive indirect aerosol effects (AIEs) may exist in some CMIP5 models (24, 25), we also analyzed the new CMIP6 multimodel ensemble simulations (table S2) that included the AIE forcing to examine its impact. Consistent with the CMIP5 model results, in the CMIP6 AA forcing–only simulations, the direct and AIE work together to generate the decadal SST changes in the North Pacific and Atlantic (fig. S7), and they also produce broadly similar atmospheric circulation changes in the Pacific and North Atlantic (figs. S4 and S7). However, the spatial pattern of aerosol loadings and the magnitudes of the responses vary between the two generations of models. For example, the different surface solar radiation responses (Fig. 5, C and D, and fig. S7, C and D) over the equatorial Pacific and south Atlantic may be due to reversed aerosol changes and distributions in South America. The Pacific Walker circulation anomaly is relatively weak in the CMIP6 models (fig. S4), the anomaly pressure center is displaced northward in CMIP6 over the Atlantic, and the SLP anomalies differ from the CMIP5 models over the Northwest Pacific (fig. S7). Although the AIE is included in these CMIP5 and CMIP6 simulations, their differences presumably result from different model parameterizations of the cloud microphysical processes and different spatial patterns in aerosol loadings (e.g., over TCP, SNP, and SSP) between the CMIP5 and CMIP6 simulations (15, 25). Given the large spread of the SST and circulation response to aerosols among individual models, further studies are needed to understand aerosols’ effects on our recent climates.


The above analyses suggest that there exist externally forced AMO-like in-phase multidecadal SST variations since 1920 not only in the North Atlantic but also across the three oceans in both model simulations and observations, although their magnitude and correlations with observed SSTs vary spatially due to varying contributions from internal variability such as IPO. The AMO-like variations result from decadal variations in external VA and AA (mainly from North America and Europe) forcing, which reduces surface solar radiation and cools SSTs during 1965–1995 relative to 1935–1960 or after the late 1990s. The aerosol forcing also causes zonally varying SST and SLP changes in the tropics that may have contributed to Pacific SST response to the aerosol forcing.

Our new finding is consistent with previous analyses of reconstructed paleo-AMO records that suggest increased AMO amplitudes since the industrial revolution (30), as the recent VA and AA decadal variations happen to be in phase with the internally generated AMO and thus contribute to about 42% to the observed AMO during 1950–2012 (fig. S2A). Our wavelet power spectral analyses of the reconstructed AMO and Niño3.4 SST indices of the last 2000 years (fig. S8) (31) show significant power around periods of 60 to 80 years for the AMO index, and this power has increased substantially since ~1700, consistent with (30). In contrast, the Niño3.4 index does not show increased power for typical IPO periods of 20 to 60 years during the past two centuries. These paleoclimate results provide another line of evidence that external forcing has enhanced the multidecadal variations in the North Atlantic (and possibly in other oceans too) since the industrial revolution, whereas the decadal-multidecadal variations in the TCP are dominated by internal variability or natural forcing.

Our study provides a global perspective on how external forcing agents (namely, VA and AA) have influenced global SSTs in models and observations. Although previous studies (10, 15) have suggested that decadal variations in aerosol forcing may have caused SST changes in the Pacific and North Atlantic, here, we show that there exist common, AMO-like in-phase multidecadal variations across the three oceans in both observations and model simulations, and it is caused by decadal changes in both VAs and AAs, not by the internally generated AMO in the North Atlantic through atmospheric teleconnection, as previously thought (8). This finding provides a new perspective for studying recent interbasin connections in SST and other climate fields on decadal to multidecadal time scales, because many of the AMO-like variations outside the North Atlantic may result from the common aerosol forcing rather than through atmospheric teleconnections linked to NASSTs as examined previously (8, 32). It also suggests that we may need to reassess the attribution of many AMO-like climate variations seen outside the North Atlantic, such as those in the TWP, Eurasia, North America, and other regions (5, 6, 8, 18, 19, 32). As we demonstrated for TWP SSTs, these AMO-like variations may be caused, at least partly, by the AMO-like aerosol forcing, rather than by the internally generated multidecadal SST variations in the North Atlantic, although it is yet to show whether and how the aerosol forcing can cause AMO-like climate variations over land areas.

Since its discovery (3335), AMO has been used to refer to the total multidecadal variations in NASSTs, which include both internally generated and externally forced components. Because of their different nature of origin, it is better to distinguish and separate those components. Some authors (36) suggested to use AMO to refer to the internal component of the NASST multidecadal variations only, while others (37) suggested to use Atlantic Multidecadal Variability (AMV) to depict the total multidecadal variations in NASSTs from both internal variability and external forcing to avoid confusion. However, how to precisely separate and define those two components requires further work.


Data and model simulations

We used four gridded monthly SST datasets, including the Extended Reconstructed SST version 5 (ERSSTv5) (38), HadISST1 (39), Kaplan Extended SST (40), and COBE-SST2 (41). These observational SST datasets share similarities in their spatial pattern of the correlation with the model-simulated AMO-like signal, although they differ substantially in their input data sources and data processing approaches. To reduce uncertainty and redundancy, we used their EM to represent the observed SSTs in this study. We focus on the period starting from 1920, because SST observations before about 1920 are sparse over the tropical and southern oceans (1). Furthermore, we exclude the Southern Ocean (south of 40°S) as the SST observations are less reliable there (1).

We analyzed the coupled climate model simulations from the CMIP5 (42), including 121 all-forcing historical (1850–2005) and Representative Concentration Pathway 4.5 (RCP4.5) future (2006–2100) simulations from 38 models (table S1). The historical forcing agents include time-varying GHGs, tropospheric and stratospheric aerosols, ozone, land use, and solar irradiance derived from observations. The effects of these individual forcing agents were investigated by analyzing single-forcing simulations forced separately by historical GHG, AA, natural (VAs plus solar irradiance) forcing, VA, and solar irradiance.

We also analyzed the large ensemble of simulations by the CESM1 (43). This ensemble includes 40 coupled simulations forced by all the historical forcing agents from 1920 to 2005 and RCP8.5 forcing from 2006 to 2100 with about 1° grid spacing. Each of the 40 simulations has the same external forcing but starts from slightly perturbed initial conditions that often lead to different and thus uncorrelated realizations of the internal variability (including IPO and AMO) among the ensemble simulations. We also used an 1800-year fully coupled preindustrial control run from the CESM1 model.

Analysis methods

We used the CMIP5 MMM or CESM1 EM as the externally forced signal, because uncorrelated internal variations among the ensemble runs were largely smoothed out in this averaging over a large number of ensemble simulations (36, 44). The all-forcing MMM was calculated by averaging over all the 121 model runs from the 38 CMIP5 models. The MMM was similar if we first calculated the EM for each model by averaging over all the available model runs for the model (which range from 1 to 17; table S1) and then averaged over the EMs over the 38 models with equal weighting. We used the arithmetic mean of the 40 members from the CESM1 large ensemble as the CESM1 EM. To separate externally forced and internally generated variations in observations, following (44), we removed the changes and variations at each grid point associated with the forced signal, which is represented by the global-mean SST time series from the CMIP5 MMM or CESM1 EM, from the observed SST time series via linear regression during 1920–2012. The regression accounts for any systematic bias in the model-simulated response to external forcing. The use of the global-mean, instead of local, SST from the MMM or CESM1 EM in deriving the forced signal further removes any remaining local internal variations in the MMM or CESM1 EM (which can be substantial for an ensemble size less than 100), as the forced signal is in phase everywhere with the same external forcing (45). After removing the forced component, the residual SST fields were considered as containing primarily internal variability.

Similarly, to remove the warming signal induced by GHG forcing only from the observed, all-forcing MMM, or CESM1 EM SSTs, we used the global-mean SST time series derived from the MMM of the CMIP5 GHG forcing-only simulations (from 1920 to 2012, 35 runs from nine models; table S1) as the GHG-forced signal to detrend the observational or modeled SST fields via linear regression at each grid point. The removal of the GHG-induced warming trend allows us to focus on the relationship between observed and model-simulated externally forced decadal-multidecadal variations. The residual (referred to as the GHG-detrended SSTs) in observations represents non-GHG components that contain internally generated variations and externally forced changes by aerosols and other forcing agents. After removing the GHG-forced component, we performed an EOF analysis of the residual SST fields from the all-forcing MMM or CESM1 EM to examine spatial and temporal characteristics of the leading modes in the forced decadal-multidecadal variations. Note that the aerosol increasing trend (gray line in Fig. 1A) should produce a cooling trend that offsets some of the GHG-induced warming trend in the all-forcing simulations, with the net warming trend being removed during the GHG-based detrending. This allows us to focus on the aerosol-induced multidecadal variations (without the cooling trend) in the GHG-trended SSTs. We applied a low-pass Lanczos filter (with 19 points and a 13-year cutoff period) to annual-mean SSTs and other variables to emphasize decadal to multidecadal variations (17). We also tested the 13-year moving averaging, and the results were qualitatively similar.

We analyzed the observed SSTs using the EM of the four widely used observational SST datasets (namely, ERSSTv5, HadISST1, Kaplan SST, and COBE-SST2) to reduce uncertainty and redundancy. To obtain the observed decadal SST variations, we removed the GHG-induced warming signal from the observed SSTs through linear regression using the global-mean SST time series derived from the CMIP5 MMM of the GHG-only simulations and then applied a low-pass Lanczos filter. To derive an estimate of internally generated variability in observations, we removed the changes and variations at each grid point associated with the forced signal, which is represented by the global-mean SST time series from the CMIP5 (CESM1) all-forcing historical MMM (EM) from the observed SST time series via linear regression. After removing the externally forced component, the residual SST fields were considered as arising from internal variability. Furthermore, we tested different methods to estimate the externally forced component in observed SSTs, and they yielded similar results (see Supplementary Text).

To examine the regional SST changes in atmospheric circulation and radiation fluxes, we performed composite analyses over the positive and negative phases of the AMO-like signal to reveal any coherent SST anomalies, and the associated atmospheric circulation and radiation changes over the global oceans. A linear correlation coefficient (r) and its significance level (P value) were calculated between two time series to quantify their association. The significance level was estimated based on a two-sided Student’s t test with an estimated effective degree of freedom to account for autocorrelation (46, 47). We also examined the statistical significance based on Monte Carlo approach (see detailed information in Supplementary Text).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP5 and CMIP6. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP5, CMIP6, and ESGF ( We also acknowledge the CESM Large Ensemble Community Project and supercomputing resources provided by NSF/CISL/Yellowstone. Funding: This study was partly funded by the National Key R&D Program of China (2016YFA0600702 and 2017YFC1502101). A.D. was supported by the U.S. National Science Foundation (grant no. OISE-1743738) and the U.S. National Oceanic and Atmospheric Administration (award no. NA15OAR4310086 and NA18OAR4310425). This study was also supported by the NUIST-UoR International Research Institute and the Innovation Research Program for College Graduates of Jiangsu Province (grant no. SJKY19_0931). Author contributions: A.D. and W.H. designed the study, contributed ideas to the analysis, and improved the manuscript. M.Q. carried out the analysis and wrote the first draft of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article