Barotropic eddy fluxes are analyzed through a geometric decomposition of the eddy stress tensor. Specifically, the geometry of the eddy variance ellipse, a two-dimensional visualization of the stress tensor describing the mean eddy shape and tilt, is used to elucidate eddy propagation and eddy feedback on the mean flow. Linear shear and jet profiles are analyzed and theoretical results are compared against fully nonlinear simulations. For flows with zero planetary vorticity gradient, analytic solutions for the eddy ellipse tilt and anisotropy are obtained that provide a direct relationship between the eddy tilt and the phase difference of a normal-mode solution. This allows a straightforward interpretation of the eddy–mean flow interaction in terms of classical stability theory: the initially unstable jet gives rise to eddies that are tilted “against the shear” and extract energy from the mean flow; once the jet stabilizes, eddies become tilted “with the shear” and return their energy to the mean flow. For a nonzero planetary vorticity gradient, ray-tracing theory is used to predict ellipse geometry and its impact on eddy propagation within a jet. An analytic solution for the eddy tilt is found for a Rossby wave on a constant background shear. The ray-tracing results broadly agree with the eddy tilt diagnosed from a fully nonlinear simulation.
The dynamics of large-scale ocean motions are strongly dependent upon the effect of the small-scale turbulent eddy field. The Gent–McWilliams parameterization (Gent and McWilliams 1990; Gent et al. 1995) is now a key ingredient in coarse resolution ocean circulation models (e.g., Fox Kemper et al. 2013) and can be interpreted as modeling the vertical flux of momentum due to eddy form stresses (Greatbatch 1998). This vertical transfer of momentum by eddies plays a fundamental role in the dynamics of the large-scale circulation; for example, it is a leading-order term in the dynamics of the Southern Ocean (e.g., Johnson and Bryden 1989; Danabasoglu et al. 1994).
While horizontal eddy momentum fluxes are less significant from a global perspective, they can play important roles in the dynamics of inertial jets. For example, they influence the dynamics of western boundary currents, where they are instrumental in transferring energy between the mean flow shear and eddies (Waterman and Jayne 2011; Waterman et al. 2011; Waterman and Jayne 2012). Horizontal eddy stress is not captured by the Gent–McWilliams parameterization, and its effects are not typically represented in coarse resolution ocean models beyond the influence of increased explicit or numerical dissipation [see Eden (2010) for an exception]. This paper focuses on the study of the geometric properties of the local horizontal eddy momentum stress, and specifically studies these properties in a series of idealized barotropic shear and jet problems. The geometric framework holds promise for elucidating the important roles of the horizontal eddy momentum fluxes in these systems, as well as for suggesting ingredients for a parameterization of their larger-scale effects.
In the mean barotropic vorticity equation, the eddies influence the mean flow through the divergence of an eddy vorticity flux vector. The vector flux, and hence the influence of the eddies on the mean flow, is conveniently visualized in the usual manner via plots of directed arrows, with the length indicating the local magnitude of the flux and the orientation indicating the local direction. In the corresponding mean momentum equation, the eddies influence the mean flow through the divergence of an eddy momentum stress tensor, whose components are related to the eddy Reynolds stresses, and that is equal to the eddy velocity covariance tensor. This stress can no longer be visualized as directed arrows—instead the natural visualization is as an oriented eddy variance ellipse (Preisendorfer 1988; Wilkin and Morrow 1994; Morrow et al. 1994), with a well-defined eddy ellipse tilt and eccentricity.
Eddy vorticity fluxes are equal to the divergence of an eddy momentum stress tensor, of which the velocity covariance (and hence the eddy Reynolds stresses) forms a part, through the Taylor–Bretherton identity (Taylor 1915; Bretherton 1966b; Plumb 1986). In the barotropic vorticity equation only the double divergence of the barotropic part of the momentum stress tensor can influence the mean dynamics directly. Specific geometric properties of the eddy velocity covariance that can influence the mean dynamics (i.e., which do not vanish under this double divergence) are identified in Waterman and Lilly (2015). Three-dimensional quasigeostrophic generalizations to the geometric eddy variance ellipse view are discussed in Hoskins et al. (1983), Marshall et al. (2012), and Maddison and Marshall (2013).
In the formulation we employ here, the average size of the eddy variance ellipse is proportional to the square root of the local eddy kinetic energy, and similarly to the square root of an average magnitude of the local eddy momentum stress.1 The direction of the major axis indicates the direction in which the corresponding component of the eddy velocity leads to the largest variance (Wilkin and Morrow 1994), and consequently the greatest stress. The ratio of the size of the ellipse in the major and minor axis directions is equal to the square root of the ratio of the variances associated with the components of the eddy velocity in each of these directions, and similarly to the ratio of square root of the magnitude of the eddy stress due to these components. These properties are illustrated in Fig. 1. The ellipse associated with the eddy velocity covariance tensor can be characterized by three key geometric properties: a magnitude, tilt, and eccentricity. Together these capture the complete structure of the local eddy momentum stress.2 Formally these geometric properties arise from a principal component analysis of the eddy velocity covariance tensor (Preisendorfer 1988).
Eddy variance ellipses and their geometric properties are useful diagnostics to characterize the eddy field and its mean flow interactions. For example, eddy variance ellipses are used to compare observational and model eddy variability in Wilkin and Morrow (1994). Southern Ocean satellite altimetry-derived eddy velocity covariance is analyzed using this approach in Morrow et al. (1994). Associations between ellipse geometry and topography, and between anisotropy and the mean flow, are noted. A similar analysis is applied in Scott et al. (2008), and associations between the eddy velocity covariance and topography are not in general found, although some association between ellipse orientation and topographic features in the Southern Ocean is discussed. Zonal versus meridional velocity variance is analyzed over regions of the Pacific Ocean in Huang et al. (2007) for observational data and a numerical ocean model. The magnitude of the ratio is dependent upon the time window over which the velocity field is averaged, prior to computing velocity variances via further averaging. A recent study in Stewart et al. (2015) extends these earlier studies using observational data and an ocean model with a resolution of one-twelfth degree. An association between increased velocity covariance anisotropy and topography is observed and, in the numerical model, an association between topographical slope and near-bottom covariance orientation is found. Differences between near-surface and near-bottom eddy velocity covariance properties are studied via decomposition of the model velocities into barotropic and baroclinic components.
It is known from linear theory that perturbations lean with the shear in stable configurations, and against the shear in unstable configurations [e.g., Pedlosky (1987), section 7.3]. Moreover in the idealized barotropic jet configuration of Waterman and Hoskins (2013) it is observed that eddy variance ellipses orient themselves in a manner consistent with stability properties of the mean flow. The time mean eddy variance ellipses show orientations consistent with instability in the upstream unstable region of the jet, and with wave radiation in the downstream region. Similar behavior is observed in Klocker et al. (2016, unpublished manuscript) in an idealized primitive equation model of the Antarctic Circumpolar Current.
The primary purpose of this paper is to provide an explicit link between predictions from linear stability theory and resulting eddy variance ellipses. The paper is organized as follows. Section 2 reviews the local geometric properties of the eddy velocity covariance tensor: the kinetic energy, horizontal eddy momentum stress anisotropy, and horizontal eddy momentum stress tilt. These are related to eddy variance ellipse properties, and their relation to wave group propagation is discussed. In section 3 these geometric properties are calculated for a piecewise linear shear layer and a piecewise linear jet on an f plane. Analytic solutions are derived from linear theory, and in the latter case these solutions are compared against calculations from a fully nonlinear numerical model. In a further numerical calculation the geometric properties are diagnosed for a piecewise linear jet on a β plane, and the results are compared against the results for a zonally evolving barotropic jet described in Waterman and Hoskins (2013). In section 4 ray-tracing theory is used to calculate eddy variance ellipse properties for the piecewise linear jet on a β plane. The predicted spatial patterns of the geometric properties are compared against those diagnosed from the numerical model. The paper concludes in section 5.
2. Eddy momentum stress decomposition
a. Geometric decomposition and interpretation
The mean barotropic vorticity equation is
where forcing and dissipation have been neglected, ∇ is the horizontal gradient operator, and t is time. A bar signifies a mean quantity and a prime denotes a deviation from the mean.3 The absolute vorticity is
where f is the planetary vorticity and u is the nondivergent velocity with x and y components u and υ, respectively.
where the eddy momentum stress tensor has the following components:
Here K is the eddy kinetic energy:
and M and N are the eddy Reynolds stresses:
where three geometric parameters appear: the eddy kinetic energy K, an eddy momentum stress anisotropy , which is bounded (0 ≤ γ ≤ 1), and an eddy momentum stress tilt θ = 1/2atan(−N/M) [see Marshall et al. (2012) for further details].
At a given spatial location, the eddy velocity forms some general trace in a u′–υ′ plot. The geometric properties in Eq. (7) can be related to an elliptical fit to this trace, with components:
Here ϕ parameterizes the ellipse and , , and ε determine the specific ellipse geometry. If the bar corresponds to a long time average, then the eddy stress and eddy kinetic energy resulting from the observed local eddy velocity u′ are identical to those resulting from an eddy velocity that, on a u′–υ′ plot, traces along this ellipse at a constant rate of change of ϕ. This ellipse is illustrated in Fig. 2a, and its relation to the signs of the Reynolds stresses M and N is shown in Fig. 2b.
With this ellipse definition, the eddy kinetic energy K is equal to one-quarter of the mean square of the semimajor and semiminor axes, and the tilt θ is the angle of the ellipse major axis to the x axis. The ratio of the semiminor axis B to semimajor axis A, r = B/A, is related to the anisotropy parameter γ via (Hoskins et al. 1983; Marshall et al. 2012):
The ratio of the focus separation 2F to the major axis 2A, equal to the ellipse eccentricity e, is given by
Additional properties of the eddy velocity covariance tensor have previously been discussed in the literature. In Huang et al. (2007) and Scott et al. (2008) a nondimensional zonal-meridional anisotropy is defined4:
which is related to the geometric decomposition via
This is further related to the geometry of the ellipse via
where F is one-half of the ellipse focus separation.
b. Spatial structure and mean flow forcing
The mean momentum equation is
where again forcing and dissipation have been neglected, p is the pressure, and ρ0 is a constant reference density. Consider a flow oriented in the x direction (or, more generally, rotate coordinates locally so that the flow is in the x direction). Then substitution of the geometric decomposition [Eq. (7)] results in a zonal eddy momentum forcing:
which relates the local mean momentum forcing by the eddies to the geometric parameters K, γ, and θ, and to their local gradients.
Now limiting consideration to the case where zonal averaging is applied, the zonal mean zonal velocity forcing by the eddies is
For zonally periodic systems, such as the idealized systems considered in this article, this local mean momentum forcing by the eddies cannot be balanced by a mean pressure gradient, and corresponds directly to a local mean acceleration. Differing geometric parameter patterns that lead to an eastward momentum tendency are illustrated in Fig. 3, and patterns that lead to a westward momentum tendency are illustrated in Fig. 4.
In more general cases, only the nondivergent component of the eddy momentum tendency is associated with a local forcing (see e.g., Maddison et al. 2015). A more general approach is taken in Waterman and Lilly (2015), where instead spatial patterns of geometric parameters that lead to local forcing of the mean absolute vorticity are considered.
c. Wave group propagation
The intrinsic group velocity of a slowly varying wave packet under the Wentzel–Kramers–Brillouin (WKB) approximation (Buhler 2009, section 2.1.6) is related to the ellipse tilt angle θ. This provides a link between the geometric properties of the eddy velocity covariance and wave activity propagation, which will be used to develop an analysis based upon ray tracing in section 4. For further details see Hoskins et al. (1983) (and also Maddison and Marshall 2013; Waterman and Hoskins 2013).
Consider a plane wave with streamfunction perturbation:
where k is the zonal wavenumber, l is the meridional wavenumber, and ω is the angular frequency, with . It follows that this leads to an anisotropy γ = 1 and a tilt:
The intrinsic group velocity (group velocity minus mean velocity) is [Eqs. (11)–(15) in Hoskins et al. (1983); Eq. (2.60) in Maddison and Marshall (2013); and Eq. (4) in Waterman and Hoskins (2013)] as follows:
where the mean potential vorticity gradient is treated here as a column vector. Substitution of the geometric decomposition in Eq. (7) leads to [see e.g., Eq. (16) in Hoskins et al. (1983)] the following:
(noting that γ = 1 here) with
where and are the angles that the mean potential vorticity gradient and intrinsic group velocity make with the x axis, respectively.
3. Geometric decomposition for piecewise linear flows
The approach used here for framing the wave instability problem is the counterpropagating Rossby wave (CRW) perspective (following Rossby 1939; Bretherton 1966a), which is an insightful pedagogical framework for explaining the fundamentals of linearized shear flows instabilities (see also Hoskins et al. 1985; Davies and Bishop 1994). Heifetz et al. (1999) performed a case study on the CRW interaction in the simple barotropic Rayleigh model (Rayleigh 1880), as well as a modified version consisting of a jetlike flow featuring two strips of vorticity with opposite signs. The growing normal-mode solutions for these configurations, obtained by Heifetz et al. (1999), are used in the current study to diagnose and interpret the eddy fields in the geometric decomposition framework.
a. Piecewise linear shear layer on an f plane
Let the bar and prime now indicate zonal mean and perturbation, respectively. A linear perturbation evolves via the linearized absolute vorticity equation:
The PV gradient is concentrated in two δ functions at y = ±b. Vorticity perturbations away from these locations can neither grow nor decay, and hence the vorticity perturbations with time-dependent amplitude are concentrated in δ functions at y = ±b. For a given wavenumber k the corresponding total perturbation is written as
where are the phases of two waves at y = ±b. In an unstable configuration the ratio of the two wave amplitudes asymptotes to one as the waves grow (Heifetz et al. 1999), and hence here their amplitudes are both set equal to Qk(t). The perturbation streamfunction is (following Heifetz et al. 1999) as follows:
Noting that, subject to zonal averaging:
it follows that in regions where q′ = 0 the Reynolds stresses M and N are independent of y, and hence M and N take uniform values in each of the three regions y < −b, −b < y < b, and y > b. Specifically, outside of the shear layer, |y| > b, the velocity perturbations are
as required. However inside of the shear layer, |y| < b, the velocity perturbations are
for −b < y < b, where is the (time independent) phase difference between the two waves. The linearized normal-mode solution then yields an analytic expression for the geometric parameters in the eddy stress tensor decomposition within the shear layer:
where the sign in Eq. (32c) is chosen so that −π/2 ≤ θ ≤ π/2.
The Reynolds stresses vanish outside the shear layer and have constant values within, in the latter case depending on the phase difference between the two interacting waves. As illustrated in Fig. 5, |M| (|N|) is maximized when the phase difference is zero (π). Outside the shear layer, the superposition of the velocities always tends to cancel out M and N, regardless of the phase difference between the waves. Inside the shear layer, the phase difference determines how the velocities interfere, and hence the resulting values of M and N. However, these values (and hence the tilt and anisotropy) remain constant within the shear layer. This result is generalized in appendix A for a continuous vorticity profile and shown to hold between any two vorticity waves.
The eddy kinetic energy K, anisotropy γ, and eddy ellipses, corresponding to the most unstable mode, are shown in Fig. 6. The eddy ellipse tilt is constant within the shear layer. The eddy kinetic energy and anisotropy both vary inside the shear layer, with a minimum in the kinetic energy and maximum in the anisotropy at the center. However, their product, equal to the dimensional anisotropy L = γK, is constant. This corresponds, as per Eq. (14), to a constant eddy ellipse focus separation within the shear layer. Outside the shear layer the eddy “ellipses” are circular, with undefined tilt and with decreasing eddy kinetic energy as y → ±∞.
Here N is constant within the shear layer and vanishes outside. This indicates that there is no mean momentum tendency due to the eddies away from y = ±b. At y = +b, moving across the boundary from inside to outside the shear layer, the eddy kinetic energy and anisotropy both decrease (discontinuously) with an eddy tilt −π/2 < θ < 0 on the nonzero shear side. This corresponds to cases A and C of Fig. 3, and hence to a westward momentum tendency. Conversely at y = −b, moving across the boundary from outside to inside the shear layer, the eddy kinetic energy and anisotropy both increase (discontinuously) with an eddy tilt −π/2 < θ < 0 on the nonzero shear side. This corresponds to cases A and C of Fig. 4, and hence to an eastward momentum tendency. In both cases the unstable linear perturbation implies a local mean momentum tendency that opposes and decelerates the flow. In this idealized case the implied deceleration is applied over δ functions at y = ±b.
As discussed in Heifetz et al. (1999), unstable solutions are composed of two phase-locked counterpropagating Rossby waves at y = ±b, with the waves mutually amplifying each other and with a phase difference 0 < εk < π. Conversely, for stable solutions the waves mutually destroy each other, with a phase difference of −π < εk < 0. From Eq. (32c), the eddy ellipse for an unstable perturbation leans against the shear, with −π/2 < θ < 0, while the eddy ellipse for a stable perturbation leans with the shear, with 0 < θ < π/2. This behavior is as discussed in Marshall et al. (2012).
Linear stability analysis gives that the most unstable mode has wavenumber kmax with 2kmaxb ≈ 0.797 and phase difference [values from Heifetz et al. (1999)]. This yields an eddy tilt within the shear layer of θ ≈ −0.177π. Note that this differs from the value −π/4, which might be expected based upon the argument that the eddies should lean maximally according to the local shear. This maximal tilt of −π/4 would corresponds to a phase difference of ε = π/2. Linear perturbation growth requires phase locking between the perturbations at y = ±b, and hence the phase difference εk and wavenumber k are not independent for growing perturbations—the choice of one implies a value for the other. As discussed in Heifetz et al. (1999), the phase difference of the fastest growing mode arises as a compromise—in particular, increasing εk beyond (leading to an increased eddy tilt) implies an increased value of k, a decrease in the amplitude of each wave as seen by the other, and a net decrease in the perturbation growth rate. This behavior is a manifestation of the nonnormality of the system (Heifetz and Methven 2005).
b. Piecewise linear jet on an f plane
Following a similar procedure as for the piecewise linear shear layer, one can find the analytic linear normal-mode solutions for a piecewise linear jet, which consists of two constant and opposite sign vorticity strips (e.g., Heifetz et al. 1999).
The background profile (Fig. 7a) is given by
In this case a nonneutral linear perturbation consists of three waves at each of the interfaces y = 0, ±b, and the vorticity perturbation for a given wavenumber k is written as
For a symmetric normal-mode solution, which corresponds to equal growth rate of all three waves, one finds , , and the corresponding streamfunction perturbation is (Heifetz et al. 1999)
where is a time-independent phase difference.
On each side of the jet, the Reynolds stresses are again independent of y, and the geometric properties of the eddy ellipse are given by
where the negative sign in Eq. (36c) corresponds to the northward side of the jet with negative shear (0 < y < b) and the positive sign corresponds to the southward side of the jet with positive shear (−b < y < 0). The eddy ellipse tilt is no longer directly proportional to the normal-mode phase difference due to the interaction with the third vorticity wave, which gives an additional 0.5e−kb term in the denominator of Eq. (36c).
To examine the validity of the analytic normal-mode solution, we compare the calculated eddy ellipse tilt and anisotropy to those found from a fully nonlinear simulation (Fig. 7). Details regarding the numerical model are given in appendix B. The nondifferentiable piecewise-linear jet in Eq. (33) is approximated by a smooth tanh profile whose vorticity is given by
where Λ is the shear, the interfaces of the layers are located at y = 0, ±b (b = 75 km), and d is a parameter that measures the relative thickness of the transition regions at the interfaces. In the limit d → 0, one recovers the piecewise linear profile. Here a small value of d = 0.05b is chosen. The zonal mean flow of the approximated piecewise linear jet at the initial moment, at t = 0 days, is shown in Fig. 8a (black solid line).
A snapshot of the relative vorticity at early stages of the simulation, after t = 57 days (Fig. 7b), shows that the vorticity perturbation is localized at the three interfaces y = 0, ±b. The system picks up a mode n = 3 solution, corresponding to a wavelength of λ ≈ 290 km (the x domain size is Lx = 872 km) and a normalized wavenumber of 2kb = 3.24. Because of the quantization requirement of our periodic domain, this is slightly larger than the marginally most unstable normal mode 2kmaxb = 2.452 [given in Heifetz et al. (1999), for the case of a jet in an infinite domain]. Our normal-mode solution has a phase difference of εk ≅ 0.418π [substituting 2kb = 3.24 in Eq. (20) of Heifetz et al. (1999)], which gives from Eq. (36c) an eddy ellipse tilt of θ ≅ 0.194π on the northward side of the jet and θ ≅ −0.194π on the southward side of the jet. Hence, the eddy ellipse is tilted against the shear in each side of the jet, consistent with an unstable normal-mode solution, fluxing momentum out of the jet core and into its flanks. Note that the tilt angle found here is not far from the one found for the single shear layer in the previous section (θ ≅ −0.177π). In addition, it is similar to the value one finds for the most unstable normal-mode solution of the jet in an infinite domain case, θ ≅ −0.208π, achieved by substituting in Eq. (36c) the theoretical values for the infinite domain case (2kmaxb = 2.452 and εk ≅ 0.628π).
Figures 7c,d show the zonally averaged eddy ellipse tilt and anisotropy, respectively, from the numerical simulation (dashed red line) and the analytic normal-mode solution (blue line). The analytic solution is calculated by interpolating the streamfunction expression in Eq. (35) for the normal-mode solution onto a numerical grid and differentiating to find the velocity field. The results are in good agreement, and show that the eddy tilt is constant in each side of the jet, positive (negative) in the negative (positive) shear, with a jump toward ±π/2 at y = 0. The anisotropy structure in each of the layers is similar to that for the single shear layer case, though here it maximizes closer to the edges y = ±b rather than in the middle of the layer. Both M and N are constant within the layers, and hence only the perturbation kinetic energy K controls the meridional structure of the anisotropy γ (which is inversely proportional to K).
The constant tilts remain close to the theoretical value even at later stages of the simulation, when the waves merge into vortices (not shown). At much later stages of the development, however, the linear solution does not describe the dynamics adequately. In addition, once the mean vorticity gradient becomes nonzero, the tilt is no longer constant. This is examined in more detail in the next section in which a background planetary vorticity gradient is added. Since there are many similarities between the results obtained on the f plane and β plane, a more detailed discussion is deferred until the latter case.
c. Piecewise linear jet on a β plane
A more physically relevant case is when a background planetary vorticity gradient is included. The normal-mode solution derived in the previous section is no longer valid in this case. However, the fastest-growing mode and its growth rate are very similar to those found for the zero beta case. To describe the meridional dependence of the tilt, ray tracing will be applied in section 4, but first a qualitative description of the evolution of the flow is given.
1) Time evolution and geometric decomposition
Figure 8 shows the initial flow as well as selected times of the zonal mean velocity (Fig. 8a), mean potential vorticity (Fig. 8b), and mean potential vorticity gradient (Fig. 8c). The initial zonal mean velocity (Fig. 8a, black solid line shows t = 0 days) is identical to that used in the zero beta case (i.e., the smooth approximation of the piecewise linear jet). It satisfies the necessary conditions for barotropic instability, since the absolute vorticity gradient changes sign within the domain (Fig. 8c, black solid line). By t = 135 days ≈ 7τ [where τ = 1/(Λe−kb sinεk) ≈ 1/(6 × 10−7 s) ≈ 19 days, the inverse growth rate of the fastest-growing wavenumber from the simulation, calculated from the linear theory for the unbounded jet on the f plane given in Heifetz et al. (1999)], the mean flow has become significantly weaker (red solid line in Fig. 8a), due to the transfer of energy from the mean flow to the eddies. At this instant, the behavior of the system changes: ceases to change sign within the domain (red solid line in Fig. 8c), so the jet can no longer support the existence of linearly unstable solutions. Later on, jet sharpening can be identified (dashed blue line in Fig. 8a, t = 172 days ≈ 9τ days).
Snapshots of the total potential vorticity are shown in Fig. 9. At early times, the system behaves linearly, and the potential vorticity at the interfaces is characterized by a wavy structure, with a mode n = 3 corresponding to a wavenumber k = 3.7 × 10−6 m−1 or normalized wavenumber of 2kb = 3.24 (Fig. 9a, t = 57 days ≈ 3τ). At later times (Fig. 9b; t = 115 days ≈ 6τ), the solution is dominated by vortices, which appear to be oriented (in spatial structure) against the shear in each layer. The orientation of the vortices changes shortly afterward (Fig. 9c, t = 155 days ≈ 8τ), and become tilted with the shear (in spatial structure). This is consistent with the jet sharpening that was identified (dashed blue line in Fig. 8a), where the mean flow strengthens as a result of upgradient eddy momentum fluxes into the core of the jet. Eventually, (Fig. 9d; t = 270 days ≈ 14τ), the vortices break nonlinearly and the system equilibrates.
The transition in the behavior of the eddy–mean flow interaction can be seen by calculating the eddy variance ellipse at different locations and times of the simulation. Figure 10 shows scatterplots of u′ and υ′ at the interfaces y = ±b, 0 for three different times, as well as their corresponding ellipses. At t = 57 days ≈ 3τ (left column), the anisotropy is close to one near the interfaces, and the eddy ellipses are elongated. In addition, the tilts are consistent with the linear instability picture [i.e., positive (negative) tilt on the negative (positive) shear on each side of the jet]. By t = 115 days ≈ 6τ, anisotropy at the flanks y = ±b has become gradually smaller (hence the eddy ellipses appear more round), though they are still characterized by a tilt that implies barotropic instability. However, this changes at later times (t = 155 days ≈ 8τ), when the ellipses at y = ±b are tilted oppositely, implying that the eddies are fluxing momentum upgradient into the jet, with backscatter of eddy energy to the mean flow. Note that at y = 0 (the center of the jet), we always find a large γ, and θ = π/2 (in fact, θ = π/2 or 0 is required by symmetry). Equivalently, N ≅ 0 and M > 0, the latter implying that so that the eddy variance ellipse is elongated in the meridional direction.
Figure 11 shows the temporal development and meridional dependence of the zonally averaged eddy ellipse tilt and anisotropy, and of the Reynolds stresses M and N. Figure 11a confirms that during the initial development of the system, in the time interval t < 135 days ≈ 7τ, the eddy tilt shows the signature of instability with positive tilts on the northward side of the jet and negative tilts on the southward side of the jet. Gradually, the tilt (absolute value) on both sides of the jet approaches π/2, until at t = 135 days ≈ 7τ all the eddies within the jet are characterized by θ = π/2. The tilt then flips, consistent with the onset of stability and the momentum fluxes being directed upgradient into the core of the jet. This holds until approximately t = 172 days ≈ 9τ. After that, there is an indication for another flip in the tilt, implying that the eddies are strengthening again. However, after t = 210 days ≈ 11τ the tilt is single signed within the jet, and its sign changes periodically. This is consistent with the sign of N oscillating periodically after t = 210 days ≈ 11τ (Fig. 11c).
Figure 11b shows the time evolution of the eddy anisotropy. During the initial time development there are three distinct regions where anisotropy is close to unity, at each of the interfaces y = ±b, 0. During the transition period (t = 135 days ≈ 7τ), this picture changes drastically, with anisotropy being maximized in a broad meridional region around the jet center (the quasicircular area of anisotropy close to unity around t = 135 days ≈ 7τ). After t = 187 days ≈ 10τ, eddy anisotropy becomes significantly smaller. However, there is a secondary (though smaller) peak in anisotropy at the jet core, around t = 238 days ≈ 12.5τ. This is consistent with a similar secondary maximum observed in M (Fig. 11d).
Looking immediately to the north of the jet in the early stage of the evolution (at around t = 57 days ≈ 3τ), near the jet core the tilt decreases with latitude with π/4 < θ ≤ π/2, corresponding to panel F of Fig. 4 and a deceleration of the mean flow. The anisotropy decreases with latitude while the eddy kinetic energy increases, with 0 < θ < π/2. Their combined effect is a decrease in the value of L with latitude and a net acceleration of the mean flow. Hence, there is a competing effect between the anisotropy, which acts to accelerate the jet core, and the tilt and eddy kinetic energy, which act to decelerate the jet core. The net combined effect is a deceleration of the zonal mean flow near the jet core. Conversely, after the tilt angle flips, the tilt increases with y, while kinetic energy and anisotropy decrease with y. Hence, all three act to accelerate the flow (Figs. 3f, 3a, and 3c, respectively).
2) Intrinsic group velocity and relation to eddy ellipse tilt
During the unstable period (t = 135 days ≈ 7τ), the eddy ellipse tilt is consistent with an inward directed group velocity. Shortly after the initial time development, by t = 25 days ≈ 1.3τ), the implied group velocity direction is almost independent of y and (the meridional component of y) is negative on the northward side of the jet, and positive on the southward side of the jet. Hence, momentum is fluxed out of the jet. The direction of the implied intrinsic group velocity at t = 25 days ≈ 1.3τ is approximately on the northward and southward sides of the jet, respectively. Shortly afterward, however, the direction of the implied intrinsic group velocity is no longer independent of y; instead it increases in the northward side and decreases in the southward side of the jet, with at the jet core.
As the jet transitions to the stable regime, once the absolute vorticity gradient ceases to change sign (at t = 135 days ≈ 7τ) and the eddy ellipse is meridionally elongated (θ = π/2), all the rays become zonal (). Immediately thereafter, we enter the stable regime, with outward-radiating rays. The sign of the eddy tilt flips and correspondingly the meridional component of the implied group velocity changes to positive on the northward side of the jet and negative on the southward side of the jet. Going away from the jet core, the ray becomes refracted such that it becomes more meridional.
3) Comparison with Waterman and Hoskins (2013)
A comparison can be made between the results obtained here and those obtained by Waterman and Hoskins (2013) for the case of a zonally evolving jet in statistically steady state. In Waterman and Hoskins (2013) an unstable jet is forced at the boundary of the domain, and the geometric properties of the eddy velocity covariance tensor are defined via time mean averaging, rather than zonal averaging. The analysis shows a similar picture to that described here: eddies tilting against the shear in the upstream unstable region, and an eddy tilt consistent with wave radiation in the downstream stable region. The broad structure of M, N, K, and the corresponding eddy ellipses described in Waterman and Hoskins (2013) are remarkably similar to those found here, though, importantly, in Waterman and Hoskins (2013) the structure varies spatially along the flow direction, whereas here the structure varies temporally through the flow evolution.
There are, however, some key differences. First, in Waterman and Hoskins (2013), the location where the eddy kinetic energy maximizes occurs downstream of the location where the potential vorticity gradient ceases to change sign as a consequence of mean flow advection, whereas here they occur simultaneously in time (at t = 135 days ≈ 7τ). At the point of maximum eddy kinetic energy Waterman and Hoskins find a “bullet of the M,” implying meridional elongation of the eddies there. Downstream of the eddy kinetic energy maximum, Waterman and Hoskins find that forcing solely from M is responsible for strengthening and extending the jet, eventually forcing their time-mean recirculation gyres. A similar bullet of M is found here (Fig. 11d) and meridional elongation at the time of maximum eddy kinetic energy. However, since a zonally symmetric jet is being considered here, the zonal mean vorticity forcing from M vanishes identically, and the mean flow forcing arises solely from N. Hence, no recirculation gyres can develop in our configuration, and the eddy forcing is instead responsible only for the acceleration or deceleration of the jet.
4. Ray tracing
For a barotropic jet on a β plane during the unstable regime, the eddy ellipse tilt is no longer constant within the shear layers (see Fig. 14a), but rather it increases toward the jet core. This effect is due to the nonzero potential vorticity gradient within the shear layers, and is intrinsically related to the wave propagation there.
Here ray-tracing theory is employed to study the propagation of waves within the shear layers under the influence of the β effect and a constant shear. The analytic solution obtained from the ray equations agrees well with the fully nonlinear simulation, aiding in the prediction of the meridional structure of the eddy ellipse tilt.
a. Theoretical background
The wave activity propagates at the group velocity. In a homogeneous medium a ray (which is the path parallel to the group velocity at every point) will propagate in a straight line. In an inhomogeneous medium, however, refraction can occur. Ray-tracing theory [Whitham (1974); Lighthill (1977), also see Buhler (2009) and Salmon (1998) for overviews] gives the leading-order asymptotic description of a slowly varying wave packet in a medium that varies slowly compared to the scale of the waves [through the WKB approximation, see Hoskins and Karoly (1981); Hoskins and Ambrizzi (1993)].
For such conditions, the streamfunction can be represented locally by a plane wave
are slowly varying. Note that here a distinction is made between the local value of the angular frequency, ω(x, y, t), and the dispersion relation, Ω(k, l, x, y, t), with ω(x, y, t) = Ω[k(x, y, t), l(x, y, t), x, y, t].
Cross differentiation yields the ray equations
along rays defined by
These are equivalent to Hamilton’s equations. In the absence of explicit time dependence in Ω(k, l, x, y, t), the analog of the Hamiltonian, is conserved along the ray path.
b. Analytic ray-tracing solution
In the case considered here, away from the interfaces, the solutions can be approximated as plane wave whose dispersion relation is governed by the Rossby wave dynamics and given by
It follows that
where ± refers to y > 0 and y < 0, respectively. Hence,
where k0 and l0 are initial zonal and meridional wavenumbers, respectively.
where is positive in the northward side of the jet (y > 0) and negative in the southward side of the jet (y < 0). This expression allows us to investigate how the ray solution depends on the parameters of the problem, namely, β, k0 and the shear Λ. Note that the solution, conveniently, does not depend explicitly on l0, but rather on the initial tilt θ(y0) = θ0, which can be calculated directly from the Reynolds stresses.
For example, for an outward-radiating ray on the northward side of the jet, originating at the jet center, we would have y0 = 0, θ0 = π/2, and , giving
which is only valid for .
For nonzero β and Λ, going away from the jet core toward the flank at y = b, the tilt becomes gradually smaller, approaching θ → 0 for . This is consistent with the results of the previous section for the outward-radiating ray (cf. Fig. 12). The analytic solution also implies that for larger wavenumber, larger shear, or smaller β, the tilt decreases faster when moving away from the jet core, which means that the ray is refracted more. Since it is the shear that is responsible for the refraction of the rays (through changing the meridional wavenumber t), it is clear that a larger shear will cause larger refraction. The wavenumber dependence is also clear, since larger wavenumbers corresponds to smaller waves, which are more easily influenced by the background shear. Finally, β has a competing effect with the shear.
c. Comparison with numerical results
The validity of the analytic solution is now verified by comparing it with numerical ray-tracing results. The ray in Eqs. (44) and (45) are solved using a basic forward Euler time stepping. In Fig. 13 results of the ray tracing are plotted for different choices of initial zonal wavenumber, k0. For the inward-radiating ray, the ray is initiated just slightly below and above the upper and lower interfaces y = b, −b, respectively. The variable k0 is specified and l0 computed such that the eddy ellipse tilt (assuming plane wave solutions) is equal to the value diagnosed at that point from the simulation. That is, given k0 and y0, we find l0 such that
which gives two possible solutions for l0. The sign is chosen such that the ray is directed toward the jet (i.e., l0 < 0 in the northward side of the jet and l0 > 0 in the southward side of the jet). The possible values for k0 can be estimated using . This gives a typical value of ≈ 2.06 × 10−5 m−1 or, in normalized units, (where b = 75 km is the width of the shear layer in each side of the jet). This is similar to the wavenumber that emerged in the simulation (mode n = 3, which corresponds to 2kb ≅ 3.24).
In Fig. 13a the implied intrinsic group velocity from the ray tracing is plotted, for rays with varying zonal wavenumbers , , , , and , where the leftmost ray corresponds to the smallest wavenumber. The track for is highlighted in blue. In all cases, the initially southwestward-pointing ray ends up pointing eastward (i.e., tends toward zero as the rays approach the core of the jet), just as was found in the results from the numerical simulation, analyzed in section 3b. The meridional component of the intrinsic group velocity, however, remains negative in the northward side of the jet, implying a positive eddy ellipse tilt, and vice versa in the southward side of the jet (i.e., the momentum fluxes point out of the jet). Since the rays are plotted as a function of the zonal location x rather than time (as in Fig. 12), in Fig. 13b the time development of the highlighted ray with is plotted, in the same zonal location (the uppermost arrow corresponds to t = 0, and the ray deflects as time progresses and it propagates toward the jet core). Note the remarkable similarity with Fig. 12 for the inward-radiating ray during the unstable regime (albeit with a specific choice of k0).
For the outward-radiating ray, the ray is initiated very close to the jet core (y = 0) and hence l0 = 0 is chosen, which implies θ = π/2 and . The sign of the ray-tracing solution is chosen such that the ray is directed outward from the jet (i.e., l is positive on the northward side of the jet and negative on the southward side of the jet). In Fig. 13c, the intrinsic group velocity from the ray tracing is plotted for the outward-propagating rays for zonal wavenumbers , , , and , where the leftmost ray corresponds to the largest wavenumber. This is consistent with our earlier inspection of the analytic solution for the outward-radiating ray, where it was found that larger wavenumbers will tend to be refracted more. The track for is highlighted in blue. In all cases, increases on the northward side and decreases in the southward side toward the jets flanks, similar to that observed in the simulation for the stable outward-radiating regime. The meridional component of the intrinsic group velocity remains positive, which implies a negative eddy ellipse tilt in the northward side of the jet and vice versa on the southward side of the jet (i.e., momentum fluxes into the jet). In Fig. 13d, the time development of the highlighted ray with is plotted, in the same zonal location. This is remarkably similar to Fig. 12 for the stable outward-radiating regime (albeit, again, with a specific choice for k0).
Finally, in Figs. 14a,b the eddy ellipse tilts from the simulation and the ray tracing are compared, for the inward- (t = 57 days ≈ 3τ) and outward-radiating (t = 155 days ≈ 8τ) rays, respectively. This uses both the actual value for the tilt (black thick line), calculated from the simulation using the Reynolds stresses, as well as the estimated value (dashed thin red line) for a plane wave solution in Eq. (19) using the estimated meridional and zonal wavenumbers , , respectively, evaluated locally at every y. For the ray-tracing solution, is chosen for the inward-radiating ray and is chosen for the outward-radiating ray (the blue rays in Figs. 13a and 13c, respectively). As a check for the analytic solution from the ray-tracing equations, both the numerical (blue thin line) and the analytic solution (blue stars) are plotted for the eddy ellipse tilt, and these indeed coincide. The ray tracing agrees well with the simulation, albeit with the particular choices for the zonal wavenumbers as described above, and captures correctly the overall dependence of the tilt with y. These results validate the use of ray-tracing predictions to predict eddy tilt.
5. Summary and discussion
In this manuscript, the role of eddy Reynolds stresses in accelerating and decelerating barotropic ocean jets has been revisited. In particular, the eddy Reynolds stresses have been analyzed by exploring a geometric decomposition of the eddy stress tensor. This decomposition involves describing the eddy stress tensor in terms of an eddy variance ellipse, the geometry of which characterizes the mean eddy shape and orientation, the direction of eddy activity propagation, and the eddy forcing of the mean flow.
Idealized linear shear and jet profiles have been analyzed and analytical results have been compared against fully nonlinear simulations. For flows with zero planetary vorticity gradient, analytic solutions have been obtained that provide a direct relationship between the geometric eddy tilt and the phase difference of a normal-mode solution. This allows a straightforward interpretation in terms of familiar concepts in classical stability theory. The initially unstable jet gives rise to eddies which are tilted “against the shear”; hence, perturbations can grow and extract energy from the mean flow. However, as the jet becomes gradually weaker, at a given moment it can no longer satisfy the conditions for instability. Once the jet stabilizes, eddies become tilted “with the shear,” all perturbations decay as they return their energy to the mean flow and strengthen the jet. For cases with a nonzero planetary vorticity gradient, ray-tracing theory has been used to investigate eddy propagation within the jet, and make predictions for the evolution of the eddy wavenumber, which in turn can indicate the evolution of the eddy ellipse tilt. An analytic solution for the eddy tilt is found for the case of a linear plane Rossby wave propagating on a constant background shear. The ray-tracing solution captures the essence of the observed eddy propagation and agrees well with the eddy tilt diagnosed from a fully nonlinear simulation, subject to a choice of initial zonal wavenumber.
We propose that the geometric framework explored in this manuscript could be used as a diagnostic tool to understand the role of Reynolds stresses in maintaining and decelerating inertial jets in ocean models and observations. For example, similar ideas have already been applied to separated western boundary currents such as the Gulf Stream and Kuroshio (Waterman and Hoskins 2013) and are currently being applied to the Southern Ocean to elucidate zonal jets embedded within the Antarctic Circumpolar Current (Klocker et al. 2016, unpublished manuscript). For such problems, it would be interesting to apply the same ray-tracing methods employed here to see if they can provide similar insights for flows that vary in a greater number of dimensions.
Moreover, we propose that the approach taken here might be employed to develop a simple parameterization of eddy Reynolds stresses for ocean general circulation models that are able to (at least partially) resolve inertial jets. While it is surely impractical to contemplate solving ray equations in such models, one could imagine assuming that some of the eddy energy generated through baroclinic instability is back-scattered to the mean flow by upgradient momentum transfer, as found in some of the idealized jet profiles in this manuscript, and proposed, for example, by Marshall and Adcroft (2010) and Jansen and Held (2014). The linear theory presented here suggests that the eddy tilt slightly deviates from that which might be expected from the argument that eddies should lean maximally with the mean shear, although the general picture of unstable eddies leaning (in the sense of eddy variance ellipse tilt) against the mean shear is observed. The anisotropy and eddy kinetic energy have in some cases been observed to counteract each other’s effects. In the zero beta cases the anisotropy and eddy kinetic energy both vary within the nonzero shear regions, but their product is constant. For the jet on a β plane in the immediate jet core the eddy kinetic energy decreases with latitude, while the anisotropy increases.
The findings of this paper might be useful ingredients for a parameterization of horizontal eddy momentum fluxes. Exploiting the fact that the component of the eddy stress tensor involving the eddy Reynolds stresses is bounded by the eddy kinetic energy, assuming the momentum fluxes are directed either upgradient or downgradient, depending on the mean flow stability properties, and prescribing a typical value for the eddy anisotropy, it should be possible to develop a parameterization that is both energetically consistent and rooted in the underlying geometry of the eddy dynamics. However, even at the simplest level there remain many questions to be addressed. For example, should such a parameterization be applied to the depth-integrated flow? Furthermore it will be necessary to model the formation, propagation and dissipation of eddy kinetic energy, as discussed by Eden and Greatbatch (2008), Marshall and Adcroft (2010), and Jansen et al. (2015). Despite these challenges, we believe this approach holds some promise and is worth exploring further.
The authors gratefully acknowledge useful discussions with Laure Zanna and helpful comments by Julian Mak. The authors also thank Stephanie Waterman and Malte Jansen for comments that significantly improved the article. TT acknowledges support by the Israeli Science Foundation through a grant to Yohai Kaspi (ISF 1310/12). JRM and DPM acknowledge the support of the Natural Environment Research Council Grant NE/L005166/1. The numerical code used in this article was developed from an original code provided by Pavel S. Berloff (see also Karabasov et al. 2009).
Generalization of the Relations for M and N between a Pair of Vorticity Waves
The two isolated δ-function vorticity waves can be generalized to the case where the vorticity field is continuous. Consider a zonal Fourier component of the vorticity anomaly,
inducing a streamfunction anomaly of the form
where the Green’s function is given by
After some algebra, we obtain
Hence, the values of M and N at some location y can be regarded as resulting from the continuum of infinite number of pairs of vorticity waves sandwiching y from below (y′ contributions in the integral) and from above (y″ contributions in the integral).
Numerical Model Description
For this study, we use the Parallel Quasi-Geostrophic Model (PEQUOD), a finite-difference code for solving quasigeostrophic equations in a rectangular domain, configured in a one-layer barotropic zonally periodic configuration. The numerical method implemented in PEQUOD and used here incorporates the Compact Accurately Boundary Adjusting high-Resolution Technique (CABARET) for advection of relative potential vorticity, combined with integration of the advection of planetary vorticity. Further details regarding CABARET can be found in Karabasov and Goloviznin (2009) and Karabasov et al. (2009). The potential vorticity inversion is performed using a fast Poisson solver using a customized version of FFTPACK.
The simulations are conducted in a rectangular domain −L ≤ x ≤ L, −0.5L ≤ y ≤ 0.5L where L =436 km, with a no-slip boundary condition applied to the perturbation from the background profile, at all lateral boundaries. The model is integrated using a grid of nx =256 and ny = 129 nodes in the zonal and meridional directions, respectively, corresponding to 3.4-km resolution and with a time step size of 52 s. The model is then run for 287 days. The eddies are defined as deviations from the zonal mean, and the corresponding eddy quantities such as eddy kinetic energy, Reynolds stresses, and the eddy variance ellipse parameters calculated accordingly.
The flow is initialized with the zonally symmetric jet U0, and perturbed with small random noise. The numerical model essentially solves
where q = ∇2ψ + βy is the absolute vorticity, ν = 106 m2 s−1 is a Laplacian viscosity coefficient parameter, and r = 5.2 × 10−7 s−1 is the relaxation time scale toward the background equilibrium flow qeq = q0. The jet strength is U0max = 0.24 m s−1, and the jet half-width is b = 75 km. This corresponds to a shear of Λ = U0max/b ≅ 3.2 × 10−6 s−1. For the case where the planetary vorticity gradient is nonzero, we use β = 2 × 10−11 m−1 s−1, which corresponds to a nondimensional β (which measures the relative importance of the planetary vorticity gradient relative to the mean shear) of . In midlatitudes, a typical ocean jet has U ~ 0.5 m s−1 at the surface (e.g., Sheen et al. 2013), whereas the barotropic jet used here has a value more typical of the depth mean over the main thermocline.
Analytic Solution for the Ray-Tracing Equations
where k0, l0 are the initial wavenumbers.
where yc is an integration constant. Similarly, using
These are the analytic solutions describing the ray path, following the intrinsic group velocity.
At t = 0, we find
Finally, solving for θ gives
which is Eq. (47), where is positive in the northward side of the jet (y > 0) and negative in the southward side of the jet (y < 0). The two possible ± solutions correspond to the two possible solutions for l.
In Wilkin and Morrow (1994) and Morrow et al. (1994), the semimajor axes of the ellipse scale with the principal velocity variances. In Waterman and Lilly (2015) the semimajor axes are set equal to the square root of the principal velocity variances. In this article, the semimajor axes are instead set equal to the square root of twice the principal velocity variances, in which case the ellipse yields a fit to observed values of eddy zonal velocity u′ and eddy meridional velocity υ′, on a u′–υ′ plot, as discussed in section 2 and shown in Fig. 10.
A crucial point to note is that that the eddy variance ellipse does not relate directly to the spatial structure of the eddy field. A locally highly isotropic eddy variance ellipse does not imply that the eddies are physically circular—rather it implies that, at that location, there is no component of the eddy velocity which preferentially leads to local eddy stress.
The averaging operator is a Reynolds operator, commutes with the partial derivatives ∂/∂x, ∂/∂y, ∂/∂t and satisfies a Cauchy–Schwartz inequality such that , with .