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 pn(z), are defined by the Sturm–Liouville eigenproblem

 
formula

with homogeneous Neumann boundary conditions at the bottom (z = z) and top (z = z+) surfaces of the domain:

 
formula

In (1), N(z) is the buoyancy frequency and f0 is the Coriolis parameter. The eigenvalue κn in (1) is the deformation wavenumber of the nth mode. With normalization, the modes satisfy the orthogonality condition

 
formula

where is the depth. The barotropic mode is p0 = 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 {pn} (Gill 1982; Pedlosky 1987; Vallis 2006; Ferrari and Wunsch 2010; LaCasce 2012). In fact, the set {pn} is mathematically complete and can be used to represent any field with a finite square integral:

 
formula

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 {pn} converges in L2(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 {pn}.

Despite the rigorous assurance of completeness in the previous paragraph, the utility of {pn} 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 {pn} 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 {pn}. We show that that both the inversion problem and evolutionary dynamics can be handled using {pn} 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 {pn} are overstated; even with nonzero surface buoyancy, the Galerkin expansion of the streamfunction ψ in terms of {pn} 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 {pn} 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 {pn} 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:

 
formula

The variable ϑ is related to the buoyancy by b = N2ϑ/f0. The QG potential vorticity (QGPV) equation is

 
formula

where the potential vorticity is

 
formula

with

 
formula

Also in (6), the Jacobian is .

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

 
formula

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

 
formula

The total energy is ρ0E, where ρ0 is a reference density.

If β = 0 then there are many quadratic potential enstrophy invariants: the volume integral of q2A(z), with A(z) an arbitrary function of the vertical coordinate, is conserved. The choice reduces to conservation of the surface integral of q2 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 q2. Multiplying the QGPV equation [(6)] by q, and integrating by parts, we obtain

 
formula

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:

 
formula

With β ≠ 0, the surface contributions in (12) are required to form a conserved quadratic quantity involving q2. 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:

 
formula

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

 
formula

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

 
formula

where the coefficients in the sum above are

 
formula

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:

 
formula

with coefficients

 
formula

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

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

 
formula

with boundary conditions at :

 
formula

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)pn(z) and integrating over the depth. Noting the intermediate result

 
formula

we obtain

 
formula

where Δn is the nth-mode Helmholtz operator

 
formula

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:

 
formula

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 ψ,

 
formula

where

 
formula

are finite approximations to distributions δ(zz±) at the surfaces. Of course, these surface δ distributions do not satisfy the L2 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 L2 convergence condition in (4), then

 
formula

as → ∞. Thus, in that limit,

 
formula

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 p0 = 1 and, for n ≥ 1,

 
formula

with κn = .

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

 
formula

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

 
formula

(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.

Fig. 1.

Nondimensional base state for the Eady problem using various truncations for the series (31). In (b), is the number of baroclinic modes. (a) Zonal velocity, although the truncation has zero slope at the boundaries there are no Gibbs oscillations. (b) Meridional PV gradient associated with the truncated series (33). (c) As in (b), but with an expanded abscissa. As increases, the PV gradient distributionally converges to two Brethertonian delta functions at the boundaries. The insets represent close-ups of the figures in regions enclosed by gray rectangles.

Fig. 1.

Nondimensional base state for the Eady problem using various truncations for the series (31). In (b), is the number of baroclinic modes. (a) Zonal velocity, although the truncation has zero slope at the boundaries there are no Gibbs oscillations. (b) Meridional PV gradient associated with the truncated series (33). (c) As in (b), but with an expanded abscissa. As increases, the PV gradient distributionally converges to two Brethertonian delta functions at the boundaries. The insets represent close-ups of the figures in regions enclosed by gray rectangles.

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

 
formula
 
formula

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:

 
formula

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

 
formula

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 :

 
formula

where

 
formula

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

 
formula

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

 
formula

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

 
formula

But the analog of the exact potential enstrophy [(12)] is not conserved by A (nor by B and C below). Finally, with the surface equations in (38), approximation A also conserves surface buoyancy variance as in (13).

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

 
formula

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

 
formula

with inhomogeneous Neumann boundary conditions

 
formula

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

 
formula

with homogeneous Neumann boundary conditions

 
formula

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

 
formula

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

 
formula

where

 
formula

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

 
formula

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:

 
formula

With (49) and (50), we have + 3 evolution equations for the + 3 fields and .

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:

 
formula

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

 
formula

where Ξnms is defined in (37), and we recall the relation between and from (22):

 
formula

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:

 
formula

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:

 
formula

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

 
formula

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:

 
formula

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).

Fig. 2.

Growth rate for the Eady problem as a function of the zonal wavenumber (l = 0) using approximations A, B (exact), and C with various numbers of baroclinic modes .

Fig. 2.

Growth rate for the Eady problem as a function of the zonal wavenumber (l = 0) using approximations A, B (exact), and C with various numbers of baroclinic modes .

Fig. 3.

Structure of κ = 8 unstable mode for the Eady problem obtained using approximation A and = 64. Streamfunction is the black curves, and PV is the colors. The streamfunction slightly tilts westward as z increases. One can see the unphysical critical layer associated with the fast-oscillating base-state PV. The critical level zc is the depth where the unstable wave speed matches the velocity of the base state. Only the top quarter of the domain is shown.

Fig. 3.

Structure of κ = 8 unstable mode for the Eady problem obtained using approximation A and = 64. Streamfunction is the black curves, and PV is the colors. The streamfunction slightly tilts westward as z increases. One can see the unphysical critical layer associated with the fast-oscillating base-state PV. The critical level zc is the depth where the unstable wave speed matches the velocity of the base state. Only the top quarter of the domain is shown.

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.

Assuming , we obtain the solution to the surface streamfunction inversion problems (42)(43):

 
formula

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 ϑ±:

 
formula

The nondimensional linearized boundary conditions (50) are

 
formula

where cB = ωB/k. Using the boundary conditions (60) in (59), we obtain an eigenvalue problem

 
formula

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

 
formula

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:

 
formula

We linearize the boundary equations [(54)] about the base states (33) and (63) to obtain

 
formula

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

 
formula

where the matrix is defined in  appendix C. It is straightforward to show that cC converges to the exact eigenspeed, that is, cCcB 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.

Fig. 4.

Absolute error as a function of number of baroclinic modes for the growth rates of the Eady problem. The solid lines show the error at the exact fastest-growing mode (κ ≈ 1.6). The dashed line is the approximation A error at κ = 8.

Fig. 4.

Absolute error as a function of number of baroclinic modes for the growth rates of the Eady problem. The solid lines show the error at the exact fastest-growing mode (κ ≈ 1.6). The dashed line is the approximation A error at κ = 8.

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

 
formula

and

 
formula

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).

Fig. 5.

Growth rate for the Green problem with as a function of the zonal wavenumber (l = 0) using approximations A, B, and C with various numbers of baroclinic modes . The black line is a finite-difference solution with 1000 vertical levels.

Fig. 5.

Growth rate for the Green problem with as a function of the zonal wavenumber (l = 0) using approximations A, B, and C with various numbers of baroclinic modes . The black line is a finite-difference solution with 1000 vertical levels.

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

 
formula

Assuming , the + 1 interior equations [(49)] linearized about (68) are

 
formula

where

 
formula

The boundary conditions of (50), linearized about (68), are

 
formula

and

 
formula

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:

 
formula

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).

Fig. 6.

Absolute error as a function of number of baroclinic modes for the growth rates of the Green problem. The solid line represents the error at the exact fastest-growing mode (κ ≈ 1.9). The dashed line is the error at κ = 8.

Fig. 6.

Absolute error as a function of number of baroclinic modes for the growth rates of the Green problem. The solid line represents the error at the exact fastest-growing mode (κ ≈ 1.9). The dashed line is the error at κ = 8.

Fig. 7.

Wave structure of the κ = 8 unstable mode for the Green problem with solved using a second-order finite-difference scheme with 1000 vertical levels. Streamfunction (black contours) and potential vorticity (colors). The streamfunction slightly tilts westward as z increases. The potential vorticity is confined to a small region, the critical layer. The critical level zc is the depth where the unstable wave speed matches the velocity of the base state. Only the bottom quarter of the domain is shown.

Fig. 7.

Wave structure of the κ = 8 unstable mode for the Green problem with solved using a second-order finite-difference scheme with 1000 vertical levels. Streamfunction (black contours) and potential vorticity (colors). The streamfunction slightly tilts westward as z increases. The potential vorticity is confined to a small region, the critical layer. The critical level zc is the depth where the unstable wave speed matches the velocity of the base state. Only the bottom quarter of the domain is shown.

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

 
formula

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

 
formula

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:

 
formula

and

 
formula

The eigenvalues are related by and

 
formula

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

 
formula

If dS/dz = 0 at a boundary then the appropriate condition at that boundary is Pn = 0.

A special case of Theorem I from Jackson (1914) states that the expansion of a function f(z) as a series in eigenfunctions Pn 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 Pn. (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

 
formula

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)],

 
formula

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 ijk, we obtain

 
formula

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:

 
formula

where the inverse of the nth-mode Helmholtz operator in Fourier space is

 
formula

The Eady problem is the special case . We use a standard eigenvalue–eigenvector algorithm to obtain the approximate eigenspeed cA.

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):

 
formula

and

 
formula

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

 
formula

where the symmetric matrix γms is

 
formula
d. Approximation C
1) The Eady problem

The 2 × 2 eigenproblem is

 
formula

where

 
formula

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

 
formula

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

 
formula

Thus,

 
formula

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

2) The Green problem

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

 
formula

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

 
formula

and

 
formula

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

 
formula

REFERENCES

REFERENCES
Bretherton
,
F.
,
1966
:
Critical layer instability in baroclinic flows
.
Quart. J. Roy. Meteor. Soc.
,
92
,
325
334
, doi:.
Callies
,
J.
,
G.
Flierl
,
R.
Ferrari
, and
B.
Fox-Kemper
,
2015
:
The role of mixed layer instabilities in submesoscale turbulence
.
J. Fluid Mech.
, in press.
Charney
,
J. G.
,
1947
:
The dynamics of long waves in a baroclinic westerly current
.
J. Meteor.
,
4
,
136
162
, doi:.
Charney
,
J. G.
,
1971
:
Geostrophic turbulence
.
J. Atmos. Sci.
,
28
,
1087
1095
, doi:.
Ferrari
,
R.
, and
C.
Wunsch
,
2010
:
The distribution of eddy kinetic and potential energies in the global ocean
.
Tellus
,
62A
,
92
108
, doi:.
Flierl
,
G. R.
,
1978
:
Models of vertical structure and the calibration of two-layer models
.
Dyn. Atmos. Oceans
,
2
,
341
381
, doi:.
Fu
,
L.-L.
, and
G. R.
Flierl
,
1980
:
Nonlinear energy and enstrophy transfers in a realistically stratified ocean
.
Dyn. Atmos. Oceans
,
4
,
219
246
, doi:.
Gill
,
A. E.
,
1982
: Atmosphere–Ocean Dynamics. Academic Press, 662 pp.
Green
,
J.
,
1960
:
A problem in baroclinic stability
.
Quart. J. Roy. Meteor. Soc.
,
86
,
237
251
, doi:.
Hua
,
B.
, and
D.
Haidvogel
,
1986
:
Numerical simulations of the vertical structure of quasi-geostrophic turbulence
.
J. Atmos. Sci.
,
43
,
2923
2936
, doi:.
Hunter
,
J. K.
, and
B.
Nachtergaele
,
2001
: Applied Analysis. World Scientific, 456 pp.
Jackson
,
D.
,
1914
:
On the degree of convergence of Sturm-Liouville series
.
Trans. Amer. Math. Soc.
,
15
,
439
466
, doi:.
LaCasce
,
J.
,
2012
:
Surface quasigeostrophic solutions and baroclinic modes with exponential stratification
.
J. Phys. Oceanogr.
,
42
,
569
580
, doi:.
Lapeyre
,
G.
,
2009
:
What vertical mode does the altimeter reflect? On the decomposition in baroclinic modes and on a surface-trapped mode
.
J. Phys. Oceanogr.
,
39
,
2857
2874
, doi:.
Lapeyre
,
G.
, and
P.
Klein
,
2006
:
Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory
.
J. Phys. Oceanogr.
,
36
,
165
176
, doi:.
Pedlosky
,
J.
,
1987
: Geophysical Fluid Dynamics. 2nd ed. Springer-Verlag, 710 pp.
Roullet
,
G.
,
J.
McWilliams
,
X.
Capet
, and
M.
Molemaker
,
2012
:
Properties of steady geostrophic turbulence with isopycnal outcropping
.
J. Phys. Oceanogr.
,
42
,
18
38
, doi:.
Smith
,
K. S.
, and
J.
Vanneste
,
2013
:
A surface-aware projection basis for quasigeostrophic flow
.
J. Phys. Oceanogr.
,
43
,
548
562
, doi:.
Tulloch
,
R.
, and
K. S.
Smith
,
2009
:
Quasigeostrophic turbulence with explicit surface dynamics: Application to the atmospheric energy spectrum
.
J. Atmos. Sci.
,
66
,
450
467
, doi:.
Vallis
,
G. K.
,
2006
: 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.