The AMOC is a large-scale oceanic circulation that brings warm and salty water in the North Atlantic high latitudes, where it cools, releases heat to the atmosphere, and eventually sinks to the deeper ocean after a number of complex oceanographic processes. Its impacts on climate and biogeochemistry is worldwide and complex as well, driven notably by large-scale atmospheric teleconnection patterns.

The classical view of the AMOC as a global conveyor belt from Broecker (1987) has been deconstructed by physical oceanography literature (Lozier, 2010). Indeed, by construction, the AMOC is aggregating a mosaic of currents with different dynamics and drivers and therefore multiple processes that can affect them. A classical belief remains that the AMOC is fed by the North Atlantic convection sites, where deep water is formed. Those sites, mainly located in the Nordic Seas and in the subpolar gyre (SPG), are sensitive to a lot of forcing, including stochastic atmospheric variations and interactions with sea ice, thus their exact location can move in a changing climate (cf. Lique and Thomas, 2018). The paradigm of the importance of deep convection on the AMOC has been recently questioned, using less than 5 years of oceanic observations (Li et al., 2021a). Water mass transformation in the deeper ocean, through eddies (Lozier et al., 2019), internal tides (de Lavergne et al., 2017), as well as the wind forcing, clearly also plays an important role in the AMOC behavior. Thus, there are still a number of issues regarding the drivers of its energetics, its dynamics, and its latitudinal coherence. In all those aspects, observational data at different locations in the Atlantic Ocean and on a long-enough time frame remain the key for our progress in understanding.

As it is clear from the experimental design, the ocean-only HR and coupled LR twin simulations (melting vs. control) are different and cannot be rigorously directly compared. The AOGCM LR simulations comprise a large set of members and are integrated over a long time-frame, while the OGCM HR simulations are forced simulations integrated over a shorter period and include both runoff and solid ice discharge at coastal locations (no iceberg transport), both considered as freshwater input. We here explore and illustrate the main responses of the ocean to the additional freshwater in the North Atlantic in both pairs of simulations. Figures 2, 3 show the changes in SSS in the HR OGCM and LR AOGCM, respectively, estimated by comparing control and sensitivity simulations. The two sets of simulations have several features in common in terms of SSS anomalies. The strongest SSS anomalies averaged over the respective simulation duration can be found in the Baffin Bay where they reach more than 1 PSU in HR. Large anomalies are also found around the Labrador basin, along the West Greenland Current and Labrador Current pathways. The fresh anomalies then mix across the North Atlantic Current (NAC) and spread eastward across the North Atlantic Ocean along the southern edge of the subpolar gyre (Supplementary Video 1).

There are however significant differences between the two sets of simulation. In HR, the largest anomalies remain along the Greenland coast, near the freshwater source while, in LR, the shelf-basin contrasts around Greenland are much less preserved. The year-to-year development of those SSS anomalies in HR (Supplementary Video 1) actually shows that the anomalies are first found at their injection points around Greenland. After about 2 years (in 2006), we see a clear accumulation in the Baffin Bay and to a lesser extent along the Labrador Current. After around 4 years (in 2008), the SSS anomalies already invade the Labrador Sea interior and start to develop along the NAC as well. At the end of the simulation (from about 2015) those SSS anomalies reach the eastern North Atlantic and start to go northward or southward. A positive SSS anomaly also starts to develop in the Irminger Sea, most notably at the end of the simulation. The surface freshwater pathways are different in LR with the freshwater continuing along the Canary Current, forming the so-called “freshwater leakage” highlighted in Swingedouw et al. (2013) as an escape pathway for freshwater to leave the SPG. This signature is less clear in the HR simulations and might be related to the longer temporal size of the LR simulations, allowing for the freshwater leakage to develop more.

In the center of the SPG, there is a slight positive anomaly, which might reflect a dynamical adjustment from the ocean in the LR simulations (Figure 3). We hypothesize that this anomaly could be due to strengthening of the SPG due to the increase in center to shelf-basin density gradient (e.g., Born and Stocker, 2014), related to the freshening of the SPG border. This strengthening can increase the entrainment of NAC water masses in the SPG, and might explain the slight increase of SSS found in the Irminger Sea in the LR simulations. This positive SSS anomaly is also clearly reminiscent of the positive anomaly found in the Irminger Sea in the HR simulation (in all years, although with varying amplitude, Supplementary Video 1).

A crucial difference between the two sets of simulations can be found in the central Labrador Sea, where the anomaly is larger in the HR. It is 0.31 PSU when averaged over the period 2004–2016 and over the convection site defined in Figure 2. It represents almost 1.5 interannual standard deviation in the HR control simulation. It even reaches more than 0.5 PSU in 2016. In the LR runs, the difference is 0.06 PSU when averaged over the whole period in the Labrador Sea convection site defined in Figure 3, or about 20% of the interannual standard deviation of the control runs. The penetration of the freshwater into the Labrador Sea interior, where convection occurs in the control runs, seems therefore more efficient in the HR runs (about 5 times more at least). Although the differences in freshwater input between HR and LR prevent any strong conclusions on this difference, we argue that part of it might be related to the differences in resolution as will be further discussed in the conclusions of this paper.

Figure 4 shows the mixed layer depth (MLD) averaged over the 2004–2016 period in the North Atlantic Ocean in the HR simulations. The MLD distribution in the control simulation (Figure 4, upper panel) captures the convection sites in the Labrador Sea, southeast and west of Iceland, in the Irminger Current, and in the Nordic Seas. The site with deepest MLD over the period considered is the Labrador Sea. This picture is in qualitative agreement with observations (e.g., Yeager et al., 2021), but we notice that convective activity in the Nordic Seas is underestimated, probably as a result of too fresh upper waters, and consequently too strong background stratification, which arises in absence of any surface salinity restoring in these simulations. The convection in the Irminger Sea seems also a bit smaller than the one estimated from Argo floats in Rühs et al. (2021). The distribution of MLD anomalies in the lower panel of Figure 4 shows that the Labrador Sea is the primary impacted region by additional freshwater input, while the other convections sites show limited changes in MLD. In the Labrador Sea, the anomalies over the whole period exceed 300 m, which is about half of the mean state value.

The time evolution of the MLD in the Labrador Sea, shown in Figure 5, allows to see that the convection in the Labrador Sea has a strong interannual variability and can be compared with 2004–2014 observation-based estimate of convective activity from Figure 3 of Kieke and Yashayaev (2015), with some years with almost no convection activity (e.g., 2010) while others exhibit strong convective events (e.g., 2008). The strong convective event found in 2015 is also in agreement with observations (Yashayaev and Loder, 2017). A first impact of the GrIS melting on the mean MLD in the Labrador Sea can be found after about 2 years, from which the sensitivity experiment persistently shows smaller MLD than the control one. Nevertheless, it is after 10 years that the impact becomes really large, with clear differences higher than 3 standard deviations of interannual variability, notably in 2015 and 2016, where the high convective activity simulated in the control run is strongly reduced in the sensitivity experiment. Nevertheless, when standardized with the mean MLD from the control simulation, the changes seem to develop more monotonically as soon as after 1 year of freshwater release. Thus, we can conclude that the GrIS melting and Arctic runoff have a considerable impact on the Labrador Sea convective activity after less than a decade of additional freshwater input in those HR OGCM simulations. We argue, in agreement with Supplementary Video 1, that most of the freshening and therefore of the impact on MLD is coming from the Greenland vicinity, since the forcing in the Arctic has lower amplitude and slowly spreads in the subpolar region. A dedicated experiment with only GrIS melting from Bamber et al. (2018) might be necessary to fully prove this.

We now compare these results of MLD in the HR with the simulations from the LR. Figure 6 shows the differences in MLD over the 95 years of simulations averaged on the two 20-member ensembles available. Similarly to the results from Devilliers et al. (2021), we see a slight reduction of MLD in the Labrador Sea, while it is enhanced in the Irminger Sea. The changes remain however limited, about an order of magnitude smaller than the mean state MLD. The largest differences in MLD can be found in the Nordic Seas, although they also remain almost one order of magnitude lower than the mean state.

The fact that the amplitude of the MLD changes remains limited in the LR simulations is further supported by Figure 7, which shows the time evolution of the MLD in the Labrador Sea (defined in Figure 3). From this figure, it is clear that the differences are very small and not significant in view of the ensemble spread. We thus conclude, as in Devilliers et al. (2021), that in this CMIP6 OAGCM, the inclusion of observed GrIS melting has a very limited impact on deep convection activity over the last century.

We now look at the impact of freshwater input most notably from GrIS melting on the AMOC. In the LR model, as shown in Devilliers et al. (2021), the impact of GrIS melting is very limited. Here, with our larger ensemble of 20 members, we do find at 45°N a slight decrease of the AMOC of 0.20 ± 0.39 Sv (0.18 ± 0.33 Sv at 26°N) over the whole period. It is 0.13 ± 1.0 Sv (0.1 ± 0.8 Sv at 26°N) over the period 2004–2014. We also compute the AMOC changes in density space (Supplementary Table 1). For the latter period, we also find a slight weakening of 0.24 ± 1.31 Sv at 45°N. When computed in sigma0-coordinate over the OSNAP section (Lozier et al., 2017) for the period 2004–2014, we find a MOC weakening of 0.57 ± 1.68 Sv along OSNAP east and a slight enhancement of 0.22 ± 0.84 Sv along OSNAP west, none of which being significant at the 95% level.

In the HR control simulation, the AMOC over the OSNAP array (Figure 8) shows similar vertical distribution and maxima as in the observations in their overlapping period (Lozier et al., 2019). Over the OSNAP array, the mean AMOC state (Figure 8) is in broad agreement with observational estimates. The AMOC weakening in the HR simulations (Figure 9) notably due to GrIS melting is about one order of magnitude higher than in the LR, a result in agreement with the larger changes in SSS and MLD in the Labrador Sea. The differences between melting and control simulations are always negative after 3 years of integration (2007). It becomes larger than 1 Sv after 10 years of integration (2014) both at OSNAP and at 45°N. It reaches a value larger than 1.5 Sv at the end of the simulation at 45°N and above 2 Sv at OSNAP array. The response of the AMOC at 45°N exhibits a negative difference emerging clearly from around 2012–2014, or after 8–10 years of integration of the meltwater anomalies. Such a time scale is short, but remains consistent with the 5–6 years estimated in Desbruyères et al. (2019) between changes in water mass transformation and the AMOC at 45°N from observation-based indices. We also notice a slight change in the distribution of the MOC transport in density space, with a shift toward lighter density range (Figure 8), which is also in agreement with what is found in the LR simulations (Supplementary Figure 3).

In this paper, we have assessed the recent literature since the IPCC-SROCC report (Collins et al., 2019) regarding changes in the AMOC over the last few decades and the future ones. Furthermore, we have presented results from a few new simulations that focus on a key element of AMOC changes, still poorly evaluated, namely the impact of Greenland Ice Sheet (GrIS) melting, by using new estimates from Bamber et al. (2018) that have been hardly introduced in model simulations up to now.

While the uncertainty in past and future AMOC changes remains very large, as discussed in the new IPCC AR6 report in 2021 (IPCC, 2021), we highlight here that sub-grid scale processes might considerably affect attribution of possible on-going AMOC weakening in the IPCC models. Indeed, HR models (<1/10°), which are not available in CMIP6-type models, do show considerable differences with the latter. First, the mean state is strongly improved with resolution, and the AMOC is usually far more vigorous in eddy rich models (Hirschi et al., 2020), which makes it more difficult to evaluate the amplitude of potential AMOC weakening and associated heat transport in the future. Secondly, the GrIS melting is still poorly accounted for in CMIP6 models. Our new results from HR simulations suggest that this melting might already have had a strong impact in the recent past on the Labrador Sea convection (as also found in Böning et al., 2016) and the AMOC, contrary to what is simulated in one CMIP6-type model. This further questions our ability to estimate future AMOC changes and therefore leads to considerable uncertainty worldwide in the future climate projections (Bellomo et al., 2021).

A number of limitations remain concerning the new results presented here and will need to be solved in future studies. The principal issue is the difficulty to compare LR and HR simulations. Here this difficulty is amplified by the fact that freshwater input differs between the two models. The anomaly of freshwater release from Greenland is about 35 mSv in the HR runs and 20 mSv in the LR over the overlapping period (2005–2014), but the HR has additional coastal runoff in the Arctic that clearly blurs the attribution possibilities. Still, we argue here based on Supplementary Video 1 that those additional fluxes elsewhere than in Greenland vicinity play a limited role in the results presented. A simulation accounting only for Greenland ice sheet perturbation from Bamber et al. (2018) would be also welcome to entirely prove that Arctic runoff plays a limited role in the results shown.

Furthermore, the fact that the response in the interior Labrador Sea seems stronger in the HR simulations than in the LR can also be limited by the fact that we compare a single OGCM realization with an ensemble mean of AOGM simulations, which might better identify the signal out of any potential noise effect. This could be justified by the fact that OGCM might exhibit lower intrinsic variability, since the atmospheric forcing is fixed. Nevertheless, the oceanic intrinsic variability is far from being negligible at HR resolution (Leroux et al., 2018) so that an ensemble might also be needed to properly assess the impact of additional freshwater input. We argue that the amplitude of the anomalies found for the MLD in the Labrador Sea is large enough to be beyond this uncertainty. Nevertheless, this remains to be fully assessed and proven.

Another candidate (on top of differences in the experimental design between LR and HR simulations) to explain the enhanced response in the Labrador Sea in HR is most probably related to the greater ability of the HR model to spread freshwater from the coast of Greenland, where the melting is injected, toward the interior basin, where most of the convection occurs. This is done through instability of the mean flow, here the West Greenland boundary current, and subsequent formation of eddies transporting shelf water into the ocean interior. The eddy kinetic energy, a measure of the eddy activity, is indeed very high off the southwestern coast of Greenland (Figure 10). Such a process has also been recently nicely illustrated in observational data, by following GrIS meltwater with gliders using fluorescing dissolved organic matter, a proxy of those meltwaters (Hendry et al., 2021). The nature of eddies in the Labrador Sea is however complex as illustrated in Rieck et al. (2019), with different types of eddies, generated either by baroclinic or barotropic instabilities and strongly dependent on bathymetry characteristics.

Furthermore, eddies might not be the only processes that are playing a role in high-resolution models. Indeed, Schulze Chretien and Frajka-Williams (2018) recently suggested that the freshwater transport from the shelves into the central Labrador Sea in the first top 30 meters might be largely due to wind-driven Ekman transport in the West Greenland Current, a process that is strongly dependent on the representation of the regional-scale energetic winds in the atmosphere (e.g., Spall et al., 2021), and is therefore poorly represented in LR models. Also, the finer scales resolved in the HR model allow to better represent the boundary currents, which are channeling most of the freshwater transport over the shelf around Greenland. In the LR model these currents are wider so that the freshwater release from the GrIS is far more diluted in the offshore direction. As a consequence, the fine scales resolved in the HR simulation increase the coastal freshwater input from eastern Greenland to the Labrador Sea shelves, which might also contribute to enhancing the freshwater transport toward the Labrador Sea center.

An important aspect of the HR simulations presented here is the fact that no SSS restoring at all is considered. Indeed, such a restoring in an ocean-only model might partially offset the impact of the meltwater in the sensitivity experiments. Thus, even though the absence of restoring strongly impacts the mean state (Biastoch et al., 2021) and might explain the loss of vigorous convective activity in the Nordic Seas in the HR control simulation, it is something necessary to correctly assess in an ocean-only model the impact of additional surface freshwater sources. In this respect, this constitutes an improvement as compared to the work of Böning et al. (2016), in addition to the update of GrIS meltwater reconstruction used in our study. Nevertheless, it is also clear that inclusion of Arctic freshwater sources in addition to the one of Bamber strongly limits the interpretation of our results in terms of attribution of GrIS freshwater impact.

The use of a regional ocean model in the North Atlantic also means that all ocean fields are prescribed at the boundaries from the parent global simulation, thus constraining the AMOC variability to that of the parent simulation. Since the prescribed fields are kept unchanged between the control and sensitivity simulations, the AMOC at the southern boundary at 28°N is the same in both runs in HR. This therefore limits the degrees of freedom of the AMOC in this model and its sensitivity close to this boundary. The differences we have computed are located further to the north, but the impact of the GrIS might be partly damped by the southern boundary conditions, meaning that the impact of GrIS melting on the AMOC could even be larger if those boundaries conditions were relaxed.

One aspect that might minor the implications of our results is the effective role of the Labrador Sea on the AMOC, a highly debated topic in the literature. Recent work from Lozier et al. (2019) based on the new OSNAP array clearly questioned its influence on AMOC variability and mean state. The key limit of those observations is their short time frame, which might strongly limit the robustness of the findings, as highlighted in Yeager et al. (2021) using 0.1° resolution coupled model simulations of 500 years duration. Nevertheless, it remains true that the Labrador Sea is contributing to a limited extent to the mean AMOC. It was estimated to be around 1/3 of AMOC strength in the synthesis of Quadfasel and Käse (2007), a number disputed by the paper of Lozier et al. (2019) that found only a few Sverdrups of overturning due to the Labrador Sea, i.e., only slightly more than 10% of the whole AMOC strength.

Another implication of our results concerns the risk of an SPG collapse, as discussed in recent literature (e.g., Born et al., 2013; Sgubin et al., 2017) and confirmed in some CMIP6 projections (Swingedouw et al., 2021). Those models do show a rapid cooling of the SPG of up to 2–3°C in less than a decade, associated with the cessation of convective activity in the whole gyre. The exact dynamics behind this non-linear behavior in some models and not in others remain unclear, but could be related to the misrepresentation of the advective salinity positive feedback, which is crucial for the SPG instability (Born et al., 2016). It has been argued that coarse resolution models might misrepresent processes related to this feedback (Born and Stocker, 2014), notably related to eddies. The results presented here clearly support this view, but the short time frame of the simulations and the fact that it is an ocean-only simulation prevent the possibility to estimate longer-term impacts on the SPG, and the possibility of a decadal-scale abrupt change, which might have crucial implications over Europe (Sgubin et al., 2019). To gain knowledge on this issue, HR climate models correctly initialized to present-day oceanic states might be necessary. However, we are still far from such tools, given the large computing cost of HR climate models as well as the issues related to initial shock in the decadal prediction system (Bilbao et al., 2021). In this respect new avenues should be explored, not only waiting for eddy rich climate models, which might necessitate a considerable amount of computing time, and will not be able to solve uncertainty related to low frequency internal dynamics (e.g., Bonnet et al., 2021) due to the low number of members that could be produced given the computing constraints.

While, at the moment, no parameterization in CMIP climate models represents the transfer of kinetic energy from small scale to large scale, the so-called energy backscatter, some work is being developed in ocean models to solve these limitations (Hewitt et al., 2020). Two main approaches are explored. The first one is introducing stochastic eddy parameterization (Grooms and Majda, 2013; Porta Mana and Zanna, 2014) in the models to account for the impact of those eddies. The other one is building on Gent and McWilliams (1990) parameterization and injects the potential energy extracted from GM directly into the momentum equation. Sub-mesoscale processes (below a few kilometers) can also be very important notably for air-sea fluxes (e.g., mesoscale and sub mesoscale winds, oceanic fronts). The parameterization of the closure of turbulent kinetic energy can also be explored using machine learning techniques, as recently done in Bolton and Zanna (2019) and Zanna and Bolton (2020). These promising approaches might help fill the gap of energy cascades in coarse resolution models and better represent some processes that appear to be essential given the results of our study.

To conclude, we can say that we are still far from a robust estimate of the fate of future AMOC weakening. There is already a very large uncertainty in coarse ocean models like CMIP6 ones (Bellomo et al., 2021). The present study further highlights the considerable role of missing processes in CMIP6 models for the response of the AMOC. The AMOC is a very complex feature that is still not well-understood both in its time scales, large-scale energetics and drivers and consequently its response to global warming is largely unknown. The results from HR model presented here seem to further amplify the potential uncertainty, which is already large in lower resolution models, by adding new sources from poorly represented processes.

Beyond the uncertainty on greenhouse gas emission scenarios, we would like here to emphasize that the causes of AMOC uncertainty are numerous. The first one is related to uncertainty in the atmospheric changes in climate projections, which are very model dependent. A second source of uncertainty is related to large-scale internal AMOC feedbacks (Swingedouw et al., 2007b), which are also not well-quantified nor constrained. The ocean-atmosphere coupling should not be forgotten in terms of uncertainty and processes, and further complexifies the picture (Bellucci et al., 2021). Last but not least, processes driving GrIS meltwater spread, deep convection or overflow processes are not correctly resolved at coarse resolution, and even at higher resolution, they remain difficult to be correctly represented (Colombo et al., 2020). We argue here that work on parameterization of those processes, based on their understanding in HR models, observations, as well as the use of new objective approaches from machine learning, is also highly necessary and will help the reduction of missing processes uncertainty.

On a shorter time frame of development, we argue here that emergent constraint approach (Hegerl et al., 2021) as used in, e.g., Weijer et al. (2020) should be further explored concerning the AMOC fate, by considering multiple driving factors of the AMOC spread that are effectively observed, in order to reduce the range of uncertainty from CMIP6 models. Finally, we support the development of several models worldwide and not only a limited number of a HR digital twin of the Earth as recently suggested (Bauer et al., 2021). Indeed, given the large number of uncertainty sources, the increase of resolution might not be enough to cover all of them, and might lead to serious underestimation of uncertainty in the case of a reduction of the number of models developed worldwide. The improvements of parametrization using HR models seem more promising and effective in that respect. A clear move toward this difficult task is therefore required.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

DS wrote most of the first draft, while all authors improved it. DS produced Figures 1, 3, 5–7. MD provided the GrIS melting forcing on the HR grid. CH, M-NH, and A-CB designed, performed, and analyzed the HR simulations. MD and DS performed the LR simulations. JD performed the OSNAP analysis from LR simulations. All authors contributed to the design and content of the article. All authors contributed to the article and approved the submitted version.

This work received support from the Blue-Action (European Union's Horizon 2020 research and innovation program, Grant No. 727852), EUCP (European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 776613) projects, and also benefited from funding from LEFE-IMAGO DECORATING project.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.