# Cumulant expansions for atmospheric flows

###### Abstract

Atmospheric flows are governed by the equations of fluid dynamics. These equations are nonlinear, and consequently the hierarchy of cumulant equations is not closed. But because atmospheric flows are inhomogeneous and anisotropic, the nonlinearity may manifest itself only weakly through interactions of non-trivial mean fields with disturbances such as thermals or eddies. In such situations, truncations of the hierarchy of cumulant equations hold promise as a closure strategy.

Here we show how truncations at second order can be used to model and elucidate the dynamics of turbulent atmospheric flows. Two examples are considered. First, we study the growth of a dry convective boundary layer, which is heated from below, leading to turbulent upward energy transport and growth of the boundary layer. We demonstrate that a quasilinear truncation of the equations of motion, in which interactions of disturbances among each other are neglected but interactions with mean fields are taken into account, can capture the growth of the convective boundary layer. However, it does not capture important turbulent transport terms in the turbulence kinetic energy budget. Second, we study the evolution of two-dimensional large-scale waves, which are representative of waves seen in Earth’s upper atmosphere. We demonstrate that a cumulant expansion truncated at second order (CE2) can capture the evolution of such waves and their nonlinear interaction with the mean flow in some circumstances, for example, when the wave amplitude is small enough or the planetary rotation rate is large enough. However, CE2 fails to capture the flow evolution when strongly nonlinear eddy–eddy interactions that generate small-scale filaments in surf zones around critical layers become important. Higher-order closures can capture these missing interactions.

The results point to new ways in which the dynamics of turbulent boundary layers may be represented in climate models, and they illustrate different classes of nonlinear processes that can control wave dissipation and angular momentum fluxes in the upper troposphere.

^{†}

^{†}: New J. Phys.

Keywords: turbulence, closure, quasi-linear approximation, atmospheric boundary layer, atmospheric convection, large-scale atmospheric circulation, jets

## 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., 2016]. 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, 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 non-trivial mean fields, which include, for example, atmospheric jet streams or the thermal stratification of the atmosphere. Because interactions between turbulent fluctuations and non-trivial 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; Farrell and Ioannou, 1993; Randel and Held, 1991; 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 non-trivial 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; Tobias and Marston, 2013; Bakas and Ioannou, 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 Cumulant expansion 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 as a sum of a constant density in a reference state and fluctuations 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]

(momentum equation) | (1a) | ||||

(continuity equation) | (1b) | ||||

(thermodynamic equation) | (1c) |

Here, denotes the three-dimensional velocity, the pressure perturbation associated with the density perturbation , the potential of the pressure-gradient accelerations, and the buoyancy ( is an effective gravitational acceleration and is the local vertical). The reference frame rotates with the constant angular velocity of the planetary rotation, as a result of which Coriolis accelerations () appear in the momentum equation. Centrifugal accelerations are subsumed in , the effective gravitational acceleration. The terms and 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 is non-divergent. 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 non-divergence condition to eliminate the time derivative.

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

(2) |

of the flow that contains all prognostic variables in Cartesian coordinates (the superscript indicates the transpose). The set of equations (1) can then be written compactly as

(3a) | ||||

(3b) |

where the outer product () of the two vectors and is defined as

(4) |

The components of the divergence of the second-order term are

(5) |

where summation over repeated indices is implied. This notation extends the usual advection operator of a scalar field by a non-divergent flow to the advection of a vector field. The linear operator 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 contains all other nonconservative terms. We have expanded , the outer product and the 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 cumulant expansion approaches because they capture the essential nonlinearity of atmospheric flows: the conservative quadratic nonlinearity of the advection operators and .

### 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 , and decompose any scalar field into a mean and a fluctuating part,

(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 and and any constant , the Reynolds properties [Monin and Yaglom, 1971]:

(7a) | ||||

(linearity) | (7b) | |||

(7c) | ||||

(commutation with derivatives) | (7d) |

Properties (7a-c) imply that the averaging operator is a projection operator and so is idempotent . They make it possible to define the Reynolds decomposition

(8) |

For a vector quantity , the average 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 A.

### 2.3 Cumulant expansion

#### First cumulant

The first cumulant is the mean , for which the equations of motion are obtained by averaging the equations of motion (3):

(9a) | ||||

(9b) |

This involves the covariance , 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

(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 and , which need not be equal. The covariance tensor is symmetric,

(11) |

We additionally define the auxiliary covariances,

(12a) | ||||

(12b) |

The first, , contains additional information about covariation of the prognostic fields with the pressure potential . These covariances can be calculated from and because the pressure potential is a diagnostic variable in the Boussinesq approximation. The second, , represents covariation of the velocity field with other prognostic variables and is already contained in .

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

(13) |

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

(14) |

where indicates the terms obtained by interchanging and , which are necessary to ensure the symmetry (11) of the covariance tensor. The third-order term is defined as in Cartesian coordinates,

(15) |

and its divergence by

(16) |

The flux represents the transport of spatial eddy correlations by the mean flow at . The term represents generation of covariance by advection down mean-flow gradients.

Additionally, continuity (1b) implies that

(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 cumulant expansion (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

(18a) | ||||

(18b) | ||||

(18c) | ||||

(18d) |

This set of equations involves the 16 terms in and the eight covariances and . The 16 components of the second cumulant are prognostic (evoloving according eq. (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 or of the equation for that can be extracted from (18c), and in which time tendencies vanish because of the non-divergence 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, b) 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, d), 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 Quasi-linear 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 quasi-linear (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 . However it linearizes the equation for the eddies , obtained by substracting (9a) from (3a),

(19) |

by setting

(20) |

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

(21) |

the eddy–eddy interaction by its mean effect . Under the QL approximation, the Boussinesq equations (3) can then be written as

(22a) | ||||

(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

Stochastic Structural Stability Theory (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 non-zonal [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 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; Webb et al., 2013; Stevens and Bony, 2013; Vial et al., 2013; Brient et al., 2016].

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 subgrid-scale 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 [e.g., Mellor and Yamada, 1982; Beljaars, 1992]. Turbulent fluxes in convective clouds, on the other hand, are often represented as vertically non-local 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 subgrid-scale 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 non-hydrostatic and filter sound waves by neglecting dynamic density variations except where they affect buoyancy. 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 cumulant expansion are described in A.

Our anelastic LES code uses the specific entropy and the three-dimensional velocity field as prognostic variables [Pauluis, 2008]. We use a second-order central difference spatial discretization scheme with strongly stability preserving Runge-Kutta timestepping [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, subgrid-scale (SGS) motions must also be modeled, and we do so by applying a constant eddy diffusivity of 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 [Smagorinsky, 1963; Lilly, 1962], 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 in the horizontal and in the vertical, with a horizontal and vertical resolution of . Horizontal boundary conditions are doubly periodic. At the upper boundary, flow fields are linearly relaxed over a deep layer toward the horizontal mean flow, which is almost motionless (horizontal mean velocities are of the order ). The relaxation coefficient varies from at the bottom of the layer to at the top.

The initial state is stably stratified with a potential temperature that increases linearly with height, at a rate of . Here, the potential temperature is related to the specific entropy by , where is the specific heat at constant pressure, is a standard temperature, and is a standard specific entropy at the standard temperature and standard pressure .^{1}^{1}1Using the ideal-gas law, it can be verified that is the usual potential temperature with reference pressure : , with specific gas constant . However, this potential temperature is evaluated at the anelastic reference-state pressure rather than at the in situ pressure , as is required for thermodynamic consistency of the anelastic approximation [Pauluis, 2008]. The initial state is destabilized by imposing a constant sensible heat (enthalpy) flux of at the lower boundary. Normally distributed random fluctuations of the potential temperature with amplitude in the lowest 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 . 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 Eq. (6). The QL truncation (22), here with the state vector , is implemented by removing at every time step the nonlinear eddy–eddy interactions from the tendencies of all prognostic variables.^{2}^{2}2The QL truncation (22) is valid for the Boussinesq approximation. The anelastic approximation requires us to use an average weighted by the background density, as explained in A. But because averages are here performed on horizontal surfaces with constant background density, the QL approximation for the anelastic system is also given by (22). 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 (Fig. 1a,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 (Fig. 1c). 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 (Fig. 1d), 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 are distributed more symmetrically around zero in the QL simulation: updrafts and downdrafts are of similar size and strength (Fig. 1b, d). This stands in contrast to the fully nonlinear simulations, in which updrafts are faster and narrower and downdrafts slower and broader.

Figure 2a 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 (Fig. 2a). 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.

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 (Fig. 2b), 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: . Its mean budget is derived from the momentum equations and reads \linenomath

(23) |

The right-hand side contains all terms that generate, destroy, or redistribute TKE: advection by the mean flow and the eddies , shear production , buoyancy production ( denotes buoyancy fluctuations), the pressure correlation term 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 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.

The triple correlation transport term 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 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 cumulant expansions.

### 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 Jr 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 rotating at constant spin angular velocity . We further assume that the flow is two-dimensional on the sphere: , with being the radial coordinate. In that case, the momentum and continuity equations (1a, 1b) reduce to

(24a) | ||||

(24b) |

The horizontal components of the wind and of the forcing are denoted and . Because the flow is two-dimensional on the sphere, only the local vertical component (latitude ) of yields non-zero terms in (24a). Taking the curl of the momentum equation and projecting it onto the radial direction yields the two-dimensional barotropic vorticity equation,

(25) |

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

(26) |

where the Coriolis parameter represents the vorticity of solid body rotation, and is the relative vorticity in the radial direction , relative to the rotating reference frame. The advection term contains the advection of planetary vorticity , with , 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

(27a) | ||||

(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 nondimensionalised using the planetary radius as the typical length scale and the length of the day as the typical time scale. With that nondimensionalisation, the operators , and become operators on the unit sphere, and the angular velocity becomes . 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 then suggests itself. Decomposing flow fields in the nondimensional barotropic vorticity equation into zonal means and eddies yields the mean and eddy vorticity equations,

(28a) | ||||

(28b) |

Here , with , is the nondimensional meridional derivative of the Coriolis parameter, and the Jacobian

(29) |

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

(30) |

The streamfunction-vorticity relations are

(31a) | ||||

(31b) |

with likewise defined on the unit sphere.

#### Cumulants

The first cumulant is the mean vorticity , and the second cumulant is the vorticity equal-time two-point correlation:

(32a) | ||||

(32b) |

The first cumulant depends on the latitude , and the second cumulant depends on the latitudes and and, because of the statistical axisymmetry of the equations, on the longitude difference [Marston et al., 2008]. Because the flow is entirely determined by the scalar (or ), all other correlations are determined by the scalar cumulant . Hence, moments like or can be calculated from and [e.g., Marston et al., 2008; Srinivasan and Young, 2012; Marston et al., 2014].

#### 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, b; Qi and Marston, 2014] with cells. To remove enstrophy that cascades to the smallest scales, hyperviscous dissipation is included, where the term ensures that angular momentum is conserved. The hyperviscosity coefficient is chosen such that the smallest resolved mode decays with an -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 .

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 on spherical wavenumber , with the zonal wavenumbers restricted to . We choose spectral cutoffs and . Hyperviscosity is adjusted to ensure the same -folding time at the maximum wavenumber 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 non-zero 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.^{3}^{3}3The application “GCM” is available for OS X 10.9 and higher on the Apple Mac App Store at URL http://appstore.com/mac/gcm More details about the simulations and the cumulant expansions 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 with a superimposed initial disturbance (eddy) with vorticity . The zonal flow and disturbance 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

(33a) | ||||

(33b) |

To mimic Earth’s upper troposphere, we choose

(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 at a latitude of and a maximum eastward velocity of at the equator; its zero is near . This zonal flow is barotropically stable. The eddy vorticity decays meridionally away from its maximum absolute value at latitude , with characteristic meridional decay scale . Its zonal wavenumber 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 of the initial zonal flow . Dimensionally, we have , where is the typical initial zonal-mean flow vorticity. Hence, is a Rossby number measuring the mean flow vorticity to the equatorial planetary vorticity . We measure the initial maximum vorticity of the eddies relative to the mean-flow vorticity through the amplitude parameter . The quantities , , , , and are then rescaled with , , , , and , respectively. The equations of motion for the mean-flow and the eddies become

(35a) | ||||

(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 (the case we will consider) is that the advection of mean-flow vorticity by meridional velocity fluctuations is negligible compared with the -term, which is a factor larger.

The nondimensional parameters 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 , and eddy-eddy interactions of order , provided 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 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 (corresponding to the Earth-like parameters in equation (34)), we run eddy lifecycle experiments for larger-amplitude initial eddies with , and for smaller-amplitude initial eddies with . 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 ) 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, Fig. 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 Fig. 4) that continue to have the initial zonal wavenumber . Such cats’ eyes are characteristic of Rossby waves that break in “surf zones” around their critical layers^{4}^{4}4Fig. 4e shows some spurious oscillations in the critical layer region. This is the consequence of a too weak hyperdiffusion. With hyperdiffusion strong enough to remove the ripples, we observed have noticable hyperdiffusive wave absorption over the time scales considered. This obscures the wave absorption due to eddy–eddy interactions and blurs the differences between fully nonlinear and CE2 simulations. Removing the spurious ripples while showing sharp differences between fully nonlinear and CE2 simulations would have required a much higher resolution, or perhaps a different numerical advection scheme. [Stewartson, 1977; Warn and Warn, 1978; McIntyre and Palmer, 1983; Killworth and McIntyre, 1985; Held and Phillips, 1987].

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 , where , becomes very small at large times (, see Fig. 5a). At those times, most of the initial kinetic energy has been transferred to the mean zonal flow. The local zonal kinetic energy ( of the mean zonal flow increases in the core of the midlatitude jet, roughly in proportion to the decrease of the eddy kinetic energy (Fig. 5c). Total energy , with , is approximately conserved in the model, up to very small losses ( of the total after a time of 50) through subgrid-scale 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.

The transfer of to 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 ():

(36) |

Acceleration of the mean zonal flow occurs where eddy momentum fluxes converge. Multiplying the mean zonal momentum equation (36) by and integrating by parts yields the equation for the zonal kinetic energy ,

(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 . The acceleration of the mean zonal jet at its core (Fig. 5c) thus is associated with eddy momentum flux convergence (EMFC, see Fig. 5e).

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 (Fig. 5e). 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 (Fig. 5a, b). However, the oscillations are too weakly damped; large eddy kinetic energies 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 (Fig. 5f). How these oscillations arise from the perspective of wave mechanics, and why CE2 cannot capture the wave absorption in this case, is illustrated in 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 (). Cats’ eyes with rolling-up vorticity filaments do not develop in the fully nonlinear simulation (Fig. 6). Instead, eddies are sheared by the mean flow, which transfers eddy kinetic energy 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 (Fig. 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 (Fig. 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.

It is worth noting that eddies are also sheared equatorward of the cats’ eyes in the larger-amplitude simulation (Fig. 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 , 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 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 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 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 () and decrease the Rossby number from to . 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 in the fully nonlinear and the CE2 simulations. As we have already seen, eddy absorption at is not captured by CE2 because it is strongly nonlinear. However, at the smaller Rossby number , it is faithfully captured by CE2. Indeed, eddy absorption appears more weakly nonlinear for , and dominated by mean-flow shearing, as is evident on the relative vorticity maps (Fig. 9). Vorticity maps for and 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 (Fig. 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 (Fig. 8). This occurs although the orders of magnitude of the relevant terms suggested by the dimensional analysis decrease by equal factors of order 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 , shear explains eddy decay and jet acceleration, even though the nondimensionalization (35) suggests shear should be much smaller than eddy-eddy interactions.

#### 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 Fig. 10. Because CE3* is computationally expensive, the resolution has been reduced to and . For simplicity we turn off eddy damping (, see Marston et al. [2014] for definitions). The full and CE2 simulations in Fig. 10 are run at this lower resolution and are consistent with the higher-resolution runs (Fig. 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 . The parameter need not be small for absorption through weakly nonlinear shearing to occur. For example, for small , 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 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 cumulant expansions capture interactions between fluctuations (e.g., thermals) and mean fields and take non-local 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 parameterisation 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 non-local 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, 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.

## Appendix A Cumulant expansion of anelastic flow

Because the flow is non-divergent 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 , with 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-c) but, importantly, not commutativity with partial differentiation (7d). Thus, we can define an eddy mean flow decomposition

(38a) | ||||

(38b) |

with hats 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 to vary with height . In the anelastic approximation, density perturbations , the pressure potential , and the buoyancy are defined relative to this height-dependent and hydrostatically balanced reference density. The state vector is defined as in (2), but the augmented state vector is now taken to be

(39) |

We have to consider covariances involving and not because multiplying by density does not commute with derivatives. For simplicity, we assume that the linear operator is homogeneous in the reference density and satisfies . It implies that it cannot involve any spatial derivative. The anelastic equations in a reference frame rotating with angular velocity are [Vallis, 2006]:

(40a) | ||||

(40b) |

The Boussinesq equations (1) are obtained from the anelastic equations by setting to a constant. The continuity equation (40b) again reduces to a non-divergence constraint. Because the reference density appears in it, we need to use the reference-density weighted average to derive the cumulant expansion.

The equation for the first cumulant is obtained by averaging (40):

(41a) | |||

(41b) |

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 and , that gives the second-order term in the first cumulant equation (41), and that corresponds to a density weighted average when . A possible choice is

(42) |

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

(43a) | |||

(43b) |

and introduce the corresponding quantities, , , and with the density difference instead of the sum, as in