Cumulant expansions for atmospheric flows

The equations governing atmospheric flows are nonlinear. 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 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. We review how truncations at second order can be used to model and elucidate the dynamics of 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 even if it does not capture important turbulent transport terms. Second, we study the evolution of two-dimensional large-scale waves representing waves 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 nonlinear eddy--eddy interactions in surf zones 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 momentum fluxes in the troposphere.


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, 1970Orszag, , 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 threedimensional 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 andIoannou, 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 andIoannou, 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.

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 ρ(r, t) = ρ 0 + δρ(r, t) as a sum of a constant density ρ 0 in a reference state and fluctuations δρ(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] ∂u ∂t Here, u denotes the three-dimensional velocity, δp the pressure perturbation associated with the density perturbation δρ, Φ = δp/ρ 0 the potential of the pressure-gradient accelerations, and b = −gδρ/ρ 0 the buoyancy (g is an effective gravitational acceleration and k is the local vertical). The reference frame rotates with the constant angular velocity Ω of the planetary rotation, as a result of which Coriolis accelerations (2Ω×u) appear in the momentum equation. Centrifugal accelerations are subsumed in g, the effective gravitational acceleration. The terms F u and 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 u 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 of the flow that contains all prognostic variables in Cartesian coordinates (the superscript (·) T indicates the transpose). The set of equations (1) can then be written compactly as where the outer product (⊗) of the two vectors Ψ and u is defined as The components of the divergence of the second-order term ∇ · (Ψ ⊗ u) are 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 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 F contains all other nonconservative terms. We have expanded u, 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 (u · ∇)u and u · ∇b.

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 f (r, t) into a mean and a fluctuating part, 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 (r, t) and g(r, t) and any constant c, the Reynolds properties [Monin and Yaglom, 1971]: Properties (7a-c) imply that the averaging operator is a projection operator and so is idempotentf =f . They make it possible to define the Reynolds decomposition 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 Appendix A.

Cumulant expansion
First cumulant The first cumulant is the mean Ψ(r, t), for which the equations of motion are obtained by averaging the equations of motion (3): This involves the covariance ∇·(Ψ ⊗ u ), 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.
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 r 1 and r 2 , which need not be equal. The covariance tensor is symmetric, We additionally define the auxiliary covariances, The first, C Φ , contains additional information about covariation of the prognostic fields Ψ with the pressure potential Φ. These covariances can be calculated from C and Ψ because the pressure potential Φ is a diagnostic variable in the Boussinesq approximation. The second, C u , represents covariation of the velocity field with other prognostic variables and is already contained in C.
The second cumulant equation can be obtained from the equations of motion (3) by evaluating Discarding terms that are third order in fluctuating quantities, one obtains where {r 1 ←→ r 2 } indicates the terms obtained by interchanging r 1 and r 2 , which are necessary to ensure the symmetry (11) of the covariance tensor. The third-order term C(r 1 , r 2 , t) ⊗ū(r 1 , t) is defined as in Cartesian coordinates, and its divergence by The flux C(r 1 , r 2 , t) ⊗ū(r 1 , t) represents the transport of spatial eddy correlations by the mean flow at r 1 . The term −C u (r 1 , r 2 , t) · ∇ Ψ (r 2 , t) T represents generation of covariance C(r 1 , r 2 , t) by advection down mean-flow gradients. Additionally, continuity (1b) implies that 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 This set of equations involves the 16 terms in C and the eight covariances C Φ (r 1 , r 2 , t) and C Φ (r 2 , r 1 , t). The 16 components of the second cumulant C 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 r 1 or r 2 of the equation for C u (r 1 , r 2 , t) 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 . 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].

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), Hence, the QL approximation amounts to replacing in the Reynolds decomposition of the nonlinear term, the eddy-eddy interaction Ψ ⊗ u by its mean effect Ψ ⊗ u . Under the QL approximation, the Boussinesq equations (3) can then be written as 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.

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 andIoannou, 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.

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

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 Appendix A.
Our anelastic LES code uses the specific entropy s and the three-dimensional velocity field u 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 ν = 1.2 m 2 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 [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 6.4 km×6.4 km in the horizontal and 3.75 km in the vertical, with a horizontal and vertical resolution of 25 m. Horizontal boundary conditions are doubly periodic. At the upper boundary, flow fields are linearly relaxed over a 800 m deep layer toward the horizontal mean flow, which is almost motionless (horizontal mean velocities are of the order 0.1 m s −1 ). The relaxation coefficient varies from τ = 0 at the bottom of the layer to τ = (100 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 K km −1 . Here, the potential temperature is related to the specific entropy by θ =T exp[(s −s)/c p ], where c p is the specific heat at constant pressure,T is a standard temperature, ands is a standard specific entropy at the standard temperatureT and standard pressurep = 1000 hPa. 2 The initial state is destabilized by imposing a constant sensible heat (enthalpy) flux of 70.46 W m −2 at the lower boundary. Normally distributed random fluctuations of the potential temperature with amplitude 0.1 K in the lowest 200 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 m 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 Eq. (6). The QL truncation (22), here with the state vector Ψ = (u, v, w, s) T , is implemented by removing at every time step the nonlinear eddy-eddy interactions from the tendencies of all prognostic variables. 3 Herring [1963] similarly used the QL approximation to study thermal convection between two parallel horizontal plates.

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 eddyeddy interactions in the horizontal plane play a role in generating the smaller scales of the cellular patterns. The vertical velocity fluctuations w 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  Figure 3. Turbulence kinetic energy budget of the fully nonlinear (a) and QL simulations (b). The different terms denote the TKE advection A, TKE transport by eddies T , shear production S, buoyancy production B, pressure correlation term P, dissipation to subgrid-scales D and rate of change of TKE R. Also shown are the residuals, which are of the same order as the SGS diffusion term in the fully nonlinear simulation [cf. Heinze et al., 2015] but are smaller in the QL simulation, likely because the latter is smoother, thus limiting numerical diffusion.
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 2 + v 2 + w 2 .
Its mean budget is derived from the momentum equations and reads The right-hand side contains all terms that generate, destroy, or redistribute TKE: advection by the mean flow A and the eddies T , shear production S, buoyancy production B (b = gθ /θ 0 denotes buoyancy fluctuations), the pressure correlation term 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 T = −1/ρ 0 ∂ z ρ 0 w 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.
The triple correlation transport term 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 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.

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 1stand 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.

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

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, 1975Held, , 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 Ω. We further assume that the flow is twodimensional on the sphere: u · e r = 0, with e r being the radial coordinate. In that case, the momentum and continuity equations (1a, 1b) reduce to ∂v ∂t The horizontal components of the wind and of the forcing are denoted v and F v . Because the flow v is two-dimensional on the sphere, only the local vertical component Ω sin φ e r (latitude φ) of Ω yields non-zero terms in (24a). Taking the curl of the momentum equation and projecting it onto the radial direction e r yields the twodimensional barotropic vorticity equation, The flow is entirely described by this equation for the absolute vorticity where the Coriolis parameter f (φ) = 2Ω sin φ represents the vorticity of solid body rotation, and ζ = (∇ × v) · e r is the relative vorticity in the radial direction e r , relative to the rotating reference frame. The advection term v · ∇q contains the advection of planetary vorticity v · ∇f = βv, with β = a −1 ∂ φ f = 2Ωa −1 cos φ, 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 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 a as the typical length scale and the length of the day 2πΩ −1 as the typical time scale. With that nondimensionalisation, the operators ∇, ∇× and ∇ 2 become operators on the unit sphere, and the angular velocity Ω becomes 2π. 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, Here β = 2Ωcosφ, with Ω = 2π, is the nondimensional meridional derivative of the Coriolis parameter, and the Jacobian represents the advection on the unit sphere of vorticity ζ by the zonal (u) and meridional (v) flow implied by the streamfunction ψ: The streamfunction-vorticity relations are with ∇ 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 −J n (ψ , ζ ). The eddy vorticity budget involves shear by the mean flowū/ cos φ ∂ λ ζ , eddy-eddy interactions J n (ψ , ζ ) − J n (ψ , ζ ) , the β-term βv and the advection of mean-flow vorticity Cumulants The first cumulant is the mean vorticityζ(φ, t), and the second cumulant is the vorticity equal-time two-point correlation: The first cumulant depends on the latitude φ, and the second cumulant depends on the latitudes φ 1 and φ 2 and, because of the statistical axisymmetry of the equations, on the longitude difference λ 1 − λ 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 ζ (r 1 , t) ⊗ u (r 2 , t) or u (r 1 , t) ⊗ u (r 2 , t) can be calculated fromζ 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. [2008Marston et al. [ , 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,b;Qi and Marston, 2014] with 163, 842 cells. To remove enstrophy that cascades to the smallest scales, hyperviscous dissipation ν(∇ 2 + 2)∇ 6 ζ is included, where the term (∇ 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 ∆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 ≤ ≤ L on spherical wavenumber , with the zonal wavenumbers restricted to |m| ≤ min{ , M }. We choose spectral cutoffs L = 60 and M = 30. Hyperviscosity is adjusted to ensure the same e-folding time at the maximum wavenumber = 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 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. 4 More details about the simulations and the cumulant expansions can be found in Marston et al. [2014].

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 (φ) with a superimposed initial disturbance (eddy) with vorticity ζ i (φ, λ). The zonal flow U and disturbance ζ 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 To mimic Earth's upper troposphere, we choose 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 ∼33 m s −1 at a latitude of ∼40 • and a maximum eastward velocity of ∼22 m s −1 at the equator; its zero is near ∼17 • . This zonal flow is barotropically stable. The eddy vorticity ζ i decays meridionally away from its maximum absolute value ζ 0 cos φ m at latitude φ m = π/4 = 45 • , with characteristic meridional decay scale δ = π/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., 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 2Ro of the initial zonal flow U . Dimensionally, we have Ro = Λ/(2Ω), where Λ ≈ 2 max(U )/a is the typical initial zonal-mean flow vorticity. Hence, Ro is a Rossby number measuring the mean flow vorticity to the equatorial planetary vorticity 2Ω. We measure the initial maximum vorticity ζ 0 of the eddies relative to the mean-flow vorticity 2Ro through the amplitude parameter = ζ 0 cos φ m /2Ro. The quantities ζ, ζ , ψ, ψ , and t are then rescaled with 4πRo, 4π Ro, 4πRo, 4π Ro, and (4πRo) −1 , respectively. The equations of motion for the mean-flow and the eddies become 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 Ro (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 Ro −1 larger. The nondimensional parameters Ro and control the vorticity of the mean flow and of the eddies and are important for the evolution of the barotropic flow. Eddymean flow interactions are of order Ro, and eddy-eddy interactions of order Ro, provided ∇ 2 n 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 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.

Results
Varying eddy amplitudes For a fixed initial mean-flow Rossby number Ro ≈ 0.06 (corresponding to the Earth-like parameters in equation (34)), we run eddy lifecycle experiments for larger-amplitude initial eddies with ≈ 6, and for smaller-amplitude initial eddies with ≈ 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 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, 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 k i = 6. Such cats' eyes are characteristic of Rossby waves that break in "surf zones" around their critical layers 5 [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 E K = 0.5 φ e K cos φ dφ, where e K = 0.5(u 2 + v 2 ), becomes very small at large times ( 30, 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 e Z = 0.5u 2 (ū ≡ −∂ φψ 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 (Fig. 5c). Total energy E K + E Z , 5 Fig. 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. with E Z = 0.5 φ e Z cos φ dφ, is approximately conserved in the model, up to very small losses (∼0.2% 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 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 (F v = 0): Acceleration of the mean zonal flow occurs where eddy momentum fluxes u v cos φ converge. Multiplying the mean zonal momentum equation (36) byū and integrating by parts yields the equation for the zonal kinetic energy E Z , 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ū cos −1 φ. The acceleration of the mean zonal jet at its core (Fig. 5c) thus is associated with eddy momentum flux convergence ∂ φ (u v cos φ) (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 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 (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 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 ( ≈ 2). 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 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 (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 eddymean 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 eddyeddy interactions seems consistent with the order of magnitude of the terms in the vorticity equations (35). The eddy-eddy interactions appear of order Ro ≈ 0.4, compared with the β-term which is responsible for Rossby wave dynamics and is of For the smaller-amplitude experiment, the dimensional analysis of the vorticity equations (35) suggests that the eddy-eddy interactions now are of order Ro 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 Ro ≈ 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 ( ≈ 6) and decrease the Rossby number 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 Ro = 0.06 is not captured by CE2 because it is strongly nonlinear. However, at the smaller Rossby number Ro = 0.02, it is faithfully captured by CE2. Indeed, eddy absorption appears more weakly nonlinear for Ro = 0.02, and dominated by meanflow shearing, as is evident on the relative vorticity maps (Fig. 9). Vorticity maps for Ro = 0.04 and 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 (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 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 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.
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 M = 40 and L = 20. 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 Figure 8. Evolution of eddy kinetic energies e K in the fully nonlinear simulation (left column) and in a direct CE2 calculation of e K (right column) for the largeramplitude eddies ( = 6), with Rossby number decreasing from Ro = 0.06 to 0.02. Time scales are non-dimensionalized with the length of the day 2πΩ −1 and length scales with the planet radius a. 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 ; they bring similar improvements.

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 Ro. The parameter need not be small for absorption through weakly nonlinear shearing to occur. For example, for small 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].

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 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 eddyeddy interactions, such as those captured by the third order CE3* approximation .
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, 1970Orszag, , 1973Marston 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 B. Eddy absorption in the QL approximation
To illustrate more clearly why CE2 fails to capture the absorption of largeramplitude 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 Fig. 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 centre of the cats' eye (labeled X). But instead of rolling up into a smallscale filament as it does in the fully nonlinear simulation (cf. Fig. 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 meanflow 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 (cf. T = 4 and T = 17.5 of Fig. 4). Figure B1. Evolution of relative vorticity in a QL simulation of the largeramplitude eddies ( = 6) and an Earth-like setting (Ro = 0.06), for which the fully nonlinear simulation is shown in Fig. 4. White X's mark the centres 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.