## Abstract

The phase speed spectrum of ocean mesoscale eddies is fundamental to understanding turbulent baroclinic flows. Since eddy phase propagation has been shown to modulate eddy fluxes, an understanding of eddy phase speeds is also of practical importance for the development of improved eddy parameterizations for coarse resolution ocean models. However, it is not totally clear whether and how linear Rossby wave theory can be used to explain the phase speed spectra in various weakly turbulent flow regimes. Using linear analysis, theoretical constraints are identified that control the eddy phase speed in a two-layer quasigeostrophic (QG) model. These constraints are then verified in a series of nonlinear two-layer QG simulations, spanning a range of parameters with potential relevance to the ocean. In the two-layer QG model, the strength of the inverse cascade exerts an important control on the eddy phase speed. If the inverse cascade is weak, the phase speed spectrum is reasonably well approximated by the phase speed of the linearly most unstable mode. A significant inverse cascade instead leads to barotropization, which in turn leads to mean phase speeds closer to those of barotropic-mode Rossby waves. The two-layer QG results are qualitatively consistent with the observed eddy phase speed spectra in the Antarctic Circumpolar Current and may also shed light on the interpretation of phase speed spectra observed in other regions.

## 1. Introduction

Mesoscale eddies (on scales of tens to hundreds of kilometers) are ubiquitous in the ocean, and are believed to be crucial in the transport of tracers and the rectification of the mean flow (e.g., Gill et al. 1974; Johnson and Bryden 1989; Hallberg and Gnanadesikan 2006; McWilliams 2008; Waterman et al. 2011). Yet many of their fundamental properties are still poorly understood. This contribution focuses on the question of what controls the zonal propagation (i.e., phase speed) of mesoscale eddies in a two-layer model of quasigeostrophic baroclinic turbulence. Eddy phase speed has recently been highlighted as an important factor in modulating the magnitude of eddy fluxes (e.g., Marshall et al. 2006; Ferrari and Nikurashin 2010; Abernathey et al. 2010; Klocker et al. 2012a; Klocker and Abernathey 2014; Bates et al. 2014). As described by the theory of Ferrari and Nikurashin (2010), eddy propagation relative to the background mean flow suppresses the efficiency of eddy-driven mixing (see also Klocker et al. 2012a,b). Incorporating this effect into mesoscale parameterizations may improve coarse-resolution ocean models that do not resolve mesoscale fluxes (Bates et al. 2014). A deeper physical understanding of what controls eddy phase speeds in simple models such as the one studied here is an important step toward this goal.

Observationally, remote sensing of sea surface height (SSH), temperature, and color has provided the primary source of information about global mesoscale eddy characteristics. Several different methods have been used to characterize eddy propagation speeds in these datasets. Arguably the most straightforward method is to fit lines to a longitude–time Hovmöller diagram, either by eye or using a Radon transform, as first implemented by Chelton and Schlax (1996) on the 3 years of TOPEX/Poseidon SSH data available at the time. Wavenumber–frequency spectral analysis has also been used to characterize phase speeds of SSH alone (Zang and Wunsch 1999; Wunsch 2009; Wortham 2013) and of the covariance between SSH, SST, and ocean color (Hill et al. 2000; Cipollini et al. 1997, 2001; Killworth and Blundell 2004; Abernathey and Wortham 2015). Fu (2009) used space–time lag correlation to identify the eddy propagation patterns and speed. Finally, the direct tracking of individual coherent eddy features has also been implemented (Chelton et al. 2011). While some of the studies cited above focus on coherent vortices, we here define the eddy field to include any transient perturbation to the time-mean flow, consistent with the use in the atmospheric science and general turbulence literature (well described in textbooks; e.g., Pedlosky 1979; Holton 1992; McWilliams 2006). The eddy field can then be characterized by the frequency–wavenumber spectrum, and eddy phase speeds are unambiguously defined by the ratio of frequency to wavenumber.

A common feature of mesoscale turbulence is that, over most of the ocean, the observed phase propagation is westward, except for the Antarctic Circumpolar Current (ACC) region and in western boundary current regions where Doppler shifting by the background mean flow is significant. Furthermore, frequency–wavenumber analysis reveals that, at least in some regions, energy is organized along “nondispersive” lines (Wunsch 2009; Early et al. 2011) (i.e., all the energy appears to propagate with the same phase speed, regardless of wavenumber), as is necessary to maintain coherent vortices (Taylor 1938; McWilliams and Flierl 1979; Nof 1981). For coherent and nondispersive eddies, the eddy phase speed is identical to the group velocity, so theories of coherent vortices may be invoked. However, it is less clear how to interpret phase speeds for dispersive dynamics in the presence of instability, which is in fact not rare in the ocean.

Mesoscale eddies are inherently nonlinear, and it is still a matter of debate to what extent their properties can be understood by linear dynamics (Chelton et al. 2007; Wunsch 2009). Nevertheless, linear Rossby wave theory (well described in textbooks; e.g., Pedlosky 1979; Vallis 2006) is an important starting point for understanding mesoscale eddy properties, especially in regions where the flow is only weakly turbulent (Tulloch et al. 2009). It is appealing to attempt to use linear dispersion relations to describe the nonlinear mesoscale eddy phase speeds. However, as Chelton and Schlax (1996) first noted, the standard theory for freely propagating linear baroclinic Rossby waves often underestimates the observed phase speeds. Moreover, the commonly observed nondispersive frequency–wavenumber spectra are inconsistent with the linear Rossby wave dispersion relation. The reasons for the observed deviations from the linear Rossby wave theory became a matter of intense discussion (Killworth et al. 1997; Dewar 1998; Killworth and Blundell 2003, 2005, 2007; Tulloch et al. 2009). Berloff and Kamenkovich (2013a,b) analyzed the spectral properties of fully nonlinear solutions in a hierarchy of idealized flows and demonstrated that in fact a very substantial part of the flow properties are explained by the underlying linear dynamics if meridional and temporal variations in the background flow are taken into account. Klocker and Marshall (2014) recently argued that empirically eddy phase speeds over most of the ocean are reasonably well approximated by the long-wave limit of the first baroclinic mode Rossby wave phase speed, Doppler shifted by the depth averaged mean velocity . The Doppler shifting is of particular importance in the ACC region, where it explains the transition from westward to eastward phase propagation—the advection by the strong eastward mean flow here exceeds the flow-relative westward propagation speed of baroclinic Rossby waves.

This paper aims to improve our understanding of eddy phase propagation by focusing on a model of quasigeostrophic (QG) baroclinic turbulence. The characteristics of eddy phase propagation are analyzed in a fully nonlinear model and compared to linear theory. Using the insights from QG theory, analysis of SSH observations is also presented as a comparison.

One motivation for our work is the paper by Early et al. (2011), which demonstrated the important role of nonlinearity in the formation and propagation of coherent mesoscale eddies. That study examined the evolution of a reduced-gravity QG model with and without nonlinearity. When nonlinearity was present, coherent vortices formed and the wavenumber–frequency spectrum collapsed onto a nondispersive line, with the phase speed given approximately by the long-wave limit of the reduced-gravity mode’s dispersion relation. The model used by Early et al. (2011), however, did not include baroclinicity in the background state, and cannot simulate the generation of eddies from baroclinic instability. Instead eddies were seeded using quasi-random initial conditions.

In a baroclinic flow, eddies and turbulence can arise as a result of baroclinic instability, and the linear dispersion relation becomes more complex (literally). The relevance of neutral mode baroclinic Rossby waves (which represent true solutions to the linearized equations only in the absence of baroclinicity) becomes unclear in the presence of a baroclinically unstable shear flow. Here we examine the eddy phase propagation properties of a homogeneous two-layer QG model, which is arguably the simplest possible model that can generate eddies through baroclinic instability. This model may be expected to reproduce some characteristics of eddies in the ACC region and possibly the western boundary current extensions, where eddies are generated by deep-mode “Phillips-like” baroclinic instability (Tulloch et al. 2011). The model instead is less likely to allow for an adequate representation of the eddying flow in the subtropical regions, where we do not expect such deep-mode instability (Tulloch et al. 2011).

This paper is structured as follows. Section 2 starts with a description of the two-layer QG model and the relevant parameters. A linear stability analysis is performed to obtain the phase speeds of the unstable modes, and the results are compared to the phase speeds of neutral modes obtained in certain limit cases. In section 3, we conduct fully nonlinear simulations to explore the dependence of the eddy phase speeds and frequency–wavenumber spectra on various model parameters; the results are interpreted in terms of the linear theory discussed in section 2. We consider some general properties of baroclinic turbulence, such as the barotropization in the turbulent cascade. With this, we argue that, in the two-layer QG model, much of the dependence of eddy phase speeds on external parameters can be understood in terms of the linear theory. In section 4 we compare the results of the QG model to SSH observations. Conclusions are given in section 5.

## 2. Linear analysis

This section describes the phase speed predictions that can be made based on linear analysis of a two-layer QG model.

### a. The model

We use a two-layer QG model on a *β* plane with a flat bottom following the formulation in Flierl (1978). To place the model in an oceanic regime, we consider the two layers to have a tunable layer thickness ratio *δ* = *H*_{1}/*H*_{2}, where *H*_{1} and *H*_{2} are the layer thicknesses for the upper layer and lower layer, respectively; for the ocean, *δ* < 1. The two-layer model is forced by an imposed background vertical flow shear of Δ*U* = *U*_{1} − *U*_{2}, where *U*_{i} is the background zonal velocity in layer *i*, with *i* = 1 denoting the upper layer and *i* = 2 denoting the lower layer. Also, *Q*_{i} is the background potential vorticity (PV) and *q*_{i} denotes the perturbation PV.

The governing equations for the PV perturbations are

where the last term on the left-hand side of the second equation denotes a linear bottom friction, with a decay time scale *τ*_{f}. The Jacobian represents nonlinear wave–wave interaction and ssd is small-scale dissipation, which will be ignored in the linear analysis.

The perturbation streamfunctions *ψ*_{1} and *ψ*_{2} are related to perturbation PV through the inversion relation

Similarly, the background PV gradient can be related to the vertical shear and planetary vorticity gradient, *β*, via

where *F*_{1} and *F*_{2} are defined as

where *L*_{d} is the baroclinic deformation radius. The barotropic component of the streamfunction (*ψ*_{BT}) and the baroclinic component of the streamfunction (*ψ*_{BC}) can be defined (Flierl 1978; Arbic and Flierl 2004) as

In this study, we assume a positive vertical mean flow shear (i.e., eastward), which is consistent with observations over most eddy-rich regions of the world oceans.

For instability to occur, in the absence of friction, the PV gradient must change sign between the upper layer and lower layer, so as to support counterpropagating Rossby waves. Note that *Q*_{y1} is always larger than zero for a positive vertical mean flow shear, and thus stability is governed by the nondimensional criticality parameter:

If *ξ* > 1 (equivalent to *Q*_{y2} < 0), instability can occur in the inviscid limit.

More generally, in the presence of friction, we can define three independent nondimensional parameters. Normalizing length scales with *L*_{d} and time scales with *L*_{d}/Δ*U*, we obtain the following:

the layer thickness ratio

*δ*,the nondimensional bottom frictional damping rate , and

the nondimensional planetary vorticity gradient .

The nondimensional planetary vorticity gradient is related to the criticality parameter via

These three parameters collectively determine the basic properties of the two-layer QG model. We explore different flow regimes by varying these three parameters.

### b. Linear dispersion relations

In a square, doubly periodic domain, we assume a plane–wave solution:

where *k* and *l* are the zonal and meridional wavenumbers, and denotes the real part. Generally *c* ≡ *ω*/*k* is a complex number whose real part represents the zonal phase propagation speed and imaginary part represents the growth/decay rate of an unstable mode. Linear solutions are obtained by replacing terms in the linearized Eq. (1) with this plane–wave solution.

Before considering the full linear instability problem, we first review the three neutral mode solutions that arise in relevant limit cases. When the lower layer is either much deeper than the upper layer, or when bottom friction is very strong, the lower-layer flow becomes very weak (i.e., *ψ*_{2} ≪ *ψ*_{1}). The lower-layer equation in Eq. (1) then becomes a second-order equation, while the linearized upper-layer equation to first order yields the reduced gravity model (considered by Early et al. 2011):

The dispersion relation for the reduced gravity mode is

where the subscript RG indicates reduced gravity.

On the other hand, in the limit of vanishing vertical shear (Δ*U* → 0) and neglecting bottom friction, the eigenvalue problem posed by the linearized version of Eq. (1) yields two neutral modes: a barotropic mode and a baroclinic mode. The dispersion relation for the barotropic mode is

where *U*_{b} ≡ (*δU*_{1} + *U*_{2})/(1 + *δ*) is the barotropic mean flow.^{1}

The dispersion relation for the baroclinic mode is

where *L*_{d} is the deformation radius of the first baroclinic mode. The long-wave limit of the baroclinic mode has recently been argued to provide a good fit to the eddy phase propagation characteristics obtained from altimetric observations. (Tulloch et al. 2009; Klocker and Abernathey 2014; Klocker and Marshall 2014; Abernathey and Wortham 2015).

Comparing *c*_{BC} with *c*_{RG}, three major differences appear: 1) *c*_{BC} is Doppler shifted by the depth-average flow, while *c*_{RG} is Doppler shifted by the upper-layer flow; 2) *c*_{BC} feels the planetary vorticity gradient *β*, while *c*_{RG} feels the upper-layer PV gradient *Q*_{y1}; and 3) *c*_{BC} feels the first baroclinic deformation radius, while *c*_{RG} feels the “upper-layer deformation radius,” . In the long-wave limit, the reduced gravity mode dispersion relationship simplifies due to a cancellation between the effects of the upper layer flow on the Doppler shift and on the PV gradient, which yields *c*_{RG} = −*β*/*F*_{1} + *U*_{2}. In this limit the reduced gravity dispersion relationship resembles that for the baroclinic mode, but with the difference that the apparent Doppler shift is given by the lower-layer flow speed, and the relevant deformation radius is that of the upper layer.

### c. Linear instability

This subsection revisits the instability properties of the full two-layer QG model. Linearizing Eq. (1) and using Eqs. (2) and (8) yields

For nontrivial solutions, the determinant of coefficients must be zero. This provides a quadratic equation for *c*, which yields two solutions. In the absence of bottom friction, a pair of complex conjugate solutions is obtained for baroclinically unstable mean shears (*ξ* > 0), denoting one growing and one decaying mode.

We now review the impact of nonzero friction on the linear stability analysis, which affects both growth rates and phase speeds. Numerical solutions for the growth rate and phase speeds are computed using parameters roughly characteristic of the ACC: *δ* = 0.25, Δ*U* = 0.04 m s^{−1} and *L*_{d} = 15 km. The layer depth ratio of 0.25 is roughly consistent with the depth of the sign reversal in the observed extratropical zonal-mean QGPV gradient (Tulloch et al. 2011). Similarly, the vertical shear and deformation radius are in rough agreement with observations in the ACC region. We will return to a more quantitative comparison between the QG model and the observations in section 4. The meridional wavenumber *l* is set to zero to obtain the fastest growing modes (e.g., Vallis 2006). The role of nonzero meridional wavenumber *l* will be discussed in the following sections.

The upper panel of Fig. 1, shows the linear growth rates as a function of wavenumber for different values of bottom friction. When friction is absent, both a short-wave cutoff and long-wave cutoff can be identified. However, with increasing frictional strength, the short waves and some long waves are destabilized while maximum growth rates decrease, as first discussed by Holopainen (1961). The length scale of the most unstable modes also slightly decreases with increasing frictional strength. The lower panel in Fig. 1 shows the corresponding phase speeds. Unstable modes propagate eastward [ > 0] at most wavelengths. Only when friction is weak, the phase speed can be slightly westward [ < 0].

The dispersion curves for the reduced gravity mode (dashed) and the barotropic mode (dot-dashed) are also plotted in the lower panel of Fig. 1. The phase speeds of the linearly unstable modes consistently fall between these two neutral curves. Since *c*_{RG} represents the limit case where the lower layer is at rest, it is perhaps not surprising that the phase speeds of the unstable modes approach *c*_{RG} as bottom friction is increased. In fact, it appears that (10) is a useful approximation to the phase speed of the unstable modes in a realistic parameter regime.

The baroclinic dispersion curve *c*_{BC} is also plotted in the lower panel of Fig. 1 for comparison. The baroclinic mode phase speed shows a relatively weak wavenumber dependence and is less clearly related to the full linear instability problem, although the overall phase speed is close to that of the fastest growing modes. Whether the agreement between the baroclinic mode phase speed and the phase speed of the fastest growing mode is a coincidence or is driven by some physical mechanism remains unclear.

Notice that *c*_{BC} and *c*_{RG} are not always easy to distinguish in the real ocean. For example, in the subtropical gyres, the first baroclinic mode typically becomes strongly surface intensified, corresponding to an upper layer that is much shallower than the lower layer. The upper-layer deformation radius in *c*_{RG} then is well approximated by the first baroclinic deformation radius as in *c*_{BC}. Moreover, the vertical shear is typically weak, such that the upper-ocean PV gradient in *c*_{RG} approximately reduces to the planetary vorticity gradient as in *c*_{BC}, and Doppler shifting by the barotropic flow becomes indistinguishable from Doppler shifting by the upper-ocean flow. On the other hand, in other regions, such as the ACC, the differences between *c*_{BC} and *c*_{RG} should be more pronounced.

## 3. Nonlinear model

The linear analysis described above reveals the phase speeds of linear modes. However, in a nonlinear equilibrated state, it is not a priori clear which wavenumbers get energized and whether the eddies in the nonlinear flow actually follow a linear dispersion relation. This section presents frequency–wavenumber spectra and phase speed spectra from a series of nonlinear two-layer QG simulations and compares them to the predictions from linear theory.

### a. Model description and methods

We numerically solve the fully nonlinear two-layer QG model described by Eq. (1) using the open-source Python package pyqg (Abernathey et al. 2015). We use a doubly periodic domain with a horizontal resolution of 256 × 256 grid points in real space. Enstrophy is removed near the grid scale with an exponential filter that is identical to that described in Arbic and Flierl (2004).

The parameters in the control simulation are as in the linear analysis discussed in section 2:

The frictional decay rate of (20 days)^{−1} is empirically chosen to obtain eddy properties that are roughly consistent with observations in the ACC region. Determination of a realistic frictional time scale from first principles is not straightforward, as it crudely represents various routes to dissipation near the bottom boundary [see also Arbic and Flierl (2004)]. The positive mean flow vertical shear of 0.04 m s^{−1} is consistent with observations in the ACC region. The parameters in our two-layer QG model correspond to a nondimensional criticality parameter of *ξ* ≈ 2.37, which is moderately supercritical.

As case studies, we present three reference simulations. These three simulations use the above listed parameters but with different frictional strengths. In addition to the control run which uses *τ*_{f} = 20 days, the strong friction run and the weak friction run use *τ*_{f} = 10 days and *τ*_{f} = 40 days, respectively. In Fig. 2, the Hovmöller diagram of the upper-layer streamfunction in the control run suggests no clear preferential direction of propagation. Comparing all three simulations in Fig. 2 suggests that decreasing the frictional strength to (40 days)^{−1} enhances westward propagation; while increasing the frictional strength to (10 days)^{−1} favors eastward propagation.

To comprehensively examine the phase speeds in the two-layer QG model, in addition to the three reference simulations, we conduct several groups of experiments, each varying only one nondimensional parameter. Each group includes more than 20 simulations to explore the effect of variations in the nondimensional parameters over at least one order of magnitude.

In the experiments dubbed THIC, the layer thickness ratio

*δ*is varied, while and are held fixed.In the FRIC experiments, is varied, while

*δ*and are held fixed.In the BETA experiments, is varied, while

*δ*and are held fixed. The additional group BETA-hf is similar to the BETA experiments, but with strong frictional strength.

The parameters considered here do not give rise to regimes with strong jets (Rhines 1975; Maltrud and Vallis 1991). Among all simulations, the output data of the first 10 years are dropped, and streamfunction snapshots for both layers from the following 50 years are sampled every 5 days.

Analysis of the nonlinear simulations is based on spectral analysis in wavenumber and frequency space (Hayashi 1971). The streamfunction *ψ*_{i}(*x*, *y*, *t*) is Fourier-transformed into spectral space as , where *k*, *l*, and *ω* are zonal wavenumber, meridional wavenumber, and frequency, respectively. By multiplication of the three-dimensional streamfunction variance spectrum with the total wavenumber and integration over all frequencies *ω*, the total eddy kinetic energy (EKE) as a function of zonal and meridional wavenumber can be obtained:

Figure 3 shows the upper-layer total EKE as a function of zonal and meridional wavenumber, *E*_{1}(*k*, *l*), for our three reference simulations. The total EKE is nearly isotropic for all three reference simulations. Because of the isotropic behavior in our model, and the focus on zonal phase propagation, we emphasize the dependence of EKE on frequency and zonal wavenumber *k*. Rather than showing total EKE, we moreover focus on meridional velocity variance spectra (hereinafter “transverse spectra”), which can be obtained by integrating the streamfunction variance spectrum over all zonal wavenumbers and multiplying by *k*^{2}:

The advantages of using transverse energy spectra are discussed in the appendix.

### b. General results

In this section, we discuss the general characteristics of the turbulent flow fields in the fully nonlinear two-layer QG simulations, and explore the dependence of the eddy properties on the three nondimensional parameters. Linear theory highlights the different propagation behavior of the barotropic and baroclinic modes, so this decomposition is a useful starting point for our analysis of the nonlinear simulations. Using Eq. (5), we can separate the energy spectra into barotropic and baroclinic components. In Fig. 4, solid curves show the ratio of barotropic mode EKE to total EKE as a function of zonal wavenumber. This ratio is bigger toward small wavenumbers in the control experiment, indicating a trend toward enhanced barotropization at larger scales. Comparing all three reference simulations, the barotropization is enhanced when friction is weak, and vice versa. This barotropization can partially be attributed to an increasing barotropic contribution to the linearly most unstable modes. However, linear theory alone does not fully predict the strong barotropization observed at low wavenumbers. The strong barotropization is instead likely attributable to nonlinearity.

In Fig. 5, we show the parameter dependence of the ratio of barotropic EKE to total EKE by averaging over all horizontal wavenumbers and frequencies, in all simulations in THIC, FRIC, BETA, and BETA-hf. The flow becomes more barotropic when 1) the layer thickness ratio approaches unity, 2) friction is reduced, and 3) is reduced. All of these dependencies are in qualitative agreement with previous results (e.g., Arbic and Flierl 2004). Throughout most of the parameter range, ratios of barotropic EKE and total EKE in nonlinear simulations are larger than the corresponding ratios predicted by linear instability analysis, suggesting an important role of nonlinearity. In the simulations whose total EKE is dominated by the barotropic mode, we expect the dispersion relation of the barotropic mode to provide a better prediction of the eddy phase speeds. We will test this hypothesis in the next subsection.

In addition to barotropization, nonlinear eddy–eddy interactions tend to lead to an increase in the horizontal eddy scale, beyond the scale of the most unstable mode. To characterize length scales, we integrate the power spectra over all frequencies to obtain a univariate zonal wavenumber spectrum for each experiment. Furthermore, the spectral energy budget at the zonal wavenumber *k* can be written as below, following the formalisms in Jansen and Held (2014):

Here, *E*_{k}, *S*_{k}, *TA*_{k}, *TK*_{k}, *F*_{k}, and *V*_{k} represent total kinetic energy for both layers, the extraction of available potential energy (APE) from the prescribed background state, spectral transfer of APE, spectral transfer of kinetic energy, frictional dissipation by bottom drag, and small-scale dissipation, respectively.

The sum of *S*_{k} and *TA*_{k} can be interpreted as the conversion from available potential energy to eddy kinetic energy. The values of *S*_{k} and *TA*_{k} can be calculated through

where the asterisk denotes the complex conjugate, and denotes the Fourier transform, here with respect to zonal wavenumber.

In Fig. 6, we use the diagnosed spectra of net conversion from APE to EKE (top row) to indicate the scale at which kinetic energy is injected into the system. The peak in the conversion spectrum roughly coincides with the peak in the linear growth rate, which is shown in the bottom row. The departure between the peak in the conversion spectrum and the peak in the barotropic EKE spectrum gives an indication of the inverse cascade range.

The THIC experiments (first column) show that simulations with equal layer depths exhibit a slightly enhanced inverse cascade compared with simulations with a shallower upper layer. The FRIC experiments (second column) exhibit the most significant changes in the inverse cascade range, with a strong inverse cascade found in limit of weak friction, consistent with Arbic and Flierl (2004). There is no direct correspondence between the most unstable wavenumbers and the energy containing wavenumbers, and linear stability analysis only predicts the injection scale. While the large inverse cascade in the limit of weak friction is the result of strong nonlinear eddy–eddy interactions, we note that the largest eddies are near the wave-turbulence crossover (Rhines 1979), suggesting that linear dynamics are likely to remain relevant. The BETA and BETA-hf experiments show that the energy cascade range is reduced for large , which is in agreement with the general understanding that the beta effect limits the inverse energy cascade (Rhines 1979). However, for small the halting scale converges to a constant value, which again depends on the bottom friction (cf. the third and fourth columns in Fig. 6).

In summary, although all nondimensional parameters influence the strength of the inverse cascade, it is here found to be most sensitive to the frictional strength. Many theoretical arguments have been developed to predict the halting scale of the inverse energy cascade (e.g., Rhines 1979; Held and Larichev 1996; Arbic and Flierl 2004; Thompson and Young 2007; Jansen et al. 2015), and the results obtained here appear to be broadly consistent with this previous work. However, our focus is on the eddy phase speed, and a detailed comparison between inverse cascade strength and different scaling laws is beyond the scope of this contribution.

### c. Phase speed spectra

#### 1) Case studies

To characterize zonal phase speeds in the nonlinear simulations, we now consider the wavenumber–frequency spectrum, which can be compared to the linear dispersion relationships discussed in section 2. To further quantify the dominant phase speed, the wavenumber–frequency spectrum can be transformed to a wavenumber–phase speed spectrum, following Randel and Held (1991) (see also Abernathey and Wortham 2015). For easy comparison to satellite observations, which are available only for the near-surface flow, all analysis of wavenumber–frequency spectra is presented for the upper model layer. The top row in Fig. 7 shows the wavenumber–frequency transverse EKE () spectra for our three reference cases: the control run, weak friction run, and strong friction run. Since is proportional to *k*^{2} times the streamfunction variance [Eq. (16)], the spectra by construction vanish in the limit of vanishing zonal wavenumber and highlight the wavenumber range of the energy-containing eddies. To make a connection between spectra and the streamfunction spectra often considered in the oceanographic literature (e.g., Wunsch 2009), and to investigate the dispersion relation at all zonal wavenumbers, the middle row of Fig. 7 shows the wavenumber–frequency spectra normalized by the total power at each zonal wavenumber. Finally, the bottom row shows the interpolated wavenumber–phase speed spectra of (without any normalization).

To facilitate the comparison with linear theory, the linear dispersion relations discussed in section 2 are also plotted in Fig. 7. We adopt spectral moments to characterize the length scales of the turbulence, as in numerous prior studies (e.g., Rhines 1975; Scott and Wang 2005; Tulloch et al. 2011; Abernathey and Wortham 2015). To quantitatively examine the accuracy of linear predictions, at each zonal wavenumber, we compare the diagnosed mean eddy phase speeds, , with the predicted mean eddy phase speeds according to the different linear dispersion relationships, . The diagnosed mean eddy phase speed is defined as the first moment of the wavenumber–phase speed spectrum :

The predicted mean eddy phase speed for a given mode is computed as

where *C*_{X}(*k*, *l*) is given by the linear dispersion relation of the respective mode, that is, *C*_{BT}(*k*, *l*), *C*_{RG}(*k*, *l*), or *C*_{BC}(*k*, *l*). Equation (20) can be used to infer an effective meridional wavenumber, *l*_{eff}, such that , which mostly lies in between *l* = 0 and *l* = *k* but varies significantly with *k* (Fig. 7). The varying effective meridional wavenumber is not a sign of anisotropy, but simply the result of projecting the two-dimensional EKE spectrum on zonal wavenumber (see appendix). Observed and predicted mean frequencies are computed as and .

Focusing first on the normalized frequency–wavenumber spectrum of the control experiment (middle panel of the middle row of Fig. 7), the barotropic dispersion curve and the reduced-gravity curve overlap with different energy-containing regions of the power spectrum— agrees better with the low-frequency part the spectrum, whereas traces the high-frequency, high-wavenumber end of the spectrum. In comparing the three simulations, as frictional strength varies from strong to weak, the maximum of the raw frequency–wavenumber spectrum (top row of Fig. 7) shifts from positive frequency to negative frequency. Correspondingly, in the phase speed spectra (bottom row of Fig. 7) the predominant phase speed shifts from eastward to westward when friction is varied from strong to weak values. The peak of the energy spectrum generally falls roughly within the region enclosed by the two neutral dispersion curves and . The phase speed of the baroclinic mode in its long-wave limit [proposed as a predictor of eddy phase speeds by Klocker and Marshall (2014)] agrees reasonably well with the energy peak in the control run and the strong friction run, but does not capture the fast westward propagation at weak friction. The faster westward propagation at weak friction is qualitatively consistent with the increased inverse energy cascade (Fig. 6) and barotropization (Fig. 5).

To quantify the overall accuracy of a particular linear model in a simulation, a normalized root-mean-square error (RMSE) for the phase speed is defined as

where *k*_{max} is a cutoff wavenumber, chosen to reduce the numerical bias due to limited sampling rates of the model output; higher wavenumbers tend to be associated with higher frequencies, which in turn require higher sampling rates. Based on our model parameter range and output sampling rate, *k*_{max} is chosen as 0.01 cycles per kilometer (cpkm). This range contains about 80%–90% of the total transverse EKE. The values of for the different linear models are given in the bottom panels of Fig. 7. Consistent with the visual impression, in the strong friction case, the reduced-gravity mode has the smallest normalized RMSE, whereas in the weak friction case the barotropic mode has the smallest normalized RMSE.

#### 2) Parameter space investigation

To examine how phase speeds vary in the THIC, FRIC, BETA, and BETA-hf experiments, for each experiment we calculate a diagnosed mean phase speed following

Notice that Eq. (22) differs from Eq. (19) by an additional integration over all zonal wavenumbers, thus yielding a single characteristic mean phase speed for each simulation.

Correspondingly, for each experiment, we calculate a predicted mean phase speed following

In Fig. 8, the top row compares diagnosed and predicted phase speeds (plotted on top of the power spectra), while the bottom row shows the growth rate as a function of phase speed from the linear stability analysis. For easy comparison, the general structure of Fig. 8 is similar to that of Fig. 6.

We first focus on the effect of varying layer thickness ratio in the THIC experiments (left column of Fig. 8). For a shallow upper layer (small *δ*), the EKE-containing eddies propagate eastward, which is close to the linear prediction from baroclinic instability analysis (bottom row). For a relatively shallow upper layer, the mean phase speed is also reasonably well captured by the predicted reduced-gravity curve . In the equal layer depth limit (*δ* ≈ 1), the observed mean phase speed is better captured by the barotropic mode predictor, . We also see an increase in the spread of the phase speed spectrum with increased layer depth ratio, which is indicative of increased turbulence leading to the dominance of chaotic advection.

We next focus on the role of varying the nondimensional planetary vorticity gradient in the BETA and BETA-hf experiments (right two columns of Fig. 8). For strong nondimensional beta (marginal criticality), the most energetic eddies generally propagate eastward, consistent with the eastward propagation of the linearly most unstable modes. As is reduced (criticality increases), the mean eddy phase speed in BETA shifts slightly westward and then slightly back eastward, while the mean eddy phase speed in BETA-hf shifts very slightly eastward across the entire parameter space. In the BETA experiments, as is reduced, there appear to exist two parameter regimes: First, reducing the nondimensional beta leads to barotropization, which leads to more westward phase propagation as the mean phase speed moves toward the phase speed of the barotropic mode. As is further reduced, the barotropic mode phase speed itself becomes weaker and eventually turns eastward. In the BETA-hf experiments, on the other hand, high friction suppresses barotropization and the mean phase speed remains well approximated by the reduced gravity mode through the entire range of . In both the BETA and (to a lesser degree) BETA-hf experiments, we also see an increase in the spread of the phase speed spectrum with decreasing . This is consistent with theoretical predictions, as the more supercritical mean state is expected to lead to more vigorous turbulence (e.g., Held and Larichev 1996).

By and large, varying the layer thickness ratio and nondimensional beta does not change the mean phase speed substantially. However, in the FRIC experiments (second column in Fig. 8), when friction becomes weaker, the predominant phase speed shifts significantly from eastward propagation to fast westward propagation. In the strong friction regime, the mean phase speed is well captured by the predicted phase speed for the reduced gravity mode. In the weak friction regime, on the other hand, the mean phase speed is well captured by the phase speed of the barotropic neutral mode, and is far away from the predictions of linear stability analysis (cf. the second row). As shown in Figs. 5 and 6, the weak friction regime is characterized by a robust inverse energy cascade and barotropization, indicating that nonlinear eddy–eddy interactions are the main reason for the departure of the eddy phase speeds from linear stability predictions. In spite of the strong nonlinearity in the weak friction limit, the linear barotropic-mode dispersion relation provides a good predictor for the eddy phase speeds, given the energy-containing wavenumber. The baroclinic and barotropic modes were also analyzed separately (not shown) and were found to have similar spectra, indicating tight coupling between modes.

The normalized RMSE in Fig. 9 confirms the aforementioned descriptions by quantifying the skill of the respective theoretical dispersion relations throughout the entire parameter space. In the THIC and FRIC experiments, the barotropic mode provides the most accurate prediction in one limit, whereas the reduced-gravity mode provides the most accurate prediction in the other limit. In the BETA and BETA-hf series, the reduced-gravity mode provides the best predictor throughout the parameter range.

We conclude that in the two-layer model, the strength of the inverse cascade determines the eddy length scale and amount of barotropization, which in turn determines the predominant eddy phase speed. When the flow is largely concentrated in the upper layer (e.g., due to strong bottom friction), the reduced-gravity dispersion curve *c*_{RG} is more apt, whereas when the flow is largely barotropic, the barotropic dispersion curve *c*_{BT} gives a better estimate. The baroclinic mode dispersion relation was found to be a generally less useful predictor than the reduced gravity mode.

## 4. Comparison to observations

In this section, we explore spatial and temporal spectral analysis of satellite observations and compare them to the QG model results discussed in the previous sections. The point of this exercise is to test the correspondence between the physical arguments in a two-layer QG model given above and the observed eddy propagation in a qualitative way. Our main test bed is a selected region of the Antarctic Circumpolar Current in the Southern Ocean, which is largely zonally symmetric and shares some essential properties with the two-layer QG model. As a comparison, we also revisit a subtropical region considered in Wunsch (2009), which is characterized by more “nondispersive” eddies and a shallow-mode instability that is not well represented by a zonally symmetric two-layer QG model (Tulloch et al. 2009).

### a. Data

We use 22 years of gridded geostrophic velocity anomalies and absolute geostrophic velocities produced by the SSALTO/DUACS (Developing Use of Altimetry for Climate Studies) multimission altimeter processing system and distributed by AVISO. Since the focus of this study is on the extratropical region, where the assumption of geostrophy is largely valid, we expect the AVISO products to reasonably reflect the underlying flow field. The grid size of AVISO is 0.25°, which in the ACC region corresponds to around 15 km. Low-pass Lanczos filtering has been applied to the product during the gridding procedure to remove the residual noise and small-scale signals below 65 km globally. Since in this study we focus on mesoscale eddies, which typically are larger than 200 km (Stammer 1997), the 65-km cutoff scale of the filter should not bias our core results, although some caution must be used in the interpretation of the data at higher wavenumbers.

For comparison with the two-layer QG model, the ideal test bed in observations should have a homogeneous zonal background flow and relatively small topographic effect on the eddies. We choose two regions with relatively flat topography and homogeneous background flow in a subtropical region and the ACC, respectively. Both regions are large enough compared to its eddy length scales. The subtropical region is in the northeast Pacific between 22° and 32°N, and from 180° to 140°W, centered around the cross section discussed in Wunsch (2009) (hereinafter the W09 box). The region from the ACC is between 54° and 64°S in the southeast Pacific from 117.375° to 77.375°W (hereinafter the ACC box). Figure 10 shows the climatological zonal component of the absolute geostrophic velocities for the two chosen regions, averaged for the year 2012. Both these regions are characterized by roughly zonally homogeneous mean flow. The ACC box region is at the core of the ACC and has a climatological zonal surface mean flow much stronger than the zonal surface mean flow in the W09 box region. Although both regions appear to be predominantly zonally symmetric, we note that even a small zonal variation of the basic state can have significant impacts on the instability problem (Tulloch et al. 2009). This may be of importance in particular for the W09 box region, where the zonal mean shear alone is stable to deep-mode baroclinic instability.

### b. Methods

To compare the observations with the results inferred from the two-layer QG model, we need to match the observed mean flow and stratification profiles to the relevant parameters in the QG model. A recipe for this procedure was developed by Flierl (1978) based on a vertical mode decomposition. If the flow in the continuously stratified fluid is dominated by the barotropic and first baroclinic mode, we can construct an analog two-layer model that adequately reproduces the dynamics of these two modes.

The first-baroclinic-mode structure of the continuously stratified QG equations provides a constraint for the baroclinic deformation radius and layer thickness ratio in the two-layer QG analog. Following Flierl (1978), we compute the layer thickness ratio to match the baroclinic self-interaction:

where the nonlinear baroclinic self-interaction, *ξ*_{111}, is determined based on the first baroclinic mode structure, Φ_{1}(*z*), as .

Following Tulloch et al. (2009), based on the vertical buoyancy profile, *N*^{2}(*z*), we calculate the first eigenmode structures, Φ_{1}(*z*), (right panel in Fig. 11) and eigenvalues, *K*_{d} = 1/*L*_{d}, from the Sturm–Liouville equation:

where the buoyancy–frequency profile *N*^{2}(*z*) is estimated from the neutral density data in Gouretski and Koltermann (2004) at the two chosen regions (left panel in Fig. 11). The right panel of Fig. 11 shows the first baroclinic mode in the two regions. In the ACC box, it has a deep structure, and a deformation radius of 14 km. The equivalent layer thickness ratio is found to be *δ* ≈ 0.49. On the other hand, the first baroclinic mode in the W09 box is more surface-intensified, with a radius of deformation of 43 km and an equivalent layer thickness ratio *δ* ≈ 0.23.

Figure 11 shows the vertical structure of the climatological zonal-mean zonal current in the W09 and ACC boxes. The mean flow for the W09 box is taken from the ECCO state estimate (Wunsch and Heimbach 2007), while the mean flow in the ACC box is taken from the Southern Ocean state estimate (Mazloff et al. 2010). In the ACC box, almost all zonal mean kinetic energy projects onto the barotropic mode and first baroclinic mode (blue and black curves in Fig. 11). This projection corresponds to a vertical mean flow shear of 2.8 cm s^{−1} in the two-layer QG analog, on the same order as the vertical mean flow shear considered in our QG two-layer reference setup in sections 2 and 3. Therefore the baroclinic instability in the ACC Box is expected to be dominated by deep “Phillips-type” instabilities, consistent with the findings of Tulloch et al. (2011).

In the W09 box, the projection on the two-layer QG model gives a mean state with a vertical mean flow shear of 0.22 cm s^{−1}, subcritical to a zonally symmetric two-layer QG-type baroclinic instability in the absence of friction. This points to several alternative routes through which eddies in this regions are generated: 1) a nonzero meridional shear, 2) instabilities of higher vertical modes, and 3) nonlocal generation. Even though a zonally symmetric two-layer QG model in this region is insufficient to directly generate baroclinic eddies, it might still be relevant to interpret the local behavior of eddies that may have been generated through other processes (so long as the surface EKE is dominated by the barotropic and first baroclinic mode). In particular the neutral reduced-gravity mode may be relevant for predicting eddy phase speeds in regions with strongly surface intensified stratification and weak vertical shear, such as the W09 region (e.g., Early et al. 2011).

### c. Phase speed spectra

Figure 12 shows a Hovmöller diagram of surface geostrophic velocity anomalies in the two regions. Visual inspection suggests that eddies in the W09 box propagate predominantly westward at a relatively coherent phase speed. In the ACC box, eddies propagate predominantly eastward, although it is less straightforward to identify one particular phase speed directly from the Hovmöller diagram. Instead, eddies seem to propagate at a range of different phase speeds. To quantify the phase speed spectra in the two boxes, we apply the same spectral analysis as for the numerical simulations discussed in section 3. Figure 13 shows the results of the spectral analysis (analogous in structure to Fig. 7). The predicted phase speeds are plotted on top of the spectra, with the two-layer model parameters chosen as described above.

The overall length scales for the ACC box (~0.025–0.004 cpkm) are slightly larger than the most energetic scales in our control two-layer QG simulation (~0.003–0.005 cpkm). This difference is attributable to the specific choice of parameters in the simulations, although shorter wavelengths may also be underrepresented in the AVISO data, due to the spatial smoothing. In the W09 box, the normalized transverse EKE spectrum (middle panel) shows robust nondispersive behavior that is consistent with the power spectrum of sea surface elevation shown in Fig. 8 of Wunsch (2009). On the other hand, the raw transverse EKE spectrum (upper panel) shows that the energy-containing eddies mostly fall on the high-wavenumber and high-frequency limit of the nondispersive line identified by Wunsch (2009). The phase speed spectrum (lower panel) for these energy-containing eddies does indicate some systematic variation in the phase speed with eddy scale. In contrast, in the ACC box, the transverse EKE spectrum exhibits dispersive behavior across a wide range of length scales. This resembles the control case in the reference simulations, which adopts parameters typical of the ACC (cf. Fig. 7).^{2}

Consistent with the results from our two-layer QG simulations, no single dispersion curve completely describes the diagnosed spectrum in either region. Both the raw and normalized transverse EKE power spectra are bounded by the barotropic dispersion curve and reduced gravity curve. It has previously been noted that the vertical structure of ocean eddies is often approximately equivalent barotropic, and lying between the barotropic and reduced gravity limits (e.g., Wunsch 1997; Arbic and Flierl 2004; Wortham 2013; Jansen et al. 2015), and our analysis supports this notion from the perspective of phase speeds.

In the ACC region, where the two-layer QG model is dynamically appropriate, the reduced gravity mode provides the best description of the dispersive wavenumber–frequency energy spectrum. In the subtropics, long waves fall roughly on a nondispersive line, as documented in Wunsch (2009), but phase speeds of the most energetic mesoscale eddies still exhibit a notable variation with wavelength suggestive of the influence of barotropic dynamics.

Interestingly, the long-wave limit of the reduced gravity mode bisects the peak in the energy spectra in both regions. (For the W09 box region, the baroclinic dispersion relation, and its longwave limit, is very similar to the reduced gravity mode, whereas in the ACC region the two differ more substantially.) However, it remains unclear whether this result is reflective of a fundamental physical process or merely arises from the combined influence of barotropic and baroclinic modes. In both regions the longwave-limit reduced gravity phase speed provides a better estimate at the most energetic scales than the longwave-limit baroclinic phase speed, which has been used in the past to fit eddy phase speeds (e.g., Klocker and Marshall 2014).

## 5. Conclusions

Our results demonstrate that, given the degree of barotropization and resulting eddy length scale, linear theory is useful to understanding eddy phase speeds in fully nonlinear two-layer QG baroclinic turbulence. Our results are consistent with those of Berloff and Kamenkovich (2013a,b), who also find linear dynamics to be relevant in a QG model. Our analysis of observed sea surface height variability further suggests that the QG model results are also relevant for the interpretation of observed mesoscale eddy phase speeds, at least in selected regions.

In the two-layer QG model, eddy phase speeds can be understood largely in terms of the reduced gravity mode and barotropic mode, which represent two relevant limit cases. If the inverse cascade is weak, the phase speed spectrum is reasonably well approximated by the phase speed of the linearly most unstable mode, which in turn can be approximated by the reduced gravity mode dispersion relation. A significant inverse cascade instead leads to barotropization, which in turn leads to phase speeds closer to those of barotropic-mode Rossby waves.

Analysis of satellite sea surface height observations suggests that mesoscale eddy phase speeds in the ACC and subtropics are also bounded by these two theoretical dispersion relations. Although two-layer QG dynamics can provide insight into the roles of linear baroclinic and barotropic modes on the eddy phase speed, such a model is clearly overly simplistic, particularly for the subtropical region. This work therefore represents only a step toward better understanding of eddy phase speeds in the real ocean.

## Acknowledgments

Much of this work was carried out at the 2014 Geophysical Fluid Dynamics Summer School while L.W. was a visiting student at Woods Hole Oceanographic Institution. M.F.J. and R.P.A acknowledge support from the GFD program. R.P.A acknowledges additional support from NASA Award NNX14AI46G. M.F.J. acknowledges additional support from NSF Award 1536450. We thank Glenn Flierl, William Young, Rui Xin Huang, and Isaac Held for helpful comments and discussions.

### APPENDIX

#### Advantages of Transverse Energy Spectra

The main advantage of using the zonal spectrum of meridional velocity (i.e., the “transverse” velocity spectrum) is that it better picks up the dominant scales in the two-dimensional total EKE(*k*, *l*) spectrum in Fig. 3 (see also Wortham et al. 2014). As shown in Fig. A1, the transverse EKE has vanishing power toward zero zonal wavenumber *k*, while the total EKE spectrum maintains high levels at zero *k* due to contributions from higher meridional wavenumbers. Comparison of Fig. A2 to Fig. 3 suggests that the meridional-mean transverse EKE spectrum exhibits similar peak zonal scales to the two-dimensional total EKE spectrum, whereas the meridional-mean total EKE spectrum exhibits peaks at somewhat smaller wavenumbers and contains significant energy at much smaller wavenumbers. At these small wavenumbers, however, most of the total energy comes from much larger meridional wavenumbers, which makes the interpretation of *ω*–*k* spectra confusing at best.

The zonal spectrum of meridional velocity further has the advantage that it is trivially related to the streamfunction variance spectrum, as the factor *k*^{2} can be taken outside of the integral in Eq. (16). Finally, it is the meridional eddy velocity that is most important for the transport of heat and tracers in the ocean (Klocker and Marshall 2014; Abernathey and Wortham 2015).

For the satellite observations, Fig. A3 verifies the largely isotropic behavior of the two selected boxes. As for the QG results, the meridional-mean transverse EKE spectrum captures the peak scales of the two-dimensional total EKE spectrum (Figs. A4 and A5). We hence adopt transverse spectra to interpret eddy phase speeds throughout this paper.

## REFERENCES

*An Introduction to Dynamic Meteorology.*3rd ed. Academic Press, 511 pp.

*Fundamentals of Geophysical Fluid Dynamics.*Cambridge University Press, 249 pp.

*Ocean Modeling in an Eddying Regime*, T. W. Hecht and H. Hasumi, Eds., American Geophysical Union, 5–15.

*Geophysical Fluid Dynamics*. Springer Verlag, 624 pp.

*Atmospheric and Oceanic Fluid Dynamics.*Cambridge University Press, 745 pp.

## Footnotes

^{1}

Formally the solutions in Eqs. (11) and (12) are derived assuming Δ*U* = 0 and thus *U*_{b} = *U*_{1} = *U*_{2}. There is thus some arbitrariness in the formulation of the “Doppler shift” in Eqs. (11) and (12). The formulations here were chosen such that the phase speeds remain independent of any background baroclinicity that may be present.

^{2}

We reemphasize that we do not mean to suggest that the true frictional drag in the ACC box is similar to the relatively large value used in the control reference simulation (i.e., 20 day^{−1}). In the two-layer QG model, frictional strength is a parameter for tuning the strength of the inverse cascade. In fact, in the two-layer QG model, friction is the only parameter that can be modified without affecting the other two nondimensional parameters. The two-layer QG model lacks many other elements that can influence the inverse cascade and barotropization in the real ocean, such as bottom topography and deep stratification. The net effect of all these factors collectively determines the strength of the inverse cascade and barotropization.