A model of the entropy flux and Reynolds stress in turbulent convection

We propose a closure model for the transport of entropy and momentum in astrophysical turbulence, intended for application to rotating stellar convective regions. Our closure model is first presented in the Boussinesq formalism, and compared with laboratory and numerical experimental results on Rayleigh-Benard convection and Homogeneous Rayleigh-Benard convection. The predicted angular momentum transport properties of the turbulence in the slowly rotating case recover the well-known Lambda-effect, with an amplitude uniquely related to the convective heat flux. The model is then extended to the anelastic case as well as the fully compressible case. In the special case of spherical symmetry, the predicted radial heat flux is equivalent to that of mixing-length theory. For rotating stars, our model describes the coupled transport of heat and angular momentum, and provides a unified formalism in which to study both differential rotation and thermal inhomogeneities in stellar convection zones.


INTRODUCTION
Turbulent convection occurs frequently in stellar interiors and other astrophysical fluid flows. While convective motion naturally transports heat and chemical elements, the transport of angular momentum by convection in rotating bodies is a more subtle issue. It is of particular interest in the case of the Sun, where the internal pattern of rotation has been measured but remains incompletely understood. It may also play an significant role in accretion flows.
Numerical simulations of astrophysical convection are becoming increasingly powerful and capable of resolving a widening range of length and time scales. Nevertheless, a simpler, statistical description of turbulent transport is desirable in order to treat the effects of convection on the structure and evolution of stars. It almost goes without saying that such a description cannot be derived strictly from the equations of fluid dynamics but must involve some modelling or parametrization.
The mixing-length theory of turbulent transport was developed by Prandtl (1925) and applied to stellar convection by Biermann (1932). It is still the basic model used in most calculations of stellar structure and evolution, usually in the form devised by Böhm-Vitense (1958). The main purpose of mixing-length theory is to relate the convective heat flux to the superadiabatic gradient; in this context it does not usually deal with the transport of (angular) momentum that arises in the presence of shear or rotation.
A standard theoretical approach to convection in differentially rotating stars is set out in the monograph by Rüdiger (1989). Angular momentum transport is described by a Reynolds stress tensor whose components can be related to the large-scale mean flows and thermodynamical gradients. A first contribution to the Reynolds stress is typically proportional to the angular velocity gradient through a turbulent viscosity coefficient. An important additional contribution comes from the Λ-effect (named after Lebedinsky), whereby even uniformly rotating convection transports angular momentum by virtue of its anisotropy. Attempts to constrain or parameterize these quantities have been made through local numerical simulations (e.g. Käpylä, Korpi & Tuominen 2004) or theoretical models (e.g. Kitchatinov & Rüdiger 1993). Mean-field models of stellar rotation (e.g. Kitchatinov & Rüdiger 1999, Rempel 2005 have been developed which use such parameterized expressions for the Reynolds stress and heat flux. Reynolds-stress models of turbulent flows have been developed in the engineering community over several decades (e.g. Pope 2000). The exact equation governing the Reynolds c 0000 RAS stress in a turbulent fluid cannot be solved because of the well known closure problem whereby an infinite hierarchy of correlations is involved. Nevertheless, by parametrizing the difficult terms in this equation, models can be constructed that bear some fidelity to the turbulent dynamics. From a more physical point of view, what is obtained is a time-dependent constitutive equation for the turbulent fluid, which relates the turbulent stress to the local history of deformation. There is a close similarity with models of non-Newtonian fluids (Ogilvie & Proctor 2003). The advection and deformation of the turbulent stress are accurately represented since they derive from linear terms in the Reynoldsstress equation, while the nonlinear 'relaxation' effects are only modelled (as is also true for non-Newtonian fluids).
A similar approach can be applied to turbulent convection in which buoyancy forces play an essential role. The additional correlations that must be considered are the flux and the variance of entropy (or temperature, in the Boussinesq approximation). This approach offers some benefits over the conventional description in terms of a turbulent viscosity and a Λ-effect. It can be formulated in a covariant manner and is not tied to the spherical geometry of a slowly rotating star. It starts from a more fundamental description and allows phenomena such as the Λ−effect to emerge in a natural way from more elementary considerations. It may also allow a more unified approach to be taken towards problems involving astrophysical turbulence.
In this paper we explore some of the consequences of a simple dynamical model of astrophysical convection of this type. The model derives from one originally conceived for magnetohydrodynamic turbulence in accretion discs (Ogilvie 2003) and later applied to rotating shear flows without magnetic fields (Garaud & Ogilvie 2005, GO05 herafter). Our motivation is to develop and test a model that can be applied to the convective zone of the Sun, to other stars or to accretion discs. We emphasize, however, that our model is chosen to be as simple as possible for the purposes of this investigation. In contrast with some of the engineering literature, we restrict the algebraic complexity in order to retain a physical understanding of the terms in the equations. Further refinements are likely to be required in order to provide an accurate match to a wide range of data.
In comparing a closure model of astrophysical convection with experimental and numerical results, we face certain difficulties. Astrophysical convection usually takes place at very high Rayleigh number, in a highly turbulent regime. Experiments have been conducted at very high Rayleigh number but mainly for the Rayleigh-Bénard problem in which the flow is dominated by boundary layers, which may not be relevant in an astrophysical context, or by mean flows not represented in the closure model. An alternative system is provided by the homogeneous Rayleigh-Bénard problem, which has periodic boundary conditions in all directions. This model, however, has certain peculiarities of its own. These issues will be addressed in the sections that follow.
In the remainder of the paper, we develop the closure model first in the Boussinesq approximation (Section 2) and apply it to the standard Rayleigh-Bénard problem (Section 3). We then consider the homogeneous Rayleigh-Bénard system with triply periodic boundary conditions (Section 4); in this section we also introduce rotation and discuss the Λ−effect. We then adapt the model to the anelas-tic approximation for use in stars and other astrophysical flows (Section 5) and finally draw conclusions (Section 6). A number of technical details are covered in the appendices.

Basic equations
In the Boussinesq approximation (e.g. Chandrasekhar 1961) the equations governing the motion of the fluid are where we have adopted a Cartesian tensor notation. The dynamical variables are the velocity u, the density ρ, the pressure p and the temperature T . Quantities regarded as constant in the Boussinesq approximation are the reference density ρ0, the reference temperature T0, the coefficient of expansion α, the gravitational acceleration g, the kinematic viscosity ν, and the thermal diffusivity κ. A simple, static basic state is possible when the temperature is uniform and the pressure gradient balances gravity, i.e.
where p0 is a reference pressure. To examine departures from this state we define obtaining the governing equations

Fluctuations
We now adopt a standard procedure and separate the dynamical variables into mean and fluctuating parts, e.g.
where the angle brackets or the overbar are interchangeably used to denote a suitable averaging operation such as a temporal, spatial or ensemble average. The mean parts of the governing equations are where is the Reynolds tensor, representing (minus) the turbulent stress, and represents the turbulent heat flux density. The problem at hand is to determineRij andFi and thereby close the system of mean equations. We also introduce the quantity representing the temperature variance. It should be noted that all three quadratic correlations Rij , Fi and Q will be redefined when we move on to the (more relevant) anelastic system in which the reference density is non-uniform, but these definitions are convenient for the Boussinesq system. The fluctuating parts of the governing equations are From these we can obtain exact equations forRij ,Fi andQ in the form (∂t +ūj ∂j)Fi +Rij ∂jΘ +Fj∂jūi + αQgi − 1 2 (ν + κ)∂jjFi The left-hand sides of these equations represent the linear interaction ofRij ,Fi andQ with the mean velocity gradient, the mean temperature gradient and the gravitational field, as well as their diffusion by the microscopic transport coefficients. There is no difficulty in treating such terms exactly as they appear. The right-hand sides of these equations contain difficult terms of three sorts: those involving correlations with the pressure fluctuation ψ ′ , those involving triple correlations of fluctuating quantities, and dissipative terms involving the microscopic diffusivities ν and κ. These effects can all be regarded as 'non-linear'; although viscous diffusion, for example, is a linear process, when the Reynolds number is large the viscous terms can be significant only when a turbulent cascade has forced structure to appear on the dissipative scales. None of the terms on the right-hand sides of these equations can be written in terms ofRij ,Fi andQ without further knowledge of the statistical properties of the fluctuating quantities, such as the spectrum of the turbulence, which are determined by the non-linear physics of the turbulent cascade.

Proposed closure model
We therefore attempt to model the system by retaining the exact forms of the left-hand sides and proposing simple closures for the right-hand sides, i.e.
where the quantities F are non-linear tensorial functions of their arguments. The dots represent the parameters of the problem, on which the functions F may depend. A simple example of such a model is (∂t +ūj ∂j)Fi +Rij ∂jΘ +Fj∂jūi + αQgi − 1 2 (ν + κ)∂jjFi where R = Rii is the trace of the Reynolds tensor, which is twice the turbulent kinetic energy per unit mass, and C1, C2, C6 and C7 are positive dimensionless coefficients of order unity, of a universal nature. (Coefficients C3, C4 and C5 are reserved for a magnetohydrodynamic extension of the model, see Ogilvie 2003) The justification for introducing non-linear terms of the above form is similar to that used in the model of magnetorotational turbulent stresses originally introduced by Ogilvie (2003). The term involving C1 causes a dissipation of turbulent kinetic energy, and allows for the free decay of hydrodynamic turbulence. The term involving C2 redistributes energy among the components ofRij, and corresponds to the tendency of hydrodynamic turbulence to return to isotropy through the effect of the pressure-strain correlation. Both are constructed assuming that these effects occur on a timescale related to the eddy turnover time, L/R 1/2 , where L is defined as the typical scale of the largest turbulent eddies. Terms C6 and C7, related to the transport of heat, are advanced by simple analogy. The coefficients must satisfy certain conditions to ensure the realizability of the model, as discussed in Appendix A.
The terms proportional to the microscopic diffusion coefficients are introduced to allow a modelling of the correlation terms 2ν ∂ k u ′ i ∂ k u ′ j , (ν + κ) ∂ju ′ i ∂j Θ ′ and 2κ (∂iΘ ′ ) 2 at moderate Reynolds number, i.e. close to the onset of convection. In such a situation a turbulent cascade does not form and the dissipative terms are proportional to, rather than independent of, the diffusion coefficients. In a similar way, for turbulent shear flows, GO05 proposed to model the momentum diffusion term as on dimensional grounds. Indeed, it is expected that near the onset of convection, most fluid motions will be on the largest scales of the system (L). By analogy, we model the other two terms here as Therefore the dissipative term in each of equations (25)- (27) is modelled by a sum of two terms, one that is independent of the diffusivity and dominates at high Reynolds numbers, and another that is proportional to the diffusivity and dominates at moderate Reynolds numbers. This completes the justification for the form of the closure model proposed in equations (28), (29) and (30).

Model setup
We now apply the closure model to the problem of Rayleigh-Bénard convection. We consider a horizontally infinite, plane-parallel system, where the bottom plate is located at height z = 0 and the top plate at height z = h. The relative temperature of the bottom plate isΘ = ∆T while that of the top plate isΘ = 0. In this setup, we look for statistically steady and horizontally homogeneous solutions assuming that mean quantities and correlations between fluctuating quantities vary only with z. We also assume that there are no mean flows in the system. Equations (13)-(15) and (28)-(30) reduce to a set of ordinary differential equations (ODEs) which can be solved to obtain the temperature profileΘ(z) between the two plates, the profiles of the turbulent kinetic energy, R(z)/2, and the temperature variance,Q(z) (for example).
By analogy with Prandtl's mixing-length formulation (Prandtl, 1932) we set L, the size of the largest eddies, to be equal to the distance to the nearest wall, i.e. L(z) = min(z, h−z) (see GO05 for applications of the same principle to pipe flows and to Couette-Taylor flows).
It can be shown with little effort thatRxy =Rxz = Ryz = 0, as well asFx =Fy = 0. The remaining set of five second-order ODEs fully characterizes the system: where g = −gz. In the case of no-slip boundaries with fixed temperature on each plate as listed above,R,Rzz,Fz and Q are zero on both boundaries. This system of ODEs with associated boundary conditions can be solved with a two-point boundary-value solver. Typical solutions are shown in Fig. 1 for various Rayleigh numbers, defined here as We set the Prandtl number to 1 for the purposes of illustration. Note the appearance of the characteristically flat temperature profile between the two plates as Ra → ∞ and of the thin thermal boundary layers. We now study in more detail the structure of the solution.

Universal profile of convection from a wall
As in the case of shear flows past a wall (see GO05), we can derive a universal profile for convection away from a wall. Let us consider a semi-infinite domain z > 0, in which case L = z, and let F0 be the convective heat flux through the system. We define dimensionless variables via so that the system of equations (34) becomes Pr + 1 2 f ′′ = Pr + 1 2 Cνκ The boundary conditions at η = 0 are r = rzz = f = q = θ = 0. Solutions very close to the wall (η ≪ 1) satisfy: r and rzz ∝ η αν with αν (αν − 1) = Cν , f ∝ η ανκ with ανκ(ανκ − 1) = Cνκ, q ∝ η ακ with ακ(ακ − 1) = Cκ.
These simple relationships provide an ideal way of calibrating each of the three constants Cν, Cνκ and Cκ individually (see Section 3.4), by analysing the power-law behaviour of the near-wall profiles of experimental or numerical data. Solutions far away from the boundary layer can be expanded as where However, unlike r0, f1 and q0 the constant θ0 cannot be determined without a numerical calculation of the boundarylayer solution for η = O(1). The scaling laws obtained for r, f , θ and q far from the wall are expected on dimensional grounds, and recover the well-known solution of Priestley (1954). They are analogous to the universal "log-law" solutions for turbulent shear flows past a wall (e.g. Schlichting, 1979). By comparing profiles of r, f and q with laboratory or numerical experiments, one can constrain some of the unknown coefficients {Ci} (see Section 3.4).

Nusselt-Rayleigh number relationship
The heat flux through the system in Rayleigh-Bénard convection is commonly measured by the dimensionless Nusselt number which compares the total heat flux with the conductive one in the absence of convection. The universal convection-froma-wall solution calculated in the previous section can be used to derive the relationship between the Nusselt number and the Rayleigh number. Indeed, by selecting a Rayleigh number we set the relative temperature at the midpoint z = h/2 to beΘ = ∆T /2 which implies through (37) that yielding an equation for the (unknown) constant heat flux F0. In dimensionless terms, we have the implicit equation for Nu: which can be solved to find Nu(Ra). In the limit of very large Rayleigh number the mid-point of the system is very far from the boundary layer, so θ ≈ θ0 which then recovers the standard scaling law (Malkus 1954) The constant θ0 depends only on Pr and on the closure parameters {Ci}, but cannot easily be expressed analytically in terms of these parameters.

Comparison with data and estimation of the model parameters
The aim of this section is to estimate, in a rough sense, the parameters {Ci} by comparing the model predictions with numerical simulations and laboratory experiments. This approach was successfully used in GO05 on pipe flow data and Couette-Taylor data, yielding: Under the assumption that the closure parameters are universal properties of the turbulent cascade, these estimated values should also apply to the case of turbulent convection without need for re-calibration. The remaining parameters C6, C7, Cνκ and Cκ may then be independently estimated.
In the following sections, we first discuss this assumption in the light of known model limitations. We then select appropriate experimental datasets and use them to constrain the remaining parameters.

Discussion of the model limitations
As discussed by Ogilvie (2003) and GO05 the closure model proposed has two intrinsic limitations: it ignores some (but not all) of the effects of pressure-strain correlations < u ′ i ∂jψ ′ >, and assumes that the effect of all modelled terms (such as the triple-correlations in (28)-(30)) is local both in time and space. As a result, it may poorly represent strongly sheared systems or systems where the turbulent eddies exhibit a strong degree of spatial or temporal coherence.
The neglected effects of the pressure-strain correlations are not thought to be important in turbulent convection, except in the presence of strong rotation or of an externally driven strong mean shear (where the timescale of rotation and shear is comparable to that of the convection). The closure should be well-suited to model convection in stellar interiors, but maybe less so for convectively unstable accretion discs. We defer this particular case to subsequent work. However, for similar reasons these effects are also likely to be important in pipe flows or Couette-Taylor flow, which were used as a basis for calibrating the constants C1 and C2 (see GO05). Consequently, the estimates given in (46) could be somewhat biased, in particular C2 which contains information on the rate of return to isotropy. Comparing the model with turbulent convection experiments (see below) can therefore help refine the estimates for C1 and C2 using more appropriate data.
As mentioned above, the closure is also less reliable when applied to systems where the turbulence exhibits coherence over large scales or long timescales. This might pose some problems when applied to convection in a finite domain, since large-scale coherent plumes which span the whole system are commonly observed in most cases ranging from Boussinesq to fully compressible systems. Comparisons with experiments can help reveal which aspects of convective transport are adequately described by the model, and which are not.

Available experimental data
Our application of the closure model to Rayleigh-Bénard convection in Sections 3.2 and 3.3 assumes for simplicity that the system is horizontally invariant, while all laboratory and numerical experiments have a limited horizontal extent. The presence, nature and geometry of the side-walls are known to affect various properties of the turbulent convection, in particular through the generation of large-scale circulations (often called "wind"). This wind influences the overall heat transport properties by changing the nature of the boundary layers (Castaing et al. 1989;Cioni et al. 1997;Grossmann & Lohse 2000, 2001, 2002, 2004. It also induces large-scale horizontal inhomogeneities, so that the measured vertical profiles of mean quantities and higher-order moments may vary with position (Maystrenko, Resagk & Thess 2007). While our formalism can in principle be applied to finite geometries and self-consistently model the effect of large-scale flows, such an extension is beyond the scope of the present paper.
In order to minimize the effect of side-walls we restrict the model comparison to experimental setups with very large aspect ratios (defined as the ratio of the horizontal to vertical extent of the domain, and denoted as Γ). There are a few large aspect ratio, high Rayleigh number experimental studies which provide measurements of the Nusselt number. Of particular interest are results of Fünfschilling et al. (2005) for convection in water (Pr = 4.38) in a cylindrical enclosure of aspect ratio up to Γ = 6, for Ra up to a few times 10 10 . Niemela & Sreenivasan (2006) provide similar information for convection in Helium (0.7 < Pr < 8) in a cylindrical container with Γ = 4, for Rayleigh numbers between 10 8 and 10 13 . Finally, the Ilmenau barrel experiments of DuPuits, Resagk & Thess (2007) provide Nu(Ra) for convection in air (Pr = 0.7) in a cylindrical enclosure with variable aspect ratio up to 11.3, for Rayleigh numbers up to a few times 10 8 (in the case of the largest aspect ratio).
By contrast, only very few large aspect ratio experimental measurements of the boundary-layer profiles of velocity and temperature correlations (such asRij ,Fi orQ) have been reported. The largest aspect ratio experiments available (Γ = 11.3) with fully resolved boundary layer profiles are presented by DuPuits,  although the data provided is limited to the mean and rms temperature profiles.
Taking a different approach, direct numerical experiments are a powerful tool for "idealized" experiments. Horizontally periodic simulations minimize the effect of sidewalls (although retain a finite aspect ratio) and permit resolved and precise measurements of all desired mean and fluctuating quantities within the flow. The main drawback is the limited range of parameter space for which resolved simulations can be run (typically, Ra < 10 8 for large aspect ratio simulations at Pr = O(1)).
For these reasons, we use a combination of experimental data (DuPuits, Resagk & Thess 2007) and numerical simulations to calibrate the remaining model parameters. Our numerical simulations are all run for Pr = 1, in a horizontally periodic domain with aspect ratio Lx/Lz = Ly/Lz = 4, using a spectral method briefly described in Appendix B. The largest Rayleigh number achieved in this case is Ra = 2.1 × 10 7 . Figure 2 shows a typical snapshot of the results, in this parameter regime, for the temperature field for example. The results of the simulations are globally consistent with those of Hartlep (PhD thesis, 2005, Göttingen).

Near-wall profiles and estimation of Cν, Cκ and Cνκ
Very close to the wall (η ≪ 1), the closure model solutions for the normalized correlations r, rzz, f , q and θ are well approximated by power laws, as described in equation (39). These relationships can be compared with data and provide a simple way of individually estimating each of the model constants Cν, Cκ and Cνκ from laboratory or numerical experiments.
Comparisons of (39) with the experimental near-wall profile for rzz(η), f (η) and q(η) yield slopes αν close to 4 (see Fig. 3), ανκ close to 3 (see Fig. 4), and ακ close to 2 (see Fig. 5). Note that while the amplitude of the power-law observed in the near-wall profile for q(η) is seen to depend on the experiment considered, the slope ακ appears to be universal. We then adopt the following values for the constants Cν, Cνκ and Cκ: Given the experimental and model uncertainties, these values and their errorbars should be thought of as rough estimates rather than precise calibrations. It is comforting to note that this independent comparison recovers the value of Cν found by GO05. Moreover, we find that within fitting errors Cνκ ≃ (CνCκ) 1/2 . Given the quantities modelled by the associated diffusive terms (see equations (31)-(33)), this result is not entirely surprising.
On the other hand, Fig. 3 reveals an important caveat of the closure model when applied to Rayleigh-Bénard convection. The universal solution for the two horizontal stress components rxx(η) and ryy(η) can easily be deduced from rxx = ryy = 0.5(r − rzz). These horizontal stresses should therefore be identical to one another and have the same power-law dependence on η as r and rzz, close to the wall and far from the wall. However, Fig. 3 clearly shows that the numerical data is at odds with the model. We attribute the discrepancy to the presence of large-scale coherent convective plumes in the system, which span the entire domain and create strong horizontally correlated fluctuations as they crash against each boundaries. As a result, the fluid Figure 3. Comparison of the universal "convection from a wall" solution with numerical data for the dimensionless Reynolds stress components rzz, rxx and ryy. The large symbols represent rzz(η) for Ra = 2.1 × 10 6 (triangles) and 2.1 × 10 7 (diamonds). The two sets of smaller diamonds show rxx(η) and ryy(η) for the case where Ra= 2.1 × 10 7 . Note that theoretically these should be lying on the same curve -the difference can be attributed to limited statistics. In all cases Pr = 1. The dotted line shows the asymptotic solution rzz = r 0zz η 2/3 using the value of C 1 estimated by GO05, while the solid line shows a numerical integration of the full universal profile, for our estimated parameter values as listed in (46)  in the viscous sublayer is much more strongly anisotropic than predicted.
3.4.4 Far-field solution and estimation of C6 and C7.   (2007) for Ra = 8.14×10 8 for air (Pr = 0.7) in a cylindrical box at aspect ratio 11.3. The discrepancy between the numerical solutions and the experimental data is attributed to the difference between periodic side-walls and impermeable side-walls. The near-wall solution was used to fit Cκ while the far-from-the-wall data was used to provide a constraint between C 6 and C 7 . The solid line show a numerical integration of the full universal profile as in Figs. 3 and 5 for Pr = 1. rzz0η 2/3 . Note that rzz0 depends only on two numbers, the Prandtl number (which is known) and the model parameter C1. It is reassuring to see that the value of C1 estimated by GO05 from wall-bounded shear flow data adequately fits the far-from wall solution for rzz in this convection problem.
The universal profiles away from the wall listed in equation (40) can also be used in conjunction with numerical and laboratory experiments to constrain C6 and C7. These constants are unfortunately difficult to extract directly from our numerical simulations. The highest Rayleigh number available (Ra = 2.1 ×10 7 ) only has a short asymptotic (η ≫ 1) range, so that estimates of C6 and C7 from these datasets are unreliable 1 . The rms temperature data measured in various laboratory experiments at higher Rayleigh number provides a more adequate point of comparison. We use the rms temperature data of the highest aspect ratio experiments of DuPuits, Resagk & Thess (2007), for Ra = 8.14 × 10 8 . This dataset exhibits a significant asymptotic range, with a power law close to the one predicted by the closure model (q ∼ q0η −2/3 ). Fitting the data yields q0 = 0.95 ± 0.05, which provides a first constraint between C6 and C7 (see Fig. 6). Note that other datasets (from Maystrenko, Resagk & Thess, 2007, for example) are generally consistent with this estimate for q0.
A second constraint between C6 and C7 is obtained by comparing the model predictions with experimental measurements of Nu(Ra). The closure model implies that Nu = 1 + KRa 1/3 where the constant K is a function of the  (2007) are reasonably well approximated by taking K = 0.06 ± 0.003. Variations of K with Prandtl number, for the range of experiments discussed, are within the errorbars. Given that C1, C2, Cν, Cκ and Cνκ are now known, for fixed Prandtl number, fitting K provides a unique relationship between C6 and C7, as seen in Fig. 6.
By combining these two constraints, we conclude that a good fit to the data can be obtained with The values for {Ci} quoted in equations (46), (47), and (48) form from here on our selected set of parameters. These values are to be taken as indicative estimates, rather than precise calibrations. We note that the parameters derived do satisfy realizability (see Appendix A). The solid lines shown in Figs. 3, 4 and 5 are the universal boundary layer profiles calculated using these parameters, and are seen to fit all datasets (except for rxx and ryy, as discussed above) satisfactorily. Fig. 7 compares our closure model prediction for the Nu(Ra) relationship, using the estimated parameters, with various available datasets for large aspect ratio experiments (Γ 4). It also shows (as dashed lines), for comparison, strict upper bounds obtained by Plasting & Kerswell (2003) and by Ierley, Kerswell & Plasting (2006) for transport by convection at finite and infinite Prandtl numbers respectively. It is reassuring to see that the Pr → ∞ prediction from our own closure model remains below the strict upper bound for the same limit. In conclusion, our model successfully reproduces most measurable features pertaining to laboratory and numerical experiments of Rayleigh-Bénard convection, for reasonable values of the model parameters {Ci}. Furthermore, comparison of the estimated parameter values across a range of experiments in other systems (pipe flows, Couette-Taylor flows) shows that they are indeed of a universal nature, a results which can only increase confidence in our approach.

Introduction
Another system that is of interest, and possibly more relevant to astrophysical applications, consists of an unbounded layer in which there is no mean flow, while the mean temperature gradient ∇Θ is uniform and parallel to the gravitational acceleration (taken to be in the z-direction). The evolution of perturbations to this mean state can be described by the following set of Boussinesq equations: where all perturbations are triply periodic, as for example This model setup is now commonly referred to as Homogeneous Rayleigh-Bénard (HRB) convection (Borue & Orszag 1997;Lohse & Toschi 2003;Calzavarini et al. 2005;Calzavarini et al. 2006). While this system cannot be studied using laboratory experiments, it lends itself relatively easily to numerical experimentation using spectral methods in particular. The relevant dimensionless parameters are the Prandtl number Pr = ν/κ, the Rayleigh number, now defined as and the aspect ratio(s) Γ = Lx,y/Lz. The microscopic diffusivities are included in the original equations (49) to regularize the system by allowing for dissipation and irreversibility. However, note that the periodic boundary conditions forbid the formation of boundary layers, so it may be conjectured that the macroscopic statistical properties of the turbulent convection should be well defined and independent of ν and κ in the limits Ra → ∞ (Spiegel 1971). Furthermore, we may expect the turbulence to be statistically steady and homogeneous, although anisotropic. These properties have been argued to be more relevant to convection in astrophysical systems than standard Rayleigh-Bénard convection. The HRB model may therefore provide a suitable local model of convection deep inside a star or planet.
On dimensional grounds, the rms turbulent velocity, for example, must be expressible in the form where f is a dimensionless function. According to the discussion above, f should tend to a non-zero function of Γ alone in the limit Ra → ∞. It is tempting to conjecture that f also becomes independent of Γ in the limit of large aspect ratio, Γ → ∞. This would imply that the vertical length-scale Lz plays a fundamental role in determining the saturation level of the turbulent convection, presumably by limiting the size of coherent structures ('eddies'). For convection deep inside a star or planet, it is the pressure scale-height that imposes a characteristic vertical scale on the turbulence (see Section 5); in the local model, the vertical extent of the box plays an equivalent role. In practice, owing to some peculiarities of the HRB system discussed below, the role of the aspect ratio in the behaviour of the solutions is not so straightforward.

Governing equations
Applying our closure model to HRB, and noting that all statistical averages are now independent of position, we obtain the system of ODEs for the temporal evolution of the second-order correlationsRij ,Fi andQ: where we have ignored for simplicity contributions from terms including Cν, Cκ and Cνκ which do not contribute to the high-Rayleigh number dynamics of HRB convection. Note that the resulting equation forR is so that these equations consist of a main system for (R,Rzz,Fz,Q), decoupled systems for (Rxz,Fx) and (Ryz,Fy), and prognostic equations forRxx,Ryy andRxy.

Choice of L and consequences for the coefficients {Ci}
While selecting L as the distance to the wall is a natural choice for wall-bounded convection or shear flows, a different approach must be used for triply periodic flows. The largest eddy size in this case is limited by the horizontal and vertical scales in the box, so that L can be assumed to be proportional to min(Lx, Ly, Lz).
It is important to note that the selection of a different L implies a potential rescaling of the {Ci} coefficients. For example, had we selected L = z/2 in the wall-bounded case instead of L = z, then the estimated C1, C2, C6 and C7 would all be half the values quoted in Section 3.4 since these parameters enter the model in the combinations C1/L, etc. Nevertheless, the ratios of any pairs of constants within the group {C1, C2, C6, C7} should (presumably) be preserved. Following these considerations, we elect to keep the estimated values of the {Ci} given in equations (46) and (48), and calibrate instead the value of the proportionality constant δ in the expression L = δ min(Lx, Ly, Lz).

High Rayleigh number HRB convection
A search for non-trivial fixed points of the dynamical system (53) (withR > 0) reveals they are the (positive) solutions of a quartic equation. In the limit of large Ra it can be shown that there is only one positive fixed point with Rxy =Rxz =Ryz = 0, Note that, in this case, This solution represents a state of fully developed turbulent convection, which is statistically steady and homogeneous. The solution exists in the statistically axisymmetric subspace in whichRxx =Ryy andRxy =Rxz =Ryz =Fx = Fy = 0, and is stable with respect to perturbations transverse to this subspace. It has the desired properties that the vertical motion is dominant (Rzz >Rxx =Ryy), while the heat flux is purely vertical and directed down the temperature gradient. Moreover, numerical integrations suggest that, where it exists, this state is stable and universally attracting.
Defining the Nusselt number Nu as the ratio of the total to the conducted heat flux, we have, in the limit Ra ≫ Pr, This scaling recovers the "ultimate turbulence" regime, where the turbulent transport properties are independent of microscopic diffusivities (Spiegel 1971

Comparison with numerical experiments
Numerical simulations of HRB convection were first performed by Borue & Orzag (1997). More recently, Toschi & Lohse (2003) and Calzavarini et al. (2005) performed a range of Lattice-Boltzmann simulations in a cubic geometry, for various values of the Rayleigh and Prandtl numbers, and report on the first evidence for scalings consistent with the "ultimate regime" of convection, namely Nu ∝ (RaPr) 1/2 and Re ∝ (Ra/Pr) 1/2 . However, it is now recognized that the dynamics of HRB convection are more subtle than previously thought. As discussed by Calzavarini et al. (2006), simulations at unit aspect ratio show huge fluctuations in the instantaneous Nusselt and Reynolds numbers arising from the intermittent or quasi-periodic (depending on Ra) exponential growth of socalled "elevator modes". These modes are thus named because they are independent of z, and have the peculiar property of being exact nonlinear and exponentially growing solutions of the governing equations (49). The most unstable mode has a horizontal wavelength equal to the larger horizontal dimension of the box. Hence, the aspect ratio of the system directly influences the macroscopic solution.
This phenomenon has a close parallel in shearingbox studies of the magnetorotational instability. In that case, forcing by a constant velocity gradient plays the role of the constant temperature gradient, while perturbations to the background fields are also assumed to be triply periodic. This system is unstable to equivalent "channel modes", exact nonlinear and exponentially growing solutions of the equations and associated periodic boundary conditions (Goodman & Xu, 1994). In this case, it is known that the channel modes are themselves subject to secondary shearing instabilities which limit their growths. However, the existence and growth rates of shearing instabilities depend sensitively on aspect ratio: they are strongly inhibited in systems where the streamwise direction is smaller than the cross-stream directions. As a result, systems with roughly cubic geometry are dominated by the channel modes and are found to have very strongly fluctuating large-scale transport properties, but for larger aspect ratio the fluctuations are much smaller and the channel modes are inhibited (Bodo et al. 2008).
For these reasons, we performed a series of HRB simulations of various aspect ratios, in order to determine whether the same phenomenon occurs, and to provide a better point of comparison for the closure model. Appendix C provides a brief description of the numerical algorithm used, and the results are summarized in Fig. 8. We studied 5 cases, with Lx = Ly and Lx/Lz =1/2, 2/3, 9/10, 1/1 and 4/3. In the last case, the elevator modes continue growing unaffected by perturbations until the code fails, which seems to corroborate the premise that the secondary instabilities are inhibited in wider-than-tall boxes. For Γ < 1, the measured Nusselt number eventually converges to a meaningful average and is found to scale as predicted by the closure model, namely proportional to (Pr Ra) 1/2 Γ 2 . A good fit with the model predictions is found by selecting L = δLx = Lx/ √ π.
For the purpose of illustration, a snapshot of the temperature field for our largest Rayleigh number, Ra = 5 × 10 6 (with Pr = 1) and aspect ratio 1/2 is shown in Figure 9.

The effect of rotation on homogeneous turbulent convection
We now consider the effect of rotation on HRB convection, where the rotation axis lies at an angle γ from the vertical direction: Ω = (0, Ω sin γ, Ω cos γ). In this section it is more  . Volume-rendered visualization of the temperature field for Ra = 5 × 10 6 and Pr = 1 for homogeneous convection in a box of aspect ratio 1/2. Note how, even at this high Rayleigh number, the size of the dominant structures is equal to the box size.
convenient to work with dimensionless variables so we select the following scalings: where for convenienceÑ is defined asÑ 2 = −N 2 , and is positive when the fluid is convectively unstable. The convective Rossby number is then defined as Stationary solutions of the closure model far from onset of convection satisfy the following equations: In the infinite Rossby number limit (equivalently in the non-rotating limit), the solution of these equations reduces to the non-dimensional form of (55) and (56). Should all of the quantities be expanded in terms of the inverse Rossby number as (for example) then we find that (and similarly for all diagonal components ofR). Our expressions for the non-diagonal terms, to first order, recover the equivalent of the well-known Λ-effect (see Rüdiger, 1989) in the coefficientRxz: The Λ-effect, as seen in the above equation, describes how rotationally constrained turbulent motions can drive differential rotation, through the non-diagonal component of the stress-tensorRxz. As expected from dimensional analysis and geometrical arguments, its amplitude scales linearly with sin γ Ω. The other two componentsRxy andRyz only become important for more rapidly rotating systems as they are both O(Ro −3 ). Finally, a non-negligible horizontal heat flux is generated in the direction of Ω × g, namelŷ although note that when applied to stellar convection zones,  (46) and (48). The Ro −1 → 0 and Ro −1 → ∞ asymptotes satisfy equations (67) and (70) respectively.
this effect is relevant only for non-axisymmetric heat transport. The "latitudinal" heat fluxFy on the other hand is of higher order in Ro −1 .
In the opposite limit of very low Rossby number (the rapidly rotating limit) an expansion in powers of Ro reveals thatR = 2 cos 2 γ C1C6 so that the rms velocity is reduced by a factor cos γ compared with the non-rotating case. Note, however, how the expected reduction (and potential suppression) of the convective heat flux in rapidly rotating systems where gravity is aligned with the rotation axis (Chandrasekhar, 1961) so that γ = 0 is not captured by this closure model. This problem, which was identified by Miller & Garaud (2007), can presumably be attributed to the incomplete modeling of the effects of the pressure-strain correlations which are known to play an important role in the limit of rapid rotation. It is therefore likely that these effects also cause our model to yield inaccurate predictions for γ = 0 in the same limit. A full resolution of the issue must eventually involve the derivation of a better closure for the pressure-strain correlation terms. For completeness note that in this limit the model predicts that a significant heat flux is carried horizontally along ey, with amplitudeFy = tan γFz, and that whileRxy andRxz are both O(Ro). Fig. 10 shows the variation of the normalizedR as a function of both γ and Ro −1 , while Fig. 11 shows the variation of the normalized −Rxz/R as a function of both γ and Ro −1 , illustrating the dependence of the Λ-effect on both parameters as predicted by our model.  (46) and (48). Note thatRxz/R ∝ Ω for low rotation rates, and to Ω −1 for large rotation rates.

Comparison with previous second-order models
We now compare our findings with the commonly used model for convective stresses originally proposed by Rüdiger & Kitchatinov (1993) and later extended by Rüdiger et al. (2005, Ral05 hereafter). Note that the related theory of  relies on the presence of a background density stratification to explain the Λ-effect. As such it is not an appropriate point of comparison for our Boussinesq calculation. Rüdiger & Kitchatinov (1993) and Ral05 assume the presence of a "background" turbulence caused by a given (unspecified) forcing mechanism, which, in the absence of rotation, is described by an eddy turnover time τ , a mixing length l and a turbulent diffusivity νt = l 2 /τ . This background turbulence also may also have some degree of anisotropy, controlled by the parameter a defined in our notation as where the superscript (0) denotes turbulent quantities of the non-rotating system. Note how a = 0 for isotropic turbulence. Ral05 show how the presence of rotation (where the rotation axis lies at an angle γ from the vertical) modifies the background turbulence, an effect which gives rise to nondiagonal components in the stress tensor. They argue that this phenomenon is controlled by the Coriolis number Ω * defined as Their eddy turnover time τ is naturally related to L/ √ R (0) in our notation, so that, for the purpose of comparison we have where the proportionality constant is of order unity.
In the limit of slow rotation, Ral05 predict a Λ−effect through the following term: where the proportionality constant is the same as in equation (74). Meanwhile, our model when written in dimensional form and recast in terms of the anisotropy factor a yields zz .
(76) whereR (0) is a dimensionless constant which depends only on the model parameters, and is given by equation (67) with Ro −1 = 0. In the same slow-rotation limit, the other offdiagonal components of the stress tensor are O(Ro −2 ) or higher order in both our and their models.
Overall, the two formalisms agree on the dependence of the stresses on the rotation rate and on latitude, as expected on dimensional and geometrical grounds. In addition, both models explicitly demonstrate the importance of the anisotropy of the background non-rotating turbulence in controlling the amplitude of the Λ-effect. However, the dependence ofRxz on the anisotropy factor a superficially appears to be different in the two theories. We interpret this in two ways. First, note that the anisotropy factor a is a "free" parameter in the works of Ral05. In our model by contrast, there is no freedom in independently specifying a since it is a solution of the model once the system is specified (e.g. shearing flow, convective flow) and depends on the {Ci} parameters. In the HRB system for example a = −6C1/(3C1 + C2).
Secondly,Rxz is directly proportional to a in the model of Ral05 while our model reveals an additional contribution to the Λ-effect arising from the background turbulent heat flux (see equation (68) for a more explicit expression). This contribution is missing from the model of Ral05 which does not take into account the heat equation. As a result, one may superficially conclude that the Λ-effect could exist even for isotropic background turbulent convection. In practice, it is difficult to conceive of a naturally occurring isotropic turbulent system which has a non-zero vertical heat flux, so the termF (0) z is in fact also indirectly related to the anisotropy of the system, although perhaps not exactly in the same way.
Finally, we emphasize that in the limit of rapid rotation, neither theory is expected to be accurate because of the extreme induced anisotropy of the rotating turbulent motions. Nevertheless it is interesting to note that the predicted dependence of the stresses on the rotation rates now no longer agree with one another. We find thatRyz tends to a constant independent of rotation rate while Ral05 find thatRyz ∝ Ro. For the other off-diagonal componentsRxz andRxy we find a dependence on Ro, while they predict a dependence on Ro 2 .
We conclude this section by emphasizing the success of our closure model in reproducing numerical experiments of HRB convection at various aspect ratios and Rayleigh numbers. Furthermore our model predictions are exactly proportional to those of Ral05 (with a proportionality constant of order unity) for convection in a slowly rotating system. Hence we expect to recover many of the results and suc-cesses of these authors in modeling differential rotation in stars, albeit with an extended model which self-consistently includes heat transport in addition to angular momentum transport. In preparation of this future modeling endeavour, we finally turn to the next natural step of this work, namely the extension of the model to the anelastic and fully compressible equations.

THE ANELASTIC SYSTEM AND COMPRESSIBLE FLOWS
So far we have worked within the Boussinesq approximation, which is applicable only to a shallow layer of fluid whose depth is much less than the density scaleheight. In order to apply our model to stars we must first adapt it to the anelastic approximation (Ogura & Phillips 1962;Gough 1969), which is relevant to subsonic convection in a deep layer.
Here we follow the more standard derivation of the anelastic approximation where the reference state is taken to be an adiabatically stratified fluid in hydrostatic equilibrium. The reference density ρ0(r) and temperature T0(r) may vary substantially, while the specific entropy s0 is uniform. In place of equations (9)-(11) we have where the dots represent terms due to viscosity (in the equation of motion) and thermal conduction (in the thermal energy equation), while ψ is, again, a modified pressure. Viscous dissipation can also be included in the thermal energy equation, although it is usually omitted in the Boussinesq approximation. A derivation of these equations, omitting diffusive effects, is given in Appendix B.
The anelastic system is formally very similar to the Boussinesq system except for the variable density of the reference state. However, the entropy perturbation and background temperature gradient respectively play the roles taken by the temperature perturbation and αgi in the Boussinesq approximation. A very similar analysis to that carried out for the Boussinesq system leads to equations for Rij ,Fi andQ of the form (∂t +ūj ∂j)Fi +Rij ∂js +Fj∂jūi +Fi∂jūj +Q∂iT0 where the dots represent terms that require a closure model. In the anelastic system the relevant definitions of the Reynolds stressRij, fluxFi and varianceQ arē (83) Note thatRij now has the correct dimensions for a stress tensor, and thatFi is really an entropy flux density. Some additional linear terms arise in the anelastic system because ∂iūi = 0.
We apply the same closure model as for the Boussinesq system, except that the relaxation timescale which was proportional to L/R 1/2 is now proportional to L/(R/ρ0) 1/2 because of the redefinition ofRij: (∂t +ūj ∂j)Fi +Rij ∂js +Fj∂jūi +Fi∂jūj +Q∂iT0 We do not include any of the terms proportional to ν or κ here because we consider the high-Rayleigh number limit in the absence of rigid boundaries only. The question arises as to how the length-scale L should be identified for anelastic convection in a deep layer. It should probably related to the pressure scaleheight or density scaleheight, as in the stellar mixing-length theory. Indeed, numerical simulations of convection in spherical shells with a substantial density variation indicate that the convective cells are much smaller near the outer surface where the scaleheight is small; nevertheless, there may be situations in which convective plumes can span several scaleheights.
Equations (84)-(86) can then be combined with equations for the mean variables in the form In the last equation we have included the turbulent viscous heating. These equations could be applied to studying convection and meridional circulation in rotating stars. The solution can be assumed to be axisymmetric and independent of time, although for practical purposes it may be easier to evolve the equations forwards in time until a steady state is reached (if it is) rather than directly seeking such a solution.
In the absence of rotation the problem becomes spherically symmetric, the mean flow disappears, the stress becomes diagonal (although anisotropic) and we obtain the local algebraic system The solution is, by direct analogy with equations (55)- (56), where now In this situation the entropy gradients is not known in advance. However, to balance the thermal energy equation, ∂iFi = 0, which implies that r 2F r is a constant, determined by the luminosity generated by the core of the star. (This conclusion is modified if the radiative flux or any sources of energy such as turbulent viscous dissipation make an important contribution to the thermal energy equation.) Then the above equations can be solved algebraically to find ∂rs, R, etc., at each radius, assuming that a prescription for L is given. The result is equivalent to a version of mixing-length theory.
Rotation couples radial and latitudinal transport of heat and momentum and induces large-scale entropy gradients and mean flows. However, if we assume that their effects can be ignored in the overall turbulent dynamics controlling the properties of the stresses, then the local Λ−effect is easily recovered as an anelastic version of equation (68). As before, the only differences with the Boussinesq case is that (i) the two terms containingR (0) , which have their origin in the eddy turnover time, should be replaced byR (0) /ρ0 and (ii) in expressing (68) in dimensional form (see equation (61)), one must also replaceÑ 2 by (∂rT0)∂rs and dΘ/dz by ρ0∂rs, as seen above. The resulting expression then directly links the turbulent transport of angular momentum and of heat to one another. Since heat transport in this model is very similar to mixing-length theory, our formalism now provides a simple framework in which to combine models of stellar structure with models of internal stellar dynamics. Note that in practice mean flows and especially latitudinal entropy gradients could play a role in the global dynamics of the system. The whole model should therefore be solved self-consistently and globally instead of using (68). This can only be done numerically and is deferred to a subsequent paper.
It is also possible to 'import' the model of anelastic convection into the full set of equations governing the motion of a compressible fluid. The idea here is that, while the convection might be assumed to be subsonic and to obey the anelastic approximation, the mean flow need not obey these constraints. An example is convection in an accretion disc, where the accretion flow, although slow, cannot be treated in the anelastic approximation with a reference density profile. Omitting now the bars on all quantities, and neglecting self-gravitation (although it can easily be restored), we propose a system of equations consisting of the equation of mass conservation, the equation of motion, and the thermal energy equation, together with the equations of the closure model, (∂t + uj ∂j)Fi + Rij ∂js + Fj∂jui + Fi∂juj + Q∂iT The total energy is then exactly conserved in the form ∂t ρ( 1 2 u 2 + Φ + e) + 1 2 R +∂i ρ( 1 2 u 2 + Φ + h)ui + 1 2 Rui + Rij uj + T Fi = 0, where e and h are the specific internal energy and the specific enthalpy, respectively, and the gravitational potential Φ is assumed to be independent of time. The existence of this conservation law implies a certain self-consistency in the equations of the model. The terms that were added in passing to the compressible model are required to have the form that they do in order that energy be conserved. We note again that Fi is really an entropy flux density, and that T Fi is the corresponding energy flux density. The physical content of this model is that the turbulent convecting fluid behaves similarly to a complex, non-Newtonian material in which there is a dynamical constitutive equation that relates the stress tensor to the deformation history of the fluid. The above equation for ∂tRij (along with those for ∂tFi and ∂tQ) plays this role.

CONCLUSIONS AND FUTURE PROSPECTS
We have laid out the foundations of a new second-order closure model for the dynamics of turbulent convection, with future applications to stellar convective regions in mind. This model is a direct extension of the work of Ogilvie (2003) and GO05, and has similar properties. The proposed closure has a straightforward physical interpretation, and wellunderstood limitations.
Comparison with laboratory and numerical experiments reveals good overall agreement of the model predictions with known properties of rotating shear flows (GO05) and high Rayleigh-number rotating convection (this work). In particular, our model naturally reproduces the standard scaling relationships between the Rayleigh and Nusselt numbers for Rayleigh-Bénard convection and for Homogeneous Rayleigh-Bénard convection, and contains the well-known Λ-effect describing angular momentum transport in a rotating turbulent fluid.
When extended to the anelastic (or fully compressible) case, our formalism can be applied to study convection in stellar interiors. Note that the effects of Maxwell stresses can also straightforwardly be included following Ogilvie (2003) if needed. We show that the model naturally reduces to a version of mixing-length theory when applied in a onedimensional framework. In the presence of rotation it becomes a powerful tool to study within a single framework the multi-dimensional balance involving large-scale mean quantities such as the entropy profile, the meridional circulation and the differential rotation. Future work applying our closure model in a spherical shell geometry will help understand some of the trends seen in the increasingly large number of available observations of stellar differential rotation. non-negative. We work here with the anelastic version of the closure model, equations (84)-(86), although a similar argument applies to the Boussinesq system in the high-Rayleigh number limit in which the terms proportional to ν or κ are omitted.
Let Ai be any vector with appropriate dimensions, and consider the tensor (within the anelastic system) at any point in the flow. The associated quadratic form is Evidently S 0 for all vectors Xi, and therefore Sij must be a positive semi-definite tensor, for any choice of Ai, at every point in the flow. Allowing the vector Ai to vary provides us with a family of realizability conditions. On the other hand, Sij can be expanded as Sij =Rij +FiAj +FjAi +QAiAj , and therefore S =Rij XiXj + 2(F · X)(A · X) +Q(A · X) 2 = (Rij −Q −1F iFj )XiXj +Q −1 (F +QA) · X 2 . (A5) Provided thatQ > 0, a necessary and sufficient condition for S to be non-negative for all choices of Ai and Xi is that the tensor be positive semi-definite. If this condition is satisfied then TijXiXj 0 for all Xi and therefore S 0 for all Xi and Ai. On the other hand, if a vector Xi exists such that TijXiXj < 0, then S < 0 for this choice of Xi if we set Ai = −Fi/Q.
We therefore aim to show that, provided that the condition (A1) is satisfied, the model ensures that Tij remains positive semi-definite at all points in the flow if, in the initial state, Tij is positive semi-definite andQ > 0 at all points.
We apply a reduction ad absurdum. If Tij has a negative eigenvalue at any event, then the quadratic form is negative at that event, for some choice of the vector Xi. Without loss of generality, let Xi be a differentiable vector field advected according to the time-reversible equation and such that T < 0 at the event in question. Here D = ∂t + ui∂i is the Lagrangian derivative following the mean flow.
Retracing the the value of T to the initial state, following the mean flow in reverse, we deduce that T must have passed through zero with DT < 0. However, using equations (84)-(86) we find DT = −T ∂iūi + 2Q −1 (F · X)XiTij ∂js −(C1 + C2)L −1 R ρ0 When T = 0, Xi is a null eigenvector of Tij and therefore DT 0, with equality only when there is no turbulence. Therefore Tij cannot in fact develop a negative eigenvalue.
horizontal directions, whereas no de-aliasing procedure is employed along the vertical coordinate. We use the Chebyshev tau method (Peyret, 2002) to solve the ODEs arising from the implicit part of the time-stepping scheme. This method has the advantage of yielding linear systems which can be manipulated into sparse form, thus keeping the memory requirements at a manageable level. The code is parallelized using transpose-based parallel FFTs, see Stellmach & Hansen (2008) for details.

APPENDIX C: NUMERICAL ALGORITHM FOR HOMOGENEOUS RAYLEIGH-BÉNARD CONVECTION
A spectral algorithm using Fourier expansions in all three spatial directions is used to solve the governing equations (49) in the homogeneous Rayleigh-Bénard case. The primitive variables u ′ ,T ′ ,ψ ′ are used, with the pressure perturbation ψ ′ being calculated in the same way as in the classical Patterson-Orzag Algorithm (Canuto et al. 2007) which is widely used in simulations of homogeneous, isotropic turbulence. Non-linear products are evaluated on a grid in physical space and aliasing errors are avoided by using the 3/2-rule. The equations are advanced in time by a semiimplicit multistep method in which all diffusive terms are treated implicitly by a third order Backward-Differencing (BDF3) algorithm, while a third-order Adams-Bashforth (AB3) scheme is applied to the non-linear terms. This time stepping method offers a relatively large stability domain that includes a part of the imaginary axis at a comparatively low cost (Peyret, 2002). As a starting scheme, we use a second-order Runge-Kutta method. A parallelization approach similar to the one employed in our Rayleigh-Bénard convection code is used, see Stellmach & Hansen (2008) for details.

APPENDIX D: DERIVATION OF THE ANELASTIC SYSTEM
The equations governing the motion of an ideal, compressible fluid are the equation of mass conservation, ∂tρ + ∂i(ρui) = 0, We adopt a system of units in which the pressure scaleheight and the sound speed are of order unity. We introduce a small parameter ǫ such that the (imaginary) Brunt-Väisälä frequency is O(ǫ) when expressed in these units. We where τ = ǫt is a slow time variable. Note that the reference state is adiabatically stratified, and therefore s0 is independent of r. The equation of motion at leading order [O(1)] requires hydrostatic equilibrium in the reference state, Poisson's equation at O(ǫ 2 ) is not required, and the departures from the reference state are not affected by selfgravitation at leading order. When the asymptotic scalings are removed, and allowance is made for diffusive effects, equations (77)-(79) are obtained.