Role of Coherent Eddies in Potential Vorticity Transport in Two-layer Quasigeostrophic Turbulence

The transport by materially coherent eddies is studied in a two-layer quasigeostrophic model of geophysical turbulence. The coherent eddies are identified by closed contours of the Lagrangian-averaged vorticity deviation obtained from Lagrangian particles advected by the flow. A series of flow regimes with different bottom friction strengths are considered---it is found that coherent eddies become more prevalent and longer-lasting as the strength of the bottom drag increases. These coherent eddies, with average core radius close to the deformation radius, propagate zonally with speeds close to the long baroclinic Rossby wave speed and meridionally with a preference for cyclones to propagate poleward and anticyclones to propagate equatorward. The meridional propagation preference of the coherent eddies gives rise to a systematic upgradient potential vorticity (PV) transport, which is in the opposite direction as the background PV transport and not captured by standard Lagrangian diffusivity estimates. The upgradient PV transport by coherent eddy cores is less than 10% of the total PV transport, but the PV transport by the periphery flow induced by the PV inside coherent eddies is significant and downgradient. This clarifies the distinct roles of the trapping and stirring effect of coherent eddies in PV transport in geophysical turbulence.


Introduction
Ocean mesoscale eddies play an important role in the transport and mixing in the ocean, and can affect wide range of oceanographic phenomena, including the stratification and overturning circulation in the Southern Ocean [2][3][4][5][6], eastern boundary currents [7][8][9][10][11][12][13][14], nutrient cycling in the upper ocean [15], and ocean uptake of carbon dioxide [16,17]. In the classical view, an "eddy" is typically defined as a deviation from a suitably defined mean state and is thought to contribute to tracer transport primarily by stirring background tracer gradients. However, with the advent of satellite observations and high-resolution ocean models, the observed eddies are increasingly considered as individual coherent structures [18][19][20]. Coherent eddies can trap fluid inside their cores and physically move water and tracers over potentially long distances. The importance of this transport mechanism to ocean circulation is currently uncertain. This study investigates the contribution of coherent eddies to tracer transport in a two-layer quasigeostrophic model of geophysical turbulence.
Coherent eddies (vortices) are a characteristic phenomenon in geophysical turbulence. Freely decaying barotropic turbulence tends to organize itself into a finite number of isolated coherent vortices starting from random initial conditions [21]. These vortices are distinct objects with closed circulations which maintain their identity for times much longer than eddy turnover time [22]. Although the formation of these coherent vortices is believed to be related to the inverse energy cascade, the specific formation mechanisms are still uncertain. Studies have shown that coherent vortices are important in determining the mixing and dispersion properties of turbulent flows [22,23]. On the one hand, the strong PV gradients around coherent vortices work as material barriers-impermeable to fluid inside vortex cores-so fluids or passive tracers are trapped inside the vortices for times comparable to the vortex lifetime. This leads to nonlocal and subdiffusive mixing [22,24]. On the other hand, in flows characterized as "vortex-dominated", Lagrangian particle dispersion [25] and Eulerian eddy fluxes [26] outside the vortex cores can also be dominated by the vortex dynamics. Despite their recognized importance, the flow regimes that favor the formation of coherent vortices and their transport characteristics remain open questions.
The majority of modern ocean models-including the oceanic components of climate models-are of insufficient resolution to resolve mesoscale eddies. The effects of eddies must therefore be parameterized in terms of large-scale quantities that are resolved by the model. The eddy flux is typically represented as a diffusive process where the flux of a tracer is related to the large-scale tracer gradient by an eddy diffusivity. However, a diffusive parameterization does not fully account for the effects of coherent eddies [22,25]. Bracco et al. [23] showed that the existence of coherent vortices leads to non-Gaussian velocity distributions, and Berloff et al. [24] revealed the nondiffusive properties of the tracer transport by mesoscale eddies in oceanic gyres, which casts doubt on using horizontally isotropic eddy diffusion for eddy parameterization. Studies by Pasquero et al. [25] and Berloff and McWilliams [27] further showed that parameterizations using stochastic models to account for effects of coherent vortices can better simulate the turbulent dispersion. It is therefore necessary to understand the dynamics and transport properties of coherent eddies in order to appropriately parameterize their effects.
Estimates of the transport by coherent eddies in ocean are quite uncertain. A factor which contribute to this uncertainty is the lack of a widely-accepted and rigorous definition of a "coherent vortex". Chelton et al. [18] assumed eddies exceeding a threshold value of a nonlinearity parameter were coherent and found that most of the mesoscale motions observable from satellite were coherent. Zhang et al. [19] defined a "coherent eddy" as a closed potential vorticity (PV) contour and estimated that the zonal volume transport due to these coherent eddies was 30-40 Sv in the subtropics. Another study by Dong et al. [28] detected coherent eddies based on the geometry of the velocity vectors and showed that the meridional heat and salt transport due to the individual coherent eddies is significant. However, the eddies in these studies were detected by Eulerian methods that can be shown to be non-objective; that is, observers in different reference frames disagree about the boundary or even the existence of the coherent eddy [29].
The particle advection studies of Beron-Vera et al. [30], Haller and Beron-Vera [31], and Wang et al. [32] have shown that eddies identified by Eulerian methods are very leaky-the defined boundary does not effectively trap the fluid inside it. Haller [29] developed an objective definition of a coherent eddy as a rotationally-coherent Lagrangian vortex (RCLV) and, in Haller et al. [1], showed that the cores of RCLVs are maxima of the Lagrangian averaged vorticity deviation (LAVD)-vorticity averaged along the trajectories of Lagrangian particles. The LAVD method reliably detects RCLVs and is significantly less computationally expensive than the previous objective coherent eddy detection methods [33]. Applying this method, Abernathey and Haller [34] showed that the RCLVs in eastern Pacific are generally smaller and shorter-lived than the eddies detected by Chelton et al. [18]. As a result, the zonal volume transport by these RCLVs was estimated to be one order smaller than Zhang et al. [19]'s estimate, and their meridional transport was negligible. Thus, the previous estimates using Eulerian methods generally overestimated the importance of coherent eddies, which causes the uncertainty about their role in tracer transport. This uncertainty extends to the fundamental question of whether coherent eddy transport is a significant player in oceanic tracer transport.
The aim of this study is to quantify the contribution of coherent eddies to tracer transport in the classical setting of two-layer quasigeostrophic (QG) turbulence. The advantage of the QG system is its simplicity, generality, and controllability which allows us to consider a range of oceanographically relevant regimes and obtain results which are reasonably independent of the vagaries of individual observational and modeling systems. Following Abernathey and Haller [34], we apply the LAVD method to detect coherent eddies. Instead of estimating the volume flux due to coherent eddies we focus on the PV flux, which is a crucial dynamical quantity in geophysical flows. PV is a materially conservative tracer, and its flux can affect the evolution of the flow. The meridional PV flux-the sum of the heat and vorticity flux-plays an important role in both driving the mean flow and maintaining the eddy enstrophy budget. Quantifying the contribution of coherent eddies to the meridional PV transport sheds light on the dynamical role of coherent eddies in the eddy-mean flow interaction. The QG model we used here mimics the dynamics in Southern Ocean, and several different values of bottom drag are used to simulate different flow regimes.
This paper is structured as follows. In section 2, we introduce the two-layer QG model, the Lagrangian diagnostics and the detection method of coherent eddies. Section 3 provides the statistics of coherent eddies in a series of different friction cases with implications for their transport properties. In section 4, we present the meridional PV transport due to coherent eddies and their contribution to the total PV transport. In section 5, we estimate the Lagrangian diffusivity due to coherent eddies, which quantifies the transport of generic passive tracers. Section 6 is the discussion and conclusions.

Model description
This study uses the Python package pyqg [35] to simulate two-layer QG turbulence. The governing equations of the model are the forced-dissipative PV evolution equations in the two layers: where the PV in layer n(n = 1, 2) is The flow is forced by a background vertical shear, ∆U = U 1 − U 2 . The Jacobian is and where k d is the inverse of the Rossby deformation radius, L d , and δ = H 1 /H 2 is the ratio of the thickness of two layers. We use the same parameters as Wang et al. [36]: L = 1200 km, L d = 15 km, H 1 = 800 m, δ = 0.25, U 1 = 0.04 m/s, U 2 = 0 and β = 1.3 × 10 −11 m −1 s −1 .
The model is run in a doubly periodic domain with a horizontal resolution of 512 × 512 grid points. A series of different values of r ek spanning an order of magnitude are used by this study. Following Wang et al. [36], r −1 ek = 10, 20, and 40 days are taken as representative values of the drag parameter and referred to as the "strong friction", "control" and "weak friction" cases, respectively.
Each case was run for 60 years with time step of 0.3125 hours, and the last 50 years were sampled daily to calculate the diagnostics. Figure 1 gives the eddy kinetic energy (EKE) time series in the upper layer over 25 years; here "eddy" refers to a deviation from the time-independent, zonal-mean flow. The standard deviation of the EKE is 4%-8% of its mean value, and all three experiments show episodes where the EKE grows or decays by 20%-30%. This implies that relatively long time averages will be necessary to obtain robust statistics.

Lagrangian particle advection and Lagrangian diffusivity
Once the flow is equilibrated, the upper layer is seeded with Lagrangian particles whose positions, x, evolve according to: where u is the velocity. The particles are initialized on a rectangular grid with a uniform spacing of 1.2 km for a total of 1,048,576 particles. They are advected for 90 (for RCLV detection) or 360 (to estimate Taylor diffusivity) days with their positions, velocities, vorticities and PV saved daily. The RCLV detection process was repeated for one hundred continuous and nonoverlapping intervals.
Particle trajectories are used to calculate the Lagrangian diffusivity. According to Taylor [37], for homogeneous and stationary turbulent flow, the Eulerian diffusivity is consistent with the single-particle Lagrangian diffusivity which can be estimated as where y i (t) is the meridional position of a particle released at y i0 at time t = 0 and N is the total number of particles.

Identification of coherent eddies
We use the Lagrangian Averaged Vorticity Deviation (LAVD) technique of Haller et al. [1] to detect coherent eddies, which are defined as a nested family of fluid tubes which rotate collectively around a common core. The cores of coherent eddies are identified as the maxima of LAVD, defined as is the vorticity of the particle located at position x 0 at time t 0 , and ζ is the domain-averaged vorticity (zero in the present case). The integrand is the vorticity deviation-the absolute difference of vorticity from its domain average-so the LAVD field is the vorticity deviation averaged along the Lagrangian trajectories over the time interval (t 0 , t 1 ). Closed contours of LAVD surrounding a maximum thus represent rings of particles rotating around the core of an eddy at the same average rate over the specified time interval. The outmost convex LAVD contour is defined as the material boundary of the eddy which is expected to trap all the fluid parcels inside during that time interval. Figure 2a and 2c show example sets of LAVD contours around a local maximum where the LAVD field is plotted at the initial positions of particles at time t 0 . The outmost LAVD contour is determined by a threshold of a Coherency Index (CI). CI measures the variation of spatial compactness of particles initially inserted within an LAVD contour and is defined as where δ 2 (t) = |X(t) − X(t) | 2 is the variance of particle positions inside an LAVD contour at time t, X(t) is the position of the particle at time t, t ∈ (t 0 , t 1 ), indicates the average over all the particles inside the contour, and || is the standard Euclidean distance. The Coherency Index compares the maximum of δ 2 (t) to R 2 /2, which is the variance of the same number of particles in a disk with radius of R. Since a disk is the most compact form of a particle cloud, any stretching or filamenting of the particle cloud will make CI negative. The definition of CI in this study is different from the one defined by Tarshish et al. [38], which measures the change of variance of particle positions at the final time relative to their initial variance. This metric tends to favor large, spatially extended initial particle clouds that subsequently collapse into compact coherent eddies. In this study, the normalization factor R 2 /2 doesn't depend on time or initial conditions and the particle spread is checked at each time step. A detailed description of the advantages of this CI will be discussed in a forthcoming paper.
Here we choose the threshold for CI as −0.75, which means we allow a maximum deviation from a disk of 75%. Figure 2b and 2d show the variation of CI on a nested set of LAVD contours for two example eddies. The value of CI typically drops rapidly below the threshold value at certain contour level, which indicates that the enclosed eddy is no longer fully coherent. Similar patterns are found in many examples, which suggests that this value is an approximate critical value below which the contour encloses incoherent particles.
In summary, the algorithm for detecting the coherent eddies within the LAVD field is the following: 1. Identify all the LAVD maxima with a minimum separation of 20 pixels (minimum separation suggested by Tarshish et al. [38]). 2. Search from each maximum in LAVD for the outermost LAVD contour satisfying CI> −0.75 using the bisection method with an initial contour level of 0.36 times the LAVD value of the maximum. 3. Remove eddies containing fewer than 200 particles (approximately equivalent to an area of 274 km 2 ). This algorithm is implemented by a Python package, floater modified from Abernathey [39].

Occurrence frequency of coherent eddies
The method outlined in the previous section was used to detect 30-day, 60-day and 90-day coherent eddies in the first 30-, 60-and 90-day periods of each of the 100 90-day intervals. The shortest interval are at least 5 times the eddy turnover time, estimated by 2π/ √ Z, where Z is the spatial mean eddy enstrophy in the upper layer. The turnover times for strong, control and weak friction cases are 5.4, 4.6 and 3.5 days, respectively. An example set of 30-day coherent eddies is shown in Figure 3, the particles inserted initially within the outer boundary of the eddies remain together in a compact shape after 30 days, with a few exceptions where eddies lose a small number of particles by the end of the 30-day period. The number of coherent eddies in the 30-, 60-and 90-day intervals are shown in Figure 4a as a function of friction parameter r ek nondimensionalized by intrinsic advection timescale L d /∆U. The coherent eddies are less numerous for smaller friction, while the coherent eddy frequency appears to plateau at large values of the friction parameter. A possible interpretation is that the flow is more energetic for weak friction, so eddy interactions are more frequent and involve stronger strain fields; both of these effects would tend to make the eddies lose coherence more quickly. Shorter-lived eddies are systematically more prevalent than longer-lived eddies. Note that each 30-day period is included in the corresponding 60-and 90-day periods, so the 60-and 90-day eddies are subsets of the 30-day eddies. In Figure 4b the eddy lifetimes are normalized by the eddy turnover time in each friction case, and the number of coherent eddies versus their normalized lifetimes is given. Both an exponential fit and a power law fit for this relationship have a p-value much smaller than 0.05, but the former one has a slightly smaller residual so is adopted by this study. The number of coherent eddies decreases exponentially with an exponent of about −0.1 as the normalized lifetime increases, implying that an eddy has a nearly 10% probability of losing coherence over each turnover time. The decay rate of coherent eddies seems to be independent on friction and may be related to the nature of the turbulent flow itself. Since the chance of an eddy loosing coherence is the same for every time interval, this implies that coherent eddies do not 'age'; that is, a recently-formed coherent eddy is as likely to loose coherence as a mature eddy.

Eddy radius distribution
Eddy radius is estimated by r = √ A/π, where A is the area enclosed by the outer boundary of the eddy (i.e., the red curves in Figure 3). The distribution of the eddy radii is given in Figure 5. The coherent eddies are generally small, with median radii slightly larger than the Rossby deformation radius L d (15 km). Note that the radii shown here is the radii of the coherent core of the eddy, which is smaller than the eddy itself. The median radius becomes slightly larger as the friction increases. An alternative eddy radial scale is the e-folding scale of its PV, which is estimated to be around 1.2 times the radius of the coherent core identified by the LAVD method.

Radial structure
The radial structure of the vorticity of isolated vortices in a two-dimensional viscous flow has been shown to have an analytical form where ω n = ω/ω 0 and r n = r/r 0 are vorticity anomaly and radial distance normalized by the amplitude ω 0 at eddy center and radius r 0 of the corresponding eddies. The form given in (10) is a similarity solution for the viscous decay of initially compact eddies with zero net vorticity; it was derived by Taylor [40], observed in laboratory experiments [41], and has been adopted to describe the radial structure of pressure anomaly-instead of vorticity-for ocean mesoscale eddies observed in satellite and Argo observations [42].
We compare the radial structure of the coherent eddies to (10), using PV instead of vorticity because the PV is the QG analog of vorticity in two-dimensional flow. The PV and radial distance r are normalized by q 0 and r 0 , respectively, where q 0 , is the PV anomaly at the center of the eddy, and the scale parameter r 0 is the e-folding scale of the PV distribution of the eddy. The radial distribution of each eddy is obtained by azimuthally averaging in contiguous radial windows of length of 0.1 r 0 . The average normalized radial PV structure for all the eddies is shown in Figure 6. The radial PV structure has a negative lobe similar to the theoretical function (10) but is steeper near the eddy center, and the minimum is shifted away from the origin. The negative lobe indicating the opposite-sign PV rings around the eddies becomes smaller as the friction increases. Weak friction average theory Figure 6. Mean radial structure of PV in coherent eddies. Panels from top to bottom are from the strong, control and weak friction cases. The orange solid curve is the PV normalized by its maximum value in the eddy center and averaged over all the coherent eddies. The blue shading spans the 20th-80th percentiles of the normalized PV at each radial distance normalized by the e-folding scale of the PV distribution. The red dashed line is the theoretical radial structure given by (10) and is normalized by its own e-folding scale.

Eddy zonal propagation velocity
The mean zonal propagation velocity of coherent eddies in three friction cases is estimated by their zonal displacements over the corresponding time interval. Studies have shown that the eddy zonal phase speeds estimated from satellite observations are predicted well by the linear baroclinic Rossby wave theory [43,44]. In the long-wave limit and taking account of the Doppler shift by the mean flow [36,44], the baroclinic Rossby wave speed is estimated as where U b = (δU 1 + U 2 )/(1 + δ) is the depth-averaged mean flow. Figure 7 shows the zonal propagation velocities of coherent eddies are within one standard deviation of c BC , although the velocities are systematically larger than c BC for all but the two smallest values of the friction parameter. This may be due to the tendency of the eddies to concentrate in the faster upper layer when friction is strong [36].

Eddy meridional displacement
The meridional displacement (Y) of the 30-and 60-day coherent eddies versus their temporal and spatial mean PV anomalies are shown in Figure 8. The left and right cluster of high-density bins represent the anticyclones and cyclones. The meridional displacement is correlated with the sign of the PV anomaly, with cyclones (anticyclones) tending to drift northward (southward). A two-sample t test shows that the means of meridional displacement of cyclones and anticyclones are significantly different (p 0.05). In all three cases approximately 60% of 30-day cyclones propagate poleward and 60% of 30-day anticyclones propagate equatorward. The identical but reversed percentages in all the three cases illustrates the symmetry between cyclones and anticyclones in QG system. The magnitude of meridional displacement of coherent eddies has a negative linear relationship with their magnitude of PV, implying that stronger The preference for cyclones and anticyclones to drift in opposite directions has been observed in a global investigation of ocean mesoscale eddies [18], which indicates that this is a universal phenomenon for ocean eddies. Indeed, McWilliams and Flierl [43] demonstrated that the secondary circulation induced by the leading and trailing streamfunction anomalies of an anticyclone acts to advect it southward. They showed that nonlinearity is essential for this self-advection process, and beta effect acts to induce the initial streamfunction dispersion. Similar mechanism is also discussed by Cushman-Roisin and Beckers [46] with clear schematics. A more fundamental interpretation is that a vortex has tendency to return to a rest latitude where the ambient PV matches that of the eddy [43,47]. Since this meridional drift reduces the eddy's PV anomaly, the drift leads to decay of the eddy and may lead to loss of coherence. Whether this effect is a significant decay mechanism still requires quantitative assessment.

Transport by trapping
The Lagrangian trajectories of the coherent eddies can be used to estimate their PV transport. The total Lagrangian advective PV flux in the meridional direction is where v and q are the velocity and PV anomalies for each Lagrangian particle and the overbar is an ensemble average over all the Lagrangian particles during each 30-, 60-or 90-day period. The PV flux due to the coherent eddies is where A is a masking function that is 1 for particles inside coherent eddies and 0 outside. A is the fractional area occupied by coherent eddies, so Q c is the mean flux per coherent particle due to coherent eddies-the total PV flux due to coherent eddies is obtained by multiplying (13) byĀ. The PV flux per particle due to incoherent motions is where 1 − A is fraction of the domain outside of coherent eddies. Figure 9 shows the statistics of coherent PV flux, Q c , calculated in each of the 30-, 60-and 90-day periods and the corresponding time averaged incoherent PV flux, Q inc , which is close to the total PV flux since the area occupied by coherent eddies is small. The coherent PV flux Q c is systematically positive (upgradient), while the mean incoherent PV flux Q inc is always negative (downgradient). This means that the PV transport inside coherent eddies is systematically opposite to the transport in the background flow, implicating that coherent eddies play a special dynamical role. The total meridional PV flux by coherent eddies is AQ c , so the fraction of the total flux due to coherent eddies is AQ c /Q t . Figure 10a shows the median of the ratio AQ c /Q t is systematically less than 10% (upper limits are less than 20 %) even in very high friction cases. Thus, the PV transport by coherent eddies is small compared to the total PV transport, but is systematically in the opposite direction.
The upgradient coherent PV flux is due to the systematically poleward translation of cyclones and equatorward translation of anticyclones (cf. section 3). Since the total PV transport must be downgradient in statistical equilibrium, the downgradient incoherent PV transport over the whole domain must be much larger to compensate for the upgradient transport due to coherent eddies. This sets a dynamical constraint on the contribution of PV transport due to coherent eddies. However, were flows less zonally uniform, advection of enstrophy by mean flow might become important [48]. Holland and Rhines [49] showed that strong advection of enstrophy by gyres leads to both up-and down-gradient PV fluxes. Thus, the upgradient PV transport by coherent eddies might be more significant in regions with strong nonuniform mean flows. 30-day Q inc 60-day Q inc 90-day Q inc 30-day Q t 60-day Q t 90-day Q t Figure 9. Net meridional PV flux (per coherent particle) due to coherent eddies as a function of the nondimensional friction parameter r ek L d /∆U. Blue, orange and green dashed curves correspond to the median of PV flux for 30-, 60-and 90-day coherent eddies, respectively. Error bars span the 20th-80th percentiles over 100 time intervals. The dash-dot red, purple and grey lines indicate the averaged incoherent meridional PV flux over the 100 30-, 60-and 90-day intervals, respectively.

Transport by stirring
The small contribution of PV transport inside coherent eddies does not necessarily mean that coherent eddies are not important in meridional PV transport. In addition to trapping fluid inside their cores, coherent eddies also induce the flow in the far field. Transport by these flows can also be thought of as "due to" the coherent eddies. Studies of two-dimensional turbulence have proposed a two-component fluid view [23,50] for modeling flow with coherent vortices immersed in background turbulence. The main idea for this view is to decompose the Eulerian flow into two dynamical components: one is generated by coherent vortices and the other is generated by the background vorticity field [25]. This allows us to estimate the transport by the flow induced by coherent eddies including the far field.
The flow induced by the coherent eddies is calculated using piecewise PV inversion (PPVI) [51][52][53]. PPVI has been widely used in diagnosing the effect of portions of the PV in geophysical systems, such as effect of a PV anomaly at the tropopause on the cyclogenesis [54]. The PV inversion process is linear in QG, so the total flow is the sum of flows induced by all the portions of the PV. Specifically, the streamfunction for the total flow is decomposed as where ψ c and ψ inc are the flow generated by the coherent eddies and the background PV field, respectively. The former ψ c is calculated by the inversion of the PV field q c which is nonzero only in the upper layer inside coherent eddies. The incoherent flow, ψ inc , is derived from the remaining PV field, including the lower layer. An example of the flow field ψ c is shown in Figure 11-the PV inside coherent eddies induce the flow throughout the whole domain.
The meridional PV flux due to ψ c is v c q, where q is the total (coherent plus incoherent) PV field. The PV transport by the flow induced by coherent eddies is 10-60% of the total PV transport, with stronger friction associated with higher ratios (Figure 10b). In addition, this transport is systematically downgradient so the net PV transport due to both the stirring and trapping effect of coherent eddies is downgradient. The meridional PV transport due to coherent eddies is therefore primarily due to stirring rather than the drift of their cores, but these two components are actually correlated with each other. More than 80% of the induced PV transport is localized within 3 radii of the coherent cores; in this region, more than 50% of the EKE is induced by the PV inside the eddy cores.

Tracer transport by coherent eddies
The PV distribution is highly correlated with the eddies-indeed, the eddies are PV extrema almost by definition. The role of coherent eddies in the transport of a "generic" passive tracer (i.e., a tracer whose initial distribution is independent of the distribution of the eddies) is also of interest.
The single-particle Lagrangian diffusivity is calculated by (7) in a 360-day period. The time series of the Lagrangian diffusivity in the three friction cases are shown in Figure 12 and compared to the Eulerian meridional PV diffusivity where an overbar indicates an area and time average. The Eulerian diffusivity becomes larger as the friction weakens. The Lagrangian diffusivity of all the three simulations increases to a maximum within the first 10 days, and then decays and asymptotes to a constant value which is close to the domain-mean Eulerian diffusivity. We can use a small modification of the Taylor diffusivity (7) to calculate the coherent diffusivity; that is, the diffusivity due to coherent eddies. The modified diffusivity is where A i is a masking function which is unity for particles inside coherent eddies and zero otherwise.  The coherent diffusivity increases to a maximum within the first 1-2 days, then drops to a local minimum near the 5th day, and then increases again ( Figure 13). The period of this diffusivity oscillation is on the same order of the eddy turnover time, so may be due to the swirling motion of the eddies. After this oscillation, the coherent diffusivity converges to the Eulerian diffusivity much faster than the Lagrangian diffusivity averaged over all the particles. However, this does not mean the coherent diffusivity represents the total Lagrangian diffusivity. The contribution of coherent eddies to the total diffusivity is AK C /K, which is the same as the 'coherent diffusivity' defined in Abernathey and Haller [34]. Indeed, the fraction of Lagrangian diffusivity due to coherent eddies is less than 2% for the three friction cases shown in Figure 13. The estimate for 30-day eddies is slightly larger than the estimate by Abernathey and Haller [34] but is still not a significant contribution to the total diffusivity. The smallness of the coherent diffusivity is due to the small fractional area, A, occupied by coherent eddies.  The coherent diffusivity is mostly positive, which is inconsistent with the upgradient coherent PV flux discussed in section 4. This implies that the PV transport by the coherent cores is not a diffusion process and isn't captured merely by the dispersion of Lagrangian particles in the cores. Taylor's [37] derivation assumes that the initial tracer gradient is uniform, and that the particle motion is independent of the tracer distribution. This is clearly not the case for PV transport by coherent eddies since the eddies themselves are the PV extrema and their motion depends on the sign of their PV. Taylor's theory therefore is not applicable to PV transport by coherent eddies since its assumptions are violated.

Discussion and conclusions
This study uses the objective Lagrangian coherent eddy detection method of Haller et al. [1] to identify the coherent eddies in two-layer quasigeostrophic turbulence. Materially coherent eddies are detected for a wide range of parameters, with stronger friction associated with a greater prevalence of coherent eddies. The coherent eddies take the form of a small (near the deformation radius) materially coherent core imbedded in a larger area of PV anomaly characterized by incoherent particle motion.
The cyclones and anticyclones have distinct preferences for meridional propagation. Cyclones systematically propagate poleward and anticyclones propagate equatorward, which gives rise to upgradient PV transport. This upgradient PV transport is the opposite of the general PV transport by the background flow and is not captured by the Lagrangian diffusivity of Taylor [37]. This implicates a special dynamical role played by the coherent eddies, which is distinct from the diffusion processes. The PV transport due to the drift of coherent eddy cores is small compared to the total PV transport, but these coherent cores can induce the flow in the far field which gives rise to significant downgradient PV transport in the periphery of the eddy cores. The former component of the transport is due to trapping, while the latter is due to stirring. This study shows that the stirring transport by coherent eddies is more likely to be parameterized as a diffusion process, while the trapping transport-although itself small-might determine the evolution of eddies themselves. How the trapping and stirring transport correlate with each other and how they are correlated with the formation and decay of coherent eddies needs further study.
The meridional propagation tendency of coherent eddies relies on the β-effect, and PV is itself an active tracer, which makes the PV transport special. It would also be interesting to investigate the role of coherent eddies in the transport of tracers like heat and biochemical constituents which have been shown to correlate with the signs of PV of mesoscale eddies in the ocean. Thompson and Young [26] investigated the eddy heat flux in a two-mode, f -plane QG model and proposed that the downgradient heat flux is due to the systematic propagation of hot anticyclones poleward and cold cyclones equatorward, which is opposite to our cases. However, Dong et al. [28] showed that the heat transport by mesoscale eddies in the ocean is significantly upgradient near the tropics. It would be worthwhile to investigate this process in the β-plane QG model.
The flows considered in this study are statistically homogeneous, while the transport and mixing in ocean gyres are anisotropic and inhomogeneous [24]. Previous studies on oceanic gyres have shown that tracer transport by eddies can differ from homogeneous case in many ways, with subdiffusive single-particle dispersion due to coherent structures [24], mixing nonlocality within jets [55] and the transport barriers in eastward jets. Whether the trapping or stirring transport of coherent eddies is significant in such systems and how coherent eddies interact with the zonal jets or meridional western boundary currents would be interesting to explore. Better understanding of the universal transport properties of coherent eddies will improve our ability to diagnose the turbulent flow structures and parameterize the eddy fluxes in coarse-resolution climate models.