Natural variability has dominated Atlantic Meridional Overturning Circulation since 1900

There is debate about slowing of the Atlantic Meridional Overturning Circulation (AMOC), a key component of the global climate system. Some focus is on the sea surface temperature (SST) slightly cooling in parts of the subpolar North Atlantic despite widespread ocean warming. Atlantic SST is influenced by the AMOC, especially on decadal timescales and beyond. The local cooling could thus reflect AMOC slowing and diminishing heat transport, consistent with climate model responses to rising atmospheric greenhouse gas concentrations. Here we show from Atlantic SST the prevalence of natural AMOC variability since 1900. This is consistent with historical climate model simulations for 1900–2014 predicting on average AMOC slowing of about 1 Sv at 30° N after 1980, which is within the range of internal multidecadal variability derived from the models’ preindustrial control runs. These results highlight the importance of systematic and sustained in-situ monitoring systems that can detect and attribute with high confidence an anthropogenic AMOC signal. The Atlantic Meridional Overturning Circulation (AMOC) is predicted to slow with climate change. Sea surface temperature data and climate model analysis show that since 1900 natural variability has been dominant in AMOC changes; anthropogenic forcing is not yet reliably detectable by this method.

radiative forcing (ERF; Methods and Figs. 2a and 3a), however, strongly increased from the 1970s onward. We note that the strong year-to-year variability in the SST of the NAWH is mostly due to atmospheric heat-flux forcing linked to the winter NAO 57 . This short-term variability is not addressed here, as the focus of this study is on the long-term variability.

Spatial structure of long-term SST change
Empirical orthogonal function (EOF) analysis is applied to the observed annual SST anomalies (blue boxes in Figs. 2b and 3b; Methods), which is a variance-maximizing multivariate statistical technique 58 . The two leading modes, EOF1 and EOF2, are dominated by long-term variability and are of interest here. EOF1 accounts for 52.2% of the total Atlantic SST variability. The time evolution (principal component) of EOF1, PC1 (Fig. 2a), is governed by an accelerating upward trend with substantial superimposed less-energetic interannual-to-multidecadal variability. PC1 follows the long-term evolution of the ERF well, with a correlation coefficient of 0.83 when the Shared Socioeconomic Pathway 5-8.5 (SSP5-8.5; Methods) scenario is considered during 2015-2019. We suggest that EOF1 describes the externally forced Atlantic SST variability, in particular the effects of global warming. The prominent, short-lived, downward spikes in the ERF during the second half of the twentieth century, however, are not well captured by PC1 and have been attributed to explosive volcanic eruptions 59 .
The pattern of EOF1, shown as the local regression coefficients upon PC1 (Fig. 2b), exhibits positive anomalies over the Atlantic basin, with the exception of the NAWH. PC1 accounts for up to 80% of the SST variability in the tropical Atlantic but hardly any in the NAWH. In most other regions of the Atlantic, the explained variances amount to at least 50%. EOF1 explains a relatively large fraction of the SST variability over most of the global ocean. We note the warming hole in the subpolar Southern Ocean ( Fig. 2b and Extended Data Fig. 1a); there, stronger westerly surface winds (associated with stratospheric ozone loss 60,61 ) or Antarctic meltwater 62 may have offset GHG-forced warming. Moreover, in climate models, the internal long-term variability is particularly strong in both the Southern Ocean and NA [63][64][65][66][67] , hindering detection of externally forced signals.
The second most energetic mode (EOF2; Fig. 3) accounts for 12.1% of the total Atlantic SST variability. The corresponding time series, PC2 (Fig. 3a), has no obvious relationship with the ERF and exhibits pronounced decadal-to-multidecadal variability but no trend. EOF2 is picking up the signal of the Atlantic interhemispheric dipole, which has been linked to long-term AMOC variability in climate models 29,30,44,[68][69][70] . Centres with opposite sign are observed in the subpolar NA and in the South Atlantic (Fig. 3b), where EOF2 typically accounts for up to 60% and 20% of the local SST variability, respectively. There is a highly significant correlation amounting to 0.7 between PC2 and an Atlantic interhemispheric dipole index as defined by Latif et al. 19 . EOF2 explains considerably less variance in the tropical Atlantic where EOF1 explains the most. The salient feature of PC2 (Fig. 3a) is the marked decline during 1930-1970. We note that the multidecadal decline in PC2 was preceded by a multidecadal decline in the winter NAO index (Extended Data Fig. 2). However, the largest significant correlation, which is negative, is observed with no lag (Extended Data Fig. 3). We hypothesize that EOF2 describes the internal, AMOC-related, SST variability in the Atlantic.

Historical climate model simulations
The observed NA-SST variability is within the range of the historical simulations with state-of-the-art climate models from the Coupled Model Intercomparison Project phase 6 (CMIP6; Methods and Fig. 1a-c), employing estimates of external forcing from 1900 to 2014 (ref. 71 ). A measure of the externally forced variability is given by the ensemble mean, that is, by averaging over all simulations. The ensemble-mean SST variability is much weaker than the observed variability in the NAWH (Fig. 1a) and subpolar NA (Fig.  1b), suggesting the prevalence of internal SST variability in these two regions. The observed SST variability is also underestimated by the ensemble mean (but to a lesser extent) when averaging over the whole NA (Fig. 1c), supporting a major role of external forcing in driving tropical and subtropical Atlantic SST. This is consistent with EOF1, which is closely linked to the ERF (Fig. 2a), explaining the most variance in the tropics and subtropics. The ensemble-mean NA SSTs lead the observed SSTs by about a decade in all three index regions ( Fig. 1a-c). We note that historical external forcing is subject to considerable uncertainties. Aerosol-cloud interaction is the largest source of uncertainty in estimating historical anthropogenic radiative forcing, and representing this interaction in climate models poses a major challenge 72,73 . Finally, in the regions of the three SST indices, the 40-year linear SST trend during 1975-2014 lies within the trend distribution derived from the preindustrial control integrations of the models (Extended Data Fig. 4). An EOF analysis is applied jointly to the ensemble-mean Atlantic SST and Atlantic meridional overturning streamfunction (MOC). The SST pattern (Fig. 4a) of the leading mode, EOF1 mod , accounting for 53.6% of the joint variance is similar to that of EOF1 (Fig. 2b), with major warming over most of the Atlantic and a warming hole in the subpolar NA south of Greenland. Variances in SST explained by EOF1 mod are generally large, exhibiting a maximum in the subtropical NA with values in localized regions exceeding 90%, and a minimum south of Greenland with values of less than 10% (Fig. 4a). The MOC pattern of EOF1 mod exhibits negative loadings north of the Equator, indicating weaker overturning, that are centred near 40° N in the depth range of 1,000-2,000 m (Fig. 4b). Positive loadings are observed south of the Equator. With regard to MOC, the explained variances are largest in the centre of the negative loadings, where they amount to about 80% (Fig. 4b). Regions of large explained variances indicate a high model consistency.
In climate models, MOC anomalies linked to changes in North Atlantic Deep Water formation first appear in the subpolar NA and then propagate southward 74 . Fully developed AMOC slowing is characterized by basin-wide negative streamfunction anomalies in the models 13 . EOF1 mod is thus consistent with the initial stage of AMOC slowing. Because EOF1 mod is the leading mode of the externally forced variability in SST and MOC in the models, PC1 mod (Fig.  4c) is expected to be significantly correlated with the ERF. The correlation coefficient between the two time series amounts to 0.79. PC1 mod is governed by multidecadal variability until about 1980, and it exhibits a sustained upward trend thereafter. Since 1980, there is The dots indicate significance at the 95% level. The blue box indicates the region over which the EOF analysis was performed, and the black box marks the region over which the NAWH index ( Fig. 1a) is defined. The dots indicate significance at the 95% level. The blue box indicates the region over which the EOF analysis was performed, and the black box marks the region over which the NAWH index (Fig. 1a) is defined.
a reduction in the AMOC transport at 30° N that amounts to about 1 Sv (10 6 m 3 s −1 ), which is assumed to be within the range of natural variability. This is supported by the distribution of non-overlapping 40-year MOC trends calculated from the preindustrial control integrations of the CMIP6 models (Fig. 4d), providing estimates of the internal MOC variability. The ensemble-mean reduction in ocean heat transport at 50° N is about 0.03 PW. We compare the AMOC index at 26.5° N reconstructed from EOF1 mod with the RAPID data 75 (Extended Data Fig. 1f). The 2004-2014 trend from RAPID amounts to −4.1 Sv and that from the reconstruction to −0.5 Sv. When assuming that the CMIP6-ensemble mean realistically represents the externally forced MOC variability, the trend from RAPID must contain a strong contribution from long-term internal variability. This is consistent with Weijer et al. 11 , who analysed another subset of the CMIP6 models, reporting that linear trends over the entire 1850-2014 historical period are generally neutral.

Discussion
Observed SSTs and a large ensemble of historical simulations with state-of-the-art climate models suggest the prevalence of internal AMOC variability since the beginning of the twentieth century. Observations and individual model runs show comparable SST variability in the NAWH region. However, the models' ensemble-mean signal is much smaller, indicative of the prevalence of internal variability. Further, most of the SST cooling in the subpolar NA, which has been attributed to anthropogenic AMOC slowing 21 , occurred during 1930-1970, when the radiative forcing did not exhibit a major upward trend. We conclude that the anthropogenic signal in the AMOC cannot be reliably estimated from observed SST. A linear and direct relationship between radiative forcing and AMOC may not exist. Further, the relevant physical processes could be shared across EOF modes, or a mode could represent more than one process. A relatively stable AMOC and associated northward heat transport during the past decades is also supported by ocean syntheses combining ocean general circulation models and data 76,77 , hindcasts with ocean general circulation models forced by observed atmospheric boundary conditions 78 and instrumental measurements of key AMOC components 9,22,79-81 . Neither of these datasets suggest major AMOC slowing since 1980, and neither of the AMOC indices from Rahmstorf et al. 20 or Caesar et al. 21 show an overall AMOC decline since 1980.
An important remaining issue is the question of how well the externally forced part of the Atlantic SSTs can be determined. Standard EOF analysis has been used above to identify the forced component. We additionally used two other methods. First, we applied principal oscillation pattern (POP) analysis (Methods) to the observed SSTs, again only over the Atlantic Ocean. The results of the POP analysis (Extended Data Fig. 5) are similar to those of the EOF analysis (Fig. 2), in that the time series of the leading POP mode (POP1) is highly correlated with the ERF (r = 0.85; Extended Data Fig. 5a) and POP1 (Extended Data Fig. 5b) is similar to EOF1 ( Fig. 2b). Second, we derived the externally forced component from the historical simulations with the CMIP6 models by applying signal-to-noise ratio (S/N) maximizing EOF analysis (Methods) to the modelled SSTs over the global ocean (Extended Data Fig. 6). As expected, the time series of the leading mode (S/N EOF1) is highly correlated with the ERF (r = 0.9, Extended Data Fig. 6a). The associated regression pattern (Extended Data Fig. 6b) is similar to that linked to both EOF1 (Fig. 2b) and POP1 (Extended Data Fig. 5b), indicating little sensitivity to the choice of the statistical method in determining the externally forced SST variability. We subtracted the variability associated with S/N-EOF1 from the observed SSTs to estimate the internal SST variability, and standard EOF analysis was applied to the residuals over the Atlantic. The second most energetic EOF mode (EOF2 res ; Extended Data Fig.  7), accounting for 19.2% of the variability in the residuals, is multidecadal (Extended Data Fig. 7a). EOF2 res is the interhemispheric dipole, with explained variances up to 40% in the NAWH and in the centre of the southern pole (Extended Data Fig. 7b). EOF2 res is similar to EOF2 (Fig. 3), which is reassuring and supports the notion that EOF2 actually reflects internal SST variability.
The leading EOF of the residuals over the Atlantic (EOF1 res ; Extended Data Fig. 8) accounts for 27.1% of the variance. Its time series (PC1 res ) is characterized by a marked centennial change with a minimum in the 1940s (Extended Data Fig. 8a). Over the Atlantic, the regression coefficients upon PC1 res are largest north of the Equator (Extended Data Fig. 8b) and reminiscent of the "horseshoe" pattern linked to the positive phase of the NAO. PC1 res and the annual-mean NAO index exhibit a statistically significant correlation of 0.4.
A number of factors, internal and external, can influence the AMOC. Menary et al. 38 report that the multi-model mean AMOC strengthened by approximately 10% from 1850 to 1985 in historical simulations with CMIP5 and CMIP6 models, which has been attributed to aerosol forcing. This is consistent with our study in that the AMOC slowing largely takes place after 1980. Analysis of single-forcing experiments by Menary et al. 38 reveal that the forced response of the AMOC is a balance of opposing contributions from aerosols, increasing AMOC, and GHGs, decreasing AMOC. Besides surface heat flux, surface freshwater input from Greenland ice melt is thought to be an important factor causing anthropogenic AMOC slowing. Meltwater forcing is not considered in the historical simulations with the CMIP6 models, which may cause underestimation of AMOC slowing. A high-resolution ocean general circulation model study finds that the present meltwater forcing from the west Greenland Ice Sheet is not sufficiently large to drive significant reductions in North Atlantic Deep Water formation and thus AMOC strength 82 . Climate model sensitivity to external forcing, however, is a long-standing issue, also with regard to the AMOC's sensitivity to freshwater forcing 83 . In summary, our results reinforce the need for systematic and sustained in-situ AMOC observation systems to detect with high confidence externally forced AMOC slowing 9,84,75 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-022-01342-4.
Mixed-layer depth climatology with 2° × 2° resolution is from de Boyer Montégut et al. 87 and based on a density threshold of 0.03 kg m −3 .
Annual satellite sea-level data for 1993-2019 are from Copernicus (https://cds. climate.copernicus.eu). In the EOF analysis, the Sea Level Anomaly product was used, which is defined relative to the mean of 1993-2012 (ref. 88 ).
Net ERF. Net ERF (expressed in W m −2 ) is the globally averaged net downward radiative flux at the top of the atmosphere after allowing for atmospheric temperature, water vapour and clouds to adjust but with surface temperature or a portion of surface conditions unchanged 1 . The ERF data are from Smith et al. 89  EOF analysis. EOF analysis 58 is a multivariate statistical method. The EOFs are the eigenvectors of the covariance matrix, and they are sorted in descending order by the explained variance. Time series, termed principal components (PCs), are obtained by projecting the original data onto the EOFs. In the joint EOF analysis of the ensemble-mean Atlantic SST and MOC, each variable was normalized by its own field sum of standard deviations, making each field have identical total variance. POP analysis. POP analysis is a multivariate statistical method and is defined as the normal modes of a linear dynamical representation of the data in terms of a first-order autoregressive-vector process with residual noise forcing [91][92][93] . For practical purposes, the original process is usually reduced into the subspace of leading EOFs. Here we use the first 10 EOFs accounting for 90.8% of the total variance.

SST-regression maps.
We show global maps of local linear regression coefficients of SST upon the principal components (PCs) that have been normalized by their respective standard deviation. An F test is used to assess the significance of the regression coefficients.
S/N-maximizing EOFs. S/N-maximizing EOF analysis refers to a method of identifying the fingerprints of external forcing in an ensemble of forced climate model experiments. The method allows us to distinguish between the response to prescribed external forcing common to all ensemble members and internal variability, which is temporally uncorrelated between ensemble members. We follow the formulation of Venzke et al. 94 . The leading EOF has the largest S/N, and the corresponding PC represents the time evolution of the most dominant forced response.
Significance of correlations. The Student's t-test and Monte Carlo simulation based on nonparametric random phase are applied to assess the statistical significance of the correlation coefficients 95 . All correlation coefficients mentioned in the main text are statistically significant at the 95% confidence level.

Code availability
The figures were generated using MATLAB and m_map (https://www.eoas.ubc. ca/~rich/map.html), where the basemap data (GSHHG, https://www.ngdc.noaa. gov/mgg/shorelines/gshhs.html) are used under the GNU Lesser General Public license. All the codes used in the data processing and visualization are available via Figshare at https://doi.org/10.6084/m9.figshare.19318004. Other code is available upon request to the corresponding author J.S.