Eddy Diffusivity Estimates from Lagrangian Trajectories Simulated with Ocean Models and Surface Drifter Data—A Case Study for the Greater Agulhas System

TheLagrangiananalysisofsetsofparticlesadvectedwiththeﬂowﬁeldsofoceanmodelsisusedtostudyconnectivity, that is, exchange pathways, time scales, and volume transports, between distinct oceanic regions. One important factor inﬂuencing the dispersion of ﬂuid particles and, hence, connectivity is the Lagrangian eddy diffusivity, which quantiﬁes the inﬂuence of turbulent processes on the rate of particle dispersal. Because of spatial and temporal discretization, turbulence is not fully resolved in modeled velocities, and the concept of eddy diffusivity is used to parameterize the impact of unresolved processes. However, the relations between observation- and model-based Lagrangian eddy dif-fusivityestimates,aswellaseddyparameterizations,arenotclear.Thisstudypresentsananalysisofthespatiallyvariable near-surface lateral eddy diffusivity estimates obtained from Lagrangian trajectories simulated with 5-day mean velocities from an eddy-resolving ocean model (INALT01) for the Agulhas system. INALT01 features diffusive regimes for dynamically different regions, some of which exhibit strong suppression of eddy mixing by mean ﬂow, and it is consistent with the pattern and magnitude of drifter-based eddy diffusivity estimates. Using monthly mean velocities decreases the estimated diffusivities less than eddy kinetic energy, supporting the idea that large and persistent eddy features dominate eddy diffusivities. For a noneddying ocean model (ORCA05), Lagrangian eddy diffusivities are greatly reduced, particularly when the Gent and McWilliams parameterization of mesoscale eddies is employed.


Motivation
Over the past decades, the Lagrangian analysis of fluid motion by following floating instruments has been used to investigate ocean general circulation patterns (e.g., Davis 1991a,b;Poulain 2001;Lumpkin and Johnson 2013;Lumpkin et al. 2017). Additionally, an increasing number of Lagrangian analyses are performed by advecting virtual fluid particles with the simulated flow fields of ocean models (van Sebille et al. 2018). They are employed in large-scale oceanography to study the sources, fates, and transformations of water masses Lique et al. 2010;Koszalka et al. 2013a,b;Gary et al. 2014;Durgadoo et al. 2017) and are particularly suited to quantify connectivity between different oceanic sites, that is, preferential linking pathways (e.g., Rühs et al. 2013;van Sebille et al. 2013 and associated time scales (e.g., Blanke et al. 2002;van Sebille et al. 2011;Koszalka et al. 2013a;Rühs et al. 2013), as well as volume, freshwater, and heat transports (e.g., Blanke et al. 2001;Biastoch et al. 2008b;Döös et al. 2012).
One important factor influencing large-scale (;1000 km) connectivity is the average rate at which particles disperse. Turbulent processes, such as mesoscale (10-100 km) eddies and jets, cause fluid particles to disperse quickly and increase the rate of mass, momentum, and tracer spreading, leading to accelerated mixing (LaCasce 2008). Turbulent flow is often described via Reynolds decomposition in terms of a long-term (or slowly varying) mean velocity and the residual eddy component. Analogously, the concept of Lagrangian eddy diffusivity is used to quantify the rate of dispersal related to the cumulative effect of eddies. Depending on the definition of the mean flow, the residual eddy component may include not only mesoscale, but also seasonal-tointerannual (Rieck et al. 2015;Laurindo et al. 2017;Uchida et al. 2017) or smaller-scale (if resolved but not low-pass filtered) variability, or may specifically refer to processes not resolved by a certain flow field (Rypina et al. 2016). In this study, we focus primarily on mesoscale eddy variability.
Owing to spatial and temporal discretization, turbulent processes are not fully resolved in modeled velocity fields, but are parameterized instead. Simulated advective Lagrangian trajectories capture the resolved turbulence but only implicitly include the effect of subgrid-scale parameterizations acting on the tracer and momentum equations, that is, by altered large-scale circulation patterns and along-track changes of temperature and salinity. This led to the notion that dispersal of advective Lagrangian trajectories is not sufficiently diffusive compared to particle dispersal in the real ocean. To circumvent this issue, Lagrangian diffusion parameterizations were introduced (Griffa 1996; Berloff and McWilliams 2002;Monti and Leuzzi 2010;Döös et al. 2011). These add an additional stochastic component to the advective particle displacements (or velocities or accelerations) and have already been employed in regional ocean circulation studies (De Dominicis et al. 2012;Koszalka et al. 2013b;Rypina et al. 2016).
However, compared to the vast number of large-scale Lagrangian applications with ocean models, there are few comprehensive validations of eddy diffusivities associated with the simulated particle dispersal. Numerous studies presented eddy diffusivity estimates derived from drifter data (e.g., Krauss and Böning 1987;Davis 1991a,b;Swenson and Niiler 1996;Poulain 2001;Lumpkin and Flament 2001;Bauer et al. 2002;Zhurbas and Oh 2004;Sallée et al. 2008;Koszalka et al. 2011;Zhurbas et al. 2014;Peng et al. 2015). Likewise, several studies addressed the diffusivity estimation using Lagrangian trajectories simulated with velocity fields from eddy-resolving ocean models (e.g., McClean et al. 2002;Koszalka and LaCasce 2010;Griesel et al. 2010Griesel et al. , 2014Chen et al. 2014;Wolfram et al. 2015). Yet there are only a handful of publications aiming at a quantitative comparison of eddy diffusivity estimates from drifter data and simulated trajectories (De Dominicis et al. 2012;Rypina et al. 2012Rypina et al. , 2016. Additionally, the relations between Lagrangian eddy diffusivity estimates and the optimal choice of diffusivities to be used in stochastic Lagrangian parameterizations and in the Eulerian diffusion parameterizations employed in OGCM tracer equations are not well understood (van Sebille et al. 2018). Unresolved issues concern the difficulty to unequivocally define the mean flow and residual eddy component , the spatial variability of eddy diffusivities associated with different turbulence regimes (Berloff and McWilliams 2002;Koszalka et al. 2011), and the sensitivity of model-based eddy diffusivity estimates to the temporal and spatial model resolution (Keating et al. 2011;Wolfram et al. 2015;van Sebille et al. 2018).
In this study, we address this gap by jointly assessing lateral near-surface eddy diffusivity estimates obtained from real drifter data and trajectories simulated with the velocity output from ocean general circulation models (OGCMs) at varying horizontal and temporal resolutions for the greater Agulhas system. The greater Agulhas system, located around the southern tip of Africa, is known for its vigorous eddy activity and its importance for interbasin exchange of heat, salt, and momentum between the Indian and Atlantic Oceans (e.g., Beal et al. 2011). It features different dynamic regimes in a confined region ( Fig. 1): the Agulhas Current (AC), a strong but stable western boundary current in the Indian Ocean; the Agulhas Retroflection (AR) into the eastward-flowing Agulhas Return Current (ARC); and associated shedding of Agulhas eddies (Lutjeharms 2006) that travel into the eastern South Atlantic Gyre (eSAG). Hence, it constitutes a good test region for evaluating spatially variable eddy diffusivity characteristics. In particular, we address the following questions: This manuscript is structured as follows: In section 2, we describe the employed OGCMs, the offline-performed Lagrangian experiments, the applied method for eddy diffusivity estimation, and the primary observational reference dataset. In section 3, we present and discuss the results: eddy diffusivity estimates obtained from trajectories simulated with 5-day mean velocity output of an eddy-resolving OGCM (section 3a), a comparison of those simulation-based estimates to the observationbased estimates (section 3b), and the sensitivity of the simulation-based diffusivity estimates to the temporal and lateral OGCM output resolution, including the impact of an applied GM parameterization (section 3c). In section 4, we further compare our diffusivity estimates with eddy diffusivity estimates based on different methods and discuss possible implications of our results for eddy parameterization approaches. Section 5 lists our conclusions.

Data and methods
One major goal of this study is to compare eddy diffusivity estimates obtained from trajectories simulated with ocean models to those obtained from observed surface drifter data. Therefore, we calculated 2D Lagrangian trajectories resembling surface drifter tracks from the simulated flow fields of an eddy-resolving and a noneddying global OGCM configuration, INALT01 and ORCA05, respectively. Subsequently, we derived eddy diffusivities from those simulated trajectories as well as from the Global Drifter Program (GDP; Lumpkin and Pazos 2007) data following the method employed by Zhurbas et al. (2014). Zhurbas et al. (2014) already presented global maps of drifter-derived eddy diffusivities, and our updated version of their estimates for  Long-term mean (1996Long-term mean ( -2006 Eulerian velocity speed (color shading) and direction (vectors). Velocity vectors for speeds .50 cm s 21 are displayed thick and at half-length compared to vectors for speeds ,50 cm s 21 . (c) Long-term mean (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) Eulerian EKE. Contours are displayed for 50, 100, 300, 700, and 1300 cm 2 s 22 . In all plots, the red dashed frame surrounds the region in which virtual fluid particles were released. The black dashed and solid frames enclose the 58 3 58 and 28 3 28 bins used for the diffusivity calculations and plotting, respectively, for the AC, AR, ARC, and eSAG.

JANUARY 2018
R Ü H S E T A L .
the Agulhas region serves as the primary observational reference in this study.

a. GDP data
The GDP array of satellite-tracked surface drifting buoys is the highest spatial resolution ocean velocity data currently available, and it is the cornerstone of our knowledge of near-surface submesoscale and mesoscale turbulence and turbulent diffusive regimes . Each drifter consists of a surface buoy with a transmitter and a temperature sensor, and a subsurface drogue, centered at 15-m depth measuring near-surface mixed layer currents. For this study, we used the latitude, longitude, and drifter velocity components at 6-h intervals obtained by objective interpolation updated through 2012 (http://www.aoml.noaa.gov/phod/dac/gdp. html; Hansen and Poulain 1996). Only trajectory segments with the drogue attached were considered (Lumpkin et al. 2013), and the trajectory data were lowpass filtered prior to the diffusivity calculations to remove variability in drifter positions and velocities with periods smaller than 2 days mostly caused by inertial oscillations.
b. Lagrangian trajectories from ocean model simulations

1) OCEAN MODEL SIMULATIONS
The global ocean/sea ice model configurations INALT01 and ORCA05 were developed under the DRAKKAR framework (Barnier et al. 2007(Barnier et al. , 2014. They were formulated with the Nucleus for European Modelling of the Ocean (NEMO, version 3.1.1; Madec 2008) and implemented on a horizontal tripolar Arakawa C-grid (Mesinger and Arakawa 1976), which is Mercator-type south of 208N. Both configurations employed have 46 vertical levels, with grid spacing increasing from 6 m at the surface to 250 m at depth and partially filled bottom cells.
ORCA05 is a well-validated global configuration (e.g., Danabasoglu et al. 2014) with a nominal horizontal resolution of 0.58 (;43 km in the Agulhas region), which realistically represents the mean flow and interannualto-decadal variability of the large-scale circulation (Biastoch et al. 2008a). It classifies as a noneddying OGCM configuration, since to fully capture the mesoscale variability in the Agulhas region, horizontal resolutions of 0.258 and finer are needed (Hallberg 2013).
INALT01 is a two-way nested model configuration that is based on the global ORCA05 configuration described above, but it is regionally refined between 508S-88N and 708W-708E to its nominal horizontal resolution of 0.18 (;9 km in Agulhas region) to capture the complex mesoscale dynamics of the greater Agulhas system . Figure 1a shows a snapshot of the daily mean current speed at 15-m depth, illustrating the ability of the model to represent, for example, Mozambique Channel eddies and Agulhas rings. INALT01 has been employed by a variety of studies investigating the local dynamics (Cronin et al. 2013;Loveday et al. 2014) and large-scale impact (Lübbecke et al. 2015;Biastoch et al. 2015) of the greater Agulhas system.
Both ocean/sea ice models were spun up for 20 years before the actual simulations were performed using interannually varying (years 1948-2007) or climatological atmospheric forcing fields from the Coordinated Ocean-Ice Reference Experiments version 2b (CORE v2b; Large and Yeager 2009;Griffies et al. 2009). The evolution of tracers was simulated using a Laplacian isoneutral diffusion operator and the total variance dissipation (TVD; Zalesak 1979) advection scheme. The momentum equations were formulated using a bi-Laplacian lateral diffusion operator and the energyand enstrophy-conserving (EEN; Arakawa and Hsu 1990) advection scheme. Diffusivity and viscosity coefficients vary horizontally according to the local grid size and are specified via their maximum values A ht0 and A hm0 (Table 1), respectively. Further details of the INALT01 and ORCA05 simulations employed in this study are described in Durgadoo et al. (2013). Table 1 contains a summary of the model simulations providing velocity data used for six Lagrangian experiments performed for this study. Note that 1) one simulation with the eddying INALT01 configuration yielded three different Lagrangian experiments that differ only in the temporal resolution of the velocity fields used for the trajectory integration (SIMeddy-1d, SIMeddy-5d, and SIMeddy-1m); and 2) for the noneddying ORCA05 configurations, two simulations and respective Lagrangian experiments were analyzed, one with (SIMpareddy-5d) and one without (SIMnoeddy-5d) GM parameterization of baroclinic eddies (Gent and McWilliams 1990). The GM parameterization mimics the impact of baroclinic eddies on tracer fluxes by adding an extra term to the tracer equation, which represents eddy-induced advection (Gent et al. 1995;Gent 2011); it is determined by the GM coefficient, which is computed from the growth of the baroclinic instability as described in Treguier et al. (1997), and is allowed to vary horizontally and temporally with a cap of 1000 m 2 s 21 .

2) LAGRANGIAN TRAJECTORY INTEGRATION
The trajectory integration was performed using the offline Lagrangian community tool ARIANE (version 2.2.6; Blanke and Raynaud 1997). Lagrangian particles were released homogenously over the greater Agulhas system ( Fig. 1) at ;15-m depth, every 0.58 in both latitudinal and longitudinal directions, every 30 days for 10 years (beginning 1200 UTC 15 January 1996), yielding nearly 900 000 particles for each Lagrangian experiment. Subsequently, particles were advected forward in time for 60 days with the modeled horizontal velocities (for SIMpareddy-5d, the Eulerian velocity combined with the eddy-induced velocity from the GM parameterization was used), whereby particle positions were stored daily. As drogued surface drifters' behavior should be mimicked, particles were kept at the same depth over the whole integration period. No additional stochastic Lagrangian parameterization was employed.
The resulting trajectories represent pathways determined by the resolved Eulerian flow, or, in the case of SIMpareddy-5d, by the resolved Eulerian flow combined with the mesoscale eddy-induced velocity from the GM parameterization. It has been shown that including the mesoscale eddy-induced transport yields a more realistic picture of the mean overturning circulation, particularly in the Southern Ocean, where the eddy-induced meridional velocity counteracts the Eulerian meridional velocity (e.g., Drijfhout et al. 2003). However, while parameterizing the advective effect of baroclinic eddies where they are not resolved by the model, GM also suppresses the explicit generation of physical as well as numerical mesoscale variability (Hallberg 2013) by flattening isopycnals. Thus, it is worth comparing eddy diffusivity estimates from Lagrangian experiments in noneddying ocean model simulations with and without GM parameterization.

c. Lagrangian eddy diffusivity estimation
Lagrangian eddy diffusivities can be estimated from both single-particle and particle-pair/cluster statistics (LaCasce 2008). While particle-pair/cluster statistics are additionally suited to infer the spatial spectra of mesoscale currents, they also require a simultaneous deployment of a large number of particles in a confined area. Because GDP data presently fulfill this requirement only in a few regions and thus cannot be used to infer spatial variability of eddy diffusivity estimates via particle-pair or cluster statistics, we employed singleparticle statistics.

PARTICLE DIFFUSIVITY ESTIMATION
The concept of Lagrangian diffusivity estimation was introduced by Taylor (1922), who determined scalar single-particle diffusivities by integrating the ensemble-mean Lagrangian velocity autocorrelation function: whereby d(t) 5 x(t) 2 x(t 0 ) 5 Ð t t 0 v L (t)dt represents the particle displacement, v L (t) 5 ›x/›t is the Lagrangian velocity, and hÁi L indicates Lagrangian averaging, that is, averaging over an ensemble of particles at a certain time lag t after their release at t 0 . The Taylor approach was formulated in the idealized context of statistically homogeneous, stationary, and isotropic flows. Davis (1991b) extended the Taylor framework for single-particle diffusivity estimation over inhomogeneous and anisotropic flows by introducing the eddy diffusivity tensor. The diffusivity is then calculated for an ensemble of particles originating at a certain location, or within a specified geographical region where the flow can be assumed as locally homogenous, for each tensor component k jk (x, t) individually [ j, k are the indices for the horizontal dimensions in Cartesian space with j, k 2 (1, 2)]. Instead of the absolute, only the residual velocities v 0 L (tjx, t o ) and displacements d 0 (tjx, t o ) for each particle (passing through x at time t 0 ) are considered. The residual velocity is defined as the departure from the local Eulerian mean velocity; the residual displacement is defined as the total displacements minus that due to the Eulerian mean velocity. This is why, strictly speaking, the eddy diffusivity is not a pure Lagrangian, but a mixed Eulerian-Lagrangian statistic (LaCasce 2008).
Single-particle eddy diffusivity estimation thus involves the following steps: 1) binning of particle trajectories on a longitudinal-latitudinal grid, 2) estimation of the mean flow and subsequent calculation of residual Lagrangian velocities and/or displacements, and 3) estimation of the ensemble-mean eddy diffusivity tensor components for each bin (either in fixed coordinate or along-and across-flow directions).
There are various approaches to estimate the components of the single-particle eddy diffusivity tensor k jk (x, t), all deriving from Taylor's and Davis's diffusivity concepts but using different procedures to reduce potential biases caused by the inhomogeneity and nonstationarity of flows, that is, different binning techniques, mean flow definitions, and diffusivity estimators. Lagrangian averaging is performed using geographical bins (e.g., Swenson and Niiler 1996;Sallée et al. 2008), which may additionally overlap and/or be rotated with respect to the velocity variance axis (e.g., Lumpkin and Garzoli 2005;Peng et al. 2015). Alternatively, trajectories may be clustered according to the nearest neighbor distance (Koszalka and LaCasce 2010). The mean flow is estimated either by simple averaging or using advanced techniques such as spline fitting or Gauss-Markov decomposition (e.g., Bauer et al. 2002;Laurindo et al. 2017). Finally, diffusivity is estimated via the half-growth rate of the residual dispersion tensor (e.g., de Verdiere 1983;Oh et al. 2000;Rypina et al. 2012Rypina et al. , 2016, the integral of the residual velocity autocorrelation (e.g., Koszalka et al. 2011;Peng et al. 2015), or the residual velocity displacement tensor (e.g., Davis 1991b; Swenson and Niiler 1996;Oh et al. 2000).
In this study, we follow the approach for diffusivity estimation used by Zhurbas et al. (2014). This approach was first introduced and thoroughly tested by random-flight simulations of Lagrangian trajectories in a sheared flow (Oh et al. 2000) and was later employed in an improved version accounting for the mean flow suppression of eddy diffusivities by Zhurbas et al. (2014). Assuming isotropy in diffusivities/eddy statistics, the approach yields one time lag-and coordinate-dependent scalar lateral eddy diffusivity K(x, t), which is defined as the semisum of the minor principal component of the symmetric part of the Davis diffusivity tensor k davis p2 (x, t) and the halfgrowth rate of the minor principal component of the single-particle dispersion tensor k disp p2 (x, t), thus combining two frequently used approaches for diffusivity estimation: The Davis diffusivity tensor k jk (x, t) is defined as whereby the notation d 0 k (t 0 2 tjx, t 0 ) represents the kth component of the residual displacement for a particle passing through x at time t 0 , obtained from following its trajectory backward in time for the period [t 0 2 t, t 0 ]; and y 0 Lj (t 0 jx, t 0 ) is the residual velocity of that particle at time t 0 .
To avoid rotational eddy fluxes, which are nondiffusive, only the symmetric part of k jk (x, t) is considered, which is here referred to as k davis jk (x, t): The half-growth rate of the single-particle dispersion whereby the notation d 0 k (t 0 1 tjx, t 0 ) represents the kth component of the residual displacement of a particle passing through x at time t 0 , obtained from following its trajectory forward in time for the period [t 0 , t 0 1 t]. The dispersion tensor is symmetric by construction and does not include rotational and advective eddy fluxes.
If diffusivities are estimated from ensembles of trajectories passing through a fixed position x at different times, k davis jk (x, t) yields the true time lag-dependent diffusivity, whereas k disp jk (x, t) can be biased by the shear effect. However, the limited spatial-temporal resolution of surface drifter data only allows for estimating diffusivities from ensembles of particles passing through a finite vicinity of x, in which case both measures can be positively biased by shear. Because Oh et al. (2000) showed that the minor principal components of k disp jk (x, t) and k davis jk (x, t), representing across-flow diffusivities (cf. Peng et al. 2015), are less biased by shear flow, these were chosen as (the only) representative estimates for lateral eddy diffusivities. They were obtained according to Consistently with Zhurbas et al. (2014), we used this approach to estimate lateral eddy diffusivities from all Lagrangian experiments and surface drifter data in overlapping 58 3 58 bins with a 28 offset in longitudinal and latitudinal direction. Residual velocities were calculated relative to climatological monthly mean currents. That means for each bin for each month of the year, all trajectories passing that bin in the respective month were selected. Then, all trajectory positions within the bin were considered as the origin of a pseudotrajectory (overlapping pseudotracks were removed), and individual backward and forward displacements d(tjx, t 0 ), as well as the Lagrangian velocity at the pseudo-origin v L (t 0 jx, t 0 ), were calculated (for the simulated trajectories, Lagrangian velocities were obtained by central time-differencing of the discrete displacements). Ensemble averaging yielded the mean displacement hd(t)i L and the climatological monthly mean velocity hv L (t 0 )i L that were used to derive the residual velocities and displacements as From these, for each set of pseudotrajectories starting in the same month, k davis jk (x, t), k disp jk (x, t), and K(x, t) were estimated; the final diffusivity estimates were obtained as averages of the respective 12 individual estimates.
Further details of the calculation, including smoothing details and small modifications with respect to the original Zhurbas et al. (2014) approach, are discussed in the appendix.
It is important to note the time-lag dependence of diffusivity estimates. Again, following Zhurbas et al. (2014), we distinguish the maximum diffusivity K max (x), defined as the local maximum of K(x, t) within the time lag interval 1 # t # 20 days, and the asymptotic diffusivity K inf (x), defined as the local mean value of K(x, t) for the time lag interval 15 # t # 20 days (Fig. 2). The usage of larger time lags was avoided because the sampling error and potential biases caused by spatial and/or temporal inhomogeneity of the residual velocity field increase with the time lag (Davis 1991b).
From a physical or parameterization point of view, K inf is the sought-after estimate of diffusivity because at those large time lags, the mean residual distance traveled by a fluid parcel in a certain time interval is approximately proportional to the square root of that time interval, just as for a diffusive process in which fluid parcels undergo random walks (cf. discussion of the Lagrangian integral time scale below). Earlier studies used K max (e.g., Oh et al. 2000;Lumpkin and Flament 2001;Sallée et al. 2008), but recently it has been shown that K max can largely overestimate the ''true'' diffusivity due to a suppression of mixing in areas where eddies and mean flow propagate at different speeds (Griesel et al. 2010;Klocker et al. 2012a,b;Wolfram and Ringler 2017a).
In this work, we purposefully assess both K inf and K max because the relation of the two quantifies the strength of diffusivity suppression by mean flow (Klocker et al. 2012a) and thus can be used to diagnose how well the ocean models represent these effects (i.e., the relationship between the mean flow and the eddy propagation speed along the mean flow).

d. Pseudo-Eulerian EKE, Lagrangian eddy time, and length scales
In addition to Lagrangian eddy diffusivities, the pseudo-Eulerian mean eddy kinetic energy (EKE) and Lagrangian eddy integral time T L and length L L scales were calculated in overlapping 58 3 58 bins.
To be consistent with the definition of the eddy diffusivity, for which we only considered the minor principal component (i.e., the across-flow component), the pseudo-Eulerian mean EKE is defined as the minor principal component of the Lagrangian mean residual velocity covariance matrix at zero time lag hy 0 Mixing, and thus, K inf , can generally be represented as scaling with EKE and the time scale T L over which mixing occurs, or equivalently, as scaling with the distance L L a particle would travel with the characteristic root-mean-square eddy velocity y 0 (which equals ffiffiffiffiffiffiffiffiffiffiffi EKE p ) before it mixes with its surroundings: . Analogously, we estimated T L and L L using our asymptotic diffusivity estimate K inf and pseudo-Eulerian mean EKE as Quantities T L and L L represent ''memory'' scales over which Lagrangian residual velocities stay strongly correlated; at time and space scales much larger than T L and L L , respectively, residual single-particle dispersion resembles diffusive spreading.

Results and discussion
a. Eddy diffusivity estimates derived from trajectories simulated with an eddying ocean model In this section, we present and discuss near-surface lateral eddy diffusivity estimates obtained from residual velocities and displacements of Lagrangian trajectories simulated with the 5-day mean velocity output of the eddy resolving INALT01 hindcast experiment (Lagrangian experiment SIMeddy-5d). Figure 2 shows time lag-dependent, ensemble-mean single-particle residual dispersion and eddy diffusivity estimates for four selected 58 3 58 bins representing different dynamic regimes (Fig. 1). Table 2 lists the respective pseudo-Eulerian and Lagrangian statistics. The bins are located in the eSAG, the ARC, the AR, and the AC.

ESTIMATES
The eSAG is an open ocean region characterized by weak mean flow and low EKE. Its dispersion and diffusivity curves (Figs. 2a,e) show the characteristic asymptotic behavior described in classical turbulence theory (LaCasce 2008): at time lags larger than ;20 days, s p2 (t) grows approximately linear with time (meaning the mean residual displacement of a fluid parcel in a certain time interval is approximately proportional to the square root of that time interval), typical for the diffusive regime. Consequently, across-flow eddy diffusivities asymptotically approach a nearly constant value at those time lags, and K max ' K inf ' 1.5 3 10 3 m 2 s 21 . Because the curves for s yy (t) and s p2 (t) are hardly distinguishable, and even s xx (t) shows a similar behavior, we conclude that the spreading is approximately isotropic. In this particular case, with weak influence of mean flow and nearly isotropic particle spreading, the diffusivity estimates k disp p2 (t), k davis p2 (t), and K(t) yield nearly the same results, with only a small spread toward larger time lags.
The ARC is a strong alongfront current with pronounced meridional excursions, supposedly related to changes in bottom topography, and vigorous mesoscale variability (Lutjeharms and Ansorge 2001). In the ARC region, dispersion and diffusivity curves (Figs. 2b,f) differ substantially from those described for the eSAG. They can be considered typical of a region where eddies propagate westward relative to an eastward mean flow (Boebel et al. 2003;Chelton et al. 2011). On one hand, the instability of the ARC introduces high EKE, leading to an overall increased across-flow eddy diffusivity. On the other hand, the across-flow eddy diffusivity is suppressed by the mean flow, as shown conceptually for the ACC by Klocker et al. (2012a). The obtained across-flow eddy dispersion estimate s p2 (t), which is nearly identical to s yy (t), increases quickly at short time lags but then levels off before finally increasing again approximately linearly with time. This yields across-flow eddy diffusivity curves to peak (and slightly oscillate) at short time lags before reaching their asymptotic values; the maximum diffusivity strongly overestimates the asymptotic diffusivity, K max ' 12.2 3 10 3 and K inf ' 3.0 3 10 3 m 2 s 21 . Still, the diffusivity estimates k disp p2 (t) and k davis p2 (t) yield nearly the same results, introducing only a small spread around K(t) toward larger time lags.
The dispersion and diffusivity curves for the AR (Figs. 2c,g) show a very similar behavior to that described for the ARC, but with even higher diffusivity TABLE 2. Statistics for four selected 58 3 58 bins centered on the listed coordinates: pseudo-Eulerian mean speed and mean EKE, as well as asymptotic K inf and maximum K max near-surface lateral eddy diffusivity estimates, and Lagrangian integral time T L and length L L scales obtained from OBS and Lagrangian experiment SIMeddy-5d.

Region
Lagrangian experiment Speed (cm s 21 ) EKE (cm 2 s 22 ) Diffusivity estimates (10 3 m 2 s 21 ) values, related to a regional maximum in EKE (cf. Fig. 1d). The across-flow eddy diffusivity curves peak at short time lags before reaching their asymptotic values. Thus, as for the ARC, the maximum diffusivity overestimates the asymptotic diffusivity, K max ' 21.8 3 10 3 and K inf ' 11.2 3 10 3 m 2 s 21 . However, the result for K inf should be interpreted with caution, because even though k disp p2 (t) and k davis p2 (t) show a similar temporal evolution, their spread around the combined estimate K(t) increases after they reach their maximum value. Their respective asymptotic diffusivity estimates K disp inf and K davis inf are 8.1 3 10 3 and 14.4 3 10 3 m 2 s 21 , respectively. The AC, heading southwestward, is one of the strongest western boundary currents of the World Ocean. However, in contrast to other western boundary currents, such as the Kuroshio or the Gulf Stream, it is remarkably stable (Lutjeharms 2006(Lutjeharms , 2007 and has lower EKE. The dispersion and diffusivity curves obtained for the AC (Figs. 2d,h) most resemble those of the eSAG, without strong suppression of eddy mixing by mean flow. This indicates that local eddies propagate at approximately the same speed and in the same direction of the mean flow, which is indeed the case for Natal pulses and Mozambique eddies traveling within and at the border of the AC (Schouten et al. 2002). In contrast to open ocean currents like the AR/ARC, where eddies can propagate westward separately from the mean flow, westward drift is impossible in the AC due to topographic constraints. However, at larger time lags, s p2 (t) does not increase linearly. Consequently, k disp p2 (t) does not show an asymptotic behavior but steadily increases with increasing time lag, suggesting that the mean flow was not successfully removed. Moreover, k davis p2 (t) shows a different behavior: it reaches its maximum at a time lag of ;5 days and afterward slowly decreases with increasing time lag. Neither the asymptotic nor the maximum diffusivity can be determined unambiguously.
In general, the results for the eSAG and ARC (and, with some limitations, also those for the AR) suggest that even in a complex eddying flow system, as represented by the greater Agulhas system, with strong eddymean flow interaction, asymptotic diffusive regimes can be found. However, the results for the AC also highlight the limitations of a generalized binning method to quantify lateral eddy dispersion and diffusivity. Because of the spatial and probably also the temporal inhomogeneity of the residual velocity field, the eddy dispersion does not reach a diffusive regime everywhere, and the derived eddy diffusivity estimates are sensitive to the applied method.
We note that the time lag interval chosen for the calculation of asymptotic diffusivity estimates adopted from Zhurbas et al. (2014) for their global analysis of surface drifters is not perfectly suited for the simulated trajectories in the greater Agulhas system because the diffusivities do not yet converge in all cases. In particular, for the eSAG bin, the Lagrangian integral time scale is ;11.6 days (Table 2), indicating that the diffusive regime is only reached at substantially larger time lags (section 2d; LaCasce 2008). Indeed, convergence seems to be reached at time lags between 25 and 30 days. Yet, as changing the averaging interval does not impact the estimates substantially, we keep them as in Zhurbas et al. (2014) for the sake of comparison. Figures 3a, 3d, and 3g show the spatial pattern of asymptotic eddy diffusivity estimates for the whole greater Agulhas system. The individual estimates K davis inf and K disp inf show the same general pattern and magnitude as the combined estimate K inf : the highest eddy diffusivities are found around the AR region (cf. Fig. 2g), an area with strong unstable mean flow and high EKE (Figs. 1c, 5c). Relatively high eddy diffusivities also occur in regions with weaker background flow but still high EKE, such as in the Cape Basin around the major pathway for Agulhas rings (Schouten et al. 2000;Dencausse et al. 2010), and in the region west and south of Madagascar, where Mozambique Channel eddies and southeast Madagascar eddies propagate (Schouten et al. 2002). In the strong but relatively stable Agulhas Current, eddy diffusivities are lower due to high mean kinetic energy but comparatively low EKE. In the vicinity of the ARC (cf. Fig. 2f), diffusivities are further decreased despite locally high EKE (Figs. 1c, 5c). Lowest diffusivities occur around the low-energy eSAG region.

2) SPATIAL PATTERN OF EDDY DIFFUSIVITY ESTIMATES
The spatial pattern of the maximum diffusivity estimates, as exemplarily shown for the combined estimate K max in Fig. 4a, is similar to that of the asymptotic diffusivity estimates. However, despite a general increase in magnitude, the spatial pattern of K max better resembles that of EKE. In particular, the region around the ARC shows strongly elevated diffusivities . The pronounced differences between K max and K inf in this region can be interpreted as the imprint of eddy mixing suppression by mean flow, as already discussed exemplarily for the ARC bin (section 2a). The relation K max /K inf displayed in Fig. 4d highlights all areas where K max greatly overestimates K inf and, thus, where eddy mixing suppression is strong. The region around the ARC clearly stands out, but eddy mixing suppression is also apparent in the northern part of the Antarctic Circumpolar Current (ACC). This fits well to the theory of eddy mixing because both regions feature an eastward mean flow and the possibility for westward propagation of eddies.
Averaging over all spatial bins yields diffusivities of ;3.3 3 10 3 , ;3.5 3 10 3 , and ;3.4 3 10 3 m 2 s 21 for K davis inf , K disp inf , and K inf , respectively. Even though the spatially averaged K disp inf constitutes a slightly higher estimate than the spatially averaged K davis inf , local differences between K davis inf and K disp inf vary in sign and magnitude. The largest discrepancies are found in the AR and AC regions, where diffusivities differ by up to more than 5.0 3 10 3 m 2 s 21 (Table 2). Furthermore, local discrepancies in K inf are higher than those in K max (not shown), highlighting the difficulty in capturing the convergence of the diffusivity estimates. Even though the approach introduced by Oh et al. (2000) and formulated in its refined form by Zhurbas et al. (2014) seems to be applicable in most regions, the combined diffusivity estimate K(x, t) should be interpreted with caution for areas with big discrepancies between the two individual diffusivity estimates.

b. Comparison of eddy diffusivity estimates derived from simulated trajectories and drifter data
In this section, we compare diffusivity estimates obtained from Lagrangian experiment SIMeddy-5d with those obtained from drifter data (OBS). A comparison with other observation-based eddy diffusivity estimates calculated with alternative techniques is given in section 4.
For nearly the whole greater Agulhas system, pseudo-Eulerian mean EKE (Figs. 5c,d,f and Tables 2, 3) is lower in SIMeddy-5d than in OBS, indicating that the on-average slightly lower diffusivity values in SIMeddy-5d are mainly related to a weaker mesoscale variability. Averaging over all spatial bins yields 197 and 307 cm 2 s 22 for SIMeddy-5d and OBS, respectively. Possible reasons for this discrepancy are the following: First, there are known weaknesses of the model configuration. The EKE in the open ocean is too low in most model simulations due to unresolved or underresolved processes, and it is further decreased in ocean model simulations forced with FIG. 3. Spatial pattern of asymptotic near-surface eddy diffusivity estimates. Shown are results obtained from Lagrangian experiment SIMeddy-5d (top panels) and OBS (middle panels), as well as the relative difference between simulations and observations (bottom panels). Asymptotic eddy diffusivity estimates from (left) the half-growth rate of the minor principal component of the single-particle dispersion tensor K disp inf ; (center) the minor principal component of the symmetric part of the Davis diffusivity tensor K davis inf ; and (right) the combined lateral eddy diffusivity estimate K inf . The relative difference is defined as (SIMeddy-5d 2 OBS)/OBS; respective contours are displayed at a distance of 0.5.

JANUARY 2018
R Ü H S E T A L .
relative winds, such as those employed in this study, due to enhanced surface drag (Eden and Dietze 2009). A limited representation of the depth-dependent Ekman drift in view of a too-scarce vertical grid spacing (6 m at the surface) of the model may have further contributed to lowered EKE values (this follows from the fact that the kinetic energy of a layer with velocity shear always exceeds the kinetic energy of a layer of the same thickness and momentum, but with uniform velocity). Second, the discrepancy between OBS and SIM could be partially related to the nature of surface drifters, which, despite being drogued at 15-m depth, do not perfectly follow the local current, but experience some additional windinduced drift (Lumpkin and Pazos 2007;Poulain et al. 2009). Finally, the different spatiotemporal coverage of drifter data and simulated trajectories matters (see online supplementary material). However, keeping in mind that one-to-one comparisons of drifters and simulated trajectories are of limited use for nonassimilative models (cf. van Sebille et al. 2009), a much better agreement between SIMeddy-5d and OBS seems unlikely, even if the sampling of simulated trajectories would be adjusted to fit the one of drifter data. The choice of the overall time period also seems to be of minor importance because diffusivity estimates obtained from Lagrangian experiment SIMeddy5d hardly differ from Lagrangian experiment SIMeddy-5d-clim (Figs. 6a, 7 and Table 3), which is identical to SIMeddy-5d, except for the fact that the OGCM was forced with climatological atmospheric fields instead of interannually varying ones. This indicates only a minor impact of interannual-to-decadal variability on the diffusivity estimates-at least in the region of interest.
Despite lower mean diffusivities in SIMeddy-5d compared to OBS, we cannot conclude that trajectories simulated with SIMeddy-5d are generally not sufficiently diffusive because locally, the difference between SIMeddy-5d and OBS strongly varies in magnitude and sign. FIG. 4. Spatial pattern of maximum eddy diffusivity estimates and the relation of maximum to asymptotic eddy diffusivity estimates. Shown are results for the combined lateral eddy diffusivity estimates obtained from Lagrangian experiment SIMeddy-5d (top panels), and OBS (middle panels), as well as the relative difference between simulations and observations (bottom panels). (left) Maximum eddy diffusivity estimate K max , included here as an approximate measure for the unsuppressed diffusivities. (right) Relation of K max to the sought-after asymptotic diffusivity estimate K inf , quantifying the strength of eddy diffusivity suppression (color shading) and pseudo-Eulerian mean velocity (vectors; cf. Fig. 6). The relative difference is defined as (SIMeddy-5d 2 OBS)/OBS, respective contours are displayed at a distance of 0.5.
Because the general pattern of the difference between SIMeddy-5d and OBS is robust with respect to the method employed for diffusivity estimation (Figs. 3c,f,i), a closer investigation of this pattern is justified. One can distinguish regions where 1) K max and K inf are lower in SIMeddy-5d than in OBS, such as in the eSAG ( Table 2) or south of Madagascar; 2) K max and K inf are higher in SIMeddy-5d than in OBS, such as in a broad region west of the retroflection; 3) K max is lower but K inf is higher in SIMeddy-5d than in OBS, such as northwest of the retroflection; and 4) K inf is lower but K max is higher in SIMeddy-5d than in OBS, such as in parts of the ARC (Table 2) and at the southern boundary of the study region. In the third and fourth cases, the local relation K max /K inf is substantially altered (Figs. 4d-f) TABLE 3. Summary of spatial mean pseudo-Eulerian and Lagrangian statistics for the greater Agulhas system obtained from OBS and the six Lagrangian experiments performed in this study: near-surface pseudo-Eulerian mean speed and EKE, as well as asymptotic K inf and maximum K max near-surface lateral eddy diffusivity, and Lagrangian integral time T L and length L L scales. The spatial mean has been calculated by averaging over all 416 spatial bins.
Lagrangian experiment Speed (cm s 21 ) EKE (cm 2 s 22 ) Diffusivity estimates (10 3 m 2 s 21 ) differs (e.g., in the region around the ARC, eddy mixing suppression is higher in SIMeddy-5d than in OBS). This is most likely related to the fact that the local pseudo-Eulerian mean flow sampled by SIMeddy-5d clearly features a coherent eastward current, whereas in OBS, the mean flow pattern is less clear (Figs. 5a,b). We conclude that differences in the near-surface lateral eddy diffusivity estimates between SIMeddy-5d and OBS cannot be directly linked to differences in EKE, but they are also influenced by differences in Lagrangian eddy length and time scales related to the suppression of eddy mixing. Even though local differences between the diffusivity estimates from SIMeddy-5d and OBS are strongly nonuniform and can be relatively high, the frequency distributions of K inf are remarkably similar ( Fig. 7; Table 3). Moreover, differences between SIMeddy-5d and OBS seem to fall in the uncertainty range. Differences between the two methods to estimate lateral diffusivities in SIMeddy-5d are of the same magnitude as differences between SIM and OBS (cf. section 3a; Table 2); differences between various observational approaches estimating near-surface lateral eddy diffusivity-for example, in the ACC region-are of comparable magnitude or even bigger (Klocker et al. 2012b and references therein). Thus, we conclude that the Lagrangian trajectories simulated with the 5-day mean output of the high-resolution ocean model INALT01 (SIMeddy-5d) capture Lagrangian eddy diffusivity characteristics of real drifter data fairly well, without the need for additional random-walk diffusion.
c. Sensitivity of eddy diffusivity estimates derived from simulated trajectories to lateral and temporal resolution The spatial OGCM resolution determines to which degree eddies are explicitly resolved during the model simulation. The temporal model output resolution further restricts the scale of processes captured by offline simulated advective Lagrangian trajectories. In this section, we quantify the sensitivity of the Lagrangian eddy diffusivity estimates to the temporal (5-day mean vs daily and monthly mean) model output resolution and spatial (1/108 vs 1/28) OGCM resolution, whereby we focus on K inf .

RESOLUTION
In general, studies of oceanic features from the mesoscale to basin scales require a temporal model output FIG. 6. Sensitivity of the asymptotic eddy diffusivity estimates K inf to OGCM choices. Shown are results for the combined lateral eddy diffusivity estimates obtained from the Lagrangian experiments performed with the eddy-resolving OGCM configuration INALT01 (a) forced with climatological atmospheric fields (SIMeddy-5d-clim), and (b) forced with interannually varying atmospheric fields, but using monthly mean velocities (SIMeddy-1m). Also shown are the results from the Lagrangian experiments performed with the noneddying OGCM configuration ORCA05 (c) without (SIMnoeddy-5d), and (d) with GM parameterization of baroclinic eddies (SIMpareddy-5d). resolution no longer than a few days. Using offline Lagrangian analyses within a high-resolution ocean model (1/108), Qin et al. (2014) showed that the connectivity along six major currents does not significantly change for temporal resolutions between 3-and 9-day means. However, degradation in flow characteristics occurs at lower temporal resolutions, resulting in changed transit times and transports between selected upstream and downstream sections.
The comparison of our diffusivity estimates from SIMeddy-5d and SIMeddy-1d shows that 5-day mean output is sufficient to capture the cumulative effect of eddies on particle trajectories in the greater Agulhas system; using daily mean output does not substantially change the spatial distribution (not shown) and magnitude of K inf and the Lagrangian integral time and length scales (Table 3).
SIMeddy-1m yields the same spatial pattern of K inf as SIMeddy-5d (Fig. 6b), but with a generally reduced magnitude (exception: eddy diffusivities in the northern part of the ACC and in near-coastal parts of the AR are increased compared to SIMeddy-5d), resulting in lower spatial mean values and a slightly smaller spatial variability ( Fig. 7; Table 3). Even though eddy diffusivities are reduced, the spatial mean K inf still reaches ;3.0 3 10 3 m 2 s 21 (compared to ;3.4 3 10 3 m 2 s 21 in SIMeddy-5d). This suggests that monthly mean velocity fields do capture a substantial amount of eddy variability-at least in the greater Agulhas system. Yet, EKE is reduced by 26% (spatial average of 145 cm 2 s 22 in SIMeddy-1m compared to 197 cm 2 s 22 in SIMeddy-5d), and average T L and L L are increased (Table 3). The disproportional reduction of K inf by only 12% again highlights that variability and changes in K inf cannot be directly explained by changes in EKE. The results further support the idea that larger and more persistent eddy features are the dominant factors determining eddy diffusivities (Wolfram et al. 2015).
These results agree with Qin et al. (2014), who showed that the connectivity transports in the AC and AR revealed no significant changes for monthly mean fields compared to 3-day mean fields (in contrast to other investigated regions). Likewise, Biastoch et al. (2015) and Cheng et al. (2016) reported interannual-to-decadal variability in Agulhas leakage transport to be captured by monthly mean data. Qin et al. (2014) also noted that the mean transit time for the AC does change significantly for monthly mean fields compared to 3-day mean fields. This can be explained by the fact that in the greater Agulhas system, connectivity transports, as well as eddy diffusivities, are dominated by the largest and most persistent mesoscale features, which are associated with the largest transports (and constitute some of the largest eddies of the World Ocean), while transit time distributions are prone to be influenced by eddy variability at all scales.

2) SENSITIVITY TO LATERAL OGCM RESOLUTION
By the criteria of Hallberg (2013), the ORCA05 model configuration is noneddying in the extended Agulhas Current system, whereas INALT01 is eddy resolving (cf. section 2a). We already showed that SIMeddy-5d (as well as SIMeddy-1d and, with slight restrictions, SIMeddy-1m) captures the general features of observed eddy diffusivities. In principle, eddy diffusivities (as well as EKEs) obtained from noneddying simulations should approach zero. In reality, the scale separation between eddies and mean flow is imperfect. Depending over which time and space scales the mean is defined, some lower-frequency variance remains in the FIG. 7. Box plots of the spatial distributions of (a) the asymptotic eddy diffusivity estimates K inf and (b) pseudo-Eulerian mean EKE. Shown are results for the main Lagrangian experiments employed in this study, whereby K inf was inferred from the combined lateral eddy diffusivity estimate. On each box, the central mark indicates the median calculated out of the 416 overlapping 58 3 58 spatial ocean bins for the greater Agulhas system; the left and right edges of the box indicate the 25th and 75th percentiles, respectively; and the whiskers extend to the minimum and maximum values. The dashed vertical lines indicate the mean values (cf. Table 3). residual velocities and may lead to nonzero eddy diffusivity estimates, even in so-called noneddying ocean models. Furthermore, numerical noise introduces artificial high-frequency variability. Our analysis for SIMnoeddy-5d yields a spatial mean pseudo-Eulerian EKE of 53 cm 2 s 22 and corresponding K inf of ;1.4 3 10 3 m 2 s 21 (Table 3). Hence, K inf in SIMnoeddy-5d (Fig. 6c) is, on average, 2.0 3 10 3 m 2 s 21 weaker than in SIMeddy-5d (cf. section 3a) but still shows the rough spatial pattern from SIMeddy-5d and does not reach values close to zero. Again, the reduction in EKE (;73%) is larger than that in K inf (;58%). The reduction of estimated lateral eddy diffusivities is consistent with the reduction of estimated isopycnal eddy diffusivities of ;50%-60% reported by Wolfram et al. (2015), based on simulations with an idealized ocean model at 32-km resolution compared to simulations at 8-km resolution (their Table 2).
Using the noneddying OGCM configuration with online GM parameterization (SIMpareddy-5d) yields spatial mean pseudo-Eulerian EKE and eddy diffusivity estimates close to zero (Fig. 6d). This is because while parameterizing the advective effect of mesoscale eddies where they are not resolved by the model, GM also suppresses the explicit generation of physical as well as numerical mesoscale variability (Hallberg 2013;cf. section 2b).
In summary, Lagrangian trajectories simulated with the output from our coarse-resolution OGCM configuration do not sufficiently capture the effect of mesoscale eddies on particle dispersal, which results in too-low eddy diffusivity values compared to estimates based on observations and eddy-resolving model simulations. Lagrangian eddy diffusivity estimates are further reduced by the employed GM parameterization. Nevertheless, for Lagrangian particle simulations, the combined Eulerian and eddy-induced velocities from a model simulation with GM may be favored over velocities from coarse-resolution model simulations without GM because they are supposed to represent large-scale circulation patterns and the advective part of particle dispersal more realistically (cf. section 2b). To capture the effect of mesoscale turbulence, one could fill the gap between coarse-resolution simulations and observations/high-resolution simulations by stochastic Lagrangian parameterizations (cf. section 4).

Comparison with other eddy diffusivity estimates and implications for eddy parameterizations
a. Comparison with other eddy diffusivity estimates for the greater Agulhas system Peng et al. (2015) estimated Lagrangian eddy diffusivity for the Indian Ocean from surface drifter data by using a binning technique similar to Zhurbas et al. (2014) but deriving diffusivities from the autocorrelation of residual velocities, which were calculated using the Gauss-Markov method (Lumpkin and Johnson 2013). They concluded that their diffusivity estimates (with seasonal effects removed) generally agree with the asymptotic estimates of Zhurbas et al. (2014), which can be confirmed by a closer inspection of Fig. 7 in Peng et al. (2015) and our Figs. 3g and 3h. In the region south of Madagascar, our estimates for K inf between 2.0 3 10 3 and 12.0 3 10 3 m 2 s 21 correspond well to their asymptotic minor principal component diffusivity estimate k 2 of 11.5 3 10 3 m 2 s 21 for region D7 (their Table 2).
Our results also match the Osborn-Cox diffusivity calculated from the evolution of passive tracers (simulated with surface velocities derived from AVISO sea surface height data) by Abernathey and Marshall (2013), who reported diffusivity values between 1.0 and 10.0 3 10 3 m 2 s 21 for the greater Agulhas system (their Fig. 2). The overall agreement in magnitude and spatial pattern of our particle-based and their tracer-based diffusivity estimates is encouraging. It supports the applicability of the results of Klocker et al. (2012b) and Wolfram and Ringler (2017b), showing that particleand tracer-based diffusivities are similar in a simple zonal channel flow, to the complex real oceanic circulation of the greater Agulhas system.

b. Implications for stochastic parameterizations in offline trajectory calculations
To derive more diffusive pathways capturing the effect of unresolved turbulence even from coarseresolution OGCMs, stochastic Lagrangian parameterizations can be employed in offline trajectory calculations (cf. section 1). However, it is still an open question how to robustly specify an appropriate stochastic model and how to fit the associated parameters with respect to the spatial and temporal resolutions of any given Eulerian ocean model output.
Though we did not implement a stochastic parameterization here, our study reveals important aspects to be considered: 1) because asymptotic diffusive regimes could be identified for different dynamic regimes in the greater Agulhas system, stochastic Lagrangian parameterizations, as described in Griffa et al. (1995), may indeed be appropriate to mimic the effect of mesoscale turbulence in this region; 2) stochastic Lagrangian parameterizations should account for the spatial variability of the diffusivity parameter, which does not necessarily scale with variability in EKE; and 3) the choice of parameters in the stochastic Lagrangian parameterizations should not only account for the horizontal (and temporal) resolution of the underlying OGCM, but should also consider possible subgrid-scale parameterizations, such as GM.
Based on these considerations, we would like to emphasize that the application of a stochastic Lagrangian parameterization implies changing from a pure Lagrangian analysis of the Eulerian OGCM model output to a Lagrangian modeling approach of particle dispersion, which may break consistency with the underlying OGCM physics and thus-depending on the scientific question-may not always be desired.

c. Relation between Lagrangian eddy diffusivity estimates and diffusivity coefficients of OGCM diffusion parameterizations
The Lagrangian eddy diffusivity is a quantitative diagnostic for the cumulative effect of eddies on the Lagrangian dispersal, and, as such, it is used in this paper to assess simulated particle dispersal. However, other studies estimating Lagrangian eddy diffusivities from high-resolution OGCM output or drifter data are partially motivated by the desire to determine more realistic parameter values for Eulerian diffusion parameterizations in coarse-resolution OGCMs. It is, therefore, worthwhile to compare our Lagrangian diffusivity estimates with the diffusivity coefficients used in those parameterizations, such as the Redi coefficient A ht (Redi 1982). The coefficient A ht should quantify the effect of scales not resolved in the OGCM (i.e., variability occurring below 0.18 and 0.58 for INALT01 and ORCA05 experiments, respectively) on tracer mixing rates. Important for this context, Lagrangian eddy diffusivities derived from the respective OGCM represent the effect of the (fully or partially) resolved mesoscale on tracer mixing rates.
Following this reasoning, one would expect the difference in the Lagrangian eddy diffusivity estimates derived from a high-resolution OGCM experiment (e.g., SIMeddy-5d) and a coarse-resolution OGCM experiment (e.g., SIMnoeddy-5d) to be comparable in magnitude to the respective difference in the employed diffusivity coefficients. Both should be representative for the diffusive effect of mesoscale variability at horizontal scales between 0.18 and 0.58. This is, however, not the case: the difference in the Lagrangian eddy diffusivities exhibit a pronounced spatial variability at magnitudes O(10 3 ) m 2 s 21 . Because A ht does not include spatial variability apart from its adaptation to the changing grid sizes, the difference in A ht also does not feature spatial variability. Even more strikingly, the difference in A ht only reaches magnitudes O(10 2 ) m 2 s 21 .
This indicates that either the model coefficients are not chosen optimally, or Eulerian diffusivity coefficients and Lagrangian eddy diffusivity estimates are two measures that cannot be easily compared.
One could argue that Lagrangian eddy diffusivity estimates derived from drifter data quantify the rate of lateral dispersal in the surface mixed layer, whereas most modern ocean models, as those employed here, use isopycnal diffusivities for their subgrid-scale parameterizations. Yet, mesoscale eddies mix tracers along isopycnals and horizontally at the sea surface (Treguier et al. 1997;Abernathey et al. 2013), so the lateral eddy diffusivities estimated from GDP data and our simulated drifter trajectories can be treated as eddy diffusivity in the surface mixed layer. However, model diffusivities have to be tuned not only with respect to physical, but also to numerical considerations. Profound simplifications in the OGCM diffusion parameterization may further inhibit a direct relation between Eulerian diffusivity coefficients and Lagrangian eddy diffusivity estimates. Also to consider, as already stated by Rypina et al. (2016), Lagrangian diffusivity estimates constitute nonlocal measures, while the diffusivity coefficient for the subgrid-scale parameterization should represent the local effect of eddies. Finally, keeping in mind the challenges of Lagrangian eddy diffusivity estimation, such as removing the mean flow effect and capturing truly asymptotic regimes, Lagrangian eddy diffusivity estimates may not exclusively quantify the impact of mesoscale eddies but may be additionally influenced by methodological choices. Therefore, tuning coarse-resolution OGCM Eulerian eddy diffusivity coefficients based on Lagrangian eddy diffusivities estimated from high-resolution models or observations remains an outstanding issue that needs to be further addressed in more dedicated studies.

Conclusions
In this study, we assessed spatially variable nearsurface lateral (mesoscale) eddy diffusivity estimates obtained from both simulated Lagrangian trajectories and real drifter data for the greater Agulhas system following the approach of Zhurbas et al. (2014). Answering the three questions outlined in section 1, we showed that in this region, 1) Using 5-day mean velocity fields from the eddyresolving ocean model configuration INALT01, asymptotic diffusive regimes could be identified for dynamically different sites, some of which exhibit strong suppression of eddy mixing by mean flow.
2) The Lagrangian model-based eddy diffusivity estimates agree in pattern and magnitude with the observation-based estimates from Zhurbas et al. (2014) and are also comparable to the tracer-based eddy diffusivity estimates provided by Abernathey and Marshall (2013). 3a) Using monthly mean velocity output decreases the EKE but decreases eddy diffusivity estimates disproportionally less, supporting the idea that larger and more persistent eddy features are the dominant factors determining eddy diffusivities (cf. Wolfram et al. 2015). 3b) Using 5-day mean velocity fields from the noneddying ocean model configuration ORCA05 greatly reduces Lagrangian lateral eddy diffusivities; if a GM parameterization is employed, diffusivities further reduce to values close to zero.
These results suggest that when employing a stochastic Lagrangian parameterization to derive more diffusive trajectories from coarse-resolution model output-particularly for large-scale applications-one should consider the spatial variability of the diffusivity parameter, which does not necessarily scale with variability in EKE, as well as possible subgrid-scale parameterizations, such as GM. tory simulations were performed at High Performance Computing Centers in Stuttgart (HLRS), Hannover (HLRN) and at the Christian-Albrechts-Universität zu Kiel. OGCM code and data are available upon request from the corresponding author. For reproducibility of the main result figures in this study, respective data and scripts are available at http://data.geomar.de. The project received funding by the Cluster of Excellence 80 ''The Future Ocean'' within the framework of the Excellence Initiative by the Deutsche Forschungsgemeinschaft (DFG) on behalf of the German federal and state governments (Siren Rühs, Grant CP1412), by the German Federal Ministry of Education and Research (BMBF) (Arne Biastoch, Grant 03F0750A of the SPACES-AGULHAS project), and by the Helmholtz Association and the GEOMAR Helmholtz Centre for Ocean Research Kiel (Jonathan V. Durgadoo, Grants IV014 and GH018). Victor Zhurbas was supported by the Russian Science Foundation (Grant 14-50-00095) and the Russian Foundation for Basic Research (Grant 18-05-00278). The authors further wish to acknowledge the DRAKKAR group for support in model development and two anonymous reviewers for their insightful comments that helped to improve this manuscript.

Details of the Eddy Diffusivity Estimation-Smoothing Choices and Deviation from the Original Zhurbas et al. (2014) Approach
We followed the approach by Zhurbas et al. (2014) in all fundamental steps summarized in Fig. A1, which were already introduced by Oh et al. (2000), including the smoothing approaches adjusted to the spatiotemporal coverage of drifter data. Because of the increase of the sampling error with increasing time lag (Davis 1991b), a time-dependent finite differencing was applied to obtain k disp jk (x, t), and the tensor components k jk (x, t) and s jk (x, t) were low-pass filtered with a cosine filter and time-dependent filter window length. The time increment dt was chosen to equal the filter window [t 2 0.2t, t 1 0.2t]. Further details can be found in Oh et al. (2000), Zhurbas and Oh (2004), and Zhurbas et al. (2014). Note that simulated drifter trajectories were processed identically to the drifter data for comparability reasons. For a pure model study with a very large number of FIG. A1. Lateral eddy diffusivity estimation following Oh et al. (2000). trajectories, less smoothing and nonoverlapping bins may be sufficient and preferred, as indicated by Fig. A2; the asymptotic diffusivity estimates inferred from SIMeddy-5d hardly change when removing the time-dependent filtering and differencing and also show only small changes when changing to nonoverlapping 28 3 28 bins.
We deviated from the original Zhurbas et al. (2014) approach by slightly altering the quality control of drifter data to allow for speeds until 3 m s 21 (instead of removing data with speeds .2 m s 21 ) because in the region of interest, velocities reaching (exceeding) 2 m s 21 were reported from observations (and are apparent in the modeled velocity fields). We were further able to increase the number of considered pseudotracks for large time lags by extending the area over which trajectories were sampled for a certain bin (also necessary due to the very high local current velocities and, thus, large dispersion due to mean currents). Finally, we fixed a slight inconsistency in the estimation of K max that previously only considered time lags until 10 days (occasionally leading to K inf . K max ).