Recent subsurface North Atlantic cooling trend in context of Atlantic decadal-to-multidecadal variability

Abstract The spatiotemporal structure of the recent decadal subsurface cooling trend in the North Atlantic Ocean is analyzed in the context of the phase reversal of Atlantic multidecadal variability. A vertically integrated ocean heat content (HC) Atlantic Multidecadal Oscillation index (AMO-HC) definition is proposed in order to capture the thermal state of the ocean and not just that of the surface as in the canonical AMO SST-based indices. The AMO-HC (5–657 m) index (defined over the area, 0°N–60°N, 5°W–75°W) indicates that: (1) The traditional surface AMO index lags the heat content (and subsurface temperatures) with the leading time being latitude-dependent. (2) The North Atlantic subsurface was in a warming trend since the mid-1980s to the mid-2000s, a feature that was also present at the surface with a lag of 3 years. (3) The North Atlantic subsurface is in a cooling trend since the mid-2000s with significant implications for predicting future North Atlantic climate. The spatial structure of decadal trends in upper-ocean heat content (5–657 m) in the North Atlantic prior to and after 2006 suggests a link with variability of the Gulf Stream–Subpolar Gyre system. The Gulf Stream leads variability in upper-ocean heat content over the Subpolar region by ∼13 years, and the lead increases from east to west from the Iceland Basin to the Irminger Sea to the Labrador Sea. The similarity between the structure of decadal mean anomalies, their change and trends in upper-ocean heat content and salinity in the North Atlantic and anomalies associated with a Gulf Stream index is striking. In this scheme, the displacements of the Gulf Stream–North Atlantic Current systems and their interactions with the Subpolar Gyre, as part of the meridional overturning circulation, have a decisive role for imposing decadal variability in both the Ocean and hydroclimate over the neighbouring continents.


Introduction
The state of the North Atlantic Ocean is a key element in understanding regional and global climate variability and change. Reversals of SST anomalies over the North Atlantic have decadal-to-multidecadal variability, which, in turn, imposes warming and cooling decadal trends when transitioning from one phase to the other. Basinwide trends have traditionally been identified via the Atlantic Multidecadal Oscillation (AMO) index (Enfield et al., 2001;Guan and Nigam, 2009), an SST-based index which integrates the vagaries of the atmosphere and as such tends toward the low-frequencies (e.g. Frankignoul and Hasselmann, 1977;Deser et al., 2010). However, the characteristic features of the Subpolar Gyre (SPG) and Subtropical Gyre, warm and cold currents (Reverdin et al. 2003), including the Gulf Stream (GS) system, and their advective properties make the region more than a passive responder to atmospheric variability (Bjerknes, 1964;Bryan and Stouffer, 1991;Wang et al., 2004;Reverdin, 2010;Zhang et al., 2016;Delworth et al., 2017;Zhang, 2017;Nigam et al., 2018).
Modelling and empirical analyses indicate that heat transports by the Gulf Stream-North Atlantic Current (GS-NAC) system are important in order to have lowfrequency variability in the North Atlantic (Sutton and Allen, 1997;Delworth et al., 2017;Zhang, 2017). Atlantic multidecadal variability requires these transports to be properly accounted for (Drews and Greatbatch, 2016).
Not surprisingly the largest year-to-decadal variability occurs along this current system (Buckley and Marshall, 2016). Extreme climate episodes in the North Atlantic at decadal scales are characterized by an out-of-phase relationship of the oceanic conditions between the Subpolar Gyre and GS regions (H€ akkinen, 2001;Zhang, 2008;H€ akkinen and Rhines, 2009;H€ akkinen et al., 2013;Zhang and Zhang, 2015). As such, cooling trends in the upper North Atlantic Ocean such as those after the 1970s (Rahmstorf et al., 2015) and after 2005 have recently been related to a weakening of the Atlantic Meridional Overturning Circulation (AMOC) by Robson et al. (2016). The latter investigators found that the density anomalies became negative due to the pronounced 1995-2005 warming of the subpolar gyre, which, in turn, decreased the strength of the formation of North Atlantic deep waters and ultimately lead to a negative anomaly in the meridional heat transport and hence a cooling. This buoyancydriven AMOC response is likely to also have triggered the cooling of the 1960s (Robson et al., 2014). Most importantly, Robson et al. (2016) concluded by examining trends in atmosphere-ocean surface fluxes and SSTs that the warming and cooling trends observed before and after 2005 were of oceanic rather than of atmospheric nature. There is, however, yet much to learn not only about this newly reported decadal cooling trend of the Subpolar North Atlantic (SPNA) but also about how regional decadal trends are produced, regarding the role of propagating subsurface anomalies as well as their link with the GS, North Atlantic current (NAC) and SPG system.
In the present paper, the surface and subsurface thermal structures of the North Atlantic basin, especially its decadal-to-multidecadal variability are closely examined to provide some insights on the role of the GS-NAC-SPG system in the generation of the recent decadal cooling. The thermal state of the North Atlantic surface waters is characterized here by using the areal average of the seasonal SST anomalies in the 75 W-5 W and Eq-60 N domain, viz. the AMO index. This index is, typically, linearly detrended and smoothed to highlight decadal-multidecadal variability; the spatiotemporal evolution of the AMO, identified as a mode of SST variability in the Atlantic, was characterized in Guan and Nigam (2009), and its subsurface structure explored in Kavvada et al. (2013). There is an interplay in upper-ocean heat content anomalies between the SPG and GS regions which changes with time before and after the mature phase. In particular, under the mature positive phase of the AMO, the subpolar region warms while the GS path region, along the North American coast, cools ( Fig. 7 in Kavvada et al., 2013;Wang and Zhang, 2013;Zhang and Zhang, 2015). Similarly to the surface AMO index, and following Kavvada (2014), the subsurface thermal state is characterized using the same areal averaging but of linearly detrended ocean heat content (HC) over various column depths; hereafter, the AMO-HC index. The ocean heat content is an informative variable reflecting the influence of both ocean-atmosphere fluxes and ocean dynamics (H€ akkinen et al., 2015Cheng et al., 2017).
The paper is divided as follows. Section 2 gives information on the data sets and methods used. Section 3 has three subsections for: (1) analyzing the HC-based AMO indices and their robustness as well as the relationship between the AMO and subsurface temperatures; (2) examining the spatial structure of decadal trends before and after 2006 in the context of decadal trends and anomalies, and their temporal changes in structure with what takes place in the SPNA; (3) investigating variability in the SPNA via HC-based indices and the relationship between the subsurface structures of this ocean region and GS variability. Finally, Section 4 provides some concluding remarks.

Data sets and methods
The present study is based on the analysis of different data sets at seasonal resolution for the period 1950-2015. Unless noted differently, anomalies are calculated with respect to the long-term 1950-2015 climatology. Analysis techniques in the present study use standard statistical tools frequently employed in climate studies such as lead/ lag correlations and regressions, and detection/extraction of linear trends via the least squares method. Statistical significance of the regressions and correlations, as well as of the linear trends, is assessed by a 2-tailed Student's t-test at the 95% confidence level using an effective sample size that accounts for serial correlation (Appendix A in Metz, 1991;Santer et al., 2000); the significant regressed anomalies, correlations and trends are stippled. The robustness of the results is also investigated by using different ocean-data sets. Sea surface temperatures (SST) are from two sources: (1) the UK's Met Office HadISST, version 1.1, data set (Rayner et al., 2003), (2) NOAA's ERv3b data set (Smith et al., 2008). Subsurface ocean temperatures used to calculate HC come from three different sources: (1) The UK Met Office EN.4.2.0 quality-controlled objective analysis data set (Gouretski and Reseghetti, 2010;Good et al., 2013), (2) The Japanese JAMSTEC/Ishii and Kimoto, version 6.12, data set (Ishii and Kimoto, 2009), (3) NOAA's NODC data set (Boyer et al., 2005) identified as NODC. The EN4.2.0 data set has 42 depth levels from 5 to 5000 m for the period 1900-present, the Ishii subsurface has 24 levels from 0 to 1500m, and the subsurface NODC data set is given at 26 vertical levels from 0 to 700 m for the period 1955 À 2010. To characterize the regional variability, the following indices are constructed: (1) The traditional AMO index, based on area-averaged SST anomalies over the Atlantic Ocean north of the Equator (75 W-5 W, Eq-60 N), which are linearly detrended via the leastsquares method and smoothed with the LOESS filter (Cleveland and Loader, 1996) using 25% of the data to highlight decadal-to-multidecadal variability.
(2) Alternative AMO indices based on linearly detrended vertically integrated HC anomalies (AMO-HC, Kavvada, 2014) at several depths. The AMO and AMO-HC indices, as defined above (detrended and smoothed), are used to calculate lead/lag correlations. AMO and AMO-HC indices with the linear trend included are also calculated and contrasted in some figures.
Temperature and salinity are used to calculate potential density over the subpolar region from the EN4.2.0 data set following the UNESCO formula (e.g. Gill, 1982, appendix three). By using climatological values of one of the two variables, one can assess the effect of the other variable for generating density anomalies. That is, by using climatological salinity (temperature) while using the anomalous temperatures (salinity) in the formula, one can assess the effect of the anomalous temperature (salinity) on the anomalous density.
In addition to the HC/SST indices mentioned above, a Gulf Stream (GS) index is also used to characterize subtropical North Atlantic variability. This index (kindly provided by T. Joyce, personal communication) is based on an EOF analysis of the 15 C isotherm at 200 m depth at selected locations between 75 W-50 W and 33 N-43 N; this isotherm represents an isotherm in the centre of the strong horizontal gradient of the Gulf Stream, just north of the maximum flow at the surface, and is considered a convenient marker for the northern 'wall' of the stream (Joyce et al., 2000). The GS index was provided as a smoothed standardized index at seasonal resolution for the period 1954-2012 and is also smoothed with the LOESS filter using 25% of the data. Because subsurface data are sparse at depths beyond 2000 m, the present analysis is largely focused on the upper layers of the ocean.
Also, two other data sets are used to explore the role of atmosphere-ocean fluxes in the recent trends. One is the UK Met Office Hadley Centre's mean sea level pressure data set (HadSLP2 - Allan and Ansell, 2006) and available on a 5 Â 5 grid from 1850 to the present. The other is the Woods Hole Oceanographic Institution's Objectively Analyzed air-sea Fluxes (OAFlux) for the Global Oceans (Yu et al., 2008) data set and available on a 1 Â 1 grid from 1958 to September 2016; sensible and latent heat fluxes as well wind speeds at 10 m are used in the present work.
A 20-year mean of dynamic topography (Pujol et al., 2016) data from the AVISO is used to approximately characterize the mean position of the SPG, subtropical gyre and the Gulf Steam current by displaying the absolute dynamic topography isolines at À0.4, 0.4, and À0.1 m, respectively. This product was produced by Ssalto/Duacs and distributed by Aviso, with support from the French Centre National D'etudes Spatiales (http:// www.aviso.altimetry.fr/duacs/).

A heat content AMO index and current cooling trend
The AMO index (Fig. 1a) exhibits decadal-scale pulses, but shows a broad cooling trend from the 1950s to the mid-1970s, followed by a warming trend up to the mid-2000s, intercepted by interannual and decadal fluctuations. Comparison of the index using two different data sets (HadISST vs. ERv3b) shows negligible differences between them. The AMO index is also shown prior to removal of the 1950-2015 linear trend, for reference; in any case these indices show a cooling trend after 2008-2009. The evolution of undetrended AMO's surface (SST) and subsurface (HC) thermal state is shown in Fig. 1b. The AMO-HC index is shown for various depths, and a comparison of the surface and subsurface thermal evolution suggests that during the period of a notable surface cooling (from 1960 to the mid-1970s), the subsurface North Atlantic was already shifting to a warming phase since 1970, however, both the surface and subsurface were in phase during this period preceding the most pronounced warming in the 1990s. There is a relative out-ofphase relationship between the surface and the subsurface North Atlantic in the second half of the twentieth century until 1998: surface/subsurface normalized anomalies were warmer/cooler than the subsurface/surface normalized anomalies in the period $1950-71, while conditions reversed in the period in $1971-98. The subsurface warming that occurred from the late-1980s to the mid-2000s is followed by cooling, which is observed to start early in the subsurface layers than at the surface with a lead time of $2-3 years. Further analysis of the temporal lead-lag between AMO (surface) and subsurface temperatures ( Fig. 1c) indicates that near-surface temperature variations at depths of $50-200 m in the subpolar latitudes (45 N-60 N) take place almost simultaneously ($0-1 years) with the AMO, while subsurface temperature variations at depths $300-500 m in the subtropical latitudes (30 N-45 N) appear to lead the surface ones from the AMO. This is consistent with the analysis done by Wang et al. (2010) indicating that temperature and salinity anomalies at subpolar latitudes lag those from subtropical latitudes by 8-9 years and this happens because oceanic advection. In essence, this picture suggests that the classical surface (basin-wide) definition of the AMO follows closely the near-surface variability at subpolar latitudes, but masks out the apparent leading influence of the variability of the subsurface layers in the subtropics.
Low frequency variability at the subsurface in the North Atlantic, which seems to be leading surface variability, is captured by the HC-based AMO indices discussed above. A closer look at the structure associated with the Atlantic Multidecadal variability in the subsurface and its temporal evolution along the 35 W meridian is shown in Fig. 2 by regressing the new smoothed AMO-HC (5-657 m) index (dark blue line in Fig. 1b) on temperature and salinity detrended anomalies. The lead-lag figure shows upward propagation of cold/fresh subsurface anomalies from below 400 m in the 40 N-50 N band (perhaps along the climatological isotherms/isohalines), and deepening of the warm/saline near-surface (at $200 m depths and shallower) anomalies to $1200-1600 m in the 50 N-60 N band, which were present at the mature stage of the AMO-HC (t panel). Changes along the GS-NAC path are evident through the warm SST anomalies, and the temperature and salinity anomalies along the 40 N-50 N band from the depth-latitude profiles that, eventually, propagate along the NAC-SPG system; nearsurface anomalies change sign 8 years after the mature phase of the AMO-HC (t þ 8 yrs panelalso, Fig. 4.6 for HC (5-315 m) anomalies in Kavvada, 2014). Thus, it is reasonable to think that the cooling and warming trends at the surface and the subsurface are induced by ocean circulation (i.e. AMOC) and reflected by the Atlantic multidecadal variability (i.e. AMO). The structure of the trends and their subtropical to subpolar links will be explored in the next two subsections, but first the robustness of the basin-wide trends will be assessed.
Trends in the AMO-HC index are similar when obtained from three different data sets as Fig. 3 shows, which displays the AMO-HC (0-700 m) indices all relative to the 1981-2010 climatology. As ocean observations are sparse and less reliable at greater depths, especially in the early record, the variations in heat content of the upper ocean (0-700 m) are further analyzed using several data sets (EN4.2.0, NODC and Ishii whose objective analyses differ in details as well in their choices of temporal resolution). All data sets are in agreement, especially in the representation of the impressive subsurface warming from the early-1970s to the mid-2000s, and the notable cooling thereafter.
The AMO-HC indices readily identify the warming and cooling periods but do not provide any information on the regions within the North Atlantic basin where such trends are focused. The geographical distribution of the linear trend in upper HC for the period 1985-2005 (when the indices show a general warming) from the EN4.2.0 (5-657 m) and Ishii (0-700 m) data sets is shown in Fig. 4a, b. The largest warming trends from both data sets are in agreement, and they both highlight an almost general warming over the North Atlantic with the steepest trends in the SPG region. Also, consistent in both datasets is the presence of a cooling trend along the US east coast/GS path and regionally confined cooling trends to the east of the Flemish Cap at the western margin of the subtropical-subpolar gyre boundary. This latter region at the western margin of the subtropical and subpolar gyres is a dynamically complicated region (Marshall et al., 2001;Zhang, 2008;McCarthy et al., 2015;Buckley and Marshall, 2016). The latter paper referred to this as the transition zone, and argued to be a pacemaker for decadal AMOC variability.

Structure of decadal warming/cooling trends
Decadal trends in both upper HC and averaged salinity (5-657 m from EN4.2.0) show, before and after the subsurface maximum in 2006 (Fig. 5), distinctive features. The most pronounced warming and salinification is seen along the mean circulation during the warming period 1996-2005 (H€ akkinen et al., 2013). During the cooling period 2006-2015, however, it is the central SPG that shows the largest cooling and freshening. It is also apparent that the trends along the GS path close to the North American coast are much warmer and more saline during the case of the cooling rather than during the warming period, and that the NAC is farther to the northeast in the case of the 2006-2015 cooling period, than in the case of the 1996-2005 warming period (see e.g. H€ akkinen and Rhines, 2009;H€ akkinen et al., 2011). The latter study finds that during warm periods, the wind-driven subtropical and subpolar gyres are weaker than normal, the subpolar density surfaces relax, which, in turn, lead to a westward shift of the subpolar front that permits warmer waters from the subtropics to invade the eastern SPNA (Hatun et al., 2005). Conversely, the cooling period, as observed in Fig. 5 (right column), is associated with an eastward shift of the subpolar front as shown recently by Zunino et al. (2017), their Fig. 8 (see also Bersch et al., 2007; note that this study only used data until 2004 and hence the implied eastward shift of the subpolar front during the recent cooling period is analogous to that of the early 1990s). In fact, variability in SST, upper-ocean HC and salt content extends further to the northeast for the cooling period than for the warming period ( Supplementary Fig. 1S). The trend rates and spatial extent of the subsurface post-2006 North Atlantic cooling trend are larger than those of the pre-2006 warming trend. Note, however, that the warming trend in the  (Fig. 4). Similarly, the recent cooling trend is expected to last several more years (e.g. Yeager and Robson, 2017).
The structure of the cooling trend (Fig. 5) reflects the central role played by ocean dynamics rather than that caused by heat fluxes (Robson et al., 2016;Zhang et al., 2016;Delworth et al., 2017). The opposite-sign regions between the generalized cooling and freshening trends in the SPG region and warming and saline trends along the GS region are known to be a distinctive pattern of the AMOC, and the evolution of the reversed dipole pattern is thoroughly discussed in Zhang and Zhang (2015). The latter study shows that due to the existence of interior pathways, a warm anomaly propagating southward along the path of the deep western boundary current from the subpolar gyre will undergo different advection speeds (faster between 55 N and 34 N and slower south of there), which induces an anomalous heat convergence in the subpolar gyre (warming), but a heat divergence south of it (cooling). The time integral of this heat convergence/ divergence is argued to give rise to this dipole. Note, however, that first the warming is enhanced in the SPG and then followed by a cooling in the GS region some years after. Zhang and Zhang (2015) also point out that  the HC dipole between the SPG and GS regions is absent when there is no southward propagation of anomalies along the western boundary current.
Analysis of trends at the atmosphere-ocean interface suggests that heat fluxes cannot be considered the only source of these trends as there is a mismatch between the structure of the trends in ocean-atmosphere heat fluxes and those in SSTs (Fig. 6). Trends in mean sea level pressure (MSLP) have an effect on the Icelandic low and Azores high or, in other words, on the mean northeasterly and westerly winds to the north and south of Iceland, respectively, as well as the mean easterly trade winds in the tropics. During the 1996-2005 warming trend years enhanced surface-wind trends along the coasts of Iceland, Greenland, Newfoundland, United States and in the central mid-Atlantic enclose weakened wind trends along the path of the NAC. This pattern drives negative surface flux trends around the rim of the SPG and positive  surface flux trends in the centre of this gyre; the trends in surface fluxes should be producing a similar pattern of positive and negative SST trends if the latter were forced by the atmosphere. However, the structure of observed trends in SSTs shows a broad warming trend along the rim of the subpolar gyre and a confined cooling trend at its centre. During the 2006-2015 cooling period enhanced surface wind trends form a horse-shoe pattern appearing along the US and Newfoundland coasts that extends to Iceland and then southwestward along the Iberian Peninsula and western Africa until the tropics; the enhanced wind trends surround a region of weak wind trends in the centre of the basin from the subtropics to the midlatitudes around 45 N. The associated heat-flux trends follow this structure and suggest a cooling of the ocean over the Labrador basin, the SPG and the subtropical latitudes off the African coasts and a warming trend along the NAC and centre of the midlatitudes. However, trends in SSTs show a distinctive dipole pattern very much like the one seen in the HC trends suggesting that the atmosphere is not the main forcing of these ocean trends. Similar results are obtained if fields from NCEP/ NCAR reanalysis (Kalnay et al., 1996) are used. In fact, these results chime with those by Robson et al. (2016 -Fig . 2) regarding the oceanic nature of the observed warming and cooling trends before and after 2005.
While it is important to analyze the structure of the decadal trends, it is also illustrative to analyze the mean decadal anomalies and the evolution of both decadal anomalies and trends in order to elucidate an underlying mechanism. The structure of the warming and decadal trends before and after 2006 shows some features over the SPG and the GS-NAC system that resemble the evolution of upper-ocean HC associated with the AMO  (d) 1965-1974 and 1964-1973, (e) 1996-2005 and 1995-2004 and (f) 2006-2015 and 2005-2014 in units of 10 8 J m À2 . Red/blue shading corresponds to positive/negative anomalies and their change. Solid black lines mark the climatological annual mean position of the subpolar and subtropical gyres using the À0.4 m and þ0.4 m absolute dynamic topography values, respectively, and the dashed black line tracks the climatological annual-mean position of the North Atlantic current through the À0.1 m topographic contour from AVISO altimetry. discussed in Kavvada et al. (2013), where contrasting HC anomalies between the GS-NAC and SPG regions, associated with the AMO, change signs within a decade. Hence, it should be taken into consideration that decadal variations in the GS-NAC-SPG system may induce decadal trends. The basin-wide AMO and AMO-HC time series in Fig. 1 show decadal cooling trends only in two periods: the current one after 2006, already described, and an early one in the 1960s; the cooling trend in the 1960s is more pronounced at the surface (AMO) than in the subsurface layers (AMO-HC) and ceases later in the mid-70s. Figure 7 contrasts 10-yr mean anomalies and their tendencies from one decade to the next one after one year in the upper HC for the 1965-1974, 1996-2005 and the 2006-2015 decades in order to highlight the regions changing the most from one decade to the immediate one in the smoothed (decadal) data. Independent of the sign (or the decade), the mean anomalies of one sign occupy almost the entire basin with maximum values along the GS region and have anomalies of the opposite sign in the SPNA region with varied extents (Fig. 7, upper row). It is significant to note the following two things. Firstly, these anomalies in the SPNA region occupy different areas  Fig. 1b. (b) Profiles of lead/lag correlations between smoothed detrended surface SPNA index and smoothed detrended subsurface temperature anomalies at the 45 N-60 N and 30 N-45 N latitude bands; the vertical axis is depth in m, horizontal axis is lead/lag in years, and red/blue shading corresponds to positive/negative correlations; the grey stippling denotes significant grid points at the 95% confidence level based on a 2-tailed Student's t-test. Fig. 9. Area-averaged wind anomalies and profiles of subsurface water property anomalies in the subpolar region (60 W-10 W, 45 N-60 N) from 1950 to 2015. (a) Wind-speed anomalies at 850 mb from the NCEP/NCAR reanalysis, (b) temperature (shaded) and salinity (contours) anomalies, (c) density (shaded) and salinity forced density anomalies (contours) and (d) density (shaded) and temperature forced density anomalies (contours). Salinity/temperature-forced density anomalies are calculated considering climatological temperature/salinity in the UNESCO formula. Wind speed is in ms À1 , temperature in K, salinity with no units, and densities are in kgm À3 . Red/blue shading and continuous/dash lines correspond to positive/negative anomalies.
within the region but all show anomalies just to the east of the Flemish Cap around 45 W, 45 N, that is, a region where the GS current in its northeastward path meets the southward side of the SPG. Secondly, the similarities between the mean anomalies for the 1996-2005 warming period and those for the 2006-2015 cooling period are notable. The 2006-2015 anomalies of one sign over the Reykjanes Ridge and anomalies of the opposite sign in the Iceland basin appear to have developed from the anomalies in the 1996-2005 decade: the cold central-gyre anomaly and the anomaly to the east of the Flemish cap. It has been argued by Nigam et al. (2018) that such structures in HC may be the result of a southward intrusion of subpolar waters through the Newfoundland Basin, which promotes a detachment of GS anomalies in its northeastern section that eventually cross and invade the subpolar region via the NAC pathways. Along this line, it is seen that the tendency in the decadal pattern (Fig. 7, lower row) highlights the regions of largest change from the previous decade and include the GS-NAC-SPG regions, which are the same regions where the largest trends are found (Fig. 5, upper row for trends before and after 2006). Mean anomalies in vertically averaged salinity and their tendencies (Supplementary Fig. 2S) highlight that anomalously cold and warm anomalies are accompanied by fresh and saline anomalies, respectively, which reinforces the role for ocean circulation as the main driver of these large-scale heat anomalies.

Variability in the subpolar North Atlantic and the Gulf stream
The previous two sections have shown the prominence that the subpolar (via the SPG and associated currents) and subtropical (via the GS) regions have in shaping the observed decadal anomalies and, in consequence, the decadal trends. Because near-surface temperatures in the subpolar latitudes occur (almost) simultaneously with the AMO (Fig. 1b), attention is focused on this region. A series of indices are created for area-averaged SST and heatcontent anomalies within the SPG region (60 W-10 W, 45 N-65 N). These indices are highly coherent and show maximum amplitude in 2006-2007 before shifting to a cooling trend (Fig. 8a), which is steeper in the subpolar region than in the whole basin, suggesting the SPG plays a prominent role in leading the basin-wide trends. The SPG surface thermal variations lag those from the subsurface subtropics by $5 years (Fig. 8b, right panel), and propagate with time to the subsurface layers and reach the 1000 m column around 5 years later (Fig. 8b, left  panel). In particular, the upper ocean layer of the SPNA region was under cold and fresh conditions during the period 1970-1995, and has been under anomalous warm and saline conditions since 1995, reaching maximum anomalies and penetration around 2006 (Fig. 9). A reversal to cold and fresh conditions start to appear in the region in 2015, suggesting a return to a cold phase seen in the 90s, if not more pronounced this time (see e.g. Yashayaev and Loder, 2016). It is thus important to note that decadal cooling was already underway since mid-2000s (Zunino et al., 2017, their Fig. 9), which makes the so called 'cold blob' (Duchez et al., 2016;Josey et al., 2018) only part of decade-long cooling trend that intensified during winter 2014 due to extreme surface heat loss (Zunino et al., 2017;Josey et al., 2018).
It is apparent that temperature and salinity anomalies deepen with time, and that the cold and warm upperocean anomalies co-vary with positive and negative wind speed anomalies, respectively, however the warming in the 2000s is not commensurable to the weak reduction in wind speed. It is important to point out that changes in HC can result from changes in air-sea heat fluxes and the convergence/divergence of ocean heat transports (see e.g. Barrier et al., 2014), with the latter a result of not just AMOC-related changes, but also wind-driven. An important aspect to note in Fig. 9 is that the potential density anomalies in the subpolar region are largely responsive to the temperature anomalies, except around 1972 when the salinity effect (Great salinity anomaly effects, Dickson et al., 1988) on density dominated over the temperature effect.
Anomalous conditions over the SPNA region (as over the entire basin) seem to be preceded by anomalous conditions at the subsurface in the subtropics (Fig. 8b, right  panel). On the other hand, variability on this subtropical region is largely confined to the western basin along the domain of the GS current ( Supplementary Fig. 1S), so it is reasonable to expect that an index for the northward excursions of the GS current (Fig. 10a) leads the subsurface variability of the SPNA region as a whole (by $13 years); but even more revelatory is the fact that the lead increases from the northeast Atlantic at $9 years, and to the central and western SPNA at 13 and 14 years, respectively (Fig. 10b, c). Similar lag years with the GS index are obtained if instead of using the area-averaged HC over the three different basins (Iceland Basin, Irminger Sea and Labrador Sea) along the mean circulation, specific points within them are used. In fact, regressions of the smoothed GS index on SST, upper HC and averaged salinity show how in a span of about 12 years warm and saline anomalies from the GS reach the subpolar region (Fig. 11). The propagating signal along the GS-NAC-SPNA path with an 8-14-year time scale agree with the findings of Sutton and Allen (1997) of northward advection of SST anomalies in North Atlantic, with the modelling findings of Teng et al. (2011) of advection of subsurface temperatures, and recently supported by Årthun et al. (2017). The latter investigators demonstrated that poleward propagating SST anomalies from the northern extension of the Gulf Stream all the way to Nordic Seas are independent of the prevailing atmospheric circulation and predominantly controlled by ocean dynamics. Note, however, that the spatial progression of anomalies in Fig. 11 appear to be better captured by upper-ocean HC rather than by SSTs, which is likely due to the subsurface nature of the GS index (Joyce and Zhang, 2010). Figure 11 further shows that while the GS is advecting anomalies in a northeasterly direction, anomalies of opposite sign are being displaced westwards by the northern branch of the SPG to eventually turn south in front of the Newfoundland coast; the southward displacement of anomalies stops around the Flemish Cap (45 W, 45 N) where it cuts off the GS current and leaves clipped anomalies to continue their northeastwards displacement toward the Icelandic and Irminger basins (Nigam et al., 2018). Similarities between regressed SST anomalies from the GS index at t À 4 years (Fig. 11, panel t þ 4 yrs) and those from the AMO-HC (5-657 m), assuming its negative phase, at t þ 10 years (Fig. 2, panel t þ 10 yrs) are worth noting, suggesting a link between these phenomena. However, the mechanisms through which the GS and the AMO are connected is not a goal in the current paper, and the reader is referred to the recent study by Nigam et al. (2018).  -50 W, 33 N-43 N); the SPNA index used in (a) is area-averaged over the blue rectangle (60 W-10 W, 45 N-60 N); a western SPNA index is defined by area-averaging over the red rectangle (60 W-44 W, 45 N-60 N); a central SPNA index is defined by area-averaging over the green rectangle (44 W-27 W, 45 N-60 N); an eastern SPNA index is defined by area-averaging over the light blue rectangle (27 W-10 W, 45 N-60 N). (c) Lead/lag correlations between smoothed detrended GS index and SPNA HC-based indices. Note that the GS leads the subpolar indices and that the farther to the west the SPNA index is defined, the larger is the lag. Fig. 11. Lead-lag regressions on detrended anomalies of the normalized smoothed (LOESS25%) GS index for the period 1954-2012. Panels in the left column show regressed SST anomalies in units of 10 À1 K, those in the central column display upper-ocean HC (5-657 m) anomalies in units of 10 8 J m À2 , and panels in the right column show vertically averaged (5-657 m) salinity anomalies with no units at 2-year intervals from t À 6 years (top row) to simultaneous (central row) to t þ 6 years (bottom row). Time runs down as indicated by the long red arrow to the left. Red/blue shading represents positive/negative anomalies. Solid black lines mark the climatological annual mean position of the subpolar and subtropical gyres using the À0.4 m and þ0.4 m absolute dynamic topography values, respectively, and the dashed black line tracks the climatological annual-mean position of the North Atlantic current through the À0.1m topographic contour from AVISO altimetry. The grey stippling denotes significant grid points at the 95% confidence level based on a 2-tailed Student's t-test.

Discussion and conclusions
The subsurface spatiotemporal structure of the North Atlantic Ocean is analyzed in the context of the recent cooling trend. In this analysis a vertically integrated (5-657 m) ocean heat content AMO index has been proposed in an effort to capture the thermal state of the ocean and not just that of the sea surface as in the canonical AMO SST-based indices, where both are defined over the same region. Regressions of the traditional AMO SST index and our new AMO-HC (5-657 m) index on subsurface temperatures and salinity anomalies highlights decadal interactions not only between the tropical, subtropical and subpolar latitudes but also between subsurface and surface anomalies with the potential to provoke decadal trends in the North Atlantic Ocean. The use of this HC-based index has been revealing as it indicates that: (1) The North Atlantic subsurface had a warming trend since the mid-1980s to the mid-2000s, a feature that was also present at the surface with an evident lag of $3 years with respect to the AMO-HC (5-657 m) index. (2) The North Atlantic subsurface shows a cooling trend since the mid-2000s, a feature that is leading the ocean surface and one with significant implications for predicting future North Atlantic climate (Yeager and Robson, 2017). (3) The traditional surface AMO index (SST-based) lags the subsurface indices (HC-based) with the leading time being latitude dependent: deeper waters centred around 350 m in the subtropics (30 N-45 N) lead by $5 years, while the surface waters centred around 50-100 m depths in the midlatitudes (45 N-60 N) lead by $0-1 year.
The spatial structure of the trends in the North Atlantic after 2006 suggests a link with decadal AMOC variability (Jackson et al., 2016;Robson et al., 2016), but other processes related to gyre-circulation changes could also have played a role (Bersch et al., 2007;H€ akkinen et al., 2011;Piecuch et al., 2017). In particular, Piecuch et al. (2017) suggest a role for mid-latitude wind stress curl changes at the subpolar-subtropical gyre boundary for the development of the recent decadal cooling of the SPNA. However, Robson et al. (2016) argue for changes of the Labrador Sea density and hence the overturning. Jackson et al. (2016) showed using an ocean reanalysis that the AMOC started a decadal decline around mid-2000s and therefore part of an internal, not forced, variability since it is a recovery from the earlier strengthening. It, thus, appears that both the gyre and overturning circulation likely had an effect on the observed decadal climate reversal to a cooling trend (Robson et al., 2016;Piecuch et al., 2017).
The structure of the 20-year warming trend before 2006 is similar to the 10-year warming trend before 2006, albeit weaker (cf. Figs. 4 and 5). In both cases warming trends surround cold trends along the GS-NAC meeting point with the cold trends reaching higher northern latitudes in the recent 10-year period than in the 20-year period before 2006. The incipient cooling trend along the GS-NAC system before 2006 evolves to the cooling trend covering the NAC-SPG region in the following 10 years after 2006. Analysis of surface heat fluxes only partially explains the observed SST warming and cooling trends before and after 2006, suggesting that ocean dynamical processes are the primary source and hence critical to the observed trends (Robson et al., 2016;Zhang et al., 2016;Piecuch et al., 2017;Zhang, 2017).
Mean decadal anomalies, their change and trends all suggest a link with variability generated by the GS-NAC-SPG system. It is not simple to locate where conditions start but at some point in the evolution of the mean HC and salinity decadal anomalies these are of one sign along the GS region, and of opposite sign in the eastern portion of the SPG region. While the anomalies along the GS move northeastwards, anomalies from the SPG move southwards along the North American coasts and the western boundary. The southward displacement of anomalies eventually reaches the region where the anomalies of opposite sign are driven by the GS current and 'clip them off' to follow the NAC. This southward extent of anomalies of one sign and the cutoff anomalies of the opposite sign sets up conditions to reverse the sign of the anomalies in both the GS and SPNA regions; an interaction that adds to the list of complex processes that are at the core of the North Atlantic decadal variability (Zhang, 2008;Joyce and Zhang, 2010;Zhang and Zhang, 2015;Delworth et al., 2017;Nigam et al., 2018).
Trends, variability as well as the near-surface subpolar and subsurface subtropical connection captured by the basin-wide (AMO and AMO-HC) indices are also identified by SPNA (SST-and HC-based) indices. The warming and cooling trends distinguished from the whole-basin indices (AMO and AMO-HC) are also found in the SPNA indices (SST-and HC-based) but with steeper values. Surface anomalies in the SPNA region as a whole penetrates into the deeper layers with time, however, water properties (i.e. density) are not homogeneous as the eastern and western SPG is more sensitive to changes in temperature and salinity, respectively.
On the other hand, and as in the case of the whole basin characterized with the AMO index, surface variability in the SPNA region is led by variability in the subsurface subtropics (30 N-45 N), and the main source of variability in the Atlantic Ocean in this latitudinal band comes from the Gulf Stream itself. The GS in fact leads variability in the upper HC (5-657 m) over the SPNA region as a whole by $13 years, but the lead increases from east to west from $9 years (over the Iceland Basin) to $13 years (over the Irminger Sea) to $14 years (over the Labrador Sea). This connection is reflected in regressions of a GS index on SST, upper-ocean HC and salt content that present coherent anomalies of one sign along the GS region and anomalies of the opposite sign over the SPNA region that evolve and change sign along the GS-NAC-SPG system on decadal time scales.
A better understanding of the decadal variations in the Atlantic is of paramount importance for the societies inhabiting the continents around it. Warming of the North Atlantic waters has been linked to multiyear droughts over the central United States (e.g. , anomalous rainfall over western Europe and the Sahel (e.g. Sutton and Hodson, 2005;Nigam and Ruiz-Barradas, 2016) and increased hurricane activity (e.g. Goldenberg et al., 2001; while the opposite regional climate conditions, including heat waves over Europe (Duchez et al., 2016;Josey et al. 2018), are observed during colder episodes of the North Atlantic waters.