Robust estimates for the decadal evolution of Agulhas leakage from the 1960s to the 2010s

Agulhas leakage, the transport of warm and salty waters from the Indian Ocean into the South Atlantic, has been suggested to increase under anthropogenic climate change, due to strengthening Southern Hemisphere westerly winds. The resulting enhanced salt transport into the South Atlantic may counteract the projected weakening of the Atlantic overturning circulation through warming and ice melting. Here we combine existing and new observation-and model-based Agulhas leakage estimates to robustly quantify its decadal evolution since the 1960s. We ﬁ nd that Agulhas leakage very likely increased between the mid-1960s and mid-1980s, in agreement with strengthening winds. Our models further suggest that increased leakage was related to enhanced transport outside eddies and coincided with strengthened Atlantic overturning circulation. Yet, it appears unlikely that Agulhas leakage substantially increased since the 1990s, despite continuously strengthening winds. Our results stress the need to better understand decadal leakage variability to detect and predict anthropogenic trends.

S outh of Africa, at the retroflection of the Agulhas Current, relatively warm and salty water leaks from in the Indian Ocean into the South Atlantic, in the form of large anticyclonic eddies (Agulhas rings), smaller cyclonic eddies, and direct inflow [1][2][3][4] . This so-called Agulhas leakage contributes to connect the subtropical gyres of the Southern Hemisphere to one large supergyre 5 and additionally constitutes a choke point for the surface branch of the global overturning circulation 6 .
Changes in Agulhas leakage can impact regional to global climate variability through various processes on interannual to millennial timescales 7 . In particular, previous studies suggest that during anthropogenic climate change Agulhas leakage has been increasing [8][9][10][11] as a result of strengthening Southern Hemisphere westerly winds, in turn caused by increasing anthropogenic greenhouse-gas concentrations combined with ozone depletion [12][13][14] . The resulting enhanced transport of salt into the South Atlantic may counteract the projected weakening of the Atlantic Meridional Overturning Circulation (AMOC) due to warming and ice sheet melting 7,12,15,16 . However, its turbulent and intermittent nature makes Agulhas leakage difficult to observe and to simulate within ocean and climate models. Consequently, estimates for its past and future evolution are sparse and individual estimates are associated with considerable uncertainties.
Here, we jointly analyze a set of already established as well as new observation-and model-based estimates for Agulhas leakage variability to robustly quantify it's (sub-)decadal evolution since the 1960s. We start with a review and update of the three major approaches to estimate Agulhas leakage variability, which are: (1) Lagrangian Agulhas leakage estimation based on the analysis of virtual fluid particle trajectories simulated with ocean models 17 (AL LA ) (2) Agulhas leakage reconstruction based on its imprint on sea surface temperature 10 (AL SST ), and (3) Agulhas leakage reconstruction based on its link to sea surface height (SSH) variability 18 (AL SSH ). We further introduce a new Lagrangian approach for the reconstruction of Agulhas leakage variability from SSH data (AL LA-SSH ), which is particularly suited for the comparison with model-based AL LA estimates. We then jointly analyze AL LA estimates from two hindcast simulations with the eddy-rich ocean model INALT20 11 (Fig. 1): one under the longestablished CORE forcing covering the period 1958-2009 11 (SIM CORE ) and another one under the newer JRA55-do forcing covering the period 1958-2018 17 (SIM JRA ). These are compared to AL LA estimates inferred from the reanalysis product BRAN2020 19 (REA BRAN ), AL SST reconstructions based on the Hadley Center Sea Ice and Sea Surface Temperature data set 20 (HadISST), as well as AL SSH and AL LA-SSH reconstructions based on altimeter products.
We find that simulated and observation-based estimates for Agulhas leakage variability feature a robust increase between the mid-1960s and mid-1980s but yield ambiguous results for the temporal evolution since the 1990s. The increase in Agulhas leakage between the mid-1960s and 1990s is at least partially due to an increased interocean exchange outside eddies via a shorter and more direct path including the so-called Good Hope Jet. It can be linked to the strengthening of Southern Hemisphere winds and coincided with an increase in the AMOC at 34°S, which propagated into the North Atlantic within one to two decades. For the 1990s and 2000s, ocean model simulations under JRA55-do forcing as well as observation-based Fig. 1 The eddy-rich ocean model INALT20 is used to estimate Agulhas leakage (AL) via Lagrangian particle tracking (AL LA ), to update an existing method to reconstruct AL from sea surface temperature (AL SST ), and to establish a new Lagrangian method to reconstruct AL from sea surface height (AL LA-SSH ). a 5-day mean snapshot of simulated upper ocean (134 m) current speed. b Release (AC) and sampling (GH) sections for the Lagrangian simulations performed to estimate AL LA , exemplary leakage and no-leakage trajectories are displayed in red and gray lines, respectively. c 5-day mean snapshot of simulated near surface (16 m)  reconstructions reveal no increase in Agulhas leakage, despite continuously strengthening Southern Hemisphere westerlies. Only ocean model simulations under CORE forcing show a continuous or even accelerated increase in Agulhas leakage after the mid-1980s, which appears to be linked to inaccuracies in the atmospheric forcing product.
Our results stress the need to better understand the nature and drivers of decadal variability of Agulhas leakage, to better detect and predict long-term trends and their impact on the AMOC.

Results and discussion
How can Agulhas leakage variability be estimated? Lagrangian Agulhas leakage estimation based on 3D velocity output from ocean model simulations (AL LA , Fig. 1b): This method consists of (i) releasing virtual fluid particles in the Agulhas Current, each representing a fraction of the current volume transport, (ii) advecting the particles with the simulated 3D velocities, and (iii) sampling and summing up the transport of all particles that enter and remain in the South Atlantic. However, the exact Lagrangian experiment design varies among studies (different applied Lagrangian tools, release and sampling strategies, and reference time for transport calculations) and can impact the estimated mean transport and interannual variability of Agulhas leakage. We here follow the method suggested by Schmidt et al. (2021) 17 , which allows for robust estimates of (sub-)decadal leakage variability: we employ the Parcels tool 21,22 , release particles every 5 days over the whole width and depth of the Agulhas Current at 32°S, and estimate annual Agulhas leakage transport as the cumulative transport of all particles that are released in the same year and are crossing the approximated Good Hope section at least once within 4 years after their release (see Fig. 1b for the exact definition of the release and sampling sections as well as exemplary simulated particle trajectories). The resulting annual leakage transport timeseries for SIM JRA , SIM CORE , and REA BRAN are displayed by the thin black, blue, and cyan lines, respectively, in Fig. 2a.
Agulhas leakage reconstructions based on sea surface temperature observations (AL SST , Fig. 1c): Biastoch et al. 10 introduced a method to reconstruct Agulhas leakage variability from sea surface temperature (SST) variability. They analyzed the relation between Agulhas leakage (AL LA ) and year-to-year SST variability in a set of high-resolution models and found that with increasing leakage, the regional SST difference (ΔSST) between the South Atlantic (SA, average between 10-20°E and 34-40°S) and the generally warmer Indian Ocean (IO, average between 25-35°E and 36-40°S) decreases, since the SST in the SA increases and that in the IO decreases (see Fig. 1c for the location of the utilized SA and IO areas). Subsequently, they performed a regression of the simulated AL LA onto the simulated ΔSST, and applied this regression to the timeseries of observed ΔSST to retrieve an observation-based reconstruction of the Agulhas leakage timeseries (AL SST ). We repeated this analysis in SIM JRA and SIM CORE (Supplementary Fig. 1). The obtained regressions can explain the major part of the simulated year-to-year leakage variability, as the AL SST obtained from applying the regression to the simulated ΔSST correlate significantly (two-sided students t-test with 95% confidence interval) with the original AL LA , with r values of 0.77 and 0.79 for SIM JRA and SIM CORE , respectively. When applying the regressions from SIM JRA and SIM CORE to the timeseries of ΔSST retrieved from HadISST data, we obtain observation-based AL SST reconstructions (the thin light-red line in Fig. 2b shows annual AL SST transport anomalies based on the regression from SIM JRA ). Notably, the two reconstructions based on the regressions from SIM JRA and SIM CORE are hardly distinguishable ( Supplementary Fig. 4a). This emphasizes that the AL SST reconstruction is much more sensitive to the observed ΔSST difference than to the exact coefficients obtained from the Simulated and reconstructed estimates for Agulhas leakage (AL) reveal that AL most likely increased between the mid-1960s and mid-1980s, but yield ambiguous results for the evolution of AL since the 1990s and the absolute transport values. a Absolute AL transport estimates. b Annual (thin lines), as well as sub-decadally filtered (thick lines, 11 yr Hamming window) AL anomalies referenced to the 1964-2005 mean (or the long-term mean for all timeseries not covering the full period). c Moving centered 21year linear AL trends (filled circles mark values that are significant according to twosided student's t-test and 95% confidence level). AL estimates are (i) inferred via Lagrangian particle tracking using the full 3D velocity fields (AL LA ) simulated in INALT20 under JRA55-do forcing (SIM JRA, SIM JRAo ) and CORE forcing (SIM CORE ), and obtained from the ocean reanalysis BRAN2020 (REA BRAN ); (ii) reconstructed from sea surface height following Le Bars et al. (2014) 18 using the difference of Agulhas Current and Agulhas Return Current strength calculated from geostrophic velocities (AL SSH ) along Topex/Poseidon satellite track T122; (iii) newly obtained from 2D Lagrangian particle tracking using geostrophic velocities from AVISO (AL LA-SSH , details in Supplementary Fig. 2); (iv) reconstructed following Biastoch et al. 10 from sea surface temperature (AL SST ) using HadISST data; and (v) inferred from observed drifter and float trajectories (AL Richardson , AL Daher ; the transparent green shading represents the error estimate provided by Daher et al. 32 ). model regressions, which speaks for the suitability of the method to retrieve a rather independent observation-based estimate for Agulhas leakage variability (the mean transport is, however, determined by the mean leakage transport of the model simulation used to retrieve the regression, so that the AL SST estimate cannot be used to infer the mean leakage strength).
Agulhas leakage reconstructions based on sea surface height observations (AL SSH and new AL LA-SSH , Fig. 1d): Le Bars et al. 18 suggested to reconstruct Agulhas leakage variability from sea surface height (SSH) observations. More precisely, they determined Agulhas leakage as the difference between the Agulhas Current and Agulhas Return Current transports, each calculated from geostrophic velocities inferred from SSH across pre-selected satellite tracks (AL SSH estimate for track 122 is displayed by yellow line in Fig. 2a). We here propose a complementary Lagrangian method for the reconstruction of year-to-year Agulhas leakage variability from the full 2D geostrophic velocity field inferred from SSH (AL LA-SSH ). It again first relies on establishing a relation between Agulhas leakage and SSH in an ocean model ( Supplementary Fig. 2). To do so, the Lagrangian leakage estimation as explained above is performed once with the full 3D field of simulated velocities (AL LA ), and another time with the 2D geostrophic velocities inferred from simulated SSH (AL LA,geo ). Subsequently, AL LA is regressed onto AL LA,geo . As for AL SST , the resulting AL LA-SSH reconstruction can explain the major part of the simulated year-to-year leakage variability, as the timeseries obtained from applying the regression to AL LA,geo derived from the model simulations themselves correlate significantly with the original AL LA , with r values of 0.86 and 0.79 for SIM JRA and SIM CORE , respectively. This implies that the dominant part of year-to-year Agulhas leakage variability can be linked to variability in the geostrophic components of the oceanic flow field related to SSH variability. Subsequently, the regressions from SIM JRA and SIM CORE are applied to the AL LA,geo timeseries retrieved from a third Lagrangian experiment, for which particles were tracked with the 2D geostrophic velocities inferred from AVISO data, and yield observation-based AL LA-SSH reconstructions (the thin purple line in Fig. 2b shows annual AL LA-SSH transport anomalies based on the regression from SIM JRA ). Again, as for AL SST , the two reconstructions based on regressions from SIM JRA and SIM CORE are hardly distinguishable, supporting the suitability of the method to infer Agulhas leakage variability (while the mean transport cannot be inferred). Hence, our proposed method yields a second largely independent observation-based Agulhas leakage estimate that is optimally suited for comparisons with Lagrangian leakage estimates from ocean model simulations. Additionally, this method is potentially very valuable for analyzing (climate) model output for which the full 3D velocity field is not stored or at least not at a sufficient temporal resolution.
We here focus on these three main methods to estimate Agulhas leakage variability, since they provide the most comprehensive and independent set of estimates. For the analysis of the representation of Agulhas leakage in ocean models additional methods have been employed, such as an Eulerian retroflection index, based on the ratio of the Agulhas Current and Agulhas Return Current transports 23,24 , or passive tracer 25,26 and thermohaline threshold [27][28][29] methods. While these approaches have been useful to address specific research questions, they are often model-specific and hence less flexible in their application. There also exist attempts to link the strength of Agulhas leakage to the location of the Agulhas retroflection, as inferred from SSH data 30 . Yet, such an approach directly assumes that specific retroflection dynamics control the variability of Agulhas leakage; the method for the AL LA-SSH reconstruction described above is less dependent on the exact mechanism controlling leakage variability. Consequently, whenever possible, we recommend to use the main methods of leakage estimation described above.
How did Agulhas leakage vary on (sub-)decadal timescales since 1958? We can now jointly analyze the main model-and observation-based estimates to identify potential robust atmospherically-forced signals in the temporal evolution of Agulhas Leakage since 1958.
First of all, we want to emphasize that in this and the following sections, we focus on (sub-)decadal Agulhas leakage variability, as we do not expect any robustness in the observation-and modelbased estimates for the mean leakage strength and interannual variability. Observation-based estimates for the mean leakage strength are still poorly constrained 18,31,32 and model-based estimates are sensitive to the representation of (sub-)mesoscale processes 33 and the large-scale circulation pattern 14 (see Fig. 2a and discussion in the method section). Moreover, though several previous studies suggested that interannual Agulhas leakage variability can be linked to atmospheric wind forcing, we argue, that interannual Agulhas leakage variability is mainly intrinsic oceanic variability, so that its exact timing can hardly be reproduced by any climate or ocean model simulation without data assimilation. This is strongly supported by the fact that a second hindcast simulation with INALT20 under JRA55-do atmospheric forcing (SIM JRAo ), which differs in its initialization and parameter set-up from SIM JRA , yields a substantially different interannual variability while the (sub-)decadal variability is largely reproduced (Fig. 2a, Supplementary Fig. 3). Additionally, the interannual leakage variability estimated from a given velocity field via the Lagrangian methodology has been shown to be sensitive to the details of the experimental set-up, in contrast to the estimated (sub-)decadal variability that can be diagnosed robustly 17 . To better quantify (sub-)decadal variability, in addition to the annual timeseries, we provide 11yr-filtered timeseries (Hamming filter), as well as moving centered 21 yr trends (the significance of trends was tested using a two-sided student's t-test with 95% confidence level).
Second, we want to highlight that all Agulhas leakage estimates, model-as well as observation-based, show an increase in Agulhas leakage between the mid-1960s and mid-1980s ( Fig. 2b), with significant decadal trends of 0.1 to more than 0.2 Sv per year (Fig. 2c). However, they yield ambiguous results regarding the temporal evolution of Agulhas leakage since the 1990s. On the one hand, SIM CORE (blue lines in all panels of Fig. 2) features an accelerated increase in Agulhas leakage in the 1990s and early 2000s, with even higher significant decadal trends than in the 1960s to 1980s; and other previously published model simulations under CORE forcing do so as well [10][11][12]34 . On the other hand, Agulhas leakage estimates from SIM JRA (black lines in Fig. 2) and REA BRAN (cyan lines in Fig. 2), as well as reconstructions based on observed SST (light-red lines in Fig. 2) and observed SSH (purple and yellow lines in Fig. 2) rather show a leveling and feature no significant positive trends in the 1990s and 2000s. Interestingly, when reexamining one of the earliest studies postulating an increase of Agulhas leakage since the 1980s 8 , we find that also their analysis already indicated that most of the increase occurred before the 1990s. We conclude that it is very likely that Agulhas leakage increased between the mid-1960s and mid-1980s, but unlikely that Agulhas leakage substantially increased in the 1990s and 2000s. Given the very good agreement between the temporal evolution of the observation-based AL SST index and the AL LA from SIM JRA , as well as the lack of a significant trend of the AL SSH and AL LA-SSH index in the 1990s and 2000s, the strong increase in the 1990s in AL LA from SIM CORE and other model simulations under CORE forcing should be considered to be to some degree artificial and related to inaccuracies in the CORE atmospheric forcing product.
Our results fit to the hypothesis that Agulhas leakage has already been increasing due to anthropogenic climate change, but that this increase is still comparatively small and masked by decadal variability 10 . The hypothesis is further supported by evaluating the AL SST reconstruction based on the full period of available HadISST data ( Supplementary Fig. 4), which yields a significant long-term linear trend of 0.02 Sv per year between 1870 and 2020, but much stronger alternating decadal trends between −0.08 Sv per year in the 1940s (long-term minimum) and +0. 17 Sv per year in the 1970s (long-term maximum). Therefore, it is crucial to better understand the nature and drivers of decadal Agulhas leakage variability, to be able to detect, quantify, and finally predict any anthropogenic trends.
Are decadal trends in Agulhas leakage reflected in changing Indo-Atlantic transport via Agulhas Rings? A fundamental question regarding the nature of decadal trends in Agulhas leakage is whether they can be found in all the individual contributions to the Indo-Atlantic interocean exchange, consisting of transport via Agulhas Rings, which are large anticyclonic eddies, cyclonic eddies, and direct inflow. Early studies on Agulhas leakage relied on the assumption that Agulhas leakage mainly occurs via Agulhas Rings 4 , implying that Agulhas leakage variability stems predominantly from changes in their number and/or size and hence can be related to changes in the eddy kinetic energy in the Agulhas retroflection region 8,35 . Yet, there also have been conflicting reports of decoupled Agulhas leakage and variability in eddy kinetic energy 14,25 . Moreover, recent studies find a non-negligible or even dominant background leakage component associated with the mean flow, which appears to largely determine the variability in Agulhas leakage 29,34,36 . We here make use of our model simulations SIM JRA and SIM CORE to analyze the potential relative role played by eddies vs direct inflow contributions for (sub-)decadal leakage variability over the past decades.
There exist different approaches to estimate the eddy vs direct inflow contributions of Agulhas leakage. One of the major difficulties for such estimations arises from the fact that the occurrence of Agulhas leakage can hardly be attributed to one distinct instantaneous process, but rather comprises all the processes that contribute to the transport of fluid particles from the Agulhas Retroflection region through the Cape Basin towards the Good Hope section. Hence, we here make use of our simulated Lagrangian trajectories together with a Eulerian eddydetection algorithm to analyze for how long each leakage particle is transported within anticyclonic eddies, cyclonic eddies, or outside eddies during its transit between 20°E (the geographical border between the Indian and Atlantic Ocean) and the approximated Good Hope line (beyond which most of the leakage particles remain in the South Atlantic and do not return into the Indian Ocean).
Averaging over all leakage trajectories and over time (leakage reference years 1958 to 2005 and 1958 to 2014 in SIM CORE and SIM JRA , respectively), we find that leakage particles spend only a bit more than a third of their transit within anticyclonic eddies (35% in SIM JRA , 34% in SIM CORE ) and most of the time outside eddies (59% in SIM JRA and 62% in SIM CORE ). This solidifies the view that a major part of the Indo-Atlantic exchange occurs outside Agulhas Rings. Likewise, decadal leakage variability appears to be mainly driven by transport outside of eddies. In fact, decadal trends in Agulhas leakage and in the average relative time leakage particles spent in anticyclonic eddies are of opposite sign (Fig. 3a, b). An increase in leakage coincides with an increase of the relative time spent outside of eddies, but a reduction of the relative time spent in anticyclonic eddies. Interestingly, this is true for both model simulations over the whole time, i.e., for the diagnosed robust increase in Agulhas leakage between the mid-1960s to mid-1980s, as well as for the potentially overestimated leakage increase in the 1990s in SIM CORE . Hence, at least for the considered time period and within our model simulations, decadal trends in Agulhas leakage cannot simply be attributed toand estimated viachanges in the transport via Agulhas rings. This also suggests that changes in Agulhas leakage do not necessarily covary with changes in eddy kinetic energy, helping to explain previous controversial results in that matter.
An increase in leakage and the reduction in the relative time spent in anticyclonic eddies goes further along with a reduction of the mean transit time between 20°E and the Good Hope line (Fig. 3c). Moreover, maps of particle crossing probabilities indicate that particles follow less often the classical eddy corridor 37 , but more often a northern path close to the African coast, as well as a more southern path (Fig. 3d-f, Supplementary  Fig. 5). Together, these results suggest that the increase in leakage between the mid-1960s and mid-1980s was at least partially the result of an increased interocean exchange via a shorter and more direct path, potentially involving the so-called Good Hope Jet.
Schubert et al. 33 showed that the resolution of submesoscale flows in an ocean model lead to a substantial increase in the mean Agulhas leakage transport, in particular, due to a strengthening of this shorter leakage pathway. The results of the present study suggest that the role of submesoscale flows on the (sub-)decadal leakage variability is however rather small, given the good agreement between the decadal variability of the AL SST estimates based on SIM JRA and observations (where the effect of the submesoscale on the larger scale circulation is included in the observations but not in the simulation). To further validate this hypothesis, high-resolution altimeter products that resolve submesoscale processes, as well as long submesoscale-permitting simulations of the Agulhas current system are necessary, which are unfortunately not yet affordable.
How do decadal trends in Agulhas leakage relate to larger-scale changes in Southern Hemisphere winds and oceanic circulation patterns? Previous studies suggest that the decadal variability as well as anthropogenic long-term trends of Agulhas leakage are mainly driven by changes in the strength of the Southern Hemisphere westerlies [12][13][14]34,38,39 and can be linked to the Southern Annular Mode (SAM) 10,34 .
Also our model simulations SIM JRA and SIM CORE feature a general tendency of a higher Agulhas leakage corresponding to a higher SAM index, here calculated from sea level pressure of the atmospheric forcing fields following Marshall (2003) 40 (Fig. 4a). However, while the detrended sub-decadal anomalies of Agulhas leakage and SAM covary to some degree in both simulations (Fig. 4b vs Fig. 2b; relatively high but not yet significant correlations with r = 0.52 and r = 0.71 for SIM JRA and SIM CORE , respectively), their decadal trends cannot necessarily be related over the whole considered time period (Fig. 4d vs Fig. 2c). The SAM index shows only minor differences between the two reanalysis products used to force SIM JRA and SIM CORE and features a long-term increase since 1958, in agreement with the station-based index (the positive trend in SAM has already been widely discussed 41 ). Until the mid-1980s, the decadal increase in SAM appears representative for the increase in AL for the two simulations as well as for the observation-based estimates. Yet, since the mid-1980s, the decadal evolution between SAM and Agulhas leakage diverges. In particular, the relative increase in Agulhas leakage in SIM CORE appears to be larger than suggested by SAM, and the leveling in Agulhas leakage in SIM JRA as well as in observation-based estimates does not fit to the continuing increase in SAM. Previously, a pause in the increase in Southern Hemisphere westerlies as expressed by SAM from the 2000s onwards was identified 42 , and hypothesized to be the cause for the lack of a significant trend in the AL SSH estimate 43 . However, the observation-based AL SST estimate (and the AL LA estimate from SIM JRA ) strongly suggests that the leveling of Agulhas leakage already started a decade earlier in the 1990s, and hence cannot be solely explained by changes in SAM, which still showed substantial positive tendencies during that period.
The missing significant correlation between (sub-)decadal variability of Agulhas leakage and SAM may be due to the fact that Agulhas leakage is particularly sensitive to regional changes in the strength of westerly wind stress over the southern Indian Ocean, which can deviate from globally-averaged changes represented by SAM 38 . Therefore, we further compared the temporal evolution of Agulhas leakage to the maximum zonalmean westerly wind stress averaged over the Indian Ocean only (20°E-110°E; Supplementary Fig. 6b, e). In some aspects, we indeed find a slightly better agreement than for the comparison of Agulhas leakage and SAM. In particular, for SIM CORE the westerly wind stress increase over the Indian Ocean is larger in the early 1990s than between the mid-1960s and mid-1980s, while for SIM JRA the positive decadal trends remain rather constant. This could help explain the accelerated increase in Agulhas leakage in the 1990s simulated with SIM CORE and suggests that the CORE forcing featuresat least regionallya too strong increase in the westerly wind strength in the 1990s (in addition to the already known overestimated increase in SAM in the 1970s related to biases in the NCEP-NCAR data from which the CORE forcing is derived 44 ). However, there remain other important aspects of the temporal evolution of Agulhas leakage that we still cannot explain. Most importantly, in the 1990s, the relative stable positive decadal trends in regional westerly wind stress still do not align well with the reduced positive or even opposing negative tendencies in Agulhas leakage as simulated with SIM JRA and reconstructed from observations. A more complete explanation can be found in the proposed mechanism for how changes in the strength of westerly winds impact Agulhas leakage, namely, via changes in the Indian ocean wind stress curl (Fig. 4c, e) at the latitudes of the Agulhas retroflection region between 35°S and 45°S 14,34 . An increase in the strength of westerly winds evokes an increase in the regional wind stress curl and corresponding increase in leakage. However, the observed and simulated increase in SAM also goes along with a general poleward shift of the westerlies, which was shown to have the potential to weaken the regional wind stress curl and Agulhas leakage 14 rather than strengthening it 9 . Moreover, the rate of change in the position of the westerlies decreases with time, whereby the decrease is more pronounced in SIM CORE than in SIM JRA (Supplementary Fig. 6c, f). In fact, SIM CORE features a stronger poleward shift than SIM JRA between the mid-1960s and mid-1980s, but a weaker poleward or even equatorward shift in the 1990s. This suggests, that in the 1990s, the change in the position of the westerlies acted to enhance the increase in Agulhas leakage due to changes in the intensity of the westerlies in SIM CORE, but weakened it in SIM JRA . To conclude, the impact of decadal trends in SH winds on Agulhas leakage are best derived directly from the wind stress curl.
It is also important to note that there may be competing drivers of decadal leakage variability related to larger-scale oceanic circulation changes on (sub-)decadal timescales as an additional response to changes in the Southern Hemisphere westerlies (such as the strength of the Antarctic Circumpolar Current), as well as to oceanic circulation changes due to concurrent wind changes at different latitudes (such as the strength of the Agulhas Current) or changes in the thermohaline forcing (such as the AMOC).
Previous studies investigated a potential relation between Agulhas leakage and the strength of the Agulhas Current, yet, with partially conflicting results. In agreement with refs. 14,25,38 but in contrast to ref. 45 , we find no robust relationship between the temporal variability in Agulhas leakage and the strength of the Agulhas Current in SIM JRA and SIM CORE (Fig. 5b, d). While (detrended) sub-decadal anomalies appear to feature some covariability (relatively high, but not statistically significant correlations of r = 0.68 and r = 0.59 for SIM JRA and SIM CORE , respectively), the decadal trends of Agulhas leakage and Agulhas Current transport have the same sign in SIM JRA and opposite signs in SIM CORE . Hence, decadal and longer-term fluctuations in Agulhas leakage mainly represent a changing proportion of the transport of the Agulhas Current flowing into the South Atlantic (Fig. 5a, c, e). This is even true on interannual timescales, as highlighted by very high and significant correlations between the annual timeseries of (detrended) Agulhas leakage to Agulhas Current transport ratio (AL/AC) and AL LA (r = 0.96 and r = 0.97 for SIM JRA and SIM CORE , respectively). Nevertheless, the differences in the simulated temporal evolution of the Agulhas Current may help explain why, from the mid-1960s to mid-1980s, SIM JRA and SIM CORE show a somewhat disproportional response to the increase in regional westerly wind stress: in SIM JRA , an increasing Agulhas Current adds to the impact of increasing regional westerly winds; in SIM CORE , a decreasing Agulhas Current weakens the impact of increasing regional westerly winds.
Several sensitivity studies with ocean models showed that idealized changes in Agulhas leakage can impact the AMOC 12,15,16,39,46 . However, detecting and quantifying the impact of the actual increase in Agulhas leakage since the 1960s has been proven difficult. Direct AMOC observations are still too sparse to investigate decadal changes. But even in realistic hindcast simulations with ocean models the multitude of other factors potentially impacting AMOC variability makes detecting Agulhas leakage driven changes a challenge. So far, the impact of changes in Agulhas leakage has only been inferred indirectly from co-varying changes in the Atlantic meridional heat transport 47,48 and the Atlantic Multidecadal Oscillation 10 , as well as from regional changes in temperature and salinity patterns in agreement with the advective pathways and timescales of Agulhas leakage along the upper limb of the AMOC 49 . Here, we analyze the rate of change in the decadal anomalies of Agulhas leakage and in the latitudinal dependent AMOC strength calculated in density coordinates (AMOC σ ), as simulated with SIM JRA and SIM CORE (Fig. 6). These analyzes show that the simulated decadal increase in Agulhas leakage between the mid-1960s and mid-1980s coincided with an increase in the AMOC σ at 34°S of comparable magnitude, which propagated into the North Atlantic within one to two decades. Notably, the increase reached 26°N in the early to mid-1980s, thereby helping to build up the positive AMOC σ anomaly in the late 1980s and early 1990s, which traditionally has been linked to changes in the remote circulation and deep water formation in the subpolar North Atlantic 50,51 (this impact is also visible in our analyzes as southward propagating positive tendencies originating around 50°N in the late 1970s to late 1980s and arriving at 26°N around 1 decade later). To our knowledge, a coinciding decadal increase in Agulhas leakage and the AMOC strength in the South Atlantic with subsequent northward propagation of AMOC changes has not been detected within hindcast simulations before. We argue that this may partially be related to the fact that it has been common to study the propagation of AMOC changes based on Hovmöller diagrams of decadal anomalies (instead of the rate of change in the decadal anomalies used in this study), which are less suited to detect corresponding changes in the AMOC and Agulhas leakage that only occur over part of the considered time period. We conclude that the South Atlantic may not only be important for future AMOC changes, but may already have modulated basin-wide AMOC variability over the past decades. This underlines the importance of sustained efforts to monitor the AMOC not only in the North Atlantic but also in the South Atlantic, e.g., across the SAMBA array 52 . Moreover, it stresses the need for a more comprehensive assessment of the resolutiondependent representation of Agulhas leakage and its impact on the AMOC in state-of-the-art climate model simulations, as they have the potential to alter regional to global future climate projections.

Methods
Ocean model simulations. INALT20 is a global eddy-permitting (1/4°) ocean/seaice model configuration with an eddy-rich (1/20°) nest in the South Atlantic, the corresponding part of the Southern Ocean and the western Indian Ocean. It is based on the "Nucleus for European Modeling of the Ocean" (NEMO version 3.6 53 ) with "Océan Parallélisé" 54 (OPA) as ocean component and "Louvain-la-Neuve Ice Model" (LIM2 55 ) as sea-ice component. It was developed at GEOMAR as a member of the INALT family 11 , and -together with its predecessorshas been widely used to study the dynamics of the greater Agulhas region and their impact on regional to large-scale ocean circulation and climate 12,14,49,[56][57][58] .
We From the model simulations, the 5-day mean output has been analyzed. The 3D velocity fields have been employed for the Lagrangian AL LA estimation (further details below). For the AL SST reconstruction, annual mean near-surface temperature taken at the vertical model level 3 (~16 m depth) has been used instead of real SST, to allow for a better imprint of ocean dynamics 10 . Sea surface height (SSH) anomalies were used to calculate geostrophic velocities for the AL LA-SSH reconstruction using cdftools, a fortran package for analysis and diagnostics on NEMO ocean model output. The strength of the AMOC was computed at each latitude and for each time step as the maximum of the annual mean AMOC streamfunction in density coordinates: AMOC σ ðt; yÞ ¼ max σ Ψðt; y; σÞ; Ψðt; y; σÞ ¼ Vðt; x; y; σ 0 Þdx dσ 0 ð1Þ where σ is the potential density referenced to the surface, V is the annual mean meridional transport, and x w and x e are the western and eastern boundaries of the Atlantic. 79 evenly spaced density bins between 1023 and 1031 kg m −3 have been used, whereby only densities larger than 1026 kg m −3 were considered for determining the AMOC strength (to avoid accidentally sampling the strength of the subtropical cells).
The representation of the general circulation and mesoscale dynamics in INALT20 has already been thoroughly evaluated and discussed 11,51 . We here want to highlight that the mean spatial pattern of simulated SSH variance, a measure of mesoscale activity, compares very well to that inferred from satellite altimetry data and captures all major regions of elevated eddy activity, including the Mozambique Channel, the Agulhas retroflection region, as well as the Brazil-Malvinas confluence zone. Moreover, the observed horizontal and vertical flow structure of the Agulhas Current along the ACT array 65 is largely reproduced in the model; and the simulated mean transport falls well in the observed range. However, as already noticed before for this ocean model as well as for a climate model at similar horizontal resolution 32 , the long-term mean Agulhas leakage transports are lower than the estimates obtained from surface drifters and floats. 2 Sv, AL/AC: 19.7 ± 2.3%) remains significantly lower than the observed values. As recently suggested 33 and already discussed above, this may be partially due to the missing representation of the submesoscale in the model simulations. In addition, the simulated mean Agulhas leakage transport seems to be impacted by the initialization strategy, the presence of tides, and the representation of the large-scale salinity and temperature fields, as SIM JRAo produces a long-term mean Agulhas leakage transport (12.9 ± 2.4 Sv) that is more than 3 Sv higher than that of SIM JRA . This highlights that an accurate representation of the mean leakage transport in ocean model simulations remains a challenge. Consequently, for this study, we focus not on the mean transport, but on the variability of Agulhas leakage.
Ocean reanalysis. The ocean reanalysis product BRAN2020 (here referred to as REA BRAN ) is a near-global ocean model constrained by ocean observations via data assimilation 19 . The ocean model is the Ocean Forecasting Australian Model, version 3 (OFAM3 66 ) with a horizontal resolution of 1/10°based on the Modular Ocean Model version 5 and forced with the JRA55-do atmospheric dataset 61 . The data assimilation system is a two-step, multiscale optimal interpolation method that uses satellite sea surface temperature and sea level anomaly as well as in situ temperature and salinity from e.g., Argo floats, surface drifting buoys, moorings, and ship-borne surveys. BRAN2020 is available at daily resolution from 1993 to present, but the vertical velocity field is only provided from 1998. The Agulhas Current system is realistically simulated by BRAN2020, but the representation of Agulhas leakage has not been analyzed in detail 67 . It is noteworthy that, in contrast to the non-assimilative ocean models, BRAN2020 features a mean leakage transport that falls within the lower range of the observational estimates (16.9 ± 1.2 Sv).
Lagrangian analyses. The Lagrangian leakage estimation was performed offline by making use of the software Parcels 21,22 and following the recipe published by Schmidt et al. 17 , with the specifications outlined above in the section "How can Agulhas leakage variability be estimated?". For the AL LA estimation, particles were advected with the simulated 5-day mean 3D velocity field, for the new AL LA-SSH estimation with the 2D horizontal geostrophic velocities obtained from 5-day means of simulated and observed SSH anomalies. The positions of the resulting Lagrangian trajectories (lon, lat, depth) were stored daily.
Leakage trajectories were subsequently further analyzed in terms of how long they needed to transit between 20°E (South of Africa) and the Good Hope section, and how much of this transit time they spent within anticyclonic eddies, cyclonic eddies, or outside eddies. The associated eddy detection method is described below. Leakage particle pathways were visualized by means of crossing probabilities obtained from spatial binning of particle positions. The bins span the whole depth range and have a width of 0.5°in latitude and 0.75°in longitude, which is larger than the 95th percentile of horizontal distance traveled per day by all particles. The number of particles weighted by their transport is counted in each bin at each time step and then summed up for each release year. This is normalized by dividing by the temporal average for each release year of the total transport in the whole region with the transport of particles captured repeatedly as long as they remain in the region.
Eddy detection. We detected eddies in the simulated 5-day mean SSH anomaly field, using an established algorithm 68 . A python implementation of this algorithm was obtained from https://github.com/ecjoliver/eddyTracking and further modified to include NEMO model output and improve the required computational time.
The algorithm searches closed contours in spatially high-pass filtered SSH anomaly fields (Gaussian filter to remove wavelength scales larger than 20°). A closed contour is classified as an eddy, if the following criteria are met: (i) the contour contains more than 100 and less than 5000 grid points, and a local maximum (anticyclonic eddy) or minimum (cyclonic eddy); (ii) the eddy amplitude, that is the difference between the local extreme value and the SSH anomaly averaged along the perimeter of the last closed contour of the eddy, exceeds 5 cm; (iii) the eddy size, measured by the maximum distance between any pair of points within the contour, must be smaller than 300 km. The Gaussian filter scale was adopted from the original algorithm. The other parameters (minimum eddy amplitude and maximum eddy size) were chosen in an explorative manner by a visual verification of structures identified/not identified as eddies in a number of snapshots, but generally fit to the characteristics of mesoscale eddies as described in the literature 69 .
A 3D temporarily varying eddy mask was calculated by, for each time-step, tagging all model grid cells within a closed contour classified as a boundary of a cyclonic or anticyclonic eddy with 1 and −1, respectively, and setting all other grid cells to 0.
Calculation of the Southern Annular Mode (SAM) index. We calculated the annual SAM index from sea level pressure (SLP) obtained from observational stations and from the two reanalysis products used to force the model, CORE and JRA55-do, following Marshall (2003) and the guidelines on this website: https:// climatedataguide.ucar.edu/climate-data/marshall-southern-annular-mode-samindex-station-based. The calculation consists of the following steps: (1) selecting SLP values at 40°S and at 65°S, (2) for each latitude calculating zonal mean SLPs, performing annual averaging, subtracting the long-term (1971-2000) mean and dividing by the standard deviation based on the annual anomalies from this mean, (3) subtracting values at 65°S from those at 40°S.
Received: 6 May 2022; Accepted: 24 November 2022; Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.