Interannual to decadal sea level variability in the subpolar North Atlantic: The role of propagating signals

. The gyre-scale, dynamic sea surface height (SSH) variability signifies the spatial redistribution of heat and freshwater in the ocean, influencing the ocean circulation, weather, climate, sea level, and ecosystems. It is known that the 10 first empirical orthogonal function (EOF) mode of the interannual SSH variability in the North Atlantic exhibits a tripole gyre pattern, with the subtropical gyre varying out-of-phase with both the subpolar gyre and the tropics, influenced by the low-frequency North Atlantic Oscillation. Here, we show that the first EOF mode explains the majority (60-90%) of the interannual SSH variance in the Labrador and Irminger Seas, whereas the second EOF mode is more influential in the northeastern part of the subpolar North Atlantic (SPNA), explaining up to 60-80% of the regional interannual SSH variability. We find that the 15 two leading modes do not represent physically independent phenomena. On the contrary, they evolve as a quadrature pair associated with a propagation of SSH anomalies from the eastern to the western SPNA. This is confirmed by the complex EOF analysis, which can detect propagating (as opposed to stationary) signals. The analysis shows that it takes about 2 years for sea level signals to propagate from the Iceland Basin to the Labrador Sea, and it takes 7-10 years for the entire cycle of the North Atlantic SSH tripole to complete. The observed westward propagation of SSH anomalies is linked to shifting wind 20 forcing patterns and to the cyclonic pattern of the mean ocean circulation in the SPNA. The analysis of regional surface buoyancy fluxes in combination with the upper-ocean temperature and salinity changes suggests a time-dependent dominance of either air-sea heat fluxes or advection in driving the observed SSH tendencies, while the contribution of surface freshwater fluxes (precipitation and evaporation) is negligible. We demonstrate that the most recent cooling and freshening observed in the SPNA since about 2010 was mostly driven by advection associated with the North Atlantic Current. The results of this 25 study indicate that signal propagation is an important component of the North Atlantic SSH tripole, as it applies to the SPNA.


Introduction
Ocean and atmosphere dynamics induce regional sea level changes with amplitudes that are often several times greater than the global mean sea level rise of about 3.5 mm yr -1 (Meyssignac & Cazenave, 2012;Cazenave et al., 2014;Volkov et al., 2017;Chafik et al., 2019). On the timescales from seasonal to longer, this dynamic sea level variability (after the global mean sea 30 level has been subtracted) is mainly due to changes in the density of water column, driven by surface buoyancy fluxes and the advection of heat and freshwater by ocean currents (Gill & Niiler, 1973;Ferry et al., 2000;Volkov & van Aken, 2003;Cabanes et al., 2006). The associated changes in the temperature and salinity of the water column determine the thermosteric and halosteric components of sea level variability, respectively, which often offset each other. Away from polar regions, the dynamic sea level variability generally follows the thermosteric changes and, therefore, represents a good proxy for the upper-35 ocean heat content variability (e.g., Chambers et al., 1998).
It has been shown that on interannual to decadal timescales, the variability of sea surface temperatures (SST), sea surface height (SSH), and oceanic heat content in the North Atlantic exhibits a tripole pattern with an out-of-phase relationship between the subtropical gyre and both the subpolar gyre and the tropics (Tourre et al., 1999;Watanabe & Kimoto, 2000;Marshall et al., 2001;Häkkinen, 2001;Volkov et al., 2019aVolkov et al., , 2019b. The tripole is obtained through the widely used Empirical Orthogonal 40 Functions (EOF) decomposition of the respective variables. Specifically, the North Atlantic SSH tripole is the leading EOF mode of the interannual dynamic SSH variability (Volkov et al., 2019b). Earlier studies suggested that both the SST tripole and the SSH tripole represent the ocean's response to the atmosphere-ocean heat flux (Cayan, 1992;Häkkinen, 2001).
However, more recent research demonstrated that the upper ocean heat content and, consequently, SSH change in the North Atlantic could be dominated by the oceanic heat transport divergence (Häkkinen et al., 2015;Piecuch et al., 2017). The latter 45 is modulated by the Atlantic meridional overturning circulation (AMOC). For example, an abrupt reduction in the AMOC observed at 26.5ºN in 2008-2009 led to a cold SST anomaly and a decrease of the thermosteric sea level (and, hence, the upper ocean heat content) in the entire subtropical gyre in 2009(Josey et al., 2018Volkov et al., 2019a). Using observations and an ocean model, Volkov et al. (2019b) showed that the tripole-related SSH variability between 26ºN and 45ºN (the subtropical band of the tripole) is correlated with the AMOC and meridional heat transport at 26.5ºN, and with the low-50 frequency North Atlantic Oscillation (NAO).
In this study, we focus on the subpolar North Atlantic (SPNA; Fig. 1), which is one of the key regions for deep water formation and, as such, it plays an important role in driving the AMOC. Warming, salinification, spindown, and contraction of the North Atlantic subpolar gyre were observed in the 1990s to mid-2000s (Häkkinen and Rhines, 2004;Sarafanov et al., 2008;Holliday et al., 2008). This warming was associated with the positive phase of the Atlantic Multidecadal Oscillation (e.g., Enfield et al., 55 2001), partly driven by the strengthening of the North Atlantic Current (NAC) and the AMOC (Msadek et al., 2014). A cooling tendency in the SPNA started in 2006, and it ended with an exceptionally cold anomaly in 2015, termed the 'cold blob' (Ruiz- Barradas et al., 2018;Chafik et al., 2019). Several studies have attributed the 'cold blob' observed in the SPNA to the dramatic wintertime heat loss event in 2014-15 (Josey et al., 2018;Grist et al., 2016;Duchez et al., 2016). Bryden et al. (2020), however, suggested that the 'cold blob' could also be linked to the 2008-2009 reduction of the AMOC observed at 26.5ºN. Piecuch et 60 al. (2017) showed, however, that horizontal gyre circulations provide a greater contribution to heat divergences in the SPNA than vertical overturning circulations.
Along with the cooling anomaly, a slowdown of the warm and salty NAC and redistribution of the fresher Arctic waters led to an unprecedented freshening of the SPNA in 2012-2016 (Holliday et al., 2020) that rapidly communicated with the deep waters (Chafik and Holliday, 2022). The 2013-2015 cold event in the SPNA also intensified deep ocean convection in the 65 Irminger and Labrador basins (de Jong & de Steur, 2016;Yashayaev & Loder, 2016). While the reduced AMOC may have influenced the intense cooling in the subpolar gyre around 2015, the resulting intensification of deep convection is expected to eventually increase the strength of the AMOC (Frajka-Williams et al., 2017). A recent analysis showed that the [2006][2007][2008][2009][2010][2011][2012][2013][2014][2015] cooling trend in the SPNA may have reversed in 2016 with a large-scale warming in the central and eastern parts of the domain due to an enhanced advection of warm and saline waters from the subtropical gyre (Desbruyères et al., 2021). 70 Since large-scale ocean circulation is an important driver for the observed low-frequency changes in the North Atlantic, propagation of signals across the basin and between the subtropical and subpolar gyres is an intrinsic element of the variability.
While westward propagating Rossby waves and eddies provide an effective mechanism for energy transfer in the tropics and within the subtropical gyres, advection by ocean currents plays a primary role for the meridional and zonal transports along western boundary currents and at high latitudes. For example, Fu (2004) identified a propagation of the interannual sea level 75 signal from the subtropical region to the eastern end of the Gulf Stream extension, demonstrating that the SSH variability is not just a steric response to the heat flux forcing, but also involves a dynamic response. Using a simple thermodynamic model, Dong and Kelly (2004) showed that advection by geostrophic currents can largely contribute to interannual variations of the upper ocean heat content in the Gulf Stream region. Desbruyères et al. (2021) hypothesized that the warming signal observed in the eastern SPNA since 2016 would progressively propagate westward into the Irminger and Labrador Basins. 80 In areas where signal propagation is important, analyses based on the EOF decomposition may be misleading, because each individual EOF pattern represents a standing oscillation. Propagating signals are often represented as a linear combination of more than one EOF pattern, such as a pair, with one mode leading the other by 90° of phase in space and time (e.g., Roundy, 2015). In this case, treating consecutive EOF modes as independent physical processes is not appropriate. In this paper, we show that signal propagation is an intrinsic element of the low-frequency dynamic SSH variability in the SPNA. Therefore, 85 the main objectives of this study are (i) to revise the definition of the North Atlantic SSH tripole by accounting for signal propagation, and (ii) to explore the fidelity of the tripole in describing sea level variability in the SPNA. Nearly 30 years of quasi-global, regular satellite altimetry measurements allow us not only to resolve the interannual variability, but also gain a preliminary insight into decadal changes. Therefore, we specifically focus on interannual to decadal time scales. In addition, we discuss how the leading modes of the variability in the SPNA are related to atmospheric wind and buoyancy forcing and 90 to advection by ocean currents.

Satellite altimetry measurements
Satellite altimetry has provided accurate, nearly global, and sustained observations of sea level since the launch of the Topex/Poseidon mission in August 1992. We use the monthly and daily maps of SSH anomalies for the time period from 95 January 1993 to December 2020 processed and distributed by the Copernicus Marine and Environment Monitoring Service (CMEMS; http://marine.copernicus.eu). The daily maps are used only for the computation of eddy propagation velocities (Section 3.3). The maps are produced by the optimal interpolation of measurements from all the altimeter missions available at a given time. Prior to mapping, the along-track altimetry records are routinely corrected for instrumental noise, orbit determination error, atmospheric refraction, sea state bias, static and dynamic atmospheric pressure effects, and tides (Pujol et 100 al. 2016). The SSH anomalies produced by CMEMS are computed with respect to a twenty-year (1993-2012) mean sea surface, but we center them around the record-long mean. We also subtract the global mean sea level from SSH anomaly time series at each grid point to focus on local dynamic sea level variability not related to global changes.

Hydrographic data
To relate sea level variability to subsurface processes, we use the EN4 monthly gridded profiles of temperature and salinity 105 (version EN.4.2.2 with mechanical/expendable bathythermograph (MBT/XBT) profile data bias adjustments from Gouretski and Reseghetti, 2010) for the period from January 1993 to December 2020 (Good et al., 2013). The profiles are based on the objective analysis of hydrographic observations (e.g., XBT, MBT, bottle, CTD, and Argo). The number of Argo profiling floats in the ocean was growing from a very sparse array of 1000 profiling floats in 2004 to a global array of more than 3000 instruments from late 2007 to the present. This means that maps in 2004-2007 and especially before 2004 are less accurate 110 than maps in 2008-2020. The EN4 temperature and salinity profiles are used to calculate the steric SSH anomalies (SSHST; due to changes in density) by integrating in-situ density anomalies with respect to the time mean and over the upper 1,000 m, which is the depth interval occupied by the NAC. The thermosteric (SSHT; due to temperature changes only) and halosteric (SSHS; due to salinity changes only) contributions to SSH are calculated by integrating in-situ density anomalies with respect to the time mean values of salinity and temperature, respectively. The obtained SSHT and SSHS are proportional to the upper 115 1,000 m heat and freshwater contents, respectively. Similar to the processing of altimetry data, the global mean SSHST, SSHT, and SSHS are subtracted from the respective time series at each grid point. To analyze the role of the mean ocean circulation, the individual Argo temperature and salinity profiles and float trajectories from the U.S. Argo Data Assembly Center (https://www.aoml.noaa.gov/argo/#argodata) are used to compute the time-mean adjusted geostrophic velocities at 1000-dbar (Schmid, 2014; see Section 3.3). 120 5

Atmospheric data
The observed changes in SSH are analyzed jointly with atmospheric forcing fields provided by ERA5 climate reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) and distributed by Copernicus (https://climate.copernicus.eu/climate-reanalysis). Specifically, we use the monthly averaged fields of sea level pressure (SLP), 10-m wind velocities, surface wind stress, and surface heat (shortwave and thermal radiation, the sensible and latent heat 125 fluxes) and freshwater (precipitation, evaporation) fluxes from January 1979 to December 2020 (Hersbach et al., 2019). In addition, we also use the monthly station-based NAO index, based on the normalized SLP difference between Lisbon (Portugal) and Reykjavik (Iceland) from the Climate Analysis Section of the National Center for Atmospheric Research (Hurrell et al., 2003;https://climatedataguide.ucar.edu).

Methods 130
To focus on interannual and longer time scales, the seasonal cycle is computed by fitting both the annual and semi-annual harmonics in a least squares sense and subtracted from all data fields and time series. The time series are further low-pass filtered with a "Lowess" filter with a cutoff period of 18 months.

EOF analysis
Empirical orthogonal function (EOF) analysis is widely used to reduce the dimensionality of datasets and to extract the leading 135 (most influential) modes of variability. Individual EOF modes, if found meaningful, are often assumed to describe different physical processes. Here, we use the conventional EOF analysis (Navarra & Simoncini, 2010) to identify the leading stationary modes of the low-frequency SSH variability in the North Atlantic. The EOF analysis is applied to both the satellite altimetry (SSH) and the EN4-derived (SSHST, SSHT, and SSHS) data. Because of the scarcity of in situ observations prior to the advent of Argo, only the EN4 data starting from 2004 was used for the EOF analysis. 140 For each mode j, the spatial pattern (map) is represented as a regression map (EOFj) obtained by projecting SSH data onto the standardized (divided by standard deviation) principal component (PCj) time series. Thus, the regression coefficients are in centimeters (local change of sea level) per 1 standard deviation change of the associated PCj. The SSH fields can be reconstructed using a limited number of modes, N: (2).

6
The conventional EOF analysis is not an effective method to identify propagating modes of variability. Therefore, the detection of propagating modes of the low-frequency SSH variability in the North Atlantic is carried out using the Complex EOF (CEOF) 150 analysis, which yields the spatial and temporal amplitude and phase information (Navarra & Simoncini, 2010). In essence, the CEOF analysis is based on the notion that a propagating signal should contain signals that are orthogonal to each other. In practice, the CEOF analysis of a data field is similar to the conventional EOF analysis, but applied to the field augmented in a manner such that propagating signals within it may be detected (Navarra & Simoncini, 2010). The augmentation of the data field is achieved through forming complex time series: where the real part is simply the original data field (SSH) and the imaginary part is its Hilbert transform (SSH H ), representing a filtering operation upon SSH(x, t) in which the amplitude of each spectral component is unchanged but each component's phase is advanced by /2 (Horel, 1984). Unlike the conventional EOF analysis, the CEOF analysis yields complex eigenvectors for the spatial patterns (maps) and for their principal components (time evolution). For simplicity, henceforth, we refer to the 160 spatial patterns as CEOFs (CEOFj) and to the complex principal components as CPCs (CPCj). Similar to the conventional EOF analysis, the real and imaginary parts of CEOFs are presented by regressing SSH onto the standardized real and imaginary parts of CPCs. The time-space progression of the CEOFj mode is obtained by multiplying CEOFj(x) and CPCj(t) by a rotation matrix whose argument may vary from 0-360°. The real and imaginary parts of CEOFs and CPCs are also used to obtain the spatial (Φ) and temporal (θ) phases of the CEOFs:

Forcing mechanisms
The statistical relationship between SSH and atmospheric circulation is established by regressing SLP and 10-m wind velocity fields on PCj and CPC1. It is reasonable to assume that on interannual-to-decadal time scales the SSH changes are mostly steric 170 (e.g., Volkov and van Aken, 2003) and, therefore, the mass-related changes can be neglected. The variations of SSHST are related to heat and freshwater fluxes as follows (e.g., Cabanes et al., 2006): where #5= = +! + =! + +, + C, is the net surface heat flux (positive into the ocean) equal to the sum of shortwave radiation (QSR), thermal (longwave) radiation (QTR), the sensible heat flux (QSH), and the latent heat flux (QLH); AB = − 175 is the surface freshwater flux (positive into the ocean) equal to the sum of precipitation (P) and evaporation (E); ɑ and β are the coefficients of thermal expansion and haline contraction averaged over the upper 1000 m, respectively; ⍴0 is the reference density; cp is the specific heat of seawater; Sa is the salinity anomaly relative to a multi-year (2004-2020) mean value averaged 7 over the upper 1000 m; ADV denotes the SSHST change due to the advection of density anomalies; and the overbar indicates the climatological averages over the entire ERA5 record . 180

Ocean circulation and eddy propagation
The potential impact of oceanic advection in the SPNA is qualitatively analyzed using (i) the mean 1000-dbar velocities obtained from Argo data and (ii) eddy propagation velocities estimated from daily satellite altimetry maps, computed for the time period of 2005-2019.
The time-mean gridded adjusted geostrophic velocities at 1000-dbar are derived from Argo and altimetry using the 185 climatological velocity field from Argo trajectories as the reference velocity following the methodology of Schmid (2014).
This method uses Argo dynamic height profiles and SSH from altimetry to derive synthetic dynamic height profiles on a 0.5° grid. These profiles are then used to derive the horizontal geostrophic velocities, followed by the barotropic adjustment, based on the climatology of the velocity derived from float trajectories at about 1000-dbar (the target parking depth of Argo floats).
Therefore, the obtained velocities at 1000-dbar are mainly constrained by float trajectories and have very little temporal 190 variability. Because the vertical extent of ocean eddies can reach 1000-dbar, eddies that are involved in horizontal property transports tend to follow the mean flow field at about the same depth.
To supplement the 1000-dbar velocities, we also compute eddy propagation velocities for the same time interval using a spacetime lagged correlation analysis of daily SSH fields as described in Volkov et al. (2013) and Fu (2006). A velocity vector at a particular grid point is obtained as follows. First, the correlations between SSH anomalies at this grid point and at neighboring 195 grid points are computed at various time lags. Second, the location of the maximum correlation is recorded, and a velocity vector is computed using the distance between the two grid points and the time lag. Finally, an average velocity vector weighted by correlation coefficients is computed from velocity vectors at various time lags. To obtain the results characteristic for eddy temporal and spatial scales, the time lags were limited to less than 70 days and the horizontal dimensions of the area for computing the correlations between the neighboring grid points were set to about 200 km. Because SSH is an integral quantity 200 characteristic for the full-depth thermohaline properties, the obtained velocities are representative of property transports.

The imprint of SSH tripole in the subpolar North Atlantic
The North Atlantic SSH tripole is defined as the leading EOF mode (EOF1) of the low-pass filtered dynamic SSH (Volkov et al., 2019b). It is characterized by the subtropical band (~20°-45°N) varying out-of-phase with both the tropical North Atlantic 205 (south of ~20°N) and the subpolar gyre (~45°-65°N) (Fig. 2a). In 1993-2020, the tripole mode explained 27.2% of the interannual to decadal dynamic SSH variability in the North Atlantic. The time evolution of the tripole is shown by PC1 (solid blue curve in Figure 3). In 1993-2010, the tripole was characterized by a general reduction of sea level in the subtropical band 8 and a sea level rise in both the subpolar gyre and in the tropics. These tendencies reversed abruptly in 2011-2015. During the latter time interval, the subtropical gyre warmed considerably, which led to an accelerated sea level rise along the U. S. 210 southeast coast with the rates of up to five times the global mean (Domingues et al., 2018;Volkov et al., 2019b), while a strong cooling occurred in the subpolar North Atlantic in 2013-2015 (Chafik et al., 2019).
According to EOF1 (Fig. 2a) and PC1 (blue curve in Fig. 3), the subtropical and subpolar gyres have been in warm and cold states, respectively, since 2015. However, Desbruyères et al. (2021) reported an upper ocean warming in the eastern SPNA that started in 2015. While this warming is apparently not part of the leading EOF mode, it appears to be consistent with the 215 time evolution (PC2; red curve in Fig. 3) of the second EOF mode (EOF2; Fig. 2b). The EOF2 mode explains 13.4% of the variance and exhibits a pattern that suggests changes in the strength and likely meridional shifts of the Gulf Stream to the east of Cape Hatteras. This mode is also responsible for decadal changes in the eastern and northeastern SPNA: a generally positive SSH tendency in 1993-2003, a SSH decrease in 2004-2014, and a SSH increase since 2015 (red curve in Fig. 3), consistent with the observations of Desbruyères et al. (2021). 220 Sea level changes associated with EOF1 and EOF2 are mostly steric, i.e., determined by density variations. This is confirmed by the spatial structures (Figs. 2 c,d) and temporal evolutions (PCs; dotted curves in Fig. 3) of EOF1 and EOF2 of the low-pass filtered SSHST accounting for 45.2% and 16.8% of the variance, respectively. These two leading modes are similar to those of the low-pass filtered SSH (Figs. 2 a,b; solid curves in Fig. 3). The correlation between PC1 of SSH and PC1 of SSHST is 0.94, and the correlation between PC2 of SSH and PC2 of SSHST is 0.88, significant at 95% confidence. Furthermore, the steric sea 225 level changes are governed by changes in the upper 1000-m ocean heat content. The spatial patterns and the signs of EOF1 and EOF2 of the low-pass filtered SSHST are determined by the thermosteric sea level variability; the EOF1 and EOF2 of the lowpass filtered SSHT are nearly identical to the EOF1 and EOF2 of the low-pass filtered SSHST (compare Figs. 4 a,b and 2 c,d).
The correlation between the PC1 of SSHST and SSHT is 0.99, and the correlation between the PC2 of SSHST and SSHT is 0.97 The percentage of the local SSH variance explained by the two leading EOF modes in the SPNA exhibits the following patterns.
The North Atlantic SSH tripole, characterized by the EOF1/PC1 of the low-pass filtered SSH, explains the majority (60-90%) 235 of the interannual SSH variance in the Labrador Sea and in the Irminger Basins (western SPNA; Fig. 5a). The EOF2 mode explains a substantial amount of the interannual SSH variance (60-80%) in the northeastern SPNA, east of Greenland (Fig.   5b). The two leading modes together (EOF1+EOF2; Fig. 5c) explain most of the interannual SSH variability in the SPNA north of 52°N, except the southern part of the Iceland Basin and the Rockall Trough, where the eddy variability associated with the NAC branches is relatively strong. It is interesting to note that while the EOF1 depicts the interannual SSH variability over the 240 deep parts of the Labrador Sea and Irminger Basin, the EOF2 depicts the interannual SSH variability over the shallower waters of the Rockall Plateau, parts of European continental shelf and slope, Reykjanes Ridge, Iceland shelf, and along the Irminger and East Greenland Currents. Chafik et al. (2019) discussed the variability of dynamic SSH in 1993-2016 averaged over 5°-45°W and 55°-65°N in the SPNA 245 (green rectangle in Fig. 1). Here, we present both the updated SSH time series (black curve in Fig. 6) and the SSHR time series reconstructed with EOF1 (blue curve in Fig. 6) and EOF2 (red curve in Fig. 6), averaged over the same region. When disregarding shorter term, year-to-year variations, SSH increased by about 6 cm in 1993-2005, decreased by the same amount in 2006-2015, and then increased again by about 5 cm in 2016-2020. As the SSHR time series show, neither the EOF1 nor the EOF2 alone can adequately explain the SSH variability in the SPNA. However, their sum (dotted curve in Fig. 6) is sufficient 250 to reasonably reconstruct SSH in the region. The SSHR computed with the two leading EOFs matches the observed SSH well, with a correlation of 0.96, meaning that SSHR explains about 92% of the SSH variance. Individually, the EOF1 and EOF2 modes explain 44% and 48% of the SSH variance, respectively.

Interannual SSH variability in the subpolar North Atlantic
In order to illustrate the relative contribution of temperature and salinity changes to the interannual-to-decadal variability of SSH in the SPNA, we compute the averages of SSHST, SSHT, and SSHS over two regions (outlined by orange rectangles in 255 The time series of SSHST closely matches those of SSH, with the correlation between them above 0.95 and the root-meansquared differences of 0.4 in the eastern SPNA and 0.8 cm in the western SPNA (compare black and blue curves in Fig. 7).
This means that the variability of SSH in the SPNA is mostly steric in nature. The remaining difference between SSH and 260 SSHST can be attributed to density changes at depths greater than 1000-m and to errors in data; the contribution of barotropic signals is expected to be small at the time scales considered. In 2005-2015, SSHST decreased by about 6 cm in the eastern SPNA and 5 cm in the western SPNA (blue curves in Fig. 7). The decrease of SSHST was caused by cooling and the associated reduction of SSHT by about 12 cm in the eastern SPNA and 6 cm in the western SPNA (red curves in Fig. 7). The reduction of SSHT in both the eastern SPNA and western SPNA was compensated by freshening induced increase of SSHS by about 5 265 cm and 1 cm, respectively (green curves in Fig. 7). It is interesting to note that while the sign of SSHST anomalies is determined by SSHT, the contribution of SSHS is substantial and comparable to the contribution of SSHT, especially in the eastern SPNA.

Propagation of SSH anomalies
We have shown that SSH in the SPNA is adequately reconstructed by the two leading EOF modes. If the observed phase difference between these statistical modes is due to signal propagation, then both EOF1 and EOF2 describe the same physical 270 process at different stages of its evolution. Indeed, the application of the CEOF analysis presented below shows that EOF1 and EOF2 resemble the real and imaginary parts of CEOF1, respectively. Recall that the real and imaginary parts of CEOFj can be interpreted as the propagating signal at two stages separated by a quarter cycle. Displayed in Fig. 8 is the time evolution of the spatial pattern of SSH reconstructed with CEOF1 mode for one full cycle at phase stages separated by 45°. Note that SSH patterns at phases ±180° and 90° are almost identical to EOF1 and EOF2, respectively (Figs. 2 a,b). Likewise, illustrating the 275 temporal evolution of CEOF1 mode, the real part of CPC1 when rotated by 180° (Fig. 9, dashed blue curve) corresponds to PC1 (Fig. 3, blue curve) and the imaginary part of CPC1 (Fig. 9, red curve) corresponds to PC2 (Fig. 3, red curve).
The reconstruction of the CEOF1 mode (Fig. 8) shows that the tripole-related sea-level variations in the tropical and subtropical bands resemble a standing wave. Negative and positive SSH in the tropical and subtropical bands at phase 0° are gradually replaced by positive and negative SSH at phase 180°, respectively, without a notable sign of signal propagation. On the other 280 hand, a signal propagation is evident in the subpolar band of the tripole. In the first half of the cycle at phases -180° and -135°, As demonstrated by the real and imaginary CPC1 time series (blue and red curves in Fig. 9, respectively) and the temporal phase (dotted black curve in Fig. 9), there were almost three full cycles of the decadal gyre-scale sea-level changes that 290 exhibited signal propagation in the SPNA in 1993-2020. One full cycle is associated with the temporal phase change from −180° to 180°. In the northeastern SPNA, the first cycle started with a low SSH in 1993 (red curve) and it progressively reached high SSH in 1996, corresponding to phases −90° and 90° in Fig. 8

Relationship between sea level and wind forcing
The maximum correlation between the PC1 (solid blue curve in Fig. 3), showing the time evolution of the North Atlantic SSH tripole, and the low-pass filtered NAO index (shaded area in Fig. 3) is -0.73, with the NAO leading by 10 months (the 95% 305 significance level for correlation is about 0.45). The lag probably indicates the oceanic adjustment time to a variable wind forcing. Regression of SLP and 10-m winds on the PC1 displays a familiar NAO dipole pattern with a cyclonic (negative) anomaly in the subtropical high and an anticyclonic (positive) anomaly in the subpolar low SLP centers (Fig. 10a). The weaker/stronger subtropical high and subpolar low associated with weaker/stronger westerly winds in the midlatitude North Atlantic lead to lower/higher sea levels in the subtropical North Atlantic and higher/lower sea levels in the SPNA. 310 Interestingly, while correlation between the PC2 and the NAO is not significant, regression of SLP and 10-m winds on the PC2 (Fig. 10b) also exhibits a dipole pattern similar to the NAO, but with somewhat shifted pressure centers. This suggests that both the EOF1 and the EOF2 of the low-pass filtered SSH are possibly driven by the same atmospheric process, but are associated with different phases of its evolution or different regimes. To support this argument, we also regressed SLP and 10m winds on the real parts of CPC1 rotated every 45° between ±180° thus yielding the full-cycle evolution of wind forcing 315 patterns associated with the leading CEOF mode (Fig. 11). The obtained patterns illustrate that the distribution of SLP and 10m wind anomalies, as the subtropical high and the subpolar low centers change their position and intensity, is clearly related to the SSH patterns associated with the CEOF1 (compare Figs. 8 and 11). The positive (negative) SLP and anticyclonic (cyclonic) wind anomalies drive the near-surface Ekman convergence (divergence) and lead to the upper-ocean warming (cooling) and, consequently, higher (lower) sea levels. It should be noted that the SLP and wind anomaly patterns associated 320 with PC1 and PC2 (Fig. 10) are identical to wind forcing patterns associated with CPC1 at phases ±180° and 90° (or at phases 90° and 0° if the sign is reversed) (Fig. 11).
The westward propagation of negative/positive SSH from the northwest European shelf towards the Labrador Sea clearly follows similar shifts in the position of SLP and wind anomalies. Indeed, the center of the positive (anticyclonic) SLP anomaly associated with EOF2 or with CEOF1 at phase 90° is located near eastern Iceland, and the anticyclonic atmospheric circulation 325 pattern covers the area of positive SSH anomalies in the eastern-northeastern part of the SPNA (compare Figs. 10b, 11, and 2b). The positive SLP anomaly pattern then intensifies and shifts westward towards the eastern coast of Greenland (Figs. 10a and 11 at phases 135° and 180°). This shift corresponds to the observed propagation of SSH anomalies towards the Irminger Basin and the Labrador Sea (Fig. 11 at phases 90°, 135°, and 180°). The established statistical relationship suggests that SSH and the upper-ocean heat content anomalies in the SPNA result from the oceanic adjustment to persistent, lasting longer than 330 a year, wind forcing. On the other hand, we acknowledge that oceanic feedback to atmospheric forcing is also possible on these time scales.
The analysis presented here indicates that the NAO, which is correlated with the PC1 only, is not a sufficient proxy of wind forcing over the SPNA. This is mostly because the definition of the NAO is limited by an assumption that it is a standing oscillation pattern with the fixed subtropical high and the subpolar low-pressure centers. We have shown that the observed 335 propagation of the dynamic SSH anomalies in the SPNA is linked to the NAO-like atmospheric pressure and wind patterns that change their location and intensity. The importance of accounting for the position and intensity of atmospheric center of action for determining causal relationships within the coupled ocean-atmosphere system in the North Atlantic has been reported earlier (e.g., Hameed and Piontkovski, 2004;Hameed et al., 2021). Our results provide further evidence that the attribution of oceanic signals to the conventional NAO index, which generalizes and often over-simplifies the atmospheric variability, is not 340 always appropriate.

Relationship between sea level and surface buoyancy forcing
To investigate the impact of surface buoyancy forcing on the interannual-to-decadal changes of SSH, we select four time intervals characteristic for the main tendencies reflected by PC1 and PC2 in Fig. 3 (by the real and imaginary CPC1 in Fig. 9 Because the SSH changes driven by surface freshwater fluxes (P+E) appear to be three orders of magnitude smaller than those driven by surface heat fluxes, they are not considered in the following. The impact of surface freshwater fluxes on the regional 350 SSHS changes (Figs. 12-15, d) is also negligible, thus, suggesting that these changes are mainly driven by the advection of freshwater. The overall increase of sea level in the Labrador Sea and Irminger Basin in 1994-2010 associated with the North Atlantic SSH tripole and depicted by the PC1 (blue curve in Fig. 3) and the real part of CPC1 (dotted blue curve in Fig. 9) is well explained by the SSHST change (Figs. 12 a,b). The SSHST change exceeds the total SSH change only over the northern part of the Iceland 355 Basin and over the Rockall Plateau and Rockall Trough. The differences between the changes of SSH and SSHST in 1994-2010 can be due to the lack of in-situ temperature measurements prior to the start of the widespread use of Argo floats in the region. The spatial pattern and the sign of the observed SSHST change are determined by the SSHT change (Fig. 12c), which is partly balanced by the SSHS change (Fig. 12d). The net surface heat flux anomalies (QNET) can explain the SSH/SSHST increase in the Labrador Sea and over the Rockall Plateau and the SSH/SSHST decrease in the southern part of the Iceland Basin (Fig.  360 12e). The shortwave (QSR) and thermal (QTR) radiation anomalies largely compensated each other (Figs. 12f,g), and the largest contribution to QNET anomalies in 1994-2010 came from the sensible (QSH) and latent (QLH) heat flux anomalies. The observed warming in the Labrador Sea in 1994-2010 was mainly caused by the QSH anomaly. This implies that there were tendencies for the ocean to be colder than the air aloft and for the air temperature just above the surface to be increasing upward. This pattern is typical during periods with an anomalous heat flux into the ocean (positive QSH anomaly). The secondary contribution 365 to the Labrador Sea warming came from the positive QLH anomaly, also suggesting that the air temperature was higher than 13 the ocean temperature and the humidity was sufficient to cause condensation of water vapor on the surface of the ocean. This resulted in a loss of heat from air into the ocean (positive QLH anomaly). In the northern part of the Irminger Basin, the increase of SSH/SSHST in 1994 was associated with the concurrent heat loss to the atmosphere (Fig. 12e), which must have been compensated by an increased oceanic heat advection into the region. 370 The 2004-2014 period, depicted by the PC2 (red curve in Fig. 3) and the imaginary part of CPC1 (red curve in Fig. 9), was characterized by a decade-long decrease of SSH over most parts of the northeastern SPNA (Fig. 2b, solid red curves in Figs. 3 and 9), in particular in the Iceland and Irminger Basins (Fig. 13a). This decrease observed by satellite altimetry corresponds well to the decrease of SSHST (Fig. 13b), which was mainly determined by ocean cooling reflected in the SSHT change (Fig.   13c). The latter was partly compensated by a freshening induced SSHS increase in the eastern SPNA (Fig. 13d). The QNET 375 anomaly (Fig. 13e) can explain the 2004-2014 cooling and SSH decrease only in the eastern-northeastern SPNA. In the interior of the Labrador Sea, while the SSH was decreasing, the atmosphere was warming the ocean (Fig. 13e), which implies an increased oceanic heat transport out of the basin. The contributions of QSR and QTR anomalies to the QNET anomaly are rather small and compensate for each other (Figs. 13f,g). It is interesting to note that unlike during the 1994-2010 period (Figs. 12f,g) the contributions of QSH and QLH anomalies are geographically distinct (Figs. 13f,g). The QSH anomaly makes the largest 380 contribution to the positive QNET anomaly in the western SPNA (Fig. 13h), meaning that in this region the ocean temperature relative to the air temperature was colder than usual. This led to an anomalous heat flux from the atmosphere into the ocean.
In contrast, in the eastern-northeastern SPNA, the QLH anomaly was the largest contributor to the negative QNET anomaly (Fig.   13i). In the latter case, the ocean temperature relative to the air was warmer than usual, thus favoring evaporation and the associated upper-ocean cooling. 385 In 2011-2015, a tripole-related decrease of SSH (PC1 and the real part of CPC1 in Figs. 3 and 9) occurred over most parts of the SPNA, with the largest anomalies in the Irminger Basin and the Labrador Sea (Fig. 14a). This decrease was mainly steric in nature (Fig. 14b), and it was associated with strong cooling (Fig. 14c) and freshening (Fig. 14d). It is interesting to note that in most regions the QNET anomaly did not provide the largest contribution to the observed SSH change in 2011-2015 (Fig.   14d). Only its contributions over the northwest European shelf, in the Norwegian Sea, and along the East Greenland and 390 Labrador Currents were substantial. This means that the oceanic advection of heat and freshwater was the major driver for the 2011-2015 decrease of SSH in most parts of the SPNA. The distributions of ΔSSHT (Fig. 14c) and ΔSSHS (Fig. 14d) suggest that advection was mainly associated with the NAC. The structure of the QNET anomalies was mostly determined by the nearly equal contributions from QSH and QLH anomalies with the strongest impacts along the East Greenland and Labrador Currents (Figs. 14h,i). The QSR and QTR anomalies were somewhat smaller and largely compensating each other (Figs. 14f,g). We recall 395 that an exceptionally cold anomaly, termed the 'cold blob', occurred in the SPNA in 2015 (Ruiz-Barradas et al., 2018;Chafik et al., 2019). The analysis presented here shows that while the onset of this strong cooling started in 2004 in the eastern SPNA (red curve in Fig. 7a) and was partly driven by the negative QNET anomalies in 2004-2014 (Fig. 13e), the advection of colder and fresher water was apparently the main driver for the cold and fresh anomaly in -2015. This agrees with 14 a recent study by Holliday et al. (2020), who attributed the unprecedented freshening in the eastern SPNA in 2012-2016 to 400 large-scale changes in ocean circulation driven by atmospheric forcing.
During the following 2015-2019 period, SSH was rising in the northern and northeastern parts of the SPNA, particularly along the bottom topographic features associated with the Rockall Plateau and Reykjanes Ridge and along the Greenland shelf (Fig.   15a). Mostly, the observed SSH increase adequately compares to the SSHST increase (Fig. 15b). What is interesting to note is that this increase was not mainly determined by the SSHT change as in the previous time intervals. The upper-ocean warming 405 was leading to the SSH rise only in the Iceland Basin, over the Rockall Plateau, and further south upstream of the NAC (Fig.   15c). Warming in these areas was in part driven by the positive QNET anomalies (Fig. 15e) and by the advection of warmer and saltier water by the NAC (Figs. 15c,d). The positive QNET anomalies were also observed over the Reykjanes Ridge and in the Irminger Basin. However, these anomalies did not lead to the upper-ocean warming and they were accompanied by a negative SSHT tendency (compare Figs. 15e and 15c), which is suggestive of the oceanic heat flux out of these areas. An anomalous 410 heat loss to the atmosphere along the East Greenland Current (Fig. 15e), which was, however, associated with the local SSH increase, also suggests the dominant role of oceanic heat advection within the area influenced by the current. The largest contribution to the observed SSH rise in the interior of the Irminger Basin and Labrador Sea in 2015-2019 was provided by the upper-ocean freshening and the associated SSHS rise (Fig. 15d). Because the impact of surface freshwater flux (P+E; not shown) is negligible, the observed freshening could only be driven by the advection of freshwater and/or the continental runoff, 415 including the meltwater from Greenland glaciers. Similar to the previous time intervals, the contribution of the QSR and QTR anomalies to the QNET anomaly was small and balancing each other (Figs. 15f,g). The QNET anomaly was mainly driven by the QSH and QLH anomalies, with the largest contribution from the latter (Figs. 15h,i). Apparently, there was a tendency for the ocean to be colder than the atmosphere, thus, favoring anomalous heat flux from the air into the ocean.

The role of the large-scale ocean circulation in the SPNA 420
It is well known that persistent wind forcing leads to changes in the upper ocean heat content through Ekman dynamics and the adjustment of the large-scale geostrophic circulation. In Section 4.4, we showed that the observed westward propagation of sea level anomalies in the SPNA is related to the shifts of NAO-like surface pressure and wind anomaly patterns. In the previous section, we demonstrated that while some of the observed regional tendencies of sea level can be explained by surface heat fluxes, there are regional tendencies that can only be accounted for by the advection of heat and freshwater. For example, 425 the strong decrease of SSH/SSHST in 2011-2015 was mainly caused by the advection of colder and fresher water masses by the NAC (Fig. 14). In this section, we present the large-scale ocean circulation in the SPNA deduced from Argo trajectories at 1000-bar depth (Fig. 16a) and from eddy propagation velocities calculated from satellite altimetry data (Fig. 16b), and we qualitatively analyze how ocean circulation could contribute to the observed westward propagation of SSH anomalies during the observational period. 430

15
The time-mean circulation at 1000-dbar depth shows a distinct cyclonic pattern constrained by bottom topography (Fig. 16a).
This pattern is similar to the schematic upper-ocean circulation shown in Fig. 1 based on previous studies (e.g., Schmitz and McCartney, 1993;Schott and Brandt, 2007). The velocities derived from Argo trajectories well correspond to those derived from hydrographic sections and satellite altimetry (e.g., Sarafanov et al., 2012). The flow associated with the NAC, once it reaches the northern part of the Iceland Basin, recirculates southwestward along the eastern flank of the Reykjanes Ridge with 435 the speeds of about 5 cm s -1 . Upon crossing the Reykjanes Ridge, it follows its western flank, finally joining the system of western boundary currents composed of the East and West Greenland Currents and the Labrador Current, where the speeds often exceed 15 cm s -1 . The eddy propagation pattern and speeds (Fig. 16b) are similar to the 1000-dbar velocities, meaning that eddies follow the mean SPNA circulation pathways, steered by bottom topography. Because the eddy propagation velocities (Fig. 16b) are directly derived from the SSH anomalies, they are representative for property transports in the upper 440 ocean. Therefore, heat and freshwater signals advected to or generated in the Iceland Basin are transferred first to the Irminger Basin and then to the Labrador Sea. Assuming the average eddy propagation speed of 5 cm s -1 (Fig. 16b), it takes about 1.5 years for a SSH signal to propagate from the northern part of the Iceland Basin (~18°W, 62°N) to Cape Farewell and about 2 years to reach the northern part of the Labrador Sea at 63°N following the mean circulation pathway along the eastern and western flanks of the Reykjanes Ridge and the Greenland continental shelf. This is the same time scale as the one we identified 445 earlier using the CEOF analysis (see Section 4.3). Therefore, it is reasonable to suggest that the advection of heat and freshwater by the mean ocean circulation in the SPNA is a possible mechanism for the observed east-to-west propagation of SSH anomalies.
To further illustrate the role of advection, we present the time-depth diagrams of potential temperature and salinity anomalies from EN4 at three selected locations in the Iceland Basin (60°N, 20°W), in the Irminger Basin (60°N, 35°W), and in the 450 Labrador Sea (58°N, 50°W) (Fig. 17). It should be noted that the EN4 data is more reliable after the Argo array achieved its

Discussion and Conclusions
This study presents a reconsideration of the interannual-to-decadal SSH variability in the North Atlantic, the leading mode of which exhibits a tripole pattern (Volkov et al., 2019b). The tripole was originally detected by the conventional EOF analysis 465 (Figs. 2 a,c; blue curve in Fig. 3), which is effective at depicting only standing oscillations. The sign of the tripole is mainly determined by SSHT (Fig. 4a), partly balanced by a sizable contribution from SSHS (Fig. 4d). The analysis presented here demonstrates that the first EOF mode alone does not adequately represent the interannual-to-decadal sea level variability in the SPNA. We show that the first mode explains the majority (60-90%) of the interannual SSH variance only in the Irminger Basin and in the Labrador Sea (Fig. 5a), while the second EOF mode accounts for 60-80% of the interannual SSH variance in 470 the northeastern SPNA (Fig. 5b).
Furthermore, we demonstrate that the two modes do not represent two distinct physical processes. Instead, they belong to the same process and arise due to the general east-to-west propagation of SSH anomalies (Fig. 8). The CEOF analysis, which is designed to detect propagating (as opposed to standing) signals, yields the real and imaginary parts of the leading CEOF mode that correspond well to the first and the second EOF modes, respectively. This suggests that the two leading EOF modes evolve 475 as a quadrature pair associated with a propagation of SSH anomalies. Based on the CEOF analysis, there were almost three full cycles of the tripole-related SSH changes that exhibited westward propagation of SSH anomalies in the SPNA in 1993-2020 (Fig. 9). The reconstruction of the leading CEOF mode at different phases of one full cycle shows that SSH anomalies first appear along the northwestern European shelf and then gradually propagate westward (Fig. 8). It takes about 2 years for a signal to travel from the Iceland Basin to the Labrador Sea, and it takes 7-10 years for one full cycle to complete. 480 Because of the signal propagation, the concept of the North Atlantic SSH tripole introduced in Volkov et al. (2019) needs to be reconsidered, at least with respect to its application in the SPNA. While the standing oscillation of the interannual-todecadal SSH anomalies is a reasonable approximation of the variability in the subtropical and tropical North Atlantic, signal propagation needs to be accounted for in the SPNA. Therefore, the tripole in the SPNA needs to be either based on the leading CEOF mode or on the first two EOF modes combined. It is necessary to mention that in order to describe the variability of 485 SSH and ocean circulation in the SPNA, several authors have also used the subpolar gyre index, which, like the tripole, is based on an EOF decomposition of SSH fields (Häkkinen & Rhines, 2004;Hátún et al., 2005;Berx & Payne, 2017;Foukal & Lozier, 2017). The main difference between the subpolar gyre index and the tripole is that the global mean sea level is not subtracted from SSH fields prior to the computation of the former. Therefore, the subpolar gyre index exhibits a trend characteristic of the global mean sea level rise. The regional dynamic changes are then represented by higher modes, which 490 led Hátún & Chafik (2018) to justly conclude that PC2 in their calculation is a better metric for a gyre index than PC1. Our results imply that with the signal propagation, if the global mean sea level is not subtracted prior to the EOF analysis, even the third mode needs to be considered.
To expand on the earlier work, we analyzed what mechanisms were responsible for the observed interannual-to-decadal SSH changes in the SPNA in 1993-2020, and what mechanisms could be responsible for the observed signal propagation. It has 495 been documented that the North Atlantic SSH tripole is correlated with the NAO: stronger/weaker than average mid-latitude westerly winds associated with positive/negative NAO phases lead to divergence/convergence and lower/higher sea levels in the SPNA (Volkov et al., 2019). We find that since the two leading EOF modes of the low-pass filtered dynamic SSH depict the same physical process, the evolution of the second EOF mode is also related to an NAO-like dipole SLP pattern, but with shifted atmospheric pressure centers (Fig. 10). The definition of the NAO implies that it is a stationary standing oscillation 500 pattern. However, the subtropical high-and the subpolar low-pressure centers in the North Atlantic change both their intensity and position (e.g., Hameed and Piontkovski, 2004;Hameed et al., 2021). Therefore, we conclude that both the first and the second EOF modes may reflect oceanic response to the NAO-like persistent atmospheric forcing at different phases of its evolution. This conclusion is supported by the space and time evolution of SLP and wind anomaly patterns associated with the first CEOF mode of the low-frequency dynamic SSH (Figure 11). The observed propagation of SSH anomalies from the eastern 505 boundary towards the Labrador Sea corresponds to the westward shifts in atmospheric pressure and wind anomaly patterns (Figs. 8 and 11).
The role of surface buoyancy forcing over the SPNA in driving the interannual-to-decadal changes of SSH, SSHT, and SSHS is investigated over several characteristic time intervals : 1994: , 2004: -2014: -2015. The impact of surface freshwater fluxes is found to be negligible in all periods, so that any changes in SSHS are mainly driven by 510 the advection of freshwater. Advection was apparently the main driver for the strong upper-ocean freshening observed in the eastern SPNA after 2010 (Fig. 7a, Fig. 14d). We find that the net surface heat flux anomalies, mainly caused by the sensible and latent heat flux anomalies, can fully or partly explain some regional tendencies in SSHT. For example, QNET anomalies drove the increase of SSHT in the Labrador Sea in 1994 and in the Iceland Basin in 2015, and the decrease of SSHT in the northeastern SPNA in 2004-2014 and in the Labrador Sea in -2015. However, there are regions and time periods when changes in SSHT can only be explained by advection, the contribution of which cannot not be directly estimated from observations. These are the regions and periods, in which changes in SSHT are either much larger or have a sign opposite to the changes implied by the QNET anomalies. A prominent example is a strong cooling in the SPNA in 2011-2015, known as the 'cold blob', which coincided with contemporary freshening (Figs. 14 c,d). This period was apparently characterized by the advection of colder and fresher water masses into the region, consistent 520 with findings of Holliday et al. (2020).
In addition to shifting atmospheric pressure patterns, the observed westward propagation of SSH anomalies could be caused by the mean ocean circulation in the SPNA. It appears that while the overall propagation is westward, SSH anomalies associated with EOF1 and EOF2 first spread over the shallower areas in the east-northeast, including the currents along the eastern and western flanks of the Reykjanes Ridge and the East Greenland Current (EOF2; Figs. 2b and 5b), and then they 525 reach the deeper parts of the Irminger Basin and Labrador Sea (EOF1; Figs. 2a and 5a). The horizontal transfer of signals from the currents to the interior basins may be carried by eddies generated by the boundary currents. The likely role of ocean currents is qualitatively assessed using the climatological velocities at 1000-dbar obtained from Argo trajectories (Fig. 16a) and eddy propagation velocities estimated from satellite altimetry measurements (Fig. 16b). Both velocities depict the cyclonic circulation in the SPNA, constrained by bottom topography. The eddy propagation velocities signify the propagation of SSH 530 anomalies and, therefore, they are characteristic for the depth-integrated temperature and salinity (that define density and, consequently SSH) transports. We find that the time required for SSH anomalies to propagate from the Iceland Basin to the Labrador Sea (1-2 years) is consistent with the time implied by the velocity estimates and by the upper 2000-m temperature and salinity anomalies (Fig. 17). This means the observed westward propagation of SSH anomalies is partly due to the mean direction of ocean currents in the SPNA. Any anomaly generated locally by atmospheric forcing or advected from another 535 region is ultimately carried towards the Labrador Sea by ocean currents.
It should be noted that, due to geostrophy, both the SSH and the general ocean circulation are linked, and both adjust to persistent atmospheric forcing. For example, an increase of SSH along the European coast starts when the negative (cyclonic) SLP anomaly is centered over the eastern coast of Greenland and the atmospheric circulation near the eastern boundary is likely to cause downwelling (phase 0° in Figs. 8 and 11). As the cyclonic SLP anomaly weakens and moves towards the 540 Labrador Sea (phase 45° in Figs. 11), the subpolar gyre weakens and contracts and the positive SSH anomalies near the eastern boundary expand westward (phase 45° in Fig. 8). It has been reported that this situation can facilitate inter-gyre exchange (Häkkinen et al., 2011;Piecuch et al., 2017). Specifically, in response to a weakening of the subtropical high and subpolar low pressure centers, the subtropical and the subpolar gyres weaken, sea level decreases in the subtropical gyre and increases in the subpolar gyre, the subpolar front moves westward, and the eastern boundary region in the SPNA widens, entraining more 545 warm and saline waters from the subtropical gyre. Consequently, positive SSH anomalies emerge first near the eastern boundary of the SPNA and then expand westward as the subpolar gyre continues to weaken (phases 45° to 180° in Fig. 8). The opposite occurs when the subtropical and the subpolar gyres strengthen (phases -135° to 0° in Fig. 8). As demonstrated by the CPC1 (Fig. 9), the local maximum SSH anomalies occurred in the eastern SPNA around 1996SPNA around , 2004SPNA around , and 2009, and they reached the western SPNA 1-2 years later. The most recent increase of SSH in the eastern SPNA since 2014 and in the western 550 SPNA since 2016, that remains present in 2020, represents a recovery from an exceptional cooling and freshening that occurred in the SPNA in 2012-2016. This means that the recent conditions are favorable again for inter-gyre exchange.
Overall, we conclude that the observed interannual-to-decadal variability of SSH, including the westward propagation of SSH anomalies, is the result of a complex interplay between the local wind and surface buoyancy forcing, and the advection of properties by mean ocean currents. The relative contribution of each forcing term to the variability is space and time dependent 555 and, therefore, difficult to assess with available observations. The observed east-to-west propagation of SSH anomalies in the SPNA suggests the potential predictability of SSH changes and conditions favoring deep convection events in the region. This study puts the interannual-to-decadal changes of SSH in the SPNA in a broader context of the gyre-scale SSH variability in the entire North Atlantic. A caveat of this study is that it is based on a rather short observational record, which covers only a few cycles of the tripole-related variability. More observations are required to explore the persistence of the SSH signal 560 propagation in the SPNA and to assess the characteristic time scales. In the meantime, it is possible to use longer SST records in a follow-on study. Further research is also needed to investigate the details of the relationship between the shifting wind forcing patterns and the tripole evolution in the SPNA, potential oceanic feedbacks on the atmospheric circulation, and to relate the tripole-related changes to inter-gyre exchange and to the AMOC.