Cumulant expansions for atmospheric flows

Farid Ait-Chaalal; Tapio Schneider; Bettina Meyer; J B Marston

doi:10.1088/1367-2630/18/2/025019

1. Introduction

Atmospheric flows shape Earth's climate and are governed by the equations of fluid dynamics, the Navier–Stokes equations augmented by the Coriolis force and thermodynamic equations (e.g., Ooyama 2001, Vallis 2006, Pauluis 2008), and equations for the microphysical processes describing, for example, the formation and re-evaporation of cloud droplets (Pruppacher et al 1998). They span an enormous range of length scales, from the micrometers of droplet formation to the planetary scale. Temporal variations range from microseconds at the smallest scales to tens of years on the largest scales (e.g., Klein 2010). Atmospheric processes are tightly coupled across all of these scales. For example, cloud droplets scatter sunlight and absorb infrared radiation, thereby affecting Earth's radiative balance globally; conversely, planetary-scale dynamics affect where and how clouds form. Current climate models cannot resolve all relevant scales. They resort to the direct simulation of dynamics on scales of tens of kilometers and larger, while representing smaller-scale processes such as turbulence in clouds and boundary layers semi-empirically (e.g., Beljaars 1992, Garratt 1994, Smith 1997, Lappen and Randall 2001, Soares 2004, Siebesma et al 2007). However, the larger-scale dynamics of weather systems, with timescales of minutes, are simulated explicitly even when only their longer-term statistics—the climate—is ultimately of interest.

Two scientific objectives would be beneficial to achieve. First, it would be desirable to obtain more accurate and more physically motivated models of the interactions between the larger scales that can currently be resolved in climate models and the smaller scales that cannot be resolved. Inaccuracies in how these interactions, in particular in clouds and boundary layers, are represented in climate models are the largest source of uncertainties in climate projections (e.g., Stephens 2005, Bony et al 2006, Soden and Held 2006, Webb et al 2013, Stevens and Bony 2013, Vial et al 2013, Brient et al 2015). Improving the representation of such interactions would have an enormous societal benefit. Second, it would be desirable to devise ways of calculating climate statistics more directly, rather than through the direct simulation of weather systems and accumulation of their statistics, as is currently done. This may in the long run lead to faster climate simulations. In the shorter term, it may lead to insight into how climate is maintained and how it varies on timescales of seasons to millennia.

Both objectives require the development of closure models for the hierarchy of statistical moment or cumulant equations associated with the equations of fluid dynamics. This hierarchy is in principle infinite because of the quadratic nonlinearity of the Navier–Stokes equations. Numerous ways of closing the hierarchy of moment or cumulant equations in a variety of circumstances have been proposed (see, e.g., Frisch 1995, Pope 2000, Lesieur 2008). Many of them concern flows that are assumed statistically homogeneous and isotropic, as an idealized benchmark from which the development of closures for more realistic applications can proceed (e.g., Orszag 1970, Orszag et al 1973). However, closures for homogeneous and isotropic turbulence often may not be easier to obtain than closures for more realistic flows: mean fields in homogeneous and isotropic turbulence can, without loss of generality, be taken to be zero; only higher-order statistics of the flows are of interest. Isotropic turbulence cannot equilibrate with any imposed driving and dissipation through interaction with mean flows; rather, it must equilibrate through nonlinear interactions across scales. By contrast, turbulence in the atmosphere usually interacts strongly with nontrivial mean fields, which include, for example, atmospheric jet streams or the thermal stratification of the atmosphere. Because interactions between turbulent fluctuations and nontrivial mean fields have the potential to be important, many atmospheric flows may be less strongly nonlinear than the oft-studied prototype problems of three-dimensional turbulence (e.g., Pedlosky 1970, Farrell 1987, Randel and Held 1991, Farrell and Ioannou 1993, Schneider and Walker 2006). Moreover, already the mean fields (e.g., mean temperatures and winds) are of primary interest for understanding climate, though, of course, higher moments (e.g., temperature extremes) also remain important to understand.

Because the nonlinearity of turbulent interactions in many atmospheric flows may be limited, truncating the hierarchy of moment or cumulant equations at a low order has potential to be successful. Here we explore the feasibility of truncations at second order—that is, neglecting third-order nonlinearities in second-order covariance equations—in two prototype problems of atmospheric flows. The first is a turbulent convective boundary layer, with scales of motion on the order of meters to a kilometer. The second is a model of large-scale turbulence in the upper atmosphere, with scales of motion on the order of hundreds to thousands of kilometers. These two problems involve disparate phenomenologies and force balances. For example, the boundary layer can be taken to be unaffected by the planetary rotation, whereas the planetary rotation and Coriolis forces are fundamental for the large-scale turbulence in the upper atmosphere. Yet the problems share that turbulent fluctuations interact strongly with a nontrivial mean state—a thermal stratification in the first case, and an atmospheric jet in the second case. Because of the strength of this interaction, truncations of moment equations at second order already achieve some success in capturing the statistics of these flows. It is essential in these truncations that nonlocal and anisotropic covariation of turbulent quantities (e.g., in waves or convective plumes) are retained. Such nonlocal truncation of the moment or cumulant equations at second order is known as second-order cumulant expansion (CE2) (Marston et al 2008, Marston 2012, Srinivasan and Young 2012, Tobias and Marston 2013) or stochastic structural stability theory (S3T) (Farrell and Ioannou 2003, 2007, Bakas and Ioannou 2014, Constantinou et al 2014a). CE2 and S3T differ in that S3T attempts to represent missing eddy–eddy interactions, whereas CE2 sets them to zero.

CE2 is a realizable closure in that its equations are the exact moment equations of a realizable system, the quasi-linear (QL) system that corresponds to the original equations of motion. The QL approximation retains the interaction of turbulent fluctuations with a mean flow but neglects the interactions of turbulent fluctuations among each others. CE2 and S3T were successful in explaining some aspects of zonal jet dynamics in rotating flows (e.g., Farrell and Ioannou 2003, 2007, Srinivasan and Young 2012, Bakas and Ioannou 2013, Tobias and Marston 2013), without relying on eddy–eddy interactions and inverse cascades (Rhines 1975, Vallis and Maltrud 1993). QL approximations of large-scale atmospheric dynamics, sometimes with added damping and stochastic forcing, were partially successful in reproducing aspects of the atmospheric climate and its variability (e.g., Whitaker and Sardeshmukh 1998, Zhang and Held 1999, DelSole 2001, O'Gorman and Schneider 2007). At smaller scales, QL approximations also capture sheared stably stratified flows when the dynamics involve the linear excitation and absorption of internal gravity waves (Orr 1907, Lindzen 1988, Bakas and Ioannou 2007). They also reproduce aspects of thermal convection, such as the dependence of the heat flux on the Rayleigh number (e.g., Malkus 1954, Herring 1963, Toomre et al 1977, Busse 1978, Niemela et al 2000).

In what follows, we derive the CE2 closure for Boussinesq flows, present fully nonlinear and QL simulations of a dry convective boundary layer using large-eddy simulations (LES), and study the evolution of a large-scale wave disturbance on a zonal jet representative of the upper troposphere. The results will demonstrate the potential and limitations of CE2 approaches.

2. CE2 of Boussinesq flow

2.1. Boussinesq flow

Atmospheric flows have low Mach number, so sound waves are generally unimportant for the dynamics. It is therefore common to study atmospheric dynamics with approximations to the Navier–Stokes equations that filter sound waves. The simplest such approximation, which ignores all density variations except where they affect the buoyancy of air masses, is the Boussinesq approximation (Boussinesq 1872). The Boussinesq equations are obtained by expressing density $\rho ({\bf{r}},t)={\rho }_{0}+\delta \rho ({\bf{r}},t)$ as a sum of a constant density ${\rho }_{0}$ in a reference state and fluctuations $\delta \rho ({\bf{r}},t)$ about it, assuming pressure variations in the reference state are hydrostatically balanced, and retaining only the leading-order terms in density and pressure fluctuations in an expansion of the Navier–Stokes equations. The resulting equations are (e.g., Vallis 2006)

$\begin{eqnarray}&&\displaystyle \frac{\partial {\bf{u}}}{\partial t}+({\bf{u}}\cdot {\rm{\nabla }}){\bf{u}}+2{\boldsymbol{\Omega }}\times {\bf{u}}=-{\rm{\nabla }}{\rm{\Phi }}+b{\bf{k}}+{{ \mathcal F }}_{u}\qquad {\rm{(momentum}}\;{\rm{equation)}},\end{eqnarray} \tag{ 1a }$

$\begin{eqnarray}&&\quad {\rm{\nabla }}\cdot {\bf{u}}=0\qquad {\rm{(continuity}}\;{\rm{equation)}},\end{eqnarray} \tag{ 1b }$

$\begin{eqnarray}&&\quad \displaystyle \frac{\partial b}{\partial t}+{\bf{u}}\cdot {\rm{\nabla }}b={{ \mathcal F }}_{b}\qquad {\rm{(thermodynamic}}\;{\rm{equation)}}.\end{eqnarray} \tag{ 1c }$

Here, ${\bf{u}}$ denotes the three-dimensional velocity, $\delta p$ the pressure perturbation associated with the density perturbation $\delta \rho$ , ${\rm{\Phi }}=\delta p/{\rho }_{0}$ the potential of the pressure-gradient accelerations, and $b=-g\delta \rho /{\rho }_{0}$ the buoyancy (g is an effective gravitational acceleration and ${\bf{k}}$ is the local vertical). The reference frame rotates with the constant angular velocity ${\boldsymbol{\Omega }}$ of the planetary rotation, as a result of which Coriolis accelerations ( $2{\boldsymbol{\Omega }}\times {\bf{u}}$ ) appear in the momentum equation. Centrifugal accelerations are subsumed in g, the effective gravitational acceleration. The terms ${{ \mathcal F }}_{{\bf{u}}}$ and ${{ \mathcal F }}_{b}$ on the right-hand side represent dissipation and forcing terms (for example, friction and diabatic heating). The continuity equation reduces to the condition that the flow ${\bf{u}}$ is nondivergent. Sound waves are filtered from these equations because no time derivative appears in the continuity equation. In effect, the speed of sound is taken to be infinite, so that pressure adjusts instantaneously across the flow domain: it can be determined from a Poisson equation obtained by taking the divergence of the momentum equation and using the nondivergence condition to eliminate the time derivative.

To write the equations in a synthetic way, we introduce the state vector

$\begin{eqnarray}&&{\boldsymbol{\Psi }}={(u,v,w,b)}^{{\rm{T}}}\end{eqnarray} \tag{ 2 }$

of the flow that contains all prognostic variables in Cartesian coordinates (the superscript ${(\cdot )}^{{\rm{T}}}$ indicates the transpose). The set of equations (1) can then be written compactly as

$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{\Psi }}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{\Psi }}\otimes {\bf{u}})={\bf{L}}\;{\boldsymbol{\Psi }}-{\rm{\nabla }}{\rm{\Phi }}+{\bf{F}}\end{eqnarray} \tag{ 3a }$

$\begin{eqnarray}&&\quad {\rm{\nabla }}\cdot {\bf{u}}=0,\end{eqnarray} \tag{ 3b }$

where the outer product ( $\otimes$ ) of the two vectors ${\boldsymbol{\Psi }}$ and ${\bf{u}}$ is defined as

$\begin{eqnarray}{\boldsymbol{\Psi }}\otimes {\bf{u}}=\left[\begin{array}{ccc}{{\rm{\Psi }}}_{1}{u}_{1} & \cdots & {{\rm{\Psi }}}_{1}{u}_{3}\\ \vdots & \ddots & \vdots \\ {{\rm{\Psi }}}_{4}{u}_{1} & \cdots & {{\rm{\Psi }}}_{4}{u}_{3}\end{array}\right]={({{\rm{\Psi }}}_{i}{u}_{j})}_{\displaystyle \genfrac{}{}{0em}{}{1\leqslant i\leqslant 4}{1\leqslant j\leqslant 3}}.\end{eqnarray} \tag{ 4 }$

The components of the divergence of the second-order term ${\rm{\nabla }}\cdot ({\boldsymbol{\Psi }}\otimes {\bf{u}})$ are

$\begin{eqnarray}&&{[{\rm{\nabla }}\cdot ({\boldsymbol{\Psi }}\otimes {\bf{u}})]}_{i}=\displaystyle \frac{\partial }{\partial {r}_{j}}[{{\rm{\Psi }}}_{i}{u}_{j}],\end{eqnarray} \tag{ 5 }$

where summation over repeated indices is implied. This notation extends the usual advection operator of a scalar field by a nondivergent flow to the advection of a vector field. The linear operator ${\bf{L}}$ contains accelerations owing to the Coriolis force, buoyancy, as well as nonconservative terms (e.g., friction or diabatic heating) that are linear in the state vector. The vector ${\bf{F}}$ contains all other nonconservative terms. We have expanded ${\bf{u}}$ , the outer product and the ${\rm{\nabla }}\cdot$ in Cartesian coordinates for simplicity, however the formalism is coordinate independent.

Despite their relative simplicity, the Boussinesq equations are commonly used to study turbulence in boundary layers, in which density variations are weak. They also underlie classical conceptual models of large-scale atmospheric dynamics (e.g., quasigeostrophic models), in which they become quantitatively inaccurate but still qualitatively capture important atmospheric phenomena such as Rossby waves and baroclinic instability. The Boussinesq equations are well suited for our exposition of CE2 approaches because they capture the essential nonlinearity of atmospheric flows: the conservative quadratic nonlinearity of the advection operators $({\bf{u}}\cdot {\rm{\nabla }}){\bf{u}}$ and ${\bf{u}}\cdot {\rm{\nabla }}b$ .

2.2. Averaging operator

Our interest is not in individual details of the atmospheric flows under consideration but in their statistics, including mean values and higher moments. Therefore, we define an averaging operator, denoted with an overline $\bar{(\cdot )}$ , and decompose any scalar field $f({\bf{r}},t)$ into a mean and a fluctuating part

$\begin{eqnarray}&&f({\bf{r}},t)=\bar{f}({\bf{r}},t)+{f}^{\prime }({\bf{r}},t).\end{eqnarray} \tag{ 6 }$

The mean is in general a function of space and time. The fluctuating part is commonly referred to as an eddy. The averaging operator is taken to satisfy, for all scalar fields $f({\bf{r}},t)$ and $g({\bf{r}},t)$ and any constant c, the Reynolds properties (Monin and Yaglom 1971):

$\begin{eqnarray}&&\bar{c}=c,\end{eqnarray} \tag{ 7a }$

$\begin{eqnarray}&&\overline{{cf}+g}=c\bar{f}+\bar{g}\qquad \qquad {\rm{(linearity)}},\end{eqnarray} \tag{ 7b }$

$\begin{eqnarray}&&\bar{\bar{f}g}=\bar{f}\bar{g},\end{eqnarray} \tag{ 7c }$

$\begin{eqnarray}&&\bar{\partial f}=\partial \bar{f}\qquad {\rm{(commutation}}\;{\rm{with}}\;{\rm{derivatives)}}.\end{eqnarray} \tag{ 7d }$

Properties (7a)–(7c) imply that the averaging operator is a projection operator and so is idempotent $\bar{\bar{f}}=\bar{f}$ . They make it possible to define the Reynolds decomposition

$\begin{eqnarray}&&\bar{{fg}}=\bar{f}\bar{g}+\bar{{f}^{\prime }{g}^{\prime }}.\end{eqnarray} \tag{ 8 }$

For a vector quantity ${\boldsymbol{\Psi }}$ , the average $\bar{{\boldsymbol{\Psi }}}$ is the component-wise average.

The choice of average is unspecified as long as it satisfies the Reynolds properties. Conceptually, ensemble averages are statistically meaningful, and they naturally satisfy the Reynolds properties. In practice, however, they can be difficult to obtain. Averages over sufficiently long times in statistically stationary (or slowly varying) flows or over sufficiently large regions in statistically homogeneous flows are more commonly used in practice, and also approximately satisfy the Reynolds properties. In concrete calculations, we choose the averages that are natural given the statistical symmetries of the problem under consideration. For example, in flows that are statistically invariant under translations along a spatial coordinate (e.g., along latitude circles), an average along that spatial coordinate suggests itself.

For more general flow equations with a variable density, the averaging operator above has to be replaced by a density-weighted average, to obtain consistent equations of motion for the statistical moments that resemble the Boussinesq moment equations formally. An example for the anelastic approximation is provided in appendix A.

2.3. Cumulant expansion

First cumulant

The first cumulant is the mean $\bar{{\boldsymbol{\Psi }}}({\bf{r}},t)$ , for which the equations of motion are obtained by averaging the equations of motion (3):

$\begin{eqnarray}&&\displaystyle \frac{\partial \bar{{\boldsymbol{\Psi }}}}{\partial t}+{\rm{\nabla }}\cdot (\bar{{\boldsymbol{\Psi }}}\otimes \bar{{\bf{u}}})=-{\rm{\nabla }}\cdot (\bar{{{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime }})+{\bf{L}}\;\bar{{\boldsymbol{\Psi }}}-{\rm{\nabla }}\bar{\phi }+\bar{{\bf{F}}},\end{eqnarray} \tag{ 9a }$

$\begin{eqnarray}&&{\rm{\nabla }}\cdot \bar{{\bf{u}}}=0.\end{eqnarray} \tag{ 9b }$

This involves the covariance ${\rm{\nabla }}\cdot (\bar{{{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime }})$ , which arises from the quadratic nonlinearity of the equations of motion. It represents eddy fluxes, for example, arising from advection of momentum fluctuations by the eddies (fluctuations) themselves. Because the equation for the mean involves a covariance, it is not closed.

Second cumulant

The second cumulant is the second central moment, or the covariance

$\begin{eqnarray}&&{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\bar{{{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{1},t)\otimes {{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{2},t)}.\end{eqnarray} \tag{ 10 }$

We only consider the prognostic variables here, because the diagnostic variables (and hence their statistics) can be obtained from them. We also only consider equal-time cumulants, that is, equal-time covariances between prognostic variables at the two points ${{\bf{r}}}_{1}$ and ${{\bf{r}}}_{2}$ , which need not be equal. The covariance tensor is symmetric

$\begin{eqnarray}&&{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)={{\bf{C}}}^{{\rm{T}}}({{\bf{r}}}_{2},{{\bf{r}}}_{1},t).\end{eqnarray} \tag{ 11 }$

We additionally define the auxiliary covariances

$\begin{eqnarray}&&{{\bf{C}}}^{{\rm{\Phi }}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\bar{{{\rm{\Phi }}}^{\prime }({{\bf{r}}}_{1},t){{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{2},t)},\end{eqnarray} \tag{ 12a }$

$\begin{eqnarray}&&{{\bf{C}}}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\bar{{{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{1},t)\otimes {{\bf{u}}}^{\prime }({{\bf{r}}}_{2},t)}.\end{eqnarray} \tag{ 12b }$

The first, ${{\bf{C}}}^{{\rm{\Phi }}}$ , contains additional information about covariation of the prognostic fields ${\boldsymbol{\Psi }}$ with the pressure potential Φ. These covariances can be calculated from ${\bf{C}}$ and $\bar{{\rm{\Psi }}}$ because the pressure potential Φ is a diagnostic variable in the Boussinesq approximation. The second, ${{\bf{C}}}^{u}$ , represents covariation of the velocity field with other prognostic variables and is already contained in ${\bf{C}}$ .

The second cumulant equation can be obtained from the equations of motion (3) by evaluating

$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\bar{{{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{1},t)\otimes \displaystyle \frac{\partial {{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{2},t)}{\partial t}}+\bar{\displaystyle \frac{\partial {{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{1},t)}{\partial t}\otimes {{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{2},t)}.\end{eqnarray} \tag{ 13 }$

Discarding terms that are third order in fluctuating quantities, one obtains

$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+{{\rm{\nabla }}}_{{{\bf{r}}}_{1}}\cdot [{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\otimes \bar{{\bf{u}}}({{\bf{r}}}_{1},t)]=-{{\bf{C}}}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t){[{\rm{\nabla }}\bar{{\rm{\Psi }}}({{\bf{r}}}_{2},t)]}^{{\rm{T}}}\\ &&\qquad +{\bf{L}}\;{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+{{\rm{\nabla }}}_{{{\bf{r}}}_{1}}{{\bf{C}}}^{{\boldsymbol{\phi }}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+\overline{{\bf{F}}^{\prime} ({{\bf{r}}}_{1},{\bf{t}})\otimes {\boldsymbol{\Psi }}^{\prime} ({{\bf{r}}}_{2},{\bf{t}})}+\{{{\bf{r}}}_{1}\;\longleftrightarrow \;{{\bf{r}}}_{2}\},\end{eqnarray} \tag{ 14 }$

where $\{{{\bf{r}}}_{1}\;\longleftrightarrow \;{{\bf{r}}}_{2}\}$ indicates the terms obtained by interchanging ${{\bf{r}}}_{1}$ and ${{\bf{r}}}_{2}$ , which are necessary to ensure the symmetry (11) of the covariance tensor. The third-order term ${\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\otimes \bar{{\bf{u}}}({{\bf{r}}}_{1},t)$ is defined as in Cartesian coordinates

$\begin{eqnarray}&&{[{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\otimes \bar{{\bf{u}}}({{\bf{r}}}_{1},t)]}_{p,q,i}={C}_{p,q}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t){\bar{u}}_{i}({{\bf{r}}}_{1},t),\end{eqnarray} \tag{ 15 }$

and its divergence by

$\begin{eqnarray}&&{\{{{\rm{\nabla }}}_{{{\bf{r}}}_{1}}\cdot [{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\otimes \bar{{\bf{u}}}({{\bf{r}}}_{1},t)]\}}_{p,q}=\displaystyle \frac{\partial }{\partial {r}_{1i}}[{C}_{p,q}{\bar{u}}_{i}].\end{eqnarray} \tag{ 16 }$

The flux ${\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\otimes \bar{{\bf{u}}}({{\bf{r}}}_{1},t)$ represents the transport of spatial eddy correlations by the mean flow at ${{\bf{r}}}_{1}$ . The term $-{{\bf{C}}}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\cdot {\rm{\nabla }}{[\bar{{\rm{\Psi }}}({{\bf{r}}}_{2},t)]}^{{\rm{T}}}$ represents generation of covariance ${\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)$ by advection down mean-flow gradients.

Additionally, continuity (1b) implies that

$\begin{eqnarray}&&{{\rm{\nabla }}}_{{{\bf{r}}}_{2}}\cdot {{\bf{C}}}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)={{\rm{\nabla }}}_{{{\bf{r}}}_{1}}\cdot {{\bf{C}}}^{u}({{\bf{r}}}_{2},{{\bf{r}}}_{1},t)=0.\end{eqnarray} \tag{ 17 }$

The set of equations for the first and second cumulants in this form is closed because the third cumulants, which would ordinarily appear in the second cumulant equations owing to the quadratic nonlinearity of the equations of motion, have been discarded. This second-order truncation of the otherwise infinite hierarchy of cumulant equations is referred to as a second-order CE2. In that CE2 assumes the first and second cumulants suffice to describe the flow statistics, it makes a normal approximation to the hierarchy of equations for the flow statistics.

CE2 equations

To summarize, the CE2 equations, with the second cumulant substituted in (9a), are given by

$\begin{eqnarray}&&\displaystyle \frac{\partial \bar{{\boldsymbol{\Psi }}}({\bf{r}},t)}{\partial t}+{\rm{\nabla }}\cdot [\bar{{\boldsymbol{\Psi }}}({\bf{r}},t)\otimes \bar{{\bf{u}}}({\bf{r}},t)]=-{\rm{\nabla }}\cdot {{\bf{C}}}^{u}({\bf{r}},{\bf{r}},t)\\ &&\qquad \quad +{\bf{L}}\;\bar{{\boldsymbol{\Psi }}}({\bf{r}},t)-{\rm{\nabla }}\bar{{\rm{\Phi }}}({\bf{r}},t)+\bar{{\bf{F}}}({\bf{r}},t),\end{eqnarray} \tag{ 18a }$

$\begin{eqnarray}&&{\rm{\nabla }}\cdot \bar{{\bf{u}}}=0,\end{eqnarray} \tag{ 18b }$

$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+{{\rm{\nabla }}}_{{{\bf{r}}}_{1}}\cdot [{\bf{C}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\otimes \bar{{\bf{u}}}({{\bf{r}}}_{1},t)]=-{{\bf{C}}}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t){[{\rm{\nabla }}\bar{{\rm{\Psi }}}({{\bf{r}}}_{2},t)]}^{{\rm{T}}}\\ &&\qquad +{\bf{LC}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+{{\rm{\nabla }}}_{{{\bf{r}}}_{1}}{{\bf{C}}}^{{\boldsymbol{\phi }}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+\bar{{{\bf{F}}}^{\prime }({{\bf{r}}}_{1},{\bf{t}})\otimes {{\boldsymbol{\Psi }}}^{\prime }({{\bf{r}}}_{2},{\bf{t}})}+\{{{\bf{r}}}_{1}\;\longleftrightarrow \;{{\bf{r}}}_{2}\},\end{eqnarray} \tag{ 18c }$

$\begin{eqnarray}&&{{\rm{\nabla }}}_{{{\bf{r}}}_{2}}\cdot {{\bf{C}}}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)={{\rm{\nabla }}}_{{{\bf{r}}}_{1}}\cdot {{\bf{C}}}^{u}({{\bf{r}}}_{2},{{\bf{r}}}_{1},t)=0.\end{eqnarray} \tag{ 18d }$

This set of equations involves the 16 terms in ${\bf{C}}$ and the eight covariances ${{\bf{C}}}^{{\rm{\Phi }}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)$ and ${{\bf{C}}}^{{\rm{\Phi }}}({{\bf{r}}}_{2},{{\bf{r}}}_{1},t)$ . The 16 components of the second cumulant ${\bf{C}}$ are prognostic (evoloving according equation (18c)), and the other 8 covariances components involve diagnostic correlations between the pressure potential and each of the other fields. The 8 corresponding diagnostic Poisson equations are obtained by taking the divergence with respect to ${{\bf{r}}}_{1}$ or ${{\bf{r}}}_{2}$ of the equation for ${{\bf{C}}}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)$ that can be extracted from (18c), and in which time tendencies vanish because of the nondivergence constraint (18d).

Physical properties

CE2 is a realizable approximation, in the sense that the implied probability distribution functions are positive (Marston et al 2014). The first cumulant equations (18a), (18b) are unchanged from the fully nonlinear system. Therefore, first-order invariants (mass and momentum) are conserved by the CE2 system in the absence of nonconservative effects. The third cumulants, which are neglected in the second cumulant equations (18c), (18d), redistribute second-order invariants (e.g., energy) among scales but do not generate or dissipate these invariants. Therefore, second-order invariants such as energy are likewise conserved by the CE2 system in the absence of nonconservative effects (see Marston et al 2014).

2.4. QL approximation

The CE2 equations can also be obtained by directly approximating the original equations of motion (3), making what has come to be called the QL approximation (O'Gorman and Schneider 2007, Srinivasan and Young 2012, Constantinou et al 2014b, Ait-Chaalal and Schneider 2015). The QL approximation keeps nonlinear terms in the equation (9) for the mean $\bar{{\boldsymbol{\Psi }}}$ . However it linearizes the equation for the eddies ${{\boldsymbol{\Psi }}}^{\prime }$ , obtained by substracting (9a) from (3a),

$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{\Psi }}}^{\prime }}{\partial t}+{\rm{\nabla }}\cdot (\bar{{\boldsymbol{\Psi }}}\otimes {{\bf{u}}}^{\prime })+{\rm{\nabla }}\cdot ({{\boldsymbol{\Psi }}}^{\prime }\otimes \bar{{\bf{u}}})+[{\rm{\nabla }}\cdot ({{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime })-{\rm{\nabla }}\cdot (\bar{{{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime }})]={\bf{L}}\;{{\boldsymbol{\Psi }}}^{\prime }-{\rm{\nabla }}{\rm{\Phi }}+{{\bf{F}}}^{\prime },\end{eqnarray} \tag{ 19 }$

by setting

$\begin{eqnarray}&&{\rm{\nabla }}\cdot ({{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime })={\rm{\nabla }}\cdot (\bar{{{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime }}).\end{eqnarray} \tag{ 20 }$

Hence, the QL approximation amounts to replacing in the Reynolds decomposition of the nonlinear term

$\begin{eqnarray}&&{\boldsymbol{\Psi }}\otimes {\bf{u}}=\bar{{\boldsymbol{\Psi }}}\otimes \bar{{\bf{u}}}+\bar{{\boldsymbol{\Psi }}}\otimes {{\bf{u}}}^{\prime }+{{\boldsymbol{\Psi }}}^{\prime }\otimes \bar{{\bf{u}}}+{{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime },\end{eqnarray} \tag{ 21 }$

the eddy–eddy interaction ${{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime }$ by its mean effect $\bar{{{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime }}$ . Under the QL approximation, the Boussinesq equations (3) can then be written as

$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{\Psi }}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{\Psi }}\otimes {\bf{u}})+{\rm{\nabla }}\cdot (\bar{{{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime }}-{{\boldsymbol{\Psi }}}^{\prime }\otimes {{\bf{u}}}^{\prime })={\bf{L}}\;{\boldsymbol{\Psi }}-{\rm{\nabla }}{\rm{\Phi }}+{\bf{F}},\end{eqnarray} \tag{ 22a }$

$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\bf{u}}=0.\end{eqnarray} \tag{ 22b }$

Because the QL equations retain as the only nonlinear interaction the interaction between eddies and mean fields, the corresponding cumulant equations are closed at second order, meaning no third-order terms appear in the second-order equations. The first two cumulant equations are exactly the CE2 equations (18). This makes it possible to simulate flows whose statistics satisfy the CE2 equations (18) simply by simulating the QL equations (22).

The QL truncation does not necessarily imply that eddy–eddy interactions and third-order correlations are equal to zero. However their evolution is decoupled from that of lower-order cumulants. This has to be kept in mind when interpreting instantaneous fields and statistics of flows simulated by the QL equations.

2.5. Differences between CE2 and S3T

S3T is a second-order statistical approach to turbulent flows that is closely related to CE2. CE2 and S3T are sometimes presented as being equivalent in the literature because they share a similar mathematical formalism.

However, CE2 and S3T differ in that S3T includes a small-scale stochastic forcing that is white in time and represents the scattering by missing eddy–eddy interactions. The resulting additional energy injection is balanced by large-scale linear damping. The stochastic forcing allows one to define rigorously an ensemble average over its realizations and permits a semi-analytical treatment of the second-order equations, whose solutions depend on the existence and statistics of the stochastic forcing.

By contrast, CE2 uses the same forcing in the truncated equations as in the fully nonlinear equations, without attempting to parameterize eddy–eddy interactions. This choice is made for two reasons. First, a stochastic noise is not necessarily a realistic model for eddy–eddy interactions because these interactions do not inject energy but only redistribute it among scales. Second, forcing the flow at small scales might appear unnatural for a wide range of planetary flows, which are forced on larger scales. For example, large-scale motion in Earth's atmosphere is essentially driven by the planetary-scale radiative imbalance between the equator and the poles, rather than by energy injection at small scales.

Another difference is that S3T is generally applied to flows for which the linear operator representing eddy–mean flow interactions does not have unstable modes (the mathematical development of S3T relies on non-normal growth and decay of stable modes, excited by the stochastic noise). In this framework, S3T has been used in barotropic rotating flows to study instabilities of zonal flows (Bakas and Ioannou 2013, Parker and Krommes 2014) and the dynamics of zonal (Farrell and Ioannou 2003, 2007, Constantinou et al 2014a) and nonzonal (Bakas and Ioannou 2014) coherent structures. Nevertheless, the simulation of unstable flows at second order is possible and well-defined mathematically. Hence, CE2 (or its QL analogue) has been applied to unstable flows in barotropic (Marston et al 2008) and baroclinic (O'Gorman and Schneider 2007, Ait-Chaalal and Schneider 2015) settings. Evaluating the ability of CE2 in capturing unstable flows is at the cornerstone of the present paper.

3. Dry convective boundary layer

The dynamics of boundary layers and clouds involve flow scales of order meters to a few kilometers. In addition, condensate formation and evaporation take place at microscales. By contrast, typical climate models have a horizontal resolution of order 100 km. Therefore, the dynamics of boundary layers and clouds are subgrid-scale (SGS) processes in climate models, which must be represented parametrically in terms of the resolved large-scale dynamics (e.g., Beljaars 1992, Smith 1997). Uncertainties about these parameterization schemes are the dominant contributor to uncertainties in climate change projections (e.g., Stephens 2005, Bony et al 2006, Soden and Held 2006, Stevens and Bony 2013, Webb et al 2013, Vial et al 2013, Brient et al 2015).

Most current parameterization schemes for the dynamics of clouds and boundary layers truncate the hierarchy of moment (or cumulant) equations at first order and represent the second-order SGS fluxes appearing in the first-order equations semi-empirically. For example, turbulent fluxes in boundary layers are often represented as diffusive fluxes down mean gradients, with diffusivities estimated from approximate spatially local second-order equations for turbulence kinetic energy (TKE) (e.g., Mellor and Yamada 1982, Beljaars 1992). Turbulent fluxes in convective clouds, on the other hand, are often represented as vertically nonlocal entraining plumes that extend deep into the boundary layer. The parameters determining the mass fluxes in the plumes are estimated from energy equations (e.g., Arakawa and Schubert 1974, Gregory 1997).

CE2 may offer a path toward improved and unified representations of SGS turbulent fluxes (e.g. in boundary layers and in convective clouds) in climate models, because CE2 explicitly retains spatial nonlocality and interactions between fluctuations (plumes or eddies) and the environment (mean field). Here we compare a QL simulation and a fully nonlinear simulation of the simplest convective boundary layer, a dry convective boundary layer, and demonstrate the potential and limitations of CE2 to represent the statistics of its dynamics. We conducted a large-eddy simulation (LES) of a dry convective boundary layer growing into a stable background stratification as a heat flux is imposed at the surface (Soares 2004). We then compared this fully nonlinear simulation with a simulation in which the equations were replaced by the corresponding QL equations (22).

3.1. Large-eddy simulations

Setup

The LES code solves the anelastic equations and is described in Pressel et al (2015). Like the Boussinesq approximation, the anelastic equations are nonhydrostatic and filter sound waves by neglecting dynamic density variations except where they affect buoyancy. But in contrast to the Boussinesq approximation, the background state depends on the vertical coordinate. Hence, the anelastic equations can capture the dynamics of flows with substantial stratification and so are better suited to study atmospheric convection. The anelastic approximation and its CE2 are described in appendix A.

Our anelastic LES code uses the specific entropy s and the three-dimensional velocity field ${\bf{u}}$ as prognostic variables (Pauluis 2008). We use a second-order central difference spatial discretization scheme with strongly stability preserving Runge–Kutta time-stepping (Spiteri and Ruuth 2002). The time step is dynamically adjusted to ensure a Courant number near 0.3. Because LES merely resolves the most energetic scales of the flow, SGS motions must also be modeled, and we do so by applying a constant eddy diffusivity of $\nu =1.2\;{{\rm{m}}}^{2}\;{{\rm{s}}}^{-1}$ throughout the domain. We choose a constant diffusivity to avoid the nonlinearities that appear in the computation of the eddy viscosity by more advanced SGS schemes such as the Smagorinsky–Lilly closure (Lilly 1962, Smagorinsky 1963), which would need to be linearized in a QL simulation; constant diffusion as an SGS closure allows a more direct comparison of fully nonlinear and QL simulations, notwithstanding that it is an inferior SGS closure. The domain extends $6.4\;\mathrm{km}\times 6.4\;\mathrm{km}$ in the horizontal and $3.75\;\mathrm{km}$ in the vertical, with a horizontal and vertical resolution of $25\;{\rm{m}}$ . Horizontal boundary conditions are doubly periodic. At the upper boundary, flow fields are linearly relaxed over a $800\;{\rm{m}}$ deep layer toward the horizontal mean flow, which is almost motionless (horizontal mean velocities are of the order $0.1\;{\rm{m}}\;{{\rm{s}}}^{-1}$ ). The relaxation coefficient varies from $\tau =0$ at the bottom of the layer to $\tau ={(100{\rm{s}})}^{-1}$ at the top.

The initial state is stably stratified with a potential temperature θ that increases linearly with height, at a rate of $2\;\;{\rm{K}}\;{\mathrm{km}}^{-1}$ . Here, the potential temperature is related to the specific entropy by $\theta =\tilde{T}\mathrm{exp}[(s-\tilde{s})/{c}_{{\rm{p}}}]$ , where c_p is the specific heat at constant pressure, $\tilde{T}$ is a standard temperature, and $\tilde{s}$ is a standard specific entropy at the standard temperature $\tilde{T}$ and standard pressure⁴ $\tilde{p}=1000\;\mathrm{hPa}$ . The initial state is destabilized by imposing a constant sensible heat (enthalpy) flux of $70.46\;{\mathrm{Wm}}^{-2}$ at the lower boundary. Normally distributed random fluctuations of the potential temperature with amplitude $0.1\;{\rm{K}}$ in the lowest $200\;{\rm{m}}$ break the horizontal homogeneity of the initial state and allow the generation of turbulent motions. The initial flow is uniform, horizontal, and has a speed of $0.01\;{\rm{m}}\;{{\rm{s}}}^{-1}$ . Together with the drag at the lower boundary, this allows for turbulent momentum fluxes to develop (Soares 2004). We ran simulations for up to 12 simulated hours, over which a dry convective boundary layer forms and grows as a result of the heating at the bottom.

Because of the statistical horizontal homogeneity of the flow, we use the horizontal average to define mean fields and eddies, as in equation (6). The QL truncation (22), here with the state vector ${\boldsymbol{\Psi }}={(u,v,w,s)}^{{\rm{T}}}$ , is implemented by removing at every time step the nonlinear eddy–eddy interactions from the tendencies of all prognostic variables⁵ . Herring (1963) similarly used the QL approximation to study thermal convection between two parallel horizontal plates.

3.2. Results

Figure 1 compares vertical and horizontal cross sections of the vertical velocity field in the fully nonlinear and in the QL simulations. The vertical cross sections (figures 1(a), (b)) show that upward motion mainly occurs in vertically coherent updrafts, as is well known (e.g., Schmidt and Schumann 1989). In the QL simulation, updrafts are more coherent because small-scale structures are missing while the larger scales are well represented. The horizontal cross section in the fully nonlinear simulation reveals that the updrafts are organized into polygonal convective cells (figure 1(c)). Such cellular patterns are well known to arise near the onset of thermal instability of a fluid heated from below (Chandrasekhar 1961, chapter 2). The QL simulation does not capture these horizontal correlation patterns (figure 1(d)), probably because eddy–eddy interactions in the horizontal plane play a role in generating the smaller scales of the cellular patterns. The vertical velocity fluctuations ${w}^{\prime }$ are distributed more symmetrically around zero in the QL simulation: updrafts and downdrafts are of similar size and strength (figures 1(b), (d)). This stands in contrast to the fully nonlinear simulations, in which updrafts are faster and narrower and downdrafts slower and broader.

Figure 2(a) shows the evolution of the vertical profile of the potential temperature in the full and QL simulations. Because of the constant heating at the lower boundary and the absence of any thermal relaxation in the fluid interior, the boundary layer continually deepens and does not reach a statistically steady state (figure 2(a)). The boundary layer is well mixed with homogenized potential temperature below its top. This indicates that the BL is close to the critical state for thermal instability. The QL simulation captures quite accurately the growth rate of the boundary layer and the mixing of potential temperature below its top. This is because the growth of the boundary layer mainly arises through interactions of fluctuations with the mean flow, which are retained in the QL simulation.

**Figure 2.** (a) Vertical profiles of potential temperature at indicated times in the fully nonlinear (solid lines) and QL (dashed lines) simulations. (b) Corresponding vertical profiles of the potential temperature variance $\bar{{{\theta }^{\prime }}^{2}}$ .
Download figure:
Standard image High-resolution image

Above the top of the boundary layer, defined as the altitude below which potential temperature is well mixed, the QL and the fully nonlinear simulations exhibit important differences. In the QL simulation, the mean potential temperature profile is identical to the initial profile. By contrast, the fully nonlinear simulation shows a layer of strong stability associated with convective overshoots of thermal updrafts into the free atmosphere and the downward entrainment of warmer free-atmospheric air into the boundary layer (e.g., de Roode et al 2004). These overshoots are missing in the QL simulation, as they are in first-order diffusive closure schemes for convective boundary layers (e.g., Stull 1988). The difference at the top of the boundary layer is emphasized in the profile of temperature variance (figure 2(b)), which shows a strong peak above the top of the boundary layer for the fully nonlinear simulation but vanishing variance throughout the mixed layer and aloft for the QL simulation. Near the surface, the mean potential temperature is reduced in the QL simulation, which reflects a more efficient upward transport of heat into the interior of the boundary layer. This probably can be ascribed to the more intense vertical velocities and the strengthened correlations between buoyancy and vertical velocity fluctuations in the more coherent updrafts of the QL flow.

An inability of diffusive closures and the QL simulation to represent the turbulent transport of turbulence kinetic energy (TKE) lies behind the failure of the QL simulation and diffusive closures to represent convective overshoots. We take TKE to be the kinetic energy of the resolved scales of the LES: $e=0.5\;({{u}^{\prime }}^{2}+{{v}^{\prime }}^{2}+{{w}^{\prime }}^{2})$ . Its mean budget is derived from the momentum equations and reads

$\begin{eqnarray}{\partial }_{t}\bar{e} & = & -\underset{{ \mathcal A }}{\underbrace{\bar{w}\;{\partial }_{z}\bar{e}}}-\underset{{ \mathcal T }}{\underbrace{\displaystyle \frac{1}{{\rho }_{0}}{\partial }_{z}({\rho }_{0}\overline{w^{\prime} e})}}-\underset{{ \mathcal S }}{\underbrace{[\overline{w^{\prime} u^{\prime} }{\partial }_{z}\bar{u}+\overline{w^{\prime} v^{\prime} }{\partial }_{z}\bar{v}+\overline{w^{\prime} w^{\prime} }{\partial }_{z}\bar{w}]}}\\ & & +\underset{{ \mathcal B }}{\underbrace{\overline{w^{\prime} b^{\prime} }}}-\underset{{ \mathcal P }}{\underbrace{\displaystyle \frac{1}{{\rho }_{0}}{\partial }_{z}(\overline{w^{\prime} p^{\prime} })}}+\epsilon .\end{eqnarray} \tag{ 23 }$

The right-hand side contains all terms that generate, destroy, or redistribute TKE: advection by the mean flow ${ \mathcal A }$ and the eddies ${ \mathcal T }$ , shear production ${ \mathcal S }$ , buoyancy production ${ \mathcal B }$ ( ${b}^{\prime }=g{\theta }^{\prime }/{\theta }_{0}$ denotes buoyancy fluctuations), the pressure correlation term ${ \mathcal P }$ and dissipation . The dissipation denotes the energy flux from the resolved scales to sub-grid scales, which in our case is parameterized as viscous dissipation of the resolved fields with constant eddy viscosity ν. All but the turbulent transport term ${ \mathcal T }=-1/{\rho }_{0}{\partial }_{z}({\rho }_{0}\bar{{w}^{\prime }e})$ are of second order in fluctuations and hence are retained in a QL truncation.

The TKE budgets for the fully nonlinear and QL simulations are shown in figure 3. The dominance of the buoyancy production term relative to the shear production term in both plots indicates that the flow is thermally driven. In the fully nonlinear simulation, buoyancy production is positive throughout the well-mixed part of the boundary layer (i.e., upward buoyancy flux), zero at the top of the boundary layer, and negative aloft (i.e., downward buoyancy flux). The negative flux is related to the overshooting thermals, which trigger a downward entrainment flux of free atmospheric air into the boundary layer. This negative flux is missing in the QL simulation, in which the buoyancy flux is zero above the top of the boundary layer. However, the buoyancy flux is well captured in the interior of the boundary layer.

**Figure 3.** Turbulence kinetic energy budget of the fully nonlinear (a) and QL simulations (b). The different terms denote the TKE advection ${ \mathcal A }$ , TKE transport by eddies ${ \mathcal T }$ , shear production ${ \mathcal S }$ , buoyancy production ${ \mathcal B }$ , pressure correlation term ${ \mathcal P }$ , dissipation to subgrid-scales ${ \mathcal D }$ and rate of change of TKE ${ \mathcal R }$ . Also shown are the residuals, which are of the same order as the SGS diffusion term in the fully nonlinear simulation (see Heinze *et al* 2015) but are smaller in the QL simulation, likely because the latter is smoother, thus limiting numerical diffusion.
Download figure:
Standard image High-resolution image

The triple correlation transport term ${ \mathcal T }$ represents the vertical transport of TKE by eddies. In the fully nonlinear simulation, eddies transport TKE from lower levels with high TKE to higher levels with low TKE. This transport seems to be crucial for the growth of TKE in the upper part of the boundary layer and especially across the top of the boundary layer. The neglect of this term in the QL dynamics is responsible for the missing overshoots of thermals across the boundary layer top and the missing associated negative buoyancy flux. The transport term ${ \mathcal T }$ can be evaluated in the QL simulations and actually is nonzero. However, it decouples from the second-order dynamics and so does not affect the TKE budget.

3.3. Implications

These results show that a QL simulation can capture important aspects of the evolution of a dry convective boundary layer, such as its well-mixed nature and its growth over time. They suggest that a corresponding CE2 closure would also capture the relevant first-order statistics and, therefore, has promise as a nonlocal second-order closure in climate models. Deficiencies such as the missing convective overshoots at the top of the boundary layer may be remedied, for example, by adding a linear, diffusive transport term in the prognostic equations of the second-order moments (e.g. Mellor 1973). The resulting parameterization scheme would solve directly for the 1st- and 2nd-order statistics in every grid box of the large-scale model.

The applicability of such a scheme depends on whether the corresponding CE2 equations can be solved efficiently—much more efficiently than an explicit QL simulation can be run. This will require a simplified representation of horizontal covariances, for example, by assuming approximate statistical symmetries such as horizontal homogeneity and isotropy of fluctuations. But the important nonlocal effects of vertical covariances need to be retained.

The eddy–eddy interactions neglected in the QL simulations may be even more critical for moist convection than they are for the dry convective boundary layer (Firl and Randall 2015). Hence, more sophisticated approaches may be needed to expand the field of application of CE2 closure schemes to moist boundary layer dynamics. Exploring to what degree the CE2 approximation and its extensions can capture the dynamics of clouds and moist boundary layers promises to be a fruitful area of study.

4. Large-scale eddy decay on the rotating sphere

While the preceding example concerned the applicability of CE2 to atmospheric dynamics on scales of meters to kilometers, we now turn to a prototype problem for atmospheric dynamics on scales of hundreds to thousands of kilometers. Eddies on such large scales are essentially the well-known weather systems. They are generated by baroclinic instability and are fundamental for maintaining Earth's climate because they are responsible for the bulk of the atmospheric transport of energy, water vapor, and angular momentum. Through these transports, they shape the distribution of temperature, precipitation, and winds at the surface (Peixoto and Oort 1992). The fundamental balances governing such large-scale eddies are different than those in the boundary layer. The Coriolis force due to the planetary rotation and the average stable stratification become of primary importance, leading to flows that are more two-dimensional in character than boundary-layer flows. A two-dimensional (latitude-longitude) model suffices to illustrate some of the issues one encounters if one wants to develop a closure for the large-scale dynamics of the atmosphere based on CE2s.

4.1. Barotropic model for the upper troposphere

Large-scale eddies in Earth's atmosphere are generated near the surface in midlatitudes, propagate upward through the troposphere, and propagate meridionally in the upper troposphere (Simmons and Hoskins 1978, Edmon et al 1980, Held and Hoskins 1985, Thorncroft et al 1993). Their meridional transport and eventual dissipation by wave breaking in latitude bands away from their generation latitudes is what gives rise to their meridional angular momentum transport: large-scale eddies in rapidly rotating atmospheres transport angular momentum from their dissipation latitudes into their generation latitudes, that is, in the opposite direction to their meridional propagation (Kuo 1951, Held 1975, 1999). This angular momentum transport ultimately shapes the strength and distribution of surface winds, with easterlies in the tropics, westerlies in midlatitudes, and weak easterlies again in polar latitudes (see Schneider 2006, for a review). To understand the strength and distribution of surface winds, it is therefore necessary to understand the meridional propagation and dissipation of large-scale eddies, which are concentrated in the upper troposphere (Ait-Chaalal and Schneider 2015). The simplest model that captures these processes is the barotropic model—a model of a two-dimensional fluid layer on a sphere, thought to represent a layer in the upper troposphere (e.g., Held and Phillips 1987).

Equations of motion

The equations governing barotropic flow can be derived from the Boussinesq equations (1). Consider a Boussinesq flow on a sphere of radius a rotating at constant spin angular velocity ${\boldsymbol{\Omega }}$ . We further assume that the flow is two-dimensional on the sphere: ${\bf{u}}\cdot {{\bf{e}}}_{r}=0$ , with ${{\bf{e}}}_{r}$ being the radial coordinate. In that case, the momentum and continuity equations (1a), (1b) reduce to

$\begin{eqnarray}&&\displaystyle \frac{\partial {\bf{v}}}{\partial t}+{\bf{v}}\cdot {\rm{\nabla }}{\bf{v}}+2{\boldsymbol{\Omega }}\times {\bf{v}}=-{\rm{\nabla }}{\rm{\Phi }}+{{ \mathcal F }}_{{\bf{v}}},\end{eqnarray} \tag{ 24a }$

$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\bf{v}}=0.\end{eqnarray} \tag{ 24b }$

The horizontal components of the wind and of the forcing are denoted ${\bf{v}}$ and ${{ \mathcal F }}_{{\bf{v}}}$ . Because the flow ${\bf{v}}$ is two-dimensional on the sphere, only the local vertical component ${\rm{\Omega }}\mathrm{sin}\phi \;{{\bf{e}}}_{r}$ (latitude ϕ) of ${\boldsymbol{\Omega }}$ yields nonzero terms in (24a). Taking the curl of the momentum equation and projecting it onto the radial direction ${{\bf{e}}}_{r}$ yields the two-dimensional barotropic vorticity equation

$\begin{eqnarray}&&\displaystyle \frac{\partial q}{\partial t}+{\bf{v}}\cdot {\rm{\nabla }}q=({\rm{\nabla }}\times {{ \mathcal F }}_{{\bf{v}}})\cdot {{\bf{e}}}_{r}.\end{eqnarray} \tag{ 25 }$

The flow is entirely described by this equation for the absolute vorticity

$\begin{eqnarray}&&q=f+\zeta ,\end{eqnarray} \tag{ 26 }$

where the Coriolis parameter $f(\phi )=2{\rm{\Omega }}\;\mathrm{sin}\phi$ represents the vorticity of solid body rotation, and $\zeta =({\rm{\nabla }}\times {\bf{v}})\cdot {{\bf{e}}}_{r}$ is the relative vorticity in the radial direction ${{\bf{e}}}_{r}$ , relative to the rotating reference frame. The advection term ${\bf{v}}\cdot {\rm{\nabla }}q$ contains the advection of planetary vorticity ${\bf{v}}\cdot {\rm{\nabla }}f=\beta v$ , with $\beta ={a}^{-1}{\partial }_{\phi }f=2{\rm{\Omega }}{a}^{-1}\mathrm{cos}\phi$ , which arises from the curl of the Coriolis force. This term is commonly referred to as the β-term. The vorticity equation (25) contains the entire dynamics of the flow because an incompressible two-dimensional flow is described by a streamfunction ψ, defined such that

$\begin{eqnarray}&&{\bf{v}}={\rm{\nabla }}\times (\psi \;{{\bf{e}}}_{r}),\end{eqnarray} \tag{ 27a }$

$\begin{eqnarray}&&\zeta ={{\rm{\nabla }}}^{2}\psi .\end{eqnarray} \tag{ 27b }$

That is, if the relative vorticity ζ is known, the advecting velocity (27a) can be determined from the streamfunction, which is the solution of a Poisson equation (27b). Thus, the equations of motion (25) and (27), supplemented by appropriate boundary conditions, specify the dynamics completely.

The equations of motion (25)–(27) can be nondimensionalized using the planetary radius a as the typical length scale and the length of the day $2\pi {{\rm{\Omega }}}^{-1}$ as the typical time scale. With that nondimensionalization, the operators ${\rm{\nabla }}$ , ${\rm{\nabla }}\times$ and ${{\rm{\nabla }}}^{2}$ become operators on the unit sphere, and the angular velocity Ω becomes $2\pi$ . Throughout the rest of the paper, we will use there nondimensionalized quantities, unless otherwise specified.

Eddy-mean flow decomposition

We consider situations in which the boundary conditions of the problem are statistically symmetric under rotations around the planet's spin axis, so that the flow statistics (but not the instantaneous flow itself) can be expected to be axisymmetric. A zonal average $\bar{(\cdot )}$ then suggests itself. Decomposing flow fields in the nondimensional barotropic vorticity equation into zonal means $\bar{(\cdot )}$ and eddies ${(\cdot )}^{\prime }=(\cdot )-\bar{(\cdot )}$ yields the mean and eddy vorticity equations

$\begin{eqnarray}&&\displaystyle \frac{\partial \bar{\zeta }}{\partial t}=-\bar{{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })},\end{eqnarray} \tag{ 28a }$

$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial t}-\displaystyle \frac{\bar{u}}{\mathrm{cos}\phi }\displaystyle \frac{\partial }{\partial \lambda }\right){\zeta }^{\prime }=-[{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })-\bar{{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })}]-{v}^{\prime }\left(\beta +\displaystyle \frac{\partial \bar{\zeta }}{\partial \phi }\right).\end{eqnarray} \tag{ 28b }$

Here $\beta =2{\rm{\Omega }}\mathrm{cos}\phi$ , with ${\rm{\Omega }}=2\pi$ , is the nondimensional meridional derivative of the Coriolis parameter, and the Jacobian

$\begin{eqnarray}&&{J}_{n}(\psi ,\zeta )=\displaystyle \frac{1}{\mathrm{cos}\phi }\left(\displaystyle \frac{\partial \psi }{\partial \lambda }\displaystyle \frac{\partial \zeta }{\partial \phi }-\displaystyle \frac{\partial \psi }{\partial \phi }\displaystyle \frac{\partial \zeta }{\partial \lambda }\right)\end{eqnarray} \tag{ 29 }$

represents the advection on the unit sphere of vorticity ζ by the zonal (u) and meridional (v) flow implied by the streamfunction ψ:

$\begin{eqnarray}&&u=-\displaystyle \frac{\partial \psi }{\partial \phi },\qquad v=\displaystyle \frac{1}{\mathrm{cos}\phi }\displaystyle \frac{\partial \psi }{\partial \lambda }.\end{eqnarray} \tag{ 30 }$

The streamfunction-vorticity relations are

$\begin{eqnarray}&&\bar{\zeta }={{\rm{\nabla }}}_{n}^{2}\bar{\psi },\end{eqnarray} \tag{ 31a }$

$\begin{eqnarray}&&{\zeta }^{\prime }={{\rm{\nabla }}}_{n}^{2}{\psi }^{\prime },\end{eqnarray} \tag{ 31b }$

with ${{\rm{\nabla }}}_{n}$ likewise defined on the unit sphere.

Equation (28) is essentially (19) for the 2D barotropic vorticity equation. The mean flow evolves in time due to vorticity fluxes $-\bar{{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })}$ . The eddy vorticity budget involves shear by the mean flow $\bar{u}/\mathrm{cos}\phi \;{\partial }_{\lambda }{\zeta }^{\prime }$ , eddy–eddy interactions $[{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })-\bar{{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })}]$ , the β-term $\beta {v}^{\prime }$ and the advection of mean-flow vorticity ${\partial }_{\phi }\bar{\zeta }{v}^{\prime }$ .

Cumulants

The first cumulant is the mean vorticity $\bar{\zeta }(\phi ,t)$ , and the second cumulant is the vorticity equal-time two-point correlation:

$\begin{eqnarray}&&\bar{\zeta }({\bf{r}},t)=\bar{\zeta }(\phi ,t),\end{eqnarray} \tag{ 32a }$

$\begin{eqnarray}&&c({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\bar{{\zeta }^{\prime }({{\bf{r}}}_{1},t)\;{\zeta }^{\prime }({{\bf{r}}}_{2},t)}=C({\phi }_{1},{\phi }_{2},{\lambda }_{1}-{\lambda }_{2},t).\end{eqnarray} \tag{ 32b }$

The first cumulant depends on the latitude ϕ, and the second cumulant depends on the latitudes ${\phi }_{1}$ and ${\phi }_{2}$ and, because of the statistical axisymmetry of the equations, on the longitude difference ${\lambda }_{1}-{\lambda }_{2}$ (Marston et al 2008). Because the flow is entirely determined by the scalar q (or ζ), all other correlations are determined by the scalar cumulant c. Hence, moments like $\bar{{\zeta }^{\prime }({{\bf{r}}}_{1},t)\otimes {{\bf{u}}}^{\prime }({{\bf{r}}}_{2},t)}$ or $\bar{{{\bf{u}}}^{\prime }({{\bf{r}}}_{1},t)\otimes {{\bf{u}}}^{\prime }({{\bf{r}}}_{2},t)}$ can be calculated from $\bar{\zeta }$ and c (e.g., Marston et al 2008, Srinivasan and Young 2012, Marston et al 2014).

The CE2 equations for this problem are given in Marston et al (2008, 2014). They are of the form (18a) and (18c), with vorticity fluxes appearing as the essential eddy terms The equivalent QL system is (28) where the eddy–eddy interactions are neglected.

Numerical implementation

We simulate a barotropic flow on a sphere, specified by the equations of motion (25) and (27) on a spherical geodesic grid (Heikes and Randall 1995a, 1995b, Qi and Marston 2014) with $163,842$ cells. To remove enstrophy that cascades to the smallest scales, hyperviscous dissipation $\nu ({{\rm{\nabla }}}^{2}+2){{\rm{\nabla }}}^{6}\zeta$ is included, where the term $({{\rm{\nabla }}}^{2}+2)$ ensures that angular momentum is conserved. The hyperviscosity coefficient ν is chosen such that the smallest resolved mode decays with an e-folding time of 5. The vorticity is stepped forward in time by a second-order leapfrog scheme using the Robert–Asselin–Williams filter (Williams 2009). The time step is fixed at ${\rm{\Delta }}t=0.01$ .

Explicit time integration of the cumulant equations is carried out in spectral space using a 4th-order Runge–Kutta algorithm with an adaptive time step. We adopt the spectral truncation $0\leqslant {\ell }\leqslant L$ on spherical wavenumber ${\ell }$ , with the zonal wavenumbers restricted to $| m| \leqslant \mathrm{min}\{{\ell },M\}$ . We choose spectral cutoffs L = 60 and M = 30. Hyperviscosity is adjusted to ensure the same e-folding time at the maximum wavenumber ${\ell }=L$ as on the smallest resolved spatial scales on the geodesic grid.

To verify that the spectral cumulant simulation has sufficient resolution and can be compared to the geodesic grid model, a simulation of the fluid is also performed in a pure spectral calculation with the same numerical methods and resolution as for the cumulant equations. The agreement between the spectral and geodesic models is excellent. QL simulations are performed in spectral space by removing the triads that correspond to the interaction of two eddies, each with nonzero zonal wavenumber.

A program that implements fully nonlinear simulations on the spherical geodesic grid, and the nonlinear, QL, and CE2 simulations in spectral space, is freely available⁶ . More details about the simulations and CE2 can be found in Marston et al (2014).

4.2. Eddy lifecycle simulations

Setup

To illustrate situations when CE2 and QL approaches succeed or fail at capturing barotropic flow dynamics, we simulate the evolution of an initial zonal flow $U(\phi )$ with a superimposed initial disturbance (eddy) with vorticity ${\zeta }_{i}(\phi ,\lambda )$ . The zonal flow U and disturbance ${\zeta }_{i}$ mimic the upper-tropospheric jet stream and disturbances that may originate, for example, from lower-tropospheric baroclinic instability. The setup is inspired by Held and Phillips (1987) and uses

$\begin{eqnarray}&&U(\phi )=A\mathrm{cos}\phi -B{\mathrm{cos}}^{2}\phi +C{\mathrm{sin}}^{2}\phi {\mathrm{cos}}^{4}\phi -D{\mathrm{cos}}^{6}\phi ,\end{eqnarray} \tag{ 33a }$

$\begin{eqnarray}&&{\zeta }_{i}(\phi ,\lambda )={\zeta }_{0}\mathrm{cos}\phi \mathrm{exp}[-{(\phi -{\phi }_{m})}^{2}/{\delta }^{2})]\mathrm{cos}({k}_{i}\lambda ).\end{eqnarray} \tag{ 33b }$

To mimic Earth's upper troposphere, we choose

$\begin{eqnarray}&&A=3.4\times {10}^{-1},B=4.1\times {10}^{-1},C=4.0,\;{\rm{and}}\;{\rm{}}D=2.3\times {10}^{-1}.\end{eqnarray} \tag{ 34 }$

The corresponding dimensionalized flow on a sphere of Earth's radius and rotation rate resembles the midlatitude jet in the upper troposphere. It has a maximum westward velocity of $\sim 33\;{\rm{m}}\;{{\rm{s}}}^{-1}$ at a latitude of $\sim 40^\circ$ and a maximum eastward velocity of $\sim 22\;{\rm{m}}\;{{\rm{s}}}^{-1}$ at the equator; its zero is near $\sim 17^\circ$ . This zonal flow is barotropically stable. The eddy vorticity ${\zeta }_{i}$ decays meridionally away from its maximum absolute value ${\zeta }_{0}\mathrm{cos}{\phi }_{m}$ at latitude ${\phi }_{m}=\pi /4=45^\circ$ , with characteristic meridional decay scale $\delta =\pi /9$ . Its zonal wavenumber ${k}_{i}=6$ is close to the dominant zonal wavenumber of transient eddies on Earth (e.g., Boer and Shepherd 1983, Straus and Ditlevsen 1999), which approximately coincides with the baroclinically most unstable zonal wavenumber (e.g., Simmons and Hoskins 1976, 1978, Thorncroft et al 1993, Merlis and Schneider 2009).

We let the flow evolve freely without forcing or dissipation, apart from hyperviscosity, and analyze the time evolution of the mean flow and the eddies. We compare CE2 to the statistics of fully nonlinear simulations for different choices of parameters. To identify relevant nondimensional parameters controlling the evolution of the flow and the adequacy of CE2 closures, we rescale the mean and eddy vorticity equations (28), using the relative vorticity $2\mathrm{Ro}$ of the initial zonal flow U. Dimensionally, we have $\mathrm{Ro}={\rm{\Lambda }}/(2{\rm{\Omega }})$ , where ${\rm{\Lambda }}\approx 2\;\mathrm{max}(U)/a$ is the typical initial zonal-mean flow vorticity. Hence, $\mathrm{Ro}$ is a Rossby number measuring the mean flow vorticity to the equatorial planetary vorticity $2{\rm{\Omega }}$ . We measure the initial maximum vorticity ${\zeta }_{0}$ of the eddies relative to the mean-flow vorticity $2\mathrm{Ro}$ through the amplitude parameter $\epsilon ={\zeta }_{0}\mathrm{cos}{\phi }_{m}/2\mathrm{Ro}$ . The quantities $\bar{\zeta }$ , ${\zeta }^{\prime }$ , $\bar{\psi }$ , ${\psi }^{\prime }$ , and t are then rescaled with $4\pi \mathrm{Ro}$ , $4\pi \epsilon \mathrm{Ro}$ , $4\pi \mathrm{Ro}$ , $4\pi \epsilon \mathrm{Ro}$ , and ${(4\pi \mathrm{Ro})}^{-1}$ , respectively. The equations of motion for the mean-flow and the eddies become

$\begin{eqnarray}&&\displaystyle \frac{\partial \bar{\zeta }}{\partial t}=-{\epsilon }^{2}\bar{{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })},\end{eqnarray} \tag{ 35a }$

$\begin{eqnarray}\mathrm{Ro}\left(\displaystyle \frac{\partial }{\partial t}-\displaystyle \frac{1}{\mathrm{cos}\phi }\displaystyle \frac{\partial \bar{\psi }}{\partial \phi }\displaystyle \frac{\partial }{\partial \lambda }\right){\zeta }^{\prime } & = & -\epsilon \mathrm{Ro}[{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })-\bar{{J}_{n}({\psi }^{\prime },{\zeta }^{\prime })}]-\displaystyle \frac{\partial {\psi }^{\prime }}{\partial \lambda }\left(1+\mathrm{Ro}\displaystyle \frac{1}{\mathrm{cos}\phi }\displaystyle \frac{\partial \bar{\zeta }}{\partial \phi }\right),\end{eqnarray} \tag{ 35b }$

equation (35) gives the relative amplitude of the different terms if we assume that the typical length scales of the mean flow and eddies are of the order of the radius of the planet. An immediate simplification that results for small Rossby number $\mathrm{Ro}$ (the case we will consider) is that the advection of mean-flow vorticity $\bar{\zeta }$ by meridional velocity fluctuations is negligible compared with the β-term, which is a factor ${\mathrm{Ro}}^{-1}$ larger.

The nondimensional parameters $\mathrm{Ro}$ and control the vorticity of the mean flow and of the eddies and are important for the evolution of the barotropic flow. Eddy–mean flow interactions are of order $\mathrm{Ro}$ , and eddy–eddy interactions of order $\epsilon \mathrm{Ro}$ , provided ${{\rm{\nabla }}}_{n}^{2}$ is of order one. This is the case initially; however, it does not remain true over the evolution of the flow, as small scales are generated. The two parameters $\mathrm{Ro}$ and can be varied independently in our setup. In what follows, we explore how these parameters affect the flow evolution and the adequacy of CE2 and QL approaches in capturing it.

4.3. Results

Varying eddy amplitudes

For a fixed initial mean-flow Rossby number $\mathrm{Ro}\approx 0.06$ (corresponding to the Earth-like parameters in equation (34)), we run eddy lifecycle experiments for larger-amplitude initial eddies with $\epsilon \approx 6$ , and for smaller-amplitude initial eddies with $\epsilon \approx 2$ . The expectation is that CE2 and QL approaches are more successful for the smaller-amplitude eddies, for which the nonlinear eddy–eddy interactions (of order $\epsilon \mathrm{Ro}$ ) are weaker, and this is indeed borne out in the simulations. It is instructive to see in what way they fail to capture aspects of the flow evolution for the larger-amplitude eddies.

For the larger-amplitude eddies, figure 4 shows the relative vorticity ζ at 5 instants during the evolution of the flow. It is evident that the initial disturbance quickly becomes nonlinear and develops drawn-out filaments on the equatorward flank of the zonal jet. The filaments roll-up anticyclonically within cats' eyes structures (marked by X's in figure 4) that continue to have the initial zonal wavenumber ${k}_{i}=6$ . Such cats' eyes are characteristic of Rossby waves that break in 'surf zones' around their critical layers⁷ (Stewartson 1977, Warn and Warn 1978, McIntyre and Palmer 1983, Killworth and McIntyre 1985, Held and Phillips 1987).

**Figure 4.** Evolution of relative vorticity in the fully nonlinear simulation for the larger-amplitude eddies ( $\epsilon =6$ ) for an Earth-like setting ( $\mathrm{Ro}=0.06$ ). The relative vorticity maps show the formation of cats' eyes, with vorticity filaments rolling up within them. White X's mark the centers of some cats' eyes. Time scales are nondimensionalized with the length of the day $2\pi {{\rm{\Omega }}}^{-1}$ .
Download figure:
Standard image High-resolution image

As the eddies break and eventually dissipate in the surf zone, they are absorbed by the mean flow, and their kinetic energy decays. The total eddy kinetic ${E}_{{\rm{K}}}=0.5{\displaystyle \int }_{\phi }{e}_{{\rm{K}}}\mathrm{cos}\phi \;{\rm{d}}\phi$ , where ${e}_{{\rm{K}}}=0.5(\bar{{{u}^{\prime }}^{2}}+\bar{{{v}^{\prime }}^{2}})$ , becomes very small at large times ( $\gtrsim 30$ , see figure 5(a)). At those times, most of the initial kinetic energy has been transferred to the mean zonal flow. The local zonal kinetic energy ${e}_{Z}=0.5{\bar{u}}^{2}$ ( $\bar{u}\equiv -{\partial }_{\phi }\bar{\psi }$ of the mean zonal flow increases in the core of the midlatitude jet, roughly in proportion to the decrease of the eddy kinetic energy e_K (figure 5(c)). Total energy ${E}_{{\rm{K}}}+{E}_{Z}$ , with ${E}_{Z}=0.5{\displaystyle \int }_{\phi }{e}_{Z}\mathrm{cos}\phi \;{\rm{d}}\phi$ , is approximately conserved in the model, up to very small losses ( $\sim 0.2\%$ of the total after a time of 50) through SGS dissipation. That is, although dissipation at small scales in the surf zone is essential to generate irreversibility, the kinetic energy loss associated with it is small compared to the transfer to the mean flow.

**Figure 5.** Evolution of kinetic energy and eddy momentum flux convergence (EMFC) in the fully nonlinear simulation (left column) and in a direct CE2 calculation of the statistics (right column) for an Earth-like setting ( $\mathrm{Ro}=0.06$ ) and for the larger-amplitude eddies ( $\epsilon =6$ ). (a), (d) Eddy kinetic energy e_K. (b), (e) Zonal kinetic energy e_Z. (c), (f) EMFC. Vertical blue dashed lines indicate at which times relative vorticity snapshots are shown in figure 4. Time scales and length scales are nondimensionalized with the length of the day $2\pi {{\rm{\Omega }}}^{-1}$ and the planet radius a.
Download figure:
Standard image High-resolution image

**Figure 5.** Evolution of kinetic energy and eddy momentum flux convergence (EMFC) in the fully nonlinear simulation (left column) and in a direct CE2 calculation of the statistics (right column) for an Earth-like setting ( $\mathrm{Ro}=0.06$ ) and for the larger-amplitude eddies ( $\epsilon =6$ ). (a), (d) Eddy kinetic energy e_K. (b), (e) Zonal kinetic energy e_Z. (c), (f) EMFC. Vertical blue dashed lines indicate at which times relative vorticity snapshots are shown in figure 4. Time scales and length scales are nondimensionalized with the length of the day $2\pi {{\rm{\Omega }}}^{-1}$ and the planet radius a.
Download figure:
Standard image High-resolution image

The transfer of E_K to E_Z implies acceleration of the mean zonal jet. This acceleration occurs through transfer of momentum from the eddies to the mean flow, as can be seen from the nondimensionalized zonally averaged momentum equation (24a) in the inviscid limit ( ${{ \mathcal F }}_{{\bf{v}}}=0$ ):

$\begin{eqnarray}&&\mathrm{cos}\phi \displaystyle \frac{\partial \bar{u}}{\partial t}=-\displaystyle \frac{1}{\mathrm{cos}\phi }\displaystyle \frac{\partial }{\partial \phi }[{\mathrm{cos}}^{2}\phi \;\bar{{u}^{\prime }{v}^{\prime }}].\end{eqnarray} \tag{ 36 }$

Acceleration of the mean zonal flow occurs where eddy momentum fluxes $\bar{{u}^{\prime }{v}^{\prime }}\mathrm{cos}\phi$ converge. Multiplying the mean zonal momentum equation (36) by $\bar{u}$ and integrating by parts yields the equation for the zonal kinetic energy E_Z,

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{E}_{Z}={\displaystyle \int }_{\phi }{\mathrm{cos}}^{2}\phi \;\bar{{u}^{\prime }{v}^{\prime }}\displaystyle \frac{\partial }{\partial \phi }\left(\displaystyle \frac{\bar{u}}{\mathrm{cos}\phi }\right)\;{\rm{d}}\phi =-\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{E}_{{\rm{K}}},\end{eqnarray} \tag{ 37 }$

where the right-hand side is obtained from a corresponding integral of the eddy momentum equations. This shows that transfer of kinetic energy from the eddies to the mean flow occurs through eddy momentum fluxes that are up the gradient of angular velocity $\bar{u}{\mathrm{cos}}^{-1}\phi$ . The acceleration of the mean zonal jet at its core (figure 5(c)) thus is associated with eddy momentum flux convergence ${\partial }_{\phi }(\bar{{u}^{\prime }{v}^{\prime }}\mathrm{cos}\phi )$ (EMFC, see figure 5(e)).

However, the eddy kinetic energy does not decay monotonically. Instead, it exhibits damped oscillations during which eddy momentum fluxes cause zonal angular momentum to slosh back and forth meridionally within the jet (figure 5(e)). The alternating poleward and equatorward momentum fluxes (with decreasing amplitude) on the equatorward flank of the jet are result of nonlinear processes within the surf zone. These processes have been described in an idealized analytical model of Rossby-wave breaking in critical layers, the Stewartson–Warn–Warn (SWW) theory (Stewartson 1977, Warn and Warn 1978; see also Killworth and McIntyre 1985). The oscillation on the poleward flank of the jet may be more linear, originating from the reflection of Rossby wavepackets from reflecting levels that arise because β decreases with latitude (e.g., Lorenz 2014).

Figure 5 (right column) shows the kinetic energies and EMFC obtained from a direct calculation of these statistics with CE2. CE2 captures the oscillation of kinetic energy between eddies and the mean zonal flow, with a period similar to the fully nonlinear simulation (figures 5(a), (b)). However, the oscillations are too weakly damped; large eddy kinetic energies e_K persist for a long time. The eddy absorption in the surf zone is not adequately captured by CE2 because CE2 cannot resolve the generation of small scales in the surf zone through eddy–eddy interactions. Consistently, unrealistically strong oscillations persist in the EMFC under CE2 (figure 5(f)). How these oscillations arise from the perspective of wave mechanics, and why CE2 cannot capture the wave absorption in this case, is illustrated in appendix B in a QL simulation that corresponds to the CE2 calculations shown here. The phenomenology of such oscillations has been described analytically by Haynes and McIntyre (1987) in the context of a QL truncation of the SWW model.

Eddy–eddy interactions are weaker and CE2 is more successful in capturing the flow dynamics when the amplitude of the initial perturbation is decreased by a factor 3 ( $\epsilon \approx 2$ ). Cats' eyes with rolling-up vorticity filaments do not develop in the fully nonlinear simulation (figure 6). Instead, eddies are sheared by the mean flow, which transfers eddy kinetic energy e_K to the mean flow through the Orr mechanism, which is only weakly nonlinear because it involves the interaction of disturbances with the slowly varying mean flow (e.g. Thomson 1887, Orr 1907, Farrell 1987, Bouchet et al 2013). The transfer of eddy kinetic energy to the mean flow occurs over time scales of a couple days, corresponding to the shear time scale of the mean zonal flow (figure 7). The damped oscillatory behavior seen in the larger-amplitude simulation disappears. Because eddy absorption results from the mean flow shearing the eddies—an eddy-mean flow interaction that is captured by CE2—statistics calculated directly with CE2 come in very close agreement with those from the fully nonlinear simulation (figure 7, right column). As in the nonlinear case, eddy absorption occurs through the formation of small-scale vorticity filaments. But instead of rolling up inside cat's eyes, here they stretch around the planet.

**Figure 6.** Evolution of relative vorticity in the fully nonlinear simulation for the smaller-amplitude eddies ( $\epsilon =2$ ) for an Earth-like setting ( $\mathrm{Ro}=0.06$ ). Time scales are nondimensionalized with the length of the day $2\pi {{\rm{\Omega }}}^{-1}$ .
Download figure:
Standard image High-resolution image

**Figure 6.** Evolution of relative vorticity in the fully nonlinear simulation for the smaller-amplitude eddies ( $\epsilon =2$ ) for an Earth-like setting ( $\mathrm{Ro}=0.06$ ). Time scales are nondimensionalized with the length of the day $2\pi {{\rm{\Omega }}}^{-1}$ .
Download figure:
Standard image High-resolution image

**Figure 7.** Evolution of kinetic energy and eddy momentum flux convergence (EMFC) in the fully nonlinear simulation (left column) and in a direct CE2 calculation of the statistics (right column) for an Earth-like setting ( $\mathrm{Ro}=0.06$ ) and for the larger-amplitude eddies ( $\epsilon =2$ ). (a), (d) Eddy kinetic energy e_K. (c), (f) EMFC. Vertical blue dashed lines indicate the times corresponding to relative vorticity snapshots on figure 4. Time scales are nondimensionalized with the length of the day $2\pi {{\rm{\Omega }}}^{-1}$ and length scales with the planet radius a.
Download figure:
Standard image High-resolution image

It is worth noting that eddies are also sheared equatorward of the cats' eyes in the larger-amplitude simulation (figure 4). Hence, weakly nonlinear eddy absorption also occurs in this simulation, but it is not the dominant effect responsible for eddy absorption.

Relative amplitude of terms in the potential vorticity budget

For the larger-amplitude experiment, the evident importance of the development of small scales through eddy–eddy interactions seems consistent with the order of magnitude of the terms in the vorticity equations (35). The eddy–eddy interactions appear of order $\epsilon \mathrm{Ro}\approx 0.4$ , compared with the β-term which is responsible for Rossby wave dynamics and is of order 1. Hence, the eddy–eddy interactions are not negligible compared with Rossby wave dynamics. By contrast, the interactions of the disturbance with the mean flow shear are of order $\mathrm{Ro}\approx 0.06$ and hence are much weaker.

For the smaller-amplitude experiment, the dimensional analysis of the vorticity equations (35) suggests that the eddy–eddy interactions now are of order $\epsilon \mathrm{Ro}\lesssim 0.2$ compared with the β-term. Hence, the eddy–eddy interactions become close to being negligible compared with Rossby wave dynamics, consistent with the simulation results. However, the dimensional analysis also suggests that interactions of the disturbance with the mean flow shear still are of order $\mathrm{Ro}\approx 0.06$ and hence are weaker still, albeit of the same order of magnitude as the eddy–eddy interactions. Yet the eddy–eddy interactions are inhibited in the fully nonlinear simulation, whereas the shear interactions dominate the eddy absorption, illustrating the limits of dimensional analysis in this nonlinear problem.

Equation (35) is useful to determine which parameter to vary. Nevertheless, one has to be careful when interpreting the relative amplitude of the different terms A small term can be fundamental for the dynamics. For example, the SWW theory has shown that eddy–eddy interactions, albeit weak, can determine at leading order momentum fluxes because they dominate the vorticity budget in a thin critical layer, where linear theory breaks down. This remains true in the limit where the relative amplitude of the eddy–eddy interactions goes to zero. Moreover, the relative amplitude of the terms in the vorticity budget evolves with time as small scales develop.

Varying Rossby numbers

To illustrate how variations of the mean-flow Rossby number affect the evolution of disturbances, we use larger-amplitude eddies ( $\epsilon \approx 6$ ) and decrease the Rossby number $\mathrm{Ro}$ from 0.06 to 0.02. Dimensionally, this is equivalent to weakening the initial flow while keeping the rotation rate of the planet constant, or to increasing the rotation rate while keeping the initial flow constant. Based on the dimensional analysis of the vorticity equations (35), this reduction of the Rossby number should decrease the relative magnitude both of eddy–eddy interactions and of shearing of eddies by the mean flow relative to the β-term, maintaining the relative amplitude of the eddy–eddy interactions to the β-term.

Figure 8 compares the time evolution of the eddy kinetic energy e_K in the fully nonlinear and the CE2 simulations. As we have already seen, eddy absorption at $\mathrm{Ro}=0.06$ is not captured by CE2 because it is strongly nonlinear. However, at the smaller Rossby number $\mathrm{Ro}=0.02$ , it is faithfully captured by CE2. Indeed, eddy absorption appears more weakly nonlinear for $\mathrm{Ro}=0.02$ , and dominated by mean-flow shearing, as is evident on the relative vorticity maps (figure 9). Vorticity maps for $\mathrm{Ro}=0.04$ and $\mathrm{Ro}=0.03$ show the transition from weakly to strongly nonlinear absorption: the meridional extent of the nonlinear surf zone decreases with increasing rotation, while shearing effects become more important (figure 9). Because CE2 can capture the weakly nonlinear shearing interactions but not the strongly nonlinear eddy–eddy interactions in the surf zone, it becomes gradually more adequate as the rotation rate increases (figure 8). This occurs although the orders of magnitude of the relevant terms suggested by the dimensional analysis decrease by equal factors of order $\mathrm{Ro}$ as the mean-flow Rossby number decreases. It appears that what is important here is that the magnitude of the advection of absolute vorticity, and hence of linear Rossby wave dynamics, increases relative to both of these terms. For $\mathrm{Ro}=0.02$ , shear explains eddy decay and jet acceleration, even though the nondimensionalization (35) suggests shear should be much smaller than eddy–eddy interactions.

**Figure 9.** Evolution of relative vorticity in the fully nonlinear simulation for the larger-amplitude eddies ( $\epsilon =6$ ), with Rossby number decreasing from $\mathrm{Ro}=0.06$ to 0.02. For comparison, the first row reproduces the vorticity maps already shown in figure 4. Time scales are nondimensionalized with the length of the day $2\pi {{\rm{\Omega }}}^{-1}$ .
Download figure:
Standard image High-resolution image

Higher-order closures

CE2 fails to capture eddy absorption when eddy–eddy interactions are important for the dynamics. We tested whether a higher-order closure (CE3*) captures eddy absorption more faithfully. CE3* is described in Marston et al (2014). It truncates the cumulant equations at third order, ensuring realizability by projecting out modes with (unphysical) negative energies.

The results are summarized in figure 10. Because CE3* is computationally expensive, the resolution has been reduced to M = 40 and L = 20. For simplicity we turn off eddy damping ( $\tau =\infty$ , see Marston et al 2014 for definitions). The full and CE2 simulations in figure 10 are run at this lower resolution and are consistent with the higher-resolution runs (figure 5); however, they do exhibit a faster eddy damping because of the stronger diffusion. It appears that CE3* captures the eddy absorption very accurately. We also tested other closures that take into account third-order terms (Marston et al 2014); they bring similar improvements.

4.4. Implications

The results show that direct CE2 calculation of barotropic flow statistics representative of the upper troposphere can succeed in circumstances when the dominant nonlinear interaction is between eddies and the mean flow, for example, by shearing. They fail when strongly nonlinear eddy–eddy interactions become important in surf zones around critical layers, where the roll-up of vorticity filaments leads to the generation of small scales. This is a process that cannot be captured in CE2. However, higher-order closures, which take some effects of third cumulants on second cumulants into account, can perform better in such cases—at the price of increased conceptual and computational complexity.

Weakly nonlinear eddy-mean flow interactions seem to be favored over strongly nonlinear eddy–eddy interactions when the eddy vorticity is small enough compared with the planetary vorticity. The transition from weakly to strongly nonlinear interactions predominating occurs above a critical value of eddy amplitude that is a decreasing function of the Rossby number $\mathrm{Ro}$ . The parameter need not be small for absorption through weakly nonlinear shearing to occur. For example, for small $\mathrm{Ro}$ , weakly nonlinear eddy absorption through shearing seems to be favored even when a large value of suggests that nonlinear eddy–eddy interactions should be larger than the mean-flow shearing of eddies.

When strongly nonlinear eddy–eddy interactions are favored, high eddy kinetic energies develop in the QL/CE2 approximation because momentum sloshes back and forth meridionally within the jet, without sufficient absorption. Lifecycle experiments carried out with a QL GCM have shown that this phenomenology is relevant to the upper troposphere in an Earth-like setting, highlighting the relevance of the simplified barotropic model: earth's upper troposphere appears to be in a regime in which nonlinear eddy–eddy interactions in surf zones are important for the structure of the momentum fluxes (Ait-Chaalal and Schneider 2015).

5. Conclusions

Atmospheric flows are highly anisotropic and inhomogeneous, with rich spatial structures. Turbulent closures that respect the anisotropy and inhomogeneity may enable the direct statistical simulation of Earth's atmosphere (Marston 2012). Expansion of statistics in equal-time cumulants yields equations of motion for the statistics that can already provide useful closures at second order (CE2), because the mean flows and interactions of perturbations with them are strong (Herring 1963, O'Gorman and Schneider 2007, Marston et al 2008, Tobias and Marston 2013, Marston et al 2014). CE2 solves the first two cumulant (central moment) equations of the QL approximation, in which interactions between mean flows and fluctuations around it are retained, while nonlinear eddy–eddy interactions are neglected. In section 2, we formulated CE2 for the Boussinesq equations by introducing a condensed tensorial notation. The case of the anelastic equations is presented in appendix A and involves the use of density weighted averages. We tested the relevance of CE2 to two distinct atmospheric flows involving different length and time scales and force balances: turbulent convection in the atmospheric boundary layer (section 3), and weak two-dimensional turbulence representative of the upper troposphere (section 4).

Convection in the atmospheric boundary layer links large-scale atmospheric dynamics aloft to the surface underneath, mediating the exchanges of momentum, energy, and water between the surface and the free troposphere. Motions in boundary layers and in clouds that have their roots in them have dynamical scales of meters, meaning that they need to be parameterized in GCMs. Current parameterization schemes have numerous shortcomings; our inability to represent cloud and boundary layer dynamics adequately in climate models is the largest source of uncertainty in climate change projections. Because CE2s capture interactions between fluctuations (e.g., thermals) and mean fields and take nonlocal correlations of fluctuations into account, without requiring the introduction of tunable parameters that proliferate in current parameterization schemes, they may offer a way to achieve more physically consistent parameterizations. We presented encouraging initial results, showing that a QL large-eddy simulation of a dry convective boundary layer captures the first-order statistics (e.g., mean boundary layer growth) of a corresponding fully nonlinear LES. However, it does not capture second-order statistics adequately. More work is required to investigate to what degree these results hold generally, in broader classes of boundary layer flows and in the presence of moisture effects and clouds, and how the QL results can be improved by including representations of higher-order effects, such as the turbulent transport of kinetic energy.

The potential for development of parameterization schemes based on CE2 is a promising direction for future research. But it requires overcoming both technical and theoretical challenges. On the technical level, fast numerical methods are required for the method to be competitive with other subgrid schemes. This could be achieved by dimensional reduction to capture only the most important nonlocal correlations in the second cumulant. At the theoretical level, it is not clear to what extent CE2 and possible extensions will be able to describe moist convection, or what effects vertical shear will have on its accuracy. It appears necessary to include some effects of eddy–eddy interactions, such as those captured by the third order CE3* approximation (Marston et al 2014).

At the planetary scale, how and where eddies in the upper troposphere dissipate controls the strength and direction of momentum fluxes and thus climatic features such as surface winds. A one-layer barotropic model that mimics the behavior of the upper troposphere illustrates different mechanisms through which eddy absorption by the mean flow can occur. CE2 describes eddy absorption well when it occurs through shearing of eddies by the mean flow. This happens when the vorticity that characterizes the eddies is small compared with the planetary vorticity (planetary rotation rate). When the eddy vorticity is large, CE2 is not adequate because eddy absorption results from the formation of small scales that form through eddy–eddy interactions in critical layers. A comprehensive theory that describes these weakly and strongly nonlinear absorptive regimes is lacking.

Our results suggest that, in general, higher-order closures are required for an accurate direct statistical simulation of large-scale and smaller-scale atmospheric flows. We have tested a few of them in the large-scale context and found improvement compared with CE2. Nevertheless, going beyond CE2 raises a number of questions. Higher-order closures are several orders of magnitude slower than CE2; currently, they are much slower than direct simulation of the flows, while CE2 can be faster that direct simulations. Dimensional reduction may be able to help here. More generally, once eddy–eddy interactions are taken into account, the whole hierarchy of cumulants is active and not completely described by any finite closure. Realizability of closures becomes an issue (e.g., Orszag 1970, Orszag et al 1973, Marston et al 2014), and it is known that intermittency cannot be adequately described in this way (Frisch 1995, Lesieur 2008). But the existence of anisotropic and inhomogeneous mean flows provides a starting point for a systematic exploration of statistical closures.

Acknowledgments

We thank Kyle Pressel for helping in the development of the QL LES and Joe Skitka for sharing results about the oceanic boundary layer. We are grateful for helpful discussions with Freddy Bouchet, Greg Chini, Baylor Fox-Kemper, Cesare Nardini, and Steve Tobias. The comments of two anonymous reviewers improved the manuscript significantly. This work was supported by the US National Science Foundation under grants CCF-1048575 (FAC and TS) and CCF-1048701 (FAC and JBM).

Appendix A.: CE2 of anelastic flow

Because the flow is nondivergent in the Boussinesq equations, it is possible to directly use any averaging operator satisfying the properties (7). However, for more general flows, in which the density may vary, the requirements on the averaging operator need to be modified so that second-order equations for fluctuations consistent with the conservation laws are obtained. Because density appears as a weight in all integrals of conserved quantities over the flow domain, the density weighted average ${\bar{(\cdot )}}^{*}=\bar{(\rho \cdot )}/\bar{\rho }$ , with $\bar{(\cdot )}$ satisfying the Reynolds properties (7), leads to a decomposition of the flow that is amenable to closures such as CE2. The density weighted average satisfies (7a)–(7c) but, importantly, not commutativity with partial differentiation (7d). Thus, we can define an eddy mean flow decomposition

$\begin{eqnarray}&&f={\bar{f}}^{*}+\hat{f},\end{eqnarray} \tag{ A.1a }$

$\begin{eqnarray}&&{\bar{{fg}}}^{*}={\bar{f}}^{*}{\bar{g}}^{*}+{\bar{\hat{f}\hat{g}}}^{*},\end{eqnarray} \tag{ A.1b }$

with hats $\hat{(\cdot )}=(\cdot )-{\bar{(\cdot )}}^{*}$ denoting fluctuations around the density weighted mean. This is sometimes called the Favre decomposition.

How a varying density affects cumulants and CE closures can be illustrated with the anelastic equations, an extension of the Boussinesq equations that allows the reference density ${\rho }_{0}={\rho }_{0}(z)$ to vary with height z. In the anelastic approximation, density perturbations $\delta \rho$ , the pressure potential ${\rm{\Phi }}=\delta p/{\rho }_{0}$ , and the buoyancy $b=-g\delta \rho /{\rho }_{0}$ are defined relative to this height-dependent and hydrostatically balanced reference density. The state vector ${\boldsymbol{\Psi }}$ is defined as in (2). We have to consider covariances involving ${\rm{\nabla }}{\rm{\Phi }}$ and not Φ because multiplying by density does not commute with derivatives. For simplicity, we assume that the linear operator ${ \mathcal L }$ is homogeneous in the reference density and satisfies ${\rho }_{0}{ \mathcal L }(\tilde{{\boldsymbol{\Psi }}})={ \mathcal L }({\rho }_{0}\tilde{{\boldsymbol{\Psi }}})$ . It implies that it cannot involve any spatial derivative. The anelastic equations in a reference frame rotating with angular velocity ${\boldsymbol{\Omega }}$ are (Vallis 2006):

$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{0}{\boldsymbol{\Psi }}}{\partial t}+{\rm{\nabla }}\cdot [{\rho }_{0}{\boldsymbol{\Psi }}\otimes {\bf{u}}]={\bf{L}}[{\rho }_{0}{\boldsymbol{\Psi }}]-{\rho }_{0}{\rm{\nabla }}{\rm{\Phi }}+{\bf{F}},\end{eqnarray} \tag{ A.2a }$

$\begin{eqnarray}&&{\rm{\nabla }}\cdot {{\boldsymbol{\rho }}}_{{\bf{0}}}{\bf{u}}=0.\end{eqnarray} \tag{ A.2b }$

The Boussinesq equations (1) are obtained from the anelastic equations by setting ${\rho }_{0}$ to a constant. The continuity equation (A.2b) again reduces to a nondivergence constraint. Because the reference density ${\rho }_{0}$ appears in it, we need to use the reference-density weighted average ${\bar{(\cdot )}}^{*}=\bar{({\rho }_{0}\cdot )}/{\bar{\rho }}_{0}$ to derive the CE2.

The equation for the first cumulant ${\bar{{\boldsymbol{\Psi }}({\bf{r}},t)}}^{*}$ is obtained by averaging (A.2):

$\begin{eqnarray}&&\displaystyle \frac{\partial {\bar{\rho }}_{0}{\bar{{\boldsymbol{\Psi }}}}^{*}}{\partial t}+{\rm{\nabla }}\cdot [{\bar{\rho }}_{0}{\bar{{\boldsymbol{\Psi }}}}^{*}\otimes {\bar{{\bf{u}}}}^{*}]=-{\rm{\nabla }}\cdot [{\bar{\rho }}_{0}{\bar{\hat{{\boldsymbol{\Psi }}}\otimes \hat{{\bf{u}}}}}^{*}]+{\bf{L}}[{\bar{\rho }}_{0}{\rm{}}\;{\rm{}}{\bar{{\boldsymbol{\Psi }}}}^{*}]-{\bar{\rho }}_{0}{\overline{{\rm{\nabla }}\phi }}^{*},\end{eqnarray} \tag{ A.3a }$

$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\bar{\rho }}_{0}{\rm{}}\;{\rm{}}{\bar{{\bf{u}}}}^{*}=0.\end{eqnarray} \tag{ A.3b }$

To define a second central moment or second cumulant, we need a quantity that respects the symmetry (11) of the covariance matrix under exchange of the spatial coordinates ${{\bf{r}}}_{1}$ and ${{\bf{r}}}_{2}$ , that gives the second-order term in the first cumulant equation (A.3), and that corresponds to a density weighted average when ${{\bf{r}}}_{1}={{\bf{r}}}_{2}$ . A possible choice is

$\begin{eqnarray}&&{{\bf{C}}}_{+}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\displaystyle \frac{1}{2}\bar{[{\rho }_{0}({{\bf{r}}}_{1})+{\rho }_{0}({{\bf{r}}}_{2})]\hat{{\boldsymbol{\Psi }}}({{\bf{r}}}_{1},t)\otimes \hat{{\boldsymbol{\Psi }}}({{\bf{r}}}_{2},t)}.\end{eqnarray} \tag{ A.4 }$

We also define the two auxiliary covariances corresponding to (12)

$\begin{eqnarray}&&{{\bf{C}}}_{+}^{{\rm{\Phi }}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\displaystyle \frac{1}{2}\bar{[{\rho }_{0}({{\bf{r}}}_{1})+{\rho }_{0}({{\bf{r}}}_{2})]\hat{{\rm{\Phi }}}({{\bf{r}}}_{1},t)\hat{{\boldsymbol{\Psi }}}({{\bf{r}}}_{2},t)},\end{eqnarray} \tag{ A.5a }$

$\begin{eqnarray}&&{{\bf{C}}}_{+}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\displaystyle \frac{1}{2}\bar{[{\rho }_{0}({{\bf{r}}}_{1})+{\rho }_{0}({{\bf{r}}}_{2})]\hat{{\boldsymbol{\Psi }}}({{\bf{r}}}_{2},t)\otimes \hat{{\bf{u}}}({{\bf{r}}}_{1},t)},\end{eqnarray} \tag{ A.5b }$

and introduce the corresponding quantities, ${{\bf{C}}}_{-}$ , ${{\bf{C}}}_{+}^{{\rm{\Phi }}}$ , and ${{\bf{C}}}_{-}^{u}$ with the density difference instead of the sum, as in

$\begin{eqnarray}&&{{\bf{C}}}_{-}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)=\displaystyle \frac{1}{2}\bar{[{\rho }_{0}({{\bf{r}}}_{1})-{\rho }_{0}({{\bf{r}}}_{2})]\hat{{\boldsymbol{\Psi }}}({{\bf{r}}}_{1},t)\otimes \hat{{\boldsymbol{\Psi }}}({{\bf{r}}}_{2},t)}.\end{eqnarray} \tag{ A.6 }$

Some algebra then gives the equations of motion for the second cumulant:

$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{{\bf{C}}}_{\pm }({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+{{\rm{\nabla }}}_{{{\bf{r}}}_{1}}\cdot [{{\bf{C}}}_{\pm }({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\otimes {\bar{{\bf{u}}}}^{*}({{\bf{r}}}_{1},t)]+{{\rm{\nabla }}}_{{{\bf{r}}}_{2}}\cdot [{{\bf{C}}}_{\pm }^{{\rm{T}}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\otimes {\bar{{\bf{u}}}}^{*}({{\bf{r}}}_{2},t)]\\ &&=-{{\bf{C}}}_{\pm }^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t){[{\rm{\nabla }}{\bar{{\rm{\Psi }}}}^{*}({{\bf{r}}}_{2},t)]}^{{\rm{T}}}-{\rm{\nabla }}{\bar{{\rm{\Psi }}}}^{*}({{\bf{r}}}_{2},t){{\bf{C}}}_{\pm }^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\\ &&\quad -\displaystyle \frac{1}{2}{{\bf{C}}}_{\pm }({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)[{\rm{\nabla }}\cdot {\bar{{\bf{u}}}}^{*}({{\bf{r}}}_{1},{\bf{t}})+{\rm{\nabla }}\cdot {\bar{{\bf{u}}}}^{*}({{\bf{r}}}_{2},{\bf{t}})]\\ &&\quad -\displaystyle \frac{1}{2}{{\bf{C}}}_{\mp }({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)[{\rm{\nabla }}\cdot {\bar{{\bf{u}}}}^{*}({{\bf{r}}}_{1},{\bf{t}})-{\rm{\nabla }}\cdot {\bar{{\bf{u}}}}^{*}({{\bf{r}}}_{2},{\bf{t}})]\\ &&\quad +{\bf{L}}\;{{\bf{C}}}_{\pm }({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+{{\bf{C}}}_{\pm }({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)\;{{\bf{L}}}^{{\rm{T}}}+{{\rm{\nabla }}}_{{{\bf{r}}}_{1}}{{\bf{C}}}_{\pm }^{{\rm{\Phi }}}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+{{\rm{\nabla }}}_{{{\bf{r}}}_{2}}{{\bf{C}}}_{\pm }^{{\rm{\Phi }}}({{\bf{r}}}_{2},{{\bf{r}}}_{1},t)\\ &&\quad +{\rm{third-order}}\;{\rm{term}}\end{eqnarray} \tag{ A.7a }$

$\begin{eqnarray}{{\rm{\nabla }}}_{{{\bf{r}}}_{2}}\cdot [{{\bf{C}}}_{+}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)+{{\bf{C}}}_{-}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)] & = & 0,\\ {{\rm{\nabla }}}_{{{\bf{r}}}_{1}}\cdot [{{\bf{C}}}_{+}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)-{{\bf{C}}}_{-}^{u}({{\bf{r}}}_{1},{{\bf{r}}}_{2},t)] & = & 0.\end{eqnarray} \tag{ A.7b }$

Truncation of the CE2 at second order (CE2) consists of equations (A.3) and (A.7), neglecting the third-order term. Clark and Spitz (1995) have derived the evolution equation of the single-time two-point covariance tensor for the compressible Navier–Stokes equation, yielding a more general version of (A.7).

Appendix B.: Eddy absorption in the QL approximation

To illustrate more clearly why CE2 fails to capture the absorption of larger-amplitude eddies, we perform QL simulations of the barotropic vorticity equation (24), eliminating eddy–eddy interactions as in the QL approximation (22) of the Boussinesq equations. The relative vorticity in the QL simulation is shown in figure B1 . The positive vorticity anomaly (labeled V) that detaches from a wave crest is initially sheared by the mean flow, leading to momentum fluxes that strengthen the mean flow and to a decrease of the eddy kinetic energy (Orr mechanism). It is then advected around the center of the cats' eye (labeled X). But instead of rolling up into a small-scale filament as it does in the fully nonlinear simulation (see figure 4), it moves to the western side of the eye, where it joins the wave lobe with positive vorticity west of the one from where it originated ( $T=5.9$ ). An analogous description applies to the negative vorticity anomaly inside the eye. This leads at $T=5.9$ to vorticity anomalies that have a southeast–northwest tilt of phase lines, consistent with an equatorward eddy momentum flux (reflection phase) and an increase of EKE because the mean-flow shear goes against the tilt (Orr mechanism). The vorticity map after two cycles of absorption and reflection ( $T=17.5$ ) is very similar to the initial one (T = 4). Differences arise from wave absorption occurring through weakly nonlinear processes and, to a lesser extent, from hyperviscosity. This is in sharp contrast to the full simulation, in which filamentation takes place, eventually leading to irreversibility (see T = 4 and $T=17.5$ of figure 4).

**Figure B1.** Evolution of relative vorticity in a QL simulation of the larger-amplitude eddies ( $\epsilon =6$ ) and an Earth-like setting ( $\mathrm{Ro}=0.06$ ), for which the fully nonlinear simulation is shown in figure 4. White X's mark the centers of what would become cats' eyes in the fully nonlinear simulation. The locations labeled by V follow a positive vorticity anomaly that detaches from an initial wave lobe.
Download figure:
Standard image High-resolution image

Cumulant expansions for atmospheric flows

Article metrics

Author affiliations

Dates

Abstract

1. Introduction

2. CE2 of Boussinesq flow

2.1. Boussinesq flow

2.2. Averaging operator

2.3. Cumulant expansion

First cumulant

Second cumulant

CE2 equations

Physical properties

2.4. QL approximation

2.5. Differences between CE2 and S3T

3. Dry convective boundary layer

3.1. Large-eddy simulations

Setup

3.2. Results

3.3. Implications

4. Large-scale eddy decay on the rotating sphere

4.1. Barotropic model for the upper troposphere

Equations of motion

Eddy-mean flow decomposition

Cumulants

Numerical implementation

4.2. Eddy lifecycle simulations

Setup

4.3. Results

Varying eddy amplitudes

Relative amplitude of terms in the potential vorticity budget

Varying Rossby numbers

Higher-order closures

4.4. Implications

5. Conclusions

Acknowledgments

Appendix A.: CE2 of anelastic flow

Appendix B.: Eddy absorption in the QL approximation

Footnotes

Cumulant expansions for atmospheric flows

Article metrics

Share this article

Author affiliations

Dates

Abstract

1. Introduction

2. CE2 of Boussinesq flow

2.1. Boussinesq flow

2.2. Averaging operator

2.3. Cumulant expansion

First cumulant

Second cumulant

CE2 equations

Physical properties

2.4. QL approximation

2.5. Differences between CE2 and S3T

3. Dry convective boundary layer

3.1. Large-eddy simulations

Setup

3.2. Results

3.3. Implications

4. Large-scale eddy decay on the rotating sphere

4.1. Barotropic model for the upper troposphere

Equations of motion

Eddy-mean flow decomposition

Cumulants

Numerical implementation

4.2. Eddy lifecycle simulations

Setup

4.3. Results

Varying eddy amplitudes

Relative amplitude of terms in the potential vorticity budget

Varying Rossby numbers

Higher-order closures

4.4. Implications

5. Conclusions

Acknowledgments

Appendix A.: CE2 of anelastic flow

Appendix B.: Eddy absorption in the QL approximation

Footnotes