## Abstract

This study investigates the representation of solutions of the three-dimensional quasigeostrophic (QG) equations using Galerkin series with standard vertical modes, with particular attention to the incorporation of active surface buoyancy dynamics. This study extends two existing Galerkin approaches (*A* and *B*) and develops a new Galerkin approximation (*C*). Approximation *A*, due to Flierl, represents the streamfunction as a truncated Galerkin series and defines the potential vorticity (PV) that satisfies the inversion problem exactly. Approximation *B*, due to Tulloch and Smith, represents the PV as a truncated Galerkin series and calculates the streamfunction that satisfies the inversion problem exactly. Approximation *C*, the true Galerkin approximation for the QG equations, represents both streamfunction and PV as truncated Galerkin series but does not satisfy the inversion equation exactly. The three approximations are fundamentally different unless the boundaries are isopycnal surfaces. The authors discuss the advantages and limitations of approximations *A*, *B*, and *C* in terms of mathematical rigor and conservation laws and illustrate their relative efficiency by solving linear stability problems with nonzero surface buoyancy. With moderate number of modes, *B* and *C* have superior accuracy than *A* at high wavenumbers. Because *B* lacks the conservation of energy, this study recommends approximation *C* for constructing solutions to the surface active QG equations using the Galerkin series with standard vertical modes.

## 1. Introduction

Recent interest in upper-ocean dynamics and submesoscale turbulence has focused attention on surface geostrophic dynamics and the role of surface buoyancy variations. A main issue is the representation of active surface buoyancy by finite vertical truncations of the quasigeostrophic (QG) equations. Standard multilayer (e.g., Pedlosky 1987) and modal approximations (e.g., Flierl 1978) assume that there is no variation of buoyancy on the surfaces.

Only few attempts have being made to represent both surface active and interior dynamics in the QG equations. The pioneering work by Tulloch and Smith (2009) developed a “two-mode two-surface” model that represents the surface dynamics exactly and approximates the interior dynamics using the barotropic and first baroclinic modes. The interaction of surface and interior dynamics motivated the development of a new set of vertical modes that simultaneously diagonalize energy and a linear combination of enstrophy and surface buoyancy variance (Smith and Vanneste 2013). Other studies of the interaction of surface and interior dynamics avoid vertical modes and use instead finite-difference schemes (Roullet et al. 2012) or idealize the interior potential vorticity as a delta function sheet (Callies et al. 2015).

Here, we explore the representation of surface and interior dynamics using the familiar vertical modes of physical oceanography. These “standard modes,” denoted here by *p*_{n}(*z*), are defined by the Sturm–Liouville eigenproblem

with homogeneous Neumann boundary conditions at the bottom (*z* = *z*^{−}) and top (*z* = *z*^{+}) surfaces of the domain:

In (1), *N*(*z*) is the buoyancy frequency and *f*_{0} is the Coriolis parameter. The eigenvalue *κ*_{n} in (1) is the deformation wavenumber of the *n*th mode. With normalization, the modes satisfy the orthogonality condition

where is the depth. The barotropic mode is *p*_{0} = 1 and *κ*_{0} = 0.

The modes defined by the eigenproblems (1) and (2) provide a fundamental basis for representing solutions of both the primitive and quasigeostrophic equations as a linear combination of {*p*_{n}} (Gill 1982; Pedlosky 1987; Vallis 2006; Ferrari and Wunsch 2010; LaCasce 2012). In fact, the set {*p*_{n}} is mathematically complete and can be used to represent any field with a finite square integral:

Even if the function *f*(*z*) has a nonzero derivative at *z*^{±}, or internal discontinuities, its representation as a linear combination of the basis functions {*p*_{n}} converges in *L*^{2}(*z*^{−}, *z*^{+}), that is, the integral of the squared error goes to zero as the number of basis functions increases (e.g., Hunter and Nachtergaele 2001, chapter 10). In quasigeostrophic dynamics, both the streamfunction *ψ* and the potential vorticity (PV) *q* satisfy the requirement (4), and thus both *ψ* and *q* can be effectively represented by linearly combining {*p*_{n}}.

Despite the rigorous assurance of completeness in the previous paragraph, the utility of {*p*_{n}} for problems with nonuniform surface buoyancy has been questioned by several authors (e.g., Lapeyre 2009; Roullet et al. 2012; Smith and Vanneste 2013). These authors argue that the homogeneous boundary conditions in (2) are incompatible with nonzero surface buoyancy and that representation of the streamfunction *ψ* as a linear combination of {*p*_{n}} is useless if *ψ*_{z} is nonzero on the surfaces.

The aim of this paper is to obtain a good Galerkin approximation to solutions of the QG equation with nonzero surface buoyancy using the familiar basis {*p*_{n}}. We show that that both the inversion problem and evolutionary dynamics can be handled using {*p*_{n}} to represent the streamfunction. As part of this program we revisit and extend two existing modal approximations (Flierl 1978; Tulloch and Smith 2009) and develop a new Galerkin approximation. We discuss the relative merit of the three approximations in terms of their mathematical rigor and conservation laws and illustrate their efficiency and caveats by solving linear stability problems with nonzero surface buoyancy.

Using concrete examples, we show that the concerns expressed by earlier authors regarding the suitability of the standard modes {*p*_{n}} are overstated; even with nonzero surface buoyancy, the Galerkin expansion of the streamfunction *ψ* in terms of {*p*_{n}} converges absolutely and uniformly with no Gibbs phenomena, and the same is true for the potential vorticity *q*. A modest number of terms provide a good approximation to *ψ* and *q* throughout the domain, including on the top and bottom boundaries. In other words, the surface streamfunction can be expanded in terms of {*p*_{n}} and, with enough modes, this representation can then be used to accurately calculate the advection of nonzero surface buoyancy. In section 5, we illustrate this procedure by solving the classic Eady problem using the basis {*p*_{n}} for the streamfunction.

## 2. The exact system

In this section, we summarize the basic properties of the QG system. For a detailed derivation, see Pedlosky (1987).

### a. Formulation

The streamfunction is denoted *ψ*(*x*, *y*, *z*, *t*), and we use the following notation:

The variable *ϑ* is related to the buoyancy by *b* = *N*^{2}*ϑ*/*f*_{0}. The QG potential vorticity (QGPV) equation is

where the potential vorticity is

with

Also in (6), the Jacobian is .

The boundary conditions at the top (*z* = *z*^{+}) and bottom (*z* = *z*^{−}) are that *w* = 0, or, equivalently,

Above we have used the superscripts + and − to denote evaluation at *z*^{+} and *z*^{−}, for example, *ψ*^{+} = *ψ*(*x*, *y*, *z*^{+}, *t*).

### b. Quadratic conservation laws

In the absence of sources and sinks, the exact QG system has four quadratic conservation laws: energy, potential enstrophy, and surface buoyancy variance at the two surfaces (e.g., Pedlosky 1987; Vallis 2006). Throughout, we assume horizontal periodic boundary conditions.

The well-known energy conservation law is

The total energy is *ρ*_{0}*E*, where *ρ*_{0} is a reference density.

If *β* = 0 then there are many quadratic potential enstrophy invariants: the volume integral of *q*^{2}*A*(*z*), with *A*(*z*) an arbitrary function of the vertical coordinate, is conserved. The choice reduces to conservation of the surface integral of *q*^{2} at any level .

Charney (1971) noted that, in a doubly periodic domain, nonzero *β* destroys all these quadratic potential enstrophy conservation laws, including the conservation of potential enstrophy defined simply as the volume integral of *q*^{2}. Multiplying the QGPV equation [(6)] by *q*, and integrating by parts, we obtain

The potential enstrophy equation [(11)] is the finite-depth analog of equation (13) in Charney’s paper. To make progress, Charney assumed *ϑ* = 0 at the ground. But the *β* term on the right of (11) can be eliminated by cross-multiplying the QGPV equation [(6)] evaluated at the surfaces with the boundary conditions (9) and combining with (11). Thus, in a doubly periodic domain, nonzero *β* selects a uniquely conserved potential enstrophy from the infinitude of *β* = 0 potential enstrophy conservation laws:

With *β* ≠ 0, the surface contributions in (12) are required to form a conserved quadratic quantity involving *q*^{2}. Notice that *Z* is not sign definite. To our knowledge, the conservation law in (12) is previously unremarked.

Finally, in addition to *E* and *Z*, the surface buoyancy variance is conserved on each surface:

Thus, with *β* ≠ 0, the QG model has four quadratic conservation laws: *E*, *Z*, and the buoyancy variance at the two surfaces.

## 3. Galerkin approximation using standard vertical modes

A straightforward approach is to represent the streamfunction by linearly combining the first + 1 vertical modes. The mean-square error in this approximation is

We use a script font, and context, to distinguish the truncation index in (14) from the buoyancy frequency *N*(*z*). The coefficients *a*_{0} through are determined to minimize err_{ψ}, and thus one obtains the Galerkin approximation to the exact streamfunction:

where the coefficients in the sum above are

Throughout we use the superscript ^{⌣} to denote a Galerkin coefficient defined via projection of a field onto a vertical mode.

In complete analogy with the streamfunction, one can also develop an ( + 1)-mode Galerkin approximation to the PV:

with coefficients

The construction of the Galerkin approximation above minimizes a mean-square error err_{q} defined in analogy with (14).

Now recall that the exact *ψ* and *q* are related by the elliptic “inversion problem”

with boundary conditions at :

The Galerkin approximations in (15) through (18) are defined independently of the information in (19) and (20). The relationship between the Galerkin coefficients and is obtained by multiplying (19) by (1/*h*)*p*_{n}(*z*) and integrating over the depth. Noting the intermediate result

we obtain

where Δ_{n} is the *n*th-mode Helmholtz operator

The relation in (22) is the key to a good Galerkin approximation to surface active quasigeostrophic dynamics.

Term-by-term differentiation of the series in (15) does not give the series in (17) unless *ϑ*^{±} = 0. In other words, term-by-term differentiation does not produce the correct relation [(22)] between and . Thus, the Galerkin-truncated PV and the Galerkin-truncated streamfunction do not satisfy the inversion boundary value problem exactly:

Despite (24), the truncated series and are the best least squares approximations to *ψ* and *q*.

Notice that, in analogy with the Galerkin approximations for *q* and *ψ*,

where

are finite approximations to distributions *δ*(*z* − *z*^{±}) at the surfaces. Of course, these surface *δ* distributions do not satisfy the *L*^{2} convergence condition in (4), and thus the series in (26) only converge in a distributional sense (e.g., Hunter and Nachtergaele 2001). For instance, if *f* satisfies the *L*^{2} convergence condition in (4), then

as → ∞. Thus, in that limit,

where denotes distributional convergence. The right-hand side of (28) is the Brethertonian modified potential vorticity (Bretherton 1966) with the boundary conditions incorporated as PV sheets. To illustrate (24) and (28), we present an elementary example that is relevant to our discussion of the Eady problem in section 5.

### An elementary example: The Eady basic state

As an example, consider the case with constant buoyancy frequency *N*. We use nondimensional units so that the surfaces are at *z*^{−} = −1 and *z*^{+} = 0. The standard vertical modes are *p*_{0} = 1 and, for *n* ≥ 1,

with *κ*_{n} = *nπ*.

We consider the basic state of the Eady problem with the streamfunction

and zero interior PV *q* = 0 and *β* = 0. The surface buoyancies are *ϑ*^{±} = −*y*.

The Galerkin expansion of the PV *q* = 0 is exact: , and therefore . The truncated Galerkin expansion of *ψ* follows from either (16) or (22) and is

(We assume that is odd so that the last term in the truncated series is as above.) Despite the nonzero derivative of *ψ* at the boundaries, the series in (31) is absolutely and uniformly convergent on the closed interval −1 ≤ *z* ≤ 0. The ^{−2} behavior of the series (31) ensures uniform convergence, for example, using the *M* test (Hunter and Nachtergaele 2001). There are no Gibbs oscillations and a modest number of terms provide a good approximation to the base velocity *U* (Fig. 1a). All these desirable properties are lost if we differentiate (31) with respect to *z*.

Thus, to illustrate (24) and (28), notice that if one attempts to calculate the Eady PV, namely, *q* = 0, by direct differentiation of (31), one obtains

The series (32) does not converge in a pointwise sense, and the partial sum is violently oscillatory as → ∞. However, in a distributional sense (Hunter and Nachtergaele 2001, chapter 11), the exact sum in (33) does converge to *δ* distributions on the boundaries (see Figs. 1b,c). These boundary *δ* distributions are the Brethertonian PV sheets (Bretherton 1966). To some extent, which we investigate in section 5a, the series (32) is rescued by this Brethertonian interpretation.

Of course the correct Galerkin approximation to the Eady PV *q* = 0 is the splendidly convergent series which is obtained if one uses either (18) or (22) to obtain . This seemingly trivial example illustrates potentially confusing issues that arise if one differentiates a Galerkin approximation; the standard modes provide good representations of *ψ* and *q*, even if *ψ*_{z} is nonzero on the boundaries. The problem is that differentiating the *ψ* series to produce a *q* series does not produce the Galerkin approximation to *q*.

## 4. Three approximations

In (24) we noted that the Galerkin approximations to *ψ* and *q* do not exactly satisfy the inversion relation. To address this error there are at least three different approximations one can make. The three approximations are equivalent when *ϑ*^{±} = 0. In the next three subsections, we provide a detailed description of each approximation. After testing, we recommend approximation *C* as the most reliable approximation using standard vertical modes.

### a. Approximation A

Approximation *A* uses the truncated series in (15) as a least squares Galerkin approximation to the streamfunction *ψ*. Approximation *A* does not use the Galerkin approximation for *q*. Instead, the approximate PV is defined so that the interior inversion relation is satisfied exactly:

This is the approximation introduced by Flierl (1978) in a context without surface buoyancy, and it is now regarded as the standard in physical oceanography. Note that in (34) is not the least squares approximation to the exact *q*. Instead the series is obtained using term-by-term differentiation of the series . The example surrounding (32) shows that with nonzero surface buoyancy, the approximation is strongly oscillatory in the interior of the domain and approaches the Brethertonian PV on the right of (28) as → ∞. The rapid interior oscillation of is a spurious creation of term-by-term differentiation. Later, in section 5, these spurious oscillations will produce significant errors in the solution of the Eady stability problem.

Following Flierl (1978), in approximation *A* the -mode approximate PV is defined via (34) and, using the modal representation for in (15), this is equivalent to

where Δ_{n} is the Helmholtz operator in (23). Following the appendix of Flierl (1978), one can use the Galerkin projection of the nonlinear evolution equation [(6)] on the modes _{n} to obtain + 1 evolution equations for the coefficients :

where

Note that Ξ_{nms} cannot be computed exactly except in cases with simple buoyancy frequency profiles. But it suffices to compute Ξ_{nms} to high accuracy, for example, using Gaussian quadrature.

Flierl (1978) implicitly assumed that *ϑ*^{+} = *ϑ*^{−} = 0 so that the surface terms in (22) vanish and then there is no difference between and . But in general, with nonzero surface buoyancy, we can append evolution equations for *ϑ*^{+} and *ϑ*^{−} to approximation *A*. That is, in addition to the + 1 modal equations in (36), we also have

Above we have evaluated the *ψ* series (15) at *z*^{±} to approximate *ψ*^{±} in the surface boundary conditions. This approach is not satisfactory because the resulting surface buoyancy equation [(38)] is dynamically passive, that is, *ϑ*^{+} and *ϑ*^{−} do not affect the interior evolution equations in (36).

Approximation *A* has the well-known energy conservation law

To obtain the energy analogous to *E* in (10), the modal sum above is multiplied by the depth *h*. Approximation *A* also has the potential enstrophy conservation law

### b. Approximation B

Approximation *B* begins with the observation that the exact solution of the inversion problem in (19) and (20) can be decomposed as

where *ψ*^{I}(*x*, *y*, *z*, *t*) is the “interior streamfunction” and *ψ*^{S}(*x*, *y*, *z*, *t*) is the “surface streamfunction” (Lapeyre and Klein 2006; Tulloch and Smith 2009).

The surface streamfunction *ψ*^{S}(*x*, *y*, *z*, *t*) is defined as the solution of the boundary value problem

with inhomogeneous Neumann boundary conditions

The interior streamfunction *ψ*^{I}(*x*, *y*, *z*, *t*) is defined as the solution of the boundary value problem

with homogeneous Neumann boundary conditions

Approximation *B* assumes that one can solve the surface problem in (42) and (43) without resorting to a truncated series. For instance, with constant or exponential stratifications, one can find closed-form, exact expressions for *ψ*^{S} (Tulloch and Smith 2009; LaCasce 2012). Approximation *B* requires that the two unknown Dirichlet boundary condition functions *ψ*^{S±} = *ψ*^{S}(*z*^{±}) can be obtained efficiently from specified Neumann boundary condition functions *ϑ*^{+} and *ϑ*^{−}. The Eady problem, discussed below in section 5, is a prime example in which one can obtain this Neumann-to-Dirichlet map.

Once *ψ*^{S} is in hand, the approximate streamfunction is

where the approximate interior streamfunction is obtained by solving the interior inversion problem (44) with the right-hand side replaced by the -mode Galerkin approximation defined in (17) and (18). The two-mode two-surface model of Tulloch and Smith (2009) is the case = 1. The exact solution of the approximate interior inversion problem is

where

To obtain the approximation *B* evolution equations, we introduce the streamfunction (46) into the QGPV equation [(6)] and project onto mode *n* to obtain

with Ξ_{nms} defined in (37). Approximation B assumes that the remaining integral of (49) can be evaluated exactly. This is only possible for particular models of the *N*(*z*) (e.g., constant buoyancy frequency profiles). In practice, however, it may suffice to compute the integral very accurately, for example, using Gaussian quadrature.

The evolution equations for approximation *B* are completed with the addition of buoyancy advection at the surfaces:

Approximation *B* conserves surface buoyancy variance. But the conservation laws for energy are problematic; because *ψ*^{I} is not orthogonal to *ψ*^{S} the energy, (10) is not conserved in approximation B (K. S. Smith 2014, personal communication). These nonconservative effects are quantitatively small but are nonetheless irritating. The nonorthogonality of *ψ*^{I} and *ψ*^{S} was a motivation for development of the “surface-aware” modes by Smith and Vanneste (2013).

With *β* = 0, approximation B conserves potential enstrophy:

But with *β* ≠ 0, the analog of the exact potential enstrophy (12) is not conserved.

### c. Approximation C

Approximation *C* uses truncated Galerkin approximations and for both *ψ* and *q*. The series is not obtained by differentiation of , and therefore, as indicated in (24), the inversion equation is not satisfied exactly by and . But instead, one will have true least squares approximations to both *ψ* and *q*. To our knowledge, approximation *C*, correctly accounting for the surface buoyancy boundary terms, has not been previously investigated.

Because method *C* approximates both the streamfunction and the PV by Galerkin series, the derivation of the modal equations is very straightforward compared with the calculations in appendix A of Flierl (1978); one simply substitutes the truncated Galerkin series for the streamfunction (15) and PV (17) into the QGPV equation [(6)] and then projects onto mode *n* to obtain

In approximation *C*, there are + 3 degrees of freedom: the + 3 modal amplitudes and the two surface buoyancy fields *ϑ*^{±}. The approximation *C* evolution equations are completed by advection of the surface buoyancy:

We emphasize that in approximation *C* the surface buoyancy fields *ϑ*^{±} are not passive; , , and *ϑ*^{±} are related through (53).

Approximation *C* conserves surface buoyancy variance as in (13). Total energy is also conserved:

With *β* = 0, approximation *C* has a potential enstrophy conservation law

But with *β* ≠ 0, as in *B*, approximation *C* does not conserve the analog of the exact potential enstrophy [(12)].

## 5. The Eady problem

We use classical linear stability problems with nonzero surface buoyancy to illustrate how solutions to specific problems can be constructed and to assess the relative merit and efficiency of approximations *A*, *B*, and *C*. The linear analysis does not provide the full picture of convergence of the approximate solutions. Nonetheless, in turbulence simulations forced by baroclinic instability, it is necessary (but not sufficient) to accurately capture the linear stability properties.

We use nondimensional variables so that the surfaces are at *z*^{+} = 0 and *z*^{−} = −1. The Eady exact base-state velocity is given by (30) with zero PV *q* = 0 and *β* = 0.

### a. Approximation A

While the surface fields *ϑ*^{±} are dynamically passive in approximation A, the Eady problem can still be considered because the base-state PV defined via (35) converges to *δ* distributions on the boundaries (section 3).

The base-state velocity in approximation *A* is given by the series (31) and is a good approximation to the exact base-state velocity of (30). But, according to approximation *A*, there is a nonzero interior base-state PV gradient given by the series (33). As → ∞, the PV gradient in (33) converges in a distributional sense to Brethertonian sheets at *z* = 0 and −1. But for numerical implementation of approximation *A*, we stop short of = ∞. While the PV gradient is much larger at the boundaries, there is always interior structure in the PV (Fig. 1c). We show that this spurious interior PV gradient has a strong and unpleasant effect on the approximate solution of the Eady stability problem.

To solve the Eady linear stability we linearize the interior equations in (36) about the base-state velocity in (31) and the PV gradient in (33). We assume , and so on, to obtain a ( + 1) × ( + 1) eigenproblem:

where are the coefficients of the series (33), and . The eigenproblem (57) can be recast in the matrix form , where ( appendix B) and solved with standard methods.

Figure 2 shows the growth rate of the Eady instability according to approximation *A* and compares this with the exact Eady growth rate. Approximation *A* does not do well, especially at large wavenumbers. The exact Eady growth rate has a high-wavenumber cutoff. At moderate values of , such as 3, 5, and 7, approximation *A* produces unstable “bubbles” of instability at wavenumbers greater than the high-wavenumber cutoff. The growth rates in these bubbles are comparable to the true maximum growth rate. As increases, the unstable bubbles are replaced by a long tail of unstable modes with a growth rate that slowly increases with *κ*. These spurious high-wavenumber instabilities are due to the rapidly oscillatory interior PV gradient that supports unphysical critical layers (see Fig. 3).

### b. Approximation B: The exact solution

In approximation *B*, the zero PV in the Eady problem implies . The + 1 modal equations (with *β = 0*) are trivially satisfied; there is no interior contribution . Thus, approximation *B* solves the Eady problem exactly.

where the magnitude of the wavenumber vector is . We evaluate the surface streamfunction (58) at the boundaries to find the relationship between the streamfunction at the surfaces and the boundary fields *ϑ*^{±}:

The nondimensional linearized boundary conditions (50) are

The eigenvalues of are given by the celebrated dispersion relation for the Eady problem (Pedlosky 1987; Vallis 2006):

### c. Approximation C

Approximation *C* expands both the streamfunction and the PV in standard vertical modes. Thus, in the Eady problem the PV is exactly zero, as it should be: . Thus, approximation *C* does not have the spurious critical layers that bedevil *A*. Moreover, in approximation *C*, the + 1 modal equations (with *β* = 0) in (52) are trivially satisfied, and the inversion relationship (53) provides a simple connection between the streamfunction and the fields *ϑ*^{±}. The base velocity for the Eady problem in approximation *C* is the series in (31) (the same as *A*). From the exact shear at the boundaries we obtain the exact base-state boundary variables:

Assuming , and using the inversion relationship (53), we obtain a 2 × 2 eigenproblem:

where the matrix is defined in appendix C. It is straightforward to show that *c*^{C} converges to the exact eigenspeed, that is, *c*^{C} → *c*^{B} as → ∞ ( appendix C). Figure 2 shows that approximation *C* successfully captures the structure of the Eady growth rate even with modest values of .

### d. Remarks on convergence

The crudest truncation (i.e., = 0) is stable for both approximations *A* and *C* (Fig. 2). With one baroclinic mode ( = 1) the growth rates are qualitatively consistent with the exact solution, and the results improve with = 2. With a moderate number of baroclinic modes ( > 2), approximations *A* and *C* converge rapidly to the exact growth rate at wavenumbers less than about 2.2 (see Fig. 2). But surprisingly the convergence of the growth rate at the most unstable mode (*κ* ≈ 1.6) is faster in approximation *A* (~^{−4}) than in approximation *C* (~^{−2}) (see Fig. 4). However, the convergence in approximation *C* is uniform; there are no spurious high-wavenumber instabilities.

Figure 4 also shows that the approximation *A* convergence of the growth rate to zero at *κ* = 8 is slow (~^{−1}). While the growth rate does converge to zero at a fixed wavenumber, such as *κ* = 8, we conjecture that there are always faster-growing modes at larger wavenumbers.

## 6. The Green problem

To further explore the relative merit and efficiency of approximations *A*, *B*, and *C*, we study the instability properties of a system with nonzero *β*. For simplicity, we consider a problem with Eady’s base state *ψ* = −(1 + *z*)*y* on a *β* plane. This is similar to the problem originally considered by Charney (1947) and Green (1960). The major difference is that Charney considered a vertically semi-infinite domain (Charney 1947; Pedlosky 1987), while we follow Green and consider a finite-depth domain with −1 < *z* < 0.

We obtain the exact system for this “Green problem” by linearizing the QG equations [(6)–(9)] about the base state (30) with background PV *βy*, where is the nondimensional planetary PV gradient. Assuming , we obtain

and

As a reference solution, we solve the eigenproblems (66)–(67) using a centered, second-order, finite-difference scheme with 1000 vertical levels (see Fig. 5).

The Green problem supports three classes of unstable modes, indicated in the lower-right panel ( = 128) of Fig. 5: 1) the “modified Eady modes,” which are instabilities that arise from the interaction of Eady-like edge waves, only slightly modified by *β*; 2) the “Green modes,” which are very long slowly growing modes (Vallis 2006); and 3) the high-wavenumber “Charney modes,” which are critical layer instabilities that arise from the interaction of the surface edge wave with the interior Rossby wave that is supported by nonzero *β*.

### a. Implementation of approximation A

The base state for the Green problem is the same as in the Eady problem. In approximation *A*, the *β* term adds only a diagonal term to the Eady system [(57)] (see appendix C).

### b. Implementation of approximation B

The base state is the same as in the Eady problem. The steady streamfunction and buoyancy fields that satisfy (49) and (50) exactly are

where

and

where is given by (58). We use the inversion relationship [(48)] and the Neumann-to-Dirichlet map [(59)] to recast this eigenproblem into standard form , where (see appendix C).

### c. Implementation of approximation C

Again the base state is the same as in the Eady problem. But now there are + 3 equations: the two boundary equations of Eady’s problem [(64)] plus + 1 interior equations:

We use the inversion relationship of (53) in (73) to recast this eigenproblem in the form , where is defined as in approximation *B* (see appendix C).

### d. Remarks on convergence

The most crude truncation ( = 0) is stable for approximations *A* and *C*. In contrast, the = 0 truncation in approximation *B* is qualitatively consistent with the modified Eady instabilities (see Fig. 5). With a moderate number of baroclinic modes ( = 2 or 3), approximations *A*, *B*, and *C* all resolve the modified Eady modes relatively well. At the most unstable modified Eady mode (*κ* ≈ 1.9), approximation *B* has typically the smallest error because it solves the surface problem exactly. As in the Eady problem, approximation *A* converges (~^{−4}) faster than approximations *B* and *C* (~^{−2}) at the most unstable mode, but *B* and *C* converge faster at high wavenumbers.

Approximations *A*, *B*, and *C* all converge very slowly to the high-wavenumber Charney modes (Figs. 5, 6). These modes are interior critical-layer instabilities (Pedlosky 1987), and the critical layer is confined to a small region about the steering level (i.e., the depth at which the phase speed matches the base velocity; see Fig. 7). With finite base-state shear, the critical layer is always in the interior. Thus, the problem is not that standard vertical modes are inefficient because they do not satisfy inhomogeneous boundary conditions; a low-resolution finite-difference solution also presents such bubbles in high-wavenumber growth rates (not shown). Resolution of the interior critical layer, not the surface boundary condition, is a problem for all methods at high wavenumbers. The surface-aware modes of Smith and Vanneste (2013) have similar performance to approximations *B* and *C*, but also have the same limitation—a large number of vertical modes are required to resolve interior critical layers (K. S. Smith 2015, personal communication).

For example, with < 25, at *κ* = 8, approximations are qualitatively inconsistent with the high-resolution finite-difference solution. For larger values of , the growth rate convergence for approximations *B* and *C* scales ~^{−3}. The growth rate for approximation *A* converges painfully slowly (~^{−1}). As in the Eady problem, at large wavenumbers, the growth rate for approximation *A* is qualitatively different from that of the finite-difference solution because of spurious instabilities associated with the rapidly oscillatory base-state PV gradient in (33).

## 7. Summary and conclusions

The Galerkin approximations *A*, *B*, and *C* are equivalent if there are no buoyancy variations at the surfaces. Thus, all three approximations are well suited for applications with zero surface buoyancy (Flierl 1978; Fu and Flierl 1980; Hua and Haidvogel 1986). But with nonzero surface buoyancy, the three approximations are fundamentally different. In particular, approximation *A*, originally introduced by Flierl (1978) in a context without surface buoyancy, obtains the approximate PV by differentiating the Galerkin series for the streamfunction, and consequently its approximate PV has violent oscillations in the interior. Approximation *B* represents the PV as a Galerkin series in standard modes and calculates the streamfunction that satisfies the exact inversion problem associated with the approximate PV (Tulloch and Smith 2009). The inversion relationship is split into surface and interior problems. Because the surface streamfunction projects onto the interior solution, the energy is not diagonalized, and consequently approximation *B* has small errors in energy conservation (K. S. Smith 2014, personal communication). The surface-aware modes of Smith and Vanneste (2013) correct this problem. Approximation *C* uses the Galerkin series for both streamfunction and PV but does not satisfy the inversion problem exactly. Nevertheless, the Galerkin series for *ψ* and *q* converge absolutely and uniformly, and approximation *C* provides a good finite truncation of the QG equations that represents surface buoyancy dynamics and also conserves energy.

With nonzero interior PV gradients, the convergence of all approximations is slow for the high-wavenumber Charney-type modes. The critical layer associated with these modes spans a very small fraction of the total depth (Fig. 7). To accurately resolve these near singularities at the steering level, there is no better solution than having high vertical resolution in the interior.

For problems with nonuniform surface buoyancy and nonzero interior PV gradient, we recommend approximation C for obtaining solutions to the three-dimensional QG equations using standard vertical modes.

The codes that produced the numerical results of this paper, plotting scripts, and supplementary figures are openly available online (at https://github.com/crocha700/qg_vertical_modes).

## Acknowledgments

CBR is grateful for a helpful conversation with G. R. Flierl. We thank Joe LaCasce and Shafer Smith for reviewing this paper. We had useful discussions with Shafer Smith, regarding approximation *B* and the surface-aware modes, and with Geoff Vallis, who pointed out that Green (1960) first considered the Eady+*β* problem. This research was supported by the National Science Foundation under Award OCE 1357047.

### APPENDIX A

#### Convergence of Galerkin Series in Standard Modes

Jackson (1914) gives conditions for the uniform convergence of series expansions in eigenfunctions of the Sturm–Liouville eigenproblem

defined on the interval *Z* ∈ [0, *π*] with boundary conditions

where *γ*_{0} and *γ*_{π} are real constants of arbitrary sign, and is the eigenvalue. The equations defining the standard modes (1)–(2) can be brought to this form using the following Liouville transformation:

and

The eigenvalues are related by and

The boundary condition for the standard modes [(2)] implies that the transformed modes satisfy (A2) with

If *dS*/*dz* = 0 at a boundary then the appropriate condition at that boundary is *P*_{n} = 0.

A special case of Theorem I from Jackson (1914) states that the expansion of a function *f*(*z*) as a series in eigenfunctions *P*_{n} converges absolutely and uniformly, provided that both *df*/*dz* and *d*Λ/*dz* are continuous and bounded, regardless of whether or not *f* satisfies the same boundary conditions as *P*_{n}. (The remainder of the theorem concerns the rate of convergence under stronger conditions on *ψ* and Λ.) The streamfunction, potential vorticity, and buoyancy profiles are typically assumed to be smooth in studies of QG dynamics, which implies that both *f* and Λ will satisfy the above conditions. Uniform convergence over *Z* ∈ [0, *π*] implies uniform convergence over *z* ∈ [*z*^{−}, *z*^{+}].

### APPENDIX B

#### Derivation of Conservation Laws for Approximation *C*

To obtain the conservation of energy in approximation *C*, we multiply the modal equations in (52) by , integrate over the horizontal surface, and sum on *n* to obtain

The triple sum term vanishes because is fully symmetric and the Jacobian is anti-symmetric. The term on the second line of (B1) is also zero: multiply the boundary conditions of (54) by and integrate over the horizontal surface. Thus, we obtain the energy conservation law in (55).

The analog of the exact potential enstrophy [(12)],

is only conserved with = 0.

### APPENDIX C

#### Details of the Stability Problems

##### a. The interaction tensor

Because the standard vertical modes with constant stratification are simple sinusoids [(29)], the interaction coefficients [(37)] can be computed analytically. First, we recall that Ξ_{ijk} is fully symmetric. Permuting the indices so that *i* ≥ *j* ≥ *k*, we obtain

The second line in (C1) corrects a factor of ½ missed by Hua and Haidvogel (1986).

##### b. Approximation A

Using the symmetry in Ξ_{nms}, and the inversion relation [(35)], we rewrite row *n* + 1 of the linear Green system:

where the inverse of the *n*th-mode Helmholtz operator in Fourier space is

The Eady problem is the special case . We use a standard eigenvalue–eigenvector algorithm to obtain the approximate eigenspeed *c*^{A}.

##### c. Approximation B

The Green eigenvalue problem in (69) through (72) can be recast in the standard form , where . The first and last rows of the system stem from the boundary conditions (71)–(72):

and

The (*n* + 1)th row originates from the *n*th interior equation [(69)]:

where the symmetric matrix *γ*_{ms} is

##### d. Approximation C

###### 1) The Eady problem

The 2 × 2 eigenproblem is

where

The sums (C9) become exact in the limit → ∞:

The base velocity also converges to the exact result. Using standard results for the summation of inverse squares, we obtain

Thus,

and the eigenvalues of the Eady problem using approximation *C* become exact, that is, as → ∞.

###### 2) The Green problem

The ( + 3) × ( + 3) eigenproblem is

where is defined as above in approximation *B*. The first and last rows of (C13) stem from the boundary conditions (64):

and

Row *n* + 1 originates from the *n*th modal equation [(73)]:

## REFERENCES

*Atmosphere–Ocean Dynamics.*Academic Press, 662 pp.

*Applied Analysis*. World Scientific, 456 pp.

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

*Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation*. Cambridge University Press, 769 pp.

## Footnotes

Current affiliation: Department of Applied Mathematics, University of Colorado, Boulder, Colorado.