Mechanisms of internal Atlantic multidecadal variability in HadGEM3-GC3.1 at two different resolutions

: This study broadly characterizes and compares the key processes governing internal Atlantic multidecadal variability (AMV) in two resolutions of HadGEM3-GC3.1: N216ORCA025, corresponding to ∼ 60 km in the atmosphere and 0.25 8 in the ocean, and N96ORCA1 ( ∼ 135 km in the atmosphere and 1 8 in the ocean). Both models simulate AMV with a time scale of 60 – 80 years, which is related to low-frequency ocean and atmosphere circulation changes. In both models, ocean heat transport convergence dominates polar and subpolar AMV, whereas surface heat ﬂ uxes associated with cloud changes drive subtropical AMV. However, details of the ocean circulation changes differ between the models. In N216 subpolar subsurface density anomalies propagate into the subtropics along the western boundary, consistent with the more coherent circulation changes and widespread development of SST anomalies. In contrast, N96 subsurface density anomalies persist in the subpolar latitudes for longer, so circulation anomalies and the development of SST anomalies are more centered there. The drivers of subsurface density anomalies also differ between models. In N216, the NAO is the dominant driver, while upper-ocean salinity-controlled density anomalies that originate from the Arctic appear to be the dominant driver in N96. These results further highlight that internal AMV mechanisms are model dependent and motivate further work to better understand and constrain the differences.


Introduction
The observational record of North Atlantic sea surface temperature (SST) shows pronounced multidecadal variability, with the whole basin going through decades of warm and cool phases relative to global-mean temperatures (Sutton et al. 2018).This low-frequency variability has become known as Atlantic multidecadal variability (AMV).The AMV has been linked to a wide range of important impacts, such as Northern Hemisphere summertime conditions (Ruprich-Robert et al. 2018;Sutton and Hodson 2005), modulation of El Niño-Southern Oscillation (ENSO; Ruprich-Robert et al. 2017), strength of monsoons around the world (Monerie et al. 2019;Martin et al. 2014), Greenland ice melt (Hahn et al. 2018;Straneo and Heimbach 2013), and more.Therefore, because of the potentially large socioeconomic impacts of these connections, understanding the origin of AMV and the potential to predict it over the coming decades is of great interest.
Due to a dearth of observations, particularly in the deep ocean, studies have relied on coupled global circulation models (CGCMs) to understand the processes shaping AMV.CGCMs have consistently shown that AMV-like variability can arise spontaneously due to natural internal processes (Delworth et al. 1993;Griffies and Tziperman 1995;Dong and Sutton 2005;Ba et al. 2014;Wills et al. 2019).This "internal" AMV is often linked to changes in the strength of the Atlantic meridional overturning circulation (AMOC) and associated northward heat transports on decadal time scales (Delworth et al. 1993;Ba et al. 2014;Menary et al. 2015a;Ortega et al. 2017;Hand et al. 2020).Indeed, the convergence of northward-heat-transported anomalies associated with the AMOC is often thought to be the dominant driver of decadal SST changes, especially in the subpolar North Atlantic (Zhang 2017).
However, a range of surface processes and feedbacks involving changes in clouds, winds, and surface fluxes have also been highlighted to be playing an important role in shaping the spatial pattern of AMV (Martin et al. 2014;Yuan et al. 2016;Brown et al. 2016).A key role for the atmosphere in shaping AMV is also highlighted by the fact that, in many models, persistent heat-loss anomalies over the subpolar North Atlantic, and the subsequent modulation of subsurface density anomalies and AMOC, are associated with positive winter North Atlantic Oscillation (NAO) conditions (Robson et al. 2012;Delworth et al. 2017;Ortega et al. 2017;Kim et al. 2020).Indeed, recent studies based on observational data have found evidence of dynamical coupling between the NAO, AMOC, and AMV, with the NAO leading the AMV by 15-20 years (Li et al. 2013;Sun et al. 2015Sun et al. , 2019)).
Recently, many studies have also argued that external radiative forcings have played a key role in the AMV.External forcing from natural sources such as solar variability and volcanic aerosols are thought to modulate AMV (Knudsen et al. 2014;Swingedouw et al. 2015;Wang et al. 2017).In particular, tropical volcanic eruptions can inject sulfate aerosols into the stratosphere, where they can directly warm the tropical lower stratosphere by absorbing longwave radiation.The resulting increased meridional temperature gradient leads to a positive winter NAO and subsequent AMOC strengthening (Otterå et al. 2010).Multidecadal variability in anthropogenic aerosol emissions in the twentieth century has also been proposed as an important driver of the observed AMV (Bellucci et al. 2017; Murphy et al. 2017;Bellomo et al. 2018; Undorf et al. 2018), through a direct response of North Atlantic SST to local radiative forcing, (Booth et al. 2012;Bellomo et al. 2018;Mann et al. 2020) or indirectly through ocean circulation changes driven by externally forced buoyancy changes (Menary et al. 2013(Menary et al. , 2020b)).
Although the AMV is a key feature of the observations and simulations, there are still ongoing debates regarding the details of the mechanism and the ultimate drivers.In particular, the relative importance of atmospheric forcing versus ocean dynamical changes in driving the AMV, and the role for external forcings, remain topics of ongoing debate.One key reason for this debate is the diversity of the internal AMV simulated by models.For example, some models show that subpolar buoyancy forcing related to the persistent phases of the NAO plays a significant role in driving lowfrequency AMOC variability and AMV (Delworth et al. 2017).Some models also show significant evidence for coupling between the ocean and atmosphere; for example, SST anomalies can modulate atmospheric circulation (Marshall et al. 2001;Wills et al. 2019) and SST-driven rainfall changes can induce density anomalies (Vellinga and Wu 2004;Menary et al. 2012).However, in some models decadal AMOC variability is primarily an ocean-only mode with no active role for the atmosphere (Kwon and Frankignoul 2014;Sévellec and Huck 2015).Others yet question the role of ocean circulation entirely, arguing that the AMV can be explained purely from atmospheric forcing (Clement et al. 2015).Unfortunately, observations remain too sparse and short for in-depth analysis to settle this debate.
Irrespective of the mechanism, it is clear that the internal AMV in most models does not explain the observational record.Indeed, internally simulated AMV usually have significant deficiencies in the time scale, magnitude, and spatial pattern of AMV (Cheung et al. 2017).However, the role of external forcings in driving the AMV is also not well understood, and there are significant differences between the simulated forced AMV and that observed (Zhang et al. 2013;Yan et al. 2018;Kim et al. 2018;Andrews et al. 2020).Unfortunately, it is not currently known whether the deficiencies in reproducing the observed AMV are due to uncertainty in the external forcing, in the underlying internal mechanism, or in how the internal mechanisms of AMV and external forcings interact.Hence, to better understand the observed AMV and improve predictions there is a need to improve our physical understanding of what drives differences in the simulated AMV.
Recently, it has been shown that both the spatial pattern of the internal AMV (Menary et al. 2018) and historical AMV (Andrews et al. 2020) in the Hadley Centre Global Environmental Model in the Global Coupled configuration 3.1 (HadGEM3-GC3.1),the U.K. climate model for phase 6 of the Coupled Model Intercomparison Project (CMIP6), differs when using different oceanic and atmospheric resolutions of the model.There are indeed reasons to expect that models with different resolutions may simulate AMV differently.For example, higher ocean resolution improves the representation of boundary currents such as the Gulf Stream, can reduce SST biases, and can affect regions of water-mass transformation and sinking (Danabasoglu et al. 2014;Katsman et al. 2018;Sein et al. 2018).However, the reasons why the simulated AMV in HadGEM3-GC3.1 is different for different resolutions (i.e., the dominant mechanisms and whether they are related to the resolution) remain unknown.
Given that the versions of HadGEM3-GC3.1 presented in Menary et al. (2018) are extremely similar (i.e., they use the same physical and numerical schemes), they provide an opportunity to understand how the internal AMV can be sensitive to the model used.Therefore, in this study, we aim to characterize the key mechanisms driving the internal AMV within the two versions of HadGEM3-GC3.1 and, particularly, to understand the main causes of the differences.We will first briefly describe the model and its mean state in section 2, before investigating the mechanism of internal AMV in sections 3-5.A discussion of the difference in mechanism between the two models will be included in section 6 and a summary will be presented in section 7.

a. Model description
HadGEM3-GC3.1 was developed by the Met Office and forms part of the United Kingdom's contribution to CMIP6 (Eyring et al. 2016).This model uses the NEMO ocean (Madec and NEMO Team 2016), CICE sea ice (Ridley et al. 2018b), and JULES land surface (Walters et al. 2019) submodels.Details about the development of HadGEM3 are given in Williams et al. (2018).
The two versions of the model used in this study are denoted N96ORCA1 and N216ORCA025 (known as HadGEM3-GC3.1-LLand HadGEM3-GC3.1-MM,respectively, in CMIP6 nomenclature, but referred to here as N216 and N96 hereafter).The medium-resolution N216 model (Williams et al. 2018) has horizontal resolution of ∼60 km in the atmosphere and 0.258 in the ocean.The low-resolution N96 model (Kuhlbrodt et al. 2018) has horizontal resolution of ∼135 km in the atmosphere and 18 in the ocean.The vertical resolution is the same in both models, with 85 pressure levels in the atmosphere, up to 85 km, and 75 depth levels in the ocean.
Although the two models are extremely similar and differ mainly in horizontal resolution, there are some small physical differences between N96 and N216, which arise mainly from the use of different parameter values to compensate for resolution-dependent processes.For example, the N96 atmosphere has different parameter values to compensate for lower maximum winds and limited resolution of gravity wave drag.There are also differences in the ocean submodels; in particular, ORCA1 includes a Gent-McWilliams parameterization for eddy-induced but ORCA025 does not, and ORCA1 also has a slightly lower albedo for snow on sea ice to compensate for weak Arctic sea ice bottom melt.A detailed summary of the differences between the models are included in Kuhlbrodt et al. (2018).
Here we focus on the analysis of the preindustrial (piControl) simulations made for CMIP6.The preindustrial control employs CMIP6-defined forcings which are interannually invariant forcings appropriate for the year 1850.These include well-mixed gases (including CO 2 ), ozone, solar, and natural and anthropogenic aerosols alongside unperturbed land use.Details of these forcings, as well as how they are implemented in the model, are described in Menary et al. (2018).For the N96 and N216 preindustrial control simulations used in this study, various modes of internal climate variability, including the AMV, have already been documented in Menary et al. (2018).The model data of the preindustrial control simulations used in this study are now being published on the Earth System Grid Federation (ESGF).

OBSERVATION-BASED DATASETS
Observations of SSTs were taken from Hadley Centre Sea Ice and Sea Surface Temperature (HadISST; Rayner et al. 2003).Sea level pressure (SLP) was taken from the NOAA Twentieth Century Reanalysis (20CR; Compo et al. 2011).The 1871-1950 climatologies were used for comparison to the control simulations to minimize the effects of varying levels of anthropogenic and natural forcings in the observations.However, the biases are insensitive to the climatological period (not shown).

1) COMPUTATION OF SEAWATER DENSITY
Seawater density used in this study was calculated offline by using the Thermodynamic Equation of Seawater-2010 (TEOS-10) referenced to the pressure level 0 dbar.The temperature component of the density was calculated by using the same equation used for total density, but with time-varying temperature while holding the salinity at the climatological value and vice versa for the salinity component.The "T versus S" variable, which is a measure of the relative importance of the temperature versus salinity component of total density, was calculated taking the difference in the absolute value of the anomalies in the T and S components of density.

2) DEFINITION OF INDICES
The AMV index used here is obtained by first linearly detrending the SST to remove model drift, before taking the area average over the North Atlantic (08-608N, 7.5-758W; marked by the black box in Fig. 2d).Finally, the AMV index is smoothed by a 15-yr moving average.Hereafter, decadal variability will refer to 15-yr smoothed data.Although this index is slightly different to the usual basin mean definition of AMV, which detrends last, extremely similar results are obtained when using the alternate index, as the correlation between the two indices is greater than r = 0.99 in these simulations (not shown).Without detrending first, in N216, there is a linear warming trend of 0.35 K over the 500 years in the AMV index (not shown).N96 is more stable, with a linear warming trend of 0.02 K over the 500 years.
The NAO index used in this study is defined as the difference in standardized winter [December-February (DJF)] SLP using values for the grid cells at 33.748N, 25.678W and 65.078N, 22.738W, approximating the common stationbased definition.The AMOC index used in this study is defined as the maximum AMOC strength at 458N, located at a depth of 1000 m in both models.AMOC at 458N is used here instead of at 26.58N, which is where the RAPID array is located, because the 458N latitude is more strongly connected to variability in the subpolar North Atlantic and less affected by the southward propagation speed of AMOC anomalies, making it more comparable between the models.

3) STATISTICAL SIGNIFICANCE
To assess the statistical significance of the relationship between time series (i.e., cross correlations).We applied a 15-yr block bootstrap method to construct 10 000 samples of the dependent variable.The resulting resampled time series are then smoothed by a 15-yr running mean.Correlations are then calculated between the resampled time series and the independent variable (e.g., the AMV index) to generate a distribution of correlation coefficients.From this distribution, correlation of the original time series is deemed significant if it was greater than the 95th percentile or less than the 5th percentile of the correlation distribution.The same process is used to assess statistical significance for the relationship between a time series and a 2D field (i.e., lagged regressions), only using regression instead of correlations and repeated for each grid point.
A similar method is used to assess significance for the power spectra: 100 000 synthetic AR(1) time series with random noise were generated based on the mean, variance, and estimated lag = 1 correlation of the respective North Atlantic SST time series.A power spectrum was calculated for each synthetic time series to generate a distribution of power at each frequency.A particular frequency is deemed significant if the spectral power of the original time series is greater than the 95th percentile.

c. Model mean state and biases
We now briefly describe the mean states and biases of the two models for variables that are relevant for North Atlantic climate.In both models, there is a warm bias of 1 K along the northern edge of the subpolar gyre (Figs.1a,b).However, N96 shows a cold surface bias of up to 26 K in the central subpolar gyre, which is not present in N216.The cold bias in N96 is a common issue in lower-resolution models (Danabasoglu et al. 2014), and can be attributed to a North Atlantic current, which is too weak and too zonal (Figs.1e,f).The reduction in cold bias in the higher-resolution model could be because of the better representation of boundary currents, which has been shown to improve eddy interaction boundary currents and the center of the gyre (Born et al. 2016).The DJF SLP is also biased high at high latitudes in both models, but is greater in N96 (Figs. 1c,d).In the previous study documenting climate variability in the same simulations, it was found that the NAO in N216 has higher multidecadal NAO variability than in N96 (Menary et al. 2018).
Figures 1e and 1f show that ocean currents are faster in N216 than in N96, especially the boundary currents associated with the subpolar gyre (SPG).Additionally, boundary currents such as the Gulf Stream, East Greenland Current, and the Atlantic inflow are more sharply resolved in N216 than N96, which is common for higher-resolution models (Marzocchi et al. 2015).As expected, the position and magnitude of sea surface height (SSH) gradients is in good agreement with the top 300-m ocean currents (Figs.1g,h), and so can be used as a proxy for surface currents.
Figures 1g and 1h show the AMOC streamfunction in depth space.In both models, the upper cell extends down to a depth of 3 km, and maximum transport occurs at approximately 1000 m.N216 has a maximum overturning strength of 15.0 Sv (1 Sv ≡ 10 6 m 3 s 21 ) at 26.58N, while N96 is slightly weaker, 14.3 Sv at 26.58N.This is weaker than the observed mean AMOC strength of 16.7 Sv (2004-18) obtained from the RAPID array also at 26.58N (Smeed et al. 2019).The southward flowing lower limb, at depths of 3000-5000 m, which corresponds to the North Atlantic Deep Water has a maximum of 22.5 Sv in both models, which is weaker than the observed estimation of 25 Sv (Smeed et al. 2018).Nevertheless these models do broadly capture the subpolar AMOC as measured FIG. 1. (a),(b) The SST biases of the N216 and N96 models, respectively, compared to HadISST. (c),(d) The winter sea level pressure (DJF SLP) biases of the N216 and N96 models, respectively, compared to 20CR.(e),(g),(i) The time-mean 0-300-m average ocean current speed, sea surface height (SSH) relative to the Earth's geoid, and Atlantic meridional overturning circulation (AMOC) streamfunction in depth coordinates, respectively, for the N216 model.(f),(h),(j) As in (e), (g), and (i), but for the N96 model.by the Overturning in the Subpolar North Atlantic (OSNAP) program (Menary et al. 2020a).

Characteristics of the AMV
We begin by first characterizing AMV in the models by comparing the AMV time series and spatial pattern between the observations and simulations.
Figure 2a shows that both models simulate pronounced multidecadal variability in the AMV time series.The magnitude of the variability, as measured by the standard deviation of the AMV index, is smaller than in the observations (0.08 K in N96 and 0.09 K in N216 and 0.13 K in the observations).Spectral analysis of the AMV index shows that in both models there is significant variability in the multidecadal band with peaks at periods of 85 and 67 years in N216 and 85 and 20 years in N96 (Fig. 2c).In the observations, much of the variance has period of around 70 years, with power greater than those of the models.However, this observed multidecadal variability is not statistically significant.This uncertainty of the observed variability is likely, in part, due to the short (145 years) time series compared to the time scales of interest.Both models also show significant variability at the interannual time scale, with periods of 2.2 years in both models and also 4 years in N96.Although there is general agreement between the models and observations, they are not directly comparable, because the observations contain varying levels of external forcing that are not present in the preindustrial simulations.
Figure 2b shows that the AMOC 458N index also shows pronounced multidecadal variability.Importantly, the evolution of AMOC 458N closely resembles the AMV index (i.e., The definition of the AMV and AMOC 458N indices are described in section 2. (c) The power spectra of the annual mean (i.e., unsmoothed) NASST time series for HadISST (black), N96 (red), and N216 (blue), estimated using Welch's method.Dots mark periods when the spectral power is significantly different than red noise at the p # 0.1 level based on a one-sided Monte Carlo test (see section 2 for details).(d) The regression of annual mean SST anomalies onto the 15-yr smoothed AMV index in HadISST observations in colors.SST anomalies have had the climatological mean subtracted and had been linearly detrended at each grid point before regression.Black contours show the variance explained, computed as the square of the correlation coefficient.(e),(f) As in (d), but for the N216 and N96 models, respectively.The black box in (d) highlights the region used for defining the AMV index.The gray dashed boxes in (f) marks the polar, subpolar, and subtropical North Atlantic regions analyzed in the text.
the maxima and minima of the AMOC index occur at nearly the same years as the AMV index).This close resemblance suggests that AMOC variability plays a key role in the internal AMV in both models.
Both models also broadly capture the key features of the AMV spatial pattern (Figs.2d-f).Specifically, the largest AMV anomalies are located in the subpolar North Atlantic and a "horseshoe" shape pattern extends into the subtropics.However, the observations show a much more spatially coherent, basinwide AMV than the simulations.In HadISST, the AMV index explains approximately 30% of the interannual local SST variance across the whole North Atlantic (as measured by the square of the correlation coefficient), whereas in the simulations, the AMV index only explains ∼10% in N96 and ∼20% in N216 of SST variance in the subpolar North Atlantic, with very little correlation in the subtropical and tropical North Atlantic.Relatedly, the magnitude of subtropical AMV anomalies is also weaker in the simulations compared to observations.An additional difference is that there are large anomalies in the Greenland-Iceland-Norway (GIN) Seas and Barents Sea in the simulations which are not present in the observations.

Evolution of the ocean and atmosphere associated with the AMV
To begin to understand the mechanisms governing the internal AMV, we first identify the large-scale changes associated with the AMV and how they vary through time.To that end, Fig. 3 shows the spatial pattern of SST, winter (DJF) SLP, and the AMOC streamfunction regressed onto the AMV index at different lags.Figure 4 shows the cross correlation between the AMV index and several indices of atmospheric and oceanic circulation changes.Additionally, Fig. 4 also shows the relationship between AMV and Labrador Sea subsurface density (700-2000 m, region shown by western portion of black box in Fig. 8), as it is often used as an indicator of buoyancy forced AMOC changes (Robson et al. 2016;Ortega et al. 2017).
The evolution of the spatial pattern of SST associated with the AMV appears broadly similar between N216 andN96.Warm anomalies appear first in the GIN Seas and Barents Sea, at lag = 225 (25 years before peak AMV), before reaching a maximum at lag = 210 (Fig. 3).The subpolar region begins warming at lag = 210 and reaches its maximum at lag = 0.A horseshoe-like pattern emerges in the subtropics at lag = 0, but it is relatively weak compared to observations and does not expand into the tropical North Atlantic.These weak warm subtropical anomalies persist for up to 10 years after peak AMV.From lag = 10 onward, the Polar region reverses into being anomalously cool.
However, there are some important differences in the evolution of subpolar SST between the models.For example, at lag = 210 and lag = 0, there are significant SST anomalies along the Gulf Stream extension in N216.These SST anomalies then appear to propagate along the North Atlantic current into the Eastern SPG.In N96 during the same lead time, strong SST anomalies instead appear in the Labrador Sea and along the southern coast of Greenland instead.These differences in SST evolution suggest a difference in ocean circulation variability between the models.
The cross correlations and the regression patterns of DJF SLP show that there are significant atmospheric circulation changes associated with the AMV in both models.In the atmosphere, spatial pattern of SLP shows NAO-like SLP anomalies develop over the subpolar North Atlantic at lag = 210 in both models (Fig. 3).However, there are key differences in SLP evolution between the models.In N216, the magnitude of the SLP anomalies is larger and its pattern resembles the canonical positive winter NAO more than in N96.Furthermore, negative NAO-like SLP anomalies follow at lag = 10 in N216, but there is no significant reversal in N96.These negative NAO anomalies likely contribute to the AMV phase reversal in N216.It is important to note that, at lag = 210, when the NAO-like anomalies are largest and therefore surface winds strongest, the SST continues to warm.This indicates that at the decadal time scale the NAO is likely not directly forcing the SST changes via surface cooling, which is consistent with previous studies (Delworth et al. 2017).
Figure 4 summarizes the difference in the role of atmospheric circulation changes.In N216, a positive NAO leads the AMV by 15 years and the subsequent negative NAO lags the AMV by 15 years.In contrast, in N96, there is no significant correlation between the NAO and AMV at any lag time.However, this difference in the correlation of NAO index and the AMV is largely due to the difference in SLP pattern; the NAO index does not do a good job of capturing the changes in N96.Focusing on the Icelandic low index (i.e., northern node of the NAO index) instead shows that there are significant atmospheric circulation changes leading the AMV in N96, but its magnitude is smaller than that found in N216 and there is no significant reversal following the AMV.
Both models also show significant correlation between the AMOC and the AMV.The AMOC anomalies are largest and appear first in the subpolar region (458N) at lag = 215 in N216 and lag = 215 in N96 (Fig. 3).Although there appears to be no significant AMOC anomalies north of 508N in N216 when expressed in depth space, significant AMOC anomalies do form at high latitudes when expressed in density space (not shown).The cross correlations show that maximum AMOC at 458N leads AMV by 5 years in both models (Fig. 4).However, there are differences in the evolution of AMOC anomalies in the models.The regression patterns show that in N216, AMOC anomalies associated with AMV are latitudinally coherent; that is, there are strong and significant AMOC anomalies ranging from the subpolar region (508N) to the tropics (08N) at lag = 210 and lag = 0. On the other hand, in N96, the AMOC anomalies associated with AMV are less coherent latitudinally with significant AMOC anomalies only reaching as far south as 158N at lag = 0.
Additionally, in both models, Labrador Sea subsurface density leads the AMV by 5 years and its evolution closely resembles AMOC changes.This suggests a close link between Labrador Sea subsurface density anomalies and AMOC in both models, consistent with previous studies (Robson et al. 2016).The drivers of subpolar subsurface density anomalies and their role in shaping the AMV will be explored in detail in later sections.
Having shown that there are statistically significant links between the AMV and large-scale ocean and atmospheric circulation changes.We will next investigate the relative importance of ocean transport versus surface fluxes in controlling regional ocean heat content (OHC) changes.

Processes governing regional OHC changes
Previous studies have shown that different processes drive AMV in different regions (Martin et al. 2014;Brown et al. 2016;Yuan et al. 2016).Therefore, to broadly identify the dominant processes shaping the internal AMV in different regions, we perform a cross correlation of ocean temperature averaged over the top 100 m (T100 m) with ocean heat transport (OHT) convergence and integrated total surface heat fluxes (SHF, accounting for net surface longwave and shortwave radiation alongside sensible and latent heat flux).Based on the evolution of SST in Fig. 3, we split the North Atlantic into three regions: polar (658-808N), subpolar (458-658N), and subtropical (158-458N).These regions are shown by the dashed gray box in Fig. 2f.For simplicity, we focus on T100 m as a measure of upper-ocean heat content because climatological mixed layer depths differ in the three regions (Figs.9a,b).However, the results of the analysis are broadly consistent when using SST or 0-300-m heat content, although with different lead/lags (not shown).
Figure 5 shows that the relative importance of either ocean advective processes or SHF for controlling AMV varies by region.In both the polar and subpolar regions (Figs.5a,b), both models show that in the years leading T100, OHT is positively correlated while SHF is negatively correlated with T100; i.e., OHT convergence drives upper-ocean warming, while SHF acts as a damping influence.This relationship in the subpolar region is reversed in the years following T100 in both models.However, there are also differences between the models.In the polar region, the magnitude of correlation for both OHT convergence and SHF is smaller in the N216 model than in N96, suggesting that there may be additional processes at play in N216.In contrast, surface fluxes instead play a dominant role in the subtropical region (Fig. 5c).In both models, SHF is positively correlated with T100 with a lead of 10 years, while OHT convergence is negatively correlated at the same lag.
The latitudinal variations in processes driving the internal AMV indicate that while slow ocean circulation changes are important in decadal variability, the basinwide AMV pattern can only be understood as a combination of changes in both ocean circulation and the processes governing surface fluxes.

a. The role of surface fluxes in ocean heat content changes
To highlight the point that the upper-ocean warming associated with internal AMV is more than just an ocean advective process, we investigate the mechanism behind the surfacedriven warming before exploring the details of ocean circulation changes in more details in later sections.In Fig. 6, net downward surface shortwave radiation (DWSW), total cloud fraction, and winter (DJF) sea ice concentration are regressed onto the AMV index.
Figure 6 shows that there are significant DWSW changes associated with the AMV in the polar, eastern subpolar, and subtropical regions in both models.In the polar region, positive DWSW anomalies in the Barents and GIN Seas reaches a maximum at lag = 210, before slowly weakening to lag = 0, and then disappearing by lag = 10.In the eastern subpolar and subtropical North Atlantic, DWSW anomalies appear at lag = 210, and reaches a maximum at lag = 0. Evolution of the DWSW anomalies (Fig. 6) is similar to the SST (Fig. 3) in some areas of the polar and subtropical regions.In the polar region, positive DWSW and SST anomalies in the GIN and Barents Seas appear at the same lag time.The weak warming in the eastern subtropical region also coincide with the weak positive DWSW anomalies However, the processes behind radiative changes in the subtropics are different from those in the polar region.In both models, positive DWSW anomalies at high latitude are correlated with the decrease in sea ice concentration, whereas positive DWSW anomalies at midand low latitudes are instead correlated with a decrease in total cloud concentration (Fig. 6).This covariation of DWSW and cloud is consistent with a positive SST-cloud feedback mechanism (Brown et al. 2016).
In the polar region, although increased DWSW is associated with a decrease in sea ice concentration, the loss of sea ice can also induce turbulent heat loss via increased air-sea interaction (Taylor et al. 2018).Figure 5a shows that, overall, the total SHF in the polar region damps upper OHC.FIG. 5. Simplified regional heat budget. (a) The cross correlation between the 15-yr smoothed 0-100-m ocean temperature (T100) with ocean heat transport (OHT) convergence (red) and total downward surface heat flux (SHF) averaged across the polar North Atlantic (658-808N, 808W-508E).(b),(c) As in (a), but for the subpolar (458-658N, 808W-508E) and subtropical (158-458N, 808W-508E) North Atlantic, respectively.The three regions are marked by the gray dashed box in Fig. 2e.OHT convergence is defined as the difference between the northward full-depth northward ocean heat transport across the southern boundary minus the northern boundary of the corresponding regions.Circles show correlation coefficients significant at the p # 0.1 level based on a two-sided Monte Carlo test (see section 2 for details).
Therefore, taken together, the cross correlation (Fig. 5a) and DWSW regression (Fig. 6) indicate that the turbulent heat loss is the primary driver of total SHF over the polar North Atlantic at this time.The magnitude of sea ice and the associated polar DWSW response greater in N216 than N96 (Fig. 6), which suggests that although turbulent heat loss is the primary driver in this region, radiative changes still likely have a stronger role in shaping the evolution of polar SST in N216 compared to N96.

b. The role of ocean heat transport in subpolar OHC changes
We now turn our attention to the subpolar North Atlantic, where ocean advective processes are most important.To further understand how heat is being advected into the subpolar, we examine the contribution across the northern (658N) and southern (458N) boundary.Figures 7a and 7b show the time series of these OHT components.Figure 7c shows the cross correlation between total subpolar OHT convergence and contribution across the northern or southern boundaries.
Figure 7a shows that decadal changes in subpolar OHC are largely driven by heat transport across the southern boundary in N216.Standard deviation of the decadal OHT is 0.028 PW across the 458N and 0.008 PW across 658N.Furthermore, OHT variability across the northern boundary is anticorrelated with the total OHT convergence and is, therefore, acting as a weak damping influence on subpolar heat content (Figs.7a,c).This relationship between the northern and southern boundary is consistent with the coherent AMOC changes shown previously in the N216 model (see Fig. 3); that is, the strengthened AMOC is both moving more heat into the subpolar North Atlantic across the southern boundary and moving more heat out across the northern boundary.
In contrast, subpolar OHC is driven by a combination of OHT across both the northern and southern boundary in N96.The magnitude of OHT variability across the northern boundary is much closer in magnitude with the southern boundary in N96; the standard deviation of the decadal OHT is 0.023 PW across 458N and 0.020 PW across 658N.Moreover, OHT across both the northern and southern boundaries are positively correlated with total OHC convergence near lag = 0, which indicates that OHT across both boundaries contribute to the total subpolar heat content.The cross correlations also show that the southern contribution leads total warming by 10 years and northern contribution lags by 10 years (Fig. 7c).
To understand the difference in the relative importance of OHT across the northern and southern subpolar boundary between the models, we explored how ocean circulation changes differ in the two models.To that end, Fig. 8 shows the spatial patterns of SSH anomalies (which is measure of surface geostrophic currents) and subsurface (700-2000 m) density anomalies, regressed onto the AMV at various lags.
Figure 8 shows that there is a marked difference in the evolution of SSH between the models.In N216, negative SSH anomalies, associated with a strengthened SPG, are largely confined to the western subpolar North Atlantic and the Labrador Sea.Positive SSH anomalies appear in the Gulf Stream extension at lag = 220, and then in the eastern SPG at lag = 0.In N216, the positive SSH anomalies from lag = 210 to lag = 0 resemble the positive SST anomalies shown in Fig. 3.Both these variables show positive anomalies that appear to move northward from the Gulf Stream region into the central subpolar North Atlantic.This covariability indicates that these positive SSH anomalies are, in part, a signature of thermosteric changes associated with upper OHC changes.However, in N96, the magnitude of negative SSH anomalies over the SPG are larger and more expansive than in N216, with anomalies stretching from the Labrador Sea to the south of Iceland (Fig. 8).These negative SSH anomalies appear at lag = 220 and reaches their maximum at lag = 0. On the other hand, there are no significant SSH anomalies associated with the Gulf Stream and North Atlantic Current in N96.
Figure 8 also shows that there are significant differences in the propagation of subsurface density anomalies.In N216, positive subsurface density anomalies appear at lag = 220 and propagate southward from the subpolar North Atlantic along the western boundary into the tropical North Atlantic.Density anomalies then cross the tropical Atlantic at lag = 0 and are consistent with subsurface fingerprints of AMOC adjustment (Johnson and Marshall 2004).In contrast, in N96, the subsurface density anomalies which appear in the subpolar gyre at lag = 220 expanding into the central SPNA and then broadly down the western Atlantic.These subsurface density anomalies in N96 reach the subtropics much later than in N216, consistent with the AMOC anomalies in N96, which are more confined to the high latitudes in comparison to N216 (Fig. 3).

Drivers of North Atlantic density anomalies
Having shown that ocean circulation and subsurface density changes are indeed associated with the internal AMV in both models (Figs. 4,5,and 8), we will now assess the drivers of these subsurface density anomalies.Previous studies have often linked AMOC variability to subpolar density anomalies and oceanic deep convection (Menary et al. 2015a;Robson et al. 2016;Ortega et al. 2017).However, the key regions of deep convection vary between models (Ba et al. 2014), and both ocean and atmospheric processes can lead to the formation of density anomalies.Therefore, we will first compare the deep convection regions between the models, before exploring the changes in surface density and the role of surface flux changes.

a. Deep convection and upper-ocean density
Figure 9 shows the mean and standard deviation of the March mixed layer depth (MLD), which is a measure of deep ocean convection.In both models, the regions with the largest mean MLD (Figs. 9a,b) are also the regions with the largest standard (Figs.9c,d).The standard deviation is used here to define the convective regions since we focus on the variability.There are two regions of wintertime deep convection in both models, and the locations of both these convective regions are very similar between the models.Deep convection is located in the GIN Seas (688-808N, 158W-158E, blue box in Figs.9c,d) and in the Western Subpolar Gyre (WSP, 538-658N, 628-308W, black box in Figs. 9c,d).The WSP site includes both the Labrador Sea and the Irminger Sea.These regions of large MLD variability are largely in agreement with observations from Argo profiles (Holte et al. 2017).The mean and variability of MLD in both the WSP and GIN Seas is larger in N216 than in N96, which is consistent with the stronger stratification and hence, resistance to convection, in N96 in both regions (Fig. 9e).
To understand the processes driving seawater density changes (see section 5b), we investigate the evolution of upper-ocean (0-700 m) density and their interaction with the oceanic deep convection regions.Figure 10 shows the spatial pattern of March MLD, upper-ocean density, and the degree of T versus S density control (see section 2 for details of this variable) regressed onto the AMV index.
Figure 10 shows that in both models, surface (0-700 m) dense anomalies begin to form across the subpolar North Atlantic at lag = 220 and reaches their maximum at lag = 210.This coincides with the deepening of WSP March MLD and also the origin of density anomalies at depth (see Fig. 8).This suggests that deep convection, associated with MLD deepening, may be the process by which subsurface density anomalies are formed in both models, although this may not be the only process through which deep water formation occurs (Katsman et al. 2018).
Figure 10 also shows that some of the positive surface density anomalies appear to originate from the Arctic.There are large positive density anomalies in the Arctic at lag = 240, which then appears along the eastern and southern coast of Greenland at lag = 220, crossing both the Fram Strait and Denmark Strait.These density anomalies appear to link the Arctic region to the WSP region, suggesting that Arctic density anomalies are propagating into the subpolar North Atlantic along the East Greenland Current (EGC) pathway and contributing to densification of the WSP at lag = 210.
The Arctic density anomalies are dominated by salinity in both models but the overall density anomalies are substantially larger in N96 than in N216.The transport of this salinitycontrolled density anomaly is associated with a weakening of the EGC and, hence, a reduction in Arctic freshwater export into the SPG, rather than the advection of negative freshwater anomalies via background flow (not shown).A similar mechanism was found by Jungclaus et al. (2005) in the ECHAM5/MPI-OM model, where AMOC variability is linked to the storage and release of Arctic freshwater.In their study, freshwater export from the Arctic is controlled by EGC changes, which are, in turn, driven by anomalous temperature and salinity transport from the Atlantic, hence forming a feedback loop.

b. Atmospheric forcing of subsurface density
Although we have identified an oceanic driver of density changes, we have also shown that changes in atmospheric circulation also lead the internal AMV in both models.Therefore, to investigate the role of the atmosphere, we examined the extent to which air-sea exchange of heat and freshwater in the WSP region are important drivers of subsurface ocean density.Although deep convection in the WSP region appears to vary coherently, Menary et al. (2020a) has shown that different surface water mass transformations occur in the Labrador and Irminger Seas.Therefore, we explore the role of these two different regions separately.
Figures 11a and 11b show that increased evaporation over the Labrador Sea lead subsurface density anomalies by a few years in both models.However, in the Irminger Sea, increased turbulent fluxes lag subsurface density by 10 years in N216 and 20 years in N96.This lead-lag relation suggests that surface forcing over the Labrador Sea rather than the Irminger Sea is key to driving subsurface density changes on these time scales.This does not rule out the importance of the Irminger Sea for the formation of deep density anomalies, only that atmospheric forcing is not an important driver in this region on decadal time scale.At first glance, this result appears to contrast with Menary et al. (2020a), who argued that Irminger Sea surface forcing of subsurface density dominates at annual time scale in the same simulation.However, cross correlation of the annual-mean time series, rather than the decadally smoothed is more in line with Menary et al. (2020a) (not shown), indicating that the dominant forcing of subsurface ocean density changes may be time scale dependent.
Although Labrador Sea surface cooling leads subsurface (700-2000 m) density anomalies in both models, the origin of this air-sea coupling differs between the models.Figure 11c shows that, in N216, increased surface wind speed leads subsurface density by 3 years.Furthermore, the significant relationship with surface wind coincides with turbulent cooling.SST anomalies then lag subsurface density by 15 years, consistent with a role for ocean advection as a result of the AMOC strengthening.
In contrast, surface winds over the Labrador Sea are not significantly correlated with subsurface density changes in N96 (Fig. 11d).Instead, the correlation between subsurface density and the ocean variables is higher in N96 than in N216, and largely in phase with the THF.For instance, SSS leads subsurface density by 1 year, and SST lags by 2 years.These relationships could suggest that anomalous turbulent surface cooling in N96 is actually the response to convection driven warming which is ultimately caused by salinity-controlled density anomalies (i.e., deep convection leads to warmer subsurface water being brought to the surface, Gelderloos et al. 2012).Furthermore, in N96, there is significant negative correlation between the SSS and WSP subsurface density at positive lags.This is consistent with the fresh, buoyant anomalies that appear in the Arctic following the AMV (Fig. 10) and likely contribute to the reversal of the AMV.

Discussions
We have shown that, although the two resolutions of HadGEM3-GC3.1 studied here appear to simulate relatively similar internal AMV variability, there are significant differences in the overall mechanism of AMV between the resolutions.For instance, there are major differences in the relative importance of the processes driving AMV in the models; atmospheric forcing of subsurface density and subsequent ocean circulation changes play a key role in the N216 mechanism, whereas N96 exhibits a more ocean-only salinity-driven mechanism.There are also differences in the relative role of ocean circulation changes, with N216 dominated by basinwide coherent AMOC changes and N96 having a much larger role for persistent anomalies in the subpolar gyre circulation.Furthermore, the AMV variability in N216 appears to be more "coupled" with significant persistent NAO anomalies both leading and following AMV.However, a key question is what controls these differences.
One feature that appears to explain the differences in the ocean circulation and hence SST pattern associated with the AMV is the propagation of subsurface density anomalies.In short, in N216 subpolar subsurface density anomalies propagate southward along the western boundary, reaching the subtropics much earlier in the phase evolution of AMV than in N96 (Fig. 8).This is consistent with the more coherent AMOC anomalies (Fig. 3) and the changes in the Gulf Stream circulation (Fig. 8) in N216.These circulation changes are, in turn, consistent with the subsequent development of SST anomalies in the Gulf Stream extension (Fig. 3), subpolar OHT anomalies that are dominated by changes to the AMOC to the south (Fig. 7), and the subsequent evolution of SST anomalies along the North Atlantic Current in N216 (Fig. 3).In contrast, in N96 subsurface density anomalies persist in the subpolar latitudes for longer and propagate into the subtropics along interior pathways, reaching the subtropics later than in N216 (Fig. 8).Consequently, the circulation anomalies and the development of SST anomalies are more centered in the subpolar North Atlantic with less connection to the Atlantic further south (Fig. 3).Furthermore, the negative SSH anomalies in N96 lag positive NAO anomalies by 10 years (Figs. 3   and 8), so it is unlikely that wind forcing plays a key role in the SPG evolution.Hence, SPG changes are likely to be driven by the subsurface density anomalies instead.A possible mechanism is the interaction of abyssal flow with sloping bathymetry (Yeager 2015(Yeager , 2020)), but further investigations will be needed to unravel the relevant mechanisms in N96.
Another key difference between the models is that the NAO-AMV linkage is much stronger in N216 than in N96, which has implications for the overall evolution of the AMV (Delworth et al. 2017;Sun et al. 2015).However, we still do not fully understand why this is or why there is significantly more decadal NAO variability in N216 than N96 (see Menary et al. 2018).One possible explanation for this could be a difference in ocean-driven turbulent heat fluxes, which has been argued to drive the NAO (Rodwell et al. 1999).As in the observations (e.g., Gulev et al. 2013), the simulated AMV is associated with ocean-driven THF over the North Atlantic current and subpolar North Atlantic in both N216 and N96 (not shown).However, these THF anomalies are larger in N216 and are, therefore, consistent with a stronger dynamical coupling between the NAO and AMV.We also find some evidence of remote teleconnection with the Pacific (not shown), and this relationship is also different between the models.Previous studies have also linked decadal NAO variability to sea ice changes (Yang et al. 2016) and stratospheric coupling Baldwin and Dunkerton (2001).Unfortunately, due to the range of processes, it is difficult to definitively disentangle the drivers and further analysis will be required.
The specific role of resolution in explaining these differences between the models, as opposed to mean-state differences, is unclear.Indeed, there is evidence to suggest that mean-state differences are important.For example, changes in stratification and SSTs may affect the formation of subsurface density anomalies (Yeager and Danabasoglu 2012;Menary et al. 2015b).Reduced SST biases in the North Atlantic may also contribute to the relationship between SST anomalies and the NAO (Lee et al. 2018).Furthermore, mean differences in the stratification and circulation between the models are also likely to affect the propagation of subsurface density anomalies.Indeed, the difference in AMOC meridional coherence on decadal time scales found here is consistent with the relationship found with the stratification in the Labrador Sea in CMIP5 models (Ortega et al. 2021).However, as already mentioned, there are several reasons to expect higher-resolution models to better simulate important aspects of the Atlantic system (Anstey et al. 2013;Roberts et al. 2016).Studies have shown that ocean resolution can influence both subsurface density propagation pathways (Spence et al. 2012) and model mean state (Menary et al. 2018;Caldwell et al. 2019;Jackson et al. 2020).Atmospheric resolution can influence NAO eddy feedback (Scaife et al. 2019) and the strength and position of the storm track (Zappa et al. 2013).Nevertheless, it is difficult to attribute any differences in simulated internal AMV to resolution alone.Moreover, simulations of even higher-resolution models (such as those submitted for HighResMIP) are not long enough to robustly compare their mechanisms of decadal variability.Further systematic study using more models, further idealized simulations and longer high-resolution simulations will be needed to isolate how different aspects of the simulated climate impact on AMV.
Although we have broadly characterized the key processes driving the internal AMV, there are still some unanswered questions, such as the drivers of Arctic circulation changes and realism of mechanism.For the Arctic mechanism, Jungclaus et al. (2005) argued that the ocean circulation variability controlling Arctic-Atlantic exchanges is largely an ocean-only mode.This is broadly consistent with our models, as we also do not see significant atmospheric circulation anomalies in phase with propagation of Arctic density anomalies into the Atlantic.Further analysis also show that Denmark Strait overflow transport changes lead AMOC changes (not shown), suggesting that high-latitude deep water formation may have some control over subsurface density changes farther south.Nevertheless, we still do not know what the key elements of this Arctic mechanism are or the processes involved.There are other details of the mechanism we have not looked into, such as the role of salinity imported from the subtropics via the North Atlantic Current, which can feed back onto the AMOC (Vellinga et al. 2002) and affect Nordic seas deep water formation (Glessmer et al. 2014).We also have not diagnosed the importance of vertical mixing for generating the AMV pattern, which have been previously highlighted as an important driver of regional SST changes (Yamamoto et al. 2020).Further targeted analysis and experiments, outside the scope of this study, will be required to resolve these outstanding questions.
Finally, on the realism of the mechanism, there is some similarity between the simulated AMV in both versions of HadGEM3-GC3.1 and the changes seen in the observations, especially for the AMV phase transition in the 1990.In particular, that event was led by a persistent positive NAO, and there is significant evidence for a change in the ocean circulation (Robson et al. 2016;Sutton et al. 2018).However, there are several features of the N216 mechanism which is in better agreement with the relatively short observational records, such as the stronger NAO-AMV linkage (Sun et al. 2015;Li et al. 2013;Sun et al. 2019) and pattern of SSH anomalies Desbruyères et al. (2021), such as the Gulf Stream signal (McCarthy et al. 2018).Nonetheless, there are significant differences to observations in both models.In particular, the internal AMV accounts for little of the tropical Atlantic variance in comparison with observations (Fig. 2).As the forced AMV signal in HadGEM3-GC3.1 (i.e., the models used here) is located primarily in the tropics (Andrews et al. 2020), and accounts for a similar fraction of the variance, this could suggest that the observed AMV is consistent with both forced and internal modes as suggested previously (Ting et al. 2014;Tandon and Kushner 2015;Watanabe and Tatebe 2019;Qin et al. 2020).However, we also find that there is significantly more SST variability in the subpolar North Atlantic in the piControls than observed (not shown).Therefore, it is difficult to compare the relative importance of decadal variability.Why the models overestimate this variability, and whether it is related to the signal-to-noise problems highlighted in the North Atlantic (Scaife and Smith 2018), is currently not known.Therefore, it is important to further improve our understanding of the variability in the subpolar North Atlantic as well as to continue to evaluate our models.

Summary and conclusions
This study presents an analysis of the Atlantic multidecadal variability (AMV) in 500-yr preindustrial control simulations with the HadGEM3-GC3.1 model at both N96ORCA1 (∼135 km in the atmosphere, 18 in the ocean) and N216ORCA025 (∼60 km in the atmosphere, 0.258 in the ocean) resolutions.We have characterized the variability in both models, and deduced the key processes leading to the simulated internal AMV.Our findings are summarized as the following: • Both versions of HadGEM3-GC3.1 simulate internal AMV-like variability with a 60-80-yr time scale that broadly compares with observations.The spatial pattern of the AMV is also generally well simulated in models, although it is weaker in the tropics and subtropics than observations.• Low-frequency changes in both ocean and atmosphere circulation lead the AMV in both models.AMOC anomalies lead the AMV by 5 years and NAO-like winter SLP anomalies lead AMV by 10 years in both models.However, the SLP anomalies are larger, more persistent, and resemble the canonical winter NAO more in N216.
• The main drivers of upper-ocean heat content (OHC) changes vary by region.In the subtropical Atlantic, (158-458N) surface heat flux changes (SHF) are the dominant driver in both models with changes in clouds playing an important role.In contrast, the polar (658-808N) and subpolar (458-658N) North Atlantic regions are primarily driven by ocean heat transport (OHT) convergence.• Details of the ocean circulation changes, and the impact on OHT, differ between the models.Specifically, AMOC changes in N216 reach the subtropics earlier and are more coherent.Significant SST and SSH anomalies also appear in the Gulf Stream extension and, subsequently, subpolar OHT anomalies are dominated by changes from the southern boundary.However, AMV in N96 is associated with much larger and more persistent changes in the subpolar gyre (SPG) and a less important role of OHT at the southern boundary.• Differences in the ocean circulation changes are consistent with the evolution of subsurface density anomalies.In N216, deep (700-2000 m) density anomalies propagate south along the western boundary, consistent with the latitudinally coherent AMOC anomalies in N216.However, in N96, deep density anomalies persist for longer in the subpolar North Atlantic, and propagate along an oceanic interior pathway and reach the subtropics much later in the AMV phase than in N216.• The relative importance of oceanic and atmospheric drivers of subsurface density is different between the models.In N216, surface winds associated with the positive winter NAO is the dominant driver of heat loss in the western subpolar North Atlantic and subsequent subsurface density anomalies.However, in N96 upper-ocean salinity-controlled density anomalies that originate from the Arctic and move into the subpolar North Atlantic through the East Greenland Current appear to be the dominant driver.
Taken together, although the simulation of internal AMV appears broadly similar in the two versions of HadGEM3-GC3.1, the relative importance of different processes driving and the shaping the evolution of AMV is different.These differences, and particularly the relationship between the NAO and AMV, will have implications for how AMV impacts on local weather and climate, and could also be important for the predictability of the Atlantic.Our results provide a benchmark for exploring how differences in climate models' underlying internal AMV mechanisms may affect their simulations of the historical AMV and future climate predictions.However, we do not know if the different mechanisms will be impacted differently by external forcing.Therefore, assessing the mechanisms of how the two versions of HadGEM3-GC3.1 respond to external forcing will be the topic of a further study.

FIG. 2
FIG. 2. (a) The time series of the AMV index in HadISST (black), N96 (red), and N216 (blue).(b) The time series of the AMOC 458N index for N96 (red) and N216 (blue).Horizontal lines show the time mean of the AMOC 45N indices.For (a) and (b), annual means are shown with thin lines and 15-yr running means are shown with thick lines.The definition of the AMV and AMOC 458N indices are described in section 2. (c) The power spectra of the annual mean (i.e., unsmoothed) NASST time series for HadISST (black), N96 (red), and N216 (blue), estimated using Welch's method.Dots mark periods when the spectral power is significantly different than red noise at the p # 0.1 level based on a one-sided Monte Carlo test (see section 2 for details).(d) The regression of annual mean SST anomalies onto the 15-yr smoothed AMV index in HadISST observations in colors.SST anomalies have had the climatological mean subtracted and had been linearly detrended at each grid point before regression.Black contours show the variance explained, computed as the square of the correlation coefficient.(e),(f) As in (d), but for the N216 and N96 models, respectively.The black box in (d) highlights the region used for defining the AMV index.The gray dashed boxes in (f) marks the polar, subpolar, and subtropical North Atlantic regions analyzed in the text.
FIG. 3. Regressions of sea surface temperature (SST), winter sea level pressure (DJF SLP), and the AMOC streamfunction in depth space onto the AMV index for the (left) N126 and (right) N96 models.The AMV index (top) lags the fields by 25 years (225 yrs) and then (bottom) leads by 10 years (10 yrs).Stippled regions show where the regression slope coefficient is significantly different to zero at the p # 0.1 level based on a two-sided Monte Carlo test (see section 2 for details).

FIG. 4
FIG. 4. (a)  The cross correlation between the AMV index with the 15-yr running-mean NAO index (light blue), Icelandic low index (blue), Labrador Sea 700-2000-m density (black), and AMOC 458 index (red).(b) As in (a), but for the N96 model.Positive lags show where the AMV index is leading the other variables.Note that the Icelandic low index is multiplied by 21 so the positive Icelandic low index is comparable to the positive NAO index.Circles show correlation coefficients significant at the p # 0.1 level based on a two-sided Monte Carlo test (see section 2 for details).

FIG. 6 .
FIG. 6. Regressions of net downward surface shortwave radiation (DW Surf SW), total cloud fraction, and winter (DJF) sea ice concentration onto the AMV index for the (left) N126 and (right) N96 models.The AMV index (top) lags the fields by 25 years (225 yrs, top) and then (bottom) leads by 10 years (10 yrs).Stippled regions show where the regression slope coefficient is significantly different to zero at the p # 0.1 level based on a two-sided Monte Carlo test (see section 2 for details).

FIG. 7
FIG. 7. (a) The time series of 15-yr running-mean full-depth northward ocean heat transport (OHT) across the northern subpolar boundary (multiplied by 21) (North, 658N, blue), across the southern subpolar boundary (South, 458N, red), and the total subpolar OHT convergence (Subpolar OHT conv, black) for the N216 model.(b) As in (a), but for the N96 model.(c) The cross correlation between subpolar OHT convergence with the with ocean heat transport across the northern boundary (blue) and southern boundary (red).Relationships for the N216 model are shown with solid lines and N96 is shown with dashed lines.Circles show correlation coefficients significant at the p # 0.1 level based on a two-sided Monte Carlo test (see section 2 for details).In all panels, OHT across the northern boundary is multiplied by 21 to show contribution to total subpolar OHT convergence.

FIG. 8 .
FIG. 8. Regressions of sea surface height (SSH) and subsurface (700-2000 m) seawater density on to the AMV index for the (left) N126 and (right) N96 models.The AMV index (top) lags the fields by 25 years (225 yrs) and then (bottom) leads by 10 years (10 yrs).Stippled regions show where the regression slope coefficient is significantly different to zero at the p # 0.1 level based on a two-sided Monte Carlo test (see section 2 for details).

FIG. 9
FIG. 9. (a),(b) The climatological mean of March mixed layer depth (MLD) in the N216 and N96 models, respectively.Contour lines are drawn to highlight mean MLD at lower latitudes.(c),(d) The standard deviation of March MLD in the N216 and N96 models, respectively.The blue box highlights the Greenland-Iceland-Norway (GIN) Seas (688-808N, 158W-158E) convection region and the black box highlights the Western Subpolar Gyre (WSP; 538-658N, 628-308W) convection region.Gray dashed line within the black box divides the WSP convection region into the Labrador Sea to the west and Irminger Sea to the east.(e) The area-averaged time-mean stratification of the WSP (black) and GIN Seas (blue) region in the N216 (solid) and N96 (dashed) models.

FIG. 10 .
FIG. 10.Regressions of March mixed layer depth (MLD), upper-ocean (0-700 m) seawater density and the T vs S control (see section 2) of the upper-ocean density on to the AMV index for the (left) N126 and (right) N96 models at lags.The AMV index (top) lags the fields by 40 years (240 yrs) and then (bottom) leads by 10 years (10 yrs).Note the lags chosen for this figure are different to the lags used in the other regression figures to better highlight the density anomalies and their evolution.Stippled regions show where the regression slope coefficient is significantly different to zero at the p # 0.1 level based on a two-sided Monte Carlo test (see section 2 for details).
11. (a) Cross correlation of WSP subsurface (700-2000 m) density anomalies with upward winter turbulent heat flux (THF, sensible plus latent heat flux, blue) and evaporation minus precipitation (E 2 P, red) averaged across the Labrador Sea.(b) As in (a), but THF and E 2 P are averaged across the Irminger Sea instead.Relationships for N216 are shown by solid lines and N96 in dashed lines.(c) Cross correlation of WSP subsurface density anomalies with upward winter THF (blue), winter surface wind speed (DJF sfc winds, red), annual SST (Ann SST, black), and annual SSS (Ann SSS, gray) for the N216 model.(d) As in (c), but for the N96 model.Circles show correlations coefficients significantly different to zero at the p # 0.1 level based on a two-sided Monte Carlo test (see section 2 for details).