Characteristics and robustness of Agulhas leakage estimates: an inter-comparison study of Lagrangian methods

. The inﬂow of relatively warm and salty water from the Indian Ocean into the South Atlantic via Agulhas leakage is important for the global overturning circulation and the global climate. In this study, we analyse the robustness of Agulhas leakage estimates as well as the thermohaline property modiﬁcations of Agulhas leakage south of Africa. Lagrangian experiments with both the newly developed tool Parcels and the well established tool Ariane were performed to simulate Agulhas leakage in the eddy-rich ocean-sea ice model INALT20 (1/20° horizontal resolution) forced by the JRA55-do atmospheric boundary con- 5 ditions. The average transport, its variability, trend and the transit time from the Agulhas Current to the Cape Basin of Agulhas leakage is simulated comparably with both Lagrangian tools, emphasising the robustness of our method. Different designs of the Lagrangian experiment alter in particular the total transport of Agulhas leakage by up to 2 Sv, but the variability and trend of the transport is similar across these estimates. During the transit from the Agulhas Current at 32°S to the Cape Basin, a cooling and freshening of Agulhas leakage waters occurs especially at the location of the Agulhas Retroﬂection, resulting in 10 a density increase as the thermal effect dominates. Beyond the strong air-sea exchange around South Africa, Agulhas leakage warms and saliniﬁes the water masses below the thermocline in the South Atlantic. Competing interests. The authors declare that they have no conﬂict of interest. Acknowledgements. The ocean model simulation was performed at the North-German Supercomputing Alliance (HLRN) and on the ESM partition at Jülich supercomputing centre (JUWELS); the trajectory simulations were conducted at the Christian-Albrechts-Universität 430 Kiel (NESH). We thank the Ariane and Parcels teams for support. This research has been supported by the German Federal Ministry of Education and Research (grant no. SPACES-CASISAC (03F0796A)). It has also received funding from the Initiative and Networking Fund of the Helmholtz Association through the project Earth System Modelling Capacity (ESM)". We thank René for helpful discussions and Willi Rath and Katharina Höﬂich for technical support.

2 Materials and methods 90 2.1 Ocean model simulation Output from a hindcast simulation with the eddy-rich ocean-sea ice model configuration INALT20  was used to conduct offline Lagrangian experiments. INALT20 is based on the ocean general circulation model NEMO (version 3.6; Madec and the NEMO Team, 2014) coupled to the Louvain-La-Neuve sea ice model version 2 (LIM-V2; Fichefet and Maqueda, 1997). A tripolar Arakawa C-grid with a horizontal resolution of 1/4°is used for the global host in which a nest with a 95 horizontal resolution of 1/20°between 63°S-10°N and 70°W-70°E is embedded via two-way nesting using AGRIF (Adaptive Grid Refinement in Fortran; Debreu et al., 2008). In the vertical, INALT20 consists of 46 z-levels with a layer thickness of 6 m near the surface increasing to 250 m in the deepest layers. Cells at the bottom are allowed to be partially filled for a better representation of the bathymetry.
The hindcast simulation used here mainly differs from the experiment described and evaluated by Schwarzkopf et al. (2019) in 100 the applied surface forcing. Here, atmospheric boundary conditions are given by the JRA55-do forcing data set covering the period from 1958 to 2019 (Tsujino et al., 2018). The main difference in this forcing set, compared to the previously applied COREv2 data set (Coordinated Ocean-Ice Reference Experiments data set version 2; Large and Yeager, 2009) is an increased horizontal (from 2°to 0.5°) and temporal (from 6 hourly, daily and monthly, depending on the forcing variable, to 3 hourly) resolution as well as the extension throughout the 2010s. Additionally, river runoff, which has previously been prescribed as 105 a climatological field, is now interannually varying at daily resolution. Furthermore, 11 of the major tidal components are simulated here, while the experiment presented by Schwarzkopf et al. (2019) was non-tidal. All other parameters, including the preceding spin-up, are identical to their experiment "INALT20 NS_RW". Although there are differences between these simulations, the large scale horizontal circulation and also the representation of mesoscale variability in the area of interest to the present study are qualitatively comparable.

Lagrangian experiments
For the comparison of the two Lagrangian methods, Ariane and Parcels, different sets of Lagrangian experiments were conducted and the mean transport, interannual variability and trend from 1958 to 2014 of Agulhas leakage were analysed. In these experiments, particles were released in the Agulhas Current and all particles crossing the approximated Good Hope section in the Cape Basin were referred to as Agulhas leakage (Fig. 1a). To advect the particles, the three-dimensional 5-day mean ve-115 locity fields from INALT20 were used. This temporal resolution of the velocity fields is sufficient for these offline Lagrangian experiments as flow characteristics of the Agulhas Current system do not change significantly when using up to 9-day mean velocity fields (Qin et al., 2014).
In Ariane (Blanke and Raynaud, 1997;Blanke et al., 1999), virtual particles are advected along analytically computed 120 streamlines, which are calculated from the three-dimensional, non-divergent velocity field on an Arakawa C-grid for each output time period. The velocities are linearly interpolated in space, while they are assumed to be constant in time for the given 4 https://doi.org/10.5194/os-2021-33 Preprint. Discussion started: 15 April 2021 c Author(s) 2021. CC BY 4.0 License. output frequency of the model simulation. The interpolation and integration methods employed by Ariane respect the local non-divergence of the input files and the resulting trajectories represent volume transport pathways, which may experience along-track tracer and density changes that reflect the sub-grid-scale parametrisations of the underlying ocean model. In our 125 Lagrangian experiment with Ariane (version 2.3.0_02), which is hereafter referred to as experiment A, particles with an initial maximum transport of 0.1 Sv were released automatically and continuously according to the transport in the Agulhas Current at 32°S over 1 year and advected forward in time for 4 years. The transport assigned to a particle takes the reduced volume of partially filled cells at the bottom into account. In grid cells with a transport smaller than the maximum transport per particle, 1 particle per grid cell was seeded on the v-point at the centre of the 5-day mean model output fields (Fig. 1b). If the transport 130 through a grid cell was greater than the maximum transport per particle, 8 particles were released, with 4 of them at the first quarter of the temporally-averaged model output and 4 particles after the third quarter. In our case, the former related to a release at 6 am the day before the centre of the 5-day mean output and the latter to a release at 6 pm the day after the centre of the 5-day mean output. Such a group of 4 particles was seeded regularly distributed in the x-z-plane centred around the v-point. In partially filled cells at the bottom, particles were also released in the middle of the cell in the vertical and hence 135 above particles at the same depth level in fully filled cells. A total of 57 of these 1-year release experiments with on average 186,000 particles per year were conducted for the years 1958-2014. As the experiments were run in the so called "quantitative" mode, the position, time and transport of a particle was only stored at its release as well as when the particle crossed a predefined section for the first time. All particles not leaving the domain enclosed by these sections within the integration time were considered as "lost".

140
In Parcels (Lange and Sebille, 2017;Delandmeter and van Sebille, 2019), virtual particles are advected numerically with different schemes and customisable advection methods being available. If the velocity and tracer fields are discretised on a Cgrid, the velocity fields are linearly interpolated in space and time, while the tracer fields are linearly interpolated in time, but constant over a grid cell. In our first Lagrangian experiment with Parcels (hereafter referred to as experiment P), particles were 145 3-dimensionally advected forward in time using the 4th-order Runge-Kutta scheme with Parcels version 2.2.0. Although this integration scheme is not necessarily volume conserving, it has been used to estimate volume transport pathways (van Sebille et al., 2012 since without diffusion trajectories calculated with the 4th-order Runge-Kutta scheme and an appropriate time step are very similar to trajectories of the analytical, volume conserving method (van Sebille et al., 2018). When using the 4th-order Runge-Kutta scheme, the time step dt should be small enough that particles do not skip grid cells and thereby 150 miss characteristic oceanographic features. After considering the trade-off between the accuracy of the time stepping and the computation time of the experiment, the timestep dt was chosen such that particles stay in one model grid cell for at least 2 time steps. Therefore, a time step dt of 16 minutes was calculated following dt < ds min /(v max ·2) where ds min is the smallest edge of all grid cells in the domain and v max is the largest horizontal velocity. The time step should also not be smaller than necessary to avoid the accumulation of round-off errors. Similarly to experiment A, particles were released in the Agulhas 155 Current at 32°S over 1 year for the years 1958-2014 and then advected for 4 years. Thereby, potential temperature θ and practical salinity S were sampled and together with the position of the particles stored daily.   and their position was determined based on the transport of the Agulhas Current at 32°S. Therefore, the southward transport per grid cell at 32°S over the whole water column was calculated for each time step of the model output, i.e. every 5 days, and the number of particles per grid cell and their individual transport identified accordingly such that the maximum transport per experiments is also different with the runtime of experiment P being about 7 times longer than for A, but there is potential to increase the model efficiency in Parcels with the aim of a fully parallel version (Delandmeter and van Sebille, 2019). The section crossed by a particle was calculated in a second step such that only the first section crossed by a particle leaving the domain was considered for later analysis.
For the comparison of the effect of different designs of the experiment, in particularly the release section, on the Agulhas 170 leakage transport, a second experiment was performed with Parcels (hereafter referred to as P-ACT). In P-ACT, the particles were released along the ACT section at 34°S in the Agulhas Current (Fig. 1a). The section "East", used to define the transport of the Agulhas Return Current, was thereby shifted westwards to 29°E. Otherwise, experiment P-ACT was conducted exactly the same as experiment P.

175
In all experiments, the transport across the defined sections was calculated by summing up the transport of all particles crossing the respective section and dividing by the number of release times per year. Only the section crossed first was considered here as this is automatically done in experiment A. The Agulhas leakage transport is given by the sum of the transport of all particles crossing the sections "West" and "Northwest", which approximate the Good Hope section.

180
The modifications of the thermohaline properties of Agulhas leakage can be evaluated in terms of a horizontal distribution of the density changes dρ.
where the overbar denotes a time average as follows α is the thermal expansion coefficient, β the saline contraction coefficient, κ the isentropic compressibility and p pressure. As we were only interested in the density changes due to the thermal and the haline effect, only the first and second term in Eq.
(1) were computed. First, the changes in density per day due to the thermal effect (Eq. (3)) and the haline effect (Eq. (4)) were calculated for each particle P along its trajectory. In a second step, all particles passing a certain geographical bin along their way were selected and their local density changes per day weighted by the individual particle's transport T were summed up.

190
Finally, this was divided by the cumulative sum over all particles of their individual transport multiplied by the length of each particle's trajectory in days. Note that particles remaining in a bin longer than a day or crossing a bin several times were also captured repeated times in the summation of the density changes as well as transports.
The bin size is ∆y = 2°in latitude × ∆x = 3°in longitude, which is larger than the maximum distance travelled per day by a particle of 206 km to avoid aliasing.

Agulhas leakage transport
The mean  Lagrangian transport across all defined sections and its variability with respect to the release year agree well for experiments A and P (Fig. 2), despite the two different Lagrangian tools Ariane and Parcels being used and the resulting differences in advection schemes, interpolation and release strategy. In both A and P, the majority of the Agulhas Current transport at 32°S reflects back into the Indian Ocean and constitutes the Agulhas Return Current, thereby crossing 205 the section "East" with a mean transport of 40.1 Sv in experiment A and 40.0 Sv in P. Particles with a combined transport of approximately 1 Sv accounting, for 2 % of the overall transport, do not leave the area within 4 years in both A and P. The Agulhas leakage transport for the years 1958 to 2014 is on average 9.7 Sv in experiment A and 9.9 Sv in P with a similar standard deviation in both experiments of 2.1 Sv. In general, we found a good agreement of the simulated transport across all sections in A and P with its differences between A and P being smaller than the standard deviation due to interannual variability.

210
The simulated Agulhas leakage transport is, however, considerably lower than the estimate of mean Agulhas leakage from floats and drifters of at least 15 Sv in the upper 1000 m (Richardson, 2007)  especially their strengthening (Durgadoo et al., 2013;Biastoch and Böning, 2013;Biastoch et al., 2015;Cheng et al., 2018).
Furthermore, first sensitivity studies indicate that the presence of tides in the experiment we analyse here, in contrast to the non-tidal simulation under COREv2 forcing, might contribute to a weaker mean Agulhas Leakage. Understanding the details of the impact of Southern Hemisphere westerlies and tides on Agulhas leakage and the differences between Agulhas leakage in INALT20 forced with either JRA55-do or COREv2 is beyond the scope of this and subject to a dedicated study.

235
If a daily timeseries of Agulhas leakage is created by summing up the transport of particles that cross the Good Hope section each day, the transport varies considerably on seasonal timescales, which can be attributed to the passage of Agulhas Rings (not shown). The high correlation between A and P of these daily timeseries (r = 0.85, p < 0.01) furthermore underlines the robustness of our method in terms of specific leakage events with a transport of up to 10 Sv.
The latitude at which particles in the Lagrangian experiments are released in the Agulhas Current changes the mean Agulhas 240 leakage transport, but it has a minor effect on the variability of Agulhas leakage (Cheng et al., 2016). The mean  Agulhas leakage transport in P-ACT, where particles were released along the ACT section at 34°S in the Agulhas Current, is 12.2 Sv and therefore 2.3 Sv (23 % of the Agulhas leakage transport in P) higher than in P with a release further north at 32°S (Fig. 2b). This can be partly explained by the increase in the transport of the Agulhas Current between these latitudes due to the increase in Sverdrup transport . Both time series show a comparable interannual variability with 245 a correlation of r = 0.92 (p < 0.01). A similar comparison of the release location in a climate model also revealed a good correlation between the Agulhas leakage transport time series, but an even stronger increase of the mean Agulhas leakage transport of 2.8 Sv (30 %) in the experiment with a release at ACT was found (Cheng et al., 2016). This agrees surprisingly well with our findings given that the atmospheric forcing of the ocean differs due to the intrinsic variability in a climate model leading to a different temporal evolution and variability of Agulhas leakage. The strong nonlinearity of the southern Agulhas

250
Current results in different degrees of recirculation south of ∼34°S and hence the Agulhas Current transport in different models usually agrees at 32°S, but disagrees at ACT (Biastoch et al., 2018). The lower Agulhas leakage transport in P compared to the estimate from observations by Daher et al. (2020) can hence be partly explained by the different reference sections through the Agulhas Current.

255
The transit times of waters from the Agulhas Current at 32°S to the Good Hope section vary strongly from less than a month up to several years, where a maximum of 4 years is the upper limit due to the design of the experiments (Fig. 3a).  interannual variability (r = 0.996, p < 0.01). The transport time series referenced to the release year and the year of the last 285 crossing of the Good Hope section do not correlate well (r = 0.37, p < 0.05), which is expected due to the transit time. The highest correlation between those time series occurs when the transport time series referenced to the last crossing lags 1 year (r = 0.76, p < 0.01). Although the different samplings influence the interannual variability, the decadal variability is not affected and agrees across all time series. Using the date of the Good Hope section crossing as the reference for the Agulhas leakage time series enables to capture crossings of e.g. Agulhas Rings. There is not a big difference between the Agulhas leakage time year is used as a reference. Due to the potential impact of an increasing leakage under global warming on the AMOC (Beal et al., 2011;Biastoch et al., 2015), the overall trend of the Agulhas leakage transport might be more important than the Agulhas leakage transport in a particular year. As a general recipe, we propose to use the year of release as reference for the Agulhas leakage transport time series calculated from the particles crossing the Good Hope section within 4 years. This represents the mean exchange between the Indian and Atlantic Ocean and its temporal changes on decadal timescales towards long-term 300 trends well while keeping the Lagrangian experiment as simple as possible by not introducing another control section at 20°E.
Whether all particles crossing the Good Hope section or only those with an odd number of crossings should be considered depends on the research question. If the advective effect of Agulhas leakage on the overturning circulation is of interest, only particles crossing the Good Hope section an odd number of times and thereby remaining in the South Atlantic shall be referred to as Agulhas leakage. On the other hand, all particles that cross the Good Hope section are important when analysing the effect 305 of Agulhas leakage on the thermohaline properties of water masses in the Southeast Atlantic as mixing with water masses from the Atlantic occurs for example already in the Cape Basin (Rimaud et al., 2012;Rusciano et al., 2012), which is so turbulent that particles can either stay there longer than the integration time of the Lagrangian experiment or escape and eventually join the Agulhas Return Current.

310
In the following, the thermohaline characteristics of Agulhas leakage are analysed for experiment P only as we showed that the average transport, its variability and the transit time of Agulhas leakage agree well between experiments A and P. Here, all particles crossing the Good Hope section, regardless of the number of crossings, are analysed.
On its transit into the Atlantic Ocean, an overall cooling and freshening of Agulhas leakage occurs especially at the location of the Agulhas Retroflection, but the modification of thermohaline properties differs between water masses. The θS-diagram with 26 < σ 0 < 26.8 (Arhan et al., 1999(Arhan et al., , 2011. During the transit from the Agulhas Current (Fig. 5a) towards the Good 320 Hope section (Fig. 5b) the amount of upper waters (σ 0 < 26) reduces by 0.9 Sv to 2.0 Sv, equivalent to a reduction of 30 %.
Thereby, the waters with σ 0 < 25 are mostly eroded by cooling due to very strong surface heat fluxes in this region (Walker and Mey, 1988). At the same time, the amount of central waters increases by 0.5 Sv to 5.0 Sv and the amount of AAIW (26.9 < σ 0 < 27.5) increases by 0.3 Sv to 2.6 Sv. As more than 95 % of the Agulhas leakage transport is within the upper 1200 m, deep waters denser than σ 0 = 27.5 constitute a negligible amount of Agulhas leakage in agreement with observations 325 (Gordon et al., 1987). This confirms the assumption by Richardson (2007) and Daher et al. (2020) that available floats drifting at 1000 m represent Agulhas leakage over the full depth well.
The modifications of temperature and salinity and their effect on density are shown in more detail in Fig. 6 and 7 for each water mass. The upper waters and the AAIW cool continuously during the transit, while the temperature distribution of the central waters broadens especially between 32°S and 20°E (Fig. 6a, d, g). The horizontal distribution of density changes due to northern part of the Cape Basin, but of smaller magnitude compared to the gain. The area of the strongest density changes due to the thermal effect is within the major pathway of Agulhas leakage around South Africa following the Agulhas Current and 335 the Agulhas Ring corridor (Fig. 7a). The density changes of the central waters and the AAIW have a similar pattern, but they have a smaller magnitude for the AAIW, with a density decrease due to warming along the coast of South Africa and a density increase in particular around the Agulhas Retroflection. The mean depth of the σ 0 =26 surface, which is the upper boundary of the central waters, suggests that different dynamical processes are responsible for the density changes of the central waters south of Africa. In the Agulhas Current and Agulhas Retroflection region central waters are mostly found below at least 100 340 m, which is below the mean winter mixed layer depth (indicated by the cyan line in Fig. 7f), suggesting a change in density mainly due to mixing with other water masses. To the northwest, in particular in the Cape Basin, the depth of the central waters shallows to a minimum depth of 50-100 m, which is within the mean winter mixed layer, and therefore air-sea fluxes might modify the thermohaline properties of the shallower central waters. The modifications of the salinity has a smaller effect on the density compared to the thermal effect, but is not negligible especially for the intermediate waters. There is a salinification and 345 hence density gain of the upper waters all around South Africa with freshening only in the southern part of the Cape Basin (Fig.   7e). This is in agreement with a slight shift towards higher salinities between the release and 20°E and a decrease afterwards in Fig. 6b. The central waters and the AAIW experience a freshening, which dominates in the region and partly cancels out the thermal effect on the density. The freshening of AAIW is a result of mixing with AAIW from the Atlantic in the Cape Basin (Rimaud et al., 2012;Rusciano et al., 2012). The changes in density of all water masses of Agulhas leakage combined shows i 26.9 < 0 < 27.5 Figure 6. Thermohaline property modification of Agulhas leakage waters for experiment P: Transport of (a-c) upper waters (22.5 < σ0 < 26), (d-f) central waters (26 < σ0 < 26.9) and (g-i) AAIW (26.9 < σ0 < 27.5) per potential temperature (a, d, g, 1°C bins), salinity (b, e, h, 0.1 psu bins) and σ0 (c, f, i, 0.2 kg m −3 bins) class. To assign waters to these water masses, σ0 in the Agulhas Current at 32°S is used, highlighted by the blue shading. The properties at the release at 32°S are shown as blue, at 20°E as cyan and at the first crossing of the Good Hope section as purple line. The thin lines show the transport of all water masses together including the deep waters (σ0 > 27.5).
that cooling in the Cape Basin has the strongest impact on the density (Fig. 7b). Initially Agulhas leakage waters are relatively warm and saline, which is a density compensated anomaly. As they lose their thermal signature during the transit, a positive density anomaly due to the anomalously high salinity spreads north into the South Atlantic (Weijer et al., 2002;Biastoch et al., 2009). especially their strengthening (Durgadoo et al., 2013;Biastoch and Böning, 2013;Biastoch et al., 2015;Cheng et al., 2018).  (Schubert et al., , 2020. In particular the representation of lee cyclones west of the Agulhas bank improves when submesoscale flows are resolved (Schubert et al., 2021, submitted manuscript). The lee cyclones are important for the formation of Agulhas filaments, which contribute to 40 % more Agulhas leakage, when the resolution is increased to 1/60° (Schubert et al., 2021, 390 submitted manuscript). The absence of submesoscale flows thus can partly explain the too low Agulhas leakage transport in INALT20 compared to observations. Submesoscale dynamics also influence the mode and intermediate water properties by the formation of mode waters in Agulhas Rings, which only occurs if the submesoscale is resolved, and the AAIW transformation in its regional varieties (Capuano et al., 2018).
In terms of thermohaline properties, we showed that they are modified not only in the Cape Basin, but also in other regions other. The impact of Agulhas leakage and its temporal changes can not only be interpreted based on the transport, but also the thermohaline properties of Agulhas waters and their spatial and temporal changes need to be taken into account. A substantial 405 increase in total heat and salt fluxes across the Good Hope section over the last decades is the result of a combined effect of changes in transport and water mass properties of Agulhas leakage (Loveday et al., 2015).
We have confirmed that estimating Agulhas leakage with the Lagrangian tool Parcels leads to very similar results as simulations with Ariane. Even though the numerical integration scheme of Parcels is not volume conserving by definition as the analytical method in Ariane, our results show that Parcels can be used for volume transport estimations. This opens up new 410 opportunities for future studies due to the flexibility and ongoing development of Parcels. Agulhas leakage can, for example, also be estimated with Parcels based on a variety of differently gridded products, e.g. reanalysis products, which do not have a non-divergent velocity field as required by Ariane. Velocity fields of an ocean model with a non-linear free surface can also not be used in Ariane. Another useful feature of Parcels is the ability to build a realistic, global velocity field consisting of model output with different resolutions for different domains. In the case of INALT20, particles would be advected with velocities 415 from the nest with its mesoscale resolving horizontal resolution of 1/20°in the Agulhas Current system and velocities from the global host (1/4°horizontal resolution) everywhere else. Such an experimental design can be used to for example analyse the far-field impact of Agulhas leakage on the AMOC. A disadvantage of Parcels at the moment is the longer computing time compared to Ariane, but the Parcels team is working on an increased efficiency and aims for a parallel version (Delandmeter and van Sebille, 2019).

420
Code and data availability. The Lagrangian software Parcels is available at http://oceanparcels.org/ and Ariane at http://mespages.univbrest.fr/~grima/Ariane/. For reproducibility of all results, the scripts to perform the Lagrangian experiments and the postprocessing scripts as well as all data required to produce the figures are made available through GEOMAR (https://hdl.handle.net/20.500.12085/b704e917-09dd-4a73-b6a1-ea24a549920c).
Author contributions. AB, SR and CS defined the overall research problem and methodology; FUS developed, ran and validated the ocean all Lagrangian simulations, produced all figures and prepared the paper. All coauthors discussed the analyses and contributed to the text.