Meridional diffusivity is assessed for a baroclinically unstable jet in a high-latitude idealized circumpolar current (ICC) using the Model for Prediction across Scales Ocean (MPAS-O) and the online Lagrangian in Situ Global High-Performance Particle Tracking (LIGHT) diagnostic via space–time dispersion of particle clusters over 120 monthly realizations of O(106) particles on 11 potential density surfaces. Diffusivity in the jet reaches values of O(6000) m2 s−1 and is largest near the critical layer supporting mixing suppression and critical layer theory. Values in the vicinity of the shelf break are suppressed to O(100) m2 s−1 because of the presence of westward slope front currents. Diffusivity attenuates less rapidly with depth in the jet than both eddy velocity and kinetic energy scalings would suggest. Removal of the mean flow via high-pass filtering shifts the nonlinear parameter (ratio of the eddy velocity to eddy phase speed) into the linear wave regime by increasing the eddy phase speed via the depth-mean flow. Low-pass filtering, in contrast, quantifies the effect of mean shear. Diffusivity is decomposed into mean flow shear, linear waves, and the residual nonhomogeneous turbulence components, where turbulence dominates and eddy-produced filamentation strained by background mean shear enhances mixing, accounting for ≥80% of the total diffusivity relative to mean shear [O(100) m2 s−1], linear waves [O(1000) m2 s−1], and undecomposed full diffusivity [O(6000) m2 s−1]. Diffusivity parameterizations accounting for both the nonhomogeneous turbulence residual and depth variability are needed.
Baroclinic instability produces mesoscale eddies (Smith and Marshall 2009, and references therein) that are ubiquitous in the global ocean (Chelton et al. 2011; Petersen et al. 2013). The eddies occur at the scale of the baroclinic Rossby radius of deformation (RRD), which is typically a few tens of kilometers at high latitudes (Chelton et al. 1998). These eddies have associated time scales ranging from days, corresponding to Lagrangian time scales (Balwada et al. 2016), to a few weeks, corresponding to eddy velocity decorrelation time scales (Klocker and Marshall 2014), to potentially hundreds of days for the eddy life cycle (Petersen et al. 2013). Both the space and time scales of mesoscale eddies are very short relative to ocean basin and ocean climate scales, suggesting that a separation of scales might exist. Temporal-scale separation between resolved and unresolved scales, for example, the mean and eddy flow, would imply that Reynolds averaging is reasonable (Bachman et al. 2015) and validate assumptions used to develop mixing suppression via weakly nonlinear theory (Ferrari and Nikurashin 2010; Klocker and Abernathey 2014). A clean scale separation, if it exists, would inform climate modeling efforts because mesoscale eddies strongly contribute to transport processes within the global climate system by enhancing mixing via stirring along isopycnals. For example, mixing by mesoscale eddies helps drive carbon sequestration into the deep ocean (Gnanadesikan et al. 2015), overturning circulations in the Southern Ocean (Naveira Garabato et al. 2011), and is likely a leading-order process for basin-scale ocean circulation and transport (Marshall and Speer 2012). Overall, an accurate parameterization of lateral mixing is necessary to understand the general circulation of the ocean (Fox-Kemper et al. 2013, and references therein).
The role of mesoscale eddies in the climate is complex and exhibits both advective- or diffusive-like transport characteristics (Garrett 2006). For example, mixing averaged to long time scales may result in advective-like behavior due to the along-tracer contour component of diffusivity, that is, an eddy-induced velocity (Gent and McWilliams 1990). Similarly, gradients in diffusivity may also produce advective-like fluxes (Garrett 2006). In contrast, at sufficiently large spatial and temporal scales, mixing may be parameterized by a downgradient Redi (1982) diffusivity that approximately describes the irreversible mixing corresponding to the covariance of eddy and scalar fluctuations. The time-averaged regime is typically quantified in mixing studies and presumes notions of scale separation and decorrelation (Corrsin 1975; Papanicolaou and Pironeau 1981; LaCasce 2008; Klocker et al. 2012b; Chen et al. 2015).
A decomposition of the Redi (1982) diffusivity term into slow and fast processes, to the authors’ knowledge, does not exist. Part of the challenge is that mixing is a time integral of cascading processes whereby stirring stretches material lines that, in turn, produce strong tracer gradients that are ultimately diffused across tracer contours. Mixing arising from the residual cross terms of a mean and eddy diffusivity decomposition (e.g., geostrophic turbulence by nonlinear eddies affected by the mean flow), may not be well understood because of the challenge of untangling complexity originating from the temporal integration, for example, because of chaos (Aref 1984; Pierrehumbert 1991; del Castillo-Negrete and Morrison 1993; Ngan and Shepherd 1997; Prants 2014). Strong progress in understanding interactions between eddies and the mean flow, especially with respect to mixing suppression, has been made (Bates et al. 2014, and references therein). However, the processes producing mixing suppression via eddy and mean flow interactions are still not fully understood, although work on contributing processes is occurring (e.g., understanding shear contributions to suppression) (Srinivasan and Young 2014). A deeper understanding of the effect of the role of residual, eddy, and mean flow interactions on diffusivity is needed to fundamentally understand mixing due to mesoscale eddies. Physical contributions to diffusivity via a diffusivity decomposition could quantify the relative mechanisms and roles of mixing in mesoscale eddying regimes. Ultimately, this knowledge will be needed to conduct accurate climate simulations at global and decadal scales that cannot directly resolve the mesoscale and instead rely upon parameterization.
In this paper, we hypothesize that mixing for an idealized baroclinic jet is due to the combination of interactions between mesoscale eddies and the background mean flow. Hence, we quantify interactions between mesoscale eddies and a sheared mean flow for an idealized baroclinic jet, producing a decomposition of mixing into mean shear, wave, and turbulent components via assessment of diffusivity for temporally decomposed Eulerian velocities. Short temporal scales produce initial filamentation via eddy motions as quantified using a high-pass filter. Long temporal scales are quantified using a low-pass filter, for example, and are associated with large-scale straining of filaments by time-mean shear. The resultant residual diffusivity quantifies turbulent contributions to mixing for interactions between the slow and fast Eulerian velocity time scales. The goal of the paper is to better understand mixing by quantifying the importance of the residual diffusivity following the removal of mean and eddy mixing from the full diffusivity, with evaluation within the context of mixing suppression and critical layer (MSCL) theory (e.g., Ferrari and Nikurashin 2010).
To this end we begin by reviewing the literature related to mean flow effects on diffusivity, that is, diffusivity suppression, in section 2. We design an experiment to study mixing in an idealized circumpolar current (ICC) in section 3a via application of Lagrangian in Situ High-Performance Global Particle Tracking (LIGHT), as detailed in section 3b. Lagrangian trajectories derived from integration of the full and decomposed eddy and mean Eulerian velocity fields are used to assess diffusivities and compute an associated residual diffusivity, as outlined in section 4. Results are presented for the ICC flow in section 5, highlighting the hydrodynamics and isopycnal mixing arising from the background flow, the eddies, and the residual interaction of the two. These results are used to analyze MSCL theory using the residual diffusivity decomposition in section 6, followed by concluding remarks in section 7.
2. Diffusivity suppression: MSCL theory
Diffusivity suppression (i.e., the resultant mixing arising from a velocity field containing turbulence and a background mean), can be estimated with MSCL theory (Bretherton 1966; Green 1970; Haynes 1985; Lozier and Bercovici 1992; Marshall et al. 2006; Shuckburgh et al. 2009a,b; Abernathey et al. 2010; Ferrari and Nikurashin 2010; Sallee et al. 2011; Naveira Garabato et al. 2011; Klocker et al. 2012a; Abernathey and Marshall 2013; Klocker and Abernathey 2014; Tulloch et al. 2014; Bates et al. 2014; Chen et al. 2014; Pennel and Kamenkovich 2014), which is derived from weakly nonlinear critical layer theory as opposed to a fully nonlinear critical layer theory (e.g., Benney and Bergeron 1969; Davis 1969; Warn and Warn 1978; Smith and Bodonyi 1982; Feldstein and Held 1989; Churilov and Shukhman 1996; Stastna and Lamb 2002). MSCL analytical derivations are provided by Ferrari and Nikurashin (2010) and Klocker et al. (2012a). The theory is empirically applied with depth-variable diffusivity profiles in the global ocean (Bates et al. 2014) and is used for observational estimates of global (Klocker and Abernathey 2014) and Southern Ocean isopycnal diffusivity (Tulloch et al. 2014), including via Lagrangian analysis (Griesel et al. 2010; Klocker et al. 2012a; LaCasce et al. 2014; Griesel et al. 2014, 2015; Balwada et al. 2016). MSCL theory postulates that mixing is suppressed by the reduction of mixing length via a mean flow, and diffusivity is largest when eddies propagate at the same speed as the mean flow, for example, following the discussion and derivation in Bates et al. (2014), a suppressed cross-stream diffusivity K is
where Kmax is the maximum unsuppressed diffusivity arising from the turbulent eddy speed urms, Doppler-shifted eddy phase speed cw, mean zonal velocity , mixing length Lmix (Klocker and Abernathey 2014), and Γ is an O(1) mixing efficiency parameter based on theory or fit to data, for example, Γ = 0.5 (Taylor 1915), Γ = 0.35 (Klocker and Abernathey 2014), etc. Note that suppression modifies the mixing length and, by association, the diffusivity. Consequently, rapid reduction of urms can inhibit a subsurface diffusivity maxima (Griesel et al. 2015). An alternative expression for Kmax in terms of eddy kinetic energy is
where TL is a Lagrangian time scale, and the form of eddy kinetic energy follows Ferrari and Nikurashin (2010). Note that TL is effectively a modeling parameter. Within the framework of Bates et al. (2014), it is presumed a constant corresponding to a decorrelation time scale; the vertically varying integral time scale has also been used (Wolfram et al. 2015). The two diffusivity models in (2) and (3) suggest different vertical structure if Lmix and TL are assumed to be uniform in the vertical, although for generality TL must be assumed to vary in the vertical and may be estimated via integration of the autocorrelation (Wolfram et al. 2015) or using assumptions of proportionality to urms (Bates et al. 2014), with similarity to (2). Thus, parameter fits within the MSCL framework can be performed in terms of two unknown parameters: b and Lmix or b and TL. Note that Γ can be absorbed into Lmix or TL during curve fitting.
The mixing suppression prescribed in (1) is due to the absolute phase speed for nonlinear eddies cw, the mean flow , and a constant b modulated by . Mixing is predicted to be at a maximum when , a location often referred to as the steering level, critical depth, or critical level, although this maximum does not necessarily have to occur depending upon attenuation of urms with depth (Chen et al. 2015; Griesel et al. 2015; Balwada et al. 2016). In an idealized Antarctic Circumpolar Current, which is a flow with a strong mean, mixing is suppressed at the surface and increases to a critical depth (Abernathey et al. 2013). Building on this foundation we therefore describe an idealized experiment designed to study mixing suppression and further extend understanding of mean flow effects on diffusivity in unstable, zonal baroclinic jets by evaluating the diffusivity for temporally decomposed Eulerian velocity fields to produce a diffusivity decomposition.
3. The idealized circumpolar current experimental design
a. Model configuration
Previous idealized Southern Ocean flows (Abernathey et al. 2011; Stewart and Thompson 2013; Abernathey et al. 2013; Saenz et al. 2015) motivate the ICC configuration. To quantify the role of mean flow effects on meridional isopycnal diffusivity, the ICC is configured such that baroclinic eddies evolve under the influence of a strong zonal-mean flow (Fig. 1). Additionally, the ICC is also used to understand the role of eddy–mean flow interactions on the momentum budget (Ringler et al. 2017).
The ICC is simulated using the unstructured grid Model for Prediction across Scales Ocean (MPAS-O; Ringler et al. 2013). MPAS-O is the ocean model component of the U.S. Department of Energy’s Accelerated Climate Model for Energy (ACME) and uses variable-resolution spherical Voronoi tessellations (Ju et al. 2011) that are capable of simulating ocean motions across a range of spatial and temporal scales.
The ICC simulation is conducted on a planar zonally periodic domain that is 1000 km wide with a meridional extent of 2000 km. The shelf break is 500 km from the southern boundary with a tanh half-width of 100 km. The on-shelf water depth is 500 m, and the deep water depth is 2500 m. Regular hexagons with cell centers separated by 5 km are used to tessellate the plane in a 200 zonal by 460 meridional cell pattern. The vertical resolution varies from approximately 0.6 m at the surface to 90 m at depth using 100 vertical layers; approximately half the vertical layers are in the upper 250 m with a z-star style arbitrary Lagrangian–Eulerian coordinate (Petersen et al. 2015). The baroclinic time step is 3 min with 28 split-explicit, barotropic subcycles per baroclinic time step.
A westerly wind stress of 0.2 N m−2 and easterly wind stress of −0.05 N m−2 are applied using sin2 profiles with 1600-km wavelengths centered at y = 1400 and 400 km, respectively, following the zonal Southern Ocean winds approximation of Large and Yeager (2009) as used by Stewart and Thompson (2013). A linear equation of state with linear thermal expansion coefficient of α = 0.255 kg m−3 °C−1 is employed. Interior temperature restoring within sponge layers at the north and south boundaries are used to drive an overturning circulation. The northern and southern wall interior restoring is similar to Saenz et al. (2015), but with ze = 1 km and an exponential decay e-folding length scale of Le = 80 km. The surface temperature is restored to the tanh profile of Saenz et al. (2015), but with Tb = 1.0 °C and a piston velocity of 1.93 × 10−5 m s−1. A β-plane approximation is used with f0 = −1 × 10−4 s−1 and β = 1 × 10−11 m−1 s−1 (Abernathey et al. 2011; Stewart and Thompson 2013). Additional configuration details are summarized in Ringler et al. (2017).
General MPAS-O details are given in Ringler et al. (2013) and Petersen et al. (2015). Grid-scale enstrophy is removed from the simulation by rotational and divergent hyperviscosities of 3.9 × 108 and 3.9 × 109 m4 s−1, respectively (Ringler et al. 2017). Bottom drag is parameterized by quadratic bottom drag with cd = 3 × 10−3. The surface mixed layer is parameterized via the K-profile parameterization (Large et al. 1994) via CVMix (Griffies et al. 2015) with a background viscosity of 10−4 m2 s−1 and diffusion of 5 × 10−6 m2 s−1.
The model has a constant salinity of 34 psu and uses the initial temperature profile of Saenz et al. (2015). The simulation is started from rest and is integrated under its steady forcing for 100 years at 20-km resolution, interpolated to 5-km resolution, and then integrated for another 35 years. Across the last 10 years of simulation, particles are released from their initialized state every month to develop 120 ensemble members from 1.2 × 108 particle trajectories with daily output.
b. Lagrangian in Situ Global High-Performance Particle Tracking
Lagrangian particle tracks are simulated using MPAS-O’s online LIGHT analysis member (Wolfram et al. 2015). This approach computes particle trajectories during simulation run time, ensuring that Lagrangian particles are advected with the same spatial and temporal fidelity as the Eulerian model. Particle trajectories are purely diagnostic and are integrated directly from a particular Eulerian velocity field without parameterization of diffusion. This online particle tracking approach is preferable to offline post hoc analysis because Lagrangian trajectories obtained from undersampled velocity output are less accurate and accuracy is not recoverable, even by addition of stochastic particle diffusion (Qin et al. 2014). The incremental cost for the simulation is an increase in runtime of approximately 15% for a particle time step of 3 min with daily output. The relatively small incremental cost for use of LIGHT is possible because its input, output, and computation are all fully parallel via the message passing interface (MPI), similar to the host MPAS-O dynamical core.
Particles are advected along isopycnal surfaces shown via the time-averaged contours in Fig. 2. At the beginning of each simulated month, 1 × 106 particles are seeded in a 200 × 460 × 11 pattern by placing one particle at each cell center across 11 potential density surfaces. The potential density surfaces (kg m−3) are shown in Fig. 2 with values ranging from 1028.5 to 1030.0 in increments of 0.15. Isopycnally constrained advection ensures, by construction, that diagnosed mixing is representative of quasi-adiabatic processes along interior surfaces of potential density. Advection on outcropping layers is performed using the velocity on the nearest potential density surface resulting in some adiabatic mixing with the surface layer. Horizontally, particles remain in the domain for Courant–Friedrichs–Levey numbers less than one by construction because Wachspress interpolation from cell vertices is employed (Gillette et al. 2012).
4. Analysis methods
a. Particle diagnostics
Particle diagnostics such as the mean velocity, eddy velocity, decorrelation time, and integral time scale are computed as in Wolfram et al. (2015) based on the statistics of ensembles of particle clusters. Quantities are plotted in buoyancy space because a potential density coordinate is a natural frame of reference to analyze isopycnal mixing. Furthermore, this choice is consistent with the diagnostic frame of reference because isopycnally constrained trajectories from LIGHT are used to compute Lagrangian diagnostics.
Cluster diffusivity is computed using the methods in Wolfram et al. (2015) via the time rate of change of absolute dispersion relative to the particle cluster center of mass, which is equivalent to double-particle statistics (LaCasce 2008) in the meridional j direction, that is,
where particles are grouped into clusters c and indexed with a cluster particle ID i at the start of each realization s at a time t = 0 where t corresponds to the time elapsed since the start of the realization. Note that for simplicity spatial coordinates x are in buoyancy space such that x = (x, y, ρ), where x is the zonal position, y is the meridional position, and ρ is a potential density surface. Diffusivity is plotted for each c at the cluster’s initial position xc. There are Np(c) particles for each cluster centered on xc that sample a horizontal area within a given radius of xc, and there are Ns realizations obtained via temporal sampling of the quasi-steady flow. The overbar operator
corresponds to averaging over each individual particle i in the cluster containing Np(c) particles for a particular realization s. Each cluster may have a different number of particles, that is, Np(c), depending upon the cluster initial location xc and the initial locations of the particles x(t = 0), which are the same, independent of realizations (i.e., particles are seeded at the same locations at the start of each realization). The brackets
correspond to averaging over Ns realizations that sample quasi-steady time for the simulation. Thus, κyy corresponds to diffusivity computed for positions in space using particles averaged over the local spatial vicinity and time.
The diffusivity is evaluated at an approximate decorrelation time t = Td, where the autocorrelation settles to zero, which yields an assessed diffusivity κyy(xc). The autocorrelation is given by
where the Lagrangian cluster-based eddy velocity urms is
The size of the cluster should correspond to the largest scales of mixing, for example, several multiples of the Rossby radius of deformation, to ensure statistical accuracy by including as many particles as possible, provide a reasonable time to decorrelation, and not impair the spatial fidelity of the computed spatial diffusivity field (Wolfram et al. 2015). A cluster radius of 100 km is sufficient to this end. Particles on outcropped potential density surfaces at the start of each realization are excluded from the calculation of cluster diffusivities to minimize the influence of the surface layer and, instead, emphasize the eddy-driven mixing along isopycnals.
The meridional diffusivity, κyy, is used to quantify mixing in the ICC and, thereby, measure the importance of eddies in setting the vertical and meridional structure of the ocean climate (Stewart and Thompson 2013; Abernathey et al. 2013; LaCasce et al. 2014; Tulloch et al. 2014). The eddy decorrelation time Td, where the autocorrelation R settles to zero, indicates the appropriate length of time needed after particle release before variance-based methods of assessing diffusivity are valid (LaCasce 2008; Wolfram et al. 2015). However, “full” integration via the integral diffusivity metric (LaCasce 2008; LaCasce et al. 2014),
may result in year-long integration times for real-world shorelines and bathymetry (Sallee et al. 2011) but are shorter and on the order of a month for an idealized zonal channel flow that is faster than the real world because it lacks form stresses (e.g., Abernathey et al. 2013). Assessment using year-long time scales is prohibitive because only a few realizations may be obtained over a decade of simulation. The decorrelation time scale was assessed using a 10-km grid with simulation over a half year, as shown in Fig. 3, and particle advection over 25 days was found to resolve the cluster decorrelation time scale while also yielding a large number of realizations that maintain reasonable spatial coherency (LaCasce et al. 2014; Wolfram et al. 2015). Diffusivity results for the 5-km grid assessed at times of 13, 25, and 50 days were also compared (not shown) to evaluate the sensitivity of this choice of Td and to ensure that reasonable values of diffusivity were obtained. The scaling of to κ considered below was also in agreement with results in Wolfram et al. (2015) ensuring that computed Lagrangian scales and diffusivities were consistent as expected using Td = 25 days.
b. Comparison with MSCL theory
Comparison of diffusivity with respect to MSCL theory requires consideration of the phase speed of nonlinear eddies cw used in (1). The phase speed is estimated from the linear baroclinic Rossby wave speed Doppler shifted by the depth- and time-averaged zonal-mean flow (Klocker and Abernathey 2014; Klocker and Marshall 2014):
where is the depth- and time-averaged zonal-mean flow, LD is the eigenvalue problem–computed Rossby radius of deformation (Chelton et al. 1998; Ferrari et al. 2010; Wolfram et al. 2015), β is the meridional gradient of the Coriolis parameter f, and the term involving Uyy = ∂2U/∂y2 is due to the modification of the eddy phase speed due to shear (Srinivasan and Young 2014; Klocker et al. 2016; A. Klocker et al. 2016, unpublished manuscript). Note, however, that in high latitudes β and LD are small, so cw is largely governed by (Klocker and Marshall 2014) and/or shear. However, computed Uyy for the idealized circumpolar current (not shown) are O(10−12) m s−1, an order of magnitude smaller than β = 10−11 m s−1. Hence, the critical layer that occurs where corresponds to the vertical location zc, such that .
For purposes of comparison with diagnostics in buoyancy space the potential density at the critical layer is computed using the inverse mapping between the time-averaged layer z coordinate and the time-averaged potential density. However, even independent of this mapping, the location of the critical layer based on the phase speed is an approximation that employs assumptions of linearity. In reality, the phase speed is much more complex and is formed from a relationship between β, , and the wavenumbers (Vallis 2006).
c. Decomposing the Eulerian velocity into high-pass and low-pass components
Online time filtering is employed to decompose the Eulerian velocity field into high- and low-pass components, corresponding to a decomposition of the velocity field into eddy and mean components, respectively. Typically, offline decompositions employ a running mean. However, use of a running mean is not possible for online computations. Consequently, a recursive filter is employed to decompose the velocity into its mean and eddy components. Derivation of the single-pole, recursive, low-pass exponential impulse filter employed in this study is presented in appendix A. A filter time constant of τ = 90 days, which is used to filter out time scales occurring at or faster than the decorrelation time scale of O(25) days, is chosen as the time scale for separating the low-frequency, background flow from the high-frequency eddying flow. Success of the filtering approach is indicated by the full and high-pass filtered flow having similar structure and magnitudes of urms (not shown), while the high-pass filtered flow also exhibits a zonal flow near zero (not shown). Note that the full and filtered Eulerian velocity fields are used to compute Lagrangian particle trajectories that are directly associated with these velocity fields.
d. Diffusivity decomposition from the Eulerian velocity decomposition
Enhanced diffusivity results from interactions between the mean flow and eddies that are obscured in the calculation of diffusivity derived from particle integrations of the full velocity field alone because it combines the effect of slow and fast time scales and their interactions:
where xFULL is a particle trajectory, uLOW is the mean velocity computed with the low-pass velocity, and uHIGH is the eddy component of the velocity computed with the high-pass velocity such that
Note that uLOW, which contains slow time scales, and uHIGH, which contains fast time scales, are both evaluated at the location xFULL that evolves due to uFULL, indicating interactions between the slow and fast velocity time scales.
All particles are advected from the same initial condition x(t = 0) = x0. Particle trajectories by the mean Eulerian velocity field are computed using
and advection by the eddy Eulerian velocity field is via
Particle trajectory interactions are heuristically quantified via the metric
The utility of (15) deteriorates at long time scales past the decorrelation time scale but provides some descriptive capability to understand the diffusivity decomposition at shorter time scales.
This decomposition approach allows direct comparison to MSCL theory and is distinct from the generalized Lagrangian mean (GLM) theory (Andrews and McIntyre 1978; Craik 1985; Bühler 2014, chapter 10) because the Eulerian velocity field is first filtered and then used to produce Lagrangian trajectories, whereas in GLM the unfiltered Eulerian velocity field is used to produce Lagrangian trajectories, that is, xFULL, which are then decomposed into their slow- and fast-varying components.
For notational simplicity, let the meridional diffusivity κyy be designated as κ with subscripts designating components of decomposed diffusivity hereafter. Diffusivity is computed from Lagrangian particle tracks derived via direct integration of velocity fields, for example, κFULL from (11), κLOW from (13), and κHIGH from (14). Diffusivity computed from the xDIFFX particle positions using (15) is designated by the dispersion-based diffusivity κDIFFX. The residual diffusivity κDIFFU is defined as
where κDIFFU is the residual diffusivity inspired by the linear velocity decomposition in (12).
Thus, interactions of the mean flow with eddy variability can be quantified by computing diffusivity for six cases: 1) a full diffusivity κFULL computed from Lagrangian particle trajectories advected by the full Eulerian velocity field in (11); 2) a low-pass or mean flow diffusivity κLOW for particle trajectories advected by the low-pass filtered Eulerian velocity field in (13); 3) a high-pass or eddy diffusivity κHIGH for particle trajectories advected by the high-pass filtered Eulerian velocity field in (14); 4) a residual diffusivity κDIFFU based on decomposition via (16); 5) the dispersion-based diffusivity corresponding to residual particle positions given by (15), κDIFFX; and 6) an estimate for the unsuppressed diffusivity case UNSPR, κUNSPR (see below). Note that a pure shear flow does not have a true diffusivity associated with it because particle motion is never decorrelated. However, for the purposes of our decomposition, we estimate its potential contribution as a “diffusivity” even though κLOW does not necessarily represent a physical diffusivity.
In a model that resolves the mean flow but not the eddies, the diffusivity derived from the resolved flow κFULL should be implemented, which is approximately composed of the sum of the eddy κHIGH and residual κDIFFU diffusivities. Because eddies are not represented in the model, the effect of eddies and the residual, κHIGH and κDIFFU, respectively, must be parameterized to compute the full κFULL.
The κDIFFU and κDIFFX estimates of mixing account for the asymptotic cases of homogeneous turbulence (e.g., high-pass filtered velocity) as well as mean meridional shear flow mixing (e.g., low-pass filtered velocity) because κDIFFU = 0 and xDIFFX = 0 for either uFULL = uHIGH or uFULL = uLOW. As explained in appendix B, the difference κDIFFX − κDIFFU provides physical insights into diffusivity suppression.
e. Estimation of unsuppressed diffusivity UNSPR
To place the components of the decomposed diffusivity in context, we estimate a range of unsuppressed diffusivities (case UNSPR) using a lower bound, which is constrained to the peak diffusivity observed near the critical layer location for the FULL diffusivity, and an upper bound Rhines scale estimate (Klocker and Abernathey 2014; Klocker et al. 2016).
The unsuppressed diffusivity lower bound estimate is computed as
where the effective mixing length is based on the location of the diffusivity maxima that is minimally suppressed, namely,
where the subscript MCL denotes evaluation of Lmix,eff at the location of the maximum at the critical layer.
The unsuppressed diffusivity upper bound can be estimated from (2) using the Taylor (1915) mixing efficiency Γ = 0.5 and an upper bound mixing length based on the Rhines scale (Rhines 1975; Klocker and Abernathey 2014; Klocker et al. 2016):
which delineates the boundary where the Rossby restoring force dominates, resulting in a regime dominated by Rossby waves instead of geostrophic turbulence (Sukoriansky et al. 2007; Klocker et al. 2016). The Rhines scale mixing length is approximately 30% larger than observed eddy scales (Tulloch et al. 2011; Klocker et al. 2016). An additional estimate for unsuppressed diffusivity, obtained with corresponding to observations from Klocker et al. (2016) for the strongly nonlinear regime (r ≈ 80) can also be made. Note, however, that the present model does not account for real-world form stresses due to variable bathymetry, and consequently this estimate may be too large.
a. Hydrodynamics of the ICC simulations
An overview of the zonally periodic ICC flow is shown in isometric view in Fig. 1. The plot is split into two sections where isosurfaces for the mean zonal velocity are plotted for 0 < x < 500 km. The resultant eastward-propagating ICC that approximates an Antarctic Circumpolar Current can be seen in red, approximately centered at y = 1500 km, and the westward-propagating slope front current (SFC) that approximates an Antarctic Slope Front is over the shelf break at y = 500 km.
Instantaneous and mean kinetic energy are plotted for 500 < x < 1000 km. The instantaneous kinetic energy, shown in green, reveals the presence of strong eddies, and the mean kinetic energy, shown in transparent gold, is largely coherent with the mean flow. As shown, these coherent structures are of a scale of O(100) km. The mean kinetic energy and mean flow are fairly zonally homogeneous. Additional details related to the ICC’s hydrodynamics are discussed by Ringler et al. (2017).
1) Lagrangian trajectories
Lagrangian pathlines and mean fluid speed over 7 days are shown in Fig. 4 for a single realization for a single buoyancy surface. This view of the flow demonstrates the strong eddy and mean eastward flow in the ICC as well as the smaller eddy and westward mean flow over the shelf break. Smaller eddies have turnover times on the order of a week, which correspond to the pathlength shown in the figure (e.g., in the ICC at approximately x = 500 km and y = 1700 km). Lagrangian paths in the ICC typically encounter large eddies but are not typically entrained within a coherent eddy, as evident by only a few Lagrangian trajectories advecting westward in the ICC. Instead, the majority of trajectories propagate toward the east with their paths curved by the superimposed eddy and mean flow. In contrast, eddies near the shelf break and on the shelf are much smaller than in the deep part of the channel. These trajectories lack strong, coherent zonal advection in the eastward direction. Particles within the SFC are advected toward the west, and particles on the shelf are advected primarily by eddies with an indistinct zonal propagation in regions away from the SFC. Although useful to qualitatively describe the flow, the unprocessed particle paths are unable to directly quantify mixing within ICC. To this end, the Lagrangian mean and eddy flows are quantified using cluster statistics developed from clustering of individual Lagrangian trajectories in space and over realizations to compute spatial structure (Wolfram et al. 2015) for the FULL case.
2) Zonal-mean and eddy velocity
The zonal-mean flow in Fig. 5 is characterized by a strong 0.25 m s−1 zonal eastward-flowing transport within the ICC (e.g., at y = 1500 km) with a smaller westward mean flow of 0.05 m s−1 over the shelf break within the SFC (e.g., at y = 500 km). Within the ICC the flow is coherent throughout the water column with the strongest velocity at the surface. The SFC, however, is strongest at depth at the shelf break at y = 500 km with very little mean flow on the shelf.
Like the zonal-mean velocity, the eddy velocity is greatest within the core of the ICC, as shown in Fig. 6. The spatial structure of the eddy velocity is similar to the zonal-mean flow in Fig. 5 with surface intensification in the core of the ICC at y = 1500 km with a fairly coherent vertical structure. The eddy velocity, however, attenuates more rapidly in the vertical than does the zonal-mean flow. The eddy velocity is smallest directly over the shelf break at y = 500 km and low eddy velocities occur for 0 ≤ y ≤ 1000 km. This result is consistent with the Lagrangian trajectories presented in Fig. 4, where Lagrangian trajectories are dominated by large eddies and large mean flows within the ICC from 1000 ≤ y ≤ 2000 km and by smaller eddies otherwise.
3) Dominant scales
The Rossby radius of deformation, shown in Fig. 7, is largely governed by the mean stratification and water depth. It ranges from approximately 20 km in the ICC to less than 5 km on the shelf break and onshore. Its associated linear Rossby wave speed scaling is O(10−3 to 10−5) m s−1. The estimated Doppler-shifted eddy phase speed, also shown in Fig. 7, is approximately 0.25 m s−1 within the ICC on account of the time- and depth-averaged zonal-mean flow. The critical layer for the ICC occurs at an approximate depth of 1 km on the 1029.5 kg m−3 potential density surface.
The meridional flow is approximately decorrelated after 20 days (Fig. 3), and we compute meridional diffusivities at t = 25 days. Spatial structure and magnitude for the diffusivity and its decomposed components are not particularly sensitive to this choice, for example, diffusivities at 13, 25, and 50 days are comparable in terms of magnitude and structure (not shown). However, spatial fidelity is lost at longer assessment times because particles diverge from their initial cluster locations so that a shorter time scale is preferable (Wolfram et al. 2015). The integral time scale is biased by the presence of a negative lobe in the autocorrelation (Abernathey et al. 2013) that produces a peak value at approximately 2 to 10 days. The integral time scale ranges from hours to several days (not shown), and the decorrelation time is consequently much larger than the integral time scale. Furthermore, the decorrelation scale of several weeks is also qualitatively consistent with Fig. 4, which shows particles traversing several fractions of a complete orbit, for example, near x = 600 km and y = 1700 km, over a period of 7 days.
b. Full diffusivity
Figure 8 shows the spatial structure of the complete diffusivity κFULL that is computed based on Lagrangian particles advected with the full Eulerian velocity field uFULL. Diffusivity values are plotted on a log10 scale with contours at 0.3 intervals corresponding to doubling of diffusivity magnitude. All spatial plots for diffusivity are presented with identical color maps and contour interval.
Mixing by the FULL flow in Fig. 8 is greatest within the ICC, reaching values of 6000 m2 s−1, and is weakest at the shelf break and onshore with values O(100) m2 s−1, an observation that is qualitatively corroborated by the scale of eddies represented by pathlines for the realization in Fig. 4. The diffusivity maximum occurs on a subsurface potential density surface corresponding to approximately the 1029.4 kg m−3 surface. The spatial location of the diffusivity maximum is in the vicinity of the estimated buoyancy surface for the critical layer and near the region of maximum wind stress in the ICC. The role of the eddies and the mean flow in contributing to the mixing within this region can subsequently be explored by considering diffusivities derived from particle trajectories from the low- and high-pass filtered Eulerian velocities, that is, (13) and (14).
c. Diffusivity from the low-pass filtered Eulerian velocity
Mixing by the LOW flow, as approximated by the low-pass Eulerian velocity filter, results in the κLOW shown in Fig. 9a. Values are one to two orders of magnitude smaller than for κFULL and are largely potential density independent for each meridional location. This result is to be expected because mixing by the LOW flow is largely due to the low-pass filter extracting the mean flow, which is predominantly characterized by meridionally varying zonal shear (e.g., Fig. 5) that does not enhance cross-stream mixing, for example, κyy, as it does zonal along-stream mixing κxx via classic shear dispersion (Taylor 1953; Fischer et al. 1979; Oh et al. 2000; Griesel et al. 2010, 2014), but this mixing enhancement is not considered here.
d. Diffusivity from the high-pass filtered Eulerian velocity
Figure 9b quantifies κHIGH, based on trajectories using the high-pass Eulerian filtered velocity field that is analogous to the standard Eulerian eddy velocity computed by removal of the mean. Diffusivities are highest at the surface, corresponding to the greatest eddy kinetic energy with attenuation toward depth. The slope of the diffusivity contours in buoyancy space at the ICC is aligned with the gradient of the nonlinear wave phase speed with values increasing toward the north. It is modulated by a mixing length that is dependent upon the Rossby radius of deformation LD, which decreases toward the shelf and is the dominant spatial scale for the eddies. Similar to κFULL in Fig. 8, diffusivity is small over the shelf break. Characteristic κHIGH values are more than a factor of 5 smaller than for κFULL diagnosed using the FULL velocity, which is a surprising finding in light of other studies (Klocker et al. 2012b; Abernathey and Marshall 2013) that found eddy-produced diffusivities, for example, κHIGH, to be larger. The relevance and consequence of this interesting result is explained in the discussion section.
e. Residual diffusivity
The residual diffusivity resulting from the difference between the full and combined mean and eddy flows is presented in Fig. 9c. Comparison of the residual κDIFFU in Fig. 9c to the FULL velocity in Fig. 8 suggests that κDIFFU ≈ κFULL. Note that the operation in (16) is applied prior to zonal averaging and consequently values within the figures may not strictly satisfy (16). This diffusivity measures the mixing due to the turbulent combination of eddies and the background mean. The maximum diffusivity occurs near the critical layer location within the ICC at y = 1500 km, and the diffusivity is of similar magnitude to that for the FULL flow. Diffusivity for outcropping portions of the layer in the southernmost portion of layers is very small.
f. Dispersion-based diffusivity
The dispersion-based κDIFFX is obtained by computing the diffusivity associated with the residual particle positions in (15) and is shown in Fig. 9d. It is a factor of 3 times larger in the core of the ICC at the surface relative to the critical layer and the MSCL-predicted maximum diffusivity at depth is absent. It is larger than κFULL and is intensified in the vicinity of the critical layer by 25% relative to κDIFFU.
a. Characteristics of decomposed diffusivities
The subsurface maximum diffusivity predicted by MSCL theory is evident for κFULL and residual κDIFFU in Figs. 8 and 9c. However, the maximum at the critical layer depth does not occur for κLOW, κHIGH, or the dispersion-based κDIFFX, that is, Figs. 9a, 9b, and 9d, respectively. Furthermore, the magnitude of mixing for κLOW and κHIGH in Figs. 9a and 9b is much smaller than for κFULL in Fig. 8, for example, more than 4 times smaller near the critical layer. The implication is that the mean flow effect on diffusivity is particularly important for mixing as κDIFFX is much larger than κFULL and κDIFFU and is qualitatively very similar in magnitude to κFULL, although the spatial extent of the maximum diffusivity is reduced.
2) Spatial structure
The diffusivity for κLOW in Fig. 9a is largely vertically homogeneous across potential density surfaces on account of the depth coherency of the mean flow shown in Fig. 5 within interior layers of the fluid. In contrast, κHIGH in Fig. 9b is primarily depth attenuated, corresponding to decreased eddy velocities at depth. Diffusivity for outcropping layers is largely driven by eddy activity with minimal contributions by the mean flow, as demonstrated by similar magnitudes and structure for the southernmost portion of layers in Figs. 8, 9a, and 9b. The residual κDIFFU in Fig. 9c displays the subsurface maximum as predicted by MSCL theory because it arises due to a suppressed mixing length for κFULL. In contrast, the dispersion-based κDIFFX in Fig. 9d shows a monotonic decrease from surface-intensified diffusivity corresponding to the vertical eddy kinetic energy maximum. Additionally, the spatial structure of κDIFFX is reminiscent of the eddy and mean zonal flow, for example, compare Fig. 9d to Figs. 5 and 6.
3) Mean flow effects on diffusivity
The κDIFFU mixing within the ICC contributes to approximately 80% of κFULL (Fig. 9c). Furthermore, the location of the maximum diffusivity κDIFFU is more closely aligned with the estimated location of the critical layer, shown by the thick white line, than for κFULL in Fig. 8. These observations indicate that mixing within the ICC is strongly governed by HIGH and LOW flow interactions with a subsurface maxima obtained via minimized suppression for the FULL and DIFFU cases. Additionally, the strong dispersion-based κDIFFX in Fig. 9d, particularly at the surface, indicates that the mean flow effect on diffusivity is to constrain the total amount of mixing that could occur because κFULL < κDIFFX, suggesting that diffusivity suppression serves to both constrain the mixing near the ocean surface and minimally reduce mixing for subsurface isopycnals near the critical layer, as predicted by MSCL theory.
b. Mixing dynamics of an inhomogeneous eddying flow
The ICC hydrodynamics within the ICC region at y = 1500 km can be conceptualized as a flow composed of quasi-independent eddies and a background mean flow. To illustrate the role of the turbulent eddy and mean flow in enhancing mixing, consider the clusters of particles that have evolved 25 days from their initial condition of 100-km clusters distributed at y = 500, 1000, and 1500 km, shown in Fig. 10. For the case of the FULL flow in Fig. 10a, with particle positions xFULL computed via (11), eddies serve to produce the initial filamentation of the clusters, and then the background mean flow further serves to elongate these filaments, resulting in advection from one eddy to the next. The implication is that mixing driven by stretching and straining by the combined eddy and mean flow is large, even for the critical layer where . However, shear contributions to diffusivity suppression, understood via MSCL theory, are minimal, as suggested from the negligible modification to cw via reduction of β by Uyy (Srinivasan and Young 2014; Klocker et al. 2016; A. Klocker et al. 2016, unpublished manuscript) based on the small size of Uyy relative to unmodified β presented in section 4b.
Particle dispersion is smallest for the case computed with only the LOW flow via (13) in Fig. 10b, xLOW, because the background flow with its meridional shear is not strong enough, in and of itself, to efficiently spread the particle cluster. The LOW case can be thought of as producing sheared mixing that occurs due to the mean flow. If the background flow were constant, no relative dispersion would occur because each point is advected precisely the same amount from one time to another. In this extreme case no real mixing occurs outside that caused by background diffusion because there is no filamentation to amplify background gradients to instigate downgradient mixing. However, for the HIGH case computed by (14) in Fig. 10c, xHIGH, the initial filamentation produced by the eddies is localized and not dispersed past the mixing length scale for the FULL flow. Mixing consequently occurs much more slowly than for the FULL case in Fig. 10a.
In MSCL theory the mean flow and eddy phase speed reduces the maximum diffusivity, for example, the denominator in (1). Consequently, κDIFFX, computed with (15), is designed to circumvent mixing suppression by reversing the effect of the individual contributions of the LOW and HIGH transport on the FULL. Conceptually, xDIFFX will estimate motion occurring due to residual motion following removal of the filtered eddy and mean trajectories. This residual motion strongly enhances dispersion through the combined action of eddies producing initial filamentation and further stretching and straining of these filaments via the background mean flow as demonstrated in the FULL flow of Fig. 10a and even more so for the DIFFX case Fig. 10d for xDIFFX. Dispersion-based κDIFFX computed from particles in Fig. 9d is even larger than κDIFFU computed with (16) in Fig. 9c, and κDIFFU arises from a decomposition due to differences in particle dispersion, whereas κDIFFX is derived from a decomposition of particle positions that are used to derive dispersion and diffusivity. The implication is that mixing driven by stretching and straining is larger than mixing due to the combined eddy and mean flow, even for the critical layer where . In the absence of the eddy or mean flow component, this mixing is maximized relative to the FULL flow. However, the full flow containing the dynamically consistent combined mean and eddy flow appears to constrain the amount of mixing that could occur as indicated by κDIFFX, although this reduction is minimized in the vicinity of the critical layer denoted by the white line (e.g., compare with Fig. 9c). The diffusivity decompositions demonstrate that κDIFFU is a key contributor to the subsurface diffusivity maxima in the critical layer. Regardless, following the removal of the κLOW and κHIGH contributions from κFULL to produce the residual κDIFFU, κDIFFX is still larger relative to κDIFFU, reinforcing the notion that κDIFFX estimates unsuppressed mixing.
c. Vertical structure of decomposed diffusivities in the baroclinic jet
Vertical profiles in y–z space are used to better understand the diffusivity decomposition and evaluate the applicability of scalings for parameterization. The profiles in Fig. 11 are for meridional diffusivities that are averaged over the range 1400 < y < 1600 km and then projected from buoyancy space back to depths via mean buoyancy layer depths, for example, as shown in Fig. 2. The approximate location of the critical layer is also shown for reference. A scaling estimate for the unsuppressed diffusivity is included in Fig. 11a, where κDIFFX is observed to approximate κUNSPR near the critical layer and away from the free surface. Unsuppressed diffusivity computed using is plotted via a gray line superimposed on the gray shaded region, denoting the range of the estimated unsuppressed diffusivity κUNSPR. The lower bound unsuppressed diffusivity corresponds to , and the higher bound is based on the Rhines scale.
Hatching denotes the z positions inside the ventilation-defined surface layer. As discussed in more depth in Ringler et al. (2017), the ventilation-defined surface layer is the set of buoyancy coordinates that ever reach the ocean surface. We note that the κDIFFU, κFULL, and κDIFFX diffusivities all show a local minimum at the base of the ventilation-defined surface. As described in Ringler et al. (2017), the base of the ventilation-defined surface layer corresponds to an unusual location within the ICC where mesoscale eddies are fluxing Ertel potential vorticity up the mean potential vorticity gradient.
In particular in Fig. 11a, note that the diffusivity κLOW corresponds to mixing due to shear and is very small because there is no turbulent motion to enhance mixing. In contrast, κHIGH, which corresponds to mixing by temporal high-frequency variability of the Eulerian velocity field, for example, eddies, is surprisingly much smaller than κFULL and requires further analysis.
1) Eddy diffusivity departure from previous studies
The κHIGH effectively quantifies the mixing occurring due to the eddy velocity resulting from the removal of the Eulerian time-mean velocity from the full velocity. However, previous researchers (Griesel et al. 2010, 2014; Klocker et al. 2012b; Abernathey and Marshall 2013) noted that the eddy-produced diffusivity was larger than the κFULL that corresponds to mixing suppression by the mean flow. Our results, for example, Figs. 8 and 9b, are at odds with these previous findings.
There are two primary reasons for this departure from the previous literature. First, eddy results in Klocker et al. (2012b) and Abernathey and Marshall (2013) are not generated by removal of the mean flow from a full velocity field. Instead, the results are produced using an AVISO-generated eddy velocity field where the full flow is constructed by addition of an independent mean flow, which subsequently results in mixing suppression in accordance with MSCL theory. This approach is the inverse of this study where the eddy flow is constructed by the removal of the mean from the full flow. Second, results in Griesel et al. (2010, 2014) are obtained performing the eddy decomposition in a similar fashion to this study via construction by removal of the mean from a full flow, albeit offline. They found that the full flow had the largest along-stream dispersion resulting in unconvergent along-stream diffusivities (Griesel et al. 2010, 2014) due to residual shear dispersion in the binning method (Bauer et al. 1998; Oh et al. 2000; Griesel et al. 2010; Koszalka et al. 2011). To the best of the authors’ knowledge they did not report results for the cross-stream meridional diffusivity for the case of the eddy-only flow, as considered here. Hence, a direct comparison to their work is not possible.
Our approach is distinct from both these types of analyses because it decomposes the full Eulerian velocity to obtain the Eulerian eddy velocity and does not consider the along-stream diffusivity. Fundamentally, the mixing that we are computing is derived from the combined action of both the mean and eddy components of the flow field, working together, to produce mixing. The addition of a velocity field, for example, to remove the mean flow, results in Doppler shifting of the eddies, thereby altering the nonlinear parameter r. This is predictable, as demonstrated in the next subsection by the high accuracy of the curve fit for the HIGH pass diffusivity. This approach demonstrates the cause of the significant reduction of the full diffusivity to the eddy diffusivity quantified via the HIGH case.
2) Shifted eddy phase speeds for the HIGH case
A feature of the velocity decomposition for uHIGH in (12) that is not immediately obvious is that the eddy phase speed cw is artificially shifted corresponding to removal of uLOW via high-pass Eulerian filtering. This effect is quantifiable via the nonlinear parameter r (Theiss 2004; Chelton et al. 2011; Klocker and Abernathey 2014; Klocker et al. 2016) that relates the eddy velocity urms to the non-Doppler-shifted eddy phase speed that accounts for the thermal component of the mean flow
where utherm is the mean zonal velocity obtained from the thermal wind relation corresponding to baroclinic shear and meridional mean flow shear effects are ignored, and is the long Rossby wave speed of westward eddy propagation (Stone 1978; Held and Larichev 1996; Klocker et al. 2016; A. Klocker 2016, personal communication). Thus, the mean flow is composed of components due to the depth-averaged and thermal wind component of the mean; for example, for simplicity, define for calculations used in this paper, where is the mean produced by low-pass filtering. Filtering results in an artificial Doppler shifting arising due to the velocity decomposition analogous to a Galilean transformation. The relative velocity of the eddy has changed in the Eulerian frame of reference following the removal of the mean. The eddy phase speed is not simply but also includes effects due to the depth-mean flow, that is, it is already Doppler shifted (Klocker and Marshall 2014). Equation (20) ignores the effects of the background barotropic mean flow by assuming the “non-Doppler effect” (Held 1983; Colin de Verdière and Tailleux 2005; Samelson 2010; Klocker et al. 2016). To better understand the modification of the eddy phase speed, consider a simplified flow of us = use + u0, where u0 corresponds to a constant advective flow and use is an eddy field. The term us has eddies that are propagating at some speed that includes the depth-mean Doppler shift by u0 (Klocker and Marshall 2014). Thus, us has eddies moving at in the absolute Eulerian frame of reference. Note that on its own u0 does not have an associated eddy propagation speed. Thus, the relative propagation speed of eddies under high-pass filtering that removes u0 is (us)HIGH = us − u0. The eddy phase speed is necessarily the Doppler-shifted value , where the filtering in this case is effectively a kinematic transformation that produces a Doppler shifting of the velocity field us that has absolute eddy propagation speed . In this simplified case the high-pass filtering acts like a Galilean transformation that essentially changes the frame of reference to modify the absolute speed of eddy propagation in (us)HIGH.
This phase speed Doppler shifting results in translation in the denominator of (20) by uLOW, namely,
where i is the unit vector in the zonal direction. A reasonable approximation for the ICC flow is because is small at this latitude relative to the depth-mean flow. Thus, the high-pass filtering step both retains the eddy velocity urms as well as Doppler shifts the effective eddy phase speed that ultimately produces mixing via eddy nonlinearity. This shifting is a numerical consequence of filtering because the eddy field is computed for the full flow by the primitive equations, and subsequent eddy nonlinearity is anchored to the Doppler shifting that occurs in the full flow. High-pass filtering of the Eulerian velocity field, consequently, results in a Doppler-shifting that modifies r.
Mixing consequently transitions from arising due to a geostrophic turbulence regime dominated by nonlinear eddies to a regime characterized by linear wavelike behavior (Klocker and Abernathey 2014; Klocker et al. 2016) under high-pass filtering depending upon the relative magnitudes of the eddy and mean flow velocities. Consequently, (14) does not correspond to motion for an unsuppressed velocity field but instead corresponds to motion for a wave-dominated regime if , yielding (Klocker and Abernathey 2014; Klocker et al. 2016). Thus, decomposed κHIGH is unable to provide information about the true unsuppressed diffusivity because the absolute eddy phase speed is shifted in the temporal decomposition of velocity. The HIGH case consequently corresponds to a different diffusivity regime than the FULL because it is in the wave-dominated regime , not the nonlinear regime r ≫ 1. The resultant κHIGH is consequently different than what would be expected for the true unsuppressed eddy diffusivity corresponding to turbulence without a mean flow for r ≫ 1, for example, as observed in other studies (Klocker et al. 2012b; Abernathey and Marshall 2013), and is more representative of the mixing occurring due to Rossby waves than turbulent mixing (Klocker and Abernathey 2014; Klocker et al. 2016). Removal of the mean velocity modifies r via (21), resulting in modification of the unsuppressed diffusivity, which is indirectly set by r via modification of Lmix in (2). The mean flow has been removed so that K = ΓurmsLmix. Consequently, κHIGH does not demonstrate a subsurface maxima as do κDIFFU or κFULL.
Thus, the diffusivity is set by a different mixing length scale relative to the FULL case, although the eddy velocities urms are approximately the same by construction via the high-pass filtering. Using (21) and the following relationship derived from Fig. 1 for the Southern Ocean from Klocker et al. (2016),
the diffusivity can be estimated with LD = 18.5 km and Γ = 0.35 (Klocker and Abernathey 2014). The resultant fit using the depth variable Lmix from (22) in Fig. 12 results in a Pearson correlation coefficient of 0.985 with a p value of 1.4 × 10−6 for the fit with κHIGH, where a value of 1 indicates perfect agreement, and the p value indicates the probability that random data would produce a fit of similar or higher quality (Crow et al. 1955). The close agreement between this theoretical estimate and the computed κHIGH demonstrates the applicability of the interpretation that nonlinearity in κLOW is reduced via a shift in r toward the linear regime that is typically characterized by linear Rossby waves (Klocker and Abernathey 2014; Klocker et al. 2016). The high-pass diffusivity exhibits a dependence on r and LD, with associated Lmix, as well as the vertical structure of urms.
3) Physical interpretation of diffusivity decomposition
The decomposed diffusivity can now be described in terms of its physical significance. The κLOW corresponds to mixing occurring due to the mean flow, which in this case is dominated by shear as indicated by Fig. 10b. As demonstrated in the previous section, κHIGH corresponds to mixing due to reduced nonlinearity of the eddies, for example, as would be expected for linear Rossby waves (Klocker and Abernathey 2014; Klocker et al. 2016). The term κDIFFU quantifies the turbulent diffusivity produced from decomposition of the FULL (suppressed), LOW (shear), and HIGH (linear wave) diffusivities via (16), assuming that for the ICC flow where the mean velocity is as strong as the eddy velocity and the mixing becomes dominated by waves for the high-pass filtered flow. Note that a fundamental challenge in the diffusivity decomposition is that a clear method to evaluate the unsuppressed diffusivity from the model output is at present unavailable. However, κDIFF, by construction, estimates unsuppressed mixing and mimics κUNSPR as estimated from scaling. The estimates diverge near the surface layer but are coincident at depth.
d. Quantifying eddy and mean flow interaction reduction to turbulent mixing
The role of the eddies and mean flow in constraining mixing can be understood by considering the reduction of turbulent mixing computed via the difference between the estimated unsuppressed turbulent mixing κDIFFX and the suppressed turbulence κDIFFU. As shown in appendix B, assuming a zero κLOW contribution and corresponding to very slow mean meridional advection relative to dispersion, the reduced mixing ΔκNyy = κDIFFX − κDIFFU in the meridional direction is
For ΔκNyy > 0, as observed in Fig. 11a, this implies that linear waves in the eddy field represented in the HIGH case, for example, Rossby waves (Klocker and Abernathey 2014; Klocker et al. 2016), are a restoring force that displaces particles in the opposite direction relative to turbulent mixing produced by baroclinic instability. The waves consequently constrain the maximum diffusivity that could contribute to mixing. This is conceptually consistent with the observation that linear Rossby waves oscillate particles meridionally along lines of constant latitude, whereas turbulent mixing works to disperse particles away from lines of constant latitude. Linear Rossby waves can thus serve to constrain mixing by restoring particles to lines of constant latitude in opposition to the unsuppressed turbulence.
e. Diffusivity parameterization
A generalized description of the vertical structure of diffusivity, including the location of the subsurface maxima, is well described by generalized MSCL theory via (1). The use of (1) in parameterizing the structure of diffusivity can be analyzed following Bates et al. (2014) via fitting computed diffusivities to the model in Eqs. (1), (2), and (3). Resultant nonlinear fits via a constrained nonlinear least squares algorithm in Fig. 11b for κFULL use the available structure of urms(z) and perform fits for constant Lrms and TL, respectively. The fits subsequently indicate that the structure of the diffusivity near the critical layer maximum can be well-represented in terms of both magnitude and shape by MSCL theory. Additionally, we also consider Taylor’s classic scaling (from Wolfram et al. 2015) that is analogous to (3) but with the depth variable TL(z) computed using the integral time scale and constant parameter Γ = α, namely,
where R is from (7) and is spatially averaged to yield vertical profiles as in Fig. 11 with potential densities ρ mapped to their associated mean vertical locations z. The integral in (25) is evaluated past Td to approximate the infinite integral. Use of (24) results in a fitted α = 0.47, as designated by the red dashed line in Fig. 11b. This was the best fit obtained and is very close to the theoretical Taylor (1915) value of 0.5 and well within the range computed for an idealized eddying double gyre using cluster diffusivities (Wolfram et al. 2015). This is not necessarily an easily accessible parameterization, however, because the vertical structure of (25) is largely distinct from readily obtainable quantities like urms (Bates et al. 2014).
Although urms and its vertical structure is relatively easy to approximate or compute in a resolved simulation, the structure of Lmix and TL must be parameterized in a more nuanced way. This is evident by examination of two other fits using MSCL theory in (1). One uses vertically varying urms and an assumed vertically constant Lmix via (2), as designated by the black dashed line. Another fit uses the vertical structure of and an assumed vertically constant TL via (3), designated by the black dotted line. Both approaches approximately represent the subsurface maxima occurring near to the approximate location of the critical layer designated by the gold line. However, neither approach assuming vertically constant Lmix or TL sufficiently represents κFULL at depth below the mixing layer. The depth-varying structure of urms and ultimately dominate at depth below the critical layer and attenuate too rapidly, producing an underestimate for diffusivity below the critical layer, at least under the assumption of a vertically constant Lmix or TL. Furthermore, the fitted parameter b in (1) was previously estimated by Bates et al. (2014) to be approximately 1 to 4. Assuming Γ = 0.35, the nonlinear fits yield Lmix ≈ 320 km and TL ≈ 25 days and result in very large values of b of 165 and 430, respectively. Part of this discrepancy is due to = O(0.05) cm s−1 versus urms = O(0.5) cm s−1, an order of magnitude larger. If they were both the same size this would correspond to reduction of b by two orders of magnitude because of the square, resulting in values of 1.7 and 4.3, which are comparable to results in Bates et al. (2014). Overall, however, the difference in the vertical structure implies that some additional understanding of the vertical structure of Lmix or TL is needed to provide an appropriate estimate of mixing at depth.
Mixing quantified by meridional diffusivity for the ICC is highly dependent on both the turbulent eddy flow, the mean flow, and their combined mixing as quantified using a residual diffusivity that separates diffusivity associated with the high- and low-pass filtered Eulerian velocity fields from the full diffusivity. Initial filamentation by eddy scales is stretched to larger-scale structures by mean flow straining. Mixing is due to combined stirring by the full velocity. Decomposition of diffusivity into eddy and mean constituents demonstrates that turbulent eddy and mean flow interactions ultimately account for 80% of the diffusivity. Interactions between the mean and eddy flow both enhance and suppress diffusivity, leading to a maximum at depth, as predicted by MSCL theory.
A simple thought experiment may help aid insight into important interactions. Let the mean flow velocity be a constant in the zonal direction and the eddy component be standing eddies, for example, a mean flow perturbed frozen-field case (Lumpkin et al. 2002). The mean motion is simply xMEAN = uMEANt + x0 and the eddy motion orbits within an eddy length scale Leddy of the particle initial position with mixing occurring due to perturbations transporting particles from eddy to eddy. Dispersion for the eddy case is small and dependent upon the background diffusivity providing the perturbation. In contrast, dispersion for the mean case is quadratic in time corresponding to ballistic motion in the zonal direction, zero otherwise. However, dispersion is not maximized under either case but is instead increased when transport occurs under the combined velocity as particles may rapidly advect from one eddy to the next, allowing spreading even in directions normal to the mean flow, assuming the eddies are not perfectly aligned in the zonal direction.
In the idealized simulation the removal of the mean affects the phase speed of eddies and shifts turbulence and its associated mixing from that dominated by nonlinearity to linearity as quantified via r (Klocker et al. 2016). This is directly accounted for in the diffusivity decomposition. To the best of the authors’ knowledge, a technique to recover the unsuppressed eddy velocity field from uFULL that does not alter r via modification of the eddy absolute phase speed, that is, convert (20) to (21), is an open question. Better estimates for the unsuppressed diffusivity would be available given this capability.
The estimated strength of the interaction of the combined eddy and mean to the full diffusivity demonstrates that residual eddy and mean flow interactions dominate mixing with a turbulent contribution κDIFFU accounting for 80% of the diffusivity. Dispersion plots of a realization (Fig. 10a) demonstrate that initial filamentation by eddy scales is stretched to larger-scale structures by mean flow straining. Ultimately though, mixing due to stirring by the full velocity is constrained by restoring due to linear waves arising from Doppler shifting due to the velocity field decomposition ( appendix B). At present, parameterizations by MSCL reasonably capture the effects of mixing within the vicinity of the critical layer but further progress is clearly needed, particularly to better parameterize diffusivity at depth (e.g., Fig. 11b).
The implication of coupled eddy and mean flow processes, for example, via nonhomogeneous turbulence, being responsible for the mixing in baroclinic jets is potentially profound for mixing parameterizations. First, mean flow suppression formulations demonstrating diffusivity amplification at depth appear necessary, as suggested by MSCL theory. But parameterizations based on linear baroclinic growth rates are likely unable to capture the enhancement of diffusivity because the time integral “memory” required to capture the turbulent interactions of eddies with the mean flow is lacking (Sinha and Abernathey 2016). The diffusivity based on the linear baroclinic growth rate (Bretherton 1966; Green 1970; Bates et al. 2014; Griesel et al. 2015) can be represented using (1) and (3) following appropriate choices for b and TL, which can be computed from our data via the nonlinear fit shown in Fig. 11b. Following the fit, the diffusivity is 20% too large near the critical layer but nearly a factor of 4 too small near the channel bottom, reflecting the inability of the theory to appropriately account for the vertical structure of the diffusivity. Third, linear mixing length theories ultimately may not be able account for nonlinear baroclinic growth rates unless Kmax is nonlinear beyond the structured imposed by depth variable urms, with depth variable Lmix, TL, and/or mixing efficiency Γ, that is, in (2) and (3) eddy size and mixing memory arising from interactions with the mean flow appear to be important. For example, this is why Lmix is typically scaled directly from the Rossby radius of deformation but requires additional information for its estimation, for example, evaluation via r as in (22).
Therefore, we ultimately conclude that a better understanding of the vertically variable feedback between the eddy and mean flows is necessary to build better parameterizations for mixing. Progress toward this endeavor is challenging and has been made here. However, a full understanding of parameterization will likely require more advanced decompositions potentially using both Eulerian spatial and temporal filtering. Lagrangian decomposition of xFULL via generalized Lagrangian mean theory (Andrews and McIntyre 1978; Craik 1985; Bühler 2014, chapter 10) may also be fruitful (B. Fox-Kemper 2017, personal communication). The limited applicability of linear parameterizations based on linear baroclinic growth rates or standard MSCL theory implies that parameterizations may require large-eddy simulation approaches to model nonlinear subgrid-scale interactions (Eden and Greatbatch 2008; Marshall and Adcroft 2010), as also noted by Sinha and Abernathey (2016). In general, the vertical structure of diffusivity does not appear to be simply self-similar to the stratification N2 or the first Rossby radius of deformation eigenvector (Bates et al. 2014; Wolfram et al. 2015). Instead, it requires consideration of mixing suppression over depth to represent the vertical diffusivity maxima in the ICC; additional improvements will likely also require some parameterization of the residual diffusivity κDIFFU that modifies the vertical structure of the diffusivity, particularly below the critical layer. These conclusions should be interpreted within the context of our diagnostic’s capabilities. Particle methods typically only track fluid motions. Mixing is a product of the time integral instantaneous straining and stretching, background fluid diffusivity, and the time history of a particular concentration field and its gradients as it is strained over time by a particular Eulerian velocity field. Simplistic particle statistics that lack representation of the concentration gradient cannot therefore directly quantify the time integral because they assume concentration invariance, as demonstrated via the connection made between particle dispersion and effective diffusivity in Klocker et al. (2012b). The concentration invariance assumption may be violated for certain spatiotemporal scales. The shortcoming of diffusivity metrics based on particle dispersion is that reversible and irreversible mixing are indistinguishable.
However, the particle statistics do loosely quantify the maximum mixing that may occur due to straining. Classical techniques (Taylor 1921, 1935) require a decorrelation time scale because the latent assumption, particularly for homogeneous turbulence, is that the straining eventually results in decorrelation that is representable by a diffusion operator or Redi (1982) diffusivity. In principle, this implies that diffusivity metrics relying upon particle statistics provide upper bound estimates for mixing. For example, no mixing occurs for a highly strained fluid if the passive tracer is already fully mixed. It is also possible that the conclusions derived in this paper, namely, that mixing is highly dependent upon mixing enhancement due to turbulent eddy and mean flow interactions, are exaggerated because the computed diffusivity may also be exaggerated relative to transport of a passive tracer in the real ocean. However, consistency between particle, tracer, and effective diffusivity methods has been demonstrated in certain cases (Klocker et al. 2012b; Abernathey et al. 2013), suggesting that the results in this paper are at a minimum a reasonable provisional estimate.
Further diffusivity decompositions, perhaps with Eulerian spatial and/or bandpass temporal filters or via use of generalized Lagrangian mean theory for the decomposition of Lagrangian trajectories (B. Fox-Kemper 2017, personal communication; Andrews and McIntyre 1978; Craik 1985; Bühler 2014, chapter 10), are likely needed to better explain the components of turbulent diffusivity and provide improved estimates for the unsuppressed diffusivity, ultimately working toward a diffusivity closure. Investigations with this focus are likely necessary to provide the key insights needed to develop robust parameterizations of diffusivity arising from mesoscale eddies.
This research was supported by the Office of Science, Office of Biological and Environmental Research of the U.S. Department of Energy Accelerated Climate Model for Energy (ACME) and used computational resources provided by the Los Alamos National Laboratory Institutional Computing facility. Code developments and simulations relied heavily on the work of the MPAS dynamical core development team at LANL and NCAR and in particular the contributions from the MPAS-Ocean development team at LANL. We thank Matt Rocklin for his support of the Dask (Rocklin 2015; Dask Development Team 2016) and distributed Python packages that we use to perform off-node parallel calculations, Stephan Hoyer and Joe Hamman for help with use and development of the Xarray Python package (Hoyer and Hamman 2017), Andrew Stewart for discussions of idealized Southern Ocean physics, Matt Rayson for sharing particle pathline plotting code, Luke van Roekel for contributions to the design of ICC, Mat Maltrud for discussions of mixing theory, Juan Saenz for discussions related to design of idealized Southern Oceans and transport along isopycnal surfaces, Andy Hogg and Ru Chen for their helpful comments on preliminary results leading to this work, Andreas Klocker for sharing preprints of his work, and reviewers Andreas Klocker, Dhruv Balwada, and Baylor Fox-Kemper whose comments greatly improved the quality of this manuscript.
Online Temporal Filtering
Online temporal filtering requires that the estimate for the low-pass filtered quantity at time step n + 1 is based on velocities from known values at time steps n, n − 1, and so on. A straightforward two-level filter can be derived by considering the ordinary differential equation for an impulse function, namely,
where ϕ is quantity to be filtered, ϕL is the low-pass filtered quantity, and τ is the time scale over which an instantaneous impulse is applied. The high-pass filtered quantity ϕH is given by ϕ − ϕL.
A running mean is recovered for the case where the impulse is approximated by a centered-in-time box function. However, for the case of online computations, values at times n + 1, n + 2, and so on, are unknown, and the simplest approximation to the impulse function is a backward-in-time box function resulting in a discrete equation of
that apportions the time-lagged perturbations to the mean over the time scale τ. This filter is equal to standard single-pole filters. Transfer function analysis is given in Smith (1997).
Quantifying Reduction of Nonlinear Turbulent Mixing
The dispersion-based κDIFFX seeks to approximate unsuppressed diffusivity, and the residual κDIFFU quantifies the effect of nonlinear turbulent mixing for flows with . The difference
provides physical insights into the mechanisms responsible for mixing suppression, for example, an estimate of physical processes reducing turbulent mixing. For the ICC flow, κLOW corresponds to mixing by shear, κHIGH corresponds to mixing by waves, κDIFFU corresponds to nonlinear turbulent mixing, and κDIFFX approximates κUNSPR. The following analysis assumes this physical interpretation of the diffusivity decomposition to better understand the role of interactions between different compartments in the flow in modifying the unsuppressed diffusivity estimate κDIFFX relative to the realized nonlinear turbulent diffusivity κDIFFU.
The differential definition of meridional diffusivity based on the meridional displacement y = j⋅x, assuming , is
where the operator and diffusivity can also be computed equivalently (LaCasce et al. 2014) from undifferentiated dispersion, that is,
where t is the time since the particle release. The formulation in (B3) is more amenable to analysis than (B2) and is subsequently used in the following derivation. Note, however, that both formulations produce similar diffusivities (LaCasce et al. 2014).
Expansion and algebraic simplification of (B5) yields
where ΔκNyy is essentially the diffusivity resulting from the variance . If Δκyy > 0, as is the case for the ICC flow considered above, this implies that 1) the LOW and HIGH interaction (corresponding to interactions between shear and linear waves) must be positive and larger than the correlated motions of the DIFFX with the LOW and HIGH (corresponding to unsuppressed interactions with combined shear and linear waves dynamics) and/or 2) the combined HIGH and LOW motion via shear and linear waves directly work to displace particles in the opposite direction relative to the unsuppressed estimate of homogeneous turbulence DIFFX and consequently constrain resultant turbulent nonlinearity (DIFFU) to ensure that , even to the point of superseding a potentially negative corresponding to interactions between linear waves and shear. Although (B6) is physically useful to describe mechanisms for the nonlinearity, its complexity supports multiple physical explanations for ΔκNyy > 0 and further reduction is desirable.
For further clarity and simplification for the ICC flow, assume that yLOW is very small and can be ignored, a reasonable assumption because for the ICC. Then,
If Δκyy > 0, as is the case for the ICC flow considered above, this implies that linear waves in the HIGH case, for example, Rossby waves (Klocker and Abernathey 2014; Klocker et al. 2016), directly work to displace particles in the opposite direction relative to the estimate for unsuppressed turbulent diffusivity κDIFFX and consequently constrain nonlinear turbulent mixing via a reduction of κDIFFX to κDIFFU.