# Increased risk of near term global warming due to a recent AMOC weakening

Rémy Bonnet, Didier Swingedouw, Guillaume Gastineau, Olivier Boucher, Julie Deshayes, Frédéric Hourdin, Juliette Mignot, Jérôme Servonnat, Adriana Sima

► To cite this version:

Rémy Bonnet, Didier Swingedouw, Guillaume Gastineau, Olivier Boucher, Julie Deshayes, et al. In-creased risk of near term global warming due to a recent AMOC weakening. Nature Communications, 2021, 12 (1), 10.1038/s41467-021-26370-0. hal-03405246

HAL Id: hal-03405246

https://hal.sorbonne-universite.fr/hal-03405246

Submitted on 27 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-021-26370-0

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# Increased risk of near term global warming due to a recent AMOC weakening

Rémy Bonnet<fnref n="1" />✉, Didier Swingedouw 1✉, ², Guillaume Gastineau³, Olivier Boucher<fnref n="1" />, Julie Deshayes 1✉, ³, Frédéric Hourdin⁴, Juliette Mignot³, Jérôme Servonnat⁵ & Adriana Sima⁴

Some of the new generation CMIP6 models are characterised by a strong temperature increase in response to increasing greenhouse gases concentration<fnref n="1" />. At first glance, these models seem less consistent with the temperature warming observed over the last decades. Here, we investigate this issue through the prism of low-frequency internal variability by comparing with observations an ensemble of 32 historical simulations performed with the IPSL-CM6A-LR model, characterized by a rather large climate sensitivity. We show that members with the smallest rates of global warming over the past 6-7 decades are also those with a large internally-driven weakening of the Atlantic Meridional Overturning Circulation (AMOC). This subset of members also matches several AMOC observational fingerprints, which are in line with such a weakening. This suggests that internal variability from the Atlantic Ocean may have dampened the magnitude of global warming over the historical era. Taking into account this AMOC weakening over the past decades means that it will be harder to avoid crossing the 2 °C warming threshold.

<footnote n="1">¹Institut Pierre-Simon Laplace, Sorbonne Université/CNRS, Paris, France. ²Environnements et Paléoenvironnements Océaniques et Continentaux, Université de Bordeaux/CNRS, Bordeaux, France. ³Laboratoire d'Océanographie et du Climat: Expérimentations et Approches Numériques, Institut Pierre-Simon Laplace, Sorbonne Université/CNRS/IRD/MNHN, Paris, France. ⁴Laboratoire de Météorologie Dynamique, Institut Pierre-Simon Laplace, Sorbonne Université/CNRS/Ecole Normale Supérieure/Ecole Polytechnique, Paris, France. ⁵Laboratoire des Sciences du Climat et de l'Environnement, Institut Pierre-Simon Laplace, CEA/CNRS/UVSQ, Gif-sur-Yvette, France. ✉email: remy.bonnet@ipsl.fr</footnote>
<figure><img src="image_4.png" /></figure>

**Fig. 2** Relationship between GSAT and AMOC over the historical period. **a** Scatter plot of Global near-Surface Air Temperature (GSAT) (K decade⁻¹) versus Atlantic Meridional Overturning Circulation (AMOC) (Sv decade⁻¹) trends calculated over the 1940-2016 period from the IPSL ensemble of extended historical simulations (IPSL-EHS; filled circles) and the observations (dot-dashed lines), with HadCRUT4-CW⁴¹,⁴² for the temperature and the Caesar index¹⁷ as a proxy for the AMOC evolution, with the related uncertainties (gray). The black solid line represents the least square regression between these two variables in IPSL-EHS, with a determination coefficient r² = 0.82 (p < 0.1, see Method). The color scale represents the RMSE between simulated and observed (HadCRUT4-CW⁴¹,⁴²) annual GSAT anomalies over the period 1900-2018. **b** Time evolution of the low-pass filtered AMOC strength anomaly from IPSL-EHS (gray), the AMOC forced signal (or ensemble mean, red), and the subset of 6 members labelled in **a** (light blue), with the subset average (dark blue). The anomaly is calculated for each member with respect to its 1850-2018 average. A Lanczos low-pass filter with a cutoff period of 11 years is used.

<figure><img src="image_5.png" /></figure>

**Fig. 3** Evaluation of Atlantic Meridional Overturning Circulation fingerprints. **a** Time evolution of the Caesar index computed as the low-pass filtered sea surface temperature anomaly (K) averaged over the North Atlantic Subpolar Gyre from November to May from the ERSSSTv5 observational dataset⁵⁴ (blue), the 32 historical members (gray) with the ensemble mean (black) and member #14 (green). **b** Same as **a** but for the Atlantic Multidecadal Variability index (see Methods). **c** Same as **a** but for the Ocean Heat Content (OHC) anomaly (10⁹ J m⁻²) between 0 and 700 m averaged over the Newfoundland region minus the OHC averaged over the North Atlantic Subpolar Gyre region from the IAP observational dataset⁵⁵,⁵⁶ (blue). **d** Same as **a** but for the interhemispheric near-surface air temperature difference (K), with the HadCRUT4-CW observational dataset⁴¹,⁴² (blue). A Lanczos low-pass filter with a cutoff period of 11 years is used in **c** and **d** and the anomalies are computed over the displayed period.

rest of the ensemble (Supplementary Fig. 5). Over the 1940-2018 period, member #14 has an averaged SST trend of 0.06 K decade⁻¹ over the North Atlantic, slightly weaker than the observed trend of 0.065 K decade⁻¹ (Supplementary Fig. 5). The pattern is overall close to the observations, with a spatial correlation coefficient of 0.73. The subset of members previously identified are also among the members with average trends over the North Atlantic closest

to the observations. Due to its strong AMOC weakening (Fig. 2a), member #25 shows the smallest SST trend in the North Atlantic with a value of 0.042 K decade⁻¹.

**Further lines of evidence from AMOC fingerprints.** In order to strengthen the hypothesis that a weakening of the AMOC in the
real world since the middle of the 20th century may have masked a fraction of the anthropogenic global warming, we evaluate four different observable AMOC fingerprints against our ensemble members. A first SST fingerprint is calculated following the Caesar index¹⁷. Then, we use the Atlantic Multidecadal Variability²⁴⁻²⁶ index (AMV index, see Methods), as basin-wide low-frequency variations in the North Atlantic SST are related to the AMOC low-frequency variability in many climate models²⁷. Whether low-frequency variations in AMOC and AMV are the results of naturally-occurring internal variability, external forcings, or a mix of both is still debated²⁷⁻²⁹. This is also the case in the IPSL-EHS, with significant Pearson correlations ranging from 0.56-0.97 between the AMOC and the AMV index over the 1850-2018 period in the various members (see Methods). We calculate a third index based on the upper (0-700 m) Ocean Heat Content (OHC) averaged offshore Newfoundland minus the average over the North Atlantic Subpolar Gyre (ΔOHC). The upper OHC is indeed significantly correlated with the AMOC in IPSL-CM6A-LR over these two regions, with a negative correlation offshore Newfoundland and a positive correlation over the North Atlantic Subpolar Gyre (Supplementary Fig. 6a). This correlation dipole is consistent with observations and other climate model simulations³⁰. As no significant correlations are found between the upper OHC and the AMOC forced signal, these correlations are mainly due to the AMOC internal variability (Supplementary Fig. 6b). Finally, the fourth index is defined as the difference between the Northern and the Southern Hemisphere anomalies of near-surface air temperature (ΔITA), which might also be a good indicator of the AMOC low-frequency variability, given its impact on interhemispheric temperature³¹. The time series of these four observation-based indices are included within the range of IPSL-EHS (Fig. 3), which gives confidence in the capacity of the model to represent climate dynamics accurately.

Furthermore, a good consistency is found between the six members identified above as best fitting GSAT over the observed period and the observations for the four AMOC fingerprints. In member #14 for example, Pearson correlations amount to 0.48 for the Caesar index (Fig. 3a), 0.49 for the AMV index (Fig. 3b), 0.46 for the ΔOHC index (Fig. 3c) and 0.74 for the ΔITA index (Fig. 3d) and are significant for the Caesar index, the ΔOHC index and ΔITA index (see Method). Notwithstanding, the decrease of −0.13 K decade⁻¹ of the Caesar index in member #14 since the 1940s, consistent with a weakening AMOC, is close to the observed trend of −0.11 K decade⁻¹. The decrease in the observed ΔOHC index of −0.45 10⁹ J m⁻² decade⁻¹ is also close to observations in member #14, with a decrease of about −0.42 × 10⁹ J m⁻² decade⁻¹. The pattern of the observed North Atlantic upper OHC trend calculated over the 1940-2020 period is consistent with a weakening AMOC, with significant positive trends over Newfoundland and significant negative trends over the Subpolar Gyre (Supplementary Fig. 6c). This pattern is no longer visible considering only the last decades (1980-2020), a period more influenced by the external forcing and characterized by strong warming of the upper OHC over the polar region (Supplementary Fig. 6d). The fact that these six members whose AMV is overall in phase with observations are those with a GSAT evolution most consistent with observations underscores the important role of the AMV in the GSAT variability at least in this model. This is consistent with a recent study³² showing the main role of the AMV in the internal variability of GSAT on multi-decadal time scales. This good consistency is reinforced by the fact that member #14 has low RMSEs for these four AMOC fingerprints.

Member #14 is not the only one with overall low RMSEs for these four AMOC fingerprints. Members #29, #4, #5, #25, and #30, identified previously, also have overall the lowest RMSEs in

comparison to the rest of the ensemble (Supplementary Fig. 6). Notwithstanding, the consistency between AMOC fingerprints in the subset of members identified in Fig. 2a and in observations reinforces the hypothesis that the AMOC may have been weakening since the second half of the 20th century, thereby masking a fraction of the anthropogenic global warming as seen in the IPSL-EHS.

**Multi-model analysis.** We now examine how other climate models support our analysis of IPSL-EHS. A similar significant positive relationship between the AMOC and the GSAT 60-year trends is found in the piControl simulation of IPSL-CM6A-LR, as well as in ten out of thirteen other CMIP6 models (IPSL-CM6A-LR, both CNRM models, MPI-ESM1-2, CESM2-WACCM, SAM0-UNICON, CMCC-CM2, MRI-ESM2, EC-Earth3, and CIESM) (Fig. 4a). Those models also generally show significant multi-centennial AMOC variability compared to red noise at periods longer than 100 years (Fig. 4b). The models with the strongest AMOC-GSAT trends relationships are also those with the largest GSAT variance in their piControl simulations. Taken together these results highlight the role of internal variability in these models and suggest that the relationship between the AMOC and GSAT trends found in IPSL-CM6A-LR is rather widespread in CMIP6 models albeit with varying strength and range.

Only a few models with large historical ensembles are available to repeat the analysis performed with the IPSL-EHS: MPI-ESM1.1 with a low ECS, CNRM-CM6A-1, and CanESM5 with a high ECS. We investigate the consistency of our previous results with these models in Supplementary Section 2, which is summarized here. First, one can see that some members in the CNRM and the MPI ensembles reproduce the observed warming (similar GSAT trend as observations), unlike the CanESM5 ensemble. This latter model is therefore omitted. This may be due to its large climate sensitivity and/or its weak internal climate variability. Furthermore, a similar positive relationship between AMOC and GSAT trends is found in the CNRM historical simulations as in the IPSL-EHS, with, nevertheless, larger regression coefficients due to the larger internal climate variability of CNRM-CM6A-1¹¹. A positive relationship is also found in the MPI ensemble, with a surprisingly large fraction of the historical members showing larger GSAT trends than observed despite a low climate sensitivity for this model. This might be due to its weak aerosol forcing. In addition, the members which fit best the GSAT and the AMOC weakening also match best the AMOC fingerprints in the CNRM-CM6 ensemble (Supplementary Fig. 9). This is not the case for the MPI ensemble where these members are scattered throughout the other members of the ensemble for the three fingerprints, which means that the members of the MPI ensemble with a weakening AMOC do not project well onto the observed AMOC fingerprints (Supplementary Fig. 8).

In conclusion, the analysis of three other large ensembles of historical simulations does not invalidate our results that suggest that the weakening in AMOC reconstructions since the second half of the 20th century might be due to internal variability. Indeed, similar behaviour can be seen in one model (CNRM-CM6 with a high ECS), the CanESM5 model does not reproduce the IPSL-CM6A-LR features as it appears to have warming not compatible with the observations, and the MPI-ESM1.1 model, in which the members with an AMOC trend consistent with a recent reconstruction¹⁷ are not those with the best representation of the four observed AMOC fingerprints analyzed in this study.

**Implications for future global warming and discussion.** Our analysis suggests that a fraction of anthropogenic warming might
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Relationship between GSAT and AMOC in CMIP6 control simulations. a Scatter plot of Atlantic Meridional Overturning Circulation (AMOC) (Sv decade⁻¹) versus Global near-Surface Air Temperature (GSAT) (K decade⁻¹) trends calculated overall 60 year windows increasing by 10 years from the first 500 years of the pre-industrial control (piControl) simulations of thirteen CMIP6 models. The lines represent the least square regression between these two variables for each model. The regression coefficients (K Sv⁻¹) are indicated in bold beside the model's name when the regression is significant (p < 0.1, see Methods). b Smoothed power spectra of the AMOC time series from the first 500 years of the piControl simulations of the same thirteen CMIP6 models (solid lines, thick when the regression is significant), with the 95% confidence limit estimated from a red noise (dashed lines). A Lanczos low-pass filter with a cutoff period of 5 years is used.</figcaption></figure>

have been hidden by an AMOC weakening that started in the middle of the 20th century and is mainly related to internal climate variability. Indeed, the subset of members identified previously, which is characterized by a strong AMOC weakening over the 1951–1990 period in comparison to the IPSL-EHS (Fig. 5a, b) shows a lower warming over the same period relative to the full IPSL-EHS (Fig. 5c, d). The same subset of members experiences an internally-generated AMOC strengthening over the next few decades. This is associated with a larger warming rate of about 0.36 K per decade relative to the IPSL-EHS mean of 0.34 K per decade. This influence of a phase change in the AMOC is however limited, as both the AMOC and GSAT low-frequency internal variability are projected to decrease in response to external forcings³³. Taken together, these results reinforce the risk of crossing the 2 °C warming objective. Indeed, our subset of members with a low-frequency variability consistent with the observations shows larger warming than the ensemble means over the next decades. Moreover, taking into account this AMOC weakening over the historical era reinforces the risk of faster warming, as a fraction of anthropogenic global warming, which could have been hidden by this low-frequency internal variability, is expected to materialize in the coming decades. These results thus seem to be in line with a recent study³⁴ suggesting that the pattern effect could be in part related to internal climate variability. This internally-driven pattern effect could have masked part of human-induced global warming in recent decades. Nevertheless, the time period investigated here is longer than the last four decades analyzed in this other study, and our mechanism is related to the North Atlantic variability, rather than the Pacific variability, as found in the previous studies³⁵. Future work is required to quantify how accounting for North Atlantic low-frequency variability, both in the Atlantic and Pacific oceans, modifies observational constraints on future warming levels.

The realism of the multi-centennial low-frequency internal variability found in some of the CMIP6 models¹⁰,¹¹ is a crucial element of our results. Although the instrumental period is short,

previous studies suggested that the Atlantic Multidecadal Variability is underestimated in CMIP5 models, with a persistence lower than the one deduced from observations³⁶,³⁷. To gain insight on this realism, paleodata provides additional information, but it usually suffers from large uncertainties concerning quantified estimations of the variability, especially at large spatial scales. While it is usually believed that model simulations might have too low multi-centennial variability as compared to proxy records³⁸, a recent study¹¹ suggests that GSAT interdecadal variability of some CMIP6 models might be overestimated over the period 1450–1840 in comparison to the pre-industrial control simulations. Therefore, we cannot exclude that some CMIP6 models, such as IPSL-CM6A-LR have too much internal variability. On the longer time scale of the last two thousand years, the different methods to reconstruct global mean surface temperature do show considerable uncertainty in terms of the magnitude of the multi-centennial variability³⁹, while the external forcings and their impact remain poorly estimated in models. Therefore, it seems difficult at the moment to properly assess the realism of model simulations with those reconstructions. Such an evaluation of model simulations might deserve a dedicated analysis using the PAGES2K database and last millennium simulations from CMIP6, using advanced techniques like pseudo-proxy approaches⁴⁰ in order to compare model simulations and reconstructions in a coherent framework.

To conclude, we have shown here that the different indications of an AMOC weakening since the mid-20th century might be mainly of internal origin coming from multi-centennial variability of the ocean. If true, this might mean that transient climate sensitivity estimated from the observational records, especially over the last 6–7 decades may be underestimated. Thus, emergent constraint approaches that try to constrain future warming using the recent decades should fully embrace the issue of low-frequency internal variability and take into account individual ensemble members rather than ensemble means, as this might have crucial implications in terms of how different models are weighted in such studies.
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Implication for near-future AMOC and GSAT change. **a** Time evolution of the low-pass filtered Atlantic Meridional Overturning Circulation (AMOC) anomaly (Sv) relative to the 1900–2018 period for the IPSL ensemble of extended historical simulations (IPSL-EHS) (gray), the ensemble mean (red), the subset of members identified in Fig. 2 (member #14, 5, 25, 29, 4 and 30 in light blue), the subset mean (dark blue) and the mean internal variability (calculated by removing the ensemble mean for each member) of the subset of members (purple). **c** Same as **a** but for the Global near-Surface Air Temperature (GSAT) anomaly (K). **b** AMOC trends (Sv per decade) in IPSL-EHS (black), with the ensemble mean (black dot) and the full minimum-maximum range (black line) and the subset of members (light blue), with the subset mean (dark blue). **d** Same as **c** but for GSAT.</figcaption></figure>

## Methods

**Observational data.** We used the infilled Cowtan and Way$^{41,42}$, named HadCRUT4-CW in this study, observational surface temperature, which avoids the problem of missing values. This dataset is composed of blending sea surface temperature over the ocean and near-surface air temperature over the land. This metric of global mean temperature, however, warms significantly less than the GSAT. To take into account this discrepancy, we applied a factor of 1.06 to the HadCRUT4-CW surface temperature, as used in a recent study$^{10}$. This factor is similar to those of 1.05$^{43}$, 1.06$^{44}$, and 1.08$^{5}$ used in several other studies. One of the interests of this dataset is to provide an ensemble of realizations, allowing to consider the associated uncertainties. Here, the median of the dataset is considered the best estimate, and the maximum and the minimum is used for the uncertainty. Note that the HadCRUT4-CW dataset is very consistent with the Berkeley$^{45}$ observational surface air temperature, with a Pearson correlation coefficient between the annual GSAT values of 0.99. Very small differences are found with the infilled HadCRUT5 dataset when accounting for the missing values of this dataset. Therefore, we keep the HadCRUT4-CW in order to have a global coverage without missing values.

**Simulations used.** The atmospheric resolution of IPSL-CM6A-LR is 1.26° × 1.25° with 79 levels (model top at 1 Pa). The model uses a nominal resolution of 1° in the ocean and 75 levels. The IPSL-CM6A-LR ensemble$^{12}$ used in this study follows the CMIP6 protocol$^{46}$ for historical simulations for the period 1850–2014. The historical simulations initial conditions are taken from different years in a long pre-industrial simulation after it has reached a quasi-stationary state. Specifically, the simulations are started from the atmospheric, oceanic, and land surface initial conditions of the 1st January from different years. The simulations were extended until 2060 using all forcings from the SSP245 scenario$^{47}$, except for the ozone field which has been kept constant to its 2014 climatology (as this particular forcing was not available at the time of performing the extensions). This implies that these simulations do not “see” the ozone hole recovery and changes in tropospheric

ozone, in contrast to the official CMIP6 projections. However, this is seen as a minor shortcoming for the purpose of this study. All the simulations used are summarized in Supplementary Table 1.

**ECS and TCR estimates.** The Sensitivity and the Transient climate response calculated over the historical period, named as S\_hist and TCR\_hist in this study, are defined as:

$$S\_hist = \frac{F_{2*CO_2} * \Delta T}{\Delta F - \Delta Q} \qquad (1)$$

$$TCR\_hist = \frac{F_{2*CO_2} * \Delta T}{\Delta F} \qquad (2)$$

with $F_{2*CO_2} = 3.77 \text{ Wm}^{-2}$ being the radiative forcing due to a doubling of the atmospheric CO$_2$ concentration. This latter quantity is calculated by averaging the Effective Radiative Forcing (ERF) values from three different experiments performed with IPSL-CM6A-LR, a fixed SST experiment$^{48}$ and two experiments based on this method$^{49}$, using linear regressions of the top-of-atmosphere net radiative flux against surface temperature carried out over the first 20 years of the abrupt experiment, with a subtraction of the corresponding piControl quantities for both, but based on a year-to-year basis for the first case or using 20-year climatologies for the second case. We estimate the uncertainty related to using the lowest (3.50 W m$^{-2}$) and largest (3.94 W m$^{-2}$) values of the three experiments. The net ERF is calculated between the pre-industrial (1850–1879) and present-day (1999-2018) periods by averaging the ERFs of the three members of piClim-histall RFMIP experiments and is estimated at $\Delta F = 2.09 \text{Wm}^{-2}$. $\Delta T$ is the change in GSAT between the same periods, and $\Delta Q$ the rate of the total increase in Earth system heat content as diagnosed in the model, calculated over the 1999-2018 period. Here, we approximate the Earth System Heat Content by calculating $\Delta Q$ from the Ocean Heat Content increase over the historical period divided by 0.89, as the oceans store about 89% of the excess heat$^{50}$. This allows us to take into account an
30. Zhang, R. Coherent surface-subsurface fingerprint of the Atlantic meridional overturning circulation. *Geophys. Res. Lett.* **35**, 035463 (2008).

31. Marshall, J., Donohoe, A., Ferreira, D. & McGee, D. The ocean's role in setting the mean position of the Inter-Tropical Convergence Zone. *Clim. Dyn.* **42**, 1967–1979 (2014).

32. Li, Z. et al. A robust relationship between multidecadal global warming rate variations and the Atlantic Multidecadal Variability. *Clim. Dyn.* **55**, 1945–1959 (2020).

33. Cheng, J. et al. Reduced interdecadal variability of Atlantic meridional overturning circulation under global warming. *Proc. Natl Acad. Sci.* **113**, 3175–3178 (2016).

34. Zhou, C., Zelinka, M. D., Dessler, A. E. & Wang, M. Greater committed warming after accounting for the pattern effect. *Nat. Clim. Change* **11**, 132–136 (2021).

35. Dong, Y., Proistosescu, C., Armour, K. C. & Battisti, D. S. Attributing historical and future evolution of radiative feedbacks to regional warming patterns using a green's function approach: the preeminence of the western pacific. *J. Clim.* **32**, 5471–5491 (2019).

36. Cheung, A. H. et al. Comparison of low-frequency internal climate variability in CMIP5 models and observations. *J. Clim.* **30**, 4763–4776 (2017).

37. Qasmi, S., Cassou, C. & Boé, J. Teleconnection between atlantic multidecadal variability and european temperature: diversity and evaluation of the coupled model intercomparison project phase 5 models. *Geophys. Res. Lett.* **44**, 11,140–11,149 (2017).

38. Laepple, T. & Huybers, P. Ocean surface temperature variability: Large model-data differences at decadal and longer periods. *Proc. Natl Acad. Sci.* **111**, 16682–16687 (2014).

39. Neukom, R. et al. Consistent multi-decadal variability in global temperature reconstructions and simulations over the Common Era. *Nat. Geosci.* **12**, 643–649 (2019).

40. Ortega, P. et al. A model-tested North Atlantic Oscillation reconstruction for the past millennium. *Nature* **523**, 71–74 (2015).

41. Jones, P. D. et al. Hemispheric and large-scale land-surface air temperature variations: an extensive revision and an update to 2010. *J. Geophys. Res. Atmospheres* **117** (2012).

42. Cowtan, K. & Way, R. G. Coverage bias in the HadCRUT4 temperature series and its impact on recent temperature trends. *Q. J. R. Meteorol. Soc.* **140**, 1935–1944 (2014).

43. Cowtan, K. et al. Robust comparison of climate models with observations using blended land air and ocean sea surface temperatures. *Geophys. Res. Lett.* **42**, 6526–6534 (2015).

44. Gillett, N. P. et al. Constraining human contributions to observed warming since the pre-industrial period. *Nat. Clim. Change* **11**, 207–212 (2021).

45. Rohde, R., Muller, R., Jacobsen, R., Perlmutter, S. & Mosher, S. Berkeley Earth temperature averaging process. *Geoinformatics Geostat. Overv.* **01**, 1000103 (2013).

46. Eyring, V. et al. Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization. *Geosci. Model Dev.* **9**, 1937–1958 (2016).

47. Gidden, M. J. et al. Global emissions pathways under different socioeconomic scenarios for use in CMIP6: a dataset of harmonized emissions trajectories through the end of the century. *Geosci. Model Dev.* **12**, 1443–1475 (2019).

48. Forster, P. M. et al. Recommendations for diagnosing effective radiative forcing from climate models for CMIP6. *J. Geophys. Res. Atmospheres* **121**, 12,460–12,475 (2016).

49. Gregory, J. M. et al. A new method for diagnosing radiative forcing and climate sensitivity. *Geophys. Res. Lett.* **31**, 01847 (2004).

50. von Schuckmann, K. et al. Heat stored in the Earth system: where does the energy go? *Earth Syst. Sci. Data* **12**, 2013–2041 (2020).

51. Hobbs, W., Palmer, M. D. & Monselesan, D. An energy conservation analysis of ocean drift in the CMIP5 global coupled models. *J. Clim.* **29**, 1639–1653 (2016).

52. Trenberth, K. E. & Shea, D. J. Atlantic hurricanes and natural variability in 2005. *Geophys. Res. Lett.* **33**, 026894 (2006).

53. Ebisuzaki, W. A method to estimate the statistical significance of a correlation when the data are serially correlated. *J. Clim.* **10**, 2147–2153 (1997).

54. Huang, B. et al. Extended reconstructed sea surface temperature, version 5 (ERSSTv5): upgrades, validations, and intercomparisons. *J. Clim.* **30**, 8179–8205 (2017).

55. Cheng, L. & Zhu, J. Benefits of CMIP5 multimodel ensemble in reconstructing historical ocean subsurface temperature Variations. *J. Clim.* **29**, 5393–5416 (2016).

56. Cheng, L. et al. Improved estimates of ocean heat content from 1960 to 2015. *Sci. Adv.* **3**, e1601545 (2017).

## Acknowledgements

The IPSL-CM6 experiments were performed using the HPC resources of TGCC under the allocations 2017-R0040110492, 2018-R0040110492, and 2019-A0060107732 (project gencmip6) provided by GENCI (Grand Equipement National de Calcul Intensif). The IPSL-CM6 team is acknowledged for data curation of the CMIP6 model output. RB and OB were supported by the European Union's Horizon 2020 research and innovation programme under grant agreement number 820829 for the "Constraining uncertainty of multi-decadal climate projections (CONSTRAIN)" project (07/2019-06/2024). JM, GG, and DS were supported by Blue-Action (European Union's Horizon 2020 research and innovation program under grant agreement no 727852) and EUCP (European Union's Horizon 2020 research and innovation programme under grant agreement no 776613) projects. This work also benefited from the French state aid managed by the ANR under the "Investissements d'avenir" programme with the reference ANR-11-IDEX-0004-17-EURE-0006.

## Author contributions

R.B., D.S., O.B., and J.M. designed the study. R.B. performed all of the IPSL-EHS and the multi-model analysis against observations. G.G. computed the A.M.O.C. index based on the O.H.C. AS provided the radiative forcing estimates. RB wrote the paper, and O.B., D.S., J.M., J.D., G.G., F.H., and J.S. contributed to the text.

## Competing interests

The authors declare no competing interest.

## Additional information

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-021-26370-0.

Correspondence and requests for materials should be addressed to Rémy Bonnet.

Peer review information *Nature Communications* thanks Nathan Gillett, Aixue Hu and the other anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021