We investigate the relationship between Ertel potential vorticity Q and Bernoulli potential B on orthobaric density surfaces in the Antarctic Circumpolar Current (ACC), using the Southern Ocean State Estimate. Similar to the extratropical atmospheres of Earth and Mars, Q and B correlate in the ACC in a function-like manner with modest scatter. Below the near-surface, the underlying function relating Q and B appears to be nearly linear. Nondimensionalizing its slope yields “Ma,” a “Mach” number for long Rossby waves, the ratio of the local flow speed to the intrinsic long Rossby wave speed. We empirically estimate the latter using established and novel techniques that yield qualitatively consistent results. Previous work related “Ma” to the degree of homogeneity of Q and to Arnol’d’s shear stability criteria. Estimates of “Ma” for the whole ACC are notably positive, implying inhomogeneous Q, on all circumpolar buoyancy surfaces studied. Upper layers generally exhibit “Ma” slightly less than unity, suggesting that shear instability may operate within these layers. Deep layers exhibit “Ma” greater than unity, implying stability. On surfaces shallower than 1000 m just north of the ACC, the Q versus B slope varies strongly on subannual and interannual time scales, but “Ma” hovers near unity. We also study spatial variability: the ACC is speckled with hundreds of small-scale features with “Ma” near unity, whereas away from the ACC “Ma” is more commonly negative or above unity, both corresponding to stability. Maps of the time-mean “Ma” show stable regions occupy most of the Southern Ocean, except for several topographically controlled hotspots where “Ma” is always near unity.
Vorticity–streamfunction correlations provide a powerful diagnostic tool for studying steady and unsteady geophysical fluid systems, including the vertical propagation of planetary waves, starting with the classic work of Charney and Drazin (1961), and free-mode structures (Read et al. 1986; Marshall and So 1990). Free-mode structures emerge in the ocean limit of weak forcing and dissipation, and have been used to model inertial western boundary currents and their separated extensions (Fofonoff 1954; Charney 1955; Marshall and Nurser 1986; Marshall and Marshall 1992), the Antarctic Circumpolar Current (ACC) (Marshall et al. 1993; Marshall 1995), and interactions between the large-scale circulation and variable bottom topography (Bretherton and Haidvogel 1976; Marshall 1995). Such studies often employ the quasigeostrophic (QG) framework (Charney and Drazin 1961), being relatively tractable but requiring strong stratification (Vallis 2017). White (1990) brought attention to the primitive framework by showing Earth’s atmosphere exhibits function-like correlations with some scatter between the Ertel potential vorticity (PV) Q and the mass-weighted streamfunction Ψ. In such a primitive-variable framework, the Bernoulli potential B is related to Q and Ψ in the steady-state limit by Q = dB/dΨ (Gill 1982, 231–233). Consequently, one may equally well look for correlations in scatterplots of Q versus B (our focus), or Ψ versus B.
Vorticity–streamfunction studies have proven insightful for the Jovian planets. Using both QG and primitive shallow-water models applied to Voyager data, Dowling (1993) proposed that Jupiter is marginally stable with respect to Arnol’d’s (1966) second branch of shear stability in the QG framework. Observational evidence supports this claim, with some caveats (Read et al. 2006). Read et al. (2009a) reached a similar conclusion for Saturn. Taking this idea to its limit, Read et al. (2009b) estimated Saturn’s rotation rate—a challenging quantity to measure because Saturn’s magnetic field is aligned with its rotation axis—from meteorological data, by assuming that Saturn resides in a state of marginal stability.
Returning to the terrestrial planets, Du et al. (2015) used vorticity–streamfunction diagnostics to show that homogeneous PV, such as assumed by Sun and Lindzen (1994), is far from the norm throughout Earth’s extratropical atmosphere. On Mars, Dowling et al. (2017) determined that PV is nearly homogeneous only in the planet’s southern three-quarters during southern summer, but otherwise is inhomogeneous.
This paper explores similar ideas in the context of Earth’s Antarctic Circumpolar Current. We focus on the “Mach” number for Rossby waves, denoted “Ma.” This emerges as a fundamental parameter, nondimensionalizing the slope of the vorticity–streamfunction relationship. It is also intimately related to the stability properties of the background flow, as discussed theoretically in section 2. We translate these ideas into the 3D primitive setting in section 3, detailing a new empirical method that estimates, on buoyancy surfaces, the eddy phase speed and “Ma.” We use the Southern Ocean State Estimate (SOSE; Mazloff et al. 2010), in section 4 to study the relationship between Q and B, the eddy phase speed, and “Ma” in the ACC. Section 5 concludes. Our empirical method is mathematically elaborated upon in appendix A, and substantiated by empirical tests in appendixes B, C, and D.
2. Theoretical background
The stability of a shear flow is linked to horizontal gradients of PV (e.g., Charney and Stern 1962). In the absence of other constraints to stabilize the flow, eddies tend to flux PV down the mean PV gradient. This idea forms the theoretical basis for some eddy parameterizations (Green 1970; Treguier et al. 1997). The eddy PV flux is closely related to the eddy force on the mean flow via the divergence of the Eliassen–Palm flux tensor (Young 2012; Maddison and Marshall 2013) which, in an integral sense, can either accelerate or decelerate the mean flow depending on the sign of the relation between Q and Ψ (Marshall and Adcroft 2010; Marshall et al. 2012; see also Maddison et al. 2015). Moreover, large-scale gradients in PV along a buoyancy surface are related to the sign of eddy-induced advective transport along that buoyancy surface, as demonstrated both theoretically and numerically by Marshall et al. (1999). Indeed, Speer et al. (2000) note that the observed PV gradients along the Upper Circumpolar Deep Water layer have the correct sign to support the observed poleward eddy transport of this water mass across the ACC.
In the absence of external forcing to maintain a large-scale PV gradient (e.g., within a closed streamline), eddies can act to homogenize PV on a buoyancy surface, an idea that has long been attractive in both oceanographic (Rhines and Young 1982) and atmospheric (Sun and Lindzen 1994) communities. Marshall et al. (1993) first advanced this idea for the Southern Ocean; however, subsequent results suggest that PV is far from uniform, as evidenced by strong variations in the e-folding depth scale of the stratification across the Southern Ocean (Karsten and Marshall 2002; Marshall and Johnson 2017). On Jupiter and Saturn, we must reckon with inhomogeneous PV; strong alternating zonal jets render the relative vorticity dominant over the planetary vorticity, producing a total vorticity that varies greatly and alternates sign, causing the northward PV gradient also to alternate sign (Read et al. 2006, 2009a).
The degree of PV homogenization, and the shear stability properties of the flow, can be assessed by the PV change over a given distance or, more usefully, over a given change of Bernoulli potential B. If PV is nearly homogeneous, the slope of the Q versus B relationship should be near zero. However, this statement is meaningless until this slope is suitably nondimensionalized. An idea emerging from Jovian atmospheric dynamics is that the ratio of the local flow speed to the intrinsic Rossby wave speed provides this nondimensionalization (Dowling 2014, 2019). As the former is related to gradients of B and the latter to gradients of Q, this ratio is proportional to the Q versus B slope. This echoes the Mach or Froude numbers, which are the predominant tools in aeronautics or hydraulics to study acoustic waves or gravity waves, respectively. One might expect the analogous ratio for Rossby waves to be equally important, but it has not yet reached that status; this work is partly an attempt to move in that direction.
Versions of this ratio for Rossby waves have been alternatively called the Froude/Rossby number (Armi 1989), but we follow nomenclature from Dowling (2014), calling this ratio the “Mach” number and denoting it “Ma” (in quotes, understanding we are not concerned with sound waves). We also use “subsonic” and “supersonic” (again in quotes) to denote 0 < “Ma” < 1 and “Ma” > 1, respectively. A third possibility arises for Rossby waves: “Ma” < 0, corresponding to waves traveling downstream rather than upstream. This negative regime exists because Rossby waves are unidirectional (in QG, their intrinsic phase velocity is perpendicular to, and to the left of, the PV gradient), whereas sound waves are omnidirectional, so the traditional Mach number is nonnegative.
We have reviewed that shear stability and PV homogenization are related to PV gradients, which are usefully measured by the Q versus B slope, which nondimensionalizes into the “Mach” number. Therefore, it should not be surprising that Arnol’d’s (1966) shear stability theory can be re-expressed using “Ma” (Dowling 2014). This closes a thought loop, and reinvigorates shear stability theory with old ideas from other branches of fluid mechanics: shear instability can be thought of as the shock of Rossby waves that occurs when “Ma” descends through unity at a critical line (Dowling 2014).
b. Arnol’d’s stability theorems and the “Mach” number line
Arnol’d’s (1966) first and second stability theorems are arguably the preeminent shear-stability theorems. The first includes the classic theorems of Rayleigh (1880), Kuo (1949), Fjørtoft (1950), and Charney and Stern (1962) as special cases.
Arnol’d’s theorems apply in several settings, the most pertinent of which for ACC dynamics is a continuously stratified quasigeostrophic fluid (McIntyre and Shepherd 1987, their appendix B), in which the instability abated by Arnol’d’s conditions is baroclinic. However, for simplicity we follow McIntyre and Shepherd (1987) and Dowling (2014), reviewing the case of a single-layer, unforced QG flow, for which the class of instability considered by Arnol’d’s conditions is barotropic. While barotropic instability is not unimportant to ACC dynamics (Youngs et al. 2017), the main point to note here is that Arnol’d’s conditions (and the forthcoming “Mach” number conditions) apply equally well to barotropic and baroclinic settings.
The flow considered is governed by
where ψ is the streamfunction and q the QG PV. If the flow is steady, then (1) implies ∇ψ and ∇q are parallel, so ψ and q are functionally related. Considering the simple case where this relationship is single valued, if one of
holds everywhere in the domain, then Arnol’d’s first or second theorem, respectively, ensures the flow is stable (in the sense of Liapunov, meaning perturbations may grow but not unboundedly).1 In (2), is the smallest eigenvalue relating PV perturbations to flow perturbations, i.e., , where primes indicate deviations from the steady background state.
Dowling (2014) re-expressed Arnol’d’s stability theorems in terms of “Ma.” The QG intrinsic Rossby wave phase speed for the gravest mode is
while the total phase speed is , with u the zonal velocity of the steady background flow. The “Mach” number is defined as
Using ∂ψ/∂y = −u and (3) gives
Thus, Arnol’d’s first branch of shear stability is associated with negative “Ma,” and the second branch with “supersonic” “Ma.” Considering the inverse of “Ma,” the two stability branches unite within “Ma”−1 < 1. Instability is possible for “Ma”−1 ≥ 1.
Note, u is reference-frame dependent, hence so too are ψ and “Ma.” Nonetheless, (6) holds in any frame. Mathematically, this is because any reference-frame-induced velocity augmenting u in (4) also augments ∂ψ/∂y, so (5) holds. Physically, it is because shocks are reference-frame invariant phenomena. Reference frames will be discussed further in section 2f.
One deficiency of (2a) is that it is dimensional. This problem dates to Rayleigh’s original inflection-point theorem of 1880, and is perhaps more serious than has been generally appreciated. The emerging perspective of “Ma” for shear stability directly resolves this concern.
c. Physical mechanisms
The two traditional viewpoints for the physical mechanism of shear instability, and its relation to critical lines, are counterpropagating Rossby waves and overreflection. These have been largely reconciled (see Harnik and Heifetz 2007, and references therein). The former, involving waves that phase lock and amplify each other for an appreciable time, is better-suited to Mach-number style thinking.
Consider a steady and zonal flow, having zonally uniform q and ψ, and possessing a single maximum of q at some critical latitude. The Charney–Stern stability theorem, here asserting stability if the meridional PV gradient is single-signed, is violated. Yet, instability is not assured. From (3), the intrinsic Rossby wave phase speed is westward to the south of this critical latitude, and eastward to the north. If ψ is also a maximum at the critical latitude, then dq/dy and dψ/dy are both positive to the south of the critical latitude and both negative to the north, so dψ/dq > 0 everywhere, so the flow is stable by (2a). The zonal flow is in the same direction as the intrinsic phase speed of Rossby waves: “Ma” < 0. Wave disturbances swim downstream and are ripped apart from each other, so cannot coamplify. Alternatively, (2b) can be satisfied when ψ has a minimum at the critical latitude, provided dψ/dq is sufficiently negative. Now the zonal velocity is eastward to the south of the critical latitude and westward to the north, opposite to but faster than the intrinsic phase speed of Rossby waves: a “supersonic” “Ma.” Wave disturbances try to swim upstream but the flow is stronger; they too are ripped apart, and cannot coamplify. Because is the smallest eigenvalue, these are the fastest such waves around; shorter, slower waves are ripped apart even sooner. See Dowling (2014, their Fig. 5) for schematics of these scenarios. Arnol’d’s two branches of shear stability are united by the physical reality that all Rossby waves are sheared apart from one another.
The physical mechanism for shear instability can be similarly understood when “Ma” = 1: the background flow exactly opposes the waves’ intrinsic phase speed, so Rossby waves are stationary, and can coamplify.
d. The onset of shear instability
Arnol’d’s (1966) theorems provide sufficient conditions for stability, but are silent on necessary conditions; the on–off switch for shear instability has not yet been mathematically found. Mach-number type thinking provides some helpful physical insight. A fundamental tenet of fluid dynamics is that a shock must form when the flow speed drops from above to below a wave speed. For sound waves, this is a sonic boom as the (acoustic) Mach number descends through unity. For gravity waves it is a hydraulic jump or bore as the Froude number descends through unity. This follows from the general theory of characteristics applied to hyperbolic partial differential equations, i.e., those supporting waves (Whitham 1974, chapter 2). Dowling (2014) conjectures: for Rossby waves, the shock that forms when “Ma” descends through unity at a critical line is the onset of shear instability itself.
We speculate in favor of this conjecture, adding support as follows. Where 0 ≤ “Ma” < 1, long Rossby waves shear apart, but instability might arise if short waves, being intrinsically slower, could phase lock and coamplify. However, short Rossby waves are dispersive. Their different phase and group speeds means that short Rossby waves cannot phase lock for much time: wave crests that are coamplifying will soon reach nodal points in the wave packet’s envelope, killing any coamplification up to that point. Indeed, both the (acoustic) Mach number and the Froude number deal with nondispersive waves. Likewise, it is nondispersive long waves (Hide 1969, p. 846) that underpin the definition of a useful “Ma” for Rossby waves.
e. Consideration of eddy mixing parameterizations
Since “Ma” is a key nondimensional parameter for shear instability, it may arise in theories describing associated processes, such as lateral (isopycnal) mixing. Improving lateral mixing parameterizations improves dynamical balances in global simulations (e.g., Danabasoglu and Marshall 2007; Gent and Danabasoglu 2011; Stanley and Saenko 2014; Farneti et al. 2015; Poulsen et al. 2018; Saenko et al. 2018). Recent advances prescribe the lateral eddy diffusivity by mixing length theory and account for a kinematic suppression of the mixing length by a mean flow that advects tracers past eddies (Ferrari and Nikurashin 2010). Including these effects in models has improved estimates of the eddy diffusivity (Klocker et al. 2012). Bates et al. (2014), for example, express the lateral eddy diffusivity as
where Γ ~ 0.35 is a mixing “efficiency,” b1 ~ 4, and the eddy mixing length Leddy is suppressed by the factor in braces. Couching this suppression factor in terms of , expressed using the background velocity , yields
Notice that (8) is insensitive to the sign of “Ma”−1. This follows from the purely kinematic perspective of the aforementioned mechanism for the suppression of eddy diffusivity by the mean flow. But there is an alternative, dynamic perspective, in which the lateral eddy diffusivity depends on critical levels at which the Rossby wave speed matches the flow speed (Smith and Marshall 2009; Abernathey et al. 2010). The parameterization (7) treats the regime “Ma”−1 ≥ 1, where shear instability is possible, the same as “Ma”−1 ≤ −1, where shear stability is guaranteed. Though beyond the scope of this study, an intriguing possibility is that an eddy diffusivity that respects these stability branches, and in particular is sensitive to the sign of “Ma”−1, might further improve the parameterization of isopycnal mixing. Moreover, Youngs et al. (2017) present evidence that barotropic instability is also active in the ACC’s standing meanders, and submit that current eddy parameterizations, being insensitive to the background shear, cannot represent this process adequately. This may offer further support to root an eddy parameterization in terms of “Ma,” since the “Ma” concept is a local one, and is ignorant of the reservoir of energy for instability.
f. Reference frames
PV is reference-frame invariant, but the streamfunction is not. If one finds a reference frame wherein ψ and q are functionally related and satisfy (2), then the flow is stable. Physically, this reference-frame freedom exists because Rossby waves can coamplify not only when they are stationary, but also when they move together with the same total angular phase speed—equivalently, they are stationary in some common reference frame.
For a simple example, suppose the PV and flow are zonally uniform, hence q and ψ are functionally related. If q increases monotonically with latitude (dq/dy > 0), the flow is stable by Charney–Stern, but it may satisfy neither (2a) nor (2b) in the planet’s reference frame: if the zonal velocity reverses direction, then dψ/dy is positive and negative at different latitudes, so dψ/dq is not single-signed. But in a rapidly rotating reference frame aligned with the planet’s rotation axis, the zonal velocity becomes single-signed, as do dψ/dy and dψ/dq, so the flow is stable. The two choices—rapidly rotating in either direction—ensure stability by the two different branches of Arnol’d’s stability theorem.
Saturn and Jupiter present more complicated examples: relative vorticity in the banded zonal jets can overwhelm the planetary vorticity, juxtaposing regions where dq/dy < 0 with regions where dq/dy > 0. This violates the Charney–Stern stability criteria, yet the flow may still be stable by Arnol’d if the zonal velocity vanishes (dψ/dy = 0) at the critical latitudes where dq/dy = 0, in such a way that dψ/dq remains finite and single-signed. If there is only one critical latitude, there exists a reference frame in which the zonal velocity vanishes there. In the presence of many critical latitudes, something special must occur if such a reference frame exists. Dowling (1993) showed this occurs on Jupiter; Stamp and Dowling (1993) further argued that Jupiter sits on the threshold between instability and stability of Arnol’d’s second branch. On Saturn, Read et al. (2009b) noticed the zonal angular velocity is nonzero, but nearly uniformly so, at the critical latitudes in Saturn’s official reference frame. Hypothesizing that Saturn also sits on this threshold, they proposed a new reference frame—a new rotation rate for Saturn—wherein the zonal flow vanishes at the critical latitudes.
Earth is no aquaplanet, and its rotation rate is accurately known: shall we use Earth’s reference frame to evaluate “Ma”? Since Arnol’d’s theorem relies on the flow being steady (and ideal and unforced), one might seek a reference frame in which the flow is steady. (We will use such a reference frame for other purposes in section 3b.) In the realistic ocean, though, no reference frame yields steadiness everywhere. Moreover, we are not actually attempting to apply Arnol’d’s stability theorems, but merely to study their quantity of central interest, re-expressed as “Ma.” Instead, Mach-number thinking supports the use of Earth’s reference frame, as follows. Shocks are physical realities with profound consequences for the flow, and are reference-frame invariant. In aeronautics, a shock is generated not because an object exceeds the speed of sound, but because there is a boundary layer around the object where the fluid velocity differs from the fluid velocity elsewhere by more than the speed of sound. The velocity difference is reference-frame invariant; the object’s speed is not. For us, measuring the velocity in Earth’s frame provides a velocity difference between the local velocity and the velocity adjacent to the (local) seafloor, which is stationary in Earth’s frame. Indeed, the waves we study must fit into their environment, which includes bottom topography because long Rossby waves propagate over long vertical distances (Charney and Drazin 1961). Hence, we measure the velocity in (4) in Earth’s frame; this provides a useful “Ma.”
3. Empirical methods
Motivated by Arnol’d’s nonlinear stability theorems in quasigeostrophy, we now leap to the primitive equations, studying the relationship between the Ertel PV Q and the Bernoulli potential B and their nondimensionalized relationship expressed using “Ma.” This is a leap because no analog for Arnol’d’s second stability theorem has been found for the 3D primitive equations. However, it is not unreasonable to suspect Arnol’d’s stability theorems are relevant to the ACC governed by the primitive equations. We will not attempt to apply Arnol’d’s theorems: the stability criteria will undoubtedly be met in some places and violated elsewhere. Nonetheless, “Ma” is an interesting diagnostic to study, since practical evidence in the QG setting suggests instability often occurs near to where the stability criteria are violated (Waterman and Jayne 2012; Waterman and Lilly 2015).
a. Primitive equations
Suppose b is a quasi-conservative buoyancy: its material time evolution, Db/Dt = ϖ, is small and due only to similar nonconservation of salinity S and potential temperature θ. If b is monotonic with depth, it can serve as the vertical coordinate, and the hydrostatic Boussinesq equations become (Young 2012)
with subscripts denoting partial derivatives. In sequence, these equations describe the zonal momentum (u the eastward velocity), meridional momentum (υ the northward velocity), hydrostatic balance, and continuity (σ the thickness). The depth of the buoyancy surface ξ is an independent variable. The thickness is , with g the gravitational acceleration. The Montgomery potential m solves the coupled equations and , where p is the pressure and ρc is the Boussinesq reference density. Tildes on the coordinates remind us to take these partial derivatives at constant b. The PV is
with the Coriolis parameter f = 2Ω sinθ, Ω Earth’s angular rotation rate, and θ the latitude. The Bernoulli potential is B = m + K, with K = (u2 + υ2)/2 the specific kinetic energy. External forcing and viscosity are represented by (Fx, Fy). The velocity through b surfaces is equivalent to ϖ.
Because seawater is a multicomponent fluid with a nonlinear equation of state, the ocean possesses no variable that is analogous to potential temperature in a dry atmosphere (McDougall and Jackett 1988), along whose isosurfaces fluid parcels can be exchanged without experiencing buoyant restoring forces—the defining characteristic of neutral surfaces (McDougall 1987). Fortunately, serviceable approximations exist, particularly when restricted to a local region. The best-established buoyancy variable for our purposes—to minimize ϖ while ensuring m exists—in a regional study is b = −gρν/ρc, where ρν is orthobaric density (de Szoeke et al. 2000). We compute orthobaric surfaces from scratch (see Stanley 2018, section 3.2.5 for details), constructing a “virtual compressibility” that approximates the true compressibility as a piecewise linear function of pressure and in situ density, using Eulerian time-averaged data from the winters of 2005 and 2006 (which produce the densest surface waters around Antarctica in SOSE). We report in terms of σν ≡ ρν − 1000 kg m−3, and frequently drop the units. For future studies, modified topobaric surfaces (Stanley 2019a) are recommended, being highly neutral while also possessing an exact geostrophic streamfunction m (Stanley 2019b).
Suppose there is a steady () background flow that solves the inviscid, unforced (Fx = 0 = Fy), and buoyancy conserving (ϖ = 0) version of (9). Then (9d) reduces to , and there exists a streamfunction Ψ such that
The gradients of B and Ψ are aligned, so a functional relationship exists between B and Ψ. The same is true for Ψ and Q [cross differentiate (12)], and Q and B [divide (12) by Q then cross differentiate]. The latter is our focus. In general, a given isovalue of B may possess many disjoint contours, along each of which there can be a different value of Q, and thus Q is actually a multivalued function of B, and vice versa (Stanley 2019a). As a first step, this paper treats these functional relations as single valued: though eddies can cause some contours of B to close, the ACC is typically dominated by one (circumpolar) contour for any given B. It is a forgiving region in this regard.
Since atmospheres and oceans are unsteady and forced, we do not expect an exact functional relationship between the instantaneous Q and B. Nonetheless, scatterplots of Q versus B (section 4a) do reveal strong correlations, suggesting an underlying functional relationship overlain with some scatter: a background state atop which unsteady and forced dynamics act (and, ultimately, adjust the background state). By empirically assessing this functional relationship, we learn about shear stability and PV homogenization in the ACC, provided we can appropriately nondimensionalize this relationship into a “Mach” number—which requires estimating the Rossby wave phase speed.
b. Rossby wave phase speed
The embedding of Arnol’d’s stability theorems into the “Mach” number line (section 2b) required that the intrinsic long Rossby wave phase speed (which is also their group speed) is proportional to the PV gradient, as in (3). For a first pass at estimating “Ma” in the atmosphere, precursor studies (Du et al. 2015; Dowling et al. 2017) effectively used this single-layer QG result (3) by analogy as
with Γ = σ0(NH/f)2, where σ0 is a mean layer thickness that relates the Ertel and QG PV gradients by (Dowling 2014), N is the buoyancy frequency, and H the pressure scale height; together, NH/f approximates the deformation radius Ld (Chelton et al. 1998). Du et al. (2015) took both N and H from the definition of the U.S. Standard Atmosphere, whereas Dowling et al. (2017) used regional averages of the instantaneous (NH)2 field for Mars. Though good enough to begin quantifying atmospheric PV homogenization, the approximation Ld ≈ NH/f is valid only for a fluid at rest; velocity shear strongly impacts the Rossby wave phase speed in the ocean (Killworth et al. 1997). We desire a more accurate method to locate the threshold to instability at “Ma”−1 = 1 in the oceanic setting.
Rossby waves in a continuously stratified and sheared fluid are typically decomposed into an infinite set of orthogonal vertical modes with specific phase speeds. Each vertical structure and phase speed is an eigenfunction and eigenvalue of an ordinary differential equation, obtained via separation of variables from a partial differential equation for the time evolution of PV, typically in the planetary geostrophic or quasigeostrophic setting. Appendix C presents this derivation in buoyancy coordinates, then demonstrates two problems with calculating the phase speed of a given normal mode: dependency upon the choice of vertical coordinate (depth or buoyancy), and whether the bottom boundary condition specifies vanishing normal or horizontal flow.
A third, more critical, problem is that this approach determines phase speeds of vertical modes, but we require the Rossby wave phase speed as a function of depth, or more specifically, at the depth of a particular buoyancy surface. Selecting any one vertical mode is clearly incorrect where our buoyancy surface intersects the zero crossing of this mode: this point is not perturbed at all by that wave mode, though it is by other modes. The phase speed at the depth of our surface arises from nonlinear (and nontrivial) interactions between modes of varying power. Rather than pursue this line of attack, we step back, leaving the wave amplitude as a function of depth (or buoyancy), and seek the phase speed directly on a buoyancy surface.
c. Eddy phase speed
In lieu of a theory for the Rossby wave speed on a buoyancy surface, we develop and test an empirical method for a related quantity, namely, the phase speed of eddies on a buoyancy surface. By “eddy,” we refer generally to any PV anomaly; this could be a closed-core vortex, part of an idealized monochromatic Rossby wave, or anything in between.
Briefly, our method finds a reference frame wherein PV anomalies on a given buoyancy surface are closest to stationary; the velocity of this reference frame relative to Earth’s frame estimates the eddy phase speed in Earth’s frame. The method only uses information from the buoyancy surface in question, and takes inspiration from how Read et al. (2009b) obtained Saturn’s rotation period.
If applied to a hypothetical fluid exhibiting a perfectly monochromatic Rossby wave field, this method would trivially discover a reference frame in which PV anomalies are perfectly stationary, and this yields the true, and in this case unique, Rossby wave phase speed. In reality, the ocean presents a complicated and evolving wavefield, and there is no single, unique phase speed to speak of. Likewise, there is no reference frame in which PV anomalies are perfectly stationary. There is often, however, a dominant mode of variability, with a dominant propagation speed. Our method captures this by selecting the “best” reference frame (in a statistical sense to be specified shortly).
The propagation speed of this dominant mode of variability is itself closely linked to a Rossby wave phase speed. Indeed, the distinction between vortices and waves in terms of phase propagation is slight (Polito and Sato 2015). Empirically, it is well known that vortical eddies tend to propagate at the group speed (which is also the phase speed) of long baroclinic Rossby waves (Chelton and Schlax 1996), provided that speed correctly accounts for nonuniformity of the background mean flow (Killworth et al. 1997), bottom topography (Killworth and Blundell 1999), and their combination (Killworth and Blundell 2003a,b). Therefore, one might reasonably interpret the eddy phase speed found by our method as a “pseudo-Rossby wave speed”: the propagation speed it empirically estimates is closely—albeit not entirely straightforwardly—related to a Rossby wave speed.
We provide three further lines of evidence to support this interpretation and verify that this reference-frame shift method produces reasonable results. First, appendix B examines Hovmöller plots that reveal eddy propagation on orthobaric density surfaces, and directly verifies our method’s results against propagation speeds inferred from the Radon transform. Second, appendix C calculates the phase speed of the most unstable Rossby wave vertical mode using four classical methods, against which our method’s results agree about as well as the four methods’ results agree with one another. Third, in a 1.5-layer model (where only one vertical mode exists), appendix D verifies that our method compares well against a direct calculation of “Ma,” and hence also accurately estimates the Rossby wave phase speed.
For our present purposes, it is unimportant whether this dominant mode of variability is conceived as a wave or a vortex: a PV eddy (wave or vortex) that tilts upstream across a shear at a critical line will experience an amplifying self-induced circulation. While the “supersonic” condition relentlessly topples this configuration, the “subsonic” condition relentlessly sets it up (Dowling 2014).
The method considers a reference frame rotating about Earth’s axis with a period of 2π/(Ω + α). When α > 0, this observer’s rotation exceeds Earth’s; objects at rest from Earth’s perspective (e.g., topography) appear to move westward. The primitive equations in this frame are akin to those in Earth’s frame, (9), but subtly modified. See appendix A for details. A velocity v in Earth’s frame becomes , where is the eastward unit vector,
and R is Earth’s radius. Similarly, Eulerian measurements now sample variables from a location moving zonally through Earth’s latitude–longitude coordinates. The Eulerian time evolution of σ, for example, was in Earth’s frame, but becomes . Also, m is augmented by −αA, where
is the specific axial angular momentum. Thus, B becomes
Crucially, Q is reference-frame independent, as are ϖ, Fx, and Fy.
As before, a steady background state of the inviscid, unforced, and buoyancy-conserving equations in the shifted frame exhibits a functional relationship between Q and B(α). In reality, unsteadiness, forcing, and nonconservation of buoyancy cause scatter atop this underlying relationship. There may be less scatter in Q versus B(α), if α is well chosen, than in Q versus B. The α variance of B gives us room to manipulate its relationship with the α-invariant Q. What is the meaning of α?
Suppose Fx = 0 = Fy and ϖ = 0, and the flow is steady in a reference frame shifted by α. Then there are no traveling waves; anomalies (eddies) are stationary in this new reference frame, meaning c + u(α) = 0, where c is the zonal velocity of anomalies in Earth’s reference frame. Therefore, we interpret
as the eddy phase speed, measured in Earth’s reference frame. The corresponding intrinsic zonal phase speed is
In practice, we empirically determine α by minimizing the error in the functional relationship between Q and B(α). We will show that, even in Earth’s frame, Q and B are nearly linearly related in the Southern Ocean, so we fit a linear model with coefficients q0 and q1:
As discussed, Q and B(α) may be more tightly related, as some unsteadiness in Earth’s frame is simply anomalies moving zonally, so we also fit a multilinear model with coefficients b0, b1, and α:
d. “Mach” number
Having identified various practical methods for obtaining the eddy phase speed, we can now estimate “Ma.” The precursor studies (Du et al. 2015; Dowling et al. 2017) made a rough estimate2 of Γ in (13), so that (4) yields
having also used the background state relationship from (12b).
For the reference-frame shift method, we use (18) and the background state relationship in the shifted frame, to get . Invoking the chain rule gives
An alternate invocation of the chain rule gives
where q1 = dQ/dB and b1 = dB(α)/dQ are the model coefficients from (19) and (20). Thus, “Ma”−1 is the ratio of the slope of the Q versus B relation in Earth’s frame to that in the frame moving with Rossby waves/eddies. If one desires Γ, (21) and (23) yield Γ = −b1/(f + ζ). However, (23) is susceptible to L’Hôpital’s trap, where q1 → 0 and b1 → ±∞, so our computations use (22). The form (23) is useful for discussing the following three limiting cases.
First, if Q is homogeneous but B is not, then q1 = 0, and (20) shows B(α) is constant and b1 is arbitrary; thus “Ma”−1 = 0 by (23). Thinking physically, as Q becomes uniform, the restoring force for Rossby waves/eddies tends to zero, so their phase speed and hence “Ma”−1 tends to zero.
Second, if Q is inhomogeneous and a perfect function of A, while B and Q are totally independent of each other, then dB/dQ = 0 but dB(α)/dQ ≠ 0, so “Ma”−1 → ±∞ by (23). Since the Ω term in (15) dominates, A contours are basically latitude circles, while B is independent of latitude. This case involves zonal Rossby waves/eddies propagating on a meridional PV gradient while the background zonal flow is at rest.
Third, when Q is inhomogeneous and a perfect function of B, but totally independent of A, then α = 0 in (20), so B(α) = B and “Ma”−1 = 1 by (22). Thinking physically, if B contours are not perfectly zonal, then a zonal Rossby wave/eddy propagates Q anomalies zonally, off the B contour, supporting the dependence of Q with latitude and diminishing its correlation with B. When α = 0, long Rossby waves/eddies cannot transport Q anomalies this way, thus Q must be constant along B contours.3
Using controlled shallow water experiments, appendix D verifies that this reference-frame shift method produces reasonable estimates of “Ma.”
For each 5-day archived state of SOSE’s 6-yr run, we compute the pressure, layer thickness, PV, and Bernoulli potential on nine orthobaric density surfaces: σν = 27.1, 27.2, 27.3, 27.4, 27.5, 27.6, 27.7, 27.75, and 27.8. (Surfaces that are much shallower or deeper are not circumpolar, year-round.) Figure 1 shows the depths of these surfaces at 240°E on 22–26 March 2006, a representative time step in late summer when orthobaric surfaces are shallowest. Isolines of PV at 240°E appear closely aligned with orthobaric surfaces—but such inspection is insufficient to establish PV homogeneity or heterogeneity.
To make contact with the global or extratropical analyses of precursor studies (Du et al. 2015; Dowling et al. 2017), we first study orthobaric density surfaces spanning the entire ACC: section 4a qualitatively examines the structure of the Q versus B relationship, prior to quantitative regression analyses in section 4b. Then, section 4c translates these quantitative techniques into the local setting to study how the Q versus B slope, the eddy phase speed, and “Ma”−1 vary across the Southern Ocean.
a. Histograms of potential vorticity versus Bernoulli potential
We predefine a region of interest (ROI) for the ACC by combining two regions. The first region is everywhere the 2005–10 mean sea surface height satisfies , where and are evaluated adjacent to the Antarctic Peninsula and Cape Horn, and , The second region is everywhere the sea surface geostrophic zonal velocity exceeds +5 cm s−1. The final ACC ROI is the connected component of the union of these two regions that contains Drake Passage, then flood-filled to include any enclosed holes. This ROI, shown in Fig. 6c, has a southern boundary similar to that of the ACC defined by Orsi et al. (1995), namely, the southernmost circumpolar streamline; its northern boundary closely hugs the northern limit of Drake Passage, but excludes slow moving flows far north of the ACC. The results are not sensitive to small changes in this ROI definition.
Figure 2a shows Q versus B in this ROI on all orthobaric surfaces during 22–26 March 2006. All but the deepest layers exhibit high |Q| on the low Bernoulli (south) side caused by high stratification where the layer lies just below the mixed layer (Fig. 1). We will focus on the ocean interior, where forcing is small, enabling a fairly tight relationship between Q and B.
Two obvious ways to eliminate this near-mixed layer signal are to remove data below some pressure threshold, or below some Bernoulli threshold. The latter is more attractive since it does not rely on a third variable. We attempted this using the Multivariate Adaptive Regression Splines algorithm (Friedman 1993), as implemented in ARESLab (Jekabsons 2016) to automatically minimize the error in Q as a piecewise linear function of B with one free, internal knot, which gives the B threshold. This produces seemingly well-chosen thresholds, but in the final analysis, “Ma”−1 is rather sensitive to the threshold. With a pressure threshold, however, the results are much less sensitive to the chosen value, provided it is not too shallow. Charney-type instabilities are surface trapped and have a (self-chosen) depth scale typically shallower than 300 dbar, though this can reach 500 dbar in some isolated regions (Smith 2007). We remove data shallower than a threshold of 300 dbar, determined by trial and error, to generally ignore these shallow instabilities.
Figure 2b shows Q versus B histograms with this threshold applied, enabling closer scrutiny of the fairly linear relations between Q and B. There are subpopulations in the Q versus B data, most noticeably on σν = 27.2 and 27.3: for a given B in the ACC, there is one low and one high |Q| cluster, corresponding to relatively constant |Q| along B streamlines in the Atlantic and Pacific sectors, respectively. In between, there is a transition region from Kerguelen to Campbell Plateau where |Q| rises along streamlines. In the other transition region, just downstream of Drake Passage, |Q| drops rapidly along streamlines, consistent with analyses of time-mean PV along time-mean streamlines (Gille 1997; Stanley 2018).
b. Q versus B regressions spanning the Antarctic Circumpolar Current
We now perform a comprehensive, primitive vorticity–streamfunction correlation for the ACC by fitting (19) within the ROI and below the pressure threshold, for each 5-day average. Since Q and B are interdependent variables, we use geometric mean regression, minimizing the sum of the areas of the right-angled triangles with edges parallel to the B and Q axes, formed by each data point and the fitted line. Our regression methods follow Greene (2015) and are detailed in the supplemental material.
In Fig. 2b and elsewhere, we use the Nash and Sutcliffe (1970) model efficiency coefficient (NSE) to quantify the fits’ quality. The NSE is 1 − (MSE/VAR), where MSE is the mean-squared error of the predicted (modeled) variable, and VAR is the variance of the observed data. For geometric mean regression, MSE is the mean area of the aforementioned right-angled triangles, and VAR is the covariance of the two variables. Positive NSE indicates the prediction is better than the mean of the observations, and negative indicates it is worse. The NSE takes values in (−∞, 1], with 1 indicating a perfect fit.
The slope of the Q versus B relation (19) is dQ/dB = q1, shown in Fig. 3a. All layers have a positive slope—generally B increases northward, as does Q (recall Q < 0)—corresponding to “Ma”−1 > 0. Deep layers (σν ≥ 27.6) have a smaller slope, and less interannual variability. Middle layers (27.3 ≤ σν ≤ 27.7) exhibit a gradual decrease in slope through the 6 years, less so for the deeper of these layers. The seasonal cycle is evident in the slope for the shallowest layer (σν = 27.1), with growth in winter and decline in summer; the maximum occurs in late winter 2008; other late winters exhibit smaller maximums, excepting 2005 when growth continues into 2006. The second shallowest layer (σν = 27.2) shows little secular trend, perhaps counterbalancing growth in the layer above and decay in those below. There is substantial coupling between layers, both in the high-frequency variability and in the seasonal cycle in the shallow layers. In summer 2005, the five shallowest layers exhibit a very similar slope of Q versus B.
Next, we determine the reference-frame shift α by fitting the model (20), using geometric mean regression with Q and B as dependent variables, and A as an independent variable, which minimizes the sum of the areas of the right-angled triangles with edges parallel to the B(α) and Q axes. Figure 4 shows an example: σν = 27.3 during 26–30 August 2005. With α = −4.19 × 10−9 s−1 (corresponding to a westward speed of 1.53 cm s−1 at 55°S), the Q versus B(α) relation tightens considerably. With this α, “Ma”−1 = 1.48, among the largest we discovered.
Figure 3b shows c = αR cos(−55°), the eddy phase speed diagnosed at a representative ACC latitude of 55°S. The major late winter maximums of dQ/dB for σν = 27.1 are evident here as more westward eddy phase speeds. The eddy phase speed is generally westward (discussed further in section 4d) for most surfaces, except the two deepest circumpolar layers σν = 27.75 and σν = 27.8, and occasionally for other layers, particularly in summer 2005.
Having α and thus B(α) from (16), Fig. 3c shows “Ma”−1 from (22). We use orthogonal regression to fit B and B(α). This minimizes the sum of the squared orthogonal distance between each data point and the fitted line (and thus requires common units between the variables). Only the two deepest circumpolar layers have “Ma”−1 < 1 most of the time (corresponding to α > 0 most of the time). The deepest layer maintains “Ma”−1 > 0.5 and the second deepest generally maintains “Ma”−1 > 0.75, so even on these deep layers the PV should not be regarded as homogeneous (“Ma”−1 ≈ 0). For all other layers, “Ma”−1 generally hovers between 0.75 and 1.25—the range aerodynamicists call “transonic”—right around the threshold of instability and clearly indicating inhomogeneous PV. Winter 2005 on σν = 27.2 exhibits the largest “Ma”−1, just shy of 1.5. The strong secular trend of the Q versus B slope is not present in “Ma”−1. This indicates, through (21), strong secular changes in , probably owing to strong changes in the stratification of the extratropical Southern Ocean (Cerovečki et al. 2019). Subannual variability of “Ma”−1 is also evident. The source of this variability lies in the deep interior, not just the relatively shallow part of each orthobaric surface near the threshold of p = 300 dbar. This variability is substantially similar, if not somewhat enhanced, when this analysis is repeated but excludes all data shallower than p < 600 dbar (not shown). Part of this variability is likely due to barotropic dynamics adjusting the zonal flow speed. Indeed, there is some resemblance to the total ACC transport (Fig. 3d). Though σν = 27.1 exhibits the largest slope of Q versus B, the largest “Ma”−1 are found on σν = 27.2 and 27.3, whose depth on the northern ACC flank is roughly 1000 m (Fig. 1), suggestive of baroclinic instability whose energy source, the available potential energy, has a depth scale of order 1000 m in the Southern Ocean (Smith 2007).
Figure 5 provides a different view, showing “Ma”−1 as vertical profiles of σν, averaged over 30 days to avoid aliasing the high frequency variability particularly evident in the lower layers. One profile is shown for each year in late winter (29 August–27 September), when “Ma”−1 > 1 in the upper layers and “Ma”−1 < 1 in the deepest layers, generally. The layer at which “Ma”−1 crosses unity varies in time but can be quite deep, between σν = 27.6 and σν = 27.75—about −1500 m and −2500 m at 55°S (Fig. 1). The profiles’ shape also varies in time: the deepest and shallowest “Ma”−1 values do not change drastically, but the interior layers (σν = 27.2–27.5) generally exhibit “Ma”−1 decreasing from one winter to the next, though puffing out to larger “Ma”−1 in 2005 and 2009. Below σν = 27.6, “Ma”−1 decreases quite smoothly. Additional profiles reveal the seasonal cycle for one year, 2008: throughout the water column, “Ma”−1 is maximum at the winter solstice and minimum at the summer solstice, except for σν = 27.1 where “Ma”−1 is largest in late winter. At the summer solstice, “Ma”−1 hovers close to unity from σν = 27.6 and shallower. (Prior to averaging over 30 days, we found “Ma”−1 on 17–21 December 2008 hovers even closer to unity here.) The range σν = 27.2–27.7 roughly bounds Upper Circumpolar Deep Water, which is known to have relatively strong lateral PV gradients, particularly at a neutral density of 27.8 (Speer et al. 2000), roughly equal to our σν = 27.5. But strong lateral PV gradients do not necessarily nondimensionalize into large values of “Ma”−1. In late winter, there does appear to be a marked drop of “Ma”−1 from σν = 27.2 to 27.3, but σν = 27.5 does not particularly stand out here, nor in the perspective of the Q versus B slope (Fig. 3a).
A tight functional relationship between Q and B is not a priori guaranteed. Even though the ACC exhibits essentially one contour per isovalue (of B or Q) implying a functional relationship between Q and B should be essentially single-valued, significant scatter overlying this could be caused by large PV tendency or forcing. However, we find an impressive relation does exist over the whole ACC. Moreover, Q and B(α) are even more tightly related, implying their gradients are more tightly aligned. Since forcing is identical between the two reference frames, this indicates that the magnitude of the PV tendency in the shifted frame is less than that in Earth’s frame. That is, the PV tendency associated with zonal propagation of anomalies in Earth’s frame is a significant term in the instantaneous PV budget, suggesting that mechanical and buoyancy forcing of PV is less significant. This conclusion is supported by time-mean PV analyses by Gille (1997) and Stanley (2018): eddy forcing of PV is a dominant term, whereas mechanical and buoyancy forcing of PV are unimportant except very near topography. (Eddy forcing in the time-mean is related to variability of the instantaneous PV, and thus to the instantaneous PV tendency.) As a whole, the ACC exhibits “Ma”−1 > 0 for all circumpolar orthobaric surfaces. But might that arise as a mixture of homogeneous PV in some locations (e.g., deep northern regions) and strong PV gradients in other locations (e.g., shallow southern regions)?
c. Local Q versus B regressions
We now apply the above methods to localized data, producing spatial maps of the Q versus B slope, the eddy phase speed, and “Ma”−1. For each water column, we use 325 data points (fewer near land) from a rectangular window, (25/6)° zonally by (13/6)° meridionally. At 60°S, this is 300 km zonally and 240 km meridionally, substantially larger than the deformation radius (roughly 15 km). The regressions are weighted by , with d the great circle distance to the central grid point. Adjusting the weights and size of the window does not significantly impact the results.
Figure 6 maps Q, B, and the slope dQ/dB, for σν = 27.3 on 22–26 March 2006. The predominantly circumpolar nature of Q and B contours is visible, though there are plenty of eddies that manage to create closed contours, particularly of Q. In most of the ACC, dQ/dB > 0 (both Q and B increase northward), but there are many patches where this slope is negative. A notable example is upstream of Drake Passage in the ACC’s southern flank, where mode water formation has injected filaments of near-zero PV water, making patches where , and hence dQ/dB < 0. In the ocean basins north of the ACC, this slope is comparable in magnitude to that in the ACC, but roughly as often negative as positive. The flow being weaker here than in the ACC, B variations are smaller, and so too are Q variations. At the southern flank of the ACC (roughly where the surface rises above 300 dbar—not shown) the slope increases drastically in magnitude.
We now switch to use ordinary least squares (OLS) regression in the model (20), with B the dependent variable. With OLS, Q should not be the dependent variable: the results would be extremely sensitive to the input data, since B and A are often tightly correlated (multicollinearity). Since Q has more small-scale structure, Q and A are less correlated, and OLS with B as the dependent variable is more robust, albeit still less robust than geometric mean regression. The present disadvantage of geometric mean regression, though, is that it requires solving a quadratic equation, so that the error-minimizing solution can swap between roots as the input data vary, leading to discontinuities in our subsequent spatial maps. In contrast, the outputs of OLS vary continuously with its inputs.
Figure 7 maps the total and intrinsic eddy phase speeds from (17) and (18). The intrinsic eddy phase speed is generally, but not entirely, westward.4 There is a band of slightly faster eastward propagating eddies on the north flank of the ACC, especially in the Indian sector, perhaps related to a weakening of the PV gradient. From this local perspective, the total eddy phase speed seems more eastward than westward, consistent with downstream advection by the background flow.
Figure 8a maps “Ma”−1 on 22–26 March 2006. The color map is circular because “Ma”−1 → ±∞ is just “Ma” approaching zero, the boundary of Arnol’d’s first stability criterion. The ACC lights up as a region with “Ma”−1 ≈ 1, though far from uniformly so. The ACC is populated with hundreds of mesoscale eddy-sized regions having “Ma”−1 ≈ 1, sometimes enveloping a pocket of “Ma”−1 > 1. An example of this can be seen at 32°S just west of the Cape of Good Hope: an Agulhas eddy. Perhaps these features are mesoscale eddies themselves, though distinguishing such features from filamentary regions of “Ma”−1 ≈ 1 is challenging. Regardless, the smallness of these features relative to a typical Rossby wavelength may somewhat nullify them: Rossby waves/eddies might only coamplify in a small region where “Ma” ≈ 1, such that any growing waves cannot be long waves, and hence are dispersive, and hence cannot coamplify for long (cf. section 2d).
North of the ACC, the ocean basins are, broadly speaking and particularly in the west, dark regions of “Ma”−1 ≤ 0, with occasional strands of larger “Ma”−1. A long feature of “Ma”−1 ≈ 0 (homogeneous PV) appears on the northern flank of the ACC between the Cape of Good Hope and Campbell Plateau, roughly aligned with the Dynamical Subtropical Front (Graham and De Boer 2013). Evidently, a reasonably strong flow exists here, but the PV gradient is weak, perhaps eroded by eddy-induced downgradient PV fluxes. Another such feature appears in the southern flank of the ACC upstream of Drake Passage, related to mode water, as discussed. Similar features appear deeper, on the σν = 27.7 surface (Fig. 9a.)
The time-mean of “Ma”−1, shown in Fig. 8b, appears considerably cleaner than its 5-day view. As mentioned, “Ma”−1 = ±∞ are equivalent, so we form the time-mean using directional statistics. Specifically, we transform each “Ma”−1 into phase angles using the inverse tangent of its negation, double these into the range (0, 2π), then average the unit vectors in the complex plane specified by these phase angles. The resultant vector R1 specifies the time-mean “Ma”−1 merely from its phase angle (reversing the above steps, i.e., −tan[arg(R1)/2]), while the standard deviation, shown in Fig. 8c, is given by .
Certain regions remain near “Ma”−1 ≈ 1 consistently (low standard deviation). In the ACC, these are downstream of the Southwest Indian Ridge, around Kerguelen Plateau and Macquarie Ridge, a rather localized region on the downstream flank of the Pacific–Antarctic Ridge, downstream of Drake Passage, and the Mid-Atlantic Ridge—all matching the hotspots of upwelling identified by Tamsitt et al. (2017) and of topographic form stress identified by Masich et al. (2015), and are likely to be sites where mixed barotropic–baroclinic instability is active (Youngs et al. 2017).
In the ocean basins to the north, the time-mean “Ma”−1 remains near unity off the west coasts of Africa, Australia, and particularly South America, whereas the opposite coasts innocuously exhibit time-mean “Ma”−1 < 1, consistent with eddy generation in the basin’s eastern margin and destruction in the western margin (Zhai et al. 2010). All these features exist on the σν = 27.7 as well (Figs. 9b,c); a notable difference is that σν = 27.7 sports larger patches of time-mean “Ma”−1 > 1 in the ocean basins, though intermittently (large standard deviation).
d. Discussion: Eastward versus westward phase speeds
It is generally accepted that the ACC’s fast eastward flow Doppler shifts the baroclinic Rossby wave speed to be eastward (Klocker and Marshall 2014). Our local method captures this eastward phase speed (section 4c), but our ACC-wide fits (section 4b) diagnosed a westward phase speed. This discrepancy requires some discussion.
It is possible that the ACC-wide result incorporates some regions where propagation is truly westward, and these slightly outweigh other regions of eastward propagation. To test this, we first expanded the rectangular window for the local method to 20° × 10°, then further to 60° × 10°, but the diagnosed phase speed maps remain generally eastward within the ACC. Then, letting the window span all longitudes and 2° meridionally, the diagnosed phase speed varies considerably with latitude. When centered near 45°S, this phase speed matches the weakly westward result for the ACC-wide method. Farther south, it is strongly eastward in the core of the ACC. Moving northward, it is (less strongly) eastward with a peak around 40°S, and is westward north of 37°S, which is about the northernmost limit of our ROI at any longitude. Hence, the westward phase speed estimated for the ROI seems to be somewhat dependent upon the ROI’s nontrivial shape.
It is possible that our empirical methods incorporate, to some degree, the fast westward barotropic Rossby wave; as well they should, since the barotropic wave also influences Q and B on a given surface. However, this does not explain the difference between our ACC-wide and local fits, about which we speculate as follows.
The view of shear instability as coamplifying Rossby waves/eddies requires that they are stationary in some common reference frame. This frame must be fairly close to Earth’s frame: otherwise, the waves/eddies that are stationary in the common frame are moving zonally together in Earth’s frame, so rapidly that they cannot significantly coamplify before they encounter topography and are interrupted. Thus, a planetary-scale disturbance is roughly at rest in Earth’s frame: the westward intrinsic speed of the long Rossby wave roughly matches the background eastward flow. Smaller-scale eddies, however, get swept downstream (eastward) because the intrinsic speed of short Rossby waves is slower than their long cousins. Jupiter nicely illustrates this (Stamp and Dowling 1993), with its Great Red Spot nearly at rest and smaller eddies advecting past it. It appears that these ACC-wide fits have conglomerated hundreds of eddies into a planetary-scale disturbance, with its own speed—close to rest (if marginally stable), or making slow headway upstream (if slightly unstable). In contrast, the local fits capture individual eddies, with their eastward (downstream) speed.
Some theoretical models of the ACC assume that PV is homogeneous on buoyancy surfaces (Marshall et al. 1993; Marshall 1995). However, near alignment of PV and buoyancy surfaces is not enough to justifiably say that PV is homogeneous. For example, Speer et al. (2000) note the presence of large-scale PV gradients on isopycnals in hydrographic data, consistent with the eddy-induced transports of water masses across the ACC. Similarly, Karsten and Marshall (2002) and Marshall and Johnson (2017) show that, moving northward across the ACC on buoyancy surfaces, the layer thickness σ increases by at least a factor of 2, and thus the large-scale PV, f/σ, decreases (in absolute terms) by a slightly greater factor. Still, this is not a proper nondimensional measure of PV inhomogeneity: over what horizontal distance is it relevant that PV changes by a factor of 2? (Also recall that PV is arbitrary up to an additive constant.)
The “Mach” number for Rossby waves gives such a nondimensionalization. Homogenized PV corresponds to “Ma”−1 = 0. This sits at the intersection of Arnol’d’s (1966) first and second shear-stability theorems, which map to “Ma”−1 ≤ 0 and 0 ≤ “Ma”−1 < 1, respectively (Dowling 2014). Instability is possible only for “Ma”−1 ≥ 1.
Dowling (2014) conjectures that the threshold of shear instability is itself the shock of Rossby waves (“Ma”−1 = 1) at a critical level. We have further supported this conjecture by arguing that, although instability is possible by Arnol’d for any “Ma”−1 ≥ 1, when “Ma”−1 > 1 only non-long Rossby waves can phase-lock and coamplify, but their dispersive nature prevents them from doing so for much time.
By estimating “Ma”−1 on orthobaric surfaces in the ACC as a whole, we find that PV is far from homogeneous. Instead, the flow configuration generally hovers around the stability threshold, taking 0.75 < “Ma”−1 < 1.25 for all but the deepest circumpolar layers, where “Ma”−1 ≈ 0.5. The prevalence of “Ma”−1 ~ 1 (“sonic” or “choked” PV) rather than “Ma”−1 ~ 0 (homogeneous PV) suggests shear instability is occurring on interior buoyancy surfaces in the ACC. This process is eliminated from theoretical models assuming homogeneous PV.
Next, we map “Ma”−1 locally on orthobaric surfaces in the Southern Ocean. We find patches of potentially unstable “subsonic” flow (“Ma”−1 ≥ 1) are sprinkled throughout the ACC at any given time, while regions north and south of the ACC more typically inhabit the stable regimes (“Ma”−1 < 1) including some patches of nearly homogenized PV (“Ma”−1 ≈ 0). Most of these patches of “subsonic” flow appear intermittently in time, but there are several topographically constrained “hotspots” in the ACC, such as downstream of Drake Passage and Kerguelen Plateau, where “Ma”−1 hovers consistently near unity throughout SOSE’s 6-yr period. These hotspots correspond to important regions of upwelling (Tamsitt et al. 2017) and bottom form stress (Masich et al. 2015).
A process study by Youngs et al. (2017) examines these hotspots in greater detail, employing the geometric eddy framework (Marshall et al. 2012) to characterize energy exchanges associated with shear instability induced by standing meanders in the ACC. In an eddy-resolving, re-entrant channel model forced by a zonal wind stress, a standing meander forms over an idealized meridional ridge. The input to eddy kinetic energy immediately downstream of the ridge, they find, is roughly balanced between two sources: mean kinetic energy and eddy available potential energy. This is indicative of mixed barotropic–baroclinic shear instability.
Where the instability types are changing and mixing along the stream, as appears to be the case in the ACC, pairing such a Lorenz energy budget analysis with the shear-instability diagnostic “Ma”−1 could prove mutually beneficial. The “Ma”−1 diagnostic identifies when and where shear instability can be born. This instantaneous and local precision refines the averaged view of the eddy geometric framework. However, “Ma”−1 does not specify whether the eddy kinetic energy feeds from the mean kinetic energy reservoir (barotropic instability) or from the eddy available potential energy reservoir (baroclinic instability), which the geometric eddy framework directly addresses.
Our results suggest that instabilities in the ACC could arise, in part, from Rossby waves or eddies phase-locking on individual layers, in particular at certain topographically controlled hotspots. This augments the common picture of homogeneous PV on individual layers with instability arising by interacting edge waves, as in Eady (1949). This picture is consistent, qualitatively, with the findings of Marshall et al. (1999) that relatively subtle changes to the PV gradient in the interior, second layer of a three-layer baroclinic jet can dramatically affect the magnitude of baroclinic slumping throughout the fluid column, suggesting that interior PV gradients act as a source of instability in their own right.
Interior PV gradients might also play an important role in refining our theoretical understanding of “eddy saturation” of the ACC, the relative insensitivity of the ACC volume transport to surface wind forcing (e.g., Straub 1993; Munday et al. 2013). Recent theoretical models have pointed to the critical role of bottom drag and other eddy damping mechanisms (Marshall et al. 2017; Mak et al. 2017) in explaining eddy saturation, the key idea being that the ACC lies in a state of marginal stability wherein the unforced Eady growth rate precisely overcomes the eddy energy damping rate from bottom drag. However, additional sources of instability due to interior PV gradients could significantly modify this picture, in particular if eddies resulting from interior PV gradients have a more confined vertical structure and are less strongly damped by bottom drag.
Since most Q and B contours in the ACC are circumpolar, our ACC-wide analysis treated these variables as single-valued functions of one another. In general, the relationship is multivalued, with a constant but unique Q value on each disjoint contour of B at the same isovalue. To handle this function as multivalued requires a topological analysis of the connected components of level sets of B by means of the Reeb (1946) graph, analogous to how Stanley (2019a) analyzes neutral surfaces. A reference frame rotating with the ACC’s mean zonal velocity, say, reveals hundreds of eddies as closed B(α) contours. Though topologically more complicated, this view may ultimately yield deeper insights about the statistical mechanics of turbulence than can be discovered in Earth’s frame. Indeed, David et al. (2017) found that, in a barotropic channel model, the probability distribution function of Q becomes far cleaner when viewed in a Lagrangian frame of reference.
A functional relationship between potential vorticity and Bernoulli potential is a powerful concept. The mass-weighted streamfunction is tied into a triad of such functional relationships, where any of Q, B, or Ψ can be determined from the other two. This triad completely specifies a solution of the unforced, steady, primitive equations. This enables one to study Rossby waves, by linearizing the unsteady primitive equations about this background state. Classic studies of Rossby waves pick a background state at rest or with a purely zonal flow; with the Q, B, Ψ triad, more realistic background states may be available.
We thank Tomos David, Ed Doddridge, Glenn Flierl, Ryan Holmes, Chris Hughes, Matt Mazloff, Scott Springer, and three reviewers whose comments improved this manuscript. GJS acknowledges partial support from the Clarendon Fund and the Canadian Alumni Scholarship at Linacre College, Oxford, and from the Australian Research Council through Grant FL150100090. DPM acknowledges partial support of the Natural Environment Research Council, NE/R000999/1. SOSE output can be downloaded at http://sose.ucsd.edu.
Equations of Motion in a Shifted Reference Frame
Consider a reference frame rotating about Earth’s rotation axis with period 2π/(Ω + α). If the velocity measured on Earth is v, in this frame it is , where u(α) = −αr cosθ, with r the radial distance from Earth’s center.
To derive the primitive equations in Earth’s reference frame from the Navier–Stokes equations, the centrifugal and gravitational potentials are combined into the geopotential, which motivates a new geometry, the geoid. Thereafter, this nonsphericity is usually ignored (Vallis 2017). We maintain this as the geopotential, to ensure the buoyancy is α invariant, i.e., b(α) = b where b(α) is the buoyancy measured in the new frame.
The total vorticity ωa is also α invariant: ωa = 2Ω + ∇ × v = 2(Ω + α) + ∇ × v(α), where and . This follows from in spherical coordinates.
The PV is the inner product of the total vorticity and the buoyancy gradient—both α invariant. Since the inner product is a coordinate-independent geometric relationship, the PV is also α invariant, as advertised.
In the new frame, the Eulerian buoyancy equation is
where the Eulerian time-evolution of b is , using the fact that is independent of time t and longitude λ. In spherical coordinates, note that , reducing (A1) to the buoyancy equation in Earth’s frame—as it should since the Lagrangian derivative is reference-frame independent.
The continuity equation is
matching that in Earth’s frame, since u(α) is zonally uniform.
The Boussinesq, vector-invariant momentum equation is
As with b, the Eulerian time evolution of v(α) in the new frame is . Forcing, encapsulated by F, is as in Earth’s reference frame. Keeping the original geopotential from Earth’s frame introduces a horizontal centrifugal acceleration in the new frame. In Earth’s frame, the centrifugal acceleration appears in B as −Ω2(r cosθ)2/2. In the new frame, the excess centrifugal acceleration appears in B(α) as [−(Ω + α)2 + Ω2](r cosθ)2/2. Similarly, the excess specific kinetic energy is [v(α) ⋅ v(α) − v ⋅ v]/2. Adding these two excesses to B, the quadratic α terms cancel identically, yielding B(α) = B − αA, with A from (15). In spherical coordinates, one can directly verify , reducing (A3) to its version in Earth’s frame.
Next, (A1)–(A3) are transformed into the primitive equations, applying the thin-shell approximation (swapping r for Earth’s radius R in u(α) and A), the traditional approximation, and hydrostatic balance. The final transformation to buoyancy coordinates proceeds exactly as in Earth’s frame. All told, the new frame’s primitive equations match those in Earth’s frame, (9), but coordinates and α-variant variables are decorated by (α) superscripts.
Eddy Phase Speeds on Orthobaric Surfaces
The phase speed of eddies or Rossby waves is often directly estimated on a buoyancy surface from a Hovmöller diagram of PV on that surface, as shown in Fig. B1. For illustrative purposes, we have selected two regions where a relatively clear wave-like signal emerges above a complicated, evolving wavefield. Other regions, such as upstream of Kerguelen Plateau where we performed additional tests at 50°S, 15°–50°E, from 13 April to 30 October 2005, typically exhibit a more erratic wavefield, but our methods are also capable of handling these regions (not shown). In Fig. B1, the overlain red lines show the phase speed estimated from our reference-frame shift method, which successfully differentiates between the predominantly westward and eastward phase propagation of the two regions (Fig. B1a vs B1b), and generally succeeds at capturing finer-scale behavior.
For comparison, each overlain white line shows the phase speed estimated using the Radon transform of the PV in a rectangle measuring 4.17° zonally and 65 days temporally, centered on that white line. Specifically, the Radon transform rotates this rectangle of data, then integrates horizontally. The vertical integral of the absolute value of these horizontal integrals provides a measure of how well aligned the data are with the rotation angle. We carry out this procedure at every degree of rotation between 0° and 180°, selecting the angle that maximizes this measure. The estimated eddy phase speed is obtained from the tangent of this angle.
Visual inspection suggests that, generally, the two methods agree reasonably well. There are clearly regions where they diverge, such as where the Radon method estimates very fast speeds (nearly horizontal lines) presumably caused by the birth or death of wavepackets, where the Radon method obtains the incorrect sign (e.g., 10°E, 30 September 2009 in Fig. B1b and several other places), or where our method estimates very slow speeds (e.g., Fig. B1a, upper right) possibly caused by mechanical or buoyancy forcing or, more likely, nonzonal propagation of PV anomalies.
We stress that one should not expect the two methods to agree perfectly. First, note that our method is applied at each time step independently: no information from other times is used to estimate the phase speed, though we do use information from nearby latitudes. In contrast, the Radon transform uses information from nearby times, but no extra latitudinal information. Second and more fundamentally, selecting a dominant wave from a rich wavefield is not only nontrivial—there is no unique, objective answer. Indeed, the 2D wavelet transform (Chen and Chu 2017) yields, at each point, the power as a function of period and wavelength, but there are many ways to reduce this to a single measure of a dominant wave.B1 This is also true of our method and the Radon transform. For the Radon transform, we could use any norm, rather than the absolute value of the horizontal integral. Rather than selecting the angle that maximizes a measure, we could choose its expected value. Rather than applying this technique to the raw data, we could filter the data (Polito and Liu 2003), preselecting certain spectral bands. If these algorithmic choices are well chosen, the Radon technique produces reasonable estimates of the phase speed of a wave component carrying high energy density. Reasonable estimates are all that can be asked for; there is no perfect answer (Subrahmanyam et al. 2009). Figure B1 shows that our reference-frame shifting method produces reasonable estimates of the eddy phase speed.
Vertical Mode Rossby Wave Phase Speeds
a. Vertical modes in depth coordinates
Traditionally, one calculates Rossby wave phase speeds by linearizing the equations of motion about a background state, and seeking wave solutions with arbitrary vertical structure. Then, one transforms that structure into an infinite sequence of mutually orthogonal vertical modes determined by the background state and boundary conditions that dictate vanishing vertical velocity at the upper and lower boundaries. This procedure resembles the following buoyancy-coordinate derivation, so is omitted for brevity; see Killworth et al. (1997) and Killworth and Blundell (2003a) for details.
The bottom boundary condition is debated, however. An alternative, appropriate for steep bathymetry, is vanishing horizontal velocity (Tailleux and McWilliams 2001). This produces vertical modes that better represent the oceanic flow (de La Lama et al. 2016), although a bottom trapped mode must be included to form a complete basis.
b. Vertical modes in buoyancy coordinates
We now derive the Rossby wave vertical modal structure and phase speed in buoyancy coordinates. In the ACC jets, can be as large as β and hence is important for small-scale features, but the PV gradient’s large-scale features are mostly determined by variations of horizontal stratification, not relative vorticity. Thus, we consider the buoyancy-conserving planetary geostrophic (PG) equations, which are
akin to (9). Happily, the Montgomery potential m serves as a state variable, determining all other fields.
Linearizing about a background state M that is steady () and zonal () yields
where m = M + m′. For wave solutions, let to find
With the background large-scale PV, and the background zonal velocity, the zonal Rossby wave phase speed c in a continuously stratified 3D fluid is
This has the form of (13): c is Doppler shifted by the local background flow, and the intrinsic phase speed, c − U, is proportional to the PV gradient, keeping high PV to its right (westward for ) when , the configuration that is consistent with a wave-like vertical structure. In general, a solution can have sections of the fluid column with , but these regions would grow exponentially (in ), so are quickly countermanded.
At the upper boundary, (C1c) implies at the surface buoyancy. Assuming the background field also satisfies these constraints, then at , hence
Similarly, at b = bb the bottom buoyancy, again leading to at , so
which is equivalent to vanishing vertical velocity at the bottom. The alternative boundary condition, zero horizontal flow at the bottom, is
To better capture the background state, the PV gradient (which contains much high-frequency noise) can be taken from the functional relationship between Q−1 and M,
where s = dQ−1/dM is depth dependent. This can be rearranged as
so “Ma”−1 arises naturally.
c. Comparison of methods
Figures C1a–d map the Rossby wave phase speed of the most unstable mode. The results, both where such a mode exists and its phase speed where it does, are sensitive to the choice of vertical coordinate and bottom boundary condition. The buoyancy is using the carefully defined orthobaric density ρν, so we should not hope that a different buoyancy would align the results. Figures C1e and C1f estimate the eddy phase speed as in section 4c on σν = 27.1 and σν = 27.3, respectively. As these estimates use information from a single surface, they should not be expected to match the phase speeds of Figs. C1a–d that are influenced by the entire water column. We selected the most unstable mode because, as the fastest growing, its amplitude will generally be large, and hence it is a good candidate for moving PV anomalies around on orthobaric surfaces, making contact with our reference-frame shifting method. Though this connection might apply somewhat broadly, it is not guaranteed. For example, if this mode has a node at the depth of our orthobaric surface of interest, then its relevance to the movement of eddies on that surface is, at best, indirect. Nonetheless, the eddy phase speed estimated by these new methods agrees reasonably well with the Rossby wave vertical mode results.
More quantitatively, Table C1 provides Pearson’s correlation coefficient (taking values in [−1, 1]) between the phase speeds in each pair of panels in Fig. C1.C1 This coefficient is positive for all pairs of panels. Moreover, it is roughly the same between pairs in Figs. C1a–d as it is between one of Figs. C1a–d and Fig. C1e, indicating that there is about as much ambiguity in the classic eigenvalue method of evaluating the Rossby wave phase speed as there is when the frame-shift method is applied on the σν = 27.1 surface. Evidently the predicted phase speeds on the σν = 27.3 surface correlates less well, though still positively, with the vertical modes’ phase speeds, suggesting this surface is somewhat deeper than the depth at which the most unstable mode exhibits its maximum amplitude. All told, these results lend support to the methods of this paper.
Idealized Tests for the Frame-Shift Method
To help illuminate uncertainties in various estimates of “Ma,” and as a proof of concept, we test our algorithms on idealized experiments with a 1.5-layer shallow-water model, having an active upper layer overlying a constant abyssal zonal jet (cf. Dowling and Ingersoll 1989). The model equations (solved in spherical coordinates, but presented in Cartesian coordinates for simplicity) are
Here, (u, υ) is the upper-layer velocity, f the Coriolis parameter, ζ = ∂xυ − ∂yu the relative vorticity, ϕ the upper-layer geopotential, K = (u2 + υ2)/2 the specific kinetic energy, Fx and Fy the viscous forcing, h = (ϕ − ϕ2)/g′ the upper-layer thickness with g′ the reduced gravity and ϕ2 the lower-layer geopotential.
The domain is periodic over 90° in longitude, with free-slip walls at 70° and 30°S. The planetary radius is R = 6378 km, and the angular velocity is Ω = 7.272 × 10−5 s−1. The grid resolution is 512 × 256, yielding 0.176° zonally (12.6 km at 50°S) by 0.156° meridionally (17.3 km). Time is discretized using a third-order Adams–Bashforth scheme. The spatial discretization is a centered second-order scheme for the continuity equation, and the Sadourny (1975) enstrophy conserving scheme for the momentum equations. Viscosity and hyperviscosity use uniform coefficients of 10 m2 s−1 and 10 m4 s−1, respectively (matching SOSE). This is a simple, Jupiter-like configuration that allows control of “Ma” (Stamp and Dowling 1993) by carefully choosing the abyssal geopotential ϕ2 as follows.
First, the quasigeostrophic approximation reduces (D1) to material conservation of the QG PV, . Here, ψ = ϕ/f0 and ψ2 = ϕ2/f0 are the geostrophic streamfunctions for the upper and lower layers, respectively, and Ld is the deformation radius satisfying , where f0 is f at a central latitude and H is a spatial average of the upper-layer thickness h. Material conservation of QG PV is by the geostrophic flow: qt − ψyqx + ψxqy = 0. Linearizing this discovers zonal Rossby waves with intrinsic speed in the long wave limit. Using u = −ψy, the QG estimate for “Ma” is
We prescribe the initial upper-layer velocity as a zonal jet, with zonal velocity , having amplitude u0 = 0.25 m s−1, width , and center θ0 = 50°S. Next, the QG PV gradient (Ingersoll and Cuong 1981) is rearranged for the deep flow u2, then (D2) yields . This u2 profile specifies the lower-layer geopotential by gradient balance: , where is the lower-layer relative vorticity. Finally, the initial upper-layer geopotential is gradient balanced with the upper-layer jet: . The constant of integration for ϕ2 is irrelevant, but for ϕ it is determined by requiring g′H = 10 m2 s−2 in all experiments, corresponding to at 50°S. Because the model uses the primitive, not QG, equations, this initial condition yields only an approximately constant “Ma” basic state, but is sufficient to test our algorithms.
To break symmetry and excite waves, we add eight small, geostrophically balanced (cyclonic and anticyclonic) vortices by adding elliptical perturbations to the initial geopotential. Their flow speed varies, the fastest reaching 0.1 m s−1. Their shape also varies, with semimajor and semiminor axes ranging from 3° × 2° to 5° × 3°. They are scattered north and south of the jet, with centers not closer than 2° of 50°S. The results do not depend significantly on these characteristics, but all experiments have the same initial vortices, and geopotential.
We run five experiments, varying the initialization “Ma” as 0.8, 0.9, 1.0, 1.1, and 1.2. Each instance is run for 8 years, averaging the output state variables, the Bernoulli potential (area-weighted then averaged from the thickness grid to the vorticity grid), and the PV over the last 5 days.
From this output, we estimate “Ma” by three methods. The first is the QG estimate (D2), with Ψy/qy the slope of an orthogonal regression of to q. Second, the linearized PG equations yield the total and intrinsic Rossby wave speeds, and where q = f/(g′h) is the PG PV. Thus, the PG estimate for “Ma” is [cf. Dowling 1993, Eq. (12)]
Figure D1 shows these “Ma” estimates, using data from the north flank of the jet (44°–47°S) and, separately, from the south flank (52°–55°S). We exclude the jet’s core, where Q, B, and A are nearly zonally uniform, so multicollinearity degrades the empirical fits. The α estimate for “Ma” agrees fairly well with the QG and PG estimates, and captures the changes to “Ma” on the jet flanks as the initialization “Ma” is changed.
Figure D2 maps “Ma” estimated at each grid point using local data, as in section 4c, at the end of the “Ma” = 1.1 experiment. The QG results are very similar to the PG results and hence omitted. The PG method estimates “Ma” very close to 1.1 in the core of the jet. Outside the jet core, the methods give qualitatively similar “Ma” estimates, despite complex spatial structure. Note the large anticyclone centered at (2°W, 56°S), taking “Ma” ≈ 1.
The perfectly zonal background state of this idealized jet poses the hardest of problems for the α method. For nonzonal flows, such as the ACC or its representation in SOSE, we expect the α method to work better still.
Denotes content that is immediately available upon publication as open access.
Technically, Arnol’d’s theorem also requires that dψ/dq does not asymptote to 0 or .
Section 3.2.4 of Stanley (2018) discusses a technical concern with the regression method of these studies.
These criteria are for the instantaneous fields Q, B, and A. In the third case, Q and B are evolving in time yet bound together by their functional relationship, which eliminates the possibility of traveling Rossby waves. Also note, the predilection with zonality is because the rotation axis in the α-shifted frame is the same as Earth’s.
Theoretically, regions of intrinsic eastward Rossby wave/eddy speed are either due to a reversal of the meridional PV gradient, or a sign reversal of the coefficient multiplying the PV gradient, e.g., in (C5). Another cause could be numerical: is calculated as c − u, with u taken from the local grid point but c estimated from surrounding data points.
We tested this method also, but found it too sensitive to parametric choices.
We use Pearson’s correlation coefficient here, rather than the Nash–Sutcliffe efficiency, because there is no “truth” on which to base the latter.