## Abstract

This study aims to understand the ocean’s circulation, which is characterized by the presence of multiple alternating zonal jets and transient mesoscale eddies, by systematic analysis of the underlying linear dynamics of this system. For this purpose, properties of the linear normal modes such as growth rates, dispersion, spatial structure, and nonlinear self-interactions are explored for a hierarchy of idealized, vertically, and horizontally sheared flows with increasing complexity. The authors find that large-scale background vertical shear, alternating multiple zonal jets, bottom friction, and the Reynolds number have important effects on these modes. This study hypothesizes that when these effects are taken into account, the linear results can be used to predict many properties of nonlinear mesoscale eddies. This hypothesis is confirmed in Part II of this paper.

## 1. Introduction

### a. Motivations and main questions

Mesoscale oceanic eddies are found virtually everywhere in the world’s oceans (e.g., Chelton et al. 2011), and they play an important role in ocean processes by affecting material transport, large-scale stratification and currents, air–sea interaction, etc. The last decades have seen remarkable progress in understanding physics governing eddies and their effects (e.g., review by McWilliams 2013), but many fundamental issues still remain unresolved. One of such issues is the subject of this paper.

The main goal of this two-part paper is to make some progress in understanding to what extent and how properties of eddies are controlled by the underlying linear dynamics. In the linear regime, the solution consists of a mean background state and linear modes, and the interactions between the modes are ignored. The linear control, that is, whether and how much predictions of the linear dynamics can be used for characterizing the nonlinearly equilibrated eddies, is a controversial topic. From the ocean observations, it is not clear whether eddies are strongly nonlinear (Chelton et al. 2011) or governed by the underlying linear dynamics (Wunsch 2009; Tulloch et al. 2009). Direct comparison of linear and nonlinear results is not straightforward. From the theoretical perspective, solving the full linear problem for the observed mean state of the ocean circulation and stratification is unfeasible. Bypassing this problem through applying local linearizations (e.g., Eden 2011) simplifies the problem, but the relevance and accuracy of this approach is questionable because there is no clear scale separation between inhomogeneities of the mean state and the linear normal modes. In this situation, the necessary first step is to examine the degree of the linear control in idealized eddy-resolving models solved at a large Reynolds number Re. However, so far even this approach remains largely unexplored and controversial. For example, some eddy-resolving studies argue that eddies are strongly nonlinear and cannot be adequately described by linear properties (e.g., Galperin et al. 2010), whereas some other studies report significant linear control of eddies (e.g., Berloff et al. 2009b). One of the main reasons for this controversy is in the appropriate definition of the mean state, which in many ways determines the properties of the linear solution. In particular, proper linearization must take into account alterations of large-scale potential vorticity (PV) gradients induced by eddies because such alterations may result in profound and spatially nonlocal transformations of the linear normal modes. The other reason for the controversy is that a systematic analysis of the linear properties should not be restricted to the most unstable linear normal modes, and should also consider other parts of the spectrum, which can be energized by the nonlinear interactions and be important in the nonlinear equilibration. Finally, we argue that the extent of the linear control must be established beyond the usual prediction of the length and time scales, by considering many other useful properties of the linear normal modes.

Our approach is to focus on an idealized, high-resolution model of anisotropic geostrophic turbulence, to analyze its solution, and to sort out the linear control question in a systematic and comprehensive way. The model equations are solved at the largest achievable Re to make their solutions the most nonlinear and turbulent. The model is baroclinic and driven by dynamically consistent instabilities of the imposed large-scale background flow—this formulation allows us to avoid imposing artificial and poorly constrained small-scale external forcing mimicking baroclinic-eddy dynamical effects. Because of the intrinsic anisotropicity, from the meridional variation of the Coriolis parameter and the imposed zonal flow, the eddying solutions develop multiple alternating zonal jets, which substantially modify the mean state and, thus, have profound effect on the underlying linear-dynamics properties. Accounting for the effects of the alternating jets places our study in the rich context of recent vigorous research on this subject (section 1b).

Our methodology in Part I of this paper aims to understand the underlying linear dynamics in terms of the linear normal modes and their properties. This is achieved by a systematic analysis of a hierarchy of background flows, from the simplest horizontally uniform currents to those that include the alternating jets. The corresponding linear normal modes are characterized by their dispersion relationships, horizontal and vertical structure, correlations with multiple jets, and nonlinear self-interactions, which can directly affect the jets. In Part II of this paper (Berloff and Kamenkovich 2013, hereafter Part II) these linear results are used to interpret the nonlinear solutions and to establish the degree of the linear control; and we argue that several important properties of the nonlinear eddies are very significantly controlled by the underlying linear dynamics.

The plan of the paper (Part I) is the following. The background is briefly outlined in section 1b; the formulation of the dynamical model and a general description of the solutions follow in section 1c. The linear analysis is in section 2, which deals with a hierarchy of linear problems. Sections 2a and 2b deal with uniform and idealized nonuniform background flows, respectively, and in section 2c we account for the realistic alternating jets by introducing the concept of “linear ensemble-averaged spectra.” In section 2d, we analyze the nonlinear self-interactions of the linear normal modes in the presence of the jets.

### b. Background

The focus of this study is on the anisotropic geostrophic turbulence, which is a dynamically rich system of interacting multiple alternating zonal jets, eddies, and isolated vortices. In general, the anisotropic geostrophic turbulence is controlled by the large-scale PV gradients that support planetary waves and by the large-scale sources of available potential energy that feeds growing instabilities. These instabilities result in the generation of mesoscale eddies that stir and redistribute the underlying PV, and this redistribution results in generation of the multiple jets that alter the mean PV distribution and thus feed back on eddies.

Over the last decade, observations of the multiple, quasi-zonal, alternating, latent jets in the global oceans came from the sea surface altimetry measurements, float trajectories, and tracer distributions, as well as from the eddy-resolving solutions of the primitive equation general circulation models [e.g., see the literature cited in Berloff et al. (2011)]. Reviews of the relevant theoretical ideas can be found in Rhines (1994), Dritschel and McIntyre (2008), and Berloff et al. (2011). The ubiquitous jets are formed because of anisotropic upscale energy transfers (e.g., Vallis and Maltrud 1993) but not necessarily through the energy cascade (e.g., Srinivasan and Young 2012). An important open question is to what extent the nonlinear dynamics of the jets and eddies is controlled by the underlying linear-wave properties of the large-scale PV distribution? Although there are indications that this control is significant at small Re (Berloff et al. 2009b; Yoo and Lee 2010), the issue was not explored in depth, particularly for large values of Re. Some properties of the most unstable linear normal modes of the alternating jets, as well as their roles in the nonlinear dynamics, were addressed in Berloff et al. (2009b, 2011). Here, we extend these preliminary results and explore the whole spectrum of the normal modes and in a comprehensive and systematic way.

Propagation of oceanic eddies represents another important and controversial aspect of the relative importance of linear dynamics. There are indications that significantly nondispersive propagation of the eddies is largely inconsistent with linear Rossby waves; these discrepancies led some authors to conclude that nonlinear dynamics is required to interpret propagating signals (Chelton et al. 2011; Early et al. 2011). However, an alternative point of view suggests that this nondispersive propagation can be explained using linear arguments, by strong coupling between the barotropic and baroclinic modes (Wunsch 2009). This issue calls for a direct comparison of the linear and nonlinear dynamical properties—this is the main subject of our two-part paper.

### c. Dynamical model

Our model of anisotropic geostrophic turbulence was widely used in the past for studying baroclinic multiple alternating jets (e.g., Haidvogel and Held 1980; Panetta 1993; Berloff et al. 2009a). The model is configured on the *β* plane, in a flat-bottomed zonal channel, with two isopycnal layers. The forcing is represented by an imposed, uniform zonal background flow with vertical shear; because this shear is unstable, it generates mesoscale eddies. The motivation for these idealizations is to establish a simple, but physically relevant, starting point, so that more physical complexity can be systematically added later on. However, Re is large, by oceanographic standards. Hence, solutions are very nonlinear. The model domain is a zonal channel with meridional width of *L*_{y} = 3600 km, and zonal extent *L*_{x} = 2*L*_{y}. The midchannel Coriolis parameter *f*_{0} = 0.83 × 10^{−4} s^{−1}, and the background planetary vorticity gradient *β* = 2 × 10^{−11} m^{−1} s^{−1}. There are two dynamically active isopycnal layers, and their thicknesses at rest are *H*_{1} = 1 km (upper layer) and *H*_{2} = 3 km (deep layer). The corresponding stratification parameters in the governing equations are

where is the reduced gravity coefficient associated with the density jump between the isopycnal layers. We chose so that the first baroclinic Rossby deformation radius , which is the fundamental length scale for the jets and ambient eddies (Berloff et al. 2009b), is 25 km.

The governing quasigeostrophic PV equations (e.g., Pedlosky 1987) are

where the layer index starts from the top; the *x* and *y* coordinates correspond to the zonal and meridional directions, respectively; *J*( , ) is the Jacobian operator; and the terms with *ν* and *γ* are the (Newtonian) lateral-eddy viscosity and the bottom friction, respectively. The eddy viscosity in our model represents momentum transfers by unresolved small-scale ageostrophic motions, and in the real ocean it is estimated to be on the order of 1–10 m^{2} s^{−1} (Muller 1976). The bottom friction parameterizes flow damping by the bottom boundary layer; in the real ocean, its value is expected to correspond to spindown time on the order of a month [see also Berloff et al. (2011), where effects of the bottom friction parameter are discussed in detail]. Isopycnal potential vorticities *q*_{i} are related to horizontal velocity streamfunctions *ψ*_{i} through the PV inversion:

No-slip lateral-boundary conditions are used on the sidewalls, and the mass and momentum constraints are imposed following McWilliams (1977). The forcing in the governing equations is introduced through an imposed, vertically sheared, baroclinically unstable background flow:

where parameters *U*_{i} are the background-flow zonal velocities, and we consider both east- and westward background flows.

For the flow consisting of the background-flow *U*_{i} and linear perturbations *ψ*_{i}, governing equations become

where advection of the PV anomalies by the background flow can be interpreted as the external forcing, which energizes some eddies by instabilities of the background flow. The perturbation PVs *ζ*_{i} are given by

The model is solved in terms of isopycnal-layer functions, but the solutions can be conveniently expressed in terms of barotropic *φ*_{BT} and baroclinic *φ*_{BC} vertical modes using the linear transformation

which is used for various analyses.

The bottom friction coefficient (i.e., *γ*), the eddy viscosity (i.e., *ν*), and the upper-layer background-flow velocity *U* = *U*_{1} are the key variable parameters, and we examined the following 12 solutions for two *U*, two *ν*, and three *γ* (Table 1). The two background flows yield dynamically different regimes (Berloff et al. 2009a): eastward background (EB) flow with *U* = 6 cm s^{−1} and westward background (WB) flow with *U* = −3 cm s^{−1}, respectively; *U*_{2} = 0 for both cases. Both east- and westward shears are moderately supercritical. Variation of *ν* allows us to study effects of the Reynolds number

which measures the nonlinearity of the flow; and the default hypothesis is that the larger Re, the weaker is underlying linear-dynamics control. Note that we based our definition of Re on fundamentally important scale Rd_{1}, rather than on the channel width, which is larger by two orders of magnitude. We consider *ν* = 50 and 2 m^{2} s^{−1}, which correspond to moderate and large (by ocean modeling standards) values of Re, respectively. In the moderate-Re case, WB and EB regimes are characterized by Re equal to 15 and 30, respectively; and in the large-Re case, WB and EB regimes have Re equal to 375 and 750, respectively. By varying the bottom friction coefficient, we obtain flow regimes with different degrees of latency of the alternating jets (Berloff et al. 2011): strong (nonlatent) jets, moderate (latent) jets, and no stationary jets. In the EB-flow regime, the small-, medium-, and large-*γ* cases correspond to *γ* = 10^{−8}, 10^{−7}, and 10^{−6} s^{−1}, respectively; and in the WB-flow regime, the corresponding cases are for *γ* = 10^{−8}, 4 × 10^{−8}, and 10^{−7} s^{−1}, respectively (Table 1).

The model equations are solved by a high-resolution numerical method (Karabasov et al. 2009), on the horizontal grid with 2048 × 1025 points (3.7-km resolution), and numerical convergence of these solutions was confirmed in Berloff et al. (2011) by several simulations with 4096 × 2049 grid spacing. For each set of parameters the model is integrated in time until a statistically equilibrated flow regime is reached (typically, after 20–40 yr of integration), and the next 60 yr are saved on a coarse 512 × 256 grid for analysis. In Fig. 1 instantaneous snapshots of the EB- and WB-flow solutions for large Re and intermediate *γ* illustrate the complexity of the flow.

We took special care in identification of statistically equilibrated-flow regimes to avoid sharp transitions of the linear properties. In particular, after initial-flow equilibrium is achieved following the first 10–20 yr of integration, some jets tend to merge over the next 20–40 yr. We examined latitude–time (Hovmöller) diagrams of zonally averaged flow (Fig. 2) to ensure that none of the solutions is contaminated by such jet mergers. Also, visual inspection of the latitude–time diagrams suggests that the jets are robustly present in the solutions at all times and most of the variability of the jet strength and position occurs on the annual and interannual time scales (Fig. 2). The latter implies that most of the interesting time scales of eddy variability are likely to be captured accurately by linearizations with respect to zonally averaged jets.

## 2. Linear analysis

In this section, we solve a hierarchy of linear problems—with and without multiple alternating zonal jets in the background state—focusing on the dispersion properties and growth/decay rates of the linear normal modes.

### a. Uniform background flow

First, we consider the simplest case of a horizontally uniform background flow (i.e., Phillips 1954) but in the presence of bottom friction and eddy viscosity. Past studies have focused mostly on the most unstable normal modes, whereas here we are interested in properties of all normal modes. Our background-flow velocity is predominantly baroclinic (the baroclinic is 4 times larger than the barotropic velocity component for all cases). Therefore, most of the background-flow effect on the dispersion properties is associated with the altered background PV, rather than with the barotropic Doppler shift of the phase velocities.

All the details of the linearization and the notations are in the appendix.

The resulting dispersion curves (Fig. 3) show frequencies—that is, the real parts of *ω*^{(1)}(*k*) and *ω*^{(2)}(*k*)—for a discrete set of 10 values of meridional wavenumber *l*, uniformly distributed over the same range of values as those used on this figure for zonal wavenumber *k*. When the frequency is negative, we trade its sign for the negative sign of *k*. Thus, normal modes with a negative/positive *k* have west-/eastward phase speeds. We used a value of *γ* = 10^{−7} s^{−1}, which in terms of our set of solutions is the medium value for the EB case but is a large value for the WB case. Three values of the background flow are discussed here: *U* = −3 (equivalent to the nonlinear WB case), *U* = 0, and *U* = 6 cm s^{−1} (equivalent to the nonlinear EB case); these values illustrate continuous changes found over this range of *U*. We classify the normal modes in terms of their growth rates (imaginary parts of the eigenfrequencies) *ω*_{i}: unstable (0 day^{−1} < *ω*_{i}), weakly damped (−0.001 < *ω*_{i} ≤ 0 day^{−1}), moderately damped (−0.002 < *ω*_{i} ≤ −0.001 day^{−1}), and strongly damped (*ω*_{i} ≤ −0.002 day^{−1}). With *U* = 0, unstable modes do not exist and all modes propagate to the west (Figs. 3b,e). The barotropic modes are damped to a larger degree than the baroclinic modes because the baroclinic modes have a smaller amplitude in the bottom layer and are less affected by the bottom friction. In the presence of the background flow, dispersion properties of the normal modes change profoundly. In the WB case, the modes retain westward phase velocities but become unstable in the region of weakly dispersive “tongue” seen in Figs. 3a and 3d. In the EB case, the first family of the modes rotates clockwise around the origin of the dispersion diagram and most of its members become eastward propagating (Fig. 3c). Also, most of these modes, except for those with relatively large *k* and *l*, become unstable. Some of the second-family modes also become eastward propagating, but none of them become unstable (Fig. 3f). Thus, an important difference between the WB and EB regimes is that in the former one there are two types of unstable modes, whereas in the latter one there is only one type of unstable mode.

In the absence of the background flow and bottom friction, the first family of the modes is purely baroclinic, that is, it projects [see (10)] only on the baroclinic vertical mode, whereas the second family is purely barotropic, that is, it projects only on the barotropic vertical mode. A common view is that the ocean is populated by purely barotropic and baroclinic eddies that become substantially coupled by nonlinearity in nonlinear eddies. To what extent is this barotropic–baroclinic segregation affected by the presence of a background flow and bottom friction (the lateral friction has no effect on this aspect)? To address this question quantitatively, we introduce a vertical-mode mixing coefficient *μ* = *μ*(*k*, *l*), that is defined as *φ*_{BT}/*φ*_{BC} for the first family and as *φ*_{BC}/*φ*_{BT} for the second family of the modes (this particular ordering is not very important): *μ* = 0 or *μ* = ∞ implies that the normal mode is “pure,” that is, it projects on the single vertical mode.

Modes with 0.25 < *μ* < 4.0 are referred to as significantly mixed normal modes, and we examined how the fraction of the spectrum corresponding to this group depends upon parameters of the problem. We varied *U* and *γ* (Fig. 4) and found that this fraction tends to increase with |*U*|, and it is 30%–40% even for |*U*| = ~0.5 cm s^{−1}, which is a relatively small value for oceanic currents. There is noticeable asymmetry between negative and positive *U*, and the WB flow always results in stronger vertical-mode mixing. Bottom friction significantly enhances vertical-mode mixing, even in the absence of the background flow, and with nonzero *γ* the second family of normal modes has systematically larger values of *μ* than the first family. Overall, the analysis of *μ* suggests that the presence of a background flow or bottom friction imposes strong coupling between the barotropic and baroclinic modes, even in the absence of nonlinear coupling. A more detailed picture of vertical-mode mixing is given by plotting *μ*(*k*, *l*) (Fig. 5) for the set of dispersion curves shown in Fig. 3. In this plot, we sort out the modes into strongly baroclinic (*μ* < 0.4), weakly baroclinic (0.4 < *μ* < 0.8), well-mixed (0.8 < *μ* < 1.25), weakly barotropic (1.25 < *μ* < 2.5), and strongly barotropic (2.5 < *μ*). With *U* = 0 and *γ* = 10^{−7} s^{−1}, the first family is strongly baroclinic, and the second family is mostly strongly barotropic, except for the moderately barotropic modes characterized by relatively large values of *l* and *k*. In the WB case, there is substantial diversity of vertical structures of the modes within both families (Figs. 5a,d). In the EB case similar diversity is found only for the second family, whereas the first family remains weakly baroclinic. The most unstable WB modes are the following: weakly and strongly baroclinic in the first family, and well-mixed and weakly barotropic in the second family; and the most unstable EB modes are weakly baroclinic. The fast westward-propagating modes, which are descendants of purely barotropic modes in the *U* = *γ* = 0 case, also become substantially mixed, except for the fastest ones.

In the two-layer system (characterized by 2 vertical degrees of freedom) two quantities are needed for the description of the vertical structure. Above, we discussed ratio of the amplitudes of the barotropic and baroclinic modes, and now we turn our attention to their phase relationship. For this purpose, we calculated correlation between the barotropic and baroclinic components of each mode:

This correlation quantifies coupling between the vertical modes, and, thus, serves as a convenient second measure of vertical-mode mixing. A large positive number suggests surface-intensified or bottom-compensated flow; such modes can be expected to be less affected by the bottom friction. In the EB-flow case, *C*_{vert} is mostly positive and large, and the only significant concentration of its negative values corresponds to the eastward-propagating normal modes with short *k*. In the WB-flow case, *C*_{vert} is mostly weakly negative, but in some parts of the *k*–*ω* space neighboring dispersion curves have *C*_{vert} of either sign.

We find that distribution of *μ*(*k*, *ω*) is rather intermittent, in the sense that neighboring branches of dispersion curves can be characterized by very different values of *μ*. Hence, provided significant linear control, estimation of *μ* from *k*–*ω*-filtered time series of the nonlinear solution will unavoidably mix up contributions from the normal modes with different vertical structures. On the contrary, distribution of *C*_{vert}(*k*, *ω*) is substantially more uniform. Therefore, it is more applicable for *k*–*ω* spectral diagnostics of the nonlinear solutions (Part II).

The effect of eddy viscosity is weak (in the considered range of wavenumbers) and limited to enhanced damping, particularly at large wavenumbers. The bottom friction stabilizes the second family but further destabilizes the unstable tongue in the WB case (cf. Figs. 3 and 6). The effect on the first family of modes depends on the background-flow direction: in the WB case it is stabilizing, but in the EB case it is destabilizing. We find that the largest considered value of *γ* = 10^{−6} s^{−1} changes only growth rates and makes stable modes even more damped and unstable modes even faster growing; and the corresponding changes of the shapes of the dispersion curves are minimal (not shown).

In this section we focused on the linear analysis of uniform background flows and showed that (i) in the WB case all unstable modes are westward propagating, whereas in the EB case there are unstable modes propagating in either direction; (ii) the unstable modes occupy clearly defined regions of the *k*–*ω* space, and in the WB case they form the characteristic nearly nondispersive tongue; (iii) most of the modes are significantly mixed in vertical, even for weak background shears, and the bottom friction enhances this mixing; and (iv) large bottom friction increases the contrast between the growth rates of the normal modes, by damping stable modes and destabilizing unstable modes even more, and by making stable/unstable modes more barotropic/baroclinic.

### b. Nonuniform background flow

In this section, we consider a more complex background flow with idealized multiple alternating zonal jets on top of the uniform WB and EB flows. The linearized problem is formulated in the appendix. To examine the effects of the jets, we varied their amplitude from zero to the values observed in the nonlinear solutions. The numerical eigenvalue problem solver does not discriminate between the first and second families of the normal modes (section 2a), because the analytical tractability is lost. The new normal modes are presented for the background flows corresponding to the WB- and EB-flow regimes, and for the situations with weak (*U*_{1} = 0.5 and 1.0 cm s^{−1} for the WB and EB cases, respectively) and moderate (*U*_{1} = 3.0 and 6.0 cm s^{−1} for the WB and EB cases, respectively) jets (Fig. 7). As in the case of the uniform flow, the modes are sorted out by their growth rates into the following: unstable (0 day^{−1} < *ω*_{i}), weakly damped (−0.0005 < *ω*_{i} ≤ 0 day^{−1}), moderately damped (−0.001 < *ω*_{i} ≤ −0.0005 day^{−1}), and strongly damped (*ω*_{i} ≤ −0.001 day^{−1}). There are several new effects because of the presence of the jets. First, there is clustering of the dispersion curves; that is, in some *k*–*ω* regions there are many normal modes and the dispersion curves overlap and cluster together. Second, most of the normal modes are meridionally localized on either east- or westward jets, and we discuss this in more detail in section 2c. The localization was reported previously for the most unstable modes (Berloff et al. 2009b), but here we observe it for most of the modes. Third, there is significant stabilization of many dispersion curves, and this is indicated by reduction of the unstable (red) curves in Figs. 7b and 7d relative to Figs. 7a and 7c. We conjecture that the underlying stabilization mechanism is similar to the one proposed by James (1987), who argues that a barotropic shear can suppress baroclinic instability by localizing the critical disturbances. Finally, the jets result in a broader range of phase speeds of the normal modes, and this is due to the fact that the modes localized on the east-/westward jets propagate faster to the east/west, both from the barotropic Doppler shift and the altered PV gradients. For example, appearance of eastward-propagating modes in the WB case is only due to the jets.

We studied the effects of bottom friction by changing *γ* from 10^{−8} s^{−1} (Figs. 7b,d) to 10^{−7} and 10^{−6} s^{−1} (Fig. 8). Only strong bottom friction has an effect by noticeably reducing the number of unstable and marginally stable normal modes, by wiping out some of the westward-propagating nearly barotropic modes (except for the fastest ones), and by making most of the normal modes nearly baroclinic rather than vertically mixed. Overall, we argue that the indirect effect of the bottom friction associated with the weakening of the jets (Berloff et al. 2011) is more important than the direct friction-induced modifications of the normal modes, which are obtained by linearizing around the pronounced jets.

### c. Linear ensemble-averaged spectra

Up to now, all linear analyses of section 2 focused on idealized-flow configurations. We now proceed with linear analysis of the time-dependent background flow diagnosed directly from the nonlinear solutions. The solution is linearized around a slowly varying mean (zonally uniform) state, assuming that the time scale of linear modes is much faster than the evolution of the mean state. As demonstrated below, this approach leads to results most closely corresponding to the nonlinear results (Part II). To follow this approach, we developed the concept of linear ensemble-averaged spectrum, which is the *k*–*ω* spectrum worked out from an ensemble of dispersion relationships. Each dispersion relationship is calculated as in the previous section, but for the mean state taken to be the zonal average of the full nonlinear-model flow, in 2-day snapshots. The *k*–*ω* domain is binned into 400 × 100 bins, and the normal-mode growth rates corresponding to the ensemble of dispersion curves are mapped into these bins. If at least one dispersion curve passes through a bin, then the largest growth rate is assigned to this bin, otherwise the bin is left “white,” that is, without information. This technique provides an upper bound on the growth rates; we have also tried averaging growth rates in each bin, but the outcome of this is qualitatively similar.

We focus the presentation on four families of solutions calculated for moderate (EB: Fig. 9; WB: Fig. 10) and large (EB: Fig. 11; WB: Fig. 12) Re. The figures show the corresponding linear ensemble-averaged *k*–*ω* spectra, and the main features of these plots are the following. First, parts of the spectra populated by the normal modes are much broader at large Re, due to the appearance of faster eastward- and westward-propagating modes meridionally localized on the east- and westward jets, respectively. Second, there are well-defined parts of the spectra corresponding to growing modes, outlined by quadrangles for future use in the nonlinear analysis (Part II). Third, around the unstable (shown in red) parts of spectra, there are much broader regions corresponding to nearly neutral (or, weakly damped) modes (shown in green). These modes can be easily energized by eddy–eddy nonlinear interactions; as shown in Part II, they play an important role in the nonlinear spectra. Fourth, at moderate Re some of the individual dispersion curves are visible in the ensemble spectra, but at large Re they are smeared out owing to more substantial variability of the background jets. Finally, we notice weakly damped westward-propagating and fast-growing eastward-propagating modes in the EB case with strong bottom friction, and no eastward-propagating modes in the WB case—all of this is in accord with findings of sections 2a and 2b.

### d. Eddy forcing of the normal modes

Analyses of the previous sections identified important effects of the jets on properties of the linear normal modes. Our next step is to explore feedbacks of the modes on the underlying jets by using the concept of internally generated eddy forcing (i.e., nonlinear self-interaction of a normal mode). In the nonlinear regime, normal modes can interact with each other through resonant triad interactions (e.g., Pedlosky 1987). In particular, self-interactions of normal modes can project on purely zonal (*k* = 0) stationary (*ω* = 0) modes corresponding to the jets. In this case, zonal jets can be maintained/resisted by eddy forcings that positively/negatively correlate with the jets. We have analyzed the normal-mode eddy forcings for all flow regimes considered here, but because the general outcomes are qualitatively similar, the discussion below is limited only to the large-Re EB case with moderate bottom friction. The background flow corresponds to Fig. 1a, and the corresponding linear ensemble-averaged spectrum is shown in Fig. 11b.

First, we characterize meridional positioning of the normal modes relative to the jets and find that most of the normal modes straddle (i.e., are meridionally localized on) either east- or westward jets. For this purpose, we calculate spatial correlations between these modes and the background flow over those latitude intervals where the mode amplitude is significant. More precisely, we normalize the meridional structure of each mode by its maximum amplitude and consider latitude intervals on which the normalized amplitude is greater than 0.02 in absolute value. Over the set of these intervals, we correlate the meridional profile of each mode with the jets in each layer and average over two layers. This correlation is applied only over the region of significant amplitude because the normal modes tend to be localized on individual jets; therefore, correlating them over the entire domain would bias correlations toward small values. (In the nonlinear analysis of Part II such conditioning of the eddy amplitude is not needed, because at a given wavenumber and frequency all jets are straddled by the eddies.) The correlations between the normal modes and jets show that most of the modes are localized on either east- or westward jets, thus, resulting in clear zonation of the *k*–*ω* spectrum in terms of this property (Fig. 13). The largest correlations (Fig. 13) correspond to the *k*–*ω* regions with the largest growth rates of the normal modes (section 2c), thus, indicating that potentially the most active normal modes are also the most meridionally localized on the jets. Short waves (large *k*) straddling eastward jets propagate eastward; longer waves can move in either direction. The modes straddling the westward jets tend to have larger zonal wavenumbers than the ones straddling the eastward jets and tend to propagate in both directions.

We next analyze the spatial patterns of the eddy forcings associated with the normal modes and project them on the jets. The eddy forcing is convergence of the eddy PV flux, which consists of fluxes of relative vorticity

and isopycnal stretching (or buoyancy anomaly)

The time- (indicated by overbar) and zonally averaged (indicated by angular brackets) eddy forcing

can be split into divergence of the relative vorticity flux [i.e., Reynolds stress (RS) forcing] and divergence of the buoyancy flux [i.e., form stress (FS) forcing]. Next, we define barotropic *F*_{BT} and baroclinic *F*_{BC} eddy forcings by projecting the governing equations on the vertical barotropic and baroclinic modes (Berloff et al. 2009a). The barotropic-eddy forcing is a sum of Reynolds stresses produced by barotropic–barotropic and baroclinic–baroclinic nonlinear interactions, and we look at these components individually. The baroclinic-eddy forcing is a sum of the Reynolds (barotropic–baroclinic and baroclinic–baroclinic) and form stresses, and we also look at these components individually.

For each normal mode, we calculated spatial correlations between *F*_{BT}(*y*), *F*_{BC}(*y*), and their components, and the corresponding meridional PV anomalies associated with the barotropic and baroclinic components of the jets, respectively. These correlations measure efficiencies of the corresponding eddy forcings and their components for maintaining (positive correlations) or resisting (negative correlations) the jets. As in the algorithm described above, each correlation is only calculated over the set of intervals on which normalized amplitude of the normal mode is significant (i.e., larger than 0.02; varying this number by a factor of 2 makes no noticeably changes in the outcome).

The main result of the above analysis is an understanding of how the *k*–*ω* spectrum is partitioned in terms of the fundamental properties of the eddy forcing. In particular, Figs. 14a and 14b show which eigenmodes have barotropic-eddy forcing (i.e., *F*_{BT}) with positive/negative feedback on the barotropic component of the jets; and Figs. 14c and 14d show the same relationship, but for the baroclinic-eddy forcing (i.e., *F*_{BC}) and the baroclinic component of the jets. From these plots, we find that most of the unstable eastward-propagating modes that straddle the eastward jets act strongly to maintain the barotropic jet but very few of those modes maintain the baroclinic jets. In contrast, most of the unstable westward-propagating modes that straddle the westward jets weakly resist the barotropic jets, but strongly maintain the baroclinic jets. There is, however, a well-pronounced “wedge” of the westward-propagating modes, characterized by relatively fast phase speeds and low wavenumbers, which strongly resists the baroclinic jets (Fig. 14d). As for the high-frequency, low-wavenumber, westward-propagating, near-barotropic modes, they have weak correlations with the eastward jets and tend to maintain both the barotropic and baroclinic jets. Overall, most of the normal modes tend to support either barotropic or baroclinic (or both) jets. The jet-resisting modes, most of which resist either barotropic or baroclinic jets but not both of them, tend to populate lower frequencies and wavenumbers.

The effects of the normal-mode RS and FS forcings can be different. Textbook theory on the baroclinic-eddying jets (e.g., Panetta 1993) states that in the EB flows, the jets are maintained by the RS and resisted by the FS, and the former process prevails. This RS effect is often described in terms of “negative viscosity” phenomenon, whereas the FS effect results from the baroclinic instability of the jets. For the WB flows, however, Berloff et al. (2009a,b) showed that this scenario is not valid, and the RS and FS eddy forcings exchange their roles, so that the RS resists the jets and the FS maintains them, and showed that this effect can be traced to the properties of the most unstable normal modes. Here, we extend these preliminary results to the whole *k*–*ω* spectrum of the normal modes, by calculating all modal FS and RS eddy-forcing components and their correlations with the underlying baroclinic jets (Fig. 15). The outcome is that we found clear partitioning of the spectrum into several dynamically distinct zones corresponding to different classes of the normal modes. We found that the whole family of the most unstable eastward-propagating modes maintains the baroclinic jets by the RS and resists them by the FS, both in the EB- and WB-flow regimes. On the other hand, the family of the most unstable westward-propagating modes maintains the jets by the FS and weakly resists them by the RS. In general, the normal modes with lower frequencies have jet-maintaining FS, whereas the eastward-propagating modes with higher frequencies tend to have jet-resisting FS.

## 3. Conclusions

The main goal of this paper is to make further progress in understanding to what extent and how properties of oceanic eddies are controlled by the underlying linear dynamics. It is still not clear whether the observed eddies are strongly nonlinear (e.g., Chelton et al. 2011) or significantly controlled by the underlying linear dynamics (e.g., Wunsch 2009). There is a similar controversy on the theoretical side, with idealized numerical-modeling studies arguing either for (Berloff et al. 2011) or against (Galperin et al. 2010) significant linear-dynamics control exerted on eddies. There are two main reasons for this controversy. First, the linear dynamics is governed by the interactions between transient motions (eddies) and steady or slowly evolving background (mean) state, and a dynamically consistent definition of the mean state is essential. In particular, we demonstrate that alterations of the large-scale potential vorticity gradients induced by eddies must be taken into account, and that the realistic mean state can dramatically modify horizontal and vertical structure as well as dispersive behavior of the linear solutions. Second, we argue that the level of the linear control must be established far beyond the traditional linear-stability prediction of the eddy length and time scales, by considering not only many other properties of the linear normal modes, but also all parts of their full spectrum that can be energized by nonlinear interactions and feedback on the mean state through mode self-interactions.

Our approach is the following. We consider solutions of an idealized, turbulent nonlinear model of mesoscale eddies generated by unstable, vertically sheared zonal flows at moderate and large Reynolds numbers (i.e., Re). Zonally averaged instantaneous solutions contain multiple alternating zonal jets that exist because of the action of eddies and can strongly modify eddies, as well as the underlying linear properties. Properties of the jets are used in a hierarchy of linear problems, which we solve and analyze. These problems start from the uniform background flow, then incorporate alterations of this flow by eddies, in the form of simple idealized jets, and then the jets observed in the nonlinear solutions. In each case the corresponding linear normal modes and their zonal wavenumber/frequency dispersion relations are found and referred to as “linear spectra.” The resulting linear spectra, as well as other properties of the normal modes, are used for interpreting the nonlinear spectral-analysis results of Part II.

The linear normal-mode properties are characterized by the dispersion relationships, horizontal patterns, vertical structure, correlations with the multiple jets, as well as by the nonlinear modal self-interactions and their projections on the jets. We explored how all of these properties depend on the mean state and frictional parameters. First, we found that in the presence of even weak uniform background flow, the vertical structure of most linear normal modes ceases to be purely barotropic or baroclinic; and instead it becomes “mixed.” This result strongly suggests that flow decomposition into the barotropic and baroclinic eddies is dynamically inconsistent in most parts of the global ocean. This implies that considering separate barotropic and baroclinic spectral energy transfers, and even cascades (e.g., Scott and Arbic 2007), may need to be upgraded to account for coupling between the vertical modes. Second, we find that in the presence of the westward (vertically sheared) background flow, the unstable linear normal modes have a tendency to be nearly nondispersive and westward propagating. In the eastward background flow, distribution of the unstable modes is more complicated, but because the westward flows are a prominent part of the wind-driven gyres, it is tempting to hypothesize that the observed nondispersive westward propagation of most of the eddies (Wunsch 2009; Chelton et al. 2011) can be explained by the linear control. We found that, for the parameters of the problem, the unstable and marginally stable linear normal modes occupy large regions of the spectral domain, and in Part II we show that these are the regions of significant spectral power of the nonlinear eddies. Next, we identified spectral regions with maximum growth rates, under the hypothesis that in the nonlinear solutions these regions correspond to the maximum energy input from the background flow to the eddies; this hypothesis is confirmed in Part II.

We explored different flow regimes in order to understand the roles of the vertical and meridional shears in the linear control and in shaping the normal modes. For this purpose, we considered flows with (i) manifest jets (low friction), (ii) moderately latent jets (medium friction), and (iii) very latent transient jets [high friction; see a discussion of these flow regimes in Berloff et al. (2011)]. We also studied less and more turbulent- and nonlinear-flow regimes, by considering intermediate and large, by the oceanographic standards, values of Re. The Reynolds number effect is mostly due to the fact that the jets become stronger, when the eddy viscosity becomes smaller. We found that the jets meridionally localize most of the normal mode amplitudes on either east- or westward jets. This process strongly modifies the spatial structure of the normal modes and changes their dispersion properties. Higher Re also expands those parts of the spectrum that are occupied by the unstable and marginally stable normal modes.

Our linear approach effectively uncouples nonlinear dynamical processes, which is the important advantage in such a complex system. In particular, in addition to the analysis of the influence of the mean state on eddies, we analyzed the eddy feedback on the mean state by finding correlations between the normal modes and the jets, as well as correlations between the eddy forcing (i.e., nonlinear self-interaction of the mode) and the jets. We found the following important properties: (i) the modes straddling the eastward jets tend to have shorter zonal wavenumbers and larger frequencies; (ii) the modes with eddy forcing that maintains/resists the barotropic jets tend to have faster/slower phase speeds; (iii) most of the modes have eddy forcing that maintains, rather than resists, the baroclinic jets; (iv) the modes that have the eddy forcing resisting the baroclinic jets via the form stress tend to be fast and eastward propagating, whereas the slower eastward-propagating and most of the westward-propagating modes tend to have the form stress acting in the opposite way; and (v) the modes that have the eddy forcing resisting the barotropic jets via the Reynolds stress tend to have both short zonal wavenumbers and low frequencies, whereas the modes with relatively high frequencies, and even more so those propagating to the east, tend to have the Reynolds stress acting in the opposite way. We hypothesize that in the nonlinear solutions, the spectral power accumulated in the modes of certain type defines the roles of the baroclinic and barotropic-eddy forcings, and of the form and Reynolds stresses. In the classical problem with “negative viscosity” effect—an eddy-driven midlatitude eastward jet—the mean jet is maintained by the eddy Reynolds stresses and resisted by the eddy form stresses; however, our results suggest that this mechanism is not generic, and there can be other alternatives [see also Berloff et al. (2009a)]. It will be interesting to find simple and clear physical arguments explaining spectral zonation of the normal modes in terms of their eddy-forcing properties, but this is beyond the scope of the present study.

## Acknowledgments

PB acknowledges hospitality of the Department of Applied Mathematics and Theoretical Physics at the University of Cambridge, where most of this work has been done. Funding was also provided for PB by NSF Grant OCE 0845150, and for IK by NSF Grant OCE 0842834 and 1154923. Numerous discussions with Tom Farrar are graciously appreciated.

### APPENDIX

#### Linearization of the Governing Equations

In the case of horizontally uniform, zonal background flow, with the velocities *U*_{1} and *U*_{2} in the upper and lower layer, respectively, the linearization of the governing equations around the background flow yields:

These equations are Fourier transformed in space and time, and their solutions are sought in the full spectral form:

where *k* and *l* are zonal and meridional wavenumbers, respectively, the corresponding eigenfrequency is *ω*, and the corresponding eigenvector is . The Fourier-transformed (A1) and (A2) are as follows:

where *κ* = *k*^{2} + *l*^{2}. For existence of a solution, the determinant of this system of linear equations must be zero. This condition yields the following quadratic equation for *ω*:

For each pair of *k* and *l*, there is unique solution of (A6) given by the pair of eigenvalues, [*ω*^{(1)}, *ω*^{(2)}]. Substitution of *ω*^{(n)} into (A4) and (A5) yields the solution pair, and , for *n* = 1, 2:

We refer to the first (*n* = 1) and second (*n* = 2) pair of solutions, as the first and second families of the eigenmodes. For consistency with the nonlinear model, in the following analysis we set *U*_{1} = *U* and *U*_{2} = 0.

In the presence of multiple alternating zonal jets on top of the uniform background flow, the background-flow-velocity streamfunction is

where *u*_{i}(*y*) is either idealized or dynamically consistent zonal velocity profile. For the latter, we used ensembles of zonally averaged solutions of the nonlinear model (section 1c). For the idealized velocity profile, we assume meridional periodicity; hence, the problem is treated as doubly periodic and there are no zonal walls. For simplicity, we used a sinusoid with the period specified by *N* and the amplitude specified by :

We used *N* = 8 for the EB case and *N* = 7 for the WB case, as those are approximate fits to the nonlinear solutions with medium values of the bottom friction (Fig. 1). To reduce the number of parameters in this problem, we also fixed ratio , and treated *U*_{1} as the only control parameter for the background-flow jets.

The perturbation streamfunction is Fourier transformed only in time and *x*:

and the Fourier-transformed linearized equations are as follows:

In principle, the eigenvalue problem can be solved analytically by looking for solutions in terms of the meridional Fourier series, but we have not done this. The above equations were discretized with second-order finite differences, the no-slip boundary conditions were applied on the sidewalls, and the resulting generalized eigenproblem was solved numerically. We tested that the outcome is not sensitive to further refinement of the spatial grid. The solutions were obtained in terms of the eigenmodes for each zonal wavenumber and complex frequencies (eigenvalues).

## REFERENCES

*J. Phys. Oceanogr.,*

*Geophysical Fluid Dynamics*. 2nd ed. Springer-Verlag, 710 pp.