Limits on the Rate of Conversion of Potential to Kinetic Energy in Quasigeostrophic Turbulence

Abstract

Flow configurations that maximize the instantaneous rate of conversion from potential to kinetic energy are sought using a combination of analytical and numerical methods. A hydrostatic model is briefly investigated, but the presence of unrealistic ageostrophic flow configurations renders the results unrealistic. In the quasigeostrophic (QG) model, flow configurations that locally optimize the conversion rate are found, but it remains unclear if these flow configurations produce the global maximum conversion rate. The difficulty is associated with the fact that in the QG model, the vertical velocity is a quadratic function of the QG streamfunction, which renders the conversion rate a cubic function of the QG streamfunction. For these locally maximal conversion rates, the rate of conversion depends on the horizontal length scale of the flow: for scales larger than the deformation radius, the maximal rates are small and decrease as the horizontal scale increases; for scales smaller than the deformation radius, the maximal conversion rate rises until it becomes comparable to the maximal rate at which potential energy can be extracted from the mean flow.

1. Introduction

The transfer of large-scale available potential energy to mesoscale eddy energy by baroclinic instability is one of the most important energy pathways in the dynamics of the global oceans [1]. In the Lorenz energy cycle [2], nonlinear baroclinic instability processes transfer large-scale potential energy to eddy potential energy (EPE); EPE is then converted to eddy kinetic energy (EKE); and EKE is either returned to the large-scale flow or dissipated through a wide range of mechanisms. Theories of geophysical macro-turbulence concern themselves, inter alia, with predicting the length scales at which the EPE production and conversion to EKE occur. An early theory proposed by Salmon [3] has become a touchstone of the modern understanding [4], having been revised and expanded by a wide range of authors. This early theory, and many successors, is formulated in the context of a highly idealized two-layer quasigeostrophic model whose vertical structure can be formulated in terms of two ‘modes’: a barotropic mode that is depth-independent, and a baroclinic mode that is antisymmetric in the vertical. (The term ‘mode’ here simply means elements of a basis, not the normal modes of a mechanical system, cf. [5,6].) These modes diagonalize the energy, and Salmon’s theory was framed in terms of energy in the modes rather than in terms of kinetic and potential energy. The barotropic mode has no potential energy, while the baroclinic mode has both potential and kinetic energy, so the connection of the modal energy perspective to the perspective based on kinetic and potential energy is murky. Nevertheless, the modern theory of geophysical turbulence (cf. [4] Section 9.3) inherits the following picture from Salmon’s theory: EPE is generated at scales larger than the deformation radius (defined below) and then cascades downscale; conversion from EPE to EKE occurs primarily at scales near the deformation radius; EKE is then dissipated by a range of processes, especially by frictional interaction with the bottom boundary at scales somewhat larger than the deformation radius.

The goal of the present investigation is to study the process of conversion from EPE to EKE from a mathematical rather than a physical perspective. Rather than studying the dynamical processes responsible for conversion, the goal is to investigate the mathematical limits on the rate of conversion per unit energy. The approach is based on considering all possible flow configurations at a fixed energy level and then seeking a configuration that will maximize the rate of potential to kinetic energy conversion. These flow configurations need not be steady solutions of the governing equations, or even states that might realistically be produced by the dynamics, so the rates obtained in this way should be understood as a limit on what a realistic flow can achieve, rather than as a prediction of what a realistic flow will achieve. The method is closely related to methods used to provide upper bounds on the heat flux in Rayleigh–Bénard convection, e.g., [7,8,9].

Barham and Grooms studied the mathematical limits on the rate by which eddies can extract potential energy from the mean flow in a fluid model incorporating only the hydrostatic approximation [10], and in a quasigeostrophic model [11]. They found that the eddies are able to extract potential energy from the mean flow at a rate that is independent of the length scale of the eddies. This was something of a null result, in the sense that the bounds on what the flow can do are not closely related to what the flow actually does: a wide range of investigations have found that the EPE is primarily generated at scales larger than the deformation radius [12,13,14,15,16,17,18,19,20]. The results of the present investigation make a closer connection to the phenomenology of geostrophic turbulence than the results of [10,11]: it is shown that conversion from EPE to EKE is much less efficient than the generation of EPE at scales large compared to the deformation radius. While the results do not provide a rigorous mathematical explanation for why EPE is typically generated at scales larger than the deformation radius (there are many dynamical theories to predict this, e.g., [3,14,15]), it does connect to the downscale cascade of EPE and to the fact that conversion from potential to kinetic energy is observed to be weak at scales larger than the deformation radius.

A limitation of the current results is that the maximal rates of conversion obtained through the analysis are not proven here to be truly maximal. The conversion rate is a function of the flow configuration, and there are an infinite number of flow configurations that are stationary points of the function, i.e., these flow configurations correspond to local extrema or saddle points of the function. The analysis here identifies an infinite set of these stationary points that can be analyzed using Fourier methods and linear eigenvalue theory. Whether there are other flow configurations that could result in higher conversion rates remains an open question. The difficulty is related to the fact that conversion, in a quasigeostrophic approximation, is a cubic function of the flow state, which is more difficult to analyze than the production of EPE (cf. [10,11]), which is a quadratic function of the flow state.

The paper is organized as follows. In Section 2, conversion is studied in the context of the hydrostatic Eady problem, similar to the development in [10]. As in [10], the results are complicated by the presence of unrealistic ageostrophic flow configurations, so, as in [11], the quasigeostrophic (QG) version of the problem is studied next. The QG problem is studied analytically in Section 3 and numerically in Section 4. Conclusions are offered in Section 5.

2. The Hydrostatic Eady Problem

The linear perturbation equations in the nondimensional Eady problem with hydrostatic and Boussinesq approximations are [10,21]

$$\begin{array}{cc}\hfill \left({\partial }_t+z{\partial }_x\right){\mathit{u}}^{\prime }+ϵw^{\prime }\widehat{\mathit{x}}+{ϵ}^{-1}{(\widehat{\mathit{z}}\times {\mathit{u}}^{\prime })}_h& =-{ϵ}^{-1}{\nabla }_h{p}^{\prime }\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\partial }_z{p}^{\prime }& =b^{\prime }\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \left({\partial }_t+z{\partial }_x\right)b^{\prime }-v^{\prime }+w^{\prime }& =0\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\nabla }_h·{\mathit{u}}^{\prime }+ϵ{\partial }_z{w}^{\prime }& =0\hfill \end{array}$$

(4)

where z is the coordinate along the axis of rotation and gravity and x is the coordinate along the direction of mean flow. The domain is periodic in the x and y directions and $z\in [0,1]$. The subscript ${}_h$ denotes the horizontal component of a vector, e.g., ${\nabla }_h=({\partial }_x,{\partial }_y)$.

The buoyancy frequency is N and the dimensional background velocity is $\overline{u}=\mathrm{\Lambda }z$. The Richardson number is ${ϵ}^{-2}$, where $ϵ=\mathrm{\Lambda }/N$. The time scale for nondimensionalization is $N/\left(f\mathrm{\Lambda }\right)$, where f is the Coriolis parameter; the vertical length scale is H; the horizontal length scale is the deformation radius $NH/f$. The horizontal velocity scale is $H\mathrm{\Lambda }$; the vertical velocity scale is $Hf{ϵ}^2$; the buoyancy scale is $HN\mathrm{\Lambda }$.

The perturbation kinetic energy equation is obtained by taking the dot product of (1) with ${\mathit{u}}^{\prime }=(u^{\prime },v^{\prime })$ and integrating over the domain (with several integrations by parts)

$$\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V{\left(u^{\prime }\right)}^2+{\left(v^{\prime }\right)}^2\mathrm{d}V={\int }_V{w}^{\prime }{b}^{\prime }-ϵw^{\prime }{u}^{\prime }\mathrm{d}V$$

(5)

where ${\int }_V$ denotes an integral over the physical domain V. The perturbation available potential energy equation is obtained by multiplying (3) by $b^{\prime }$ and integrating over the domain

$$\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V{\left(b^{\prime }\right)}^2\mathrm{d}V={\int }_V{v}^{\prime }{b}^{\prime }-w^{\prime }{b}^{\prime }\mathrm{d}V.$$

(6)

The term of interest here is the conversion from potential to kinetic energy

$$C={\int }_V{w}^{\prime }{b}^{\prime }\mathrm{d}V.$$

(7)

Because this is a hydrostatic model, the vertical velocity $w^{\prime }$ is obtained from the horizontal velocity using

$$w^{\prime }=-{ϵ}^{-1}{\int }_0^z{\partial }_x{u}^{\prime }+{\partial }_y{v}^{\prime }\mathrm{d}s$$

(8)

where the variable s is a stand-in for the vertical coordinate.

The goal is to obtain a configuration of $(u^{\prime },v^{\prime },b^{\prime })$ that maximizes the conversion C at a fixed energy level $E_0$. We therefore define the Lagrangian

$$I[u^{\prime },v^{\prime },b^{\prime },\lambda ]=C-\lambda (E-E_0)$$

(9)

where $\lambda$ is the Lagrange multiplier and the total energy is

$$E=\frac{1}{2}{\int }_V{\left(u^{\prime }\right)}^2+{\left(v^{\prime }\right)}^2+{\left(b^{\prime }\right)}^2\mathrm{d}V.$$

(10)

The Lagrangian is a quadratic function of the buoyancy and velocity and the energy level $E_0$ can be scaled out of the problem, so, without loss of generality, let $E_0=1$.As usual, the Euler–Lagrange equations are derived by finding conditions that describe stationary points of the Lagrangian.

The Euler–Lagrange equations for this constrained optimization problem are obtained as follows. We first consider the Fréchet derivative of the energy, which is simply

$$\begin{array}{c}\hfill \mathrm{d}E={\int }_V{u}^{\prime }\delta u+v^{\prime }\delta v+b^{\prime }\delta b\mathrm{d}V.\end{array}$$

(11)

To derive the Fréchet derivative of the shear production, start from the expression

$$\mathrm{d}C={\int }_V{b}^{\prime }\delta w+w^{\prime }\delta b\mathrm{d}V.$$

(12)

To proceed, we need the following simple integration by parts identity, also used in [10]

$${\int }_0^{1}g\left(z\right){\int }_0^{z}h\left(s\right)\mathrm{d}s\mathrm{d}z=\left({\int }_0^{1}g\left(z\right)\mathrm{d}z\right)\left({\int }_0^{1}h\left(z\right)\mathrm{d}z\right)-{\int }_0^{1}h\left(z\right){\int }_0^{z}g\left(s\right)\mathrm{d}s\mathrm{d}z$$

(13)

which is valid for integrable functions g and h. This expression enables the following manipulation of the first term in (12):

$$\begin{array}{c}ϵ{\int }_0^1{b}^{\prime }\delta w\mathrm{d}z=-{\int }_0^1{b}^{\prime }{\int }_0^z\left({\partial }_x\delta u+{\partial }_y\delta v\right)\mathrm{d}s\mathrm{d}z\\ =-\left({\int }_0^1{b}^{\prime }\mathrm{d}z\right)\left({\int }_0^1\left({\partial }_x\delta u+{\partial }_y\delta v\right)\mathrm{d}z\right)+{\int }_0^1\left({\partial }_x\delta u+{\partial }_y\delta v\right){\int }_0^z{b}^{\prime }\mathrm{d}s\mathrm{d}z.\end{array}$$

(14)

The fact that $\delta w=0$ at both the upper and lower boundaries sets ${\int }_0^1{\partial }_x\delta u+{\partial }_y\delta v\mathrm{d}z=0$. The Fréchet derivative of the conversion is therefore

$$\mathrm{d}C=-{ϵ}^{-1}{\int }_V\left[\delta u\left({\int }_0^z{\partial }_x{b}^{\prime }\right)+\delta v\left({\int }_0^z{\partial }_y{b}^{\prime }\right)+\delta b\left({\int }_0^z{\nabla }_h·{\mathit{u}}^{\prime }\right)\right]\mathrm{d}V.$$

(15)

Configurations of $(u^{\prime },v^{\prime },b^{\prime })$ that are stationary with respect to the conversion rate subject to the condition of fixed energy satisfy

$$\mathrm{d}I=\mathrm{d}C-\lambda \mathrm{d}E=0$$

(16)

for all $(\delta u,\delta v,\delta b)$. The Euler–Lagrange equations are therefore

$$\begin{array}{c}{\int }_0^z{\partial }_x{b}^{\prime }\mathrm{d}s=ϵ\lambda u^{\prime }\end{array}$$

(17)

$$\begin{array}{c}{\int }_0^z{\partial }_y{b}^{\prime }\mathrm{d}s=ϵ\lambda v^{\prime }\end{array}$$

(18)

$$\begin{array}{c}{\int }_0^z{\nabla }_h·{\mathit{u}}^{\prime }\mathrm{d}s=ϵ\lambda b^{\prime }.\end{array}$$

(19)

Perturbation fields $(u^{\prime },v^{\prime },b^{\prime })$ that satisfy these equations are associated with conversion rates $\lambda$. Note that horizontally incompressible flow with $b^{\prime }=\lambda =0$ is a stationary point of the Lagrangian, but an uninteresting one since it has no conversion.

A system of partial differential equations for the optimal flow configurations can be obtained by taking the partial derivative of these equations with respect to z and then condensing to a single equation

$${\nabla }_h^2{b}^{\prime }={ϵ}^2{\lambda }^2{\partial }_z^2{b}^{\prime }.$$

(20)

Noting that (19) implies that, for these optimal configurations $b^{\prime }=0$ on the top and bottom boundaries, we can expand solutions as

$$b^{\prime }={\widehat{b}}_{\mathit{k},n}{e}^{\mathrm{i}\mathit{k}·\mathit{x}}sin\left(n\pi z\right)+\mathrm{c}.\mathrm{c}.$$

(21)

where ‘c.c.’ denotes the complex conjugate and the amplitude of ${\widehat{b}}_{\mathit{k},n}$ is determined by the condition of unit total energy. The conversion rate $\lambda$ satisfies

$$\lambda =\pm \frac{\left|\mathit{k}\right|}{ϵn\pi }.$$

(22)

The conversion rate exhibits an ultraviolet catastrophe where the conversion rate approaches infinity at small scales. Ultraviolet catastrophes can call into question the well-posedness of a system of equations; in this case, this is not a concern. Conversion is not a component of the total energy budget, since the contributions to the kinetic and potential energy budgets cancel on adding, so this ultraviolet catastrophe does not cause unbounded growth of energy at small scales. Instead, it means that certain small-scale flow configurations will almost instantaneously convert their potential energy to kinetic energy, which is of course not problematic for the dynamics.

The flow configurations associated with these high conversion rates have divergent horizontal velocities: (19) implies

$${\nabla }_h·{\mathit{u}}^{\prime }=ϵ\lambda {\partial }_z{b}^{\prime }=\pm \frac{\left|\mathit{k}\right|}{n\pi }{\partial }_z{b}^{\prime }.$$

(23)

As $\left|\mathit{k}\right|\to \infty$, the magnitude of the divergence grows without bound, which causes the magnitude of the vertical velocity to also grow without bound per (8). These flow configurations are mathematically permissible—vertical velocity is not included in the conserved energy for this model, so infinite vertical velocity is not precluded by the constraint of finite energy—but physically unrealistic for balanced ocean dynamics, which have small vertical velocities. Thus, as in [11], the next section pursues the question in the context of a quasigeostrophic model. An alternative approach, not taken here, would be to analyze the conversion in the non-hydrostatic model; this would include the vertical velocity as part of the kinetic energy, which would prevent it from growing unboundedly in the progression to small scales.

3. Conversion in the Quasigeostrophic Model

The nondimensional eddy vorticity and buoyancy evolution equations in an inviscid quasigeostrophic (QG) approximation in the presence of a zonal mean flow are

$$\begin{array}{cc}\hfill {\partial }_t\omega +J[\psi ,\omega ]+\overline{u}\left(z\right){\partial }_x\omega +\beta {\partial }_x\psi -{\partial }_z{w}^{\prime }& =0\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\partial }_t{b}^{\prime }+J[\psi ,b^{\prime }]+\overline{u}\left(z\right){\partial }_x{b}^{\prime }-v^{\prime }{\partial }_z\overline{u}\left(z\right)+w^{\prime }{N}^2\left(z\right)& =0\hfill \end{array}$$

(25)

where $\omega ={\nabla }_h^2\psi$ is the eddy vorticity, $b^{\prime }={\partial }_z\psi$ is the eddy buoyancy, $N\left(z\right)$ is the nondimensional buoyancy frequency, and $J[\psi ,\omega ]={\partial }_x\psi {\partial }_y\omega -{\partial }_y\psi {\partial }_x\omega ={\mathit{u}}^{\prime }·\nabla \omega$. The depth H is used to nondimensionalize z, and the deformation scale $L=N_{0}H/f$ is used to nondimensionalize x and y, where f is the Coriolis parameter at a fixed reference latitude, and $N_0$ is the maximum value of the dimensional buoyancy frequency. The time scale is $N_0^{-1}$; the scale of $\psi$ is $N_0{L}^2$ and of $w^{\prime }$ is $NL$. The nondimensional planetary vorticity gradient is $\beta ={\beta }_{0}L/N$, where ${\beta }_0$ is the meridional rate of change of the Coriolis parameter at the same reference latitude.

These two evolution equations for derivatives of $\psi$ could, in principle, be incompatible; the incompatibility is prevented by vertical velocity $w^{\prime }$, which acts to keep the two evolution equations consistent. The condition on $w^{\prime }$ that keeps the two equations consistent is obtained by applying $-{\partial }_z$ to (24), applying ${\nabla }_h^2$ to (25), and adding the results. This leads to the Omega equation [22] for $w^{\prime }$:

$$N^2\left(z\right){\nabla }_h^2{w}^{\prime }+{\partial }_z^2{w}^{\prime }={\partial }_z\left[J[\psi ,\omega ]+\beta v\right]-{\nabla }_h^2\left[J[\psi ,b^{\prime }]\right]+2\left({\partial }_z\overline{u}\left(z\right)\right){\partial }_x\omega .$$

(26)

The evolution of kinetic and available potential energy in the QG model can be obtained, respectively, by multiplying (24) by $-\psi$ and integrating over the domain, and by multiplying (25) by $b^{\prime }/N^2\left(z\right)$ and integrating over the domain. The results are

$$\begin{array}{cc}\hfill \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V{\left(u^{\prime }\right)}^2+{\left(v^{\prime }\right)}^2\mathrm{d}V& ={\int }_V{w}^{\prime }{b}^{\prime }\mathrm{d}V\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V\frac{{\left(b^{\prime }\right)}^2}{N^2\left(z\right)}\mathrm{d}V& ={\int }_V\frac{v^{\prime }{b}^{\prime }{\partial }_z\overline{u}\left(z\right)}{N^2\left(z\right)}-w^{\prime }{b}^{\prime }\mathrm{d}V\hfill \end{array}$$

(28)

where one integration by parts was performed to obtain the kinetic energy equation, using $w^{\prime }=0$ conditions on the boundaries.

3.1. Derivation of the Lagrangian

The form of the potential to kinetic conversion and of the total energy are the same as in the hydrostatic Eady problem, with the important difference that in the QG problem, the velocity and buoyancy are all derived from a single variable $\psi$. In terms of the QG streamfunction, we can write the total energy as

$$E=\frac{1}{2}{\int }_V{\left|{\nabla }_h\psi \right|}^2+\frac{1}{N^2\left(z\right)}{\left({\partial }_z\psi \right)}^2\mathrm{d}V.$$

(29)

Unlike in the hydrostatic Eady problem, conversion is a cubic function of $\psi$, although it can be split into a sum of cubic and quadratic components. For notational convenience, $\mathrm{\Omega }$ is defined to be the differential operator acting on $w^{\prime }$, i.e.,

$$\mathrm{\Omega }\left[w^{\prime }\right]=N^2\left(z\right){\nabla }_h^2{w}^{\prime }+{\partial }_z^2{w}^{\prime }.$$

(30)

The operator $\mathrm{\Omega }$ is self-adjoint when acting on the Sobolev space of functions with homogeneous Dirichlet boundary conditions whose weak derivatives of order $\le 2$ belong to $L^2\left(V\right)$, which implies that ${\mathrm{\Omega }}^{-1}$ is also self-adjoint (cf. [23], examples 10.1 and 10.4). Noting that $\mathrm{\Omega }$ is linear, we can split conversion as follows:

$$\begin{array}{cc}\hfill C& =C_3+C_2\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill C_3& ={\int }_V{\mathrm{\Omega }}^{-1}\left[{\partial }_z\left[J[\psi ,{\nabla }^2\psi ]\right]-{\nabla }_h^2\left[J[\psi ,{\partial }_z\psi ]\right]\right]{\partial }_z\psi \mathrm{d}V\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill C_2& ={\int }_V{\mathrm{\Omega }}^{-1}\left[\beta {\partial }_z{\partial }_x\psi +2\left({\partial }_z\overline{u}\left(z\right)\right){\partial }_x{\nabla }^2\psi \right]{\partial }_z\psi \mathrm{d}V\hfill \end{array}$$

(33)

where $C_3$ is cubic and $C_2$ is quadratic. The fact that ${\mathrm{\Omega }}^{-1}$ is self-adjoint and commutes with ${\partial }_x$ implies that

$${\int }_V{\mathrm{\Omega }}^{-1}\left[{\partial }_x{\partial }_z\psi \right]{\partial }_z\psi \mathrm{d}V=\frac{1}{2}{\int }_V{\partial }_x\left({\mathrm{\Omega }}^{-1}\left[{\partial }_z\psi \right]{\partial }_z\psi \right)\mathrm{d}V=0.$$

(34)

From this, we conclude that the vertical velocity generated by the $\beta$ term in the Omega equation does not lead to any conversion, and thus

$$C_2=2{\int }_V{\mathrm{\Omega }}^{-1}\left[\left({\partial }_z\overline{u}\left(z\right)\right){\partial }_x{\nabla }^2\psi \right]{\partial }_z\psi \mathrm{d}V.$$

(35)

Note that $\beta$ has a profound influence on the dynamics, and thus an indirect influence on the rate of conversion of EPE to EKE.

With this notation, the problem of interest is to maximize the rate of conversion at fixed energy level $E_0$. The Lagrangian for this constrained optimization problem is

$$I\left[\psi \right]=C_3+C_2-\lambda (E-E_0).$$

(36)

Unlike the non-QG case, the presence of a cubic term implies that the energy level $E_0$ cannot be set to unity by a rescaling of $\psi$.

3.2. Analysis of the Cubic Term

The presence of the cubic term $C_3$ implies that the Euler–Lagrange equations for this Lagrangian constitute a quadratic eigenvalue problem, unlike the linear eigenvalue problems associated with EPE production analyzed in [10,11]. However, the cubic term has a curious property: for all flow configurations $\psi$ for which $\mathrm{d}{C}_2-\lambda \mathrm{d}E=0$, we have that $C_3=\mathrm{d}{C}_3=0$. This implies that any conclusions drawn from analyzing only the quadratic part of the Lagrangian apply also to the full Lagrangian, with the only caveat being that the optimal conversion rates discovered in this way may not be globally optimal for the full problem.

The following theorem justifies the foregoing statement about the cubic part of the Lagrangian.

Theorem 1

(Univariate degeneracy of $C_3$). Any flow configuration that satisfies $\psi (x,y,z)=p(k_{x}x+k_{y}y,z)$ sets $C_3=\mathrm{d}{C}_3=0$ for any $k_x$ and $k_y$.

Proof. 

The fact that $C_3=0$ for any $\psi (x,y,z)=p(k_{x}x+k_{y}y,z)$ follows from the fact that the Jacobian advection operator is zero for two univariate arguments: $J[p(k_{x}x+k_{y}y)$, $r(k_{x}x+k_{y}y)]=0$.

To prove that $\mathrm{d}{C}_3=0$ for any $\psi (x,y,z)=p(k_{x}x+k_{y}y,z)$, begin by using the fact that $\mathrm{\Omega }$ is self-adjoint to write

$$C_3={\int }_V\left({\partial }_z\left[J[\psi ,{\nabla }^2\psi ]\right]-{\nabla }_h^2\left[J[\psi ,{\partial }_z\psi ]\right]\right){\mathrm{\Omega }}^{-1}\left[{\partial }_z\psi \right]\mathrm{d}V.$$

(37)

Integrations by parts in the vertical and horizontal now produce the following, where horizontal periodicity has been used, as well as the fact that ${\mathrm{\Omega }}^{-1}[·]$ produces a function that is zero on the upper and lower boundaries:

$$C_3=-{\int }_{V}J[\psi ,{\nabla }^2\psi ]\left({\partial }_z{\mathrm{\Omega }}^{-1}\left[{\partial }_z\psi \right]\right)+J[\psi ,{\partial }_z\psi ]\left({\nabla }_h^2{\mathrm{\Omega }}^{-1}\left[{\partial }_z\psi \right]\right)\mathrm{d}V.$$

(38)

With this expression in hand, we can consider the Fréchet derivative of $C_3$

$$\begin{array}{cc}\hfill \mathrm{d}{C}_3& =-{\int }_{V}J[\psi ,{\nabla }^2\psi ]\left({\partial }_z{\mathrm{\Omega }}^{-1}\left[{\partial }_z\delta \psi \right]\right)+J[\psi ,{\partial }_z\psi ]\left({\nabla }_h^2{\mathrm{\Omega }}^{-1}\left[{\partial }_z\delta \psi \right]\right)\mathrm{d}V\hfill \\ \hfill & -{\int }_{V}J[\delta \psi ,{\nabla }^2\psi ]\left({\partial }_z{\mathrm{\Omega }}^{-1}\left[{\partial }_z\psi \right]\right)+J[\delta \psi ,{\partial }_z\psi ]\left({\nabla }_h^2{\mathrm{\Omega }}^{-1}\left[{\partial }_z\psi \right]\right)\mathrm{d}V\hfill \\ \hfill & -{\int }_{V}J[\psi ,{\nabla }^2\delta \psi ]\left({\partial }_z{\mathrm{\Omega }}^{-1}\left[{\partial }_z\psi \right]\right)+J[\psi ,{\partial }_z\delta \psi ]\left({\nabla }_h^2{\mathrm{\Omega }}^{-1}\left[{\partial }_z\psi \right]\right)\mathrm{d}V.\hfill \end{array}$$

(39)

The first line is zero whenever $\psi (x,y,z)=p(k_{x}x+k_{y}y,z)$ because of the property of the Jacobian noted above. The remaining terms can all be expressed in the form

$${\int }_{V}qJ[r,s]\mathrm{d}V$$

(40)

where either r or s is a perturbation. (Note that q does not refer to potential vorticity in this expression.) The identity $J[r,s]=-J[s,r]$ allows us to consider s to be the perturbation without loss of generality. Integration by parts puts these terms into the form

$${\int }_{V}qJ[r,s]=-{\int }_{V}s\left({\partial }_y\left(q{\partial }_{x}r\right)-{\partial }_x\left(q{\partial }_{y}r\right)\right)\mathrm{d}V.$$

(41)

If $\psi =p(k_{x}x+k_{y}y,z)$, then, for all of terms on the second and third lines of (39), q and r are of the form $q=q(k_{x}x+k_{y}y,z)$ and $r=r(k_{x}x+k_{y}y,z)$. For functions of this form,

$${\partial }_y\left(q{\partial }_{x}r\right)-{\partial }_x\left(q{\partial }_{y}r\right)=0$$

(42)

which implies that $\mathrm{d}{C}_3=0$. □

Note that no claim is made that these are the only flow configurations for which $\mathrm{d}{C}_3=0$. A complete, rigorous analysis of the stationary points of the full Lagrangian is outside the scope of this investigation.

3.3. Euler–Lagrange Equations for the Quadratic Part

Having established the foregoing property of $C_3$, we next derive Euler–Lagrange equations for the quadratic component of the Lagrangian

$$I_2\left[\psi \right]=C_2-\lambda (E-E_0).$$

(43)

As with the hydrostatic Eady problem, the fact that this Lagrangian is quadratic implies that the energy level $E_0$ can be set to unity without loss of generality. The Fréchet derivative of the energy is

$$\mathrm{d}E=-{\int }_V\left({\nabla }_h^2\psi +{\partial }_z\left(\frac{1}{N^2\left(z\right)}{\partial }_z\psi \right)\right)\delta \psi \mathrm{d}V.$$

(44)

The Fréchet derivative of the quadratic component of conversion is

$$\mathrm{d}{C}_2=-2{\int }_V\left({\partial }_z{\mathrm{\Omega }}^{-1}\left[\left({\partial }_z\overline{u}\left(z\right)\right){\partial }_x{\nabla }_h^2\psi \right]+\left({\partial }_z\overline{u}\left(z\right)\right){\mathrm{\Omega }}^{-1}\left[{\partial }_x{\nabla }_h^2{\partial }_z\psi \right]\right)\delta \psi \mathrm{d}V.$$

(45)

Expand both $\psi$ and $\delta \psi$ as Fourier series,

$$\psi (x,y,z)=\frac{1}{\sqrt{L_x{L}_y}}\sum _{\mathit{k}}{\widehat{\psi }}_{\mathit{k}}\left(z\right)e^{\mathrm{i}\mathit{k}·\mathit{x}}$$

(46)

$$\delta \psi (x,y,z)=\frac{1}{\sqrt{L_x{L}_y}}\sum _{\mathit{k}}{\widehat{\delta \psi }}_{\mathit{k}}\left(z\right)e^{\mathrm{i}\mathit{k}·\mathit{x}}$$

(47)

where $L_x$ and $L_y$ are the nondimensional length and width of the periodic domain. With this expansion, we have the following representations of $\mathrm{d}E$ and $\mathrm{d}{C}_2$:

$$\mathrm{d}E=-\sum _{\mathit{k}}{\int }_0^1\left(-{\left|\mathit{k}\right|}^2{\widehat{\psi }}_{\mathit{k}}+{\partial }_z\left(\frac{1}{N^2\left(z\right)}{\partial }_z{\widehat{\psi }}_{\mathit{k}}\right)\right){\widehat{\delta \psi }}_{\mathit{k}}^{*}\mathrm{d}z$$

(48)

$$\mathrm{d}{C}_2=2\mathrm{i}\sum _{\mathit{k}}{k}_x{\left|\mathit{k}\right|}^2{\int }_0^1\left({\partial }_z{\mathrm{\Omega }}^{-1}\left[\left({\partial }_z\overline{u}\left(z\right)\right){\widehat{\psi }}_{\mathit{k}}\right]+\left({\partial }_z\overline{u}\left(z\right)\right){\mathrm{\Omega }}^{-1}\left[{\partial }_z{\widehat{\psi }}_{\mathit{k}}\right]\right){\widehat{\delta \psi }}_{\mathit{k}}^{*}.$$

(49)

where ${\widehat{\delta \psi }}_{\mathit{k}}^{*}$ is the complex conjugate of ${\widehat{\delta \psi }}_{\mathit{k}}$. Stationary points of the quadratic Lagrangian are thus single Fourier modes that satisfy the Euler–Lagrange equations

$$\begin{array}{c}2\mathrm{i}{k}_x{\left|\mathit{k}\right|}^2\left({\partial }_z{\mathrm{\Omega }}^{-1}\left[\left({\partial }_z\overline{u}\left(z\right)\right){\widehat{\psi }}_{\mathit{k}}\right]+\left({\partial }_z\overline{u}\left(z\right)\right){\mathrm{\Omega }}^{-1}\left[{\partial }_z{\widehat{\psi }}_{\mathit{k}}\right]\right)=\hfill \\ \hfill \lambda \left(-{\left|\mathit{k}\right|}^2{\widehat{\psi }}_{\mathit{k}}+{\partial }_z\left(\frac{1}{N^2\left(z\right)}{\partial }_z{\widehat{\psi }}_{\mathit{k}}\right)\right).\end{array}$$

(50)

Because stationary points of the quadratic part of the Lagrangian are single Fourier modes, they are of the form $\psi (x,y,z)=p(k_{x}x+k_{y}y,z)$. Theorem 1 implies that $C_3=\mathrm{d}{C}_3=0$ for these configurations, so these flow configurations are also stationary points of the full Lagrangian.

3.4. Asymptotic Analysis

Exact solutions of the Euler–Lagrange equations for the quadratic part (50) remain elusive due to the presence of the ${\mathrm{\Omega }}^{-1}$ operator. Before proceeding directly to numerical investigations of conversion, it is valuable to provide some asymptotic analysis to complement the numerics. Two limits of interest are large and small scales. To analyze these limits, we return to the Omega Equation (26), written for a single Fourier mode and omitting the passive $\beta$ term

$$\left[-N^2\left(z\right){\left|\mathit{k}\right|}^2+{\partial }_z^2\right]{\widehat{w}}_{\mathit{k}}=-2\mathrm{i}{k}_x{\left|\mathit{k}\right|}^2\left({\partial }_z\overline{u}\left(z\right)\right){\widehat{\psi }}_{\mathit{k}}.$$

(51)

Clearly, ${\widehat{w}}_{\mathit{k}}\to 0$ as $k_x\to 0$, which implies that conversion also goes to zero as $k_x\to 0$. More generally, at large scales, i.e., $\mathit{k}\ll 1$, it is clear that ${\widehat{w}}_{\mathit{k}}\sim k_x{\left|\mathit{k}\right|}^2$, so conversion must go to zero at scales much larger than the deformation radius.

For small scales, $\mathit{k}\gg 1$, it is evident that

$${\widehat{w}}_{\mathit{k}}\sim 2\mathrm{i}{k}_x\frac{{\partial }_z\overline{u}\left(z\right)}{N^2\left(z\right)}{\widehat{\psi }}_{\mathit{k}}.$$

(52)

Exponentially thin boundary layers return ${\widehat{w}}_{\mathit{k}}$ to zero on the boundary if necessary, but these layers contribute an asymptotically subdominant amount to the total conversion. The total conversion for large $\mathit{k}$ thus scales as twice the real part of

$$2\mathrm{i}{k}_x{\int }_0^1\frac{{\partial }_z\overline{u}\left(z\right)}{N^2\left(z\right)}{\widehat{\psi }}_{\mathit{k}}{\partial }_z{\widehat{\psi }}_{\mathit{k}}^{*}\mathrm{d}z.$$

(53)

The elementary inequalities $2\left|ab\right|\le (a^2+b^2)$ and $k_x^2\le {\left|\mathit{k}\right|}^2$ imply that the conversion in this limit is bounded by

$$2{\int }_0^1\frac{\left|{\partial }_z\overline{u}\left(z\right)\right|}{N\left(z\right)}\left({\left|\mathit{k}\right|}^2{\left|{\widehat{\psi }}_{\mathit{k}}\right|}^2+\frac{{\left|{\partial }_z{\widehat{\psi }}_{\mathit{k}}\right|}^2}{N^2\left(z\right)}\right)\mathrm{d}z$$

(54)

which further implies that the rate of conversion per unit energy is bounded by twice the Richardson number

$$2{∥\frac{{\partial }_z\overline{u}\left(z\right)}{N\left(z\right)}∥}_{\infty }.$$

(55)

This boundedness at small scales is in contrast to the behavior seen in the hydrostatic Eady problem in Section 2, and mirrors the behavior of EPE production rates found in [11].

4. Numerical Investigation into Conversion Rates

This section provides a numerical investigation of maximal conversion rates, applying only to the quadratic component of conversion in the presence of a zonal mean flow. Rather than discretize the Euler–Lagrange Equations (50), the analysis is framed directly in terms of the conversion for a single Fourier mode

$${\int }_V{w}^{\prime }{b}^{\prime }\mathrm{d}V=2\mathrm{Re}\left[{\int }_0^1{\widehat{w}}_{\mathit{k}}{\widehat{b}}_{\mathit{k}}^{*}\mathrm{d}z\right].$$

(56)

By applying standard equispaced, second-order, centered finite differences to the Omega equation for a single Fourier mode (51) and to ${\partial }_z{\widehat{\psi }}_{\mathit{k}}={\widehat{b}}_{\mathit{k}}$, one obtains an expression of the form

$$C_2\approx {\widehat{\mathit{\psi }}}_{\mathit{k}}^{*}\mathbf{C}{\widehat{\mathit{\psi }}}_{\mathit{k}}$$

(57)

where $\mathbf{C}$ is a Hermitian matrix and ${\widehat{\mathit{\psi }}}_{\mathit{k}}$ is a vector of values of ${\widehat{\psi }}_{\mathit{k}}$ at $N_z$ equispaced points on $z\in (0,1)$. The energy may similarly be numerically approximated using

$$E\approx {\widehat{\mathit{\psi }}}_{\mathit{k}}^{*}\mathbf{E}{\widehat{\mathit{\psi }}}_{\mathit{k}}$$

(58)

where $\mathbf{E}=({\left|\mathit{k}\right|}^2\mathbf{I}-\mathbf{L})/\left(4N_z\right)$ and $\mathbf{L}$ is a discrete approximation to the operator ${\partial }_z\left(N^{-2}\left(z\right){\partial }_z(·)\right)$. With these discrete approximations, the generalized eigenvalues $\lambda$ of

$$\mathbf{C}{\widehat{\mathit{\psi }}}_{\mathit{k}}=\lambda \mathbf{E}{\widehat{\mathit{\psi }}}_{\mathit{k}}$$

(59)

are the conversion rates per unit energy associated with flow configurations given by the generalized eigenvectors ${\widehat{\mathit{\psi }}}_{\mathit{k}}$. The fact that $\mathbf{C}$ and $\mathbf{E}$ are both Hermitian and that $\mathbf{E}$ is positive definite imply that there are exactly $N_z$ real generalized eigenvalues $\lambda$ for each Fourier mode $\mathit{k}$.

For comparison, the same treatment is applied to the rate of eddy potential energy production following [11]; the continuous expression is

$$P={\int }_V{v}^{\prime }{b}^{\prime }\frac{{\partial }_z\overline{u}\left(z\right)}{N^2\left(z\right)}\mathrm{d}V=2\mathrm{Re}\left[\mathrm{i}{k}_x{\int }_0^1{\widehat{\psi }}_{\mathit{k}}{\partial }_z{\widehat{\psi }}_{\mathit{k}}\frac{{\partial }_z\overline{u}\left(z\right)}{N^2\left(z\right)}\mathrm{d}z\right].$$

(60)

The discrete approximation takes the form

$$P\approx {\widehat{\mathit{\psi }}}_{\mathit{k}}^{*}\mathbf{P}{\widehat{\mathit{\psi }}}_{\mathit{k}}$$

(61)

and the optimal rates of production per unit energy are found as generalized eigenvalues $\lambda$ of the Hermitian generalized eigenvalue problem

$$\mathbf{P}{\widehat{\mathit{\psi }}}_{\mathit{k}}=\lambda \mathbf{E}{\widehat{\mathit{\psi }}}_{\mathit{k}}.$$

(62)

The code used to generate the $\mathbf{C}$, $\mathbf{P}$, and $\mathbf{E}$ matrices can be found in [24].

Results are computed for three mean flow profiles:

$$\begin{array}{cc}\hfill \mathrm{Eady}\overline{u}\left(z\right)& =z-\frac{1}{2}\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill \mathrm{Phillips}\overline{u}\left(z\right)& =-\frac{1}{2}cos\left(\pi z\right)\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill \mathrm{Ocean}\mathrm{Charney}\overline{u}\left(z\right)& =\frac{3}{2}{z}^2-\frac{1}{2}\hfill \end{array}$$

(65)

all with $N\left(z\right)=1$ and using $N_z=200$ points. These three profiles represent three different ways of violating the Charney–Stern–Pedlosky criteria for baroclinic stability ([4], Section 6.4.3; see also [25]). For each mean flow configuration and over a range of wavenumbers $\mathit{k}$, the optimal rates of conversion and production per unit energy are computed.

These conversion rates are shown as functions of $k_x$ and $k_y$ in Figure 1, and as a function of $k_x$ for $k_y=0$ in Figure 2. As predicted by the asymptotic theory, the conversion rates go to zero as $k_x\to 0$, which is in contrast with the production rates. For the latter, there is a singularity at $\mathit{k}=\mathbf{0}$ such that, for $k_y=0$ and $k_x\to 0$, the production rate remains nonzero, whereas for $k_y\ne 0$, the production rate goes to zero as $k_x\to 0$. This contrast is a key finding of this investigation: QG flow is able to produce EPE at scales larger than the deformation radius much more efficiently than it is able to convert EPE to EKE, at least when considering only linear conversion processes. However, this mismatch only occurs near the $k_y=0$ axis; at large scales away from this axis, the maximal rates of conversion and production are both weak.

Figure 2 compares the conversion and production rates on the same plot for $k_y=0$. For all mean flow configurations, the maximum possible rate of EPE production is larger than the maximum possible rate of linear conversion at large scales. The rate of conversion increases towards small scales, and the maximum possible rate of conversion exceeds that of production at $k_x<10$ for the Eady and Phillips mean flows; for the Charney mean flow, conversion exceeds production at much smaller scales (not shown). None of the computed conversion rates violates the asymptotic bound derived in Section 3.4.

Figure 3 shows the flow configurations $\psi$ that optimize conversion and production for all three mean flow configurations at $k_x=1$, $k_y=0$. At this wavenumber, the maximal production rate is significantly greater than the maximal conversion rate, and the flow configurations that achieve these maximal rates are markedly different: the vertical tilt of the horizontal flow is opposite for optimal conversion and optimal production.

5. Conclusions

The goal of this paper was to study the mathematical limits on the rate of conversion from eddy potential energy (EPE) to eddy kinetic energy (EKE) with reference to the ocean’s Lorenz energy cycle. The first result parallels one of the results of [10], namely that when the problem is studied using only the hydrostatic approximation, the results are clouded by the presence of highly unrealistic ageostrophic flow configurations that convert EPE to EKE with an efficiency that is unbounded in the progression to infinitesimal horizontal length scales. Thus, following [11], the quasigeostrophic version of the problem is studied next so that, somewhat ironically, the more restrictive approximation might lead to more realistic results.

In the quasigeostrophic approximation, the conversion rate is a cubic function of the eddy flow configuration because the vertical velocity is a quadratic function of the QG streamfunction through the Omega Equation (26) [22]. Nevertheless, in the presence of a mean flow, the vertical velocity can be written as the sum of two components that depend linearly and quadratically on the eddy flow state; the conversion rate can thus also be written as the sum of two components that depend quadratically and cubically on the eddy flow state. The second result of this study is the discovery that eddy flow configurations that correspond to local extrema of the quadratic part of the conversion rate also correspond to local extrema of the full conversion rate including the cubic term. This enables the use of methods based on Fourier analysis and linear eigenvalue problems to find eddy flow configurations that correspond to local maxima of the conversion rate; whether these local maxima are also global maxima remains an open question.

Studying the local maxima only, it is found that the rate of EPE to EKE conversion goes to zero as the horizontal length scale of the eddy flow configuration grows above the deformation radius, and that as the horizontal length scale decreases past the deformation radius, the rate of conversion appears to grow towards an upper bound. Putting these results together with those of [11] on the maximal rate of EPE generation leads to the following picture: at large scales, the maximal rate of EPE generation is much larger than the maximal rate of conversion from EPE to EKE. Thus, if EPE is generated near the maximal rate at large scales, it cannot be converted equally rapidly to EKE; to achieve a statistically steady state, the EPE must be transferred towards smaller scales, where it can be converted to EKE efficiently. As these results are only concerned with limits on what the flow can do, they do not predict what the flow will do. Nevertheless, the results are consistent with the phenomenological theory of QG turbulence.

This study raises some questions that remain unresolved. An obvious example is whether there are eddy flow configurations that can achieve higher conversion rates than the ones identified here. The question can be addressed using more advanced analytical methods, or using numerical methods for partial-differential-equation-constrained optimization. Barotropic shear could be introduced into the mean flow profiles; this would prevent the problem from being partially diagonalized by a Fourier basis, and would similarly require more advanced analytical methods, or using numerical methods for partial-differential-equation-constrained optimization. Aside from the cubic problem, there still remain open questions about the quadratic component of the problem. For example, one might attempt to derive a rigorous bound on the quadratic conversion rate to complement the asymptotic analysis provided in Section 3.4. Finally, the problems in the hydrostatic model could potentially be alleviated by studying conversion in the non-hydrostatic model. In the non-hydrostatic model, all three components of velocity are included in the kinetic energy, which prevents the vertical velocity from growing unboundedly at small scales.