# An estimate of the eddy diffusivity tensor from observed and simulated Lagrangian trajectories in the Benguela Upwelling System

Ria Oelerich, Birte Gülk, Julia Dräger-Dietel, Alexa Griesel

► To cite this version:

Ria Oelerich, Birte Gülk, Julia Dräger-Dietel, Alexa Griesel. An estimate of the eddy diffusivity tensor from observed and simulated Lagrangian trajectories in the Benguela Upwelling System. Ocean Science, 2025, 21 (2), pp.727-747. 10.5194/os-21-727-2025. hal-05075923

HAL Id: hal-05075923

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

Submitted on 21 May 2025

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

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

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

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

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

# An estimate of the eddy diffusivity tensor from observed and simulated Lagrangian trajectories in the Benguela Upwelling System

Ria Oelerich<fnref n="1" />, Birte Gülk<fnref n="2" />, Julia Dräger-Dietel<fnref n="3" />, and Alexa Griesel<fnref n="3" />

**Correspondence:** Ria Oelerich (riaoelerich@gmail.com)

Received: 6 September 2024 – Discussion started: 18 September 2024

Revised: 24 January 2025 – Accepted: 24 January 2025 – Published: 3 April 2025

**Abstract.** Lateral mixing of unresolved processes in ocean models is usually parameterized with a scalar diffusivity, although the observed mixing can be highly anisotropic. Estimating the full diffusivity tensor from Lagrangian dispersion observations is challenging because shear dispersion from background currents can prevent the diffusive limit from being reached. This study investigates the diffusivity tensor with Lagrangian single- and pair-particle statistics in the Benguela upwelling region, using (1) a set of Lagrangian trajectories derived from a recent drifter dataset with an hourly resolution and background currents from the Ocean Surface Current Analysis Real-time (OSCAR) and (2) another trajectory set from simulations using the 1/10° Parallel Ocean Program (POP) simulation. Theory predicts that pair-particle diffusivities, which are expected to be independent of background mean flows, are twice the single-particle diffusivities if the pair velocities are uncorrelated. In this study, it is found that, although pair-particle diffusivities are much less influenced by mean flow, they are generally significantly smaller than twice the single-particle diffusivities. Subtracting the mean flow reduces this discrepancy and improves convergence in both methods, although single-particle diffusivities remain higher. Velocity autocorrelations decay faster than pair correlations, with mean flow subtraction accelerating decorrelation, especially in the zonal direction. The pair correlation term in the diffusivity equation contributes significantly to the differences between single-particle and pair-particle diffusivities, explaining why pair-particle diffusivities are generally smaller, making them a less accurate esti-

mate in diffusive parameterizations. In both the POP simulation and the observations, convergence properties improve significantly after mean flow subtraction. Mean flow removal plays a critical role in achieving convergence in the components of the diffusivity tensor and in the major-axis component after diagonalization. The significant anisotropy in the diffusivity tensor is mainly explained by the anisotropy in the Lagrangian integral timescales, while the major-axis component of the velocity variance tensor is only about 1.2 times the minor-axis component. The motions that are not resolved by the OSCAR surface current product but are captured by the surface drifters contribute significantly to the diffusivities, accounting for 8% and 42% of the contributions to the zonal and meridional components, respectively, after mean flow subtraction. This study highlights the importance of including the full diffusivity tensor in the Benguela upwelling region in lateral mixing parameterizations.

## 1 Introduction

Upwelling systems, found in equatorial and coastal regions, are areas where cold, nutrient-rich water from the deep ocean rises to the surface. This phenomenon is primarily driven by subtropical basin-scale circulations that create a slow, Equator-ward gyre flow near the shore (Capet et al., 2008). All eastern boundary upwelling systems share specific features, such as wind-driven flows, alongshore currents, shallow thermoclines, and significant vertical and offshore nu-

<footnote n="1">¹Institute of Environmental Physics, University of Bremen, Bremen, Germany</footnote>

<footnote n="2">²LOCEAN, Sorbonne Université-CNRS-IRD-MNHN, Paris, France</footnote>

<footnote n="3">³Institut für Meereskunde, Universität Hamburg, Hamburg, Germany</footnote>
trient transport (Large and Danabasoglu, 2006). However, each system exhibits unique characteristics in terms of circulation, primary productivity, and phytoplankton biomass, influenced by their specific geographic settings (Berger and Wefer, 2002).

The focus of this study is the Benguela Upwelling System, characterized by a frontal system of strong sea surface temperature (SST) gradients. This frontal system is frequently interrupted by finger-like structures of lower SST, known as filaments, resulting from meso- and submesoscale dynamics (Hösen et al., 2016). These filaments are influenced by eddies from two distinct origins: warm-core eddies, known as Agulhas rings, produced by the interaction between the Agulhas Current and the Antarctic Circumpolar Current (Agulhas Retroflection) after passing the Agulhas Bank at the southern tip of Africa (Lutjeharms, 2006), and cold-core eddies generated from boundary current instabilities near the coast of Namibia.

The dynamics in coastal upwelling systems like the Benguela Upwelling System range from well-understood processes to complex meso- and submesoscale flow instabilities (Capet et al., 2008). Due to the intricacies of these processes, climate simulations frequently face significant challenges, resulting in notable SST biases, although an adequate resolution of atmospheric processes is important as well (Small et al., 2015). This underscores the crucial need for ongoing research to enhance the accuracy of climate models and better understand the complexities of upwelling systems.

In coupled climate models with a resolution of 0.5° (~50 km) or lower, meso- and submesoscale eddies at their smallest spatial scales are likely to remain unresolved (Richter, 2015). Eddies contribute significantly to the transport of momentum and tracers (e.g., salinity and temperature) due to the comparatively weak mean ocean currents in upwelling systems, and they play a crucial role in the oceanic transport of colder water from the coast to the interior. In most climate models, the effects of eddies on a tracer T are parameterized using the tracer conservation equation:

$$ \frac{\partial \overline{T}}{\partial t} + \overline{u} \cdot \nabla \overline{T} = \nabla \cdot \overline{u' T'} + \overline{Q}, \quad (1) $$

where the overline denotes an average over intervals in time and space, ′ denotes the deviation from that aforementioned average, and $\overline{Q}$ describes the sources and sinks of the tracer T. Using the Transformed Eulerian Mean framework, the eddy flux can be decomposed into components across and along isolines of the mean tracer, e.g., diffusive and skew-diffusive fluxes (e.g., Griesel et al., 2019). Using $\overline{u'T'} = -K\nabla\overline{T} + B \times \nabla\overline{T}$, we receive the following expression:

$$ \frac{\partial \overline{T}}{\partial t} + (\overline{u} + \nabla \times \mathbf{B}) \cdot \nabla \overline{T} = \nabla \cdot \mathbf{K} \nabla \overline{T} + \overline{Q}, \quad (2) $$

where B is the vector stream function defining the eddy-driven advection velocity $\nabla \times B$ depending on the tracer T,

and K is a 3 × 3 diffusivity tensor denoting irreversible mixing of the tracer T. This diffusivity tensor K contains both lateral and vertical mixing components and has to be estimated for parameterizing eddies in climate models. Thus, an accurate representation of eddies in climate models is crucial to improve the general understanding of involved oceanic processes and their impact on future climate predictions (e.g., Ernst et al., 2023; Huot et al., 2022; Hewitt et al., 2020). This study focuses on the lateral components of the diffusivity tensor.

Diffusivities can be quantified using both Eulerian and Lagrangian methods (Griesel et al., 2014, 2019; Klocker and Abernathey, 2014). The Lagrangian approach is based on the spreading of floats or tracers as they follow the flow (Taylor, 1953, 1921; Nakamura, 1996), whereas Eulerian diffusivities can be quantified relative to a fixed geographical location, e.g., from Eulerian eddy tracer fluxes (e.g., Eden, 2006; Griesel et al., 2014). The traditional way to estimate eddy diffusivities originates from the dispersion of Lagrangian particles (Taylor, 1921, 1953), which provides the means to test the applicability of the eddy-diffusion model, where particles are expected to spread in a diffusive manner only after some time period determined by the temporal and spatial scales of the largest eddies (e.g., LaCasce et al., 2008). There are two possible statistical methods to estimate diffusivities, referred to as single- and pair-particle statistics. Single-particle statistics measure the spreading of particles from their origin and depend strongly on the background mean flow. Pair-particle statistics measure how two particles spread apart and would not depend on the background mean flow if both particles experienced the same mean flow during their spreading (Dräger-Dietel et al., 2018; Sansòn et al., 2017). In a flow with homogeneous statistics, the pair-particle diffusivity is expected to be twice the single-particle diffusivity.

In climate models, lateral eddy diffusion is typically represented by a single scalar, $\kappa_i$, aligned along isopycnal surfaces below the mixed layer. This approach is suitable when eddy diffusion is isotropic in the horizontal or along-isopycnal directions. In the presence of strong background flows, eddy diffusion is suppressed in the cross-stream component (Ferrari and Nikurashin, 2010; Klocker et al., 2012b, a; Klocker and Abernathey, 2014; Griesel et al., 2014) and eddy mixing becomes anisotropic (Bachman et al., 2020). Shear dispersion increases the effective diffusivity in the along-stream component (Taylor, 1953; Oh et al., 2000). This shear dispersion in the real ocean results from large-scale background mean flow, but it also occurs due to the eddies themselves, as large velocity shears exist at the eddy rims. The prevailing premise is that it is the mesoscale eddy motions with timescales of weeks and spatial scales corresponding to the first Rossby radius of deformation that contribute to eddy mixing (e.g., Ni et al., 2020; Stammer, 1998). Mixing is likely to be most efficient at the edge of the eddies, where large temperature gradients and large horizontal velocity shears occur. However, we further hypothesize that sub-
mesoscale motions with spatial scales below the first internal Rossby radius of deformation and associated with timescales shorter than a few days may contribute substantially. Recently, Peng et al. (2015) and Peng et al. (2020) reported that symmetric and shear instability occurs at the southern boundary of a filament in the Benguela Upwelling System, the same one that is investigated in this study. Diffusion might become more isotropic once the effect of the shear dispersion associated with the local shear from the mesoscale eddies is subtracted.

Dispersion processes are particularly important in the Benguela Upwelling System, where filaments are elongated in an east-west direction by the action of the velocities at the eddy rims, leading to increased mixing in this direction. This effect should be captured in ocean models. However, a question remains as to whether mixing in the along-stream component can be adequately described with a diffusivity. Oh et al. (2000) illustrated how the along-stream component of the diffusivity tensor is influenced by mean flow shear and does not always converge to a constant value. Oh et al. (2000) further advised using the minor-axis component of the diffusivity tensor, which is oriented across the background mean flow, as the diffusivity estimate. The minor-axis component of the diffusivity tensor can be estimated from drifter data (Zhurbas et al., 2014; Rühs et al., 2018). Davis (1987) and Davis (1991) devised the underlying theory regarding the computation of the diffusivity tensor in principle in the presence of an inhomogeneous background mean flow where, instead of diagnosing the statistics from the absolute Lagrangian velocities and displacements, diffusivities are computed from residual velocities and displacements after the Eulerian mean has been subtracted.

While most observational studies have focused on the minor axis of the diffusivity tensor, only a few, such as Rypina et al. (2012) in the North Atlantic and Peng et al. (2015) in the Indian Ocean, have considered all components of the diffusivity tensor from observational drifter data. Griesel et al. (2010, 2014) and Chen et al. (2014, 2015) demonstrated that, in many cases, the along-stream component of the diffusivity tensor can also be determined once the Eulerian mean is subtracted at each float point. However, they used trajectories from eddying ocean models where the background mean flow is known at each grid point.

In this study, for the first time, all components of the diffusivity tensor are estimated from the Lagrangian single- and pair-particle statistics (Taylor, 1921; Davis, 1987, 1991; Lacasse et al., 2008) in the Benguela upwelling region. Additionally, for the first time, our study compares the effect of the mean flow on both single- and pair-particle diffusivities simultaneously. Furthermore, we specifically quantify the contribution not resolved by the altimeter product to the diffusivity. An observed set of drifters (Dräger-Dietel et al., 2018) and simulated trajectories from drifters deployed in a high-resolution state-of-the-art eddy-permitting global ocean model (the 1/10° Parallel Ocean Program, hereinafter re-

ferred to as the POP simulation) are analyzed in detail. Similarly to Rypina et al. (2012), the background flows from the Ocean Surface Current Analyses Real-time (OSCAR) (Bonjean and Lagerloef, 2002) are considered in the analysis of the surface drifter dataset. This study aims to address the following questions:

- What is the dependence of the components of the diffusivity tensor on time lag, and can a diffusive limit be reached for all components?

- How do the diffusivities depend on eddy-mean flow compositions considering both single- and pair-particle statistics, and can the pair-particle diffusivity emerge as an alternative to the single-particle diffusivity, as it should be independent of the mean flow, which does not resolve the motions induced by the eddies?

- What is the role of smaller-scale motions that are captured by the drifters but not by current altimeter products and that have timescales shorter than a few days for the diffusivities and anisotropy?

The study is organized as follows: Sect. 2 presents the datasets used for this study, including the observed and simulated drifter trajectories. Section 3 describes the theory and methods used to estimate all components of the diffusivity tensor and the mean flow decomposition. Section 4 examines the prevailing background mean flow and eddy kinetic energy in the Benguela Upwelling System. Section 5 evaluates all components of the diffusivity tensor for single- and pair-particle statistics and highlights the impact of the background mean flow decomposition on the eddy diffusivities. Autocorrelations and pair correlations are examined in detail, and the impact of motions not resolved by satellite altimetry and their contributions to the diffusivities are investigated. Section 7 focuses on the evaluation of anisotropy by decomposing the diffusivity tensor in major- and minor-axis components in the Benguela Upwelling System. Section 8 presents the main conclusions and suggestions for future work.

## 2 Data

The observational dataset includes 35 surface drifters. The drifters of the type SVP-I-XDGS from MetOcean, consist of a surface buoy and a subsurface drogue, mounted in a 16 in. (40.6 cm) hull and equipped with Iridium telemetry and an SST sensor (Dräger-Dietel et al., 2018). With the drogue at 15 m depth, they provide the ability to monitor currents from this depth with less influence from surface winds.

The drifters were deployed in the eastern South Atlantic in November–December 2016 by the R/V *Meteor* during cruise M132 from Walvis Bay (Namibia) to Cape Town (South Africa) (Dräger-Dietel et al., 2018). The drifters were deployed at the northern boundary (at 26.0° S, 12.58° E) and
the southern boundary (26.4° S, 12.0° E) of a filament identified by satellite altimetry (Fig. 1a). The drifters were released as triplets in a triangular shape with initial separations of 100–200 m. In total, there were four drifter release sites, with the first and third groups released at the southern boundary and the second and fourth groups released at the northern boundary of the filament. After approximately 150 d, the drifters were transferred to the Global Drifter Program (GDP), where the temporal resolution was increased from 30 to 60 min. For this study, the drifter dataset with the lower resolution, as in Elipot et al. (2016), is used for a period of 250 d.

The set of simulated drifter trajectories consists of 948 drifters released at 10 m depth in the 1/10° global Parallel Ocean Program (POP) simulation. The run was initialized from the ocean state by Maltrud et al. (2010), which was integrated for 120 years using the annually repeating normal year Coordinated Ocean-ice Reference Experiments (COREs) surface forcing (Large and Yeager, 2004). In 1983, the COREs forcing was switched to vary interannually (Large and Yeager, 2009) and the model was integrated for an additional 27 years (Chouksey et al., 2022). The simulated drifters were released on 2 November 1996 with a spacing of 1/5° in both latitude and longitude between 24 and 30° S and between 10 and 18° E, covering a filament in the Benguela Upwelling System (Fig. 1b). The drifters are advected online for each model time step by the three-dimensional velocity field with a fourth-order Runge-Kutta scheme (Griesel et al., 2010). We note that the particles do not move vertically in our analysis period of 250 d; rather, they remain in the 5–15 m depth range.

Moreover, in the analysis of the observational drifter dataset, the OSCAR product is used to investigate the surface conditions with respect to the background mean flow and the eddy kinetic energy. The OSCAR product is provided by Earth Space Research (ESR, 2009, dataset accessed in November 2020). The dataset is on a 1/3° grid with a 5 d temporal resolution. The OSCAR product offers near-surface ocean current estimates that are derived using quasi-linear and steady-flow momentum equations. The spatial velocities in this dataset are directly estimated from sea surface height (SSH), surface winds, and sea surface temperature (SST) measurements taken by various altimetric satellites and in situ instruments.

## 3 Theory and methods

This section outlines the theoretical framework and methodologies used to investigate all components of the diffusivity tensor through both single- and pair-particle statistics, focusing on diffusive behavior, mean flow dependence, and smaller-scale motions.

### 3.1 Single-particle diffusivities

In a one-dimensional, stationary flow with zero mean, the mean square displacement $r^2$ of particles satisfies the following:

$$ \frac{d\langle r^2(t) \rangle}{dt} = \left\langle 2 \frac{dr(t)}{dt} r(t) \right\rangle \quad (3) $$

$$ = \left\langle 2u(t) \int_{0}^{t} u(t')dt' \right\rangle \quad (4) $$

$$ = 2 \int_{0}^{t} \langle u(t)u(t') \rangle dt' \quad (5) $$

$$ = 2 \int_{0}^{t} R(t')dt', \quad (6) $$

where $r(t)$ is the displacement of the particles at time $t$ from their position at time $t = 0$; angle brackets ($\langle \rangle$) denote an average over all particles considered; and $R(t)$ is the Lagrangian velocity autocovariance, which is characterized by its value at the origin $R(0) = \langle u(0)^2 \rangle$ and the Lagrangian integral timescale

$$ T_{\text{int}} = \frac{1}{R(0)} \int_{0}^{\infty} R(t')dt'. \quad (7) $$

For times much longer than the integral timescale, it follows that

$$ \langle r^2(t) \rangle = 2Kt \quad \text{with} \quad K = \langle u^2 \rangle T_{\text{int}}; \quad (8) $$

i.e., in analogy to molecular diffusion, the mean square displacement grows linearly with time with constant turbulent diffusivity $K$ (Taylor, 1921; LaCasce et al., 2008). In reality, with finite times and spatially varying background flows, the integral may not necessarily converge, and the diffusivity and timescale do not converge to a constant.

In a two-dimensional flow field that includes a spatially varying background current, the diffusivity is a $2 \times 2$ symmetric tensor and the equations generalize to

$$ \kappa_{ij} = \frac{1}{2} \frac{d}{dt} \langle r'_{i}(t) r'_{j}(t) \rangle \quad (9) $$

$$ = \int_{0}^{t} \langle u'_{i}(t) u'_{j}(t') \rangle dt' \quad (10) $$

$$ = \frac{\langle r'_{i}(t) r'_{j}(t) \rangle}{2t}, \quad (11) $$

with indices $i$ and $j$ both run over the zonal (1) and meridional (2) velocity and displacement directions. If the velocity
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The drifter trajectories are shown for (a) surface drifter observations and (b) particles advected with the velocities from the POP simulation from the deployment date (red dots) to 120 d after deployment (orange dots). The background colors display the snapshots of SST for the observations and the POP simulation. For the observations, the SST data are shown for the day of 28 November 2016 and represent measurements from the Moderate Resolution Imaging Spectroradiometer (MODIS) provided by the Earth Observing System (EOS) (dataset accessed in November 2020). For the POP simulation, the SST is shown for the day of 28 November 1996.</figcaption></figure>

statistics are homogeneous and stationary, Eqs. (9)–(11) lead to the same diffusivity – for real oceanic flows, the methods can lead to different results, as the averaging ⟨⟩ is applied at different stages (LaCasce et al., 2014). In this study, the integral of the velocity autocovariance (Eq. 10) is used to calculate the single-particle diffusivities.

## 3.2 Flow decompositions

Following Davis (1987, 1991), in Eqs. (9)–(11), $u'(t)$ represents the drifter velocities at time $t$ with the spatially varying Eulerian background flow removed and $r'(t)$ denotes the displacements with displacement due to the background flow removed. In that sense, the diffusivity introduced by Davis (1987) is a mixed Eulerian–Lagrangian quantity that can directly be related to the eddy tracer fluxes in a diffusive parameterization. Hence, the diffusivity reflects the turbulent mixing due to the residual components of the flow field and depends on the flow decomposition. In this study, two different flow decompositions are tested to estimate the contribution of different flow components to the magnitude and orientation of the diffusivity tensor:

$$u'(t) = u(t) - \overline{U}(x, y)|_{r(t)}, \quad (12)$$

$$u''(t) = u(t) - U(x, y, t_E)|_{r(t)}.$$

Here, $x$ is the zonal direction, $y$ is the meridional direction, the subscript $|r(\tau)$ indicates that the Eulerian velocities $U$ are evaluated at the float positions and times, and $t_E$ is the time at which the Eulerian currents are available.

In the first decomposition, the temporal mean but spatially varying Eulerian background flow $\overline{U}$ is interpolated to each drifter location and subsequently removed from the drifter

velocities. In the observational drifter dataset, $\overline{U}$ is the temporal mean velocity of the 1/3° OSCAR surface currents interpolated to the drifter locations, whereas in the POP simulation, $\overline{U}$ is the temporal mean of the 1/10° Eulerian currents interpolated to the numerical drifter locations (Table 1). Note that the current altimeter-derived products, such as the OSCAR product, could have effective resolutions in space and time closer to about 100 km and 30 d (Ballarotta et al., 2019), respectively, and should, in those instances, be regarded as a coarse-grained velocity.

In the second decomposition, which is only applied to the observed drifter data that capture motions with an hourly temporal resolution, the time-dependent background velocities $U(x, y, t)$ from the OSCAR surface current product are interpolated spatially and temporarily to match each drifter location and time and are subsequently removed from the drifter velocities. The $U(x, y, t)$ is the spatially varying OSCAR surface current available at 5 d intervals, and it allows for the estimation of the contribution of velocity variations to the diffusivity tensor that are not captured by OSCAR. These contributions can be associated with unresolved mesoscale motions as well as submesoscale motions that contain significant ageostrophic components (e.g., North et al., 2024). In addition, we can use the velocities from the OSCAR surface current product only to estimate diffusivities from $u$ and $u'$ (Table 1).

To visualize the displacements due to different flow components, similarly to Rypina et al. (2012) and Griesel et al. (2010), trajectories are calculated using the following residual velocities:
<table><caption>**Table 1.** Summary of different flow components used for diffusivity calculations.</caption><thead><tr><th>Velocity field</th><th>Description</th></tr></thead><tbody><tr><td><b>u</b><sub>OBS</sub></td><td>Total velocity at hourly drifter locations</td></tr><tr><td><b>u</b>'<sub>OBS</sub></td><td>Total velocity at hourly drifter locations minus 1/3° OSCAR yearly mean at drifter locations</td></tr><tr><td><b>u</b>''<sub>OBS</sub></td><td>Total velocity at hourly drifter locations minus 1/3° OSCAR 5 d means at drifter locations and times</td></tr><tr><td><b>u</b><sub>OSCAR</sub></td><td>1/3° OSCAR 5 d means at drifter locations and times</td></tr><tr><td><b>u</b>'<sub>OSCAR</sub></td><td>1/3° OSCAR 5 d means at drifter locations and times minus OSCAR yearly mean at drifter locations</td></tr><tr><td><b>U</b><sub>OSCAR</sub></td><td>Eulerian 1/3° OSCAR yearly mean velocity</td></tr><tr><td><b>u</b><sub>POP</sub></td><td>Total 1/10° POP velocity at daily numerical drifter locations</td></tr><tr><td><b>u</b>'<sub>POP</sub></td><td>Total 1/10° POP velocity at daily numerical drifter locations minus 1/10° POP yearly mean at numerical drifter locations</td></tr><tr><td><b>U</b><sub>POP</sub></td><td>Eulerian 1/10° POP yearly mean velocity</td></tr></tbody></table>

$$r'(t) = \int_{0}^{t} u'(t')dt', \quad r''(t) = \int_{0}^{t} u''(t')dt'. \qquad (13)$$

These trajectories are referred to as pseudo-trajectories, as (in contrast to the trajectories obtained from advection with the Eulerian flow fields) they are artificial trajectories obtained from integration with residual velocities.

## 3.3 Pair-particle diffusivities

As an alternative to the single-particle diffusivities, pair-particle diffusivity is calculated considering two particles (m and n), their separation (s), and their separation velocity ($v = u^m - u^n$):

$$\kappa^{\text{rel}}(t) = \frac{1}{2} \frac{d\langle s^2(t) \rangle}{dt} \qquad (14)$$

$$= \langle s(t)v(t) \rangle = \int_{0}^{t} \langle v(t)v(t') \rangle dt' \qquad (15)$$

$$= 2 \int_{0}^{t} \langle u^{m}(t)u^{m}(t') \rangle dt' - 2 \int_{0}^{t} \langle u^{m}(t)u^{n}(t') \rangle dt' \qquad (16)$$

$$= 2\kappa - 2 \int_{0}^{t} \langle u^{m}(t)u^{n}(t') \rangle dt'. \qquad (17)$$

From Eqs. (16) and (17), it is seen that the pair-particle diffusivity $\kappa^{\text{rel}}$ is twice the single-particle diffusivity $\kappa$ if the individual velocities are uncorrelated. For turbulent diffusion, the average correlation is positive for pair distances smaller than the correlation length (Koszalka et al., 2009) and leads to a reduced $\kappa^{\text{rel}}(t)$, where the reduction depends on the average initial pair separation. The pair-particle diffusivity is expected to be less dependent on the background mean flow, as any influence from the same background flow experienced by both particles is automatically subtracted in their relative dispersion. However, if the particles move apart into regions

with different background flows, the relative diffusivity will include contributions from the background shear. As for the single-particle diffusivity, all tensor components of the pair-particle diffusivity are considered using the zonal and meridional particle separations $s_i = x^m - x^n$, $s_j = y^m - y^n$ and the zonal and meridional separation velocity $v_i = u_i^m - u_i^n$, $v_j = v_j^m - v_j^n$:

$$\kappa_{ij}^{\text{rel}}(t) = \frac{1}{2} \frac{d\langle s_i(t)s_j(t) \rangle}{dt} \Big|_{t \gg T_{\text{int}}} = \frac{\langle s_i(t)s_j(t) \rangle}{2t} \qquad (18)$$

Furthermore, pair-particle diffusivities are calculated with the background flows removed, similar to the single-particle diagnostics. However, for the pair-particle diagnostics, $s'(t)$ is derived from the pseudo-trajectories (Eq. 13).

## 3.4 Anisotropy

The diffusivity tensor **K**, which is symmetric by construction, can be diagonalized $\mathbf{K} = R(\alpha) \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} R(\alpha)^{\text{T}}$, where $R(\alpha) = \begin{pmatrix} \cos(\alpha) & -\sin(\alpha) \\ \sin(\alpha) & \cos(\alpha) \end{pmatrix}$ is a rotation matrix characterizing the orientation of the principal axes of the diffusivity components. The eigenvalues of the symmetric Lagrangian eddy diffusivity tensor are given by

$$\lambda_{1,2} = \frac{1}{2} \left( \kappa_{xx} + \kappa_{yy} \pm \sqrt{(\kappa_{xx} - \kappa_{yy})^2 + 4\kappa_{xy}^2} \right), \qquad (19)$$

where the angle is given by

$$\tan 2\alpha = 2\kappa_{xy} / (\kappa_{xx} - \kappa_{yy}) \qquad (20)$$

(Haigh et al., 2020), and $\lambda_1$ is the diffusivity along the major axis of diffusion and corresponds to the eigenvector $e_1 = (\cos(\alpha), \sin(\alpha))$. As described by Oh et al. (2000), it is the along-stream component that is amplified by the shear dispersion of the background mean flow and does not converge to a constant value. However, when the background mean flow is subtracted from each drifter velocity, the shear dispersion effect should be minimized even in the along-stream direction.
## 3.5 Bootstrapping and diffusive limit

Finally, the important concept of this study is the investigation with regards to the saturated behavior of the diffusivities (the convergence to a constant value $\kappa^\infty$) and how this behavior is affected by the subtraction of the background mean flow. A convergence criterion is established, requiring the standard deviation to be less than 15% of the average eddy diffusivity from 100 to 250 d. This criterion ensures consistent comparison of convergence properties and diffusivities.

Given the limited sample size, an error estimation is derived through bootstrapping. This involves subsampling 35 bootstrap samples from the original datasets and conducting the analysis for single and pair particles with 100 repetitions. From these values, the mean and standard deviation from the mean are calculated. This approach yields a more accurate standard deviation than a simple calculation of the standard deviation from the original datasets.

## 4 Background mean flow, eddy kinetic energy, and pseudo-trajectories

The annual mean surface currents in the POP simulation (Fig. 2a, b) closely resemble those in the OSCAR surface current product, with both showing a predominantly north-westward Benguela Current near Africa's western coast between 6 and 18°E. In both datasets, the mean surface currents are primarily westward across the region. Generally, the zonal mean flow at the drifter locations is stronger in the POP simulation than in the OSCAR product, particularly near the coast, while the flows are similar in magnitude further offshore. Snapshots of surface currents from both datasets, taken 200 d after drifter release (5 d interpolation for OSCAR and daily for the POP simulation; Fig. 2c, d), reveal extensive eddy activity of similar magnitudes, filling the upwelling region. These eddies, with many originating from the Agulhas Retroflection site, appear in both datasets, although the POP simulation benefits from higher spatial and temporal resolution. The eddy kinetic energy (EKE) calculated from the annual means of both products (Fig. 2e, f) shows higher EKE magnitudes along the drifter trajectories in the POP simulation compared with the OSCAR product.

The POP simulation with numerical trajectories is only available for 1996–1997, while the observational dataset was obtained in 2016–2017. Regardless, as Fig. 2e and f illustrate, the horizontal distributions of kinetic energies are similar. There is visible interannual variability for the region (Fig. 3a) that is largely driven by the variability in the Agulhas rings. While there are discrepancies between OSCAR and POP, particularly around 2004, the area-averaged energy levels in 1996 in POP (red solid and dashed lines in Fig. 3a) and in 2016 for the OSCAR product (black solid and dashed lines in Fig. 3a) are similar, and so are the total kinetic energies of POP and OSCAR averaged along the drifter loca-

tions (red and cyan diamonds in Fig. 3a). However, the total kinetic energy of the surface drifter dataset at an hourly resolution averaged over the surface drifter locations in the observations is substantially higher (black diamond in Fig. 3a). This is largely due to inertial oscillations, visible as semidiurnal oscillations in the time series of the observational drifter velocities (black lines in Fig. 3b and c) as compared to the velocities obtained from POP and OSCAR that are similar in magnitude (red and cyan lines in Fig. 3b and c) but do not capture motions on these timescales.

The pseudo-trajectories (Fig. 4) reveal various motions occurring at different scales and due to different processes. Both the observations and the POP simulation exhibit a dominant westward mean flow. This is evident from the zonal extent of the pseudo-trajectories with the mean flow subtracted (Fig. 4a, d), which is reduced by up to 1000 km in both datasets. The total motion of particles (black trajectories in Fig. 4a and c) in both datasets consists of two components: the northwestward motion driven by the mean Benguela Current, which shows a stronger zonal flow in the POP simulation and a stronger northward component in the observations (gray trajectories in Fig. 4b and d), and the motion caused by eddies, which includes a net westward component as the eddies entrain filaments of cold water offshore (red trajectories in Fig. 4b and d) as also described in Dräger-Dietel et al. (2018).

## 5 The effect of the mean flow on eddy diffusivities

### 5.1 Diffusivities

As introduced in Sect. 3, the diffusivity reflects the turbulent mixing according to residual components of the flow field and, thus, depends on the flow decomposition. In this section, the eddy diffusivities calculated from the drifter velocities before and after subtraction of the annually averaged and spatially varying mean flow are compared for the observations and the POP simulation and for the single- and pair-particle approaches. This allows for investigation of the contribution of each mean flow decomposition on the diffusivity tensor and also enables a comparison of the mean flow contributions within each of the datasets and methods, for example, the comparison of single- and pair-particle-derived diffusivities. Note that, from this point onward, the diffusivities derived from pair-particle statistics are defined as $\kappa_{xx,xy,yx,yy}^{\text{rel}*} = 0.5\kappa_{xx,xy,yx,yy}^{\text{rel}}$ (Eq. 14) to allow for a direct comparison between pair- and single-particle diffusivities. Before the effect of the mean flow on single-particle diffusivities is discussed, we note that the eddy diffusivities from the POP simulation and observations are comparable to each other, as discussed in Sect. 4.

For single-particle diffusivities without mean flow subtraction (Fig. 5, red), the $\kappa_{xx}$ component is smaller in the POP simulation ($(50.98 \pm 10.28) \times 10^3 \, \text{m}^2 \, \text{s}^{-1}$) compared
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The annual means of surface current speeds are calculated for (a) observations (OSCAR product, 2016–2017) and (b) the POP simulation (1996–1997, 15 m depth). The annual mean current speeds are superimposed with velocity vectors (black arrows) for each respective dataset. The reference vector is given in red. The 5 d interpolation (OSCAR) and daily (POP simulation) snapshots of current speeds (c, d) 200 d after the drifter release are presented to emphasize the strong daily variations in current speeds. The eddy kinetic energy (e, f), calculated from the annual means, highlights the effects of eddies on the background mean flow after the drifter deployments. Blue markers are as in Fig. 1.</figcaption></figure>

with the observations ((60.14 ± 10.15) × 10³ m² s⁻¹). Similarly, the κyy component is also smaller in the POP simulation ((1.56 ± 0.39) × 10³ m² s⁻¹) compared with the observations ((4.59 ± 0.55) × 10³ m² s⁻¹). In contrast, the off-diagonal component κxy = κyx shows larger diffusivities in the POP simulation ((3.84 ± 1.54) × 10³ m² s⁻¹) than in the observations, with the latter exhibiting negative eddy diffusivities ((−3.80 ± 4.95) × 10³ m² s⁻¹).

The mean flow subtraction from absolute drifter velocities (blue in Fig. 5) affects the component of the diffusivity tensor that is most strongly impacted by the mean flow. κxx represents the zonal component of the diffusivity tensor and demonstrates the strongest reduction in diffusivities after subtraction of the mean flow across both datasets. This is to be expected, as the mean flow in the Benguela Upwelling

System is predominantly zonal, as discussed in Sect. 4. The reductions in κxx are about 55% for the observations and 86% for the POP simulation. Thus, the zonal mean flow contribution is largest for the diffusivities of the POP simulation. After the mean flow subtraction, κxx fulfills the convergence criterion (Sect. 3) for the POP simulation and the observations.

The κxy components are reduced by 55% in the POP simulation, while they show an increase of 342% in the observations (from negative diffusivities to positive diffusivities). It follows that largest reductions in κxy are again found for the POP simulation, but the reductions are smaller (with respect to the percentage) compared with the κxx component for the POP simulation. Additionally, with respect to the κxy component, the POP simulation does not fulfill the conver-
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The time series of annual mean kinetic energies for POP (red) compared with the OSCAR product (black) averaged over the Benguela upwelling region (as shown in Fig. 1) are shown in panel (a), where the total kinetic energies (solid lines) and the annual mean background flows (dashed lines) are displayed. Additionally, the total kinetic energies averaged over all float locations for POP, the observations, and the OSCAR velocities interpolated to drifter locations (red, black, and cyan diamonds, respectively) are shown. The time series of zonal and meridional velocities averaged over all drifters are displayed in panels (b) and (c) for POP (red), the observations (black), and OSCAR velocities interpolated to drifter locations (cyan).</figcaption></figure>

gence criterion and, thus, does not show saturated behavior after the subtraction of the mean flow. In contrast, $\kappa_{xy}$ shows convergence towards a constant value after the subtraction of the mean flow.

Finally, for the meridional component $\kappa_{yy}$, the diffusivities show an increase after the subtraction of the mean flow, where the diffusivities of the observations increase by 36% and the diffusivities for the POP simulation increase by 57%. This demonstrates that the subtraction of the mean flow does not necessarily lead to a reduction in diffusivities for all tensor components as the meridional component $\kappa_{yy}$ increases for both datasets. After the subtraction of the mean flow, the convergence criterion is fulfilled for both datasets and saturates towards a constant value.

For the pair-particle diffusivities derived from the relative drifter velocities (Fig. 6), the $\kappa_{xy}^{\text{rel}*}$ component for the observations ((3.38 ± 1.83) × 10³ m² s⁻¹) is about 55% larger than for the POP simulation ((1.41 ± 0.26) × 10³ m² s⁻¹), but it still remains the component with the largest eddy diffusivities. The $\kappa_{yy}^{\text{rel}*}$ component also demonstrates almost 78% smaller eddy diffusivities for the POP simulation ((0.92 ± 0.11) × 10³ m² s⁻¹) compared with the observations ((4.31 ± 0.47) × 10³ m² s⁻¹). The $\kappa_{xy}^{\text{rel}*}$ components for the POP simulation ((0.73 ± 0.05) × 10³ m² s⁻¹) are also about 70% smaller than the components from the observations ((2.5 ± 0.91) × 10³ m² s⁻¹).

In contrast to single-particle diffusivities, the pair-particle diffusivities show significantly smaller effects on eddy diffusivities after the subtraction of the mean flow. This is to be expected, as the relative drifter velocities are, in theory, less affected by the mean flow because they are derived from the velocity difference in the pairs, from which a common mean flow is then automatically subtracted (Sect. 3). $\kappa_{xx}^{\text{rel}*}$ shows that, also for pair-particle statistics, the subtraction of the mean flow does not necessarily lead to a reduction in eddy diffusivities.

The mean flow subtraction demonstrates eddy diffusivities with an increase of 30% for the observations and an increase of 35% for the POP simulation. This indicates that, in contrast to single-particle diffusivities, which mostly decrease after the mean flow subtraction, $\kappa_{xx}^{\text{rel}*}$ increases for the observations and the POP simulation. After the mean flow subtraction, $\kappa_{xx}^{\text{rel}*}$ fulfills the convergence criterion for the POP simulation but not for the observations. For $\kappa_{xy}^{\text{rel}*}$, there is a similar increasing behavior for both datasets after the mean flow subtraction. The eddy diffusivities change with an increase of 19% for the observations and 36% for the POP simulation. The convergence criterion is fulfilled for the POP simulation but not for the observations. $\kappa_{yy}^{\text{rel}*}$ shows opposite signs and presents decreasing behavior of the diffusivities after the subtraction of the mean flow with respect to the other tensor components. The diffusivities after the mean flow subtraction for $\kappa_{yy}^{\text{rel}*}$ show a reduction of 17% for the observations and
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The pseudo-trajectories (Eq. 13) for the observations (a, b) and the POP simulation (c, d). Note that all trajectory components are shown using the Cartesian coordinate system. For the POP simulation, every fifth trajectory is shown to improve the visibility of underlying patterns. The original trajectories are marked in shades of blue (a, c), whereas the pseudo-trajectories from the integration with $\mathbf{u}'$ (Eq. 12) are marked in shades of red (a, c). The pseudo-trajectories from the integration with the OSCAR surface currents $\mathbf{u}_{\text{OSCAR}}$ and the POP simulation daily mean velocities $\mathbf{u}_{\text{POP}}$ are highlighted in shades of blue (b, d) and those with the annual mean velocities interpolated to the drifter locations $\mathbf{U}|_{r(t)}$ are shown in shades of red (b, d). All pseudo-trajectories have been constructed such that (0,0) is the origin. Note that the coloring of trajectories was chosen such that every trajectory has a slightly different color in the range of blue or red so that patterns of individual trajectories are highlighted.</figcaption></figure>

an increase of 36% for the POP simulation. The convergence criterion after the mean flow subtraction is fulfilled for both datasets. It is important to note that the impact of the mean flow on pair-particle diffusivities is smaller compared with single-particle diffusivities across all tensor components.

The comparison of single- and pair-particle diffusivities clearly highlights that, specifically for the zonal component, $\kappa_{xx}$ components are significantly larger for the single-particle approach than for the pair-particle approach. This is mostly the case before and after the subtraction of the mean flow. Before the mean flow subtraction, the $\kappa_{xx}$ components with the single-particle approach are significantly larger, by 94% for the observations and 97% for the POP simulation. After the mean flow subtraction, $\kappa_{xx}$ components (compared with $\kappa_{xx}^{\text{rel}*}$ components) are larger, by about 84% for the observations and 73% for the POP simulation. Although, $\kappa_{xx}$ and $\kappa_{xx}^{\text{rel}*}$ do not show the same overall behavior before or after the subtraction of the mean flow, it is of note that smaller differences between $\kappa_{xx}$ and $\kappa_{xx}^{\text{rel}*}$ improve comparability in all datasets. The $\kappa_{xy}$ components compared with the $\kappa_{xy}^{\text{rel}*}$ components are larger, by 166% for the observations and 189% for the POP

simulation, before the mean flow subtraction. After the mean flow subtraction, the $\kappa_{xy}$ components (compared with $\kappa_{xy}^{\text{rel}*}$) are 68% larger for the observations and 42% larger for the POP simulation. The meridional components $\kappa_{yy}$ compared with $\kappa_{yy}^{\text{rel}*}$ before the mean flow subtraction show that $\kappa_{yy}$ is 6% larger for the observation and 41% larger for the POP simulation. After the subtraction of the mean flow, the $\kappa_{yy}$ components are larger than the $\kappa_{yy}^{\text{rel}*}$ components across both datasets, with percentages of 43% for the observations and 48% for the POP simulation.

This study further reveals that the convergence criterion is predominantly satisfied for the pair-particle approach, where 8 out of 12 tensor components across all datasets before and after subtraction of the mean flow showed saturated behavior. For the single-particle approach, 6 out of 12 tensor components demonstrated saturated behavior, where the subtraction of the mean flow had a positive effect on the convergent behavior for 4 out of 12 tensor components. The slightly stronger convergent properties for the pair-particle approach are to be expected, as relative drifter velocities are less affected by the mean flow. However, the results indicate that
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Single-particle eddy diffusivities are shown without subtraction of the mean flow (red) and with subtraction of the annually averaged mean flow (blue), as shown in Fig. 2a and b. The eddy diffusivities are given for all components of the diffusivity tensor $\kappa$ for the observations (a, c, e) and the POP simulation (b, d, f) (in m$^2$s$^{-1}$). The vertical black line in each panel marks the 100th day after drifter deployment, after which the eddy diffusivities are temporarily averaged. The temporarily averaged eddy diffusivities are shown with their standard deviations below each panel in the respective colors with and without mean flow subtraction (blue and red). The shaded areas in each panel indicate the uncertainties in eddy diffusivity calculations from bootstrapping (35 subsampled drifters with 100 repetitions).</figcaption></figure>

the single-particle approach generally provides larger diffusivities before and after the mean flow removal, which might indicate that some components of the mean flow may still be affecting the data, not only across the zonal component but also across the other tensor components. In most cases, the diffusivities from the observations and the POP simulation demonstrate similar values, although the values from the

POP simulation are slightly smaller before and after the mean flow subtraction.

In summary, our results show that the impact of the mean flow is smaller on pair-particle diffusivities than on single-particle diffusivities, whereas the subtraction of the mean flow improved comparability between single- and pair-particle diffusivities for all tensor components and across dif-
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Pair-particle eddy diffusivities are shown without subtraction of the mean flow (red) and with subtraction of the annually averaged mean flow (blue), as shown in Fig. 2a and b. The eddy diffusivities are given for all components of the diffusivity tensor $\kappa$ for the observations (a, c, e) and the POP simulation (b, d, f) (in $m^2 s^{-1}$). The vertical black line and shaded areas mark properties of the plot, as in Fig. 5.</figcaption></figure>

ferent datasets. In some cases, the subtraction of the mean flow improved the convergent properties specifically for the single-particle approach. However, even after the subtraction of the mean flow, single-particle diffusivities tend to be larger than pair-particle diffusivities and are not equal (as one might expect). This indicates that there might be an effect of the mean flow varying on scales smaller than the separation distance of the pair (which should be notable especially for

longer time lags) or a contribution of the correlation term (Eq. 17).

## 5.2 Velocity correlation term analysis

To analyze the effect of the correlation term, the velocity autocorrelations (Fig. 7a, b) and the velocity pair correlations (Fig. 7c, d) are compared in detail. The correlations are normalized by their values at zero time lag. In particular, the
comparison focuses on the *xx* and *yy* correlations after the mean flow is subtracted. Furthermore, the integral timescales from the autocorrelations and pair correlations are calculated as a measure of the decorrelation timescales (Eq. 7). It should be noted that if significant negative lobes and oscillations appear in the velocity correlations after the first zero crossing, the integral timescale will be considerably smaller than the time required for the diffusivities to converge.

First, it is noted that the kinetic energy averaged along the drifter trajectories (i.e., the non-normalized velocity autocorrelation at zero time lag) is higher for the observations than for the POP simulation (cf. black and red diamonds in Fig. 3). The total zonal and meridional kinetic energy is approximately 2.5 and 3.2 times greater in the observations compared with the POP simulation. When the mean flows are subtracted, the total zonal and meridional kinetic energy are approximately 3.6 and 3.3 times larger, respectively. The increase in the ratio with mean flow subtraction, particularly for the zonal component, is consistent with the stronger mean flow observed in the POP simulation, especially in the zonal direction (Fig. 4), leaving a smaller residual in the POP simulation.

The tails of the normalized correlations (Fig. 7) reveal that both the velocity autocorrelations and pair correlations for the *xx* component do not approach zero when the mean flow is not subtracted. This behavior aligns with the observation that single-particle diffusivities continue to increase with time lag (Fig. 5a, b). When the mean flow is subtracted (red lines in Fig. 7a and b), the autocorrelations for the *xx* components first cross zero at 83 d for the observations and at 182 d for the POP simulation. The *yy* components (blue lines in Fig. 7a and b) decay to zero much faster, occurring after 10 d in the observations and after 12 d for the POP simulation, followed by a significant negative lobe within the first 40 d and subsequent oscillations around zero. The meridional components exhibit similar behavior to the autocorrelations with mean flow subtraction (not shown), with the primary difference being reduced values at zero time lag. The times to the first zero crossing are longer for the pair correlation terms. For the *xx* component, the first zero crossing occurs at 236 d for the observations and at 167 d for the POP simulation when the mean flow is subtracted. In contrast, the *yy* components reach their first zero crossing much earlier, at 10 d for the observations and at 67 d for the POP simulation.

The autocorrelation integral timescales (Fig. 7a, b) are slightly longer for the POP simulation compared with the observations and are larger in the zonal direction than in the meridional direction. In contrast, the pair correlation integral timescales (Fig. 7c, d) are significantly larger, ranging from 2 to 5 times those of the autocorrelations.

The diffusivity values (Fig. 7) indicate that the correlation term in Eq. (17) is substantial and plays a key role in explaining the differences between single- and pair-particle diffusivities. This term is predominantly positive, contribut-

ing to the generally smaller pair-particle diffusivities. For instance, the correlation term for the *xx* component, without mean flow subtraction, is $(51.43 \pm 9.65) \times 10^3 \text{ m}^2 \text{s}^{-1}$ for both the observations and the POP simulation, accounting for the over 90% higher single-particle diffusivities. When the mean flow is subtracted, this correlation contribution decreases significantly, leading to better convergence between single- and pair-particle diffusivities, although the pair-particle diffusivities themselves remain relatively unaffected by the mean flow. Even after mean flow subtraction, the contribution of the pair correlation term for the *xx* component remains $(3.75 \pm 1.27) \times 10^3 \text{ m}^2 \text{s}^{-1}$. Notably, there are instances where the correlation contribution is negative; for example, in the *xy* component of the observations without mean flow subtraction, there is a negative contribution of $(-2.38 \pm 8.25) \times 10^3 \text{ m}^2 \text{s}^{-1}$, which corresponds to cases in which pair-particle diffusivity exceeds single-particle diffusivity (not shown). It is important to note that these values are sensitive to the chosen averaging interval with large error bars.

In summary, the analysis of the correlation terms reveals that pair diffusivities are not a reliable measure of lateral mixing, primarily because pair velocities remain correlated, especially in the zonal direction, where the first zero crossings occur only after more than 150 d following mean flow subtraction. The velocity autocorrelations decay to zero more rapidly than the pair correlations after mean flow subtraction.

Given that autocorrelation integral timescales are similar in both the observations and the POP simulation, the generally larger single-particle diffusivities observed in the observations compared with the POP simulation can be attributed to the higher kinetic energies (see Eq. 7). Furthermore, oscillations in the velocity autocorrelation reflect particle motions with timescales of a few days associated with mesoscale eddies (Griesel et al., 2010). Additionally, inertial oscillations with periods of about half a day, which are only present in the observational data, appear as oscillations in Fig. 7a and c, although they are not identifiable in the figure (due to the small periods) but, rather, appear as thickened lines. As inertial oscillations are associated with circular motions, they are hypothesized to not contribute much to net mean square displacement and, thus, should not impact the diffusivities. However, the energy contained in the inertial oscillations may contribute to the larger diffusivities in the observations compared with POP, which is further discussed in Sects. 6 and 7.

## 6 Effect of motions not resolved by satellite data

This section examines the contribution of unresolved motion in the OSCAR product to diffusivities from the observations by subtracting OSCAR-derived velocities from the total observed velocities, thereby isolating the impact of these unresolved motions on diffusivity estimates. The correspond-
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Panels (a) and (b) show the velocity autocorrelations normalized by their value at zero time lag $\langle u(t)u(t') \rangle$ (black), $\langle u'(t)u'(t') \rangle$ (red), and $\langle v'(t)v'(t') \rangle$ (blue), where the angle brackets $\langle \cdot \rangle$ represent an average over all particles, as a function of time lag. Panels (c) and (d) display the corresponding normalized velocity pair correlations $\langle u_i^m(t)u_j^n(t') \rangle$, where the angle brackets $\langle \cdot \rangle$ represent an average over all particle pairs. Also given are the integral timescales $T_{u'u'}$ and $T_{v'v'}$, averaged over 100–250 d (Eq. 7), when the mean flow is subtracted. Panels (a) and (c) show the results for observations, whereas panels (b) and (d) show the results for the POP simulation. The temporarily averaged (100–250 d) eddy diffusivities associated with the pair correlation terms $\int_0^t \langle u_i^m(t)u_j^n(t') \rangle dr'$ (last term on the right-hand side of Eq. 17) are shown with their standard deviations below the lower panels for $\kappa_{xx}$ with and without mean flow subtraction (black and red) and for $\kappa_{yy}$ with mean flow subtraction (blue) (in $m^2 s^{-1}$).</figcaption></figure>

ing pseudo-trajectories (Fig. 8) display inertial oscillations as well as other motions that are not resolved by the OSCAR product. The net displacements are found to be about a quarter of the total displacements in the pseudo-trajectories (Fig. 4a).

The diffusivity tensor components for both single and pair particles (Fig. 9a–c) indicate contributions in the range of $(1.00–2.50) \times 10^3 \text{ m}^2 \text{ s}^{-1}$. For single-particle diffusivities, the xx and yy components are found to be similar within the error bars, while the xy component becomes negative and does not converge, continuously decreasing with time lag, which suggests the presence of substantial nonzero autocorrelations for this component. Significant differences are observed between single- and pair-particle diffusivities for the xx and xy components, highlighting persistent pair correlations, whereas the yy components remain consistent within the error bars.

Overall, the contribution to the xx component is approximately 8% of that with the mean flow subtracted, while

<figure><img src="image_11.png" /><figcaption>**Figure 8.** The pseudo-trajectories from the integration with the residual velocity $u''$ from the observations. All pseudo-trajectories have been constructed such that (0,0) is the origin.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Single-particle (black) and pair-particle (red) eddy diffusivities using *u''*<sub>OBS</sub> (Eq. 11). The eddy diffusivities are given for all components of the diffusivity tensor for the observations in panels (a)–(c). The vertical black line and shaded areas mark the properties of the panels, as in Fig. 5.</figcaption></figure>

the contribution to the yy component is around 42% of the
value with the mean flow subtracted. This indicates that dif-
fusivities calculated from trajectories that are advected with
altimeter-derived velocities in the Benguela upwelling region
are likely to be significantly underestimated, especially in the
zonal component.

7 Anisotropy

The diffusivity tensors are characterized in terms of their
major- and minor-axis components as well as the angle that
defines the orientation of the principal axes of the diffusivity
components (Eq. 19). In the following, the anisotropy of the
diffusivity tensor is measured by the ratio of the major- to
minor-axis components.

In both the POP simulation and the observations, the major-axis components do not converge but, rather, increase continuously, as expected (Fig. 10). However, convergence is observed when the mean flow is removed, with the satisfaction of the convergence criterion, although this is only the case for the diffusivity plateaus for large time lags. For the POP simulation, mean flow subtraction results in a decrease in anisotropy, with the ratio of major- to minor-axis diffusivities dropping from about 44 to 4. Averaged over time lags of 100–250 d, the major-axis component drastically decreases from $(50.77 \pm 10.24) \times 10^3 \text{ m}^2 \text{s}^{-1}$ to $(7.77 \pm 0.71) \times 10^3 \text{ m}^2 \text{s}^{-1}$, whereas the minor-axis component actually increases from $(1.18 \pm 0.57) \times 10^3 \text{ m}^2 \text{s}^{-1}$ to $(1.81 \pm 0.33) \times 10^3 \text{ m}^2 \text{s}^{-1}$, consistent with the increase in the yy component (Sect. 5). In contrast, the observations show an increase in anisotropy that becomes negative, due to the minor-axis component being predominantly negative for time lags greater than about 70 d, although the error bars from the bootstrapping samples are large. The absolute value of the minor-axis component, averaged over time lags from 100 to 250 d, is less than $(0.4 \pm 1.0) \times 10^3 \text{ m}^2 \text{s}^{-1}$. The negative minor-axis component would indicate a decrease in mean square displacement with time along the minor-axis component due to convergence of particles, although the error bars from the standard deviation are large.

The diffusivities can be written as the product of major-
and minor-axis kinetic energies and integral timescales, as
in Eq. (7) (Table 2). It is found that after mean flow sub-
traction, the ratio of major- to minor-axis kinetic energies is
1.24 and 1.25 for the observations and the POP simulation,
respectively. It is mainly the integral timescales that explain
the large anisotropy in diffusivities. After mean flow subtrac-
tion, the major-axis timescales are about 14 d for the observa-
tions and 12 d for the POP simulation, while the minor-axis
timescales are -0.21 and 3 d for the observations and the
POP simulation, respectively. It should be noted that, while
the small anisotropy in the kinetic energy components is sim-
ilar in the POP simulation and the observations and the major
integral timescales are similar, the kinetic energies are larger
in the observations (Table 2), explaining the larger major-axis
diffusivities compared with the POP simulation. The R(0)
values are significantly larger in the observations compared
with the POP simulation, even after subtraction of the OS-
CAR velocities (third row in Table 2). This implies that the
kinetic energy associated with inertial oscillations might con-
tribute to the diffusivities in Sect. 6, even though the integral
timescales for u'' are small. Using the pseudo-trajectories
<figure><img src="image_13.png" /><figcaption>**Figure 10.** (a) Major- and minor-axis diffusivities (Eq. 19) for the POP simulation and the observations without mean flow subtraction. (b) Major- and minor-axis diffusivities for the POP simulation and the observations with mean flow subtraction. (c) Angle (Eq. 20) for the observations and the POP simulation with and without mean flow subtraction. (d) Major- and minor-axis diffusivities for subtraction with daily mean for the observations. Shown are the means of the diagonalization of the bootstrapping samples and their standard deviations (shading). The numbers in black, red, and blue are the ratios of the major- to minor-axis diffusivities, where the diffusivities were averaged over time lags of 100–250 d.</figcaption></figure>

derived from the contributions $u''$ not captured by OSCAR, the analysis shows that convergence is not achieved and the major-axis diffusivity continuously increases with time lag, while the minor-axis diffusivity decreases progressively and becomes negative after time lags exceeding approximately 220 d. The anisotropy remains significant, with a ratio of major- to minor-axis diffusivities of around 15.

The results show that the angle of the diffusivity ellipse is close to approximately 30° in the observations when the mean flow is subtracted, while it converges to nearly zero without mean flow subtraction (Fig. 10c). Note that the angle is less steep for the POP simulation compared with the observations after mean flow subtraction. When the daily mean is subtracted, the average angle is -30°, although with large error bars (not shown). The angle is expected to be correlated with the mean flow direction, as the maximum diffusion is anticipated to be oriented in the direction of the mean flow. The mean flow has a stronger meridional component in the observations compared with the POP simulation (Fig. 2a, b), especially for longitudes east of 6° E. Further offshore, the mean flow is zonally oriented and perpendicular to the coast for the observations as well as the POP simulation.

In summary, this study highlights that the minor-axis diffusivity can become negative and that anisotropy occurs for all flow components. The net westward displacements of the drifters and the anisotropy is not only caused by the action of the mean flow but also by the action of the mesoscale eddies that entrain the cold, upwelled water westward. Thus, the mixing is suppressed across both strong mean flow and mesoscale motions.

## 8 Discussion and conclusions

In this study, two sets of drifter trajectories were analyzed: one derived from observational drifter data and another from trajectories simulated using the POP (Parallel Ocean Program) simulation. The model simulation served as a test bed with about 30 times more drifter trajectories than available in the observations. One challenge in diffusivity estimates is that the influence of shear dispersion by background currents often inhibits the attainment of a diffusive limit, which is why many studies focus on the minor-axis or cross-stream components of the diffusivity tensor. Here, the first goal was to in-
<table><caption>**Table 2.** Velocity autocorrelation at zero time lag for major- and minor-axis components and integral Lagrangian timescales $T^{\text{int}}$ (in days) for different flow velocities (rows). The products of the timescales and velocity autocorrelations are the diffusivities, as given in Fig. 10, averaged over time lags of 100–250 d. A black asterisk indicates that the integral timescales and, therefore, that the diffusivities do not converge for that component.</caption><thead><tr><th>Velocity field</th><th><i>R</i>(0)<sub>major</sub> (cm<sup>2</sup> s<sup>-2</sup>)</th><th><i>T</i><sup>int</sup><sub>major</sub> (d)</th><th><i>R</i>(0)<sub>minor</sub> (cm<sup>2</sup> s<sup>-2</sup>)</th><th><i>T</i><sup>int</sup><sub>minor</sub> (d)</th></tr></thead><tbody><tr><td><b><i>u</i></b><sub>OBS</sub></td><td>299 ± 26</td><td>22.42 ± 3.60*</td><td>228 ± 12</td><td>1.52 ± 0.42</td></tr><tr><td><b><i>u'</i></b><sub>OBS</sub></td><td>277 ± 25</td><td>13.83 ± 1.40</td><td>223 ± 12</td><td>-0.21 ± 0.57</td></tr><tr><td><b><i>u''</i></b><sub>OBS</sub></td><td>165 ± 14</td><td>3.38 ± 0.53*</td><td>153 ± 7</td><td>0.24 ± 0.21*</td></tr><tr><td><b><i>u</i></b><sub>POP</sub></td><td>119 ± 10</td><td>49.20 ± 10.75*</td><td>72 ± 7</td><td>1.91 ± 0.95</td></tr><tr><td><b><i>u'</i></b><sub>POP</sub></td><td>78 ± 9</td><td>11.56 ± 1.71</td><td>68 ± 6</td><td>3.08 ± 0.63</td></tr></tbody></table>

investigate all components of the diffusivity tensor, study their dependence on mean flow subtraction using available background currents, and assess whether diffusive limits were reached. In the observations, the 1/3° OSCAR surface current product was used, while the 1/10° Eulerian background mean flow was readily available in the model. The single-particle diffusivities were also compared with pair-particle diffusivities, which are hypothesized to be less influenced by background mean flows, making them a potential alternative for single-particle diffusivity estimates where high-quality mean flows are not available.

The results show, as expected, that single-particle diffusivities are much more sensitive to mean flow subtraction than pair-particle diffusivities. In both the POP simulation and the observations, convergence properties significantly improve after mean flow subtraction. Mean flow removal plays a critical role in achieving convergence in the *xx* and *xy* tensor components for the single-particle diffusivities. Pair-particle diffusivities, on the other hand, demonstrate more consistent convergence and are less impacted by the mean flow, likely due to the inherent subtraction of common flow components between particle pairs.

However, pair-particle diffusivities are predicted by theory to be twice the single-particle diffusivity if the pair velocities are uncorrelated, and our study demonstrates that pair-particle diffusivities are considerably smaller than that (even after mean flow subtraction). We find that velocity autocorrelations decay more rapidly than the pair correlations, with the subtraction of mean flow leading to even faster decorrelation, particularly in the zonal direction. The integration of the pair correlation term in the equation for pair-particle diffusivity significantly contributes to the observed differences between single- and pair-particle diffusivities and explains why pair-particle diffusivities are generally smaller. As demonstrated by Davis (1987, 1991), the single-particle diffusivity after mean flow subtraction appears in diffusive parameterizations, so our pair-particle estimates would likely underestimate mixing in diffusive parameterizations. Pair-particle diffusivities can be an alternative to the single-particle estimates with mean flow subtraction if the deployment locations

are such that the pair velocities are uncorrelated from the beginning or quickly decorrelate.

The diffusivity can be written as the product of kinetic energy, i.e., the velocity autocorrelation at zero time lag, and the integral timescale, and we summarize the results for all flow components in Table 3. The results show that, after mean flow subtraction, the kinetic energies for the zonal and meridional directions become more similar to each other with ratios of 1.1–1.2 for the observations, the POP simulation, and the OSCAR product. However, the zonal integral timescales are 3.6 (observations), 2.6 (POP), and 1.5 (OSCAR) times larger than the meridional ones. The reduction in timescale for the meridional direction can be related to mixing suppression by a zonal mean flow, as discussed in Klocker et al. (2012a), Klocker and Abernathey (2014), and Griesel et al. (2015). The degree of mixing suppression can be based on mixing-length theory, although eddy propagation relative to the mean flow must also be taken into account. In the presence of a zonal background flow, the meridional diffusivity can then be written as follows:

$$ \kappa_{yy} = \alpha R(0) \frac{\gamma}{\gamma^2 + k^2(\bar{U} - c)^2}, \qquad (21) $$

where $\gamma$ is a typical Lagrangian decorrelation timescale and is equal to the growth rate of unstable waves in linear instability theory (Griesel et al., 2015), $k$ is related to a typical eddy size and can be regarded as the wavenumber of maximum growth in linear instability theory (although both might differ in the presence of an inverse energy cascade), $U$ is the zonal background mean flow, $c$ is a typical phase speed of the eddies, and $\alpha = 0.35$ as diagnosed by Klocker and Abernathey (2014). It is the difference in $U$ and $c$ that leads to the mixing suppression effect; otherwise, the diffusivity is equal to the velocity autocorrelation at zero lag times the Lagrangian decay scale $\gamma$. For the POP simulation, we use the eddy sizes and translational speeds as diagnosed from an eddy-tracking algorithm (Griesel et al., 2015), which amount to 90 km for the average eddy radius and $-0.033 \text{ m s}^{-1}$ for the zonal translational speed. We note that the translational speeds are close to, but not necessarily equal to, the phase speed in Eq. (21) (Griesel et al., 2015). With $\gamma$ as the decay scale (time to first zero crossing) for the meridional direc-
tion of about 1/(12d) and using our mean flow values aver-
aged over the drifter trajectories, the timescale is about 3d,
which is close to the Lagrangian integral timescale obtained
for POP for yy (Table 3). On the other hand, if we apply
Eq. (21) to the suppression by the mean meridional flow, us-
ing the region-averaged meridional translational eddy speeds
of 0.003 m s⁻¹ and γ = 1/(182d), we arrive at a timescale
of 14d, which is indeed close to Txx from Table 3. These
timescales are also consistent with those obtained by Klocker
and Abernathey (2014) and Rühs et al. (2018). We note that
mixing suppression in the zonal direction due to the mainly
zonally oriented topographic background potential vorticity
gradient might also be relevant here, a mechanisms that was
recently discussed by Sterl et al. (2024).

We find that the diffusivities in the observed drifter dataset
are larger than those inferred from the POP trajectories, par-
ticularly after mean flow subtraction. After mean flow sub-
traction, the integral timescales are slightly larger in the POP
simulation compared with the observations, although they
are similar overall (second and fifth rows of Table 3). How-
ever, the velocity autocorrelation at zero lag after mean flow
subtraction is 2–3 times larger in the observations than in the
POP simulation and, thus, explains the larger eddy diffusivi-
ties in the observations. This is partly due to the fact that the
background flow is larger and more highly resolved in the
POP simulation than in the OSCAR surface current prod-
uct, hence leaving a smaller EKE residual, and illustrates the
importance of using high-resolution background flow com-
ponents. However, the total kinetic energy (R(0) for uOBS
in Table 3) is also larger in the observations than in the POP
simulation, while the R(0) values, when considering the con-
tribution from the OSCAR velocities alone (last two rows in
Table 3), are more comparable to the POP simulation.

This relates to the research question that addresses the role
of smaller-scale motions which are detected by drifters but
not by current altimeter products in influencing diffusivities
and anisotropy. The study highlights the significant contribu-
tion of these unresolved motions. We found that these mo-
tions contribute 8% to the xx component and 42% to the
yy component of diffusivities after mean flow subtraction.
The observed differences between single- and pair-particle
diffusivities, particularly in the xx and xy components, em-
phasize the influence of ongoing pair correlations within the
unresolved small-scale dynamics. These small-scale motions
include inertial oscillations, which are hypothesized to not
contribute much to the net diffusion, as they may only lead
to oscillations in the velocity autocorrelation, which aver-
age out in the integral over the time lag. Indeed, as Table 3
shows, the integral timescales for u''_OBS are only about 1d
for the xx and yy components. However, the velocity au-
tocorrelation at zero lag is significant and is largely due to
the kinetic energy in the inertial motions. Further studies
with more trajectories are needed to clarify the role of in-
ertial motions for mixing. The underestimation of diffusiv-
ities when using altimeter-derived velocities indicates that

small-scale motions substantially influence diffusivity val-
ues and anisotropy, with the potential for large discrepan-
cies when these motions are not accounted for. Promising im-
provements are anticipated from the recent SWOT (Surface
Water and Ocean Topography) satellite mission, which will
provide a much higher resolution in geostrophic surface ve-
locity measurements. However, as highlighted by North et al.
(2024), velocities in the Benguela upwelling filament con-
tain significant ageostrophic components at scales smaller
than about 15 km, which will not be captured by geostrophic
velocities derived from sea surface height data. A promis-
ing approach was recently suggested by Zhang and Wolfe
(2024), who showed that eddy diffusivity can be inferred
from the eddy tracks without the need for drifter deployments
or altimetry-derived geostrophic velocity measurements.

Finally, the study finds significant anisotropy (measured
by the ratio of major- to minor-axis components) in the diffu-
sivities with all flow components. Similarly to the difference
between the xx and yy components of the diffusivity tensor,
the anisotropy of the diffusivity tensor is mainly explained
by the difference in major- and minor-axis timescales, as the
anisotropy in the velocity variance ellipses is small and sim-
ilar in both the observations and the POP simulation (Ta-
ble 2). The major-axis components of the diffusivity tensor
increase with time when no mean flow is subtracted, while
the major-axis components converge in the POP simulation
and the observations after the mean flow subtraction. The
minor-axis component in the observations becomes negative
for time lags greater than about 70 d after mean flow subtrac-
tion. This highlights the complex nature of diffusivity behav-
ior, where the Lagrangian diffusivity tensor is not necessarily
positive definite, as mean square displacement can decrease
over time. This might be due to the presence of convergent
velocities associated with the filament. Negative minor-axis
components were also recently found by Haigh et al. (2020),
who diagnosed Eulerian eddy diffusivities from tracer flux-
gradient relations in an eddying North Atlantic model sim-
ulation, and by Chen et al. (2014), who also diagnosed a
significant number of negative Lagrangian cross-stream eddy
diffusivities in the Kuroshio Extension. Negative Lagrangian
diffusivities are, in principle, consistent with mean square
displacements decreasing over time and net up-gradient eddy
tracer fluxes.

Subtracting the mean flow reduces the anisotropy in the
POP simulation but increases it in the observations due to
the near-zero (negative) minor-axis component (division by
a very small number). Mean flow subtraction aligns the dif-
fusivity ellipse more closely with the mean flow direction
(orientation northwestward). Hence, there is a difference be-
tween the yy component of the diffusivity tensor and the
minor-axis component. This finding is different from what
was found for the Indian Ocean from Global Drifter Program
data in Peng et al. (2015), where the mean flow was pre-
dominantly meridional; the zonal diffusivity, in most cases,
aligned with the minor-axis component; and the merid-
<table><caption>**Table 3.** Velocity autocorrelation at zero time lag and integral Lagrangian timescales $T^{\text{int}}$ (in days) (compare Eq. 7) for different tensor components (columns) and flow velocities (rows). The products of timescales and velocity autocorrelations are the diffusivities, as given in Figs. 5 and 9. A black asterisk indicates that the integral timescales and, therefore, that the diffusivities do not converge for that component.</caption><thead><tr><th>Velocity field</th><th><i>R</i>(0)<sub>xx</sub> (cm<sup>2</sup> s<sup>-2</sup>)</th><th><i>T</i><sup>int</sup><sub>xx</sub> (d)</th><th><i>R</i>(0)<sub>yy</sub> (cm<sup>2</sup> s<sup>-2</sup>)</th><th><i>T</i><sup>int</sup><sub>yy</sub> (d)</th><th><i>R</i>(0)<sub>xy</sub> (cm<sup>2</sup> s<sup>-2</sup>)</th><th><i>T</i><sup>int</sup><sub>xy</sub> (d)</th></tr></thead><tbody><tr><td><b><i>u</i></b><sub>OBS</sub></td><td>299 ± 19</td><td>23.24 ± 4.20*</td><td>228 ± 18</td><td>2.33 ± 0.34</td><td>-7 ± 70</td><td>658.32 ± 7066.65*</td></tr><tr><td><b><i>u</i></b><sub>OBS</sub></td><td>274 ± 18</td><td>11.46 ± 1.25</td><td>226 ± 19</td><td>3.20 ± 0.45</td><td>13 ± 6</td><td>81.26 ± 41.34</td></tr><tr><td><b><i>u</i></b><sub>OBS</sub></td><td>164 ± 11</td><td>1.61 ± 0.31</td><td>154 ± 9</td><td>1.97 ± 0.32</td><td>23 ± 3</td><td>-85.10 ± -130.75*</td></tr><tr><td><b><i>u</i></b><sub>POP</sub></td><td>119 ± 9</td><td>49.444 ± 10.64*</td><td>72 ± 8</td><td>2.51 ± 0.68</td><td>-3 ± 1</td><td>-213.33 ± 174.97</td></tr><tr><td><b><i>u</i></b><sub>POP</sub></td><td>77 ± 8</td><td>10.76 ± 1.26</td><td>69 ± 8</td><td>4.14 ± 0.50</td><td>3 ± 1</td><td>74.16 ± 44.11</td></tr><tr><td><b><i>u</i></b><sub>OSCAR</sub></td><td>106 ± 11</td><td>45.94 ± 9.08*</td><td>75 ± 0.0011</td><td>5.42 ± 0.90</td><td>-0.5 ± 4</td><td>-666.83 ± 4946.82*</td></tr><tr><td><b><i>u</i></b><sub>OSCAR</sub></td><td>85 ± 10</td><td>21.35 ± 3.01</td><td>74 ± 12</td><td>14.36 ± 2.99</td><td>8 ± 3</td><td>143.23 ± 59.62</td></tr></tbody></table>

ional diffusivity aligned with the major-axis component. The anisotropy illustrates the role of mesoscale eddies in contributing to westward displacements and anisotropy while suppressing mixing across strong mean and mesoscale flows.

Overall, the minor-axis components before mean flow subtraction are of the order of 3000 m² s⁻¹ for the observations and 1000 m² s⁻¹ for the POP simulation, which is similar to what has been found previously for this component from Lagrangian observations (Zhurbas et al., 2014; Rühs et al., 2018; Peng et al., 2015). However we find that the minor-axis component is also influenced by the mean flow subtraction, with a 1.5-fold increase for the POP simulation but a drastic decrease for the observations to a small, negative value. This may depend on the temporal averaging interval used for the mean flow, which is 1 year in our case. The values of major- and minor-axis components might change if a longer averaging interval is used.

Overall, our findings emphasize the necessity of using the full 2 × 2 diffusivity tensor to represent lateral mixing in climate models, rather than simplified isotropic representations. We speculate that the current warm SST bias in climate models in the eastern boundary upwelling regions may be reduced if the offshore, predominantly zonally oriented mixing of cold upwelled water by the action of eddies entraining the cold water of the filaments is represented properly.

*Data availability.* The observed drifter trajectories are part of the Global Drifter Program (GDP, https://www.aoml.noaa.gov/phod/gdp/index.php, last access: January 2024; https://doi.org/10.25921/x46c-3620, Elipot et al., 2022). The simulated drifters adapted within the POP simulation are available from https://doi.org/10.5281/zenodo.13710842 (Griesel, 2024). The data from the Ocean Surface Current Analyses Real-time (OSCAR) product are available from http://podaac.jpl.nasa.gov/dataset/OSCAR_L4 OC_ third-deg (ESR, 2009).

*Author contributions.* RO and AG developed the concept for the manuscript. RO and AG led the writing process. RO, BG, and AG processed and analyzed the observational and simulated drifter data. The field experiment design for drifter deployments was led by JDD

and AG. BG, AG, and JDD provided valuable input and guidance for the development of this study.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

*Acknowledgements.* The authors wish to thank the science party and the crew of cruise M132 on R/V *Meteor*, for their support, and Kerstin Jochumsen, as lead scientist. This paper is a contribution to project L3 (Meso- to submesoscale turbulence in the ocean) of the Collaborative Research Center TRR 181 “Energy Transfer in Atmosphere and Ocean”, funded by the German Research Foundation (DFG; grant agreement no. 274762653). We are grateful to the editor and the reviewers for their constructive comments, suggestions, and questions that improved our paper.

*Financial support.* This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. 274762653).

*Review statement.* This paper was edited by Julian Mak and reviewed by three anonymous referees.

## References

Bachman, S., Fox-Kemper, B., and Bryan, F. O.: A diagnosis of anisotropic eddy diffusion from a high-resolution global ocean model, J. Adv. Model. Earth Sy., 12, 1–14, https://doi.org/10.1029/2019MS001904, 2020.

Ballarotta, M., Ubelmann, C., Pujol, M.-I., Taburet, G., Fournier, F., Legeais, J.-F., Faugère, Y., Delepoulle, A., Chelton, D., Dibar-
boure, G., and Picot, N.: On the resolutions of ocean altimetry maps, Ocean Sci., 15, 1091–1109, https://doi.org/10.5194/os-15-1091-2019, 2019.

Berger, W. and Wefer, G.: On the reconstruction of upwelling history: Namibia upwelling in context, Mar. Geol., 180, 3–28, https://doi.org/10.1016/S0025-3227(01)00203-1, 2002.

Bonjean, F. and Lagerloef, G.: Diagnostic model and analysis of the surface currents in the tropical Pacific Ocean, J. Phys. Oceanogr., 32, 2938–2954, https://doi.org/10.1175/1520-0485(2002)032<2938:DMAAOT>2.0.CO;2, 2002.

Capet, X., McWilliams, F., Penven, J. C., and Marchesiello, P.: Eddies in eastern boundary subtropical upwelling systems, Geoph. Monog. Series, 177, 131–147, https://doi.org/10.1029/177GM10, 2008.

Chen, R., McClean, J. L., Gille, S. T., and Griesel, A.: Isopycnal eddy diffusivities and critical layers in the Kuroshio extension from an eddying ocean model, J. Phys. Oceanogr., 44, 2191–2211, https://doi.org/10.1175/JPO-D-13-0258.1, 2014.

Chen, R., Gille, S., McClean, J., Flierl, G., and Griesel, A.: A multi-wavenumber theory for eddy diffusivities and its application to the southeast Pacific (DIMES) region, J. Phys. Oceanogr., 45, 1877–1896, https://doi.org/10.1175/JPO-D-14-0229.1, 2015.

Chouksey, A., Griesel, A., Chouksey, M., and Eden, C.: Changes in global ocean circulation due to isopycnal diffusion, J. Phys. Oceanogr., 52, 2219–2235, https://doi.org/10.1175/JPO-D-21-0205.1, 2022.

Davis, R. E.: Modeling eddy transport of passive tracers, J. Mar. Res., 45, 635–666, 1987.

Davis, R. E.: Observing the general circulation with floats, Deep-Sea Res., 38, 531–571, https://doi.org/10.1016/S0198-0149(12)00293-9, 1991.

Dräger-Dietel, J., Jochumsen, K., Griesel, A., and Badin, G.: Relative dispersion of surface drifters in the Benguela upwelling region, J. Phys. Oceanogr., 48, 2325–2341, https://doi.org/10.1175/JPO-D-18-0027.1, 2018.

Eden, C.: Thickness diffusivity in the Southern Ocean, Geophys. Res. Lett., 33, L11606, https://doi.org/10.1029/2006GL026157, 2006.

Elipot, S., Lumpkin, R., Perez, R. C., Lilly, J. M., Early, J. J., and Sykulski, A. M.: A global surface drifter data set at hourly resolution, J. Geophys. Res.-Oceans, 121, 2937–2966, https://doi.org/10.1002/2016JC011716, 2016.

Ernst, P. A., Subrahmanyam, B., Trott, C. B., and Chaigneau, A.: Characteristics of submesoscale eddy structures within mesoscale eddies in the Gulf of Mexico from 1/48° ECCO estimates, Frontiers in Marine Science, 10, 1181676, https://doi.org/10.3389/fmars.2023.1181676, 2023.

ESR: OSCAR third degree resolution ocean surface currents, NASA Physical Oceanography Distributed Active Archive Center [data set], https://doi.org/10.5067/OSCAR-03D01, 2009.

Ferrari, R. and Nikurashin, M.: Suppression of eddy mixing across jets in the Southern Ocean, J. Phys. Oceanogr., 40, 1501–1519, https://doi.org/10.1175/2010JPO4278.1, 2010.

Griesel, A.: Particle Trajectories from 1/10 degree POP simulation in the Benguela upwelling region, Zenodo [data set], https://doi.org/10.5281/zenodo.13710842, 2024.

Elipot, S., Sykulski, A., Lumpkin, R., Centurioni, L., and Pazos, M.: Hourly location, current velocity, and temperature collected from Global Drifter Program drifters world-wide, NOAA

National Centers for Environmental Information [data set], https://doi.org/10.25921/x46c-3620, 2022.

Griesel, A., Gille, S. T., Sprintall, J., McClean, J. L., LaCasce, J. H., and Maltrud, M. E.: Isopycnal diffusivities in the Antarctic circumpolar current inferred from Lagrangian floats in an eddying model, J. Geophys. Res.-Oceans, 115, C06006, https://doi.org/10.1029/2009JC005821, 2010.

Griesel, A., McClean, M. J., Gille, S. T., Sprintall, J., and Eden, C.: Eulerian and Lagrangian isopycnal eddy diffusivities in the Southern Ocean from an eddying ocean model, J. Phys. Oceanogr., 44, 2191–2211, https://doi.org/10.1175/JPO-D-13-039.1, 2014.

Griesel, A., Eden, C., Koopmann, N., and Yulaeva, E.: Comparing isopycnal eddy diffusivities in the Southern Ocean with predictions from linear theory, Ocean Model., 94, 33–45, https://doi.org/10.1016/j.ocemod.2015.08.001, 2015.

Griesel, A., Dräger-Dietel, J., and Jochumsen, K.: Diagnosing and parameterizing the effects of oceanic eddies, in: Energy Transfers in Atmosphere and Ocean, edited by: Eden, C. and Iske, A., Mathematics of Planet Earth, Springer-Verlag, 193–224, https://doi.org/10.1007/978-3-030-05704-6_6, 2019.

Haigh, M., Sun, L., Shevchenko, I., and Berloff, P.: Tracer-based estimates of eddy-induced diffusivities, Deep-Sea Res. Pt. I, 160, 103264, https://doi.org/10.1016/j.dsr.2020.103264, 2020.

Hewitt, H. T., Roberts, M., Mathiot, P., Biastoch, A., Blockley, E., Chassignet, E. P., Fox-Kemper, B., Hyder, P., Marshall, D. P., Popova, E., Treguier, A.-M., Zanna, L., Yool, A., Yu, Y., Beadling, R., Bell, M., Kuhlbrodt, T., Arsouze, T., Bellucci, A., Castruccio, F., Gan, B., Putrasahan, D., Roberts, C. D., Van Roekel, L., and Zhang, Q.: Resolving and parameterising the ocean mesoscale in Earth System models, Current Climate Change Reports, 6, 137–152, https://doi.org/10.1007/s40641-020-00164-w, 2020.

Hösen, E., Möller, J., Jochumsen, K., and Quadfasel, D.: Scales and properties of cold filaments in the Benguela upwelling system off Lüderitz, J. Geophys. Res.-Oceans, 121, 1896–1913, https://doi.org/10.1002/2015JC011411, 2016.

Huot, P.-V., Kittel, C., Fichefet, T., Jourdain, N. C., and Fétweis, X.: Effects of ocean mesoscale eddies on atmosphere-sea ice-ocean interactions off Adélie Land, East Antarctica, Clim. Dynam., 59, 41–60, https://doi.org/10.1007/s00382-021-06115-x, 2022.

Klocker, A. and Abernathey, R.: Global patterns of mesoscale eddy properties and diffusivities, J. Phys. Oceanogr., 44, 1030–1046, https://doi.org/10.1175/JPO-D-13-0159.1, 2014.

Klocker, A., Ferrari, R., and LaCasce, J. H.: Estimating suppression of eddy mixing by mean flows, J. Phys. Oceanogr., 9, 1566–1576, https://doi.org/10.1175/JPO-D-11-0205.1, 2012a.

Klocker, A., Ferrari, R., LaCasce, J. H., and Merrifield, S. T.: Recombining float-based and tracer-based estimates of eddy diffusivities, J. Mar. Res., 70, 569–602, 2012b.

Koszalka, I., Lacasce, J., and Orvik, K.: Relative dispersion in the Nordic Seas, J. Mar. Res., 67, 411–433, 2009.

LaCasce, J. H., Ferrari, R., Marshall, J., Tulloch, R., Balwada, D., and Speer, K.: Statistics from Lagrangian observations, Prog. Oceanogr., 77, 1–29, https://doi.org/10.1016/j.pocean.2008.02.002, 2008.

LaCasce, J. H., Ferrari, R., Marshall, J., Tulloch, R., Balwada, D., and Speer, K.: Float derived isopycnal diffusivities
in the DIMES experiment, J. Phys. Oceanogr., 44, 764–780, https://doi.org/10.1175/JPO-D-13-0175.1, 2014.

Large, W. G. and Danabasoglu, G.: Attribution and impacts of upper-ocean biases in CCSM3, J. Climate, 19, 2325–2346, https://doi.org/10.1175/JCLI3740.1, 2006.

Large, W. G. and Yeager, S. G.: Diurnal to decadal global forcing for ocean and sea-ice models: the datasets and flux climatologies, NCAR Technical Note TN-460+STR, 8, National Center for Atmospheric Research, https://doi.org/10.5065/D6KK98Q6, 2004.

Large, W. G. and Yeager, S. G.: The global climatology of an interannually varying air-sea flux data set, Clim. Dynam., 33, 341–364, https://doi.org/10.1007/s00382-008-0441-3, 2009.

Lutjeharms, J.: The Agulhas Current, Springer Verlag Berlin-Heidelberg, https://doi.org/10.1007/3-540-37212-1, 2006.

Maltrud, M. E., Bryan, F., and Peacock, S.: Boundary impulse response functions in a century-long eddying global ocean simulation, Environ. Fluid Mech., 10, 275–295, https://doi.org/10.1007/s10652-009-9154-3, 2010.

Nakamura, N.: Two-Dimensional Mixing, Edge Formation, and Permeability Diagnosed in an Area Coordinate, J. Atmos. Sci., 53, 1524–1537, https://doi.org/10.1175/1520-0469(1996)053<1524:TDMEFA>2.0.CO;2, 1996.

Ni, Q., Zhai, X., Wang, G., and Marshall, D. P.: Random movement of mesoscale eddies in the global ocean, J. Phys. Oceanogr., 50, 2341–2357, https://doi.org/10.1175/JPO-D-19-0192.1, 2020.

North, R. P., Dräger-Dietel, J., and Griesel, A.: Characterization of physical properties of a coastal upwelling filament with evidence of enhanced submesoscale activity and transition from balanced to unbalanced motions in the Benguela upwelling region, Ocean Sci., 20, 103–121, https://doi.org/10.5194/os-20-103-2024, 2024.

Oh, S. I., Zhurbas, V., and Park, W. S.: Estimating horizontal diffusivity in the East Sea (Sea of Japan) and the northwest Pacific from satellite-tracked drifter data, J. Geophys. Res.-Oceans, 105, 6483–6492, https://doi.org/10.1029/2000JC900002, 2000.

Peng, J.-P., Holtermann, P., and Umlauf, L.: Frontal instability and energy dissipation in a submesoscale filament, J. Phys. Oceanogr., 50, 2017–2035, https://doi.org/10.1175/JPO-D-19-0270.1, 2020.

Peng, S., Qian, Y.-K., Lumpkin, R., Li, P., Wang, D., and Du, Y.: Characteristics of the near-surface currents in the Indian Ocean as deduced from satellite-tracked surface drifters. Part II: Lagrangian statistics, J. Phys. Oceanogr., 45, 459–477, https://doi.org/10.1175/JPO-D-14-0049.1, 2015.

Richter, I.: Climate model biases in the eastern tropical oceans: causes, impacts and ways forward, WIREs Clim. Change, 6, 345–358, https://doi.org/10.1002/wcc.338, 2015.

Rühs, S., Zhurbas, V., Koszalka, I. M., Durgadoo, J. V., and Bias-toch, A.: Eddy diffusivity estimates from lagrangian trajectories simulated with ocean models and surface drifter data – a case study for the greater agulhas system, J. Phys. Oceanogr., 48, 175–196, https://doi.org/10.1175/JPO-D-17-0048.1, 2018.

Rypina, I. I., Kamenkovich, I., Berloff, P., and Pratt, L. J.: Eddy-induced particle dispersion in the near-surface North Atlantic, J. Phys. Oceanogr., 42, 2206–2228, https://doi.org/10.1175/JPO-D-11-0191.1, 2012.

Sansòn, L. Z., Pèrez-Brunius, P., and Sheinbaum, J.: Surface relative dispersion in the southwestern Gulf of Mexico, J. Phys. Oceanogr., 47, 387–403, https://doi.org/10.1175/JPO-D-16-0105.1, 2017.

Small, J. R., Churchister, E., Hedstrom, K., Kaufmann, B., and Large, W.: The Benguela upwelling system: quantifying the sensitivity to resolution and coastal wind representation in a global climate model, J. Climate, 28, 9409–9432, https://doi.org/10.1175/JCLI-D-15-0192.1, 2015.

Stammer, D.: On eddy characteristics, eddy transports, and mean flow properties, J. Phys. Oceanogr., 28, 727–739, https://doi.org/10.1175/1520-0485(1998)028<0727:OECETA>2.0.CO;2, 1998.

Sterl, M. F., LaCasce, J. H., Groeskamp, S., Nummelin, A., Isachsen, P., and Baatsen, M.: Suppression of mesoscale eddy mixing by topographic PV gradients, J. Phys. Oceanogr., 54, 1089–1103, https://doi.org/10.1175/JPO-D-23-0142.1, 2024.

Taylor, G.: Diffusion by continuous movements, P. Lond. Math. Soc., 20, 196–212, https://doi.org/10.1112/plms/s2-20.1.196, 1921.

Taylor, G.: Dispersion of soluble matter in solvent flowing slowly through a tube, P. Roy. Soc. A-Math. Phy., 219, 186–203, https://doi.org/10.1098/rspa.1953.0139, 1953.

Zhang, W. and Wolfe, C.: Inferring tracer diffusivity from coherent mesoscale eddies, J. Adv. Model. Earth Sy., 16, 1–25, https://doi.org/10.1029/2023MS004004, 2024.

Zhurbas, V., Lyzhkov, D., and Kuzmina, N.: Drifter-derived estimates of lateral eddy diffusivity in the World Ocean with emphasis on the Indian Ocean and problems of parameterization, Deep-Sea Res. Pt. I, 83, 1–11, https://doi.org/10.1016/j.dsr.2013.09.001, 2014.