The Major Role of Air-Sea Heat Fluxes in Driving Interannual Variations of Gulf Stream Transport

The Gulf Stream (GS) is central to the global redistribution of heat due to the transport of large volumes of warm water from the tropics to high latitudes and the extreme ocean heat loss to the atmosphere. This study assesses the extent to which winter surface heat fluxes and wind stress curl can drive interannual variations of the full-depth GS transport. Intensification of the has been observed immediately after a of frequent cold air outbreaks that led to a deepening of the mixed layer and subsequent steepening of meridional temperature gradients to the south of the GS. This forcing process is further investigated here using the ORCA12 hindcast (1978–2010) from the global, eddy-resolving, ocean-only Nucleus for European Modeling of the Ocean model in order to understand GS forcing mechanisms. Lagrangian analysis is also undertaken to examine the impact on the southern recirculation gyre. Results show that surface heat fluxes and wind stress curl can act in concert to effect year-to-year changes of up to 38% of the spring GS transport at a of from and southern recirculation, which can also limit the formation of deeper layers to the south of the GS

Strong heat losses, and associated buoyancy forcing, are often induced during widespread outbreaks of cold, continental air from the northeast of North America extending over the GS region (Bane & Osgood, 1989;Grossman & Betts, 1990;Joyce et al., 2009;Kelly et al., 2010;Ma et al., 2015). These losses are often the largest seen globally with extreme values sometimes exceeding 1,000 Wm −2 during an outbreak (Silverthorne & Toole, 2013). This surface heat loss initiates strong convection in the Sargasso Sea, which can lead to mixed layer depths (MLD) of up to 500 m (Buckley et al., 2014;Rossby, 1999). Re-stratification then preserves beneath the surface the 18 Degree Water formed in the early spring (Leetmaa, 1977;Mensa et al., 2013;Silverthorne & Toole, 2013). This deepening of the thermocline combined with an increase in meridional temperature gradients across the GS will have a significant influence on its baroclinic transport, given the thermal wind relation (Bane & Osgood, 1989;Kelly et al., 1996;Leetmaa, 1977;Sato & Rossby, 1995;S. Dong et al., 2007), see Section 2.2.
There is also a view that GS variability is primarily associated with large wintertime air-sea momentum fluxes (Kelly et al., 1999;Xue et al., 1995;Xue & Bane, 1997). In particular, the role of wind stress curl (WSC) variability has been found to drive changes to the transport of the GS (e.g., Anderson & Corry, 1985;De Coetlogon et al., 2006;Di Nezio et al., 2009;N. G. Hogg & Johns, 1995;Sturges & Hong, 2001). For example, a stronger (i.e., more negative) WSC in the subtropical gyre (STG) causes an intensification to the gyre, and therefore the GS (A. H. Chaudhuri et al., 2011;Hakkinen et al., 2011). WSC variability is also closely linked to the North Atlantic Oscillation with the positive phase causing a more negative WSC and stronger GS transport (DiNezio et al., 2009;Meinen et al., 2010). Yet, De Coetlogon et al. (2006) emphasized that wind forcing cannot be the sole influence on GS transport, and by implication that an important role is also played by buoyancy forcing. The barotropic transport is primarily associated with the large-scale wind driven circulation while baroclinic transport variability is related to buoyancy forcing that causes steric changes to the thermal structure of the water (Gill & Niiler, 1973;Wang & Koblinsky, 1996).
Recognizing the importance of buoyancy forcing for the wider GS system, L. V. Worthington (1972) proposed a mechanism termed "anticyclogenesis," whereby the STG is strengthened via a deepening of the thermocline on the southern flank of the GS. The clearest example of anticyclogenesis was observed by L. V. Worthington (1977) after the severe winter of 1976/1977, when cold air outbreaks led to large heat losses over the western subtropical gyre. This led to anomalously deep mixed layers and a high renewal of mode water to the south of the GS, which strengthened the meridional temperature gradients across the current. He found that this increased baroclinic GS transport (relative to 2,000 m) by around 50% compared to corresponding transport estimates after the much milder winters of 1974/1975and 1975/1976(L. V. Worthington, 1977. The effects of winter 1976-1977 were also investigated by Leetmaa (1977) who revealed a thermocline depth of ∼800 m at 34°N, 71°W, which is close to where L. V. Worthington (1977) measured an increase in transport. The extensive renewal of 18 degree water during this winter followed a series of years when its formation was anomalously low, coincident with mild winters over northeast North America (Worthington, 1977). This severe winter, therefore, caused an increased transport, compared to the previous decade, of ∼40 Sv (Worthington, 1977). In addition to the strengthening of the main current (of up to 60%) following winter cooling, Huang (1990) found an accompanying increase in both the recirculation gyres ( Figure 1) using a two-layer model. However, a study by Zheng et al. (1984) proposed that cooling over the Sargasso Sea is not the sole driving force for an increase in GS transport. They suggest that this mechanism is reinforced by the heat loss observed in the Slope Water, that is, the combination of isotherms upwelling to the north of the GS and downwelling to the south. This study will explore the extent to which winter surface heat fluxes have led to spring GS intensification over the last three decades and in particular to examine the contributions coming from the north, for example, Zheng et al. (1984), and from the south, for example, L. V. Worthington (1977) of the GS core. The analysis is conducted using the hindcast of a 1/12° global ocean model forced using surface observations over the period 1978-2010. With continuous high-resolution coverage in both time and space, multitimescale GS variability is more representative than previously available in sparse observations or coarser resolution models. This permits a more thorough investigation of GS forcing mechanisms. We begin by outlining the model and diagnostic techniques. We then present findings that highlight the relative importance of both mechanisms outlined above (L. V. Worthington [1977]; Zheng et al. [1984]), and additionally propose a mechanism that suggests westward intensification of the SRG.

Model Hindcast
We employ the high resolution, global, ocean model Nucleus for European Modeling of the Ocean (NEMO, Madec, 2015), using the eddy-resolving ORCA12 configuration setup of the DRAKKAR project (e.g., Blaker et al., 2015;Duchez et al., 2014;Marzocchi et al., 2015), which spans 1978-2010. This hindcast has a horizontal resolution of approximately 1/12° and 75 vertical levels with finer grid spacing near the surface. Sea ice is represented by the Louvain la Neuve (LIM2) sea ice model (Timmerman et al., 2005). In order to prevent numerical instability where the meridians converge in the northern hemisphere, the model is tri-polar with two North Poles (Madec & Imbard, 1996). Additionally, the bottom topography is represented as partial steps with a free-slip lateral friction condition applied at the lateral boundaries. The run is initialized using the World Ocean Atlas 2005 climatological fields Locarnini et al., 2006). It uses Drakkar Forcing Set 4.1 surface forcing fields produced by DRAKKAR (Brodeau et al., 2010), with 6-hourly mean 10 m winds, 2 m temperature, 2 m humidity, daily mean heat and radiative fluxes and monthly mean precipitation fields adjusted from ERA40 reanalysis (Uppala et al., 2005). Bulk formule are used to calculate the turbulent momentum, heat and freshwater fluxes from the surface atmospheric conditions and the sea surface temperature and sea ice, while the salinity is relaxed toward the climatology to avoid excessive drifting. The output fields are stored as 5-day means but for the purpose of this study we use monthly means, to reduce the signal of chaotic mesoscale variability.
The realism of the modeled surface velocities with the altimetric sea surface currents (1/4° resolution) distributed by AVISO (http://www.aviso.altimetry.fr/duacs/) (AVISO, 2014) is illustrated in Figure 2. There is strong qualitative agreement between the modeled and observed surface circulation on decadal timescales (Figures 2a and 2b) fields with high velocities seen in the GS region and an accurate GS separation. The model also represents realistic levels of mesoscale variability (in the form of meanders and eddies) in the monthly means (Figures 2c and 2d), and a broadly correct location and strength for the North Atlantic Current and Azores Current. This is further illustrated by Marzocchi et al. (2015).
We note that the model successfully represents GS separation from the coast throughout the hindcast, which is a common problem in models of resolution lower than 1/10° (Chassignet & Garraffo, 2001). These coarser resolution models tend to have a GS that separates from the coast too far north (e.g., Dengg et al., 1996). However, the sensitivity to subgrid scale parameterizations continue to make realistic GS separation a challenge (Chassignet & Marshall, 2008). This is discussed extensively in the review chapter by Chassignet and Marshall (2008) who considered factors such as topography, boundary conditions, viscosity, interaction with the Deep Western Boundary Current.

Diagnostics, Datasets, and Transport Calculations
From the ORCA12 hindcast, the following fields are analyzed: mean winter (January-March [JFM]) surface heat flux (Q net ), which is diagnosed from the model surface fields, winter mean WSC and early spring (April) MLD, defined as the depth where the temperature differs by 0.5°C from the surface over the domain 20-50°N, 85-30°W in the North Atlantic. Early spring is chosen as the ocean will respond to the culmination of winter forcing. Additional to the model forcing, we analyze the winter mean sea level pressure field, 2 m air temperature, and 2 m wind speed from the National Centres for Environmental Prediction and National Centre for Atmospheric Research (NCEP/NCAR) atmospheric reanalysis (2.5° × 2.5° horizontal resolution), all over the same period as the ORCA12 hindcast   (Kalnay et al., 1996).
The temperature and meridional temperature gradient profiles in April are examined along a transect at 70°W (Figure 1). The GS still has a narrow core at 70°W, which makes it ideal to measure the geostrophic transport with minimal mesoscale interference, for example, from meanders or eddies, which are commonplace farther east. It is also close to where L. V. Worthington (1977) found an intensification of the current after the severe winter of 1976/77. Using the model eliminates the need to employ a reference level for a depth of no motion as the sea surface height (SSH) field, and therefore the absolute value of the pressure, p, is known. This enables the calculation of the geostrophic velocity to the sea floor. The geostrophic transport, U, is then calculated as the total depth-integrated eastward components of velocity (to exclude westward flows from recirculations, slope water current and deeper flows) from 30°N to the coast, as follows: where H is the height of the water column and s and n are the south and north latitudinal limits respectively, which here, correspond to 30°N and 42°N. By accounting for transport across a wide range of latitude, we account for substantial lateral movements in the GS core, although this will include the eastward component of eddies and recirculation. As long as these relative contributions to eastward transport do not substantially vary through time, this method should capture GS variability with only limited bias. A similar method was employed by Richardson (1985) who used eastward-only velocities between the two westward countercurrents and estimated the mean GS transport as 93 Sv, which is comparable to that found here. As discussed previously, we examine the influence of the winter Q net on the density fields and baroclinic transport of the GS. Using the thermal wind relation, we relate horizontal density gradients to the vertical shear of velocity, and anomalies thereof-that is, a strengthening of the density gradients across the GS leads to an increase in the current shear. Integrating upwards, this translates into an intensification of the baroclinic transport.
In order to separate transport changes that are due to changes in the SSH field (barotropic and baroclinic) and in the density field (just baroclinic), the pressure gradient associated with each can be examined. Some change in the SSH will be related to the expansion/contraction of the water column due to changes in density, that is, the steric height. If the change in steric height is exactly equal to the SSH then this translates to a transport change that is, entirely driven by the density gradient. In general, part of the SSH change will not be explained by the change in steric height, for example, it could be associated with wind-driven mass convergence. In order to evaluate this, the surface steric height is calculated using 4,000 m as a reference depth of no motion. 4,000 m is used as it is the deepest depth that remains clear of bottom topography along the section until about 38°N. This assumes that density changes below this depth create negligible transport changes on the timescales of interest here. Specifically, the variation in the missing transport from 4,000 m is found to be much smaller (+/5 Sv) than the overall transport changes (+/40 60 Sv) in the model.
Throughout the study, Monte Carlo analysis is undertaken to quantify the statistical significance of anomalous results by calculating 10,000 random composite anomalies, that is, finding the anomaly between two sets of average variables, for example, Q net, MLD, over 5 random years (from 1978 to 2010). The significance of the actual composite anomaly is then determined by defining the 90% and 95% confidence intervals.

Lagrangian Particle Tracking
In order to investigate variability in GS pathways, we use the particle tracking software Ariane, an offline mass-preserving Lagrangian scheme (Blanke & Raynaud, 1997). This is used with the 5-day mean velocity fields of ORCA12, using multiple point particles (as e.g., Jacobs et al., 2019;Popova et al., 2013). The particles are advected in the model velocity field where they characterize and quantify selected features of ocean circulation from a Lagrangian perspective.
For each experiment, 400 particles are uniformly distributed within a single 1/12° grid box (i.e., 20 × 20) at a particular depth and released into the GS upstream of separation from Cape Hatteras. This is close to the maximum number that can be released successfully within a single grid box, which ensures all particles are immediately entrained into the GS. The specific release point is at 26°N 80°W ( Figure  with possible flows entering the Azores Current or the northern recirculation (Schmitz, 1996). A total of 31 experiments are conducted, one for each year from 1980 to 2010 (starting from 1980 to avoid the potential for spurious results related to model spin-up).

Geostrophic Transport Variability
Interannual variability of GS geostrophic transport in April at 70°W is evident in Figure 3a. Excluding 1978, the range in April transport exceeds 70 Sv with a maximum of 175 Sv in 2007 and a minimum of 105 Sv in 1993. This variability is of the same order as observed by L. V. Worthington (1977), who found a range of ∼40 Sv in the early-mid 1970s. The large absolute transport values in Figure 3a are due to selective use of eastward velocities-see Equation 1. Much of the eastward transport is in fact compensated by westward transport in energetic mesoscale eddies. However, the eastward-only velocities were found to be most effective in capturing GS transport variability due to the existence of strong westward counter-currents. Despite the existence of mesoscale variability, the transport at 70°W is significantly correlated (and thus representative) of similar sections from 74 to 60°W.
The five years of greatest (weakest) transports are marked in red (blue) in Figure 3 and will henceforth be referred to as strong (weak) years, see Table 1. Figure 3b shows the increase or decrease in geostrophic transport between consecutive Aprils, that is, April ( year-to-year changes of the geostrophic transport (measured in [a]) referred to as delta transport that is, ΔT = T i −T i−1 , where T is the transport. Here the red markers represent the five highest increases in transport from the previous year and the blue markers represent the five highest decreases in transport from the previous year. greatest increases/decreases in transport are +40 Sv over 1984+40 Sv over -1985 and −60 Sv (∼38%) over 2007-2008. As in Figure 3a, red (blue) markers represent years of largest increase (decrease) in transport from the previous year, which are subsequently referred to as strengthened (weakened) years, see Table 1.

Strong years Weak years Strengthened years Weakened years
As outlined in Section 2.2, two components are required to calculate GS transport; the SSH and the three-dimensional pressure gradient. The SSH slopes down from south to north, inducing a large eastward GS transport. This is opposed by the pressure gradient that is, controlled by isopycnals sloping upwards from south to north, which acts to reduce the large eastward transport relative to the surface flow. These two components are interrelated as changes in the water column's pressure, and therefore, density gradients will induce changes in its steric height. This will lead to changes in both the baroclinic and barotropic component, that is, both the steric height and SSH.
In order to establish how much GS geostrophic transport variability (in Figure 3a) is due to density changes, that is,. changes in the three-dimensional pressure field (baroclinic component), the steric height of the water column is calculated and is used as a proxy for SSH. We call this the steric sea level. The model SSH includes changes in both the barotropic component, that is, due to convergence from changes in the wind stress and other minor processes, and changes in the baroclinic component, that is, thermal expansion/ contraction due to changes in the temperature of the water column. The steric sea level is used to calculate the geostrophic transport, which is compared to the original geostrophic transport, that is, using the model SSH. If the new transport (using the steric sea level) is similar to the original transport (using model SSH) this proves that the majority of the transport change from year-to-year is due to density changes, that is, changes to the steric height of the water column.
The geostrophic GS transport is calculated at each latitude from 30 to 37.5°N along 70°W, using the proxy and real transports, in April from 1978 to 2010. For the purpose of this figure the latitudinal limit excludes shallower waters further north, where we expect GS transport to be lower. The squared correlation coefficient (r 2 ) between the two transports as a function of latitude, from 30 to 37.5°N, is shown in Figure 4a. It reveals that r 2 > 0.5 at most latitudes. This confirms that the SSH, and therefore the geostrophic transport, along 70°W is mostly explained by the steric changes associated with density with only a minor role being played by processes such as convergence/divergence of mass associated with wind stress variability that induce changes in the barotropic flow.
Additionally, the total change (from the previous year) in eastward transport (as in Figure 3b) using the steric sea level is calculated along 70°W (from 30 to 37.5°N) from 1978 to 2010. This is shown as a percentage of the original method (using model SSH) in Figure  × 100. The original transport time series for both methods is shown in Figure S1. Taking the average across the timeseries, it reveals that more than 88% of the transport changes are due to changes in the density gradient, that is, most of the SSH changes are explained by steric changes associated with density, and despite the large range, no strengthened or weakened year has a transport difference of less than 60%, further highlighting the importance of the density field during years of large changes in GS transport. Additionally, for the westward transport changes (not shown), the mean for the timeseries is more than 86%. This confirms that the majority of the interannual variability in GS transport at 70°W, in April, is due to density changes in the water column. These density changes could be due to strong winter heat losses, which occur during severe winters (L. V. Worthington, 1977), and are investigated in relation to geostrophic transport changes in the next section. Figure 5a shows the mean winter Q net over the wider GS region, which reveals a band of high surface heat loss, exceeding 400 Wm −2 , over the GS path. Figure 5b is a composite of the winter Q net for the strong years minus the weak years (defined in Section 3.1). A negative anomaly of up to −50 Wm −2 is seen to the south of the GS from ∼72 to 55°W, which indicates that a larger heat loss has occurred in this region in the winters preceding the years of strong GS transport compared to the years of, relatively, weak GS transport. It should be noted that if an anomalous Q net over the course of a given year significantly changes the thermal structure (and thus transport) across the GS, then the Q net will change the thermal structure from that existing at the beginning of the year (rather than some mean state). Consequently, as well as examining composites of strong versus weak transport years we also examine composites of the strengthened minus the weakened years. This is consistent with the analysis of (L. V. Worthington, 1977) who noted that the severe winter of 1976/1977 led to an increase in transport (50%) relative to an earlier winter.

Surface Heat Flux and Mixed Layer Depth
The Q net composite for strengthened minus weakened years ( Figure 5c) reveals a larger, that is, more negative, Q net anomaly that covers the entire western subtropical gyre region. Much of this Q net anomaly is statistically significant (90%) with the greatest heat loss, that is, in the years of greatest transport increase,  Figure 5d shows the mean April MLD, which reveals deeper mixed layers exist to the south of the GS from 65 to 45°W, which is consistent with climatological observations (de Boyer Montegut, 2004). The April MLD composite for the strong minus weak years, shown in Figure 5e, reveals patches of positive anomalies (deepening) of up to 100 m in the Sargasso Sea from ∼75 to 40°W. However, there is a lack of any coherent anomaly pattern. In contrast, a statistically significant positive anomaly of more than 200 m, corresponding to mixed layer deepening, can be seen in Figure 5f, which is a MLD composite of the strengthened minus weakened years. This deepening is centered on a region further east, that is, from 30 to 40°N, 64 to 40°W, which is the main site for mode water formation and is, therefore, susceptible to deep mixed layers as this is where the 18°C isotherm outcrops at the surface (black line in Figure 5d) during late winter/early spring (Joyce et al., 2000). The deeper mixing that is, taking place in this region in years of strengthened GS transport leads to MLDs more than 200 m deeper than in years of weakened GS transport.
Following Worthington's anticyclogenesis mechanism, deepening of the mixed layer during strengthened years is expected near 70°W. However, in strengthened years mixed layer deepening occurs farther east, as discussed above, while a weaker shoaling (∼50 m) occurs farther west. This agrees with Qiu and Huang (1995), who found that winter mixed layers were generally shallower than 200 m to the west of 60°W. A possible explanation of the MLD dipole is that there is an increased upper ocean heat content through heat advection by the GS, leading to early re-stratification. This could cause deeper mixed layers to be capped farther to the west (Dong & Kelly, 2004;S. Dong et al., 2007) and implies that the oceanic impact of severe winters can be damped by an elevated ocean heat content. The correlation between heat loss and MLD properties may thus break down; which is commonly found to occur farther west (S. Dong et al., 2007;Talley & Raymer, 1982;Warren, 1972;Yasuda & Hanawa, 1997). In order to examine this further, selected temperature profiles will be examined in the next section.

Temperature Variability
Stronger, more significant relationships have been found between MLD, Q net , and transport in strengthened/weakened years (as opposed to strong/weak years), so we hereafter examine only these years. This is due to expectations that the Q net will change the thermal structure from that at the previous end-of-winter state to the current end-of-winter state, as mentioned in the previous section. Figure 6a is a composite of the meridional temperature profile along 70°W in April for the strengthened years minus the weakened years with the corresponding mean field revealed in Figure 6c. In order to avoid spurious temperature anomalies associated with shifts in the GS path and associated mesoscale variability, the composite is calculated relative to the GS core in each month, which is defined as the latitude of the maximum surface velocity at 70°W.
In the difference between strengthened and weakened years (Figure 6a), a cold anomaly with values of up to −4°C is revealed to the north of the GS core in the Slope Water and a smaller, but still significant, negative anomaly of up to −2.5°C is seen to the south of the GS core at depths greater than 300 m. A warm anomaly of up to 1°C is also found in the upper (∼200 m) region of the Sargasso Sea between 0.4° and 3.0° south of the GS core. This implies that when the GS transport is intensified relative to the previous year, there is a large heat loss to the north (cold anomaly), which suggests that the warmer isotherms outcrop farther to the south in these years (see the mean field in Figure 6c). In addition to this, there is evidence that the shallow MLD seen to the south of the GS at 70°W (Figure 5f) may have been caused by warmer surface waters (above 200 m), preventing the GS transport from increasing by strengthened temperature gradients from the south. However, the cold anomaly seen beneath 200 m could suggest a deeper mixed layer existed in previous months but by April was capped by the warmer waters. This capping seen to the south could be due to heat advection processes as suggested by S. Dong et al. (2007). Further investigation of the evolution of the vertical temperature profiles at 34°N, 70°W between September and April (averaged over strengthened years) supports this (Figure 7a). The deepest mixed layer can be seen in February, with shoaling and re-stratification beginning in March, when warm advection is likely causing capped mixed layers to occur farther west. The MLD shoals slightly later (April) at 54°W and the advection of warmer surface water is absent (Figure 7b). Additionally, the advection of warmer waters at 70°W is not as pronounced during the weakened years (Figure 7c).
The changes in the temperature structure could be either associated with the anomalous Q net (Figure 5b Table 2). South of the GS, strengthened years on average are preceded by net cooling from surface fluxes and warming from OHT convergence. This implies that it is the surface heat fluxes rather than the OHT convergence that contribute to the low temperature to the south of the GS. However, we note one strengthened year in particular, 1999, where anomalous OHT convergence acted strongly to cool the region. In spite of this, on average, the OHT acted to warm the southern box during the strengthened years, which supports the argument that these years are sometimes characterized by warm advection that caps the mixed layer. Immediately north of the GS, anomalous OHT convergence can either act to warm or cool the top 400 m in strengthened years. However, during the strengthened years the average effect of OHT convergence is to warm this region while the average effect of the anomalous surface fluxes is to cool it. In contrast, both boxes are associated with net warming from the surface heat flux and from OHT convergence during the weakened years. Figure 6. Composited difference (strengthened minus weakened years) in (a) April temperature (°C) and (b) April meridional temperature gradient (°Cm −1 ), both calculated relative to the Gulf Stream core (here defined as the location of the maximum surface velocity) along 70°W of the ORCA12 run. Black contours denote the 90%-95% confidence intervals and straight dashed line represents the location of the core of the Gulf Stream. (c) Mean April temperature (°C), lines represent the 18°C isotherm for the total mean (black), mean strengthened years (white) and mean weakened years (magenta) and (d) mean April meridional temperature gradient (°Cm −1 ), also calculated relative to the Gulf Stream core, along 70°W in the ORCA12 hindcast for the period 1980-2010.
The composite of the differences in meridional temperature gradients along 70°W (also relative to the GS core) in April, for strengthened years minus the weakened years, is shown in Figure 6b, with the corresponding mean field displayed in Figure 6d. A positive difference is associated with increasing temperature to the north and vice versa. Figure 6d reveals that negative temperature gradients exist on the northern flank of the GS above ∼400 m and on the southern flank of the GS from ∼400 to 900 m. In Figure 6b, a negative anomaly is apparent on the northern flank of the GS to depths of ∼600 m, while a warm anomaly is seen beneath this, with the greatest difference occurring in the top 200 m. This suggests that during years of increased transport, from the previous April, there is a steepening of the meridional temperature gradients across the GS core but predominantly on the northern flank due to more southward outcropping of warm isotherms (Figure 6c). This led to an intensified GS in the upper 200 m at 70°W ( Figure S3).
This agrees with findings of Zheng et al. (1984), who observed that heat loss in the Slope Water causes an uplift of isotherms that outcrop further south ( Figure 6c) and cause an increase in the baroclinic transport. A similar situation was found by S. Dong et al. (2007), identifying that years with a large transport are associated with strong thermal wind shear across the current. They also suggested that a warming on the southern flank of the GS could correspond to a strengthened upper southern recirculation. This will be investigated in the following section and may provide a link between the deeper mixed layers found farther east and the strengthened GS seen at 70°W.
In contrast, deepening of the mixed layer during strengthened years (observed in Figure 5f) suggests a Worthington style response in the horizontal temperature gradient farther to the east, near 54°W. At this longitude, the GS no longer has a narrow core and there is greater occurrence of mesoscale activity (Krauss et al., 1990), which makes it difficult to account for latitudinal path shift when examining the temperature and temperature gradients.  visible at 54°W from about 200 to 500 m during the strengthened years compared to the weakened years ( Figure S2). This caused greater meridional temperature gradients in the GS core (0-600 m) during the strengthened years from 37 to 38°W. The strengthened temperature gradients at depth at 54°W is consistent with a Worthington style response (L. V. Worthington, 1977), which contrasts with a Zheng style (Zheng et al., 1984) response occurring at 70°W.

Lagrangian Analysis
In order to examine the variability of the SRG during strengthened and weakened years, Lagrangian analysis is undertaken using the particle tracking software Ariane (see Section 2.3). Figure 8c is the composite of the particle density (i.e., the number of particles that have passed through each 1° grid box) for the strengthened (Figure 8a) minus the weakened years ( Figure 8b) after they have traveled in the GS system for 3 months. The differences reveal that, during the strengthened years, more particles are found to the south of the main GS core between 62°W and 75°W with fewer particles found farther to the east. This suggests that during years of intensified GS transport there is a stronger westward flow on the southern flank of the GS, which indicates a stronger southern recirculation. This is consistent with the greater transport seen at 70°W and the weaker flow continuing eastward into the North Atlantic Current. The greater number of particles seen to the west of 62ºW suggests that there is a tighter southern recirculation cell with an intensified westward component during the strengthened years, which enables a stronger return flow and increases the transport at 70°W. Additionally, this potentially provides the warm advection necessary for the capping of mixed layers seen farther west (Figure 7a).
In order to quantify the potential changes in the strength of the southern recirculation, three metrics have been calculated based on the location of the particles after 3 months and the direction of flow: 1. the number of particles residing in the eastern SRG (Box 1, located from 25 to 40°N and 65 to 35°W) and in the western SRG (Box 2, located from 25°N to 40°N and 65°W-coast); 2. the number of particles residing in the total SRG, which is defined as the sum of Box 1 and Box 2; 3. the proportion of particles that flow west during the third month when residing in either the whole SRG (Box 1) or the western SRG (Box 2) is quantified. The percentage of particles that fall into Box 1 and Box 2 for each year are averaged over all years, and for the strengthened and weakened years, with the difference taken between the two-see Table 3. Nearly all of the particles are located in the SRG and more than half are found in the western SRG after 3 months during the strengthened years. In contrast to this, during the weakened years, far fewer particles released reside in the SRG or western SRG after 3 months. Specifically, during the strengthened years, 32.1% more particles are found in the SRG and 47.7% more particles are found in the western SRG than during the weakened years. Both of these anomalies are statistically significant to the 95% confidence level, which was calculated using Monte Carlo analysis, described in Section 2.2. A graphical representation of this is shown in Figure 9a, which illustrates the average location of particles after bifurcation of the GS for all 31 years, the strengthened years, and the weakened years. Of the particles not residing in the SRG, the remainder (Rem) is comprised of the North Atlantic Current, northern recirculation and those still located in Note. Errors represent 1 standard deviation with the SRG and western SRG anomalies significant to more than 2 standard deviations. Abbreviation: SRG, southern recirculation gyre.

Table 3 Percentage of Particles Found in the SRG (Box 1 and Box 2), in the Western SRG (Box 2) and the Eastern SRG (Box 1) After 3 Months, Averaged for the Strengthened Years and Weakened Years With the Difference Taken
Between the Two. Errors in italics Figure 9. The number of particles residing in the western SRG (southern recirculation gyre) eastern SRG, and the remainder (Rem) (which is comprised of the North Atlantic Current, northern recirculation and the remainder left in the main current) after 3 months (a) and the number of particles traveling west during the third month for particles residing in the SRG and western SRG (b). Both averaged over all years, strengthened years, and weakened years.
the main GS core. This further shows that more particles reside in the SRG and fewer in the remainder, after three months in the strengthened years, and the opposite in the weakened years. Figure 9b shows the number of westward-moving particles in the SRG in the third month after release. The sharp increase in westward-moving particles in the western SRG during strengthened years confirms that much of the increase in particle numbers found in the western SRG ( Figure 8) is associated with westward flow and thus with strengthened surface recirculation. This stronger recirculation augments transport at 70°W, contributing to the increase seen there. Additionally, a tighter recirculation implies a strengthening of the anticyclone, which is consistent with a Worthington style response. SRG, southern recirculation gyre.

Regional Wind Stress Curl
In addition to Q net , the wind stress curl (WSC) is also known to influence the transport of WBCs via Sverdrup balance with a greater negative WSC inducing an increased southward interior Sverdrup transport and an increased poleward return flow (GS) at the western boundary (N. G. Hogg & Johns, 1995;Sverdrup, 1947) according to: where β is the variation of the Coriolis parameter f with latitude, V is the meridional transport, ρ is the density of seawater and τ s is the wind stress. To focus on regional WSC changes we integrate westwards from 45°W e  to get anomalous (strengthened minus weakened) Sverdrup transport, ΔV: where  is the longitude, θ is the latitude and x is the distance. Figure Figure 11. JFM, January-March strengthened years, the anomalous WSC caused changes in the Sverdrup transport of up to −1.5 Sv over the subtropical gyre, which translates to a northward GS transport of +1.5 Sv. This confirms that changes in Q net are the main driver for GS transport changes during the strengthened years (increases of ∼40 Sv-see Figure 3b). Greater negative WSC anomalies may also be driving the tight western recirculation seen during the strengthened years (A. H. Chaudhuri et al., 2009).
To investigate if the two forcing anomaly fields are connected (Q net and regional WSC), we compare the Sverdrup transport (integrated WSC at the western boundary from Figure 10) from 25 to 35°N with the mean Q net , averaged over a box to the south of the GS (see Figure 10), for each year in the hindcast. This is shown in Figure 11 which highlights the strengthened and weakened years. No significant correlation between the two forcing indices exists, but two standout strengthened years are noticeable, which correspond to 1985 and 2007, and exhibit considerably negative Q net and Sverdrup transport values. Although one of the other strengthened years (2009) has been forced with near-average winter Q net and WSC values, the remaining two strengthened years appear to occur immediately after winters with large heat losses that do not coincide with a substantial response in WSC. This indicates that WSC occasionally varies in conjunction with Q net , reinforcing changes in GS transport, while Q net can act without a significant WSC response to drive changes in GS transport.

Meteorological Characteristics
The mean winter JFM sea level pressure, 2 m air temperature and 2 m wind speed patterns are investigated in order to ascertain if a common set of meteorological conditions prevails during the strengthened or weakened years, which may be strongly influencing Q net . Figure 12 shows these anomalies for the strengthened years minus the weakened years. It reveals a cooler air temperature of up to 2°C over the western North Atlantic, which is associated with the flow of polar, continental air shown by both the negative sea level pressure and southward wind anomalies.
This suggests that, during the strengthened years, more frequent and/or more intense low-pressure systems occur over this region. These cyclonic systems, known as cold air outbreaks, initiate a northerly wind that brings cold, dry air south over the GS and, consequently, causes a large Q net due to the higher temperature difference and stronger wind at the air-sea interface. The contours of anomalous Q net and MLD, from Figures 5b and 5d, are shown on this plot, which reveal that the location of greatest heat loss and MLD are located within this region of cold, dry air.
This further highlights the dominant influence of Q net over WSC in driving GS transport changes, as a stronger subtropical high would also be required to initiate greater WSC forcing south of 35°N. Although in some years in which the GS strengthens (for example, 2007, not shown), such a pattern (i.e., a stronger subtropical high) can occur, supporting the case for both roles of Q net and WSC, Figure 12 indicates that this is not typically the case.

Conclusions
We have presented results on the role of intense air-sea heat fluxes in driving fluctuations of GS transport using the high resolution ORCA12 hindcast. Significant relationships have been identified between the transport at 70°W and Q net and MLD in the western subtropical gyre, and the temperature and meridional temperature gradient with depth along 70°W. Stronger relationships are found with the change in transport compared to the previous April, as opposed to the absolute transport. Greater increases in spring transport are associated with (i) high winter heat losses throughout the western subtropical gyre, (ii) deeper early Figure 12. Mean winter JFM 2 m air temperature (°C) composite difference for the strengthened minus the weakened years with contours (thin black lines) of the winter JFM sea level pressure (mb) composite difference for the strengthened minus the weakened years. Arrows represent the vectorized winter JFM wind speed (ms −1 ) anomaly for the strengthened minus the weakened years. The green and gray contours denote the 50 Wm −2 and 100 m anomalies for the Q net and mixed layer depths composites respectively (i.e., from Figure 5). Thick black line represents the mean Gulf Stream position in April (denoted by the 0.4 ms −1 contour). JFM, January-March spring mixed layers to the south of the main GS core and east of 60°W, and (iii) stronger meridional temperature gradients to the north of the GS at 70°W.
The large heat losses in the Slope Water, causing temperature decreases of up to 4°C, that led to the strengthening of temperature gradients to the north ( Figure 6) agrees with the findings of Zheng et al. (1984). This is in contrast to that of L. V. Worthington (1977), who found an intensification of the GS occurred in conjunction with strengthening temperature gradients to the south from a deeper thermocline after the severe winter of 1976/1977. The formation of deeper mixed layers further to the west was found here to be prevented by surface warming, potentially from increased heat advection in the intensified GS. However, the strengthening of temperature gradients in the main GS core at greater depths occurs alongside deeper mixed layers farther to the east at 54°W, which is consistent with a Worthington style response.
Lagrangian analysis indicated that an increased, westward-intensified, southern recirculation tends to occur during strengthened years, which is expected to contribute to the increased transport seen at 70°W. It may also be facilitating warm advection in this region, leading to the capping of deeper mixed layers farther west, referred to above. The deeper mixed layers found farther to the east (from ∼60 to 50°W), than noted by L. V. Worthington (1977), have the potential to contribute to the increased transport seen at 70°W via this intensified recirculation.
The influence of WSC in setting GS transport fluctuations at 70°W was also investigated. Although WSC can act in conjunction with Q net to produce an intensification of the GS at this location (i.e., in 1985 and 2007), Q net has also been found to dominate changes in GS transport (i.e., in 1987 and 1999). Additionally, integrating WSC from the eastern boundary yields an anomalous Sverdrup transport of only around 2 Sv during the strengthened years. This, along with the prevailing winter meteorological conditions during the strengthened and weakened years, highlights the dominant role of buoyancy forcing.
As mentioned in Section 1, heat fluxes over the GS region may have implications for the North Atlantic storm track and winter conditions over western Europe (e.g., Joyce et al., 2009;Kwon et al., 2010). This is beyond the scope of this study but it is of interest to determine how Sea Surface Temperature anomaly patterns, and therefore widespread variable heat loss, develop in strengthened and weakened years. Specifically, this would have a substantial impact on the baroclinicity of the atmosphere and, consequently, the North Atlantic storm track. Then given the sea surface temperature tripole's influence on the North Atlantic Oscillation (e.g., Buchan et al., 2014;Grist et al., 2019;Visbeck et al., 2001), this could implicate the GS in large-scale climate variability.
To conclude, we have found that air-sea heat exchange plays a dominant role in driving interannual variations in GS transport at 70°W in the ORCA12 hindcast of the NEMO model. Early spring transport increases of up to 35% (compared to the previous year) were found to occur immediately after severe winters were experienced in North America. These increases are the result of strengthened meridional temperature gradients on the northern flank of the GS and/or an intensified southern recirculation.

Data Availability Statement
The model data used for this publication are available from http://gws-access.ceda.ac.uk/public/nemo.