The Atlantic Meridional Overturning Circulation (AMOC), a tipping component of the climate system, is projected to slowdown during the 21st century in response to increased atmospheric CO2

concentration. The rate and start of the weakening are associated with relatively large uncertainties. Observed sea surface temperature-based reconstructions indicate that AMOC has been

weakening since the mid-20th century, but its forcing factors are not fully understood. Here we provide dynamical observational evidence that the increasing atmospheric CO2 concentration

affects the North Atlantic heat fluxes and precipitation rate, and weakens AMOC, consistent with numerical simulations. The inferred weakening, starting in the late 19th century, earlier

Arctic sea ice, Greenland ice-sheet, among others). Such a tipping element, AMOC, has a quasi-global impact and played a central role in past abrupt climate changes2,3. Its fate during the

twenty first century is a topic of major scientific and socio-economic interest.

Most climate projections indicate that AMOC will suffer a centennial scale slowdown during the twenty-first century, mainly in response to intensified North Atlantic freshwater and heat

fluxes, induced by increased atmospheric CO2 concentration4. However, there are significant quantitative differences between the results of various model simulations, which imply large

uncertainties regarding the future evolution of AMOC.

Numerical integrations over the historical period simulate a modest AMOC decrease of 1.4 ± 1.4 Sverdrup (Sv) between preindustrial (1850–1900) and present day (2006–2015), which is most

pronounced during the last decades, indicating that anthropogenic warming may have already weakened it5. Direct observations indicate that the overturning weakened by 30% during the

1957–2004 period6 and that it was in a relatively weak state between 2008 and 20177. Alternatively, indirect measures of ocean circulation changes, based on historical sea surface

temperature (SST) fields, suggest that AMOC has been weakening since the mid-twenty century8. When calibrated with an ensemble of model simulations from the CMIP5 project, the weakening over

1870–2006 period was estimated to be of 3 ± 1 Sverdrup (Sv) (15%) and it was most pronounced since the mid-twenty century9. Proxy and reconstructed data suggest that the reduced AMOC

intensity during the 1975–1995 period is at the lowest level in the last millennium10. The attribution of this reconstructed trend to external or internal factors remains an open problem of

fundamental importance in climate research.

Here we investigate a potential contribution of atmospheric CO2 concentration to AMOC slowdown, based on observational and reanalysis data. First, we separate the SST-based reconstructed

long-term AMOC weakening trend and emphasize the routes through which this greenhouse gas could affect AMOC. Then we probe the associated causal chains using the Convergent Cross Mapping

(CCM) technique, a method based on the theory of dynamical systems used to identify causation in weakly coupled systems based on two timeseries11.

In practical terms, CCM causation is tested using the technique of “cross mapping”: a time delay embedding is constructed from the time series of Y and the ability to estimate the values of

X from this embedding quantifies how much information about the later has been encoded into the former variable. The accuracy of the prediction is measured using Pearson’s correlation

coefficient (ρ), between observed and predicted values. One notes the counterintuitive fact that the cross map estimate runs in the reverse direction of causality: if Y predicts X, then X

causes Y. However, one key property that distinguishes causality from mere correlation is the convergence of the cross estimation. When constructing the embedding, only a given portion of

the time series is used. Increasing this library length should improve the accuracy of the prediction, since the additional points fill in the trajectories in the attractor, resulting in

closer nearest neighbors. At some point, the information contained in the affected variable has been exhaustively harnessed and the cross map saturates to a plateau. The asymptotic increase

of the cross-map skill with library length is called convergence. The strength of the causal interaction may be linked with the rate of growth toward the convergence level with library size,

but also with the level of the cross map skill.

In order to investigate the AMOC response to this greenhouse gas, in a preliminary analysis, we aim to increase the signal-to-noise ratio by separating the centennial from multidecadal

overturning variations emphasized in previous numerical and observational studies8,12,13, with the former having the same characteristic time scale as the increasing trend of the atmospheric

CO2 concentration4. The separation is performed through two EOF analyzes performed on North Atlantic and South Atlantic annual SST anomalies from the ERSST.v5 dataset, extending over the

1854–2016 period14 (Supplementary Figs. 1–4). Similar results are obtained if the EOF method is applied to the Hadley SST fields15.

Whereas the time component of one North Atlantic mode is marked by a decreasing centennial-scale trend (hereafter Trend Mode—TM), which starts around 1890, that of the other one is dominated

by multidecadal fluctuations, typical for the Atlantic Multidecadal Oscillation (AMO) (Fig. 1b). While the TM pattern is dominated by three centers of alternating signs, disposed from SW to

NE of Greenland (Fig. 1a), a structure which is linked to AMOC changes9,10,16,17, the other mode has a typical monopolar AMO structure (Supplementary Fig. 2a), reflecting multidecadal AMOC

variations18,19. The superposition of the time components of the two North Atlantic modes is strongly correlated (0.77) with a SST-based reconstruction9 (Fig. 1e). The time series of the

corresponding modes identified based on the South Atlantic SST fields show similar centennial-scale and multidecadal variations (Fig. 1d). The South Atlantic TM pattern has a north-south

oriented dipolar pattern (Fig. 1c). The other AMOC mode is marked by multidecadal fluctuations (Fig.1d) and by monopolar structures in both hemispheres of the Atlantic basin (Supplementary

Figs. 2a and 4c), which are typical for AMO8,20. The TM(AMO) derived from North Atlantic SSTs has a 0.65 (0.49) maximum correlation with the corresponding mode obtained from South Atlantic

fields, when the former leads the later by 19 (9) years.

The pattern of the SST Trend Mode (°C) (a) and the associated time series (dark blue line), together with the AMO time series (turquoise line) (b), derived through EOF analysis of North

Atlantic sector (80°W-0°E, 0°N-80°N) annual SST fields. The corresponding pattern and time series, for the South Atlantic sector, are shown in panels c and d. The sums of the two timeseries

in panel b and that of the pair of time components in panel d are shown in panel e (indigo and magenta lines, respectively), together with an AMOC reconstruction9 (orange line). Details on

both EOF analyzes are presented in the supplement.

The superposition of the time components of the two South Atlantic modes is significantly correlated (0.60) with the SST-based reconstruction9 (Fig. 1e). Therefore, these two SST modes are

characterized by distinct spatial and temporal features and are responsible for the largest part of decadal and longer AMOC variability.

A qualitatively identical and quantitatively similar bimodal decomposition of decadal and longer-term AMOC variability over the instrumental period was emphasized in an ensemble of

historical simulations21. When these were combined with control simulations with constant radiative forcing, AMO appears associated with internal AMOC variations, whereas the centennial

scale trend was interpreted as an externally forced mode. One notes that in some of the historical simulations the AMOC centennial-scale decreasing trend starts at the end of the nineteenth

century21 (Fig. 11 in their study), as does the time component of TM, derived from observed SSTs (Fig. 1b). Furthermore, the Atlantic TM and AMO SST patterns were linked to forced and

internal variations, respectively, based on model simulations and observations22. Consequently, hereafter we consider TM as an indicator of externally forced AMOC centennial-scale

variations, separated here from the multidecadal internal fluctuations, reflected by AMO.

It was shown that the average SST anomalies over the subpolar gyre can be translated to AMOC changes, by using a calibration factor of 3.8 ± 0.5 Sv K−1, derived from CMIP5 simulations9. The

fact that the dominant center of TM coincides with the SST anomalies in this area (Fig. 1a), makes possible an estimation of the AMOC change associated with this mode. Consequently, the

linear decreasing trend of TM (Fig. 1b) translates in an AMOC weakening of 3.7 ± 1.0 Sv over the 1854–2016 period.

The North Atlantic SST dipole pattern located south of Greenland in the TM structure (Fig. 1a) was associated with the AMOC slowdown in climate model simulations with increasing atmospheric

CO2 concentration9,23,24,25. This suggests that the TM’s weakening trend, which starts in the late nineteenth century, is induced by this greenhouse gas (Fig. 1b).

Model integrations indicate that an increase of atmospheric CO2 concentration can weaken AMOC through increases of North Atlantic surface heat and freshwater fluxes4,26,27,28. Here we

investigate these simulated connections, based on observational (SST—ERSSTv5 dataset) and reanalysis (Sea Level Pressure (SLP), heat fluxes and precipitation rate—NOAA/CIRES/DOE 20th Century

Reanalysis V3 datasets) fields14,29 (sea Methods).

The regression of spring North Atlantic Ocean heat fluxes on the time series of atmospheric CO2 concentration30 (Fig. 2a) is dominated by positive values (Fig. 2c), consistent with a direct

thermodynamic influence of this greenhouse gas on ocean surface temperature. A quasi-identical pattern is obtained if a composite map is constructed as difference between average values over

the 1935–2015 and 1854–1934 periods (Supplementary Fig. 5). The regression of the North Atlantic spring precipitation rate on CO2 record is dominated by negative anomalies (not shown).

a Reconstruction of annual atmospheric CO2 concentration (ppm); b Time series of the dominant North Atlantic spring SLP mode (black line) and of the SST/AMOC Trend Mode (blue line) (c)

regression map of spring surface heat fluxes on the annual atmospheric CO2 concentration (W/m2/ppm); for positive values the flux is directed from the atmosphere into the ocean (d)

regression map of spring precipitation rate on the SLP time series shown in panel b ([kg/(m2 × s)/Pa] × 106); the dominant North Atlantic SLP mode is also shown (solid black contour lines)

(Pa); e time series of average North Atlantic ocean spring heat fluxes (orange line) and time series of the SST/AMOC Trend Mode (blue line). f Time series of average North Atlantic spring

precipitation rate (violet line) and the time series of the SST Trend Mode (blue line). In panels c and d the regions for which the regression is significant above the 95% level are marked

by red grids.

Model projections show that increasing atmospheric CO2 concentration results in a positive trend of the dominant atmospheric mode in this sector, the North Atlantic Oscillation31. The

pattern and time series of this mode are derived through EOF analysis (not shown). The regression map of the North Atlantic precipitation rate field on the NAO index is dominated by a

prominent center of positive values located in the low-pressure center of this mode, consistent with an influence from the later to the former (Fig. 2d).

A similar regression map of heat fluxes on NAO time series reveals a dipolar structure, which has a small projection on the average value of the North Atlantic sector (not shown). Therefore,

these analyzes suggest influences of increased atmospheric CO2 concentration on spring North Atlantic surface heat fluxes (directly, thermodynamically) and on precipitation rate

(indirectly, dynamically, through NAO), which were shown to weaken AMOC in model simulations4,26,27,28.

The dominant growing character of these two fluxes, associated with increasing atmospheric CO2 concentration (Fig. 2c, d), is reflected also in the upward trends of the integrated North

Atlantic (70°W-20°E, 40°N-80°N) heat flux and precipitation rate time series (Fig. 2e, f). In order to test these observed connections linking increasing CO2 with weakening AMOC, anticipated

by model simulations, we apply CCM on their associated pairs of variables.

TDCCM is first applied to two possible connected time series in order to estimate the lag for which the cross-map estimate is maximum and to infer the correct sense of the causal links, when

an unambiguous interpretation is available11. With the identified lag, CCM is used (Supplementary Fig. 9) to check for convergence and its statistical significance, based on which a

potential causal relationship can be inferred (see Methods).

The cross map skill from North Atlantic spring surface heat flux to atmospheric CO2 concentration increases with the library length and reaches a plateau around \(\rho \cong 0.8\), well

above the 95% significance level (p