A shift in the ocean circulation has warmed the subpolar North Atlantic Ocean since 2016

The Subpolar North Atlantic is known for rapid reversals of decadal temperature trends, with ramifications encompassing the large-scale meridional overturning and gyre circulations, Arctic heat and mass balances, or extreme continental weather. Here, we combine datasets derived from sustained ocean observing systems (satellite and in situ), idealized observation-based modelling (advection-diffusion of a passive tracer), and a machine learning technique (ocean profile clustering) to document and explain the most-recent and ongoing cooling-to-warming transition of the Subpolar North Atlantic. Following a gradual cooling of the region that was persisting since 2006, a surface-intensified and large-scale warming sharply emerged in 2016 following an ocean circulation shift that enhanced the northeastward penetration of warm and saline waters from the western subtropics. The long ocean memory of the Subpolar North Atlantic implies that this advection-driven warming is likely to persist in the near-future with possible implications for the Atlantic multidecadal variability and its global impacts. The subpolar North Atlantic Ocean warmed sharply after 2016 as the transport of warm saline waters from the western subtropics was enhanced - following a period of cooling between 2006 and 2016 - according to analyses of observations via idealized modelling and machine-learning techniques.

T he Subpolar North Atlantic (SPNA) plays a fundamental role in the global ocean circulation and therefore in the Earth climate system. It provides a favorable environment for surface waters to mix and sink towards the seafloor before spreading south as part of the Meridional Overturning Circulation (MOC), which transports climate-relevant properties globally (e.g., heat, carbon) 1,2 . It is also the primary gateway for warm and salty Atlantic waters to reach the Nordic Seas and Arctic Ocean 3 , and therefore exerts a key control on the mass and heat budgets of those rapidly-changing and most-vulnerable areas [4][5][6] . Finally, it likely influences atmospheric regimes (e.g., storm tracks) and the development of continental heat waves 7 , and could be able to trigger multicentennial climate variability by governing complex processes within the coupled oceanatmosphere-ice system of the North Atlantic 8 .
The climatic imprint of the SPNA is evidently tightly linked to the temperature of its upper layer 9 , which is known to respond to both intermittent atmospheric events and decadal-scale changes in ocean circulation and associated heat transport [10][11][12] . The latter has notably been found responsible for a sudden warming-tocooling reversal of the SPNA in the mid-2000s 13 , with ramifications seen as far as in the sea-level of the East Atlantic shelves and coasts 14 . From an ad hoc observational reconstruction of the meridional heat transport at 45°N derived from the theory of surface-forced water mass transformation, the onset of a new advection-driven reversal in 2016 has been predicted, with suggestive temperature predictions for the early 2020s approaching their previous 2005-2006 maximum 15 . Here, we employ new observational and statistical approaches to document this recent and ongoing warming trend in the context of decadal-scale variability and demonstrate its dominant circulation-driven mechanism, namely the increased poleward transport of warm subtropical water.

Results
Entering a new warming phase. Our analysis targets the broad eastern portion of the SPNA since this is where significant decadal-scale temperature anomalies first transit before spreading cyclonically within the gyre circulation and invading the western portion of the domain 16 . An ocean analysis product covering the period 2002-2019 and primarily built with Argo-derived hydrographic profiles 17 is used to show the evolution of temperature anomalies in the eastern SPNA (10°W-35°W, 40°N-65°N), as averaged above a range of depth levels and relative to the 2002-2019 reference period (Fig. 1). On top of rapid intra-annual changes, the record is dominated by two trend reversals in 2006 and 2016. The former interrupts a decadal warming of the region triggered in the mid-1990's by a surge in the meridional oceanic heat transport at the entrance of the SPNA 18,19 . The following 10 years were marked by a gradual cooling and punctuated by a sharp drop in 2014-2015-the well-documented SPNA cold anomaly partly driven by anomalously strong ocean heat loss to the atmosphere 10 . Since 2016, however, the temperature tendency is seen to have changed sign again suggesting that the SPNA has most likely entered a new warming phase, i.e., the main focus of this study.
The decadal temperature reversal in the SPNA (Fig. 1) is surface-intensified but is still apparent when averaging the temperature of the water column down to 2000 m depth. They have induced concomitant and significant changes in sea surface height (SSH), as illustrated here by a strong agreement on intraannual to decadal timescales between temperature anomalies and altimetry-derived SSH anomalies (dashed line in Fig. 1; r = 0.82 with 0-100 temperature time series). Overall, the ocean heat content (OHC) of the 0-2000 m layer increased by about 8.5 × 10 21 J between April 2016 and December 2019, with 80% of this heat gain found above 700 m depth. Anomalous air-sea heat fluxes during this time-span suggests that the recovery from an intermittent period of extreme ocean heat loss in 2014-2015 10 potentially contributed to about a quarter of that value (Supplementary Figure 1). Thus, this leaves ocean heat transport convergence as a dominant contributor, and the remaining of the paper focuses on the dynamical mechanism underlying this specific contribution.
Basin-scale warming and shifting currents. The horizontal pattern of upper (0-100 m) temperature changes in the northern North Atlantic from the cold 2014-2016 period to the warm 2017-2019 period shows a large area of positive anomalies with a maximum magnitude of 0.6°C and covering the Iceland Basin, the Iberian plain, and the eastern Newfoundland basin (Fig. 2a). The concomitant distribution of air-sea heat flux anomalies confirms their minor (yet non-negligible) contribution to the warming trend (Supplementary Figure 1). As expected, the distribution of SSH anomalies shows a consistent pattern (Fig. 2b), with positive anomalies covering most of the eastern SPNA, which is where the largest thermosteric contributions to SSH changes are known to occur 20 . Those recent SSH changes are particularly intense in the Iceland Basin where the northern and strongest branch of the North Atlantic Current (NAC)-the subpolar front-is observed 21 (see Fig. 2d for a schematic of the circulation features relevant to the present study).
This particular SSH distribution and the strong correlation between upper temperature and SSH tend to support the advective origin of the warming, and a first glimpse into concomitant changes in the intensity and pathway of surface currents is provided by the mean kinetic energy (MKE) calculated from the altimetry-derived geostrophic velocity field (Fig. 2c). A striking intensification of the subtropical portion of the NAC, also visible as an increased SSH cross-stream gradient (Fig. 2b), is observed south of 45°N and west of 30°W. The large MKE positive anomalies in the western subtropics can be tracked downstream within the northern branch of the NAC (circa 50°N), while negative anomalies occupy its southern portion. This overall MKE pattern is in good agreement with Lagrangian model experiments showing that stronger subtropical inflows within the NAC drive an increased (decreased) transport of its northern (southern) branch in the eastern SPNA 22 . This observed seesaw of the NAC branches in the eastern SPNA is often described as "contraction" and "expansion" of the subpolar gyre (SPG) and subtropical gyre (STG). The associated changes in the relative proportion of warm waters (i.e., from the western subtropical gyre) and cold waters (i.e., from the western subpolar gyre) in the eastern SPNA are likely to drive significant temperature changes 19,23-25 (see Fig. 2d for a schematic of the source water pathways feeding the eastern SPNA).
Passive tracing of advective origins. To evidence and quantify the potential link between the observed shift in circulation and the 2016-to-present warming of the SPNA, we first turn to a simple advection-diffusion model for a passive tracer of local concentration C: This simple model enables to extract from surface Eulerian velocities and mesoscale eddy diffusivities a Lagrangian decomposition of the eastern SPNA water mass origins (see schematic in Fig. 2d). Two separate experiments are conducted, in which tracers are kept constant (i.e., C = 1) in either a western subtropical box (Fig. 3a) or a western subpolar box (Fig. 3b). They are allowed to evolve freely from those source regions for N years within a N-year mean altimetry-derived surface geostrophic current field v and a background surface eddy field parametrized with an observation-based estimate of mesoscale eddy diffusivities κ 26 . We will refer to the final tracer concentration in these two experiments as C STG and C SPG , respectively. The experiments are run for successive and overlapping N-year windows so that the final distribution of C associated with year y is the result of N years of tracer advection within the altimetry-derived velocity field hvi ¼ 1 N R t¼y t¼yÀNþ1 vdt. Dependency to the parameter N, which represents a typical advective timescale of surface waters between the two source regions and the eastern SPNA, is evaluated through an ensemble mean and an ensemble spread (standard deviations) of model outputs with N = 2, 3, and 4 years (Methods and Supplementary Fig. 4). The choice of N was guided by the literature 27,28 and validated herein via an additional tracer experiment that quantified the mean travel-time of tracers between the two source regions and the eastern SPNA ( Supplementary Fig. 2).
Noteworthy, the non-inclusion of surface Ekman velocitieswhich should, in reality, hamper the surface exchange of C STG across the inter-gyre boundary 27,29,30 -enables the tracers to be strictly advected within a velocity field that primarily reflects the surface-intensified first baroclinic mode and the motions of the main thermocline 31 . The altimetry-derived geostrophic surface circulation used herein is hence likely to be a fair proxy of the large-scale upper circulation and its changing intensity and horizontal structure 32 . The consistency between geostrophic surface and interior pathways will be demonstrated herein with Current (NAC) that feeds the eastern SPNA with a mixture of cold subpolar waters (blue -C SPG ) and warm subtropical waters (red, C STG ). The relative proportion of these two source waters, denoted P STG in the present study and represented by the color scale here, is the key quantity behind temperature trend reversals in the eastern SPNA (green box as considered here).
a companion tracer experiment run on the three-dimensional velocity field of an ocean reanalysis (ECCOv4r4 33 ), as well as with a statistical clustering of vertical ocean properties.
The long-term mean (2003-2019) final distributions of C STG and C SPG show the general pathways of the two main source waters composing the inter-gyre flow (Fig. 3a, b). From those, the local proportion of western-origin subtropical water P STG ¼ C STG C STG þC SPG (and western-origin subpolar water P SPG ¼ 1 À P STG ) can be deduced (Fig. 3c). The Lagrangian nature of P STG provides an alternative view of circulation contours and gyre shapes usually derived directly from SSH patterns 34,35 . More importantly, it here enables a direct link with the pathways and upstream properties of both source waters advected by the NAC along the intergyre boundary. The time-mean proportion of subtropical-origin and subpolar-origin surface waters in the eastern SPNA box amounts to 60% and 40%, respectively (in line with Lagrangian decomposition of the eastern SPNA thermocline waters in a hindcast simulation 22 ). We note here that the contribution of parametrized eddy mixing to the distributions of C STG and C SPG is as expected substantially smaller than the advective contribution and has the average effect of decreasing P STG in the eastern SPNA box by about 4% ( Supplementary  Fig. 3).
Running the Lagrangian experiment with successive (and overlapping) N-year averaged velocities enables to construct a time series of P STG within the eastern SPNA and to compare it with a time series of temperature anomalies of the upper layer (Fig. 4a). The striking correlation between those two quantities, with the variability of P STG capturing both trend reversals circa 2006 and 2016, suggests a clear causal relationship (estimated here as 0.05°C/%). The respective evolutions of C STG and C SPG reveal the important role of the former in driving P STG variability, and most particularly in driving the latest and abrupt cooling-towarming reversal. There is likewise a good agreement between the horizontal patterns of temperature and P STG increase ( Fig. 2a and Fig. 3d, respectively), with positive anomalies covering most of the eastern SPNA. Negative P STG anomalies in the central and eastern subtropics further highlight the recent enhanced penetration of warm subtropical waters 36 , which, in turn, agrees with the observed intensification of the northern NAC branch (Fig. 2c). We also note a consistent northward displacement of the 50% P STG isoline-an indicative measure of the gyre expansions and contractions-from its 2015 southernmost position, with the C STG -dominated region currently expanding into the Iceland Basin ( Supplementary Fig. 5). Overall, all important changes revealed by the simple model are, as expected, associated with horizontal advection 11,19,24,25,37 (thin lines in Fig. 4a; note that because κ is time-invariant, the time-varying diffusive signal reflects changes in large-scale advective pathways only).
Following our previous discussion on the representativeness and use of surface (geostrophic) velocities to diagnose changes in interior pathways, we evaluate a companion experiment performed with the three-dimensional velocity field of the ECCOv4r4 ocean state estimate 33 . An identical model setup (with the advective scheme augmented by a vertical component) is used to infer the mean and anomalous concentration of C STG and C SPG from the surface to about 700 m depth (recall that decadal temperature signal we focus on is largely contained above 700 m depth). Our purpose here is strictly to demonstrate that the largescale geostrophic pathways are largely coherent from the surface down to the main pycnocline, using an ocean model that should best reconcile the surface altimetry fields with interior observational constraints-i.e., we do not intend to reproduce dedicated model-based analysis of those pathways, which are published elsewhere, and for which higher-resolution simulations could be preferred. We first note that the time-mean horizontal distributions of C STG and C SPG appear similar to the altimetry-derived fields but also highlight some discrepancies associated with wellknown deficiencies of low-resolution numerical models in the region (e.g., note the absence of a proper Northwest Corner at the intergyre boundary, Supplementary Fig. 6). In line with aforementioned works, the northward penetration of C STG is maximum below the Ekman layer (peaking at 100 m), while the eastward penetration of subpolar tracer is surface-intensified (i.e., Ekman drift is enhancing the eastward tracer advection from the Labrador Sea). Most importantly, the (decadal) changes in C STG and C SPG are vertically-coherent between the base of the Ekman layer and the depth of the main pycnocline (taken here to be σ 0 = 27.6), i.e., within the natural layer along which tracers can efficiently spread and mix from both source regions (Supplementary Fig. 6). Because geostrophic currents are not likely to significantly change from the surface to the base of the Ekman layer, we conclude that altimetry-derived surface flows allow to capture large-scale changes in the interior circulation (in line with ref. 31 ) and in the associated propagation of the two primary water masses that compose the upper layer of the eastern SPNA.
To better quantify the relative importance of P STG (or P SPG ) changes in driving the observed temperature signal in the eastern SPNA, we now associate those changes with the initial temperature of both source waters θ STG and θ SPG as averaged within their respective source regions (boxes in Fig. 3a, b) using the following equation: This reconstruction enables to decompose the advection-driven temperature anomalies θ 0 ADV of the eastern SPNA into three distinct contributions and to evaluate their relative importance for the heat budget of the region. In line with the tracer experiment configuration, we assume a typical timescale of N years for temperature anomalies to be advected adiabatically above the main pycnocline (σ 0 < 27.6) from their source regions to the eastern SPNA where they mix into a single water mass (see Methods). Three distinct drivers of temperature variability are thereby considered: anomalies in the relative proportion of the source waters (first term in r.h.s. of Eq. 2) and anomalies in their respective initial temperature (second and third terms). It is found that the first of these contributions largely dominates the others, confirming that shifts in the large-scale circulation projecting onto the time-mean ocean thermal structure is the primary mechanism triggering basin-scale reversals of temperature trends in the SPNA (Fig. 4b). The advection of remote temperature anomalies from both source regions remain nonnegligible and consistently tend to oppose themselves-most especially since 2013-in line with the well-known dipole pattern of temperature anomalies between the western subtropical and subpolar regions [38][39][40] .
Statistical clustering of water masses. An alternative confirmation of the importance of circulation shifts in driving the recent cooling-to-warming reversal in the eastern SPNA is provided by applying a machine learning algorithm-namely a profile classification model (PCM)-to our three-dimensional observational dataset 41 . The PCM is built upon a synthetic statistical and unsupervised classification of ocean vertical profiles into a chosen number of "clusters" or "classes" using the classical "Gaussian Mixture Model" method. Once trained, the PCM can be used to assign any given profile to its most probable class (see Methods). It is herein used to reveal coherent temporal and spatial patterns in the SPNA and to further verify whether the altimetry-derived results (i.e., surface-focused advective shifts) are representative of the ocean interior. In other words, we seek to complement the Lagrangian clustering of surface pathways with a statistical clustering of interior ocean properties, which in the present case is an indirect indicator of anomalous three-dimensional pathways (i.e., anomalous penetration of stratified subtropical waters into the vertically-homogeneous subpolar gyre).
As pointed out repeatedly in PCM analysis [41][42][43] , the determination of the most appropriate number of clusters can only be indicatively given by statistical metrics (e.g., Bayesian Information Criteria 44 ) and it is the scientific meaning of the clusters for the given purpose that must prevail. Because our present goal is to mimic the altimetry-based and objectively distinguish the typical vertical structures of (homogeneous) subpolar and (stratified) subtropical water masses, we employ a simple approach and train a 2-class PCM using the 2002-2019 temperature and salinity gridded product (see Methods). The reader is referred to Supplementary Fig. 7 for the presentation of a 3-class PCM.
The horizontal distribution of the subpolar and subtropical classes, referred to as K STG and K SPG hereafter, is shown in Fig. 5a along with their respective temperature and salinity vertical structures and associated spread in Fig. 5b. One will easily note the close resemblance with the altimetry-derived distribution of P STG previously obtained (Fig. 3c). Here, the PCM provides the missing information on the typical vertical structure of both "water masses", with subpolar waters being logically colder, fresher, and less stratified than their subtropical neighbors. This trained PCM is finally used to classify each year of the data set and to reveal the role of both classes in explaining the observed temperature variability. This role is quantified by computing the respective contribution of each class K to the total OHC of the 0-700 m eastern SPNA as: with p K being the time-varying local probability that a given θ profile belongs to the considered class K (i.e., K STG or K SPG ) 45 .
The resulting time series concur with the tracer-based monitoring of P STG : since the mid-2000s and until the mid 2010s, in the eastern SPNA the relative proportion of K STG decreased and the proportion of K SPG increased, so did their associated OHC contributions (Fig. 5c). From the mid-2010s, trends reversed and K STG and K SPG respectively started to increase and decrease their OHC contributions. One will note that this recent decrease in OHC KSPG seemingly preceded that of K SPG proportion, which is in line with the local atmospheric drivers (air-sea heat fluxes) of the subpolar cold anomaly that occurred during 2014-2015. Because of their different mean background temperature profile (Fig. 5b), OHC KSTG overcomes OHC KSPG overall and accounts for both trend reversals in 2006 and 2016. Therefore, both surface-based index (Fig. 3) and stratification-based index (Fig. 5) of source water pathways concur with one another in pointing out a significant circulation shift and increased penetration of warm subtropical waters into the eastern SPNA since 2016.

Conclusion and discussion
This study documents an ongoing large-scale hydrographic change occurring within the upper layer of the SPNA: the decadelong cooling trend prevailing since 2006 rapidly reversed in 2016 and a large-scale warming settled in the central and eastern portions of the domain. The responsible mechanism-also at play for the previous warming-to-cooling reversal in 2006-was revealed with observational datasets by solving a simple altimetry-based advective model for a passive tracer and by performing a statistical clustering of ocean in situ profiles. This combined approach provided two independent and observation-based decompositions of eastern SPNA waters into subtropical-origin and subpolar-origin waters, and both depicted the changing relative proportion of these two source waters as a predominant cause for the 2016 trend reversal. Confirming previous model-based analyses of the eastern SPNA heat budget 19,25,46 , this change was more specifically induced by a rapid increase in the northeastward transport of warm subtropical waters within the North Atlantic Current, in line with a concomitant strengthening of the MOC heat transport as detected across 45°N 15 and presumably pending at 26°N 47 .
Embedded within the large-scale cyclonic gyre circulation, this warming signal is expected to spread westward from its current eastern SPNA location and progressively invade the Irminger and Labrador Sea while propagating to greater depth via boundaryfocused mixing and sinking 48 . This warming signal may also influence the intensity of light-to-dense water mass transformation and the depth of open-ocean convective mixing in those regions, and hence the forthcoming magnitude of the Atlantic MOC and related property transports. Ongoing investigations are specifically dedicated to this important issue. The variability of the SPNA heat content is known to exert a key influence on the so-called Atlantic Multidecadal Variability (AMV)-index based on sea surface temperature-that is generally used to characterize the climate of the North Atlantic and its global-scale ramifications [49][50][51][52] . In fact, it was recently suggested that the AMV would soon transition towards a negative (i.e., cold) phase owing to the cooling of the SPNA since the mid-2000s 53 . It is too early to ascertain the low-frequency nature of the ongoing SPNA warming described herein, and whether this thermal reversal will ultimately damp the current downturn of the AMO index remains as of now an open question 54 (see Supplementary Fig. 8 for an updated AMV index). We argue, however, that the advective origin of the ongoing SPNA warming and the long ocean memory of this region 55 makes it likely to persist over the coming years. Several indicators in fact describe oceanic and atmospheric conditions similar to those of the mid-1990's when another sharp cooling-to-warming reversal triggered a decadelong warming trend 18 . Those indicators include the preceding persistence of strongly positive North Atlantic Oscillation conditions (see Supplementary Fig. 8), the occurrence of relatively deep convective activity 56,57 and intense water mass transformation in the SPNA in the preceding winters 28 , and a resulting intensification of the MOC and its meridional heat transport across the southern SPNA boundary 15 . Investigating and separating the atmospheric drivers of the ocean-focused changes and mechanisms reported herein will motivate a follow-up study.
Methods ISAS-15 temperature and salinity gridded fields. The observational in situ dataset used is the In Situ Analysis System, which provides gridded threedimensional fields of temperature and salinity during 2002-2019 58 . Originally based on the Argo network of profiling floats, ISAS is now optimally interpolating all types of vertical profile (e.g., XBTs, CTDs) and time series. The version used here is ISAS15 and is available online 17 . Data were quality-controlled in delayed mode for the years 2002-2016 and in near-real-time mode for 2017-2019, with this difference in data treatment having no impacts on the large-scale temperature anomalies investigated herein. The field of eddy diffusivities was computed by 26 from observational data and theory and is available online at https://figshare.com/ articles/Groeskamp_et_al_2020_-_mixing_diffusivities/12554555. Satellite altimetry. The dataset for the sea surface height (SSH) and associated geostrophic velocities used to derive the mean kinetic energy (MKE) is the DUACS DT2018 multimission satellite altimetry 59 . The data is on a 1/4°grid resolution and spans the period from January 1993 until December 2019. Since the delayed time data is only available until May 2019, we extended the data until December 2019 using with the near-real-time fields.
NCEP atmospheric fluxes. The atmospheric fluxes used here are derived from the 1948 to 2019 NCEP/NCAR reanalysis 60 . The net heat flux was derived by adding the shortwave, longwave, sensible and latent heat flux. For the momentum fluxes, we used the zonal and meridional wind stress to derive the wind-stress curl. The horizontal grid resolution of the data is 2.5°× 2.5°.
ECCOv4r4 ocean state estimate. The state-of-the-art ocean state estimates Estimating the Circulation and Climate of the Ocean (ECCO) are based on leastsquares fitting of the Massachusetts Institute of Technology General Circulation Model to numerous types of global-scale ocean observations during 1992-2017. ECCO is thoroughly described in the literature 33 and was shown to be a good fit to observations of global-scale and basin-scale ocean heat content, especially in analysis of the North Atlantic heat budget 11,48 . We here use the monthly interpolated field of the ECCOv4 Release 4, which were downloaded from https://ecco. jpl.nasa.gov/drive/files/Version4/Release4. Tracer advection and diffusion model. Equation 1 is solved on a subset of the ¼°c artesian AVISO grid using a forward timestep and an upstream differencing scheme for the advective terms. An explicit 2-time level scheme is used with a timestep of 1500 s. The open boundary conditions for the overall domain are set as ∂C ∂x ¼ 0 and ∂C ∂y ¼ 0. The diffusive term is built upon a recent estimate of the three-dimensional mesoscale eddy diffusivities from observation data and theory 26 . Two types of experiments are conducted in which the initial tracer concentration is set to C=0 everywhere except in a western subtropical box (C STG = 1; 60°W-75°W, 35°N-45°N) or a western subpolar box (C SPG = 1; 50°W -60°W, 55°N -65°N). We note that the present definitions of the subtropical and subpolar boxes encompass the bulk of the western boundary currents of both gyres (i.e., the Gulf Stream and the West Greenland and Labrador Currents, respectively) that connect with each other as the North Atlantic current downstream near Flemish Cap and the Grand Banks of Newfoundland. In the case of the subtropical box, we hence consider the Northern Recirculation Gyre 61 as part of the subtropical source region. The tracer is then allowed to spread freely during N years within a N-year mean surface geostrophic velocity fields obtained from altimetry data while being kept constant within the source regions. The experiments are run for successive and overlapping N-year windows so that the final distribution C associated with year y is the result of N years of tracer advection within the velocity field v h i ¼ 1 N R t¼y t¼yÀNþ1 vdt. In order to ensure the robustness of the method, the time series presented in the paper are ensemble mean of three experiments using different (yet realistic) typical surface advective time scales N = 2, N = 3, and N = 4 years. Those values were deduced from a separate tracer experiment in which tracer was added once in the western STG and western SPG boxes and advected freely froward during 10 years within the time-mean (2002-2019) velocity field. The tracer concentration was averaged within the eastern SPNA box at each time step of the model integration, providing a precise quantification of the advective timescale from both region ( Supplementary Fig. 2). The spread (i.e., standard deviation) between the results of the experiment is used as an uncertainty range in Fig. 4 (see Supplementary Fig. 4 for the individual time series from these experiments). An identical model setup was used in the ECCOv4r4 tracer experiment, with the vertical component of advection added to the advective scheme, and the tracer concentration set to unity over the full vertical extent of the source regions (C STG = 1; 60°W-75°W, 35°N-45°N, surface-bottom) or a western subpolar box (C SPG = 1; 50°W-60°W, 55°N-65°N, surface-bottom). Final tracer concentrations were eventually averaged above the main pycnocline, taken here to σ 0 = 27.6 kg m −3 .
The local proportion of western-origin subtropical water P STG ¼ C STG C STG þC SPG (or equivalently western-origin subpolar water P SPG ¼ 1 À P STG ) in the SPNA is then calculated. It is used to estimate the three primary advective drivers of temperature changes in the eastern SPNA using Eq. 2. θ STG and θ SPG represent the time-mean initial temperature of subtropical and subpolar source waters as averaged above σ 0 = 27.6 kg m −3 . This isopycnal is chosen here to represent the main pycnocline of the domain above which the two source waters spread downstream (several tests were performed with σ 0 ranging between 27.1 and 27.7 and limited impacts on the results were found). θ 0 STG and θ 0 STG represent departure from θ STG and θ SPG . Their respective time series are out-phased by N years in Eq. 2 to account for the appropriate advection timescales of anomalies from their source regions to the eastern SPNA.
Profile classification model. A PCM with K = 2 classes was trained on vertical temperature and salinity profiles from 2002 to 2019 time mean of the ISAS dataset (described above). Only profiles where depth does exceed 2000 m were included in the PCM. Each variable was first normalized by vertical levels (e.g., all-temperature values at 200 m are centered and standardized, this allows for deep structures to be preserved against the more diverse surface values) and secondly compressed using a principal component analysis (to reduce the number of vertical dimensions to a few eigen values). Then a Gaussian Mixture Model with 2 modes is fitted on reduced profiles using a classic Expectation-Maximization algorithm 62 (See 41 for full details of the PCM methodology and 43 for the presentation of a PCM with two variables). Finally, temperature and salinity profiles of the full ISAS15 time series are classified using the fitted PCM, i.e., the p K posterior probabilities of the Gaussian Mixture components (clusters) are computed at each grid points. A PCM is a fuzzy classifier because the p K can vary from 0 to 1 and their sum for each profile is 1. This means that each grid point is classified and can possibly be in a smooth intermediate position between the cluster core structure. This also means that the grid cell heat content is preserved if split into each cluster (see Eq. 3). It is important to note that the PCM training does not consider the latitude/longitude of profiles. Therefore, the cluster distribution in space naturally arises from similarities and coherence in ocean interior structures. To compute the typical vertical structures of clusters, we used the mean and standard deviation of profiles attributed to a given class. The PCM analysis was conducted using the pyXpcm software 63 .

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Code availability
The codes that support the findings of this study are available from the corresponding author upon reasonable request.