Inferring Tracer Diffusivity From Coherent Mesoscale Eddies

Abstract

Mixing along isopycnals plays an important role in the transport and uptake of oceanic tracers. Isopycnal mixing is commonly quantified by a tracer diffusivity. Previous studies have estimated the tracer diffusivity using the rate of dispersion of surface drifters, subsurface floats, or numerical particles advected by satellite-derived velocity fields. This study shows that the diffusivity can be more efficiently estimated from the dispersion of coherent mesoscale eddies. Coherent eddies are identified and tracked as the persistent sea surface height extrema in both a two-layer quasigeostrophic (QG) model and an idealized primitive equation (PE) model. The Lagrangian diffusivity is estimated using the tracks of these coherent eddies and compared to the diagnosed Eulerian diffusivity. It is found that the meridional coherent eddy diffusivity approaches a stable value within about 20–40 days in both models. In the QG model, the coherent eddy diffusivity is a good approximation to the upper-layer tracer diffusivity in a broad range of flow regimes, except for small values of bottom friction or planetary vorticity gradient, where the motions of same-sign eddies are correlated over long distances. In the PE model, the tracer diffusivity has a complicated vertical structure and the coherent eddy diffusivity is correlated with the tracer diffusivity at the e-folding depth of the energy-containing eddies where the intrinsic speed of the coherent eddies matches the rms eddy velocity. These results suggest that the oceanic tracer diffusivity at depth can be estimated from the movements of coherent mesoscale eddies, which are routinely tracked from satellite observations.

Introduction

Ocean mesoscale eddies play an important role in the transport and mixing of oceanic tracers, such as heat, carbon, and nutrients. Tracer mixing by mesoscale eddies impacts the large-scale ocean circulation (Hallberg & Gnanadesikan, 2006; Marshall & Radko, 2003, 2006; Wolfe & Cessi, 2010) and biogeochemical environment (Gnanadesikan et al., 2015; McGillicuddy Jr et al., 2003; Steinberg et al., 2019). Mesoscale eddies can contribute to tracer mixing by both stirring the background tracer distribution and by trapping tracers in their cores and moving them as the eddies drift (Frenger et al., 2015; W. Zhang et al., 2020). Here we focus on the former mechanism and refer to transport by stirring when mesoscale tracer transport or mixing is mentioned. Mesoscale tracer mixing is commonly parameterized using a combination of eddy-induced advection (Gent & McWilliams, 1990) and diffusion of tracers along isopycnals (Redi, 1982) in the coarse-resolution ocean component of climate models. Climate simulations are sensitive to the magnitude and distribution of the isopycnal tracer diffusivity (Pradal & Gnanadesikan, 2014; Sijp et al., 2006), which needs to be constrained by observational measurements.

Lagrangian methods have been used to estimate the tracer diffusivity in the ocean using surface drifters (Roach et al., 2018; Rühs et al., 2018; Rypina et al., 2012; Zhurbas et al., 2014; Zhurbas & Oh, 2003, 2004), subsurface floats (J. LaCasce et al., 2014; Balwada et al., 2016, 2021), and numerical particles advected by satellite derived flow fields (Klocker, Ferrari, LaCasce, & Merrifield, 2012; Rypina et al., 2012). According to Taylor (1922), for homogeneous and stationary turbulent flow, the eddy diffusivity can be estimated by the continuous movements of single Lagrangian particles, [Image Omitted. See PDF]where x_i(t) is the position of a particle found at x₀ at time t = 0, $${\left\langle \cdot \right\rangle }_{L}$$ indicates the Lagrangian mean, which is the average over the ensemble of particles, and $${\left\langle {\left(\mathbf{x}(t)-{\mathbf{x}}_{0}\right)}^{2}\right\rangle }_{L}$$ is the absolute dispersion of particles. Modified versions of Equation 1 have been developed by Davis (1987, 1991) to account for the inhomogeneity and anisotropy of mixing. Accurate estimates using these methods require averages over large numbers (order of hundreds) of drifters (Klocker, Ferrari, LaCasce, & Merrifield, 2012), but the spatial distribution of drifters and floats is generally sparse and many surface drifters are contaminated by wind effects (Lumpkin et al., 2013). Also, the estimates of diffusivity by these methods generally takes a long time (the order of months) to asymptote to the “true” value of diffusivity (Klocker, Ferrari, & LaCasce, 2012; Rypina et al., 2012). The long convergence timescale makes these estimates inefficient and allows errors to accumulate as the sampling error grows with time (Davis, 1991). Further, particles in inhomogeneous flow might move to a different region with a different mixing rate during the period of diffusivity calculation.

Mesoscale eddies are increasingly observed and studied as individual coherent structures that can be identified and tracked from satellite observations (Chelton et al., 2011; Dong et al., 2014; Z. Zhang et al., 2014). Coherent eddies are swirling structures that can move in the ocean over a potentially long distance (Chelton et al., 2011). Studies of vortex-dominated 2D turbulence have used Equation 1 to estimate a diffusivity from the movements of coherent eddies, with eddy displacements replacing particle displacements (Chong et al., 2020; Hansen et al., 1998; J. H. LaCasce, 2008b; Weiss et al., 1998). This is based on the observation that the movement of coherent eddies resembles that of particles (Weiss et al., 1998); that is, the motion is initially ballistic (absolute dispersion quadratic in time) and then transitions to diffusive (absolute dispersion linear in time). This evolution of the absolute dispersion is also a typical feature of Brownian motion (Chong et al., 2020). When the motion of coherent eddies becomes diffusive, a Lagrangian diffusivity can be estimated from their motion using Equation 1.

The movement of coherent ocean eddies has both systematic and chaotic features. Eddies are impacted by the β-effect, which causes them to drift systematically westward relative to the mean flow (Cushman-Roisin et al., 1990). The β-effect also leads to meridional “beta drift”: cyclonic eddies tend to propagate poleward and anticyclonic eddies tend to propagate equatorward (e.g., Holland, 1982; Nycander, 2001; R. B. Smith, 1993; Sutyrin et al., 1994; Korotaev, 1997). In addition to these systematic drifts, coherent eddies also move randomly due to eddy-eddy interactions (Ni et al., 2020; Samelson et al., 2014, 2016). The random movements of coherent eddies have been used to estimate a diffusivity by Ni et al. (2020). However, they interpreted this diffusivity to represent the spreading of eddy energy rather than the mixing of tracers.

W. Zhang et al. (2020) recently found that the Lagrangian diffusivity estimated from the dispersion of coherent eddies using Equation 1 can provide an accurate estimate of the Eulerian PV diffusivity in a two-layer QG model. W. Zhang et al. (2020) defined a coherent eddy as a “rotationaly coherent Lagrangian vortex” (RCLV) (Haller et al., 2016), which can trap particles (water parcels) inside them over a long time. It was found that particles trapped within coherent eddies have a negligible contribution to the total particle dispersion (consistent with Abernathey & Haller, 2018), but that the diffusive movement of coherent eddies themselves was representative of the diffusivity of the flow. If this finding also applies to mesoscale eddies in the ocean, then ocean tracer diffusivity can also be inferred from the dispersion of coherent mesoscale eddies. To support this application, two questions need to be addressed. First, the definition of a coherent eddy (i.e., RCLV) used in W. Zhang et al. (2020) is stricter than the commonly used definition of a coherent eddy, which is an Eulerian feature (e.g., a sea surface height or vorticity extremum) that can be identified and tracked over a long time (Chelton et al., 2007, 2011; Mason et al., 2014). Does the dispersion of Eulerian coherent eddies also provide an accurate estimate of tracer diffusivity? Second, W. Zhang et al. (2020) only qualitatively compared the diffusivity estimated from coherent eddies and the Eulerian PV diffusivity in three QG simulations with varying bottom friction. Can the results of W. Zhang et al. (2020) be generalized to broader QG regimes and more realistic 3D ocean circulations?

In this study, we compare the Lagrangian diffusivity estimated from the movement of coherent eddies to the tracer diffusivity in a two-layer QG model and a 3D primitive equation (PE) model. We find the coherent eddy diffusivity provides an accurate estimate of the upper-layer meridional tracer diffusivity in the QG model and is highly correlated with the tracer diffusivity at a depth determined by a nonlinearity parameter in the PE model. This depth is close to the e-folding vertical scale of the energy-containing eddies, which can be estimated from sea surface height and hydrography. These findings can be further used to infer and interpret the lateral tracer diffusivity in the ocean based on tracking the coherent mesoscale eddies.

This manuscript is structured as follows. In Section 2, we describe the configuration of the numerical models and analysis methods. In section3, the coherent eddy diffusivity is shown to accurately reproduce the Eulerian tracer diffusivity in the QG model except when either bottom friction or beta are small. We discuss the reason for the discrepancy for small beta and friction simulations in Section 4. The coherent eddy diffusivity is then compared to the tracer diffusivity in the PE model in Section 5. The conclusions are summarized in Section 6.

Approach

In this section, we introduce the key components of our approach: (a) the configuration of the QG and PE models that are used to examine the coherent eddy diffusivity; (b) the methods used to track coherent eddies and compute the coherent eddy diffusivity in both models; and (c) the methods to estimate the tracer diffusivity in the two models. The detailed approaches for estimating the coherent eddy and tracer diffusivities differ in the two models, due to models' different complexity.

Models Two-Layer QG Model

An idealized two-layer QG model is used to simulate geostrophic turbulence. The governing equation is the nth layer perturbation QG PV, q_n, equation: [Image Omitted. See PDF]where subscripts n = 1, 2 represent the upper and lower layers, respectively, ψ is the perturbation streamfunction, J(⋅, ⋅) is the horizontal Jacobian operator, U_n is the background zonal mean flow, ssd is the small-scale dissipation implemented by a spectral filter of enstrophy (the same as J. H. LaCasce Jr, 1996, except the exponential damping factor is 23.6), and r_ek is the bottom friction damping rate, which is only active on the lower layer (δ_mn is the Kronecker delta).

The PV of the nth layer is [Image Omitted. See PDF]where [Image Omitted. See PDF]k_d = 1/L_d, L_d is the Rossby deformation radius, and δ = H₁/H₂ is the ratio of the thicknesses of the two layers.

The background meridional PV gradient, [Image Omitted. See PDF]is due to both the planetary vorticity gradient, β, and the background vertical shear. The background flow is baroclinically unstable when the sign of β_n differs between the two layers.

We use the same parameter setting as L. Wang et al. (2016) and W. Zhang et al. (2020) for the control simulation: L_d = 15 km, H₁ = 800 m, δ = 0.25, U₁ = 0.04 m/s, U₂ = 0, r_ek = (20 days)⁻¹ and β = 1.3 × 10⁻¹¹m⁻¹s⁻¹. The nondimensional frictional rate, [Image Omitted. See PDF]and nondimensional β, [Image Omitted. See PDF]are varied by an order of magnitude to explore a wide range of flow regimes. We pick three typical simulations with r* = 0.43, 0.22, and 0.11 for displaying example results in Figures 1 and 3 below.

Figure 1. Snapshot of the upper-layer streamfunction anomaly field relative to the background zonal mean flow for a simulation with r* = 0.22 and β* = 0.073. Blue and red lines indicate the boundaries of cyclonic and anticyclonic eddies, respectively. Black dots are the centroids of the eddies. Black lines are the tracks of the eddies from the current time to when they terminate.

The model domain is doubly periodic with side length L = 1,200 km. The horizontal grid spacing is 2.3 km, which resolves the Rossby deformation radius L_d. The imposed background flow (U₁, U₂) is baroclinically unstable for all parameters used in this study. Each simulation is initialized randomly to stimulate baroclinic instability, which rapidly saturates into forced-dissipative geostrophic turbulence. The simulations are run for 30 years to forget their initial conditions and then for another 20 years with daily snapshots for analysis. Details of the model setup can be found in W. Zhang et al. (2020). The model is implemented using the Python package pyqg version 0.1.3 (Abernathey et al., 2016).

Primitive Equation Model

The PE model is an idealized configuration of the Massachusetts Institute of Technology general circulation model (MITgcm checkpoint67r; Campin et al., 2020; Marshall, Adcroft, et al., 1997; Marshall, Hill, et al., 1997), which has been used for several previous studies (Cessi & Wolfe, 2009; Cessi et al., 2010; Wolfe, 2014; Wolfe & Cessi, 2009, 2010, 2011; Wolfe et al., 2008; W. Zhang & Wolfe, 2022). The circulation simulated by this model is more complex than the QG model, as it contains multiple gyres, boundary currents, and a zonally reentrant channel flow analogous to the Antarctic Circumpolar Current.

This model is formulated in a two-hemisphere basin on an equatorial β-plane, with β = 2.3 × 10⁻¹¹ m⁻¹ s⁻¹. The model domain is approximately half the width, length, and depth of the Atlantic Ocean. The horizontal extent of the domain is 2,440 km in zonal direction and 9,760 km in meridional direction and is enclosed by vertical walls everywhere except for the southernmost eighth of the domain, where the flow is zonally reentrant. The bottom is flat with a uniform depth of 2,440 m. The model has horizontal resolution of 5.4 km and 20 vertical levels with reduced vertical grid spacing near the surface. The model is forced by zonally uniform zonal winds and a surface heat flux provided by relaxation to a zonally uniform surface temperature distribution. Both wind and temperature relaxation fields are idealizations of the forcing of the Atlantic Ocean. Dissipation is provided primarily by linear bottom drag. Since the model lacks topography, the drag coefficient must be set to a large value (damping time of 25 days) to control the transport of the circumpolar current (Gill, 1968; Tréguier & McWilliams, 1990). The details of the model setup are described by Wolfe and Cessi (2010) and W. Zhang and Wolfe (2022).

Analysis Methods Eddy Identification and Tracking

Coherent eddies are typically defined as swirling structures that can be distinguished from their surroundings (i.e., tracked) for a long time (e.g., many eddy turnover times) (Samelson, 2013). Based on this definition, we apply an algorithm to identify and track coherent eddies in both the QG and PE models. We use the eddy tracking package described by Mason et al. (2014), which provides a tracking method identical to Chelton et al. (2011) as an option. This method detects coherent mesoscale eddies as sea surface height (SSH) extrema from snapshots of SSH. The boundaries of eddies are identified as the outermost SSH contour that satisfies an area and amplitude threshold and contains no more than one local SSH extremum. Eddies are then tracked by connecting the proximal eddies within a restricted distance in successive time frames. The distance limit is determined from the local long baroclinic Rossby wave speed, and the eddies' amplitude and radius must be within a factor of 2.5 of the corresponding eddies in the last time step. Only the eddies that last longer than 30 days are used. See Mason et al. (2014) and Schlax and Chelton (2016) for details on the eddy identification and tracking algorithms.

Recent studies have used more stringent methods to identify Lagrangian coherent structures among coherent eddies (e.g., Abernathey & Haller, 2018; Beron-Vera et al., 2013; Haller et al., 2016; Y. Wang et al., 2016). Lagrangian coherent structures can trap particles over their lifetimes, while ordinary coherent eddies, such as the ones detected using the method of this section, are not necessarily Lagrangian coherent (Beron-Vera et al., 2013; Liu et al., 2019). Since we focus on the tracer mixing due to eddy stirring, we track coherent eddies as sustained Eulerian features, rather than Lagrangian coherent structures.

Coherent Eddy Diffusivity

The trajectories of coherent mesoscale eddies are used to calculate the Lagrangian diffusivity. To account for inhomogeneity and anisotropy, we use the modified version of the single-particle diffusivity tensor developed by Davis (1987, 1991), [Image Omitted. See PDF]where $${u}_{i}^{\prime }\left(t\vert \mathbf{x},{t}_{0}\right)$$ is the residual velocity of a particle at time t that was found at x at time t₀. The velocity is statistically stationary in our simulations, so $${K}_{ij}^{L}(\mathbf{x},\tau )$$ is independent of t₀. Here $${u}_{i}^{\prime }$$ is taken to be the drift velocity of the centroid of a coherent eddy to estimate a “coherent eddy diffusivity” using Equation 8. The residual velocity is calculated as $${u}_{i}^{\prime }={u}_{i}-{\left\langle {u}_{i}\right\rangle }_{L}$$, which is the deviation of the eddy velocity from the Lagrangian mean over the ensemble of eddies, which are defined differently for the QG and PE models below. The first and last 10% of the eddy trajectories are excluded for calculating the coherent eddy diffusivity. The rationale for this exclusion is given in Section 3.

Considering that the meridional drift tendency due to the β-effect is opposite between cyclones and anticyclones (McWilliams & Flierl, 1979; Nycander, 2001), the Lagrangian mean velocity is estimated separately for cyclones and anticyclones. The method, Equation 8, used in this study is different from that in Ni et al. (2020), who did not remove the Lagrangian mean velocity $${\left\langle {u}_{i}\right\rangle }_{L}$$ and only calculated the meridional diffusivity instead of the whole diffusivity tensor. We instead compute the two eigenvalues of the symmetric part of K_ij, which is necessary to reduce the bias (e.g., due to shear dispersion) of the diffusivity estimate (Griesel et al., 2014; Oh et al., 2000), and focus on analyzing the minor (i.e., second) eigenvalue of the tensor. The minor eigenvalue usually corresponds to mixing across the mean flow, is less biased by the shear dispersion, and is more relevant to eddy tracer transport since along-stream transport is typically dominated by the mean flow.

In the QG model, the Lagrangian mean velocity is the average over coherent eddies in the whole domain, since the flow is statistically homogeneous. The maximum integration time lag τ is taken as five times the eddy turnover time τ_e, [Image Omitted. See PDF]where ζ is the vorticity and $$\overline{(\cdot )}$$ is a 15-year and domain average in the upper layer. The time scale, τ, is chosen to account for the different lifetimes of coherent eddies in different simulations, as the eddy lifetime depends on the eddy turnover time, as found by W. Zhang et al. (2020). Five times the eddy turnover time is about 20–70 days in most of the QG simulations, which is comparable to the eddy lifetime and allows the Lagrangian diffusivity to asymptote to a constant value during this time period. A longer time window (e.g., seven times eddy turnover time) is also tested (not shown), which does not qualitatively change the results. The eddy tracks that are used to estimate K_ij are grouped in overlapping 100-day time windows with a 15-day interval, which is similar to the pseudo-track approach used by previous studies (e.g., Chen & Waterman, 2017; Klocker, Ferrari, LaCasce, & Merrifield, 2012; Swenson & Niiler, 1996).

In the PE model, the Lagrangian average and integration time lag τ are taken to be different from those in QG model to account for flow inhomogeneity. The domain is divided into a set of 304 × 304 km spatial bins and the mean eddy velocity is estimated as the average velocity of all coherent eddies within each bin. The Lagrangian diffusivity is then calculated using the segments of coherent mesoscale eddy trajectories that start from those bins. The size of the spatial bin (304 km) is chosen because it can contain sufficient number of coherent eddy tracks for the diffusivity calculation and also not so large that it averages away the spatial variability. The bin size evenly divides the model domain into 32 × 8 bins in the meridional and zonal directions, respectively. The integration time lag, [Image Omitted. See PDF]is chosen to ensure that most of the eddies remain in the same bin during the integration time of the diffusivity calculation. In Equation 10, Δx is half the length of the bin, |u_coh| is the speed of coherent eddies, and 30 days is the minimum lifetime of coherent eddies required by the eddy tracking method in Section 2.2.1.

Note that the Lagrangian diffusivity tensor, $${K}_{ij}^{L}$$, has four components, but only the symmetric component contributes to diffusion. The two eigenvalues of the symmetric part of K_ij give the diffusivities in the directions of the associated eigenvectors, which are often aligned parallel and perpendicular to the mean flow direction (Fox-Kemper et al., 2013; Riha & Eden, 2011), topographic gradients (Isachsen, 2011; Mechoso, 1980), or vorticity gradients (Bachman et al., 2020; J. LaCasce & Speer, 1999; K. S. Smith, 2005). The diffusivities are typically anisotropic, with the cross-stream (approximately parallel to the PV gradient) diffusivity smaller than the along-stream (approximately perpendicular to the PV gradient) diffusivity.

Eulerian Tracer Diffusivity

The Eulerian tracer diffusivity is diagnosed and compared to the coherent eddy diffusivity in the QG and PE models. In the QG model, the flow field is homogeneous and the background PV gradient is aligned with the meridional direction, so the PV diffusivity is straightforward to diagnose from the PV flux and gradient averaged over the whole domain. In the PE model, the inhomogeneity and anisotropy make it challenging to estimate the diffusivity based on a single tracer, since the tracer gradient might vanish or misalign with the tracer flux at many locations. We simulate a total of 27 different passive tracers and diagnose the tracer diffusivity using the multiple tracer inversion method of Bachman et al. (2015, 2020).

In the QG model, the PV diffusivity is diagnosed from the meridional flux-gradient relation. The coherent eddy diffusivity is compared to the upper-layer PV diffusivity since coherent eddies are detected from the upper layer streamfunction (proportional to SSH). The upper-layer PV diffusivity, κ_q, is calculated as [Image Omitted. See PDF]where $$\overline{(\cdot )}$$ indicates a 20-year and domain average, and β₁ is the upper-layer PV gradient defined in Equation 5.

In the PE model, the tracer diffusivity tensor is diagnosed by advecting multiple passive tracers, τ_α, with eddy-resolving velocity fields and inverting the course-grained flux-gradient relationship, [Image Omitted. See PDF]in a least-squares sense. In Equation 12, K is the diffusivity tensor, i and j index the horizontal spatial dimensions, α indexes the tracer number, and repeated indices imply summation. At least three different tracers with misaligned gradients are required to uniquely define K, but using more tracers (Bachman et al., 2020, suggests nine) provides a smoother estimate and reduces bias. Note that the tracer advection is done online.

The evolution equation of the α^th tracer is [Image Omitted. See PDF]where λ is the relaxation rate and $${\tau }_{\alpha }^{0}$$ is the initial condition of α^th tracer.

Following W. Zhang and Wolfe (2022), nine different initial distributions [Image Omitted. See PDF]are used, and tracers are relaxed to these nine initial distributions with relaxation time scales, λ⁻¹, of 1 year, 3 years, and 9 years, which leads to a total of 27 tracers with different distributions. Detailed reasons for choosing these initial distributions and relaxation time scales are described in W. Zhang and Wolfe (2022).

With all 27 tracers, the diffusivity tensor is solved using the Moore-Penrose pseudoinverse (Moore, 1920; Penrose, 1955)—denoted by (⋅)^†—to obtain [Image Omitted. See PDF]where $$\bar{(\cdot )}$$ indicates a 20-year and 304 km average. The 304 km coarsening scale is the same as the size of the spatial bins for the Lagrangian diffusivity estimate in Section 2.2.2, so that we can make a grid-wise comparison between the tracer and coherent eddy diffusivities. This coarsening scale is significantly larger than the mesoscale eddy scale and allows for the scale separation between the mean and eddies. Also, this coarsening scale roughly represents the largest resolved scale by models with a resolution between a half and one degree, which are typical resolutions used by many ocean climate models (Fox-Kemper et al., 2014). Note that the largest resolved scale is usually about several times the model grid spacing since hyperviscous damping is usually used to control noise at the grid scale (Danilov et al., 2019). We also tried coarsening scales of 76 and 152 km and did not find significant differences in the diagnosed tracer diffusivity (not shown).

The multiple tracer inversion method accounts for the anisotropy of eddy diffusion by diagnosing each component of a diffusivity tensor. It is shown in W. Zhang and Wolfe (2022) that the diagnosed diffusivity tensor is generic and is effective in representing the flux of an arbitrary tracer. The principal diffusivities are the eigenvalues of the symmetric part of the diffusivity tensor. The first and second eigenvalues were found to be in the zonal and meridional directions, respectively, in most regions, except where the mean flow is strong and non-zonal (e.g., the western boundary currents and southeastern part of the northern subpolar region) where the two eigenvalues align better with the along- and cross-mean flow directions (W. Zhang & Wolfe, 2022). The cross-stream diffusivity has been the focus of many oceanic studies, since the along-stream transport is mainly attributed to the mean flow (e.g., Ferrari & Nikurashin, 2010; Klocker, Ferrari, & LaCasce, 2012; Riha & Eden, 2011). Cross-stream/meridional mixing is found to have important dynamical impacts in the Southern Ocean, such as on the meridional heat transport and the overturning circulation (e.g., Bates et al., 2014; Chapman & Sallée, 2017). W. Zhang and Wolfe (2022) found that the along-stream diffusivity can be reconstructed from the cross-stream diffusivity with the suppressed mixing length formula of Ferrari and Nikurashin (2010). In this study, the cross-stream tracer diffusivity is compared with the second eigenvalue of the Lagrangian diffusivity tensor, Equation 8, estimated using coherent eddy tracks.

Coherent Eddies and Tracer Diffusivity in QG Turbulence

Coherent eddies are identified and tracked in the upper layer of the QG simulations, using the eddy tracking algorithm described in Section 2.2.1. Examples of eddy boundaries and tracks in the simulation with r* = 0.22 and β* = 0.073 are shown in Figure 1. Coherent eddies form and move randomly in QG turbulence. They also exhibit a systematic eastward drift, a net result of the intrinsic westward propagation of eddies due to beta effect and advection by the eastward mean flow.

The statistics of the temporal variation of eddy amplitude are studied following Samelson et al. (2014). Each eddy amplitude time series, A_n(t) (n = 1, 2…, N), is normalized by the time mean over the eddy lifetime, $${\widehat{A}}_{n}(t)={A}_{n}(t)/{\bar{A}}_{n}$$. Eddies with the same lifetime, T, are grouped together and their normalized amplitude is averaged, $${\widehat{A}}^{T}(t)={\sum }_{m}{\widehat{A}}_{m}^{T}(t)\ (m=1,2,\text{\ldots },M)$$, where M is the number of eddies with lifetime T.

Figure 2 shows the time series of all $${\widehat{A}}^{T}(t)$$ versus the dimensionless time, t*, where t* = t/T, with 0 ≤ t ≤ T. The normalized amplitude tends to grow during the initial stage (0 < t* < 0.1), slowly grow or decay during the mature stage (0.1 ≤ t* ≤ 0.9), and decay rapidly during the end stage (0.9 < t* < 1). This result is consistent with that found by Samelson et al. (2014) for the altimeter-tracked mesoscale eddies in the ocean. Eddies tend to experience significant merging and splitting events or interaction with submesoscale processes during their growing or decaying stages (Samelson et al., 2016; Z. Zhang & Qiu, 2018). Since the eddy behavior during the initial and final phases is atypical compared to the mature phase, we drop the first and last 10% of eddy trajectories before calculating the diffusivity below.

Figure 2. Dimensionless time series of normalized amplitude for coherent eddies. Each line is the time series of normalized amplitude of eddies with the same lifetime T, 30 days ≤ T ≤ 150 days. The dimensionless time is t* = t/T, where 0 ≤ t ≤ T, and T is the eddy lifetime. The amplitude of each eddy is normalized by the mean amplitude over its lifetime.

The Lagrangian diffusivity is calculated based on the tracks of coherent eddies over the whole domain and over 10 years using Equation 8 described in Section 2.2.2. Figures 3a–3c shows the time series of the minor (i.e., second largest magnitude) eigenvalue of the coherent eddy diffusivity tensor, $${\kappa }_{2}^{\text{coh}}$$, for three simulations with r* = 0.43, 0.22, and 0.11 (and β* = 0.073). The eigenvectors associated with $${\kappa }_{2}^{\text{coh}}$$ are nearly in the meridional direction. The coherent eddy diffusivity, $${\kappa }_{2}^{\text{coh}}$$, approaches the domain-averaged meridional PV diffusivity in the upper layer (black dashed line in Figures 3a–3c) in 20–40 days. This consistency between the diffusivity estimated from coherent eddies and the upper-layer Eulerian PV diffusivity is also found by W. Zhang et al. (2020), where coherent eddies are identified and tracked as Lagrangian coherent structures. While the PV diffusivity can become ill-defined in regions with weak PV gradients (Uchida et al., 2023), the background PV gradient is held fixed in this model and Lagrangian and Eulerian estimates of the PV diffusivity are expected to agree.

Figure 3. Top panels: upper-layer Lagrangian diffusivity calculated from coherent eddy tracks over 10 years. The blue line is the minor eigenvalue of the diffusivity tensor κ2coh ${\kappa }_{2}^{\text{coh}}$, which represents mixing in the meridional direction. Error bars are 2 times the standard error. Bottom panels: upper-layer Lagrangian diffusivity calculated from numerical particle trajectories. The orange line is the minor eigenvalue of the particle diffusivity tensor, κ2, averaged over 10 nonoverlapping time windows. The shading gives 2 times the standard error for the estimates in 10 time windows. Black dashed line indicates the Eulerian PV diffusivity, κq. The parameters for the simulations are r*= (a, d) 0.43, (b, e) 0.22, and (c, f) 0.11 and β* = 0.073.

For comparison with the coherent eddy diffusivity, the Lagrangian diffusivity Equation 8 is also calculated using the paths of a total of 1,048,576 initially uniformly spaced numerical particles advected over the whole domain. The time series of the minor eigenvalue of the particle diffusivity tensor, κ₂, are shown by the orange lines in Figures 3d–3f for the same three simulations as in Figures 3a–3c. The direction of κ₂ is also aligned with the meridional direction. The diffusivity, κ₂, estimated from particles also approaches the domain-averaged Eulerian PV diffusivity (black dashed line in Figures 3d–3f). However, it takes more than 100 days for this convergence to occur, which is more than 5 times slower than the coherent eddy diffusivity. The slow convergence of particle diffusivity is also found in some oceanic observations (Balwada et al., 2016; Klocker, Ferrari, LaCasce, & Merrifield, 2012; J. LaCasce et al., 2014; Rypina et al., 2012). The finding here implies that it is more efficient to estimate the tracer diffusivity by identifying and tracking coherent eddies than by deploying and tracking Lagrangian particles.

The slow convergence of particle diffusivity is also apparent in the two-particle relative diffusivity. In Appendix A we compute the relative diffusivity using particle and coherent eddy pairs. The relative diffusivity of particles initially increases with the particle separation, following the Richardson's 4/3 law (Richardson, 1926), until the separation becomes about 100 km. This separation is generally comparable to the energy-containing wavelength of eddies, which is 100–250 km in these simulations. Beyond the 100-km separation, the relative diffusivity fluctuates and then slowly converges to the PV diffusivity after the particle separation increases by a further 300–500 km. The relative diffusivity should behave like the single-particle diffusivity (i.e., κ₂) once the particle separation is large enough that their motions become uncorrelated (J. H. LaCasce, 2008a). The additional 300–500 km required for the relative diffusivity to converge is consistent with the average particle drift rate multiplied by the single-particle diffusivity convergence time. The relative diffusivity of coherent eddies approaches to the PV diffusivity within tens of kilometers, which is much shorter than the convergence distance of particles.

There are several possible reasons why the coherent eddy diffusivity approaches the Eulerian diffusivity faster than the particle diffusivity. First, the coherent eddy movements are due to the low-frequency component of advection by the mesoscale or large-scale flows, which dominate the mesoscale diffusivity (J. H. LaCasce, 2008a). In contrast, particles or drifters may initially spread quickly due to filaments and other small-scale processes, with their movement constrained by large-scale dynamics (e.g., the meridional PV gradient) only on longer time/length scales. In that case, the particle diffusivity increases quickly and then slowly asymptotes to the background Eulerian PV diffusivity, as is shown in Figures 3d–3f. Second, the diffusion itself has been found to be due to the chaotic movement of coherent eddies in point vortex models (Aref, 1984; Weiss et al., 1998). In this case, the diffusivity estimate based on coherent eddy tracks is directly related to the source of the diffusion. Third, particle motion in turbulent fluids have been found to transition to the diffusive regime more slowly than standard Brownian motion due to entrainment by surrounding fluid which generates long-range correlations (Chong et al., 2020; Franosch et al., 2011). If a similar effect also exists in geostrophic turbulence, it can lead to long decorrelation time scales for particle motion and a slower transition to diffusive behavior.

The correspondence between the coherent eddy and Eulerian diffusivities is examined in broad regimes of QG simulations by varying r* and β*. Since the coherent eddy diffusivity occasionally still displays some drift at the end of the integration window, a “final” value is estimated by the average over the last tenth of integration window. This final coherent eddy diffusivity (orange dashed line) and domain averaged Eulerian PV diffusivity (blue line) are compared in Figure 4. Figure 4a shows that the coherent eddy diffusivity is consistent with the Eulerian PV diffusivity over a range of frictions varying by an order of magnitude, except when the friction is very small, where the coherent diffusivity is smaller than the PV diffusivity. Figure 4b shows that the coherent eddy diffusivity matches the PV diffusivity in the regimes when β* > 0.05, while the coherent eddy diffusivity underestimates the PV diffusivity when β* < 0.05. Possible explanations are discussed in Section 4.

Figure 4. Comparison between the upper-layer PV diffusivity (blue line) and the final value of the meridional diffusivity from coherent eddy tracks (orange line) in simulations with varying (a) r* (with β* = 0.073) and (b) β* (with r* = 0.22). Error bars are two times the standard error, which is the rms of the individual standard errors over the last tenth of the integration time window.

Why Does the Coherent Eddy Diffusivity Underestimate the PV Diffusivity for Small Beta and Friction?

Coherent eddy diffusivity tends to underestimate the PV diffusivity in simulations where beta or friction are small (Figure 4). This result suggests that the mixing of tracers is not dominated by the movement of coherent eddies in these simulations. To make this observation more concrete, we examine the spatial scales responsible for the meridional PV flux, vq, by comparing the co-spectra of v and q between simulations with small and large beta (using β* = 0.01 and β* = 0.09 as an example). The co-spectrum, ⟨v, q⟩, between the meridional velocity anomaly v and PV anomaly q is [Image Omitted. See PDF]where $$\widehat{v}$$ and $$\widehat{q}$$ are the spectra of v and q, respectively. Figures 5a and 5b show ⟨v, q⟩ (blue line) in the two simulations with β* = 0.01 and β* = 0.09, respectively, in Figure 4b. The kinetic energy (KE) spectrum $$\widehat{E}=\vert \widehat{u}{\vert }^{2}+\vert \widehat{u}{\vert }^{2}$$ is also plotted (orange line) for comparison. The coherent eddy diffusivity underestimates the PV diffusivity when β* = 0.01, and matches the PV diffusivity well when β* = 0.09 (Figure 4b). The peak of ⟨v, q⟩ is at a smaller wavenumber than that of $$\widehat{E}$$ when β* = 0.01, while they overlap when β* = 0.09. This difference indicates that PV mixing is dominated by motions with scales larger than (equal to) the energy containing scale—the inverse of the peak wavenumber of $$\widehat{E}$$—when β* = 0.01 (β* = 0.09). In general, the peak of the PV flux co-spectrum is at larger scales than the peak in the energy spectrum for β* ≲ 0.07, while the two peaks overlap for larger values of β* (not shown).

Figure 5. Top panels: upper-layer spectrum of meridional PV flux (blue) and kinetic energy (orange) in the simulations with (a) β* = 0.01 and (b) β* = 0.09 in Figure 4. 2D spectra are azimuthally integrated to obtain 1D spectra. Red dashed line indicates the wavenumber below which 80% of the meridional PV flux is contained. Bottom panels: Snapshots of upper-layer streamfunction anomaly fields for the simulations shown in the upper panels. Black dashed lines are contours of spatially low-pass filtered field for the same snapshot, using a cutoff kc.

To better understand how the relative locations of the peaks of PV flux and energy spectra relate to the flow in physical (rather than spectral) space, we introduce a low-wavenumber cutoff, k_c, that satisfies [Image Omitted. See PDF]where k_max is the largest resolved wavenumber. Motions with wavenumber smaller than k_c thus account for 80% of the total PV mixing. Figures 5c and 5d show snapshots of upper-layer streamfunction anomaly for the two simulations with β* = 0.01 and β* = 0.09, respectively. Black dashed lines indicate low-pass filtered streamfunction field using the cutoff k_c. Mesoscale eddies appear as local streamfunction extrema in the shading. The eddies have similar sizes in both simulations, consistent with their similar peak wavenumbers of the energy spectrum $$\widehat{E}$$ in Figures 5a and 5b. In the β* = 0.01 simulation, the low-pass filtered contours capture large-scale structures represented by patches of positive or negative streamfunction anomalies. These patches consist of multiple same-signed eddies that are correlated through their velocity fields. These long-range correlations between multiple eddies signal that mixing is nonlocal in physical space, which is why the peak of ⟨v, q⟩ appears at smaller wavenumbers than that of $$\widehat{E}$$. In the β* = 0.09 simulation, the filtered contours predominantly capture individual coherent eddies, with fewer long-range correlations due to limited correlations between same-signed eddies. Consequently, the mixing is local in physical space and dominated by individual coherent eddies in this simulation.

Long-range correlations between eddies are also evidenced by the slope of KE spectra in Figure 6, where KE spectra of different simulations are normalized by their peak value and the corresponding horizontal wavenumbers for comparison. The slope of energy spectrum on the left of the peak (i.e., scales larger than the energy-containing scale) is shallower in the simulation with β* = 0.01 than that with β* = 0.09, indicating relatively more energy at scales larger than the energy-containing scale in the former. This result is consistent with the prevalence of large patches of streamfunction anomalies formed by long-range correlations between individual eddies in the simulation with β* = 0.01.

Figure 6. Normalized 1D upper-layer kinetic energy spectra for simulations with varying (a) r* (with β* = 0.073) and (b) β* (with r* = 0.22). The energy spectrum Eˆ=|uˆ|2‾+|vˆ|2‾ $\widehat{E}=\overline{\vert \widehat{u}{\vert }^{2}}+\overline{\vert \widehat{v}{\vert }^{2}}$, where uˆ $\widehat{u}$ and vˆ $\widehat{v}$ are the Fourier transform of the zonal and meridional velocity anomalies, and ⋅‾ $\overline{\cdot }$ indicates a 20-year average. The 1D spectrum is obtained as the azimuthal integral of the 2D spectrum. Each energy spectrum is normalized by its maximum, Eˆmax ${\widehat{E}}_{\text{max}}$. The horizontal wavenumber k is also normalized by the wavenumber corresponding to Eˆmax ${\widehat{E}}_{\text{max}}$, k0. The black dashed line indicates the slope of k3.

Overall, the energy spectrum becomes shallower than k³ at scales larger than the energy-containing scale when β* < 0.07 (Figure 6b) or when r* < 0.1 (Figure 6a). This result is consistent with the prevalence of long-range correlations between the eddies in those simulations with small beta or friction (not shown), which are similar to the results in the simulation with β* = 0.01 (Figures 5a and 5c). Geostrophic turbulence theory holds that the inverse energy cascade can be halted by friction or beta (e.g., Held & Larichev, 1996; Larichev & Held, 1995; Rhines, 1975), so smaller friction or beta leads to larger fraction of energy at large scales. According to mixing length theory, larger scale motions are more efficient in mixing due to their larger mixing lengths. In addition, a stronger inverse energy cascade tends to make eddies more barotropic (L. Wang et al., 2016). Barotropic eddies have a longer range of influence than baroclinic eddies due to their Green's functions (Bracco et al., 2004; Provenzale et al., 2008). Consequently, tracer mixing is increasingly dominated by the large eddy patches rather than individual eddies as beta or friction is made smaller, resulting in a tracer diffusivity that is larger than the coherent eddy diffusivity.

Application to Inhomogeneous 3D Ocean Circulation

The coherent eddy diffusivity and the tracer diffusivity are further compared in the PE model, described in Section 2.1.2. This model contains multiple gyres, western boundary currents, and a circumpolar current in the channel (Figure 7e) (see also W. Zhang & Wolfe, 2022). The mean and turbulent flows in this model are inhomogeneous, which is an important property of the World Ocean that idealized models typically used to study turbulence—like the QG model of the previous section—often lack. This PE model thus provides a more stringent test for the relationship between coherent eddy movement and tracer diffusivity than the QG model.

Figure 7. (a–d) Minor eigenvalue of the diffusivity tensor estimated from coherent eddy tracks in the four bins labeled (a–d) in (e) from the idealized basin circulation model. Error bars are 2 times standard error. (e) Cyclonic (blue) and anticyclonic (red) eddy trajectories in 20 years (there are sightly more cyclonic eddies than anticyclonic eddies, though many of the trajectories obscure each other here). Black lines indicate 20-year-mean SSH contours.

Coherent Eddy Diffusivity

Coherent eddies are detected and tracked using SSH snapshots at a 3-day interval in the PE model, using the method described in Section 2.2.1. The tracks of the 28 205 cyclonic and 22 678 anticyclonic eddies identified over the 20 years simulation are shown by the blue and red lines, respectively, in Figure 7e. Coherent eddy tracks cover most parts of the model domain, except for the tropics and a quiescent region in the western part of the northern subpolar gyre.

Unlike the QG model, which has a fixed Rossby deformation radius, the deformation radius in the PE model varies from about 10 km at high latitude to larger than 100 km in the tropics with a similar distribution to the World Ocean (Chelton et al., 1998; W. Zhang & Wolfe, 2022). The nondimensional planetary vorticity gradient, β*, is estimated by the Charney-Green number (Charney, 1947; Green, 1960) in the same model by W. Zhang and Wolfe (2022) (their Figure 2), which shows that β* is larger than 0.05 at most of the latitudes, except the southern part of the channel. This means that most of the regions in this model are characterized by values of β* that are in the regime where the coherent and Eulerian diffusivities agree in the QG model (Section 3). The magnitude of β* becomes larger than 0.2 (the largest value used in the QG model) in the subtropics. The regime with β* > 0.2 is never tested in the two-layer QG model because this value of β* stabilizes the two-layer system to baroclinic instability. In the PE model, we find the baroclinic instability is mainly due to that the PV gradient changes sign near the surface in the subtropics, which is a Charney-type instability (Charney, 1947; Tulloch et al., 2011) and not simulated by the two-layer QG model.

The Lagrangian diffusivity tensor is calculated using coherent eddy tracks in 304 × 304 km spatial bins, as described in Section 2.2.2. The first and last 10% of eddy tracks are excluded before calculating the diffusivity, as with the QG model. Figures 7a–7d show the minor (i.e., the second largest magnitude) eigenvalue of the Lagrangian diffusivity tensor estimated using coherent eddy tracks in the four regions shown by the black boxes in Figure 7e. The direction of the corresponding eigenvector is within 15° of the meridional direction except in the western boundary currents, where the eigenvector is about 30° from the meridional direction. The time series of coherent eddy diffusivity in this model are similar to those of the meridional coherent eddy diffusivity in the QG model. The coherent eddy diffusivity rises rapidly within the first 10–20 days and then stabilizes (with some fluctuations) within about 20–40 days in most of the regions. The “final” value of the diffusivity is taken as the average over the last three time steps (i.e., 9 days) of the integration window, described in Section 2.2.2.

Ni et al. (2020) found that the meridional diffusivity estimated from tracks of coherent eddies kept increasing with time rather than stabilizing at many locations in the ocean, which is different from what we find for the minor diffusivity in Figure 7. The reason is likely because we remove the Lagrangian mean from the eddy velocity in the diffusivity calculation in Equation 8 and estimate the minor principal component of the diffusivity tensor, both of which are necessary to reduce the bias (e.g., due to shear dispersion) of the diffusivity estimate (Griesel et al., 2014; Oh et al., 2000). We find that the coherent eddy diffusivity estimate has greater variance (i.e., larger error bars) if we do not remove the Lagrangian mean (not shown). Note that the diffusivity integration window used by Ni et al. (2020) (their Figure B1) is significantly longer than that used here, because the lifetime of eddies in our model is shorter—generally less than 100 days—than those in the ocean, where a significant fraction of eddies live longer than 16 weeks (Chelton et al., 2011). The horizontal and vertical extent of our model is about a half of the Atlantic Ocean, so the eddies in this model might be more impacted by the strong bottom friction and eddy-eddy interactions, which reduce their average lifetime. Another reason that we use relatively short integration time window is that many eddies propagate far from their initial position and experience different mixing regimes. Using too long an integration window (50 days or longer) thus leads to bias in the diffusivity estimate.

Comparison of the Coherent Eddy Diffusivity to the Tracer Diffusivity

The final coherent eddy diffusivity is compared with the local tracer diffusivity estimated using the multiple tracer inversion method, described in Section 2.2.3. The vertical profiles of the minor diffusivity in the four black boxes in Figure 7e are shown as blue lines in Figure 8. The direction of the associated eigenvector is generally in the meridional direction and is almost perpendicular to the direction of the mean flow in the upper ocean. The tracer diffusivity has a complicated vertical structure and tends to have a subsurface maximum due to the variation of eddy velocity with depth and the mean flow suppression effect, discussed in detail in W. Zhang and Wolfe (2022).

Figure 8. Vertical profiles (blue lines) of the second eigenvalue of the tracer diffusivity tensor diagnosed in the same bins as the coherent eddy diffusivities in Figure 7. Orange dashed line indicates the final second eigenvalue of the coherent eddy diffusivity tensor with shaded error of 2 times the standard error, which is the rms of individual standard errors over the last 9 days of the integration time window shown in Figure 7. Black and red dashed line indicates the depth where r = 1, and z=−hL0 $z=-{h}_{{L}_{0}}$, respectively.

The final coherent eddy diffusivity in the same region is plotted as the orange dashed line in Figure 8 with shadings of two times the averaged standard error. The tracer diffusivity has a complicated vertical structure, but the coherent eddy diffusivity provides a single estimate in each bin. This makes the comparison between the two diffusivities challenging. The reason for this difference is that the coherent eddies themselves are 3D structures and move as deep water columns. The eddy swirling velocity and the mean flow tends to decrease with depth, which causes the vertical variation of tracer diffusivity (W. Zhang & Wolfe, 2022). The coherent and tracer diffusivities are most likely to be comparable at the vertical level where the translation speed of coherent eddies is close to the eddy swirling velocity. This depth can be determined by the nonlinearity parameter, r, (Flierl, 1981) defined as [Image Omitted. See PDF]where u_rms is the rms eddy velocity estimated from a 20-year average, and c is the intrinsic speed of coherent eddies, estimated as [Image Omitted. See PDF]where u_coh is the translation velocity of coherent eddies, and $${\overline{\mathbf{U}}}^{z}$$ is the vertically averaged mean flow.

The depth where r = 1 (black dashed line in Figure 8) is the depth where the intrinsic translation speed of coherent eddies is equal to the movement speed of surrounding water parcels. At this depth, the coherent eddy diffusivity (orange dashed line in Figure 8) becomes comparable to the tracer diffusivity, though the coherent eddy diffusivity is biased low. The physical meaning of r = 1 is that it indicates the depth where eddies transition from nonlinear to linear dynamics. The depth where r = 1 is referred as the “transition depth” hereafter. If r > 1, eddies can form closed streamlines, while if r < 1, the streamlines do not close and the eddy is more wave-like (Flierl, 1981). W. Zhang and Wolfe (2022) find that the tracer diffusivity follows a mixing-length scaling (i.e., the diffusivity is equal to a depth-invariant mixing length times u_rms) where r > 1. If we express the coherent eddy diffusivity as the same mixing length times the eddy movement speed, we should expect the coherent eddy diffusivity to be comparable to the tracer diffusivity at the depth where r = 1, which is consistent with the finding here. The observation that the coherent eddy diffusivity underestimates the tracer diffusivity at this depth suggests that the mixing length for the coherent eddy diffusivity may be smaller than that for the tracer diffusivity. This discrepancy can occur when the tracer mixing is dominated by the long-range eddy correlations, as discussed in Section 4.

The coherent eddy diffusivity and the tracer diffusivity at the depth where r = 1 are compared in every 304 × 304 km bin, and their spatial distribution is shown in Figures 9a and 9b. The two diffusivities share similar spatial patterns; for example, they are strong in the western boundary current and circumpolar current regions and weak in the gyres. A quantitative comparison at all bins is given in Figure 10a. The two diffusivities are highly correlated (R² = 0.7). These results suggest that the coherent diffusivity is still meaningful in inhomogeneous 3D ocean circulations, but that it represents the diffusivity at a specific depth due to the fact that coherent mesoscale eddies generally have deep vertical extents (Frenger et al., 2015; Z. Zhang et al., 2014) and move as a whole water column. The coherent eddy diffusivity is biased low compared to the tracer diffusivity at depth r = 1 for many grid points, which might be due to the occurrence of long-range correlations between eddies, discussed in Section 4. We do observe a surface energy spectrum with a slope shallower than k³ at scales larger than the energy-containing scale in these regions (not shown), indicating the presence of long-range correlations.

Figure 9. Comparison between the distributions of the second eigenvalue of (a) coherent eddy diffusivity tensor, (b) the tracer diffusivity tensor at the level where r = 1, and (c) the tracer diffusivity tensor at the depth where z=−hL0 $z=-{h}_{{L}_{0}}$.

Figure 10. Scatter plot of the second eigenvalue of coherent eddy diffusivity tensor versus that of the tracer diffusivity tensor at the level where (a) r = 1 and (b) z=−hL0 $z=-{h}_{{L}_{0}}$ for all 304 × 304 km bins, excluding the estimates on the boundaries and negative values. The one-to-one line is dashed orange.

Estimation of the Transition Depth

The depth where r = 1 is an estimate of the vertical extent of the eddies. A direct estimate of this depth would require the vertical profile of eddy velocities, which is usually not available from ocean observations. However, the vertical structure of ocean eddies has been found to be well-described by surface quasigeostrophic (SQG) dynamics (Isern-Fontanet et al., 2008; Klein et al., 2009; Lapeyre & Klein, 2006; Qiu et al., 2016), which connects the vertical scale of eddies to their horizontal scale through the stratification. If the stratification is uniform (i.e., N = N₀) and the eddy amplitude decays to zero at infinite depth, SQG shows that the eddy streamfunction is [Image Omitted. See PDF]where $$\widehat{\psi }(k,z)$$ is the Fourier transform of the eddy streamfunction at depth z, and k is the horizontal wavenumber. The vertical decay scale of an eddy with wavenumber k is therefore |f|/(kN₀).

In varying stratification, W. Zhang et al. (2024) shows the WKB approximation to the vertical structure is [Image Omitted. See PDF]For a coherent eddy with scale L = k⁻¹, its e-folding vertical scale, h_L, can be obtained by solving [Image Omitted. See PDF]

We assume the mean horizontal scale of coherent eddies is close to the local energy containing scale, L₀, estimated following Thompson and Young (2006) and W. Zhang and Wolfe (2022) as [Image Omitted. See PDF]where $$\overline{(\cdot )}$$ denotes a 20-year mean, and η′ is the SSH anomaly relative to the mean. The vertical scale of coherent eddies is estimated as $${h}_{{L}_{0}}$$ by replacing L with L₀ in Equation 22. The vertical scale $${h}_{{L}_{0}}$$ is compared to the depth where r = 1 in Figure 11. The two depths compare well in most of the regions, except in the western boundary currents, where $${h}_{{L}_{0}}$$ is an underestimate. Outside of these regions, $${h}_{{L}_{0}}$$ is a good estimate of the transition depth and has the advantage that it can be estimated from SSH and mean stratification observations.

Figure 11. Comparison between the depth of (a) h(r = 1) and (b) hL0 ${h}_{{L}_{0}}$. (c) Zonal averages of h(r = 1) and hL0 ${h}_{{L}_{0}}$.

The tracer diffusivity at the depth $${h}_{{L}_{0}}$$ is compared with the coherent eddy diffusivity in Figures 9c and 10b. The horizontal distribution of the tracer diffusivity at depths $${h}_{{L}_{0}}$$ and r = 1 are quite similar and the coherent diffusivity is more highly correlated with tracer diffusivity at $${h}_{{L}_{0}}$$ than with the diffusivity at the depth where r = 1. This result suggests that the spatial pattern of the diffusivity at the depth $${h}_{{L}_{0}}$$ can be effectively estimated from the dispersion of coherent eddies in the ocean.

Summary and Conclusions

This study compares the coherent eddy diffusivity—the Lagrangian diffusivity estimated from the dispersion of coherent mesoscale eddies—with the tracer diffusivity in a QG model and a PE model. Coherent eddies are identified and tracked in two sets of two-layer QG simulations with varying bottom friction and β, respectively. The meridional coherent eddy diffusivity, which is the minor eigenvalue of the diffusivity tensor, is generally within 30% of the upper layer PV diffusivity except for simulations with the smallest values of friction and β. It takes 20–40 days to reach an approximately stable value, which is at least 5 times faster than the meridional Lagrangian diffusivity estimated from numerical particles inserted evenly throughout the domain. The more rapid stabilization of the coherent eddy diffusivity is likely because coherent eddies are large objects and only feel the low-frequency component of chaotic advection, which sets the diffusivity. This result suggests that coherent eddy dispersion may give an efficient estimate of the tracer diffusivity in the ocean.

The coherent eddy diffusivity underestimates the PV diffusivity for small β (β* < 0.05) and bottom friction (r* = 0.1), likely to due to prevalence of large patches of streamfunction anomaly, which play a more important role in mixing than individual coherent eddies. These large patches are formed by the coupling between multiple same-sign eddies, as a result of inverse energy cascade and barotropization in simulations with small β or bottom friction (e.g., Maltrud & Vallis, 1991; L. Wang et al., 2016; Rhines, 1975). The existence of these long-range correlations between eddies is a symptom of shallower kinetic energy spectra at scales larger than the energy-containing scale. Because these long-range correlations cover spatial scales larger than a single coherent eddy, they are more efficient in tracer mixing and play a more important role in setting the tracer diffusivity than single eddies. Consequently, the tracer diffusivity is larger than the coherent eddy diffusivity in simulations with small β and friction. We therefore expect dispersion of coherent eddies to provide a good estimate of meridional mixing where the β-effect or friction is significant, for example, in the subtropical and midlatitude oceans.

We then compare the minor eigenvalues of the diffusivity tensors estimated from coherent eddies and from tracer inversion in the PE model. The coherent eddy diffusivity is estimated in spatial bins, using the eddy tracks that pass through each bin. Estimates of the coherent eddy diffusivity stabilize in 20–40 days in most regions, which is similar to what is found in the QG model. The tracer diffusivity in the PE model has a complicated vertical structure, while the coherent eddy diffusivity provides a single estimate of the diffusivity over the whole water column impacted by the eddies. It is found that the coherent eddy diffusivity is strongly correlated with the tracer diffusivity at the “transition depth”—the depth where the translation speed of coherent eddies is equal to the rms eddy velocity. This is the depth where coherent eddies transition from nonlinear dynamics to more linear, wave-like dynamics (W. Zhang & Wolfe, 2022). The transition depth is found to be close to the e-folding vertical scale of the energy-containing eddies, which can be estimated based on SQG dynamics from SSH observations and hydrography. The tracer diffusivity at the e-folding depth of the energy-containing eddies can thus be estimated using the tracks of coherent mesoscale eddies.

The differences between the QG and PE model results are likely due to the coarse vertical resolution of the QG model (two layers). In the QG model, coherent eddies are primarily confined to the upper layer and their diffusivity is close to the upper-layer tracer diffusivity. The vertical structure of tracer diffusivity is better resolved in the PE model, and the diffusivity estimated from coherent eddies gives an estimate of the tracer diffusivity at depth—coherent eddies generally have a great vertical extent (Z. Zhang et al., 2014), which can impact the mixing over a broad range of depths. In addition, the eddy field is statistically homogeneous in the QG model, which allows us to estimate a bulk upper-layer PV diffusivity and that can be compared to the coherent eddy diffusivity estimated over the whole domain. In the PE model, lateral mixing is inhomogeneous so coherent eddy movements in each bin might be less representative of the local tracer diffusivity than those in the QG model.

A natural next step to verify the findings in this study is to compare the coherent eddy diffusivity with the tracer diffusivity in a realistic and eddy-resolving ocean model. Estimates of the full-depth tracer diffusivity in a global ocean model has been attempted by Bachman et al. (2020). It would be worth examining whether the coherent eddy diffusivity is correlated with the tracer diffusivity at the e-folding depth of the energy-containing eddies in such a more realistic model.

Mesoscale eddies have been routinely identified and tracked using satellite observations and provided as products by AVISO⁺ Altimetry (2019). These eddy tracks could be used to estimate a Lagrangian diffusivity from observations. Indeed, Ni et al. (2020) have discussed a diffusivity calculated from coherent eddy tracks, although their purpose and approach are different enough to make direct comparisons difficult. In particular, Ni et al. (2020) did not remove the systematic meridional drift of coherent eddies, which results in eddy dispersion that is often superdiffusive. In contrast, we find that the coherent eddy diffusivity in the PE model stabilizes at a constant value. Ni et al. (2020) also do not connect their eddy diffusivity to the tracer diffusivity in the ocean interior, while we find that the coherent eddy diffusivity is correlated with the tracer diffusivity at a depth that can be estimated from SSH observation and hydrography.

The coherent eddy diffusivity obtains a stable value much faster than the particle diffusivity, suggesting that using coherent eddy tracks may be more efficient than using surface drifters and subsurface floats to estimate diffusivity, though the coherent eddy diffusivity does not estimate the tracer diffusivity at the same depth as the drifters and floats do. Connecting these diffusivity estimates require knowledge of the vertical structure of eddy velocity (W. Zhang & Wolfe, 2022). Ni et al. (2023) provides an estimate of full-depth eddy kinetic energy based on composite analysis of altimeter and Argo observations. Such a data set can be used to combine the deep diffusivity estimated from coherent eddies, the surface diffusivity estimated from drifters (e.g., Zhurbas et al., 2014), and the subsurface diffusivity estimated from floats and Argo (e.g., J. LaCasce et al., 2014; Balwada et al., 2016; Roach et al., 2018) to infer the full vertical profile of the tracer diffusivity. This approach will be pursued in future work.

A AppendixRelative Diffusivity

The relative diffusivity is estimated using pairs of particles. Provided both are defined, the relative diffusivity asymptotes to the absolute (i.e., single-particle) diffusivity described by Equation 1 at sufficiently large particle separations (J. H. LaCasce, 2008a). We compute the relative meridional diffusivity using particle pairs in the QG model following LaCasce and Bower (2000) via [Image Omitted. See PDF]where y_r(t, r₀) is the meridional separation between particle pairs that have an initial separation r₀ and ⟨⋅⟩_p indicates the mean over all particle pairs. Note that Equation A1 is one-half the relative diffusivity defined by LaCasce and Bower (2000). With this factor, $${\kappa }_{p}^{r}$$ is supposed to converge to the single-particle diffusivity when particle velocities are uncorrelated (J. H. LaCasce, 2008a). Figures A1a–A1c show the upper-layer relative meridional diffusivities as a function of pair separation, $$\sqrt{{\langle {r}^{2}\rangle }_{p}}$$, for the three simulations in Figure 3. The relative diffusivity follows a 4/3 power law for separations less than about 100 km, a regime similar to the results of Richardson (1926). This separation is comparable to the energy-containing wavelength, 148 km, 175 km, and 248 km, in these three simulations with r* = 0.43, 0.22, and 0.11, respectively. Since the PV flux cospectrum generally peaks around or slightly larger than the energy-containing scale (Figure 5), the length scales associated with the particle (Lagrangian) and PV (Eulerian) diffusivities are consistent. Beyond 100 km, the relative diffusivity slowly converges to the PV diffusivity when the separation increases further by 300–500 km in these simulations. This slow convergence does not depend on the initial separation: we tried several initial separations (2.4, 20, 100, 400 km) and found that the particles always drift hundreds of kilometers apart before the relative diffusivity converges to the PV diffusivity.

Figure A1. Top panels: Upper-layer relative meridional diffusivity (blue line) calculated from particle pairs, κpr ${\kappa }_{p}^{r}$, with an initial separation of 2.4 km. The black dashed line indicates the Eulerian PV diffusivity, κq. Bottom panels: Upper-layer relative meridional diffusivity calculated from coherent eddy pairs, κcohr ${\kappa }_{\mathit{coh}}^{r}$, with an initial separation of 400 km. Blue, orange and green lines are the relative diffusivities computed using cyclonic, anticyclonic, and all eddy pairs, respectively.

The relative diffusivity is also computed using coherent eddy pairs. Coherent eddies are paired for cyclonic eddies, anticyclonic eddies, and all eddies separately with an initial separation of 400 km. This large initial separation is required to ensure a sufficiently large sample of eddy pairs, since the number of eddy pairs decreases rapidly with decreasing separation below 400 km. The relative coherent diffusivities generally approach the PV diffusivity within a 20 km increment of separation. The relative coherent diffusivity slightly underestimates the PV diffusivity in Figure A1f, which could be attributed to the smaller sample size of eddy pairs, leading to larger uncertainty of the diffusivity estimate in this simulation.

Acknowledgments

We thank Ryan Abernathey, Scott Bachman, Edmund Chang, Stephen Griffies, and Shafer Smith for helpful comments and discussions. We also thank the three anonymous reviewers for their insightful comments that improved the manuscript. This work was supported by the NSF (OCE-2048826). W.Z. was also supported by award NA18OAR4320123 from the National Oceanic and Atmospheric Administration, U.S. Department of Commerce. The statements, findings, conclusions, and recommendations are those of the author(s) and do not necessarily reflect the views of the National Oceanic and Atmospheric Administration, or the U.S. Department of Commerce. The authors would like to thank Stony Brook Research Computing and Cyberinfrastructure, and the Institute for Advanced Computational Science at Stony Brook University for access to the SeaWulf computing system, which was made possible by the NSF (OAC-1531492).

Data Availability Statement

Model configuration, analysis scripts, data files used for this study are available at https://doi.org/10.5281/zenodo.8298569 (W. Zhang, 2023).