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

Siren Rühs, Christina Schmidt, René Schubert, Tobias G Schulzki, Franziska

U Schwarzkopf, Dewi Le Bars, Arne Biastoch

► To cite this version:

Siren Rühs, Christina Schmidt, René Schubert, Tobias G Schulzki, Franziska U Schwarzkopf, et al. Robust estimates for the decadal evolution of Agulhas leakage from the 1960s to the 2010s. Communications Earth & Environment, 2022, 3, 10.1038/s43247-022-00643-y. hal-04737755

HAL Id: hal-04737755

https://hal.science/hal-04737755v1

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Siren Rühs ID 1,4✉, Christina Schmidt ID 1,5, René Schubert1,6, Tobias G. Schulzki ID 1,
Franziska U. Schwarzkopf ID 1, Dewi Le Bars ID 2 & Arne Biastoch ID 1,3

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 find 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.

<footnote n="1">¹GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany. ²Royal Netherlands Meteorological Institute, De Bilt, The Netherlands. ³Christian-Albrechts-Universität zu Kiel, Kiel, Germany. ⁴Present address: Institute for Marine and Atmospheric Research, Utrecht University, Utrecht, The Netherlands.</footnote>

<footnote n="5">⁵Present address: Climate Change Research Centre and the Australian Centre for Excellence in Antarctic Science, University of New South Wales, Sydney, NSW, Australia. ⁶Present address: Univ Brest, CNRS, IRD, Ifremer, Laboratoire d'Océanographie Physique et Spatiale (LOPS), IUEM, Brest, France. ❀email: s.ruhs@uu.nl</footnote>
South of Africa, at the retroreflection 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¹⁻⁴. This so-called Agulhas leakage contributes to connect the subtropical gyres of the Southern Hemisphere to one large supergyre⁵ and additionally constitutes a choke point for the surface branch of the global overturning circulation⁶.

Changes in Agulhas leakage can impact regional to global climate variability through various processes on interannual to millennial timescales⁷. In particular, previous studies suggest that during anthropogenic climate change Agulhas leakage has been increasing⁸⁻¹¹ as a result of strengthening Southern Hemisphere westerly winds, in turn caused by increasing anthropogenic greenhouse-gas concentrations combined with ozone depletion¹²⁻¹⁴. 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⁷,¹²,¹⁵,¹⁶. 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 its (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¹⁷ (AL<sub>LA</sub>) (2) Agulhas leakage reconstruction based on its imprint on sea surface temperature¹⁰ (AL<sub>SST</sub>), and (3) Agulhas leakage reconstruction based on its link to sea surface height (SSH) variability¹⁸ (AL<sub>SSH</sub>). We further introduce a new Lagrangian approach for the reconstruction of Agulhas leakage variability from SSH data (AL<sub>LA-SSH</sub>), which is particularly suited for the comparison with model-based AL<sub>LA</sub> estimates. We then jointly analyze AL<sub>LA</sub> estimates from two hindcast simulations with the eddy-rich ocean model INALT20¹¹ (Fig. 1): one under the long-established CORE forcing covering the period 1958–2009¹¹ (SIM<sub>CORE</sub>) and another one under the newer JRA55-do forcing covering the period 1958–2018¹⁷ (SIM<sub>JRA</sub>). These are compared to AL<sub>LA</sub> estimates inferred from the reanalysis product BRAN2020¹⁹ (REA<sub>BRAN</sub>), AL<sub>SST</sub> reconstructions based on the Hadley Center Sea Ice and Sea Surface Temperature data set²⁰ (HadISST), as well as AL<sub>SSH</sub> and AL<sub>LA-SSH</sub> 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

<figure><img src="image_3.png" /><figcaption>**Fig. 1** The eddy-rich ocean model INALT20 is used to estimate Agulhas leakage (AL) via Lagrangian particle tracking (AL<sub>LA</sub>), to update an existing method to reconstruct AL from sea surface temperature (AL<sub>SST</sub>), and to establish a new Lagrangian method to reconstruct AL from sea surface height (AL<sub>LA-SSH</sub>). 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<sub>LA</sub>, 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) temperature, boxes used to calculate the temperature difference between the South Atlantic (SA) and Indian Ocean (IO) used for the AL<sub>SST</sub> estimation are added. d 5-day mean snapshot of sea surface height (SSH) anomaly (global mean subtracted) used to calculate geostrophic velocities for the new Lagrangian AL<sub>LA-SSH</sub> estimation, satellite track employed for the old AL<sub>SSH</sub> estimate by Le Bars et al. (2014) is overlayed.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 2 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 two-sided student's t-test and 95% confidence level). AL estimates are (i) inferred via Lagrangian particle tracking using the full 3D velocity fields (AL<sub>LA</sub>) simulated in INALT20 under JRA55-do forcing (SIM<sub>JRA</sub>, SIM<sub>JRAo</sub>) and CORE forcing (SIM<sub>CORE</sub>), and obtained from the ocean reanalysis BRAN2020 (REA<sub>BRAN</sub>); (ii) reconstructed from sea surface height following Le Bars et al. (2014)<sup>18</sup> using the difference of Agulhas Current and Agulhas Return Current strength calculated from geostrophic velocities (AL<sub>SSH</sub>) along Topex/Poseidon satellite track T122; (iii) newly obtained from 2D Lagrangian particle tracking using geostrophic velocities from AVISO (AL<sub>LA-SSH</sub>, details in Supplementary Fig. 2); (iv) reconstructed following Biastoch et al.<sup>10</sup> from sea surface temperature (AL<sub>SST</sub>) using HadISST data; and (v) inferred from observed drifter and float trajectories (AL<sub>Richardson</sub>, AL<sub>Daher</sub>; the transparent green shading represents the error estimate provided by Daher et al.<sup>32</sup>).</figcaption></figure>

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<sub>LA</sub>, 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)¹⁷, which allows for robust estimates of (sub-)decadal leakage variability: we employ the Parcels tool²¹,²², 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<sub>JRA</sub>, SIM<sub>CORE</sub>, and REA<sub>BRAN</sub> are displayed by the thin black, blue, and cyan lines, respectively, in Fig. 2a.

Agulhas leakage reconstructions based on sea surface temperature observations (AL<sub>SST</sub>, Fig. 1c): Biastoch et al.<sup>10</sup> introduced a method to reconstruct Agulhas leakage variability from sea surface temperature (SST) variability. They analyzed the relation between Agulhas leakage (AL<sub>LA</sub>) 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<sub>LA</sub> 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<sub>SST</sub>). We repeated this analysis in SIM<sub>JRA</sub> and SIM<sub>CORE</sub> (Supplementary Fig. 1). The obtained regressions can explain the major part of the simulated year-to-year leakage variability, as the AL<sub>SST</sub> obtained from applying the regression to the simulated ΔSST correlate significantly (two-sided students t-test with 95% confidence interval) with the original AL<sub>LA</sub>, with *r* values of 0.77 and 0.79 for SIM<sub>JRA</sub> and SIM<sub>CORE</sub>, respectively. When applying the regressions from SIM<sub>JRA</sub> and SIM<sub>CORE</sub> to the timeseries of ΔSST retrieved from HadISST data, we obtain observation-based AL<sub>SST</sub> reconstructions (the thin light-red line in Fig. 2b shows annual AL<sub>SST</sub> transport anomalies based on the regression from SIM<sub>JRA</sub>). Notably, the two reconstructions based on the regressions from SIM<sub>JRA</sub> and SIM<sub>CORE</sub> are hardly distinguishable (Supplementary Fig. 4a). This emphasizes that the AL<sub>SST</sub> reconstruction is much more sensitive to the observed ΔSST difference than to the exact coefficients obtained from the
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.¹⁸ 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²³,²⁴, or passive tracer²⁵,²⁶ and thermohaline threshold²⁷⁻²⁹ 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³⁰. 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 model-based estimates for the mean leakage strength and interannual variability. Observation-based estimates for the mean leakage strength are still poorly constrained¹⁸,³¹,³² and model-based estimates are sensitive to the representation of (sub-)mesoscale processes³³ and the large-scale circulation pattern¹⁴ (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¹⁷. 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¹⁰⁻¹²,³⁴. 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⁸, 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<sub>LA</sub> from SIM<sub>CORE</sub> 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¹⁰. The hypothesis is further supported by evaluating the AL<sub>SST</sub> 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⁴, 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 retroreflection region⁸,³⁵. Yet, there also have been conflicting reports of decoupled Agulhas leakage and variability in eddy kinetic energy¹⁴,²⁵. 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²⁹,³⁴,³⁶. We here make use of our model simulations SIM<sub>JRA</sub> and SIM<sub>CORE</sub> 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 Retroreflection region through the Cape Basin towards the Good Hope section. Hence, we here make use of our simulated Lagrangian trajectories together with a Eulerian eddy-detection 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<sub>CORE</sub> and SIM<sub>JRA</sub>, respectively), we find that leakage particles spend only a bit more than a third of their transit within anticyclonic eddies (35% in SIM<sub>JRA</sub>, 34% in SIM<sub>CORE</sub>) and most of the time outside eddies (59% in SIM<sub>JRA</sub> and 62% in SIM<sub>CORE</sub>). 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<sub>CORE</sub>. Hence, at least for the considered time period and within our model simulations, decadal trends in Agulhas leakage cannot simply be attributed to – and estimated via – changes 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³⁷, 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.³³ 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<sub>SST</sub> estimates based on SIM<sub>JRA</sub> 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¹²⁻¹⁴,³⁴,³⁸,³⁹ and can be linked to the Southern Annular Mode (SAM)¹⁰,³⁴.

Also our model simulations SIM<sub>JRA</sub> and SIM<sub>CORE</sub> 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)⁴⁰ (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<sub>JRA</sub> and SIM<sub>CORE</sub>, 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<sub>JRA</sub> and SIM<sub>CORE</sub> 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⁴¹). Until the mid-1980s, the decadal increase in SAM appears representative for the increase in AL for the two
<figure><img src="image_5.png" /></figure>

**Fig. 3** The increase in Agulhas leakage (AL) between the mid-1960s and mid-1980s cannot directly be related to an increase in the transport via anticyclonic eddies, but is associated with an increase in the Indo-Atlantic exchange outside eddies. a Simulated (SIM<sub>JRA</sub> in black and SIM<sub>CORE</sub> in blue) temporal evolution of (non-)eddy contributions to the interocean exchange, expressed by the mean (averaged over all particles) relative time each particle spent within anticyclonic eddies or outside eddies during its transit between 20°E and GH. b Decadal mean AL vs anticyclonic eddy contribution (calculated based on the annual timeseries of AL<sub>LA</sub> from Fig. 2a and the anticyclonic contribution from panel a of this figure); decadal means for all possible averaging periods are displayed by semi-transparent dots, decadal means for subsequent periods D1 (1965-1984) to D4 (1995-2004) are marked by solid connected dots. c Mean transit time between 20°E and GH. d Mean simulated (SIM<sub>JRA</sub>) leakage pathways for period D1, displayed via crossing probabilities, and changes in the pathways between periods (e) D3 and D1, and (f) D5 and D3.

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<sub>CORE</sub> appears to be larger than suggested by SAM, and the leveling in Agulhas leakage in SIM<sub>JRA</sub> 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<sup>42</sup>, and hypothesized to be the cause for the lack of a significant trend in the AL<sub>SSH</sub> estimate<sup>43</sup>. However, the observation-based AL<sub>SST</sub> estimate (and the AL<sub>LA</sub> estimate from SIM<sub>JRA</sub>) 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<sup>38</sup>. Therefore, we further compared the

temporal evolution of Agulhas leakage to the maximum zonal-mean 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<sub>CORE</sub> 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<sub>JRA</sub> the positive decadal trends remain rather constant. This could help explain the accelerated increase in Agulhas leakage in the 1990s simulated with SIM<sub>CORE</sub> and suggests that the CORE forcing features – at least regionally – a 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<sup>44</sup>). 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<sub>JRA</sub> and reconstructed from observations.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Decadal trends in Agulhas leakage (AL) since 1958 can only partially be related to changes in the annual Southern Annular Mode (SAM) index.
Between the mid-1960s and mid-1980s the increase in AL fits to an overall increase in SAM, afterwards AL evolves disproportionally to SAM. **a** Decadal mean AL vs SAM as simulated in SIM<sub>JRA</sub> (blue) and SIM<sub>CORE</sub> (black) (AL<sub>LA</sub> from Fig. 2a and SAM from panel b of this figure), and as inferred from observations (red) (AL<sub>SST</sub> estimate based on HadISST and SAM from panel b of this figure); decadal means for all possible averaging period are displayed by semi-transparent dots, decadal means for subsequent periods D1 (1965-1984) to D4 (1995-2004) are marked by solid connected dots. Annual (thin lines) and sub-decadally filtered (thick lines, 11 yr Hamming window) (b) SAM index calculated following Marshall (2003) as the difference (40°S-65°S) in normalized zonal-mean sea level pressure from station data (OBS) and the JRA55-do (SIM<sub>JRA</sub>) and CORE (SIM<sub>CORE</sub>) atmospheric forcing data sets; and (c) average Indian Ocean wind stress curl (35-45°S, 20°E-110°E). **d**, **e** Moving centered 21 yr linear trends of annual timeseries in panels b and c (dots mark values that are significant according to student's t-test and 95% confidence level).</figcaption></figure>

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<sup>14,34</sup>. 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<sup>14</sup> rather than strengthening it<sup>9</sup>. Moreover, the rate of change in the position of the westerlies decreases with time, whereby the decrease is more pronounced in SIM<sub>CORE</sub> than in SIM<sub>JRA</sub> (Supplementary Fig. 6c, f). In fact, SIM<sub>CORE</sub> features a stronger poleward shift than SIM<sub>JRA</sub> 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<sub>CORE</sub>, but weakened it in SIM<sub>JRA</sub>. 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. <sup>14,25,38</sup> but in contrast to ref. <sup>45</sup>, we find no robust relationship between the temporal variability in Agulhas leakage and the strength of the Agulhas Current in SIM<sub>JRA</sub> and SIM<sub>CORE</sub> (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<sub>JRA</sub> and SIM<sub>CORE</sub>, respectively), the decadal trends of Agulhas leakage and Agulhas Current transport have the same sign in SIM<sub>JRA</sub> and opposite signs in SIM<sub>CORE</sub>. 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<sub>LA</sub> (r = 0.96 and r = 0.97 for SIM<sub>JRA</sub> and SIM<sub>CORE</sub>, 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<sub>JRA</sub> and SIM<sub>CORE</sub> show a somewhat disproportional response to the increase in regional westerly wind stress: in SIM<sub>JRA</sub>, an increasing Agulhas Current adds to the impact of increasing
<figure><img src="image_7.png" /></figure>

**Fig. 5** Decadal trends in Agulhas leakage (AL) are largely decoupled from changes in the Agulhas Current (AC) and mainly represent a changing fraction of AC waters entering the South Atlantic (AL/AC). a Decadal mean AL vs AL/AC as simulated in SIM<sub>JRA</sub> (blue) and SIM<sub>CORE</sub> (black) (AL<sub>LA</sub> from Fig. 2a and AL/AC from panel c of this figure); decadal means for all possible averaging periods between 1958 and 2005 (overlapping simulation period SIM<sub>CORE</sub> and SIM<sub>JRA</sub>) are displayed by semi-transparent dots, decadal means for subsequent periods D1 (1965-1984) to D4 (1995-2004) are marked by solid connected dots. Simulated annual (thin lines) and sub-decadally filtered (thick lines, 11 yr Hamming window) (b) AC and (c) AL/AC timeseries. **d, e** Moving centered 21 yr linear trends for timeseries in panels b and c of this figure (dots mark values that are significant according to student's t-test and 95% confidence level).

regional westerly winds; in SIM<sub>CORE</sub>, 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<sup>12,15,16,39,46</sup>. 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<sup>47,48</sup> and the Atlantic Multidecadal Oscillation<sup>10</sup>, 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<sup>49</sup>. 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<sub>σ</sub>), as simulated with SIM<sub>JRA</sub> and SIM<sub>CORE</sub> (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<sub>σ</sub> 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<sub>σ</sub> 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<sup>50,51</sup> (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<sup>52</sup>. Moreover, it stresses the need for a more comprehensive assessment of the resolution-dependent 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/sea-ice model configuration with an eddy-rich (1/20°) nest in the South Atlantic, the
<figure><img src="image_8.png" /><figcaption>**Fig. 6** The increase in Agulhas leakage (AL) between the mid-1960s and mid-1980s coincided with an increase in the Atlantic Meridional Overturning Circulation (AMOC), which propagated from the South into the North Atlantic. **a** Annual mean (thin lines) and decadally filtered (thick lines, 21 yr Hamming window) strength of the AMOC in density coordinates ($AMOC_σ$) at 34°S (solid lines) and at 26°N (dashed lines). **b** Annual rate of change in decadally filtered $AMOC_σ$ strength as a function of time and latitude. **c** Annual rate of change in decadally filtered AL strength, as simulated with $SIM_{JRA}$. **d-f** Same as panels a to c but for $SIM_{CORE}$.</figcaption></figure>

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 predecessors - has 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-58}$.

We here made use of three different hindcast simulations: (1) $SIM_{CORE}$ (experiment identifier INALT20.L46-KFS044$^{11}$); forced with the CORE atmospheric data set version 2$^{59,60}$ which is provided as 6-h to monthly fields on a 2° grid for the period 1958–2009, together with climatologically varying runoff; (2) $SIM_{JRA}$ (experiment identifier INALT20.L46-KFS10X$^{17,51}$) similar to $SIM_{CORE}$, but forced with the JRA55-do atmospheric data set version 1.4$^{61}$, which is provided at 3-h resolution on a 1/2° grid for the period 1958–2018, together with interannually varying daily runoff; (3) $SIM_{JRAo}$ (here newly introduced, experiment identifier INALT20.L46-KFS119): similar to $SIM_{JRA}$, but with a different initialization strategy based on the OMIP-2 protocol$^{62}$, starting from rest and initializing temperature and salinities from the World Ocean Atlas 2013 (WOA13$^{63,64}$) instead of starting from the 30-year spin-up used to initialize both $SIM_{CORE}$ and $SIM_{JRA}$, as well as slight changes in the parameter settings regarding the sea surface salinity restoring (now 50 m per yr). Further note that in $SIM_{JRA}$ 11 of the major tidal components are simulated, which were not included in $SIM_{JRAo}$ nor in $SIM_{CORE}$.

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_{\sigma}(t, y) = \max_{\sigma} \Psi(t, y, \sigma); \Psi(t, y, \sigma') = \int_{\sigma_{\text{max}}}^{\sigma} \int_{x_w}^{x_e} V(t, x, y, \sigma') dx \, d\sigma' \quad (1)$$

where $\sigma$ 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

retroreflection 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. They amount to 10.2 ± 2.4 Sv and 9.7 ± 2.1 Sv with corresponding Agulhas leakage to Agulhas Current transport ratios (AL/AC) of 16.4 ± 4.2% and 17.8 ± 3.8% for $SIM_{CORE}$ and $SIM_{JRA}$, respectively; compared to at least 15 Sv (AL/AC: 25%) and 21.3 ± 4.7 Sv (AL/AC: 27.6 ± 2.5%) for observation-based estimates considering the transport in the upper 1500 m$^{31}$ and 2000 m$^{32}$, respectively. Part of the differences may be explained by temporal variability, as the observation-based estimates are mostly representative for the 2000s and 2010s and hence do not include the period of weaker leakage values prior to the mid-1980s. Yet, even if only considering the period post 2000, the mean leakage transport in $SIM_{JRA}$ (10.4 ± 1.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⁶⁸. 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⁶⁹.

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-sam-index-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.

**Data availability**

The station-based annual SAM index can be downloaded from https://legacy.bas.ac.uk/
met/gjima/sam.html. The HadISST 1.1 data is made available through the Met Office,
Hadley Center and can be downloaded from https://www.metoffice.gov.uk/hadobs/
hadisst/data/download.html. Altimeter satellite gridded sea surface height anomalies as
well as derived geostrophic currents are provided by E.U. Copernicus Marine Service
Information: https://doi.org/10.48670/moi-00148, https://doi.org/10.48670/moi-00149.
The updated Al_SSH timeseries are publicly available on zenodo: https://doi.org/10.5281/
zenodo.7022671. Output from the BRAN2020 reanalysis is available from NCI (Australia
National Computational Infrastructure) OPeNDAP servers at https://doi.org/10.25914/
6009627c7af03. The output from the INALT20 ocean model as well as other post-
processed data necessary to reproduce all figures of this manuscript are available through
GEOMAR: https://hdl.handle.net/20.500.12085/3b8db9da-5aa0-476c-81bc-
9b5d1d53186e.

**Code availability**

The Lagrangian software Parcels is available at http://oceanparcels.org/. The code for
Lagrangian leakage estimation in ocean models can be retrieved from https://hdl.handle.
net/20.500.12085/b704e917-09dd-4a73-b6a1-ea24a549920c. The cdftools can be
obtained from https://github.com/meom-group/CDFTOOLS. The scripts necessary to
reproduce the figures of this manuscript are available through GEOMAR: https://hdl.
handle.net/20.500.12085/3b8db9da-5aa0-476c-81bc-9b5d1d53186e.

Received: 6 May 2022; Accepted: 24 November 2022;

Published online: 14 December 2022

**References**

1. Gordon, A. L. Interocean exchange of thermocline water. *J Geophys. Res.* **91**, 5037–5046 (1986).

2. Lutjeharms, J. R. E. *The Agulhas Current*. (Springer, 2006).

3. Lutjeharms, J. R. E. The exchange of water between the South Indian and South Atlantic Oceans; a review. in *The South Atlantic Present and Past Circulation* (ed. Wefer, G.) (Elsevier, 1996).

4. De Ruijter, W. P. M. et al. Indian-Atlantic inter-ocean exchange: dynamics, estimation and impact. *J. Geophys. Res.* **104**, 20820–885910 (1999).

5. Speich, S., Blanke, B. & Cai, W. Atlantic meridional overturning circulation and the Southern Hemisphere supergyre. *Geophys. Res. Lett.* **34**, L23614 (2007).

6. Broecker, W. S. The great ocean conveyor. *Oceanogr* **4**, 79–89 (1991).

7. Beal, L. M., de Ruijter, W. P. M., Biastoch, A. & Zahn, R., SCOR/WCRP/IAPSO working group, 136. On the role of the Agulhas system in ocean circulation and climate. *Nature* **472**, 429–436 (2011).

8. Rouault, M., Penven, P. & Pohl, B. Warming in the Agulhas Current system since the 1980's. *Geophys. Res. Lett.* **36**, L12602 (2009).

9. Biastoch, A., Böning, C. W., Schwarzkopf, F. U. & Lutjeharms, J. R. E. Increase in Agulhas leakage due to poleward shift of the Southern Hemisphere westerlies. *Nature* **462**, 495–498 (2009).

10. Biastoch, A. et al. Atlantic multi-decadal oscillation covaries with Agulhas leakage. *Nat. Commun.* **6**, 10082 (2015).

11. Schwarzkopf, F. U. et al. The INALT family - A set of high-resolution nests for the Agulhas Current system within global NEMO ocean/sea-ice configurations. *Geosci. Model Dev.* **12**, 3329–3355 (2019).

12. Biastoch, A. & Böning, C. W. Anthropogenic impact on Agulhas leakage. *Geophys. Res. Lett.* **40**, 1138–1143 (2013).

13. Ivanciu, I., Matthès, K., Wahl, S., Harlaß, J. & Biastoch, A. Effects of prescribed CMIP6 ozone on simulating the Southern Hemisphere atmospheric circulation response to ozone depletion. *Atmos. Chem. Phys.* **21**, 5777–5806 (2021).

14. Durgadoo, J. V., Loveday, B. R., Reason, C. J. C., Penven, P. & Biastoch, A. Agulhas leakage predominantly responds to the Southern Hemisphere westerlies. *J. Phys. Oceanogr.* **43**, 2113–2131 (2013).

15. Weijer, W., De Ruijter, W. P. M., Dijkstra, H. A. & Van Leeuwen, P. J. Impact of interbasin exchange on the Atlantic overturning circulation. *J. Phys. Oceanogr.* **29**, 2266–2284 (1999).

16. Weijer, W., De Ruijter, W. P. M., Sterl, A. & Drijfhout, S. S. Response of the Atlantic overturning circulation to South Atlantic sources of buoyancy. *Glob. Planet Change* **34**, 293–311 (2002).

17. Schmidt, C., Schwarzkopf, F. U., Rühs, S. & Biastoch, A. Characteristics and robustness of Agulhas leakage estimates: an inter-comparison study of Lagrangian methods. *Ocean Sci.* **17**, 1067–1080 (2021).

18. le Bars, D., Durgadoo, J. V., Dijkstra, H. A., Biastoch, A. & de Ruijter, W. P. M. An observed 20-year time series of Agulhas leakage. *Ocean Sci.* **10**, 601–609 (2014).

19. Chamberlain, M. A. et al. Next generation of Bluelink ocean reanalysis with multiscale data assimilation: BRAN2020. *Earth Syst. Sci. Data* **13**, 5663–5688 (2021).

20. Rayner, N. A. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. *J Geophys. Res.* **108**, 4407 (2003).

21. Lange, M. & van Sebille, E. Parcels v0.9: prototyping a Lagrangian ocean analysis framework for the petascale age. *Geosci. Model Dev.* **10**, 4175–4186 (2017).

22. Delandmeter, P. & Van Sebille, E. The Parcels v2.0 Lagrangian framework: New field interpolation schemes. *Geosci. Model Dev.* **12**, 3571–3584 (2019).

23. Dijkstra, H. A. & De Ruijter, W. P. M. On the physics of the Agulhas Current: steady retroflection regimes. *J. Phys. Oceanogr.* **31**, 2971–2985 (2001).

24. Hermes, J. C., Reason, C. J. C. & Lutjeharms, J. R. E. Modeling the variability of the greater Agulhas Current system. *J. Clim.* **20**, 3131–3146 (2007).

25. Loveday, B. R., Durgadoo, J. V., Reason, C. J. C., Biastoch, A. & Penven, P. Decoupling of the Agulhas leakage from the Agulhas Current. *J. Phys. Oceanogr.* **44**, 1776–1797 (2014).

26. Le Bars, D., De Ruijter, W. P. M. & Dijkstra, H. A. A new regime of the Agulhas Current retroflection: turbulent choking of Indian-Atlantic leakage. *J. Phys. Oceanogr.* **42**, 1158–1172 (2012).

27. Van Sebille, E., Van Leeuwen, P. J., Biastoch, A. & De Ruijter, W. P. M. Flux comparison of Eulerian and Lagrangian estimates of Agulhas leakage: A case study using a numerical model. *Deep Sea Res. 1 Oceanogr. Res. Pap.* **57**, 319–327 (2010).

28. Putrasahan, D. A., Beal, L. M., Kirtman, B. P. & Cheng, Y. A new Eulerian method to estimate 'spicy' Agulhas leakage in climate models. *Geophys. Res. Lett.* **42**, 4532–4539 (2015).

29. Holton, L. et al. Spatio-temporal characteristics of Agulhas leakage: a model inter-comparison study. *Clim. Dyn.* **48**, 2107–2121 (2017).

30. Van Sebille, E. et al. Relating Agulhas leakage to the Agulhas Current retroflection location. *Ocean Sci.* **5**, 511–521 (2009).

31. Richardson, P. L. Agulhas leakage into the Atlantic estimated with subsurface floats and surface drifters. *Deep Sea Res. Part I Oceanogr. Res. Pap.* **54**, 1361–1389 (2007).
32. Daher, H., Beal, L. M. & Schwarzkopf, F. U. A new improved estimation of Agulhas leakage using observations and simulations of Lagrangian floats and drifters. *J. Geophys. Res. Oceans* **125**, 1–17 (2020).

33. Schubert, R., Gula, J. & Biastoch, A. Submesoscale flows impact Agulhas leakage in ocean simulations. *Commun. Earth Environ.* **2**, 1–8 (2021).

34. Loveday, B. R., Penven, P. & Reason, C. J. C. Southern Annular Mode and westerly-wind-driven changes in Indian-Atlantic exchange mechanisms. *Geophys. Res. Lett.* **42**, 4912–4921 (2015).

35. Backeberg, B. C., Penven, P. & Rouault, M. Impact of intensified Indian Ocean winds on mesoscale variability in the Agulhas system. *Nat. Clim. Change* **2**, 608–612 (2012).

36. Cheng, Y., Putrasahan, D., Beal, L. & Kirtman, B. Quantifying Agulhas leakage in a high-resolution climate model. *J. Clim.* **29**, 6881–6892 (2016).

37. Garzoli, S. L. & Gordon, A. L. Origins and variability of the Benguela Current. *J. Geophys. Res. Oceans* **101**, 897–906 (1996).

38. Cheng, Y., Beal, L. M., Kirtman, B. P. & Putrasahan, D. Interannual Agulhas leakage variability and its regional climate imprints. *J. Clim.* **31**, 10105–10121 (2018).

39. Biastoch, A., Böning, C. W. & Lutjeharms, J. R. E. Agulhas leakage dynamics affects decadal variability in Atlantic overturning circulation. *Nature* **456**, 489–492 (2008).

40. Marshall, G. J. Trends in the Southern Annular Mode from observations and reanalyses. *J. Clim.* **16**, 4134–4143 (2003).

41. Fogt, R. L. & Marshall, G. J. The Southern Annular Mode: Variability, trends, and climate impacts across the Southern Hemisphere. *Wiley Interdiscip. Rev. Clim. Change* **11**, 1–24 (2020).

42. Banerjee, A., Fyfe, J. C., Polvani, L. M., Waugh, D. & Chang, K. L. A pause in Southern Hemisphere circulation trends due to the Montreal Protocol. *Nature* **579**, 544–548 (2020).

43. Ivanciu, I., Matthes, K., Biastoch, A., Wahl, S. & Harlaß, J. Twenty-first-century Southern Hemisphere impacts of ozone recovery and climate change from the stratosphere to the ocean. *Weather Clim. Dyn.* **3**, 139–171 (2022).

44. Swart, N. C., Fyfe, J. C., Gillett, N. & Marshall, G. J. Comparing trends in the Southern Annular Mode and surface westerly jet. *J. Clim.* **28**, 8840–8859 (2015).

45. Van Sebille, E., Biastoch, A., Van Leeuwen, P. J. & De Ruijter, W. P. M. A weaker Agulhas Current leads to more Agulhas leakage. *Geophys. Res. Lett.* **36**, L03601 (2009).

46. Van Sebille, E. & Van Leeuwen, P. J. Fast northward energy transfer in the Atlantic due to Agulhas Rings. *J. Phys. Oceanogr.* **37**, 2305–2315 (2007).

47. Kelly, K. A., Drushka, K., Thompson, L., Le Bars, D. & McDonagh, E. L. Impact of slowdown of Atlantic overturning circulation on heat and freshwater transports. *Geophys. Res. Lett.* **43**, 7625–7631 (2016).

48. Castellanos, P., Campos, E. J. D., Piera, J., Sato, O. T. & Silva Dias, M. A. F. Impacts of Agulhas leakage on the tropical Atlantic western boundary systems. *J. Clim.* **30**, 6645–6659 (2017).

49. Lübbecke, J. F., Durgadoo, J. V. & Biastoch, A. Contribution of increased Agulhas leakage to tropical Atlantic warming. *J. Clim.* **28**, 9697–9706 (2015).

50. Yeager, S. & Danabasoglu, G. The origins of late-twentieth-century variations in the large-scale North Atlantic circulation. *J. Clim.* **27**, 3222–3247 (2014).

51. Biastoch, A. et al. Regional imprints of changes in the Atlantic Meridional Overturning Circulation in the eddy-rich ocean model VIKING20X. *Ocean Sci.* **17**, 1177–1211 (2021).

52. Meinen, C. S. et al. Meridional overturning circulation transport variability at 34.5°S during 2009-2017: Baroclinic and Barotropic flows and the dueling influence of the boundaries. *Geophys. Res. Lett.* **45**, 4180–4188 (2018).

53. Madec, G., and the NEMO team. NEMO ocean engine. *Notes du Pôle de Modélisation de l' Institut Pierre-Simon Laplace*, **27**, 396 (2016).

54. Madec, G., Delecluse, P., Imbard, M. & Levy, C. OPA 8.1 Ocean General Circulation Model reference manual. *Notes du Pôle de Modélisation de l' Institut Pierre-Simon Laplace*, **11**, 91 (1998).

55. Fichefet, T. & Maqueda, M. A. M. Sensitivity of a global sea ice model to the treatment of ice thermodynamics and dynamics. *J. Geophys. Res.* **102**, 12609–12646 (1997).

56. Rühs, S., Durgadoo, J. V., Behrens, E. & Biastoch, A. Advective timescales and pathways of Agulhas leakage. *Geophys. Res. Lett.* **40**, 3997–4000 (2013).

57. Rühs, S., Schwarzkopf, F. U., Speich, S. & Biastoch, A. Cold vs. warm water route - sources for the upper limb of the Atlantic Meridional Overturning Circulation revisited in a high-resolution ocean model. *Ocean Sci.* **15**, 489–512 (2019).

58. Wagner, P., Rühs, S., Schwarzkopf, F. U., Koszalka, I. M. & Biastoch, A. Can Lagrangian tracking simulate tracer spreading in a high-resolution ocean general circulation model? *J. Phys. Oceanogr.* **49**, 1141–1157 (2019).

59. Griffies, S. M. et al. Coordinated ocean-ice reference experiments (COREs). *Ocean Model (Oxf)* **26**, 1–46 (2009).

60. Large, W. G. & Yeager, S. G. The global climatology of an interannually varying air-sea flux data set. * Clim. Dyn.* **33**, 341–364 (2009).

61. Tsujino, H. et al. JRA-55 based surface dataset for driving ocean-sea-ice models (JRA55-do). *Ocean Model (Oxf)* **130**, 79–139 (2018).

62. Griffies, S. M. et al. OMIP contribution to CMIP6: Experimental and diagnostic protocol for the physical component of the Ocean Model Intercomparison Project. *Geosci. Model Dev.* **9**, 3231–3296 (2016).

63. Zweng, M. M. et al. World Ocean Atlas 2013, Volume 2: Salinity. In NOAA *Atlas NESDIS* (eds. Levitus, S. & Mishonov, A.) **74**, 39 (2013). https://doi.org/10.7289/V5251G4D

64. Locarnini, R. A. et al. World Ocean Atlas 2013, Volume 1: Temperature. in *NOAA Atlas NESDIS73* (eds. Levitus, S. & Mishonov, A.), 40 p. (2013). https://doi.org/10.7289/V55X26VD

65. Beal, L. M., Elipot, S., Houk, A. & Leber, G. M. Capturing the transport variability of a western boundary jet: results from the Agulhas Current Time-Series Experiment (ACT). *J. Phys. Oceanogr.* **45**, 1302–1324 (2015).

66. Oke, P. R. et al. Evaluation of a near-global eddy-resolving ocean model. *Geosci. Model Dev.* **6**, 591–615 (2013).

67. Russo, C. S., Veitch, J., Carr, M., Fearon, G. & Whittle, C. An intercomparison of global reanalysis products for Southern Africa's major oceanographic features. *Front. Mar. Sci.* **9**, 284 (2022).

68. Chelton, D. B., Schlax, M. G. & Samelson, R. M. Global observations of nonlinear mesoscale eddies. *Prog. Oceanogr.* **91**, 167–216 (2011).

69. Chelton, D. B., Schlax, M. G., Samelson, R. M. & de Szoeke, R. A. Global observations of large oceanic eddies. *Geophys. Res. Lett.* **34**, L15606 (2007).

## Acknowledgements

This study has been funded by the German Federal Ministry of Education and Research through the SPACES-II CASISAC project (grant no. 03F0796A) and by the European Union's Horizon 2020 research and innovation programmes (grant agreement no. 818123, iATLantic). It has also received funding from the Initiative and Networking Fund of the Helmholtz Association through the project "Advanced Earth System Modeling Capacity (ESM)". R.S. acknowledges funding through the project DEEPER (ANR-19-CE01-0002-01) from the French National Agency for Research (ANR). The ocean model simulations were performed at the North German Supercomputing Alliance (HLRN) and on the ESM partition (JUWELS) at Jülich Supercomputing Center (JSC); the Lagrangian trajectory simulations as well as all postprocessing was conducted at the Christian-Albrechts- Universität zu Kiel (system NESH). BRAN data has been used in this study. BRAN is made freely available by CSIRO Bluelink and is supported by the Bluelink Partnership: a collaboration between the Australian Department of Defence, Bureau of Meteorology and CSIRO. We thank Eric Oliver for openly distributing his python implementation of the eddy tracking algorithm, Willi Rath and the Parcel teams for technical support with the Lagrangian analyses, and the GEOMAR ocean dynamics group for their contributions to the development and validation of the employed ocean model configuration.

## Author contributions

S.R. conceived the study, developed the new AL<sub>LA-SSH</sub> estimate together with C.S., F.U.S., and A.B., performed the formal analyses, produced all figures, and wrote the manuscript; C.S. performed the Lagrangian experiments and provided software for the visualization of particle pathways; R.S. provided software for the calculation of the AMOC in density coordinates and performed together with T.G.S. the eddy detection along the particle trajectories; F.U.S. developed the INALT20 model configuration, performed all hindcast simulations, and contributed to the validation of the model output; D.L.B. provided the updated AL<sub>SSH</sub> timeseries; A.B. provided scientific guidance to the analysis under the umbrella of the funded project CASISAC; all authors contributed to the writing through review and editing.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-022-00643-y.

**Correspondence** and requests for materials should be addressed to Siren Rühs.

**Peer review information** *Communications Earth & Environment* thanks Björn Backeberg, Tamaryn Morris and the other anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Regina Rodrigues, Joe Aslin, Heike Langenberg. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_9.png" /></figure>

**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/.

© The Author(s) 2022