Labrador Sea subsurface density as a precursor of multidecadal variability in the North Atlantic: a multi-model study

. The subpolar North Atlantic (SPNA) is a region with prominent decadal variability that has experienced remarkable warming and cooling trends in the last few decades. These observed trends have been preceded by slow-paced increases and decreases in the Labrador Sea density (LSD), which are thought to be a precursor of large-scale ocean circulation changes. This article analyses the interrelationships between the LSD and the wider North Atlantic across an ensemble of coupled climate model simulations. In particular, it analyses the link between subsurface density and the deep boundary density, the Atlantic Meridional Overturning Circulation (AMOC), the subpolar gyre (SPG) circulation


30
The North Atlantic Ocean is a key component in Earth's climate through, for example, its role in redistributing heat and in taking up excess heat and carbon from the atmosphere. It is also a region that has varied significantly in the past. This is particularly true for the North Atlantic subpolar gyre, that has varied significantly on multi-decadal timescales across a range of different variables (Häkkinen and Rhines, 2004;Holliday et al., 2020;Reverdin, 2010;Robson et al., 2018b). Basin-mean sea surface temperature (SST) over the North Atlantic has also been observed to vary on multi-decadal timescales (Schlesinger
in monsoon regions (Knight et al., 2006;Monerie et al., 2019;Zhang and Delworth, 2006). The North Atlantic is also expected to change significantly in the future due to the effects of climate change, and consequently produce substantial climate impacts on the surrounding regions (Sutton and Hodson, 2005;Woollings et al., 2012). On decadal timescales, it is the interaction between natural variability and externally forced changes that will shape how the Atlantic regions climate will 40 evolve. Therefore, in order to improve predictions of the North Atlantic, it is imperative that we improve our understanding of the processes that control decadal timescale changes in this region.
It has generally been thought that changes in the ocean circulation, and particularly the Atlantic Meridional Overturning Circulation (AMOC), have played a significant role in shaping multi-decadal variability across the North Atlantic (Knight et 45 al. 2005). In particular, changes in the strength of the AMOC, and its related ocean heat transports have been shown to control multi-decadal internal variability in a range of coupled climate models (Danabasoglu, 2008;Dong and Sutton, 2005;Jungclaus et al., 2005;Ortega et al., 2011Ortega et al., , 2015. Changes in AMOC and the wider ocean circulation have been used to explain the changes observed in the subpolar North Atlantic (SPNA), on decadal and longer timescales (Moat et al., 2019). In particular, the SPNA underwent a rapid warming and salinification in the mid 1990s before a decadal timescale cooling and freshening 50 started in 2005, which is consistent with decadal-to-multidecadal variability of the AMOC (Robson et al., 2012(Robson et al., , 2013(Robson et al., , 2016. The recent cooling has been linked to climate impacts over the continents, including heat waves (Duchez et al., 2016), through an effect on the position on the jet stream . A long term relative cooling of the SPNA since ~1850 has also been attributed to a centennial weakening of the AMOC (Caesar et al., 2018;Rahmstorf et al., 2015), an AMOC reduction that most CMIP6 model projections predict to continue in the future (Weijer et al., 2020). However, a lack of direct observations of 55 the strength of the AMOC or the ocean circulation more generally have hindered our ability to make a direct attribution of recent changes.
In order to understand the aforementioned changes in the SPNA on multi-decadal timescales many authors have turned to indirect measurements of the AMOC. One particular proxy of AMOC strength that has received some focus recently are 60 density anomalies at depth in the western SPNA or Labrador Sea region. In climate models, density anomalies in the western SPNA are a key predictor of density anomalies further south on the western boundary, and hence of the AMOC strength via thermal wind balance (Hodson and Sutton, 2012;Ortega et al., 2017;Robson et al., 2014Robson et al., , 2016. Observations show considerable decadal variability in subsurface density anomalies; density anomalies in the western SPNA or Labrador Sea between ~1000-2500 m increased significantly and peaked in ~1995 and subsequently declined (Robson et al., 2016; Although there is consistency across proxies of AMOC changes in the North Atlantic, there are considerable gaps in our understanding and major uncertainties to overcome. For example, the development of the subsurface density proxies has been investigated so far in just a few models (Ortega et al., 2017;Robson et al., 2014). However, there is considerable spread across climate models in the simulations of AMOC mean state and variability (Reintges et al., 2017;Zhang and Wang, 2013), and also in the latitudinal coherence of AMOC anomalies (Li et al., 2019;Roberts et al., 2020;Hirschi et al., 2020), which might 80 reflect different roles of deep density anomalies in the western SPNA on the AMOC, as well as different interplays between the subpolar and subtropical gyre contributions (Zou et al., 2020). Models also do not resolve realistically many key features of AMOC, most notably the overflows, and this affects the subsurface stratification downstream and on the western boundary (Zhang et al., 2011). There also remains significant uncertainty for other important processes. For example, it is not yet clear whether the recent changes in the SPNA are an ocean response to buoyancy forcing, or whether mechanical wind forcing has 85 shaped the recent observed changes (Robson et al, 2016;Piecuch et al. 2017). Local surface fluxes are also likely to explain a significant proportion of the recent cooling . Subsurface density anomalies are not just a proxy for the AMOC, but more generally for buoyancy forced (or thermohaline) circulation changes, including gyre changes (Ortega et al., 2017;Yeager, 2015). Finally, the AMOC variability is also thought to respond to local wind forcing on a range of timescales, especially at lower latitudes (Polo et al., 2014;Zhao and Johns, 2014), which could disrupt or "mask" the influence of 90 subsurface density anomalies as they propagate further south.
There is also considerable uncertainty in how and where subsurface density anomalies are formed in the SPNA, and how they are related to the AMOC. In observations and models, most water transformation associated with the AMOC occurs within the SPNA, and particularly in the eastern SPNA (Desbruyères et al., 2019;Grist et al., 2014;Langehaug et al., 2012).

95
However, subsurface density anomalies in the western SPNA on decadal timescales have often been linked with buoyancy forcing and changes in deep convection in the Labrador Sea or with changes in the volume of Labrador Sea Water production (Yashayaev and Loder, 2016;Yeager and Danabasoglu, 2014). Many studies have also reported that the basin-wide AMOC in ocean-only and coupled models is sensitive to heat flux or buoyancy forcing in the Labrador Sea (Kim et al., 2019;Ortega et al., 2011Ortega et al., , 2017Xu et al., 2019;Yeager and Danabasoglu, 2014). However, the link between deep convection, deep water 100 formation, and density anomalies at depth in the Labrador Sea is complex, and not fully understood (Katsman et al., 2018).
Observations suggest that very little water transformation and deep water formation actually occurs in the Labrador Sea (Pickart and Spall 2007;Lozier et al. 2019). Indeed, recently it has been shown that the Labrador Sea (i.e. OSNAP-west) played very little role in the interannual variability so far observed across the whole OSNAP line (Lozier et al., 2019).
Moreover, ocean-only models appear to significantly overestimate the amount of deep water formed within the Labrador Sea, 105 with likely implications for coupled models (Li et al., 2019). These inconsistencies raise the question of whether models are simulating the right relationships.
In this study we will address some of the above uncertainties by performing a multi-model analysis of the North Atlantic in coupled climate models. We focus on the question of how robust is the relationship between subsurface Labrador Sea density 110 anomalies and the basin-wide Atlantic Ocean circulation on decadal timescales. We also assess the question of whether Labrador Sea density can robustly induce density changes over the western continental slope and generate a geostrophic response in the meridional circulation (Bingham and Hughes, 2009; the basin-wide AMOC cell, as well as to identify the models that can produce more reliable predictions and projections of the 115 SPNA. For this, we will assess specifically the connection between subsurface density and AMOC at high and low latitudes via the western boundary. Furthermore, we will determine whether models consistently support an impact of AMOC changes on the SPNA upper ocean temperatures, and if not, investigate why. Our primary aim is to provide, for the first time in a multi-model context, a broad characterization of these relationships using consistent analysis frameworks and tools, and to document the uncertainty. The reasons for the uncertainty in the relationships will also be explored, establishing links with key 120 model climatological properties that could eventually be exploited as emergent constraints. We intentionally do not explore in detail how subsurface density anomalies are formed in these models, and leave this for further study.
The paper is organised as follows. Section 2 describes the experiments and methods. Labrador Sea density, and its link with the ocean circulation and the wider North Atlantic are explored across the multi-model ensemble in Section 3. The

125
characteristics of the intermodel spread in the previous relationships are explored in Section 4, and Section 5 presents the main conclusions of this study and discusses its implications.

Experiments and methods
Here we provide an overview and brief description of the models used in this study and provide some statistical considerations 130 for the intermodel comparison.

Experiment selection
For the multimodel analysis, we use the preindustrial control simulations (picontrol) from the fifth phase of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012). In particular, we only use those models in which 3D fields of 135 ocean temperature and salinity, as well as the streamfunctions of meridional overturning circulation and/or the barotropic circulation, were available. Twenty different models meet this condition, and their main characteristics have been summarized in  Robson et al. (2016) and Ortega et al. (2017). Note that we will assume that the present day control in HiGEM can be compared with the other preindustrial simulations due to the large uncertainty these later show in their 150 climatological biases, and so, for the sake of simplicity, we will only refer to preindustrial control experiments from now on. Figure 1 demonstrates that this assumption is reasonable, since the mean Labrador Sea stratification in HiGEM is very similar to that in the other models.
As an observationally-constrained reference, this study also includes the assimilation run from DePreSys3, a decadal 155 prediction system from the MetOffice based on GC2 (Dunstone et al., 2016). In the ocean, the assimilation is performed through a strong nudging (ten-day relaxation timescale) towards the full fields of a three-dimensional objective temperature and salinity analysis (Smith and Murphy, 2007). Since it covers a comparatively shorter period , and therefore different timescales than the control experiments, its comparison with the other simulations will be done with caution, in particular regarding the indices of the large-scale Atlantic circulation, for which other assimilation products show important 160 discrepancies (Karspeck et al., 2015), thus highlighting significant uncertainty. Statistical significance of correlation coefficients is assessed following a two-tailed Student's t-test that takes into account the 170 series' autocorrelation to correct the sample size, reducing the degrees of freedom of a series to its effective value (Bretherton et al., 1999).

Because our goal is to provide further insight into the suggested relationships established from observed trends in the North
Atlantic (e.g., Robson et al., 2016), all statistical analyses in this study exploring the relationships between variables and associated lags are based on 10-year running trends. This is analogous to the calculation of a typical 10-year running mean, but 175 computing over each 10 year period a linear trend instead and keeping the slope value. Note also that our main results remain similar if decadal running means are applied instead (not shown), as both are alternative approaches to concentrate on the lowfrequency variability. Running trends have also the particular advantage of not being sensitive to long-term drifts, which are still present (and can be important for some simulations and variables) when running means are computed. To illustrate how decadal running trends represent low-frequency variability, these have been included in Figure 2b for an index of Labrador 180 Sea density.

Labrador Sea density as an index of multi-decadal North Atlantic variability
This section explores the potential of Labrador Sea density as a proxy of the ocean circulation changes in the North Atlantic.
As in our previous studies (Ortega et al., 2017;Robson et al., 2016), the indices that we will herein define represent waters within the Labrador Sea and not those that are necessarily formed in the region (e.g. Labrador Sea Water). Since Labrador Sea 185 variability is affected by different processes (e.g. vertical mixing, Arctic-Atlantic overflows, sea ice interactions) that can be represented differently in the models, both in time and space, we characterize its variability over a relatively broad box (60°W-35°W; 50°N-65°N, blue box Figure 1a) that also includes part of the Irminger Sea region. Note that over this large area EN4 (Good et al., 2013), an objective analysis of monthly temperature and salinity 3D observations, shows the weakest density stratification in the North Atlantic (characterised in Figure 1a as the density difference between 1000m and the surface).

Labrador Sea density across models
A first indicator of potential model discrepancies is Labrador Sea stratification, which can lead to differences in the representation of deep ocean convection (i.e. weaker density stratifications will facilitate the mixing, fostering convection activity, and vice versa for stronger density stratifications). temperature. Most models present their warmest waters at the surface, and temperatures decrease sharply to minimum values around 100 m and increase again at deeper levels, reaching uniform conditions after approx 300 m. However, the location and magnitude of this temperature minimum and the two maxima are highly variable. It is important to note that the profile for one of the models, MRI-CGCM3, is noticeably different to the others, with a subsurface minimum more than 2 degrees colder than for any of the other models. In terms of salinity, the general profile is more coherent across models, with minimum 200 salinity at the surface that progressively increases with depth and attains uniform values after 500 m. Density stratification seems to be determined by salinity, as their two vertical profiles show similar features. This similarity includes exceptionally strong density and salinity stratification in MRI-CGCM3 as compared with the other models. This stratification is so strong that it precludes the occurrence of deep convection (not shown). Because of this, MRI-CGCM3 is an outlier for many of the metrics used in the paper, and has been excluded for the subsequent analyses to facilitate the interpretation of our results. We 205 also note that the profiles for the two eddy-permitting models (green and orange lines in Figure 1b

Labrador Sea density linkages with the ocean circulation
The link between PC1-LSD and other ocean circulation indices in the North Atlantic is now examined. Three indices are  there is some inter-model spread regarding the lag of maximum correlation, which ranges between 0 and 2 years (with PC1-LSD leading), although both variables are in phase for the majority of models. The AMOC at 26°N (AMOC26) is also positively related to PC1-LSD, with PC1-LSD leading AMOC26 by three years on average (Figure 4b). However, the average correlation between PC1-LSD and AMOC26 is weaker, and the spread in the magnitude and lag of the maximum correlation is larger than for AMOC45. Therefore, it appears that the link with the subtropics is weaker than for 45°N and that AMOC 285 coherence between subpolar latitudes and the subtropics in coupled models is model dependent. The reasons for the spread in the relationship between PC1-LSD and AMOC26 are explored in Section 4. A strong relationship is also found between PC1-LSD trends and those in SPGSI (Figure 4c), of similar order than for AMOC45. Thus, overall, PC1-LSD is a good proxy for https://doi.org/10.5194/esd-2020-83 Preprint. Discussion started: 5 November 2020 c Author(s) 2020. CC BY 4.0 License.
the large-scale ocean circulation in the Subpolar North Atlantic, and can also be a precursor for a fraction of the AMOC variability in the subtropical Atlantic. This result is further supported by a parallel analysis in Figure 5, looking at the 290 maximum correlation between the decadal AMOC trends and those in Labrador Sea density as a function of depth, when the latter leads the AMOC by up to 10 years. Figure 5 reveals that the strongest link between the Labrador Sea densities and the AMOC, both at 45 and 26°N, occurs in its first 1000 m, the same levels where the first EOF of LSD show the maximum loadings (Figure 2a), which confirms the appropriateness of using PC1-LSD to represent the ocean circulation. The same analysis also supports a strong link between SPGSI and LSD, although in that case the largest correlations usually happen at 295 deeper levels (between 1000 and 2000 m). Note also that the main conclusions drawn from PC1-LSD are also valid for the dLSD index: however, the inter-model differences are larger (Supplementary Figure 1).

300
Significance is assessed as in Figure 2d and indicated with a circle. For positive lags, PC1-LSD leads. b-d The same as in a but between PC1-LSD and the maximum AMOC streamfunction at 26°N after the Ekman transport is removed (AMOC26), the subpolar gyre strength index (SPGSI) and the vertically averaged top 700 m temperatures in the eastern subpolar gyre (ESPNA-T700; grey box in Figure 1a).

Figure 5: a Maximum correlation (at any lag) between the AMOC45 index (after the Ekman transport is removed) and
Labrador Sea Densities as a function of depth for all the simulations. Colored dots indicate correlations that are significant at the 95% confidence level. b-c The same as in a but between the AMOC26 index and LSD, and between the SPGSI and LSD, respectively.

310
Previous studies based on the GC2 picontrol simulation have suggested LSD to be also a potential predictor of wide-spread  (Figure 6a) show that there is a coherent relationship between both variables across models, with PC1-LSD increases (decreases) being consistently followed by ESPNA-T700 warmings (coolings).
Nevertheless, there are inter-model differences concerning the magnitude and lag of the strongest positive correlations, revealing important uncertainty in the relationship. The spread in the PC1-LSD vs ESPNA-T700 relationship is thus reminiscent of the spread found between PC1-LSD and AMOC26, which suggests that they might be related. We also note 320 significant negative correlations when ESPNA-T700 leads PC1-LSD by 2-4 years that might be explained by the opposed (and nearly concomitant) impacts that the NAO exerts on both variables (Figure 6b,c). Positive NAO phases, and associated surface buoyancy forcing (Lozier et al., 2008) lead in first instance to negative SST (Barrier et al., 2014;Lohmann et al., 2009) and an almost simultaneous cooling in ESPNA-T700 (Figure 6b). In comparison, on the western side of the SPNA, positive NAO phases contribute to reduce vertical density stratification, favoring convection and a more positive LSD index (Robson 325 et al., 2016), which in the models lags the NAO by 2-3 years (Figure 6c). The fact that correlations between NAO and ESPNA-T700 are weaker than between PC1-LSD and ESPNA-T700 suggests that the ocean might also be playing an additional role.
https://doi.org/10.5194/esd-2020-83 Preprint. Discussion started: 5 November 2020 c Author(s) 2020. CC BY 4.0 License.  Figure 1a). Correlations are based on 10-year running trends. Significance is assessed as in Figure 2d and indicated with a circle. For positive lags, PC1-LSD leads. b-c The same as in a but between the North Atlantic Oscillation (NAO; defined as the standardised difference in sea level pressure between the closest grid-points to Azores and Reykjavik) and the ESPNA-T700, and between the NAO and the PC1-LSD, respectively. In these two cases, for negative lags the NAO leads.

335
The link between PC1-LSD and the ESPNA could be explained through an influence of the first on the meridional ocean heat transport. This link is now investigated in the two eddy-permitting simulations (Figure 7) and in the five CMIP5 models for which the ocean heat transport fields are publicly available. In the two high resolution experiments and two of the CMIP5 ones the decadal trends in the meridional ocean heat transport at 45°N (OHT45) are strongly linked with those in PC1-LSD. This is 340 a similar relationship to the one previously found in Figure 4 between PC1-LSD and both the AMOC45 and SPGSI, but in this case with PC1-LSD leading with slightly longer lead time. The other CMIP5 experiments support a weaker, yet significant, link, as well as a longer lag between OHT45 and PC1-LSD. Altogether, Figure 7a confirms that PC1-LSD is a good precursor of the changes in the meridional ocean heat transport, although with some differences across models which might reflect a different representation of certain processes. The contributions of two different processes to this delay are further investigated 345 in HiGEM, for which OHT had been decomposed online at each time-step into vertical and horizontal heat transports (as in Bryan, 1969), which can be respectively interpreted as the "overturning" (i.e. characterised by the zonal mean transport) and "gyre" (i.e. characterised by variations from the zonal mean transport) components (Robson et al., 2018a). While the overturning contribution (OHT45 over ) increases in phase with the AMOC45, SPGSI and PC1-LSD changes (Figure 7b

Characteristics of the inter-model spread in subpolar to subtropical AMOC
This section investigates which particular climatological model features are linked to the large inter-model spread in the PC1-LSD vs AMOC26 relationships. The most relevant model features thus identified will improve our process understanding, and 360 can eventually be used to identify which models are most realistic and, in turn, can deliver more reliable projections of the future changes in the North Atlantic. Figure 8 shows that models that simulate a stronger and deeper climatological AMOC (both at 45°N and 26°N) tend to have a stronger correlation between PC1-LSD and the subtropics. All these linear relationships between climatological AMOC strength and depth and the PC1-LSD vs AMOC26 connectivity are significant at the 95% confidence level. These

365
climatological AMOC values (without Ekman) can be put in context with those from RAPID observations and DPS3. RAPID observational uncertainties have been considered by including the mean values over three different non-overlapping periods (i.e. 2004-2007, 2008-2012 and 2013-2016; dotted lines in Figure 8). The scatterplots show that the majority of models whose climatological AMOC26 lies within the RAPID/DPS3 climatological spread have a relatively weak link between PC1-LSD and AMOC26, although some models supporting a strong link are also included or remain close to the RAPID/DPS3 values.

370
However, caution is recommended, e.g., before defining emerging constraints, because model and observations are not directly comparable for numerous reasons: both RAPID and DPS3 cover shorter periods than the simulations and relate to different background forcing conditions (present day vs preindustrial) which might imply different mean states (Thornalley et al., 2018). Furthermore, evidence suggests RAPID calculations from mooring arrays might be underestimating the AMOC strength by ~1.5 Sv .

385
A potentially important factor behind the inter-model spread in Figure 4b is the mean density stratification in the Labrador Sea. Figure 9 suggests that, indeed, the PC1-LSD vs AMOC26 spread is partly influenced by the density stratification in this region. Models with a weaker density stratification have a stronger relationship between PC1-LSD and AMOC26. However, no link is found between the PC1-LSD vs AMOC26 relationship and temperature or salinity stratification in the Labrador Sea.

390
Models that have a weaker density stratification (here defined as the difference between the top 100 m, and the average between 500-1000 m), and thus favor deeper convection in the Labrador Sea, generally exhibit a stronger link between PC1-LSD and AMOC26. This result is robust for other stratification indices based on different depth levels (See Supplementary   Figure 2). It is also worth mentioning that all models except CanESM2 are more weakly stratified in the Labrador Sea than the observations (represented herein by the DePreSys3 assimilation run and EN4). Hence, the real link of LSD with the AMOC26 395 may not be as strong as some models suggest.
https://doi.org/10.5194/esd-2020-83 Preprint. Discussion started: 5 November 2020 c Author(s) 2020. CC BY 4.0 License. correlations are based on 10-year running trends. The correlation coefficient between the two metrics is shown in the top-left corner. The presence of an asterisk indicates that the correlation is significant at the 95% confidence level. Colors indicate the lag at which the maximum correlation between PC1-LSD and AMOC26 is obtained. The grey (green) vertical lines depict the mean stratification value in the DePreSys3 assimilation run EN4. In both cases, their overlap period is used to compute the climatology (i.e., 1960-2013). b-c The same as in a but for the Labrador Sea salinity and density, respectively.

405
Another key aspect of the LSD vs AMOC26 connectivity is the western boundary density. Indeed, boundary density is critical to the mechanism through which LSD influences the AMOC at lower latitudes. Positive (negative) LSD anomalies propagate equatorward following this boundary, and as they do so they strengthen (weaken) the zonal density gradient, triggering a thermal wind response that accelerates (decelerates) the AMOC. In the following we investigate differences in the propagation 410 of boundary densities across models, and if these differences can affect the inter-model PC1-LSD vs AMOC26 spread. Figure   10 focuses on the two high-resolution simulations, where important differences already manifest. It represents the correlations of PC1-LSD with the density fields near the western boundary at four different longitudinal transects: 57°N (cutting across the Labrador Sea), 45°N, 35°N and 26°N. In both models, the depth of the maximum correlation near the continental shelf is coherent across latitudes. However, in HiGEM these occur at deeper levels (1000 to 3000 m) compared to GC2 (1000 to 2000 415 m), and the difference is especially clear at 35°N, where the highest correlations occur at ~2000 m in HiGEM, while only at 1000 m in GC2. Similar depth differences are also found at 26°N, but with slightly weaker correlations. In addition to the difference in the depth of the maximum correlation between HiGEM and GC2, there are differences in the vertical structure between the two models. For example, at 35°N in GC2, density anomalies on the western boundary form a tripole (low correlation above and below the maximum correlation at ~1000 m), but in HiGEM the density anomalies form a dipole 420 ( Figure 10f). We note some differences in bathymetry at this latitude (which is steeper in HiGEM), which might partly explain some of the differences in terms of the density correlation structure.
https://doi.org/10.5194/esd-2020-83 Preprint. Discussion started: 5 November 2020 c Author(s) 2020. CC BY 4.0 License. 430 Figure 11 shows that the diversity in the depth of these boundary densities is even more evident when including the CMIP5 models. The depth of the maximum correlation between PC1-LSD and the western boundary density at the four latitudinal sections relates linearly (and significantly at the 95% confidence level) across models with their PC1-LSD vs AMOC26 correlation. In this case, models exhibiting deeper WBDs generally show stronger links between LSD and the subtropical AMOC. We also checked if models with stronger correlations at the WBDs also support a stronger link between the LSD and 435 the AMOC (correlations in magenta in Figure 11), but this only holds true at 57°N. This suggests that the depth along which WBDs propagate southward, and/or the vertical structure of anomalies, are the key aspects to understand and potentially narrow down the spread.

Conclusions and discussion
This article has explored, in a multi-model context, the linkages between subsurface density in the subpolar North Atlantic

465
-The connectivity between anomalies in LSD and AMOC at 26°N is sensitive to different model features, including the strength and depth of the climatological AMOC maximum, the mean density stratification in the Labrador Sea, and the depths at which the LSD propagates southward along the western boundary. Stronger LSD connectivity with the subtropics tends to occur in models with a stronger and deeper AMOC, weaker Labrador Sea stratification and western boundary density propagating at deeper levels.
-Observationally derived constraints of the model based relationships tend to suggest that the link between LSD and the subtropical AMOC is weak. This suggests that observations of AMOC via RAPID may not be representative of the basin wide buoyancy forced AMOC variability. However, caution is advised because simulations and observations are not directly comparable, and so significant uncertainty remains in constraining the relationship between LSD and subtropical AMOC.

475
-The multi-model ensemble does also support a significant lagged relationship between LSD and the upper ocean temperature in the eastern SPNA, in line with previous studies linking LSD to the recently observed changes in the North Atlantic. However, models disagree regarding the strength of the link (correlations between 0.3-0.7), and the maximum lag (3 to 10 years).

480
We have shown that, in coupled climate models at least, subsurface density anomalies in the wstern SPNA are an important predictor of the wider North Atlantic ocean circulation and upper ocean temperature in the SPNA. This importance on the ocean circulation is especially clear at the latitudes of the SPNA itself. Given the important role of the wind in driving lower latitude AMOC anomalies, and the range of processes by which wind can act on the AMOC (Duchez et al., 2014b(Duchez et al., , 2014aKanzow et al., 2010;Polo et al., 2014;Zhao and Johns, 2014) it is not surprising that the relationship between LSD and 485 AMOC at 26°N is much weaker. Nevertheless, the reasons behind the large spread in these relationships across models is not so clear.
We have tried to constrain this uncertainty by looking at a range of observed metrics that may explain the spread in the correlation strength, including the density anomalies on the western boundary, the stratification of the Labrador Sea, and the mean-AMOC strength. Overall, these constraints point to a relatively weak relationship between LSD and AMOC at 26°N on 490 decadal timescales (i.e. r ~ 0.4) in the real world. However, there are many reasons why this number is still very uncertain.
The simulations and observation-based datasets are not directly comparable, as they differ in the background radiative forcing levels, the length of the period used to compute the climatologies, and even the way some indices, like the AMOC, are computed. There is also emerging evidence that models underestimate AMOC variability on decadal timescales (Roberts et al., 2013;Yan et al., 2018). We also recognise that there is large uncertainty within the observationally derived metrics. For 495 example, the assimilation run in DePreSys3, which is used to constrain relationships clearly does not represent the mean AMOC strength at 26°N from RAPID (see Figure 8b). Models generally tend to underestimate the depth of the return flow, and this may still affect how density anomalies project on the basin-wide AMOC. It has also been argued that ocean-only models produce too much deep water in the western basin and Labrador Sea (i.e., Li et al., 2019), and recent observations even challenge the prevailing view from models that Labrador Sea convection dominates the AMOC variability (Koenigk and 500 Brodeau, 2017), suggesting that the key deep water formation occurs North East of the Labrador Sea (Lozier et al., 2019).
Therefore, further in-depth study is required to understand the uncertainty in the relationships.
Most of the models considered in this study have relatively coarse resolution, including non-eddying oceans (≥ 1°x1°), which means that they might be missing some key dynamics for the AMOC (Johnson et al., 2019) that could be important to represent realistic linkages. The current analysis, however, includes two models at eddy-permitting resolution (HadGEM3-505 GC2 and HiGEM), whose relationships lie within the spread of those in the coarser models, suggesting that climatological features (like Labrador Sea stratification) can be more important than the representation of mesoscale processes. However, it could be that higher resolution is needed (e.g. enabling meso-scale eddies in subpolar latitudes) to identify substantial https://doi.org/10.5194/esd-2020-83 Preprint. Discussion started: 5 November 2020 c Author(s) 2020. CC BY 4.0 License.
differences (Hirschi et al., 2020;Johnson et al., 2019). A recent analysis based on HadGEM3-GC3.1 (a later version of HadGEM3-GC2) configured at different horizontal resolutions has shown that long-standing model biases, including in the 510 North Atlantic, are reduced at eddy-resolving resolution (1/12° x 1/12° in the ocean), and that the strength of the AMOC, the boundary currents and the northward heat transport is higher than for the coarser resolutions (Hirschi et al., 2020;Roberts et al., 2019). High resolution coupled models also generally support the new view from OSNAP observations in which the largest fraction of AMOC variability (on sub annual to decadal timescales) originates at the eastern SPNA (Hirschi et al., 2020). Eddy-resolving resolutions have also been shown in a multi-model study (Roberts et al., 2020) to represent the AMOC 515 response at 26°N differently in future projections, leading to stronger declines than in non-eddying simulations, declines mostly associated with a weakening in the Florida Current. Roberts et al. (2020) also compares the meridional coherence of the AMOC, which does not seem to be resolution-dependent, a result that is in line with another multi-model comparison between non-eddying and eddy-permitting simulations (Li et al., 2019).
Despite the current limitations in the models considered for this study, it is important to highlight that they provide a rather 520 consistent picture of a chain of relationships in the North Atlantic that is able to explain some of the recent observed trends (Robson et al. 2016). This paper has broadly characterized this behaviour, and highlighted the uncertainty. These relationships are also consistent with the mechanisms proposed by Yeager and Robson (2017) to explain high levels of predictive skill in the SPNA on decadal timescales. Our analysis has also helped to identify specific metrics (such as LSD stratification and the depth of the boundary density) that could be used as emergent constraints for future projections, i.e. to subset the simulations 525 expected to more realistically represent the future changes in the region. Having a more realistic subpolar gyre stratification at present day conditions has been shown in CMIP5 simulations to increase the probability of a future collapse in convection (Sgubin et al., 2017), that would lead to a widespread SPG cooling. It remains to be tested if similar conclusions can be drawn from eddy-resolving simulations.

530
Code availability. The main scripts used in the analysis and other supporting information that may be useful to reproduce the results of this article are archived at the Barcelona Supercomputing Center and will be shared upon request by the corresponding author.
Author contributions. P. O., J. R. and R. S. conceived the study, which was later discussed and refined with the other coauthors. M. M. downloaded and processed the CMIP5 data, computing the main climate indices. P. O. led the analysis, and together with J. R. prepared the manuscript with contributions from all co-authors.

540
Competing interests. The authors declare that they have no conflict of interest.