Phenomenology of buoyancy-driven turbulence: recent results

In this paper, we review the recent developments in the field of buoyancy-driven turbulence. Scaling and numerical arguments show that the stably-stratified turbulence with moderate stratification has kinetic energy spectrum $E_u(k) \sim k^{-11/5}$ and the kinetic energy flux $\Pi_u(k) \sim k^{-4/5}$, which is called Bolgiano-Obukhov scaling. The energy flux for the Rayleigh-B\'{e}nard convection (RBC) however is approximately constant in the inertial range that results in Kolmorogorv's spectrum ($E_u(k) \sim k^{-5/3}$) for the kinetic energy. The phenomenology of RBC should apply to other flows where the buoyancy feeds the kinetic energy, e.g. bubbly turbulence and fully-developed Rayleigh Taylor instability. This paper also covers several models that predict the Reynolds and Nusselt numbers of RBC. Recent works show that the viscous dissipation rate of RBC scales as $\sim \mathrm{Ra}^{1.3}$, where $\mathrm{Ra}$ is the Rayleigh number.


Introduction
Gravity pervades the whole universe, and it plays a dominant role in the flow dynamics of the interiors and atmospheres of planets and stars. The gravitational force also affects stratification. Under equilibrium condition, the density profile is where ρ b , ρ t are the densities at the bottom and top layers respectively (see Fig. 1).
We denoteρ(z) as the mean density profile. With fluctuations, the local density ρ l (subscript l stands for local) is ρ l (x, y, z) =ρ + ρ(x, y, z).
The gravitational force on a unit volume is −ρ l gẑ, where −gẑ is the acceleration due to gravity. Hence the gravitational force density on the fluid is The force ρgẑ occurring due to the change in density from the local value is the buoyancy. It is along −ẑ for ρ > 0 , but alongẑ for ρ < 0. The fluid flow is described by the Navier-Stokes equation where u, p are the velocity and pressure fields respectively, µ is the dynamic viscosity of the fluid, and f u is the external force in addition to the buoyancy. Substitution of Eq. (4) in Eq. (5) yields ρ l ∂u ∂t + (u · ∇)u = −∇σ − ρgẑ + µ∇ 2 u, where σ = p + g zρ (z )dz (7) is the modified pressure. The continuity equation for the density is ∂ρ l ∂t + ∇ · (ρ l u) = ∇ · (κ∇ρ l ), where κ is the diffusivity of the density. We assume that κ is constant in space and time. We can rewrite Eq. (8) as Now we employ Oberbeck-Boussinesq (OB) approximation according to which (dρ l /dt)/ρ l ≈ 0. Hence the relative magnitude of ∇ · u is where L, U are the large length and velocity scales respectively, and Pe is the Péclet number. Hence for large Pe, which is often the case for buoyancy-driven flows, we can assume that ∇ · u = 0. Therefore, under the Oberbeck-Boussinesq approximation, Eq. (8) gets simplified. In addition, we replace ρ l of Eq. (6) with the mean density of the fluid, ρ m . Hence the governing equations for the buoyancy-driven flows are ∂u ∂t where ν = µ/ρ m is the kinematic viscosity. The assumption that ν, κ are constants in space and time is also considered to be a part of the OB approximation. Also note that the buoyancy term, which is a function of variable density, is retained in the Navier-Stokes equation since it is comparable to the other terms of the momentum equation. In the stably-stratified turbulence, the total energy decays without f u , hence, f u is employed to maintain a steady state. Note that the system is stable when heavy fluid is below the lighter fluid, or dρ/dz < 0. Such systems yield wave solution in the linear limit. On the contrary, when heavy fluid is above the lighter fluid, dρ/dz > 0 and the flow becomes unstable and convective. See Schematic diagram of Fig. 1 for an illustration.
Temperature field T induces density variation in the following manner: where α is the thermal expansion coefficient, which is assumed to be constant in space and time. Hence we can rewrite Eqs. (11,12) in terms of the temperature field. Let us consider a fluid confined between two thermally-conducting horizontal plates kept at constant temperatures, as shown in Fig. 1(b). We denote the temperatures of the bottom and top plates to be T b and T t respectively, and ∆ = T b − T t . Thermal convection is absent for small ∆. Under this condition, the temperature profile is linear as and the heat is transported by conduction. This configuration has no fluctuation, i.e., u = 0 and ρ = 0. The flow however becomes unstable and convective when ∆ exceeds a certain critical value. For such flows it is customary to write the temperature as where θ is the temperature fluctuation over the background conduction profileT . A comparison of Eqs. (3,13,15) yields substitution of which in Eqs. (11,12) yields the following set of governing equations: The above fluid configuration under OB approximation is called Rayleigh-Bénard convection (RBC). For moderate temperature difference, say 30C for water, Onerbeck-Boussinesq approximation is satisfied. The flow dynamics of RBC is described by Eqs. (17,18,19).

Non-Boussinesq flows
Oberbeck-Boussinesq approximation provides a useful simplification for the analysis of the fluid flow. Without this approximation, we would need to solve the equations for the velocity, density, and temperature fields. For an illustration, refer to the set of equations in Sameen et al. [81]. The above description, called non-Boussinsq convection, is useful in stellar convection where the temperature difference is too large for the OB approximation to be valid. This topic, however, is beyond the scope of this paper.

Nondimensionalized equations
Fluid flows are conveniently described by nondimensional equations since they capture relative strengths of various terms of the equations. Also, they help reduce the number of parameters of the system, which is quite useful for analysis, as well as for the numerical simulations and experiments. Equations (11,12) have been nondimensionalized in various ways. Here, we present two such schemes. When we use d as the length scale, κ/d as the velocity scale, d 2 /κ as the time scale, and ∆ρ = |ρ b − ρ t | as the density scale, we obtain the following nondimensional equations: where ρ → ρ/(∆ρ), and For the stably-stratified flows, S = −1, but S = 1 for RBC. Using Eqs. (16) we can write the above equation in terms of temperature field as follows: for which where ∆ is the temperature difference between the bottom and top plates, as defined earlier. Note however that for large Ra, the aforementioned nondimensional velocity becomes very large (∼ √ RaPr) [34,102] that becomes an obstacle for numerical simulations due to very small time-steps. Hence, in numerical simulations, it is customary to employ √ αg∆d as the velocity scale, which yields the following set of equations: For stably-stratified flows, researchers often employ dimensional equations, but with density converted to units of velocity by a transformation [54] where is the Brunt-Väisälä frequency. In terms of the above variables, the equations become ∂u ∂t ∂b ∂t The other important nondimensional parameters used for describing the buoyancydriven flows are Froude number Fr = u rms dN , Richardson number Ri = 1 where u rms is the rms velocity of flow. Note that the Richardson number is the ratio of the buoyancy and the nonlinearity (u · ∇)u. Another important nondimensional parameter for RBC is the Nusselt number Nu, which is the ratio of the total heat flux (convective plus conductive) and the conductive heat flux, and is computed using the following formula:

Boundary conditions
For the velocity field we employ the following set of boundary conditions: (i) No-slip: All the components of the velocity field vanish at the walls, i.e., u = 0.
(ii) Free-slip: At a wall, the normal component of the velocity field vanishes, i.e., u ·n = 0, and the gradient of the parallel components of the velocity vanishes, i.e., ∂u /∂n = 0.
(iii) Periodic: The velocity is periodic, i.e., u(x + lL xx + mL yŷ + nL zẑ ) = u(x), where l, m, n are integers, and the box is of the size L x × L y × L z .
For the temperature field, the typical boundary condition used are (i) Conducting : Uniform temperature field at the walls, i.e., θ = 0.

Exact relations
Equations (11,12) are nonlinear, and hence researchers have not been able to write down general analytic solutions for them. However, Shraiman and Siggia [85] derived the following exact relations for RBC flows: Also, in the idealized limit of ν = κ = 0, using Eqs. (32,33), we deduce that the total energy is conserved for periodic and vanishing boundary conditions. In the above, the positive sign is for the stably-stratified flow, while the negative sign for the RBC. A stablystratified flow is stable, for which the u 2 /2 and b 2 /2 terms are the the kinetic and potential energies of the system, analogous to a harmonic oscillator. In RBC, the conserved quantity is also written as [u 2 − αgθ 2 /(dT /dz)]/2dr, where θ 2 /2 is called entropy. Note that θ 2 /2 is not the thermodynamic entropy that quantifies the degree of disorder at the microscopic scales. It is convenient to describe behaviour of turbulence flows in spectral or Fourier space since it captures the scale-by-scale energy and interactions quite well. In the next subsection, we describe the definitions used for such descriptions.

Equations in Fourier space
We rewrite Eqs. (17)- (19) in the Fourier space as In the above equations, i represents two things: √ −1 in front of the k i p(k, t)/ρ m term, and i = x, y, z in u i . Note that u(k), p(k), and θ(k) are the Fourier transforms of u, p, and θ respectively. The above equations are in terms of θ, but we can easily convert them as a function of ρ.
In the Fourier space, E u (k) denotes the kinetic energy spectrum, which is the sum of the kinetic energy of all the modes in a given shell (k − 1, k]. Similarly we define the spectra for the entropy and potential energy, which are denoted by E θ and E b respectively. They are computed using the following formulas:

Linear and nonlinear regimes
The behavior of buoyancy driven flows depends on the parameters and dimensionality.
Here we present a bird's-eye view of the observed states of RBC and stably-stratified flows. rolls. For Ra just above the onset, the instability saturates due to nonlinearity leading to the "roll" solutions. At larger Ra, the nonlinearity yields patterns and chaos [4,16,22,60,66,72]. For even larger nonlinearity, spatio-temporal chaos, weak turbulence, and strong turbulence emerge [60]. In this paper we will focus on only the strong turbulence regime.

2.7.2.
Stably-stratified flow For S = −1, the linearised version of Eqs. (20,21) yields internal gravity waves whose dispersion relation is where k ⊥ = k 2 x + k 2 y is the wavenumber component perpendicular to the buoyancy direction. Clearly ω = N for k = 0. These internal gravity waves persist for weak nonlinearity and inviscid case (ν = κ = 0). Strong nonlinearity has two kinds of generic behaviour: Strong stratification (Fr 1) suppresses the flow along the buoyancy direction and yields a quasi two-dimensional (2D) stratified flow; on the other hand, moderate and weak stratification (Fr 1) yields near isotropic turbulent flows. For Fr ≈ 1, Kumar et al. [46] obtained Bolgiano-Obukhov [10,70] scaling as predicted (to be described in Sec. 3.3.1). In this paper we focus on the Fr 1 regime.

Temperature profile and related equations
In this subsection we derive the properties of temperature fluctuations. For convenience we work with nondimensional variables.
Experiments and numerical simulations of RBC reveal that the horizontally averaged temperature T m (z) remains approximately a constant (≈ 1/2) in the bulk, and its value drops sharply in the thermal boundary layers [31,84], as shown in Fig. 2. The quantitative expression for T m (z) = T xy can be approximated as where δ T is the thickness of the thermal boundary layer, and xy represents averaging over the xy planes. A horizontal averaging of Eq.
as exhibited in Fig. 2. For thin thermal boundary layers, θ m (0, 0, k z ), which is the Fourier transform of θ m (z), is dominated by the contributions from the bulk. Hence The corresponding velocity mode, u z (0, 0, k z ) = 0 because of the incompressibility condition k · u(0, 0, k z ) = k z u z (0, 0, k z ) = 0. Also, u x (0, 0, k z ) = u y (0, 0, k z ) = 0 in the absence of a mean flow along the horizontal direction. Hence for the k = (0, 0, k z ) modes, the momentum equation yields or dσ m (z)/dz = ρ 0 αgθ m , and the dynamics of the remaining set of Fourier modes is governed by the momentum equation as where Hence, the modes θ m (0, 0, k z ) and σ m (0, 0, k z ) do not couple with the velocity modes in the momentum equation, but θ res and σ res do.
Equation (52) has strong implications on the scaling of the Reynolds and Nusselt numbers, which will be discussed in Sec. 4. In addition, the set of Fourier modes θ(0, 0, k z ) of Eq. (50) yields E θ (k) ∼ k −2 . This issue will be discussed in Sec. 3.

Other related systems
Several buoyancy-driven systems can be related to RBC. Here we list some of these systems.

Rayleigh-Taylor instability (RTI)
A fluid configuration with a denser fluid above a lighter fluid is unstable. The heavier fluid falls and the lighter fluid rises. After an initial stage of RTI, the flow develops significant nonlinearity and becomes turbulent [23]. We will discuss later that the turbulence phenomenology of RTI is similar to that of RBC.

Taylor-Couette flow
Two coaxial rotating cylinders create random flow at large Taylor number. This flow has been related to RBC with significant similarities in their phenomenology. See Grossmann et al. [39] for a review of such flows. [2] performed experiments in which a flow in a vertical tube is driven by an unstable density difference across the tube. They placed a brine solution at the top and distilled water at the bottom. This system has significant similarities with RBC [2]. Note however that the above system does not have walls or boundary layers at the top and bottom; this feature helps us study the ultimate regime quite conveniently. Exchanging the top and bottom containers will lead to behaviour similar to stably-stratified flows.

Turbulent exchange flow in a vertical pipe Arakeri and coworkers
2.9.4. Bubbly flow Bubbles are introduced in a tank in which turbulence is generated by an active grid [77]. Naturally this system has certain similarities with RBC.
In the next section, we will relate the turbulence behaviour of the above systems.

Spectra and fluxes of buoyancy-driven turbulence
It is convenient to describe behaviour of turbulence flows in spectral or Fourier space since it captures the scale-by-scale energy and interactions very well. In the next subsection, we describe the definitions used for such description.

Definitions
We can derive the time-evolution equation for E u (k) using Eq. (11) as [51,98] where T u (k) is the energy transfer rate to the shell k due to nonlinear interaction, and F B (k) and F ext (k) are the energy supply rates to the shell from the buoyancy and external forcing f u respectively, i.e., For brevity we set ρ m = 1. In Eq. (54), D(k) is the viscous dissipation rate defined by The kinetic energy (KE) flux Π u (k 0 ), which is defined as the kinetic energy leaving a wavenumber sphere of radius k 0 due to nonlinear interactions, is related to the nonlinear interaction term T u (k) as Under a steady state (∂E u (k)/∂t = 0), using Eqs. (54) and (58), we deduce that or In computer simulations, the KE flux, Π u (k 0 ), is computed using the following formula [29,97], Similarly, the potential energy (PE) flux Π ρ (k 0 ) is the potential energy leaving a wavenumber sphere of radius k 0 , which is computed using where b is defined in Eq. (30). For RBC, we replace u and b by nondimensional u and θ respectively. For a more detailed description of the energy transfers, we divide the wavenumber space into a set of wavenumber shells. The energy contents of a wavenumber shell of radius k and of unit width is denoted by E(k). The shell-to-shell energy transfer rate from the velocity field of the mth shell to the velocity field of the nth shell is defined as One of the most interesting problems in the field of buoyancy driven turbulence is the scaling of energy spectra and fluxes [56,79]. In the next section, we will review some of the theoretical results obtained for the aforementioned topic.

Turbulence phenomenology
3.2.1. Classical Bolgiano-Obukhov scaling for stably-stratified turbulence (SST): For the inertial range of isotropic hydrodynamic turbulence, Kolmogorov [43] first proposed a phenomenology according to which the energy spectrum is independent of the fluid properties and nature of large-scale forcing. He showed that the one-dimensional energy spectrum E(k) = K Ko Π 2/3 u k −5/3 in the inertial range of wavenumbers, where Π u (k) is the constant energy flux, and K Ko is the Kolmogorov's constant.
Buoyancy (forcing) act at all scales, hence Kolmogorov's theory may not work for the buoyancy-driven turbulence. In this section we will describe how the buoyancy affects the energy spectra and fluxes of the buoyancy-driven flows. For stable stratification, Bolgiano [10] and Obukhov [70] argued that the KE flux Π u (k) is depleted at different length scales due to the conversion of KE to PE via buoyancy (u z ρg). Subsequently, Π u (k) decreases with k, and E u (k) is steeper than that predicted by Kolmogorov's theory; we refer to the above as BO phenomenology or scaling. According to this phenomenology, for L B l L, buoyancy is important and the buoyancy term is balanced by the nonlinear term [ρg ≈ (u · ∇)u]. Here L B is the Bolgiano scale [10] and L is the large length scale or the box size. The force balance at wavenumber k = 1/l yields According to BO phenomenology, PE has a constant flux, i.e., Π ρ ≈ ku k ρ 2 k ≈ ρ . Hence, where c i 's are constants. At smaller length scales (k > k B ), where k B = 1/L B is the Bolgiano wavenumber, Bolgiano [10] and Obukhov [70] argued that the buoyancy is relatively weak, hence Kolmogorov-Obukhov (KO) scaling is valid in this regime, i.e., where K Ba is the Batchelor's constant. A comparison of Π u (k) of Eq. (69) with that of Eq. (73) yields the crossover wavenumber k B as The BO phenomenology implicitly assumes isotropy in Fourier space, which is a tricky assumption. For BO scaling, the gravity must be strong enough to compete with the nonlinear term u · ∇u, but not too strong to make the flow quasi two-dimensional (quasi-2D). This corresponds to Fr ≈ 1 regime. A large number of earlier explorations in SST have been for Fr 1 regime, see for example, Lindborg [53], Brethouwer [13], and Bartello and Tobias [3]. SST can be broadly classified in three regimes. Note that nonlinearity is strong (Re 1) for turbulent flows.
(ii) Moderate gravity (Ri ≈ 1): Comparable strengths of gravity and nonlinearity yields nearly isotropic turbulence with BO scaling, as described earlier.

Generalization of Bolgiano-Obukhov scaling to RBC:
Using mean field theory approximation, Procaccia and Zeitak [78] argued that the Bolgiano-Obukhov scaling is applicable to convective turbulence. Later, L'vov [58] assumed that in convective turbulence, the kinetic energy is converted to the potential energy and therefore, favored BO scaling. L'vov and Falkovich [59] employed a differential model for energy and entropy fluxes in k-space and found that the BO scaling is valid for convective turbulence. Rubinstein [80] employed renormalization group analysis to RBC and observed that the renormalized viscosity ν(k) ∼ k −8/5 , E u (k) ∼ k −11/5 , and E ρ (k) ∼ k −7/5 . Based on these observations Rubinstein claimed BO scaling for RBC.
The aforementioned theories had profound influence in the field, and a large number of analytical, experimental, and numerical works attempted to verify these ideas. Recently Kumar et al. [46] showed that the BO scaling does not describe RBC turbulence since the energy supply by buoyancy in RBC is very different from that in stably stratified flow. We will provide these arguments below.

3.2.3.
A phenomenological argument based on kinetic energy flux: Kumar et al. [46] and Verma et al. [100,101] presented a phenomenological argument based on the KE flux to derive a spectral theory of buoyancy-driven turbulence. Equation (60) provides important clues on the energy spectrum and flux of the buoyancy-driven flows. Here we list three possibilities for the inertial range (k f < k < k d ), where k f is the forcing wavenumber, and k d is the dissipation wavenumber: (i) For the inertial range of hydrodynamic turbulence, F B (k) = 0 and D(k) → 0, therefore Π u (k + ∆k) ≈ Π u (k) and E u (k) ∼ k −5/3 , which is a prediction of the Kolmogorov's theory [43]. Note that F ext (k) = 0 in the inertial range.
Reprinted with permission from Kumar et al. [46].
(ii) For the stably stratified flows, as argued by Bolgiano [10] and Obukhov [70], in the k f < k < k B wavenumber band, buoyancy converts the kinetic energy of the flow to the potential energy, i.e., F B (k) < 0. Hence, Eq. (60) predicts that Π u (k) will decrease with k in this wavenumber range, as shown in Fig. 3(a). In the wavenumber range, k B < k < k d , buoyancy becomes weaker, hence Π u (k) ≈ constant.
(iii) For RBC, buoyancy feeds the kinetic energy, hence F B (k) > 0. Therefore we expect the KE flux Π u (k) to increase. Numerical simulation of Kumar et al. [46] however show that the energy supplied by buoyancy is dissipated by the viscous force, i.e., F B (k) ≈ D(k). Hence Π u (k) ≈ constant in the inertial range, and we recover Kolmogorov's spectrum for RBC. Note that L'vov [58] argued that F B (k) < 0, which is not the case.

Modeling and field theory:
Researchers [30,44,52,61,62,107] employed fieldtheoretic techniques to understand the physics of turbulent fluid. In field theory, the nonlinear terms of the equations are expanded perturbatively. Some of the popular fieldtheoretic techniques are direct interaction approximation (DIA) [44,52], renormalization group analysis [30,44,61,62,107], mean field approximation [78], etc. Field theory has been applied to buoyancy-driven flows as well.
As described in Sec. 3.2.2, that Procaccia and Zeitak [78] employed mean field approximation to convective turbulence and obtained BO scaling. Rubinstein [80] used Yakhot-Orszag's [107] renormalization group procedure and proposed an isotropic model for convective turbulence. His results are consistent with that of Procaccia and Zeitak [78]. Recently, using self-consistent field theory, Bhattacharjee [6] obtained E u (k) ∼ k −13/3 for RBC in the infinite Prandtl number limit. Bhattacharjee [5] used the global energy balance for the stratified fluid and argued that the BO scaling could be observed in stably stratified flow at high Richardson number. In addition, he also added the possibility of BO scaling for RBC in some range of Prandtl numbers.
In the next section, we will present numerical results for the stably stratified turbulence and Rayleigh-Bénard convection.

Numerical analysis of buoyancy-driven turbulence
Stably stratified turbulence: Researchers simulated the stably-stratified turbulence for the three regimes described in Subsection 3.2.1. First we discuss the results for a strong gravity that corresponds to Ri 1 or Fr 1. Such configurations are observed in some regimes of plantery and stellar atmospheres. Strong gravity makes such flows quasi-2D with dual scaling, k −3 and k −5/3 . In this regime, Lindborg [53], Brethouwer [13], and Bartello and Tobias [3] showed that the spectra of the horizontal KE and PE follow k −5/3 ⊥ scaling, while the energy spectrum of the vertical velocity follows k −3 . Vallgren et al. [95] included rotation in their simulation and obtained KE spectra k −3 and k −5/3 for two different wavenumber bands.
For weak stratification (Ri 1), Kumar et al. [46] performed a 3D stably-stratified turbulence simulation and reported Kolmogorov's spectrum for the kinetic energy as expected. Kumar et al. also studied the moderate stratification regime and reported BO scaling, which will be described below. In this paper we focus on the results for Fr ≈ 1 since they are recently observed. Normalized E ρ (k) Kumar et al. [46] simulated stably stratified flows in a cubical box of size (2π) 3 with periodic boundary conditions at all the walls. They forced the small wavenumber modes randomly to achieve a steady state. The parameters of their simulations are Ra = 5×10 3 and Pr = 1 that yields Ri = 0.01 and Fr = 10. Figure 4(a) exhibits the normalized KE spectra-E u (k)k 11/5 for the BO scaling, and E u (k)k 5/3 for the KO scaling. The numerical data fits better with the BO scaling for than the KO scaling, thus confirming the BO phenomenology for the stably-stratified turbulence when Fr ≈ 1. This is also verified by the PE spectrum as shown in Fig. 4(b) in which E ρ (k)k 7/5 provides a better fit to the data than E ρ (k)k 5/3 .   [46] also computed the energy supply rate by buoyancy, F B (k), and the viscous dissipation spectrum, D(k), which are illustrated in Fig. 6. Note that F B (k) < 0, as argued in BO phenomenology. The Bolgiano wavenumber k B of Eq. (75) is approximately 8.5, which is only 3 to 4 times smaller than k d , wavenumber where the dissipation range starts. Therefore Kumar et al. [46] did not observe a definitive crossover from k −11/5 to k −5/3 in their simulations.
The aforementioned observations demonstrate applicability of the BO scaling for SST with a moderate stratification.

Rayleigh-Bénard Convection:
A large number of numerical simulations have been performed with an aim to identify which among the two, BO or KO, scaling is applicable to RBC. Grossmann and Lohse [33] using simulation for Pr = 1 under Fourier-Weierstrass approximation and reported Kolmogorov's scaling. Based on periodic boundary condition, Borue and Orszag [11] andŠkandera et al. [88] reported KO scaling for the velocity and temperature fields. Kerr [42] reported the horizontal spectrum as a function of horizontal wavenumber and observed Kolmogorov's spectrum. Verzicco and Camussi [103], and Camussi and Verzicco [19] showed BO scaling using the frequency spectrum of real space probe data. Kaczorowski and Xia [41] reported KO scaling for the longitudinal velocity structure functions, but BO scaling for the temperature structure functions in the centre of a cubical cell. Kumar et al. [46] computed E u (k) and Π u (k), and showed Kolmogorov-like behaviour for RBC, i.e., E u (k) ∼ k −5/3 and Π u (k) ∼ const. In this paper we present the above quantities for 4096 3 resolution and very high Ra that unambiguously demonstrates KO scaling for RBC. We also report the shell-to-shell energy transfers and the ring spectrum for RBC that show close resemblance with the hydrodynamic turbulence. We performed RBC simulations in a unit box with 4096 3 grid for Pr = 1 and Ra = 1.1 × 10 11 . For the velocity field, we employed the free-slip boundary condition at the top and bottom plates, and periodic boundary condition at the side walls. The temperature field satisfies conducting boundary condition at the top and bottom plates, and the periodic boundary condition at the side walls. We computed the spectra and fluxes of the KE and the entropy (θ 2 /2) using the steady state data. Figure 7(a) exhibits the KE spectra normalized with k 11/5 and k 5/3 . The plots indicate that in the wavenumber band 15 < k < 600 (inertial range), the shaded region of the figure, the KO scaling fits better than the BO scaling.
We exhibit the KE and entropy fluxes in Fig 7(b). We observe that the kinetic energy flux Π u (k) remains constant in the inertial range, a band where E u (k) ∼ k −5/3 . Thus we claim that the convective turbulence exhibits Kolmogorov's power law in the inertial range. We also computed F B (k), Π u (k), and dΠ u (k)/dk as further tests. According to Fig. 8(a) F B (k) > 0 in the inertial range, consistent with the discussion of Sec. 3.2.3 and Fig. 3(b), and it approximately balances D(k). Therefore, (59)]. The constancy of Π u (k) yields E u (k) ∼ k −5/3 , consistent with the energy spectrum plots of Fig. 7(a). Fig. 8    shows that [dΠ u (k)/dk]/Π u (k) 1 in the inertial range consistent with the constant Π u (k). Interestingly, We also compute the shell-to-shell energy transfers [Eq. (63)] using the steadystate data of our simulation. We divide the Fourier space into 40 concentric shells; the inner and outer radii of the nth shell are k n−1 and k n respectively with k n = {0, 2, 4, 8, 8 × 2 s(n−3) , ..., 6432}, where s = (1/35) log 2 (804). The radii of the inertialrange shells are binned logarithmically due to the power law physics of RBC in the inertial range. In Fig. 9(a) we exhibit the shell-to-shell energy transfers with the indices of the x, y axes representing the receiver and giver shells respectively. The plot indicates that mth shell gives energy to (m + 1)th shell, and it receives energy from the (m − 1)th shell. Thus the energy transfer in RBC is local and forward, very similar to hydrodynamic turbulence. This result is consistent with the energy spectrum and flux studies described earlier.
Convective flows are expected to be anisotropic due to buoyancy; hence it is important to quantify anisotropy using quantities that are dependent on the polar angle, the angle betweenẑ and k. For the same, we divide a wavenumber shell into rings [67]. The energy contents of the rings are called ring spectrum E(k, β), where β represents the sector index for the polar angles (for details see Nath et al. [67]). The ring spectrum E(k, β), depicted in Fig. 9(b), shows that the flow is nearly isotropic, again similar to hydrodynamic turbulence. These results clearly demonstrate that the turbulent convection for Pr = 1 has a very similar behavior as hydrodynamic turbulence. The temperature fluctuation however exhibit a unique behaviour. As illustrated in Figure 10, we observe dual branches for the entropy spectrum (E θ (k)). The upper branch varies as k −2 because θ(0, 0, k z ) ≈ −1/(πk), as discussed in Sec. 2.8. The lower branch shows neither KO (k −5/3 ) nor BO (k −7/5 ) spectrum. Note that both the branches of entropy spectrum generate a constant entropy flux Π θ (k) (see Fig. 7(b)), and the modes θ(0, 0, k z ) also participate in energy transfers.

Experimental results
For stably-stratified flows, there are not many laboratory experiments to verify BO phenomenology. However, scientists have measured the KE spectrum of the Earth's atmosphere and relate it to the theoretical predictions. Most notably Gage and Nastrom [32] observed a combination of k −3 and k −5/3 energy spectra. Some researchers attribute the k −3 spectrum at lower wavenumbers to the two-dimensionalization of the flow, while k −5/3 spectrum at larger wavenumbers to the forward cascade of kinetic energy; yet these issues are still unresolved. These features are expected to arise for Fr 1. There are a significant number of laboratory experiments on RBC, with some favouring the BO scaling [24,110], while some others in support of the KO scaling [27]. In most convective experiments, the velocity field, u z (r, t), and/or the temperature field, T (r, t), are probed near the lateral walls of the container. For such experiments, the Taylor's hypothesis [49,83,93] is invoked to relate the frequency power spectrum E(f ) to the one-dimensional wavenumber spectrum E(k); this connection is under debate due to the absence of any constant mean velocity field. Researchers [50,92,108, 109] employ 2D particle image velocimetry (PIV) for high-resolution visualization and computation of an approximate energy spectrum under the assumption of homogeneity and isotropy, which is not strictly valid in convection [67]. In summary, on the experimental front, there is no convergence on which of the two scaling, BO of KO, is valid. For details refer to the review papers [1,56].

Turbulence in thermal boundary layer and in two dimensions
A burning question is whether KO or BO scaling are applicable to the boundary layers of RBC. The flux arguments of Sec. 3.2.3 provide some insights into the dynamics of boundary layers. Here, typically u z u ⊥ , hence the flow is quasi-2D, and we expect an inverse cascade of KE. Using Π u (k) < 0, F B (k) > 0, and dΠ u (k)/dk ≈ F B (k), we may argue that Π u (k) may increase with k as shown in Fig. 11. An application of scaling arguments of Sec. 3.2.1 may yield spectra and fluxes according to Eqs. (67)(68)(69)(70), i.e., Bolgiano-Obukhov scaling for k < k B . For k > k B , the KE spectrum may exhibit a mixture of k −5/3 (regime of inverse cascade of energy) and k −3 (regime of forward cascade of enstrophy) depending on where the effective forcing band lies in relation to k B . Thus, in the boundary layer, RBC may exhibit BO scaling, and it needs to be investigated carefully using numerical simulations and experiments. The aforementioned scaling arguments may work for 2D RBC (xz plane in which the buoyancy is along the z direction), as well as in quasi 2D RBC (when L x L y ). Toh and Suzuki [94] simulated 2D RBC and reported E u (k) ∼ k −11/5 and Π u (k) ∼ −k −4/5 in line with the above arguments. Calzavarini et al. [18] also reported similar results in their structure function computations.

Turbulence in Rayleigh-Taylor instability (RTI)
Chertkov [23] proposed that a fully-developed 3D RTI will exhibit Kolmogorov's spectrum due to the Rayleigh-Taylor pumping at large scales. Boffetta et al. [9] observed this behaviour in their numerical simulations. Chertkov [23] however does not take into account the buoyancy at all scales (see Sec. 3.2.3). In a quasi-2D box (L y L x ), Boffetta et al. [8] show coexistence of BO and KO scaling (k −11/5 and k −5/3 ), consistent with the arguments of Sec. 3.5.

Turbulence in miscellaneous systems
Scientists have studied spectra of the velocity and the scalar field in other buoyancydriven systems. Pawar and Arakeri [76] performed experiment on the vertical tube described in Sec. 2.9.3. They observed that the velocity field exhibits k −5/3 spectrum, while the scalar spectrum is closer to k −7/5 .
Prakash et al. [77] studied the energy spectrum of the bubbly turbulence using an experiment. For the velocity field, they reported k −5/3 energy spectrum for k < 1/b, and k −3 for k > 1/b where b is the bubble size. They argue that that the large and intermediate scales exhibit k −5/3 spectrum due to the standard Kolmgorov's argument. For k > 1/b, Prakash et al. [77] explained the k −3 energy spectrum by invoking equipartition between the energy dissipation and energy feed by the buoyancy. For this system it may be interesting to investigate the energy spectrum using the flux arguments.
The turbulent Taylor-Couette flow [39] may exhibit spectral behaviour similar to RBC since both the systems are unstable with similar energetics (see Secs. 3.2.3 and 3.3.2). We believe that the Non-Bouusinesq convective flows may also exhibit Kolmogorov-like spectrum for weak compressibility since here too the thermal plumes feed the kinetic energy, as in RBC.

Turbulence in small and large Prandtl number RBC
In Sec. 3.2.3 we derived the spectra and fluxes of the velocity and temperature for RBC when Pr ∼ 1. These arguments are not applicable to RBC with extreme Prandtl numbers. However, we can easily deduce the spectrum for very small and very large Pr's as follows. These computations have been first reported in [65] and [75] respectively. In RBC with zero or small Prandtl numbers, thermal diffusivity κ → 0 that leads to u z (k) ∼ θ(k)/(κk 2 ) [65]. Hence, the buoyancy, which is proportional to θ(k), is dominant at small wavenumbers. Therefore, the assumption of the Kolmogorov's phenomenology that the forcing acts at large length scales is valid, and we expect the Kolmogorov's phenomenology for the hydrodynamic turbulence to be applicable to RBC with Pr → 0. Mishra and Verma [65] verified the above phenomenology using numerical simulations.
In the limit of infinite Prandtl number (ν → ∞), the momentum equation is linear [75]. However if the Péclet number is large, the temperature equation is nonlinear and it yields an approximate constant entropy flux. Using scaling arguments, Pandey et al. [75] derived that for infinite and large Pr, E u (k) ∼ k −13/3 . They also verified the above scaling using numerical simulations.

Simulation of turbulent convection in a periodic box and shell model
Borue and Orszag [11],Škandera et al. [88], Lohse and Toschi [55], and Calzavarani et al. [17] simulated turbulent thermal convection in a periodic box. They simulated Eqs. (17,18) under a gradient dT /dz. In the absence of boundary layers, the velocity and temperature fields exhibit k −5/3 spectra [11,88]. In addition, the Nusselt number Nu ∼ Ra 1/2 [17,55], which is expected in the ultimate regime when the effects of boundary layers are negligible. Note that the temperature spectrum for the periodic box is very different from that with conducting walls that exhibit dual spectra. It is important to note that turbulent thermal convection in a periodic box is numerically unstable; the system exhibits steady behaviour for carefully chosen set of initial conditions.
Direct numerical simulation of turbulent systems is quite demanding due a large number of interacting Fourier modes. Therefore, scientists often use shell models, which are based on much fewer number of modes. Brandenburg [12], Lozhkin and Frick [57], Mingshun and Shida [63], Ching and Cheng [26], and Kumar and Verma [47,48] constructed shell models for buoyancy-driven turbulence. The advantage of the shell model of Kumar and Verma [47] is that it describes both turbulent stably-stratified and convective flows using a single set of equations. It also enables flux computation of the kinetic energy and density. Kumar and Verma [47] showed that the results of the shell model are consistent with the DNS results described earlier.

Concluding remarks on the energy spectrum
We summarise the important results of this section in the following. (ii) Turbulence in RBC has significant similarities with hydrodynamic turbulence. For example, the KE flux is nearly constant in the inertial range; the shell-to-shell energy transfer is local and forward; the ring spectrum exhibits a near isotropy in Fourier space. The constant KE flux is due to the near cancellation between the KE supply by buoyancy and the viscous dissipation rate.
(iii) Under nearly isotropic conditions (when Froude number is of the order of unity), the stably-stratified turbulence exhibits Bolgiano-Obukhov scaling.
(iv) The temperature fluctuations exhibits dual spectra, with the upper branch scaling as k −2 . In Sec. 2.8 we discussed the origin of k −2 spectrum in terms of the structures of the boundary layers and the bulk.

Modelling of large-scale quantities of RBC
In this section we quantify the large-scale quantities of RBC, namely the Nusselt and Reynolds numbers. Many researchers have worked on this problem; for details and references, refer to the review articles [1,7,56,86]. Despite complexities of the flow, RBC exhibits certain universal behaviour; in the turbulent limit, Pe ∼ √ RaPr, but in the viscous regime, Pe ∼ Ra 3/5 [34,73]. The Nusselt number however scale as Nu ∼ Ra β with β ranging from 0.27 to 0.33. Researchers have attempted to explain the above behaviour. For brevity, in this review we only discuss the models of Grossmann and Lohse (GL) [34][35][36][37][38] and that of Pandey et al. [73].

Grossmann-Lohse model
Grossmann and Lohse (GL) [34][35][36][37][38]90] derived the formulas for Nu(Ra, Pr) and Re(Ra, Pr) by exploiting the fact that the global viscous dissipation rate, u , and thermal dissipation rate, T , get contributions from the bulk and boundary layers, i.e., where BL and bulk denote the boundary layer and the bulk respectively. They invoked the exact relations of Shraiman and Siggia [85] for the global viscous and thermal dissipation rates (see Eqs. (38,39)), and estimated the aforementioned contributions of the boundary layers and the bulk to u and T in various Ra-Pr regimes. For Pr ≈ 1 and very large Ra they used u,bulk = U 3 /d and T,bulk = U ∆ 2 /d, but for extreme Prandtl numbers, these estimates get altered by the boundary layer widths.
Using the above formulae, GL computed the Nusselt and Reynold numbers as a function of Ra and Pr that agree with presently available experimental and numerical simulation results quite well [1].

An alternate derivation of Péclet number
Recently Pandey et al. [73] and Pandey & Verma [74] provided an alternate derivation of Péclet number. Note that Pe = RePr. Pandey et al. [73] analysed the rms values of various terms of the momentum equation, which are exhibited in the schematic diagram of Fig. 12. Under statistical steady state ( ∂u/∂t ≈ 0), Pandey et al. observed that in the turbulent regime, the acceleration u · ∇u is primarily provided by the pressure gradient −∇σ, and the buoyancy and viscous terms are relatively small. The above features are consistent with similarities between the turbulence in RBC and hydrodynamics (see Sec. 3.3.2). However, in the viscous regime (Re 1), −∇σ is small, and the buoyancy and viscous terms cancel each other resulting in a very small acceleration of the fluid.
Dimensional analysis of the momentum equation yields where c i 's are dimensionless coefficients defined as Pandey et al. [73] observed c i 's to be functions of Ra and Pr that yields interesting and nontrivial scaling relations. It is important to contrast this behaviour with free turbulence (without walls) where c i 's are constants. Multiplication of Eq. (80) with d 3 /κ 2 yields Turbulent Viscous (a) (b) Figure 12. The relative strengths of the forces acting on a fluid parcel. In the turbulent regime, the acceleration u·∇u is provided primarily by the pressure gradient.
In the viscous regime, the buoyancy and the viscous force dominate the pressure gradient, and they balance each other. Reprinted with permission from Pandey and Verma [74].
where Pe = U d/κ is the Péclet number. The solution of the above equation is using which Pe can be computed as a function of Ra and Pr.
In the turbulent regime, the viscous term of Eq. (82) can be ignored, hence This limit is applicable when The scaling for the viscous regime is obtained by equating the buoyancy and viscous terms of the momentum equation that yields Pandey et al. [73] computed c i 's using the RBC simulation data for Pr = 1, 6.8, 10 2 , 10 3 and Ra from 10 6 to 5 × 10 8 . These simulations were performed for no-slip boundary condition at all the walls using a finite volume solver OpenFoam [71]. They reported the following functional form for c i 's   [106]. Reprinted with permission from Pandey and Verma [74].
The errors in the above exponents are 0.01, except for the Ra exponent of c 4 that has error of the order of 0.10. In Fig. 13 which is independent of Pr, consistent with the results of Silano et al. [87], Horn et al. [40], and Pandey et al. [75]. For the turbulent regime, Eq. (84) yields For mercury (Pr ≈ 0.025) as an experimental fluid, Cioni et al. [28] observed that Re ∼ Ra 0.424 , which is close to the predicted exponent of 0.38 discussed above. The range of Rayleigh numbers in the experiment of Cioni et al. [28] is 5 × 10 6 ≤ Ra ≤ 5 × 10 9 that is consistent with the turbulent regime estimated above (Ra 10 6 Pr). The aforementioned results are in general agreement with those of Grossmann and Lohse [34][35][36][37][38].

Scaling of Nusselt number and dissipation rates
The Nusselt number, a measure of the convective heat transport [1,25,105], is quantified as where V stands for a volume average, u z = u z d/κ, θ res = θ res /∆, and C uθres is the normalized correlation function between the vertical velocity and the residual temperature fluctuation [102]: The deviation of the exponent from 1/2 in the ultimate regime [45] is due to the nontrivial scaling of C uθres , u z , and θ res . We observe that C uθres , and the rms values of u z and θ res scale with Ra in such a way that Nu ∼ Ra 0.32 ; without these corrections, Nu ∼ Ra 1/2 in the turbulent regime. For details refer to Pandey et al. [73] and Pandey and Verma [74]. In hydrodynamic turbulence, the viscous dissipation rate u ≈ U 3 /d. However this is not the case in RBC, primarily due to walls or boundary layers. Using numerical data, Verma et al. [102] and Pandey et al. [73] have shown that u ∼ Ra 1.32 or See Fig. 14 for illustration for Pr = 1. One of the exact relations of Shraiman and Siggia [85] yields Substitution of Pe ∼ Ra 0.51 and u ∼ (U 3 /d)Ra −0.21 yields Nu ∼ Ra 0.32 . These arguments show that the reduction of the viscous dissipation rate could be a reason for the deviation of the observed scaling Nu ∼ Ra 0.32 from Nu ∼ Ra 1/2 corresponding to the ultimate regime. We conclude this section with a remark that the walls or the boundary layers affect the scaling of Péclet and Nusselt numbers significantly.  and Ahlers [15], Xi and Xia [104], and Sugiyama et al. [91] observed that the vertical velocity near the lateral wall switches sign randomly. Deciphering how the reversals take place is an interesting puzzle, and it is not yet fully solved. In this section we briefly describe the present status of the field. It is believed that the flow reversals are caused by the nonlinear interaction among the large-scale structures of the flow. For a closed cartesian box, these structures can be conveniently described by the small-wavenumber Fourier modes [20,21]. This description is useful even for no-slip boundary conditions since the flow structures inside the boundary layers contribute to the large wavenumber modes. For a cylindrical geometry, partial information about the flow structures can be obtained by computing the azimuthal Fourier modes corresponding to the velocity field measured at various angles near the later walls [15,64,104]. Here we summarise the main results on the properties of flow reversals.
(i) During a reversal, the amplitude of the most dominant large-scale mode vanishes, while that of the secondary mode rises sharply. Chandra and Verma [20,21] reported that during a reversal in a unit cartesian box, the Fourier mode (1, 1) vanishes, while the mode (2, 2), corresponding to the corner rolls, become the most dominant mode [20,21]. See Fig. 1 of Chandra and Verma [20]. This numerical result is consistent with the experimental results of Sugiyama et al. [91].
(ii) The nature of dominant structures depends on the box geometry and boundary conditions. For example, for a box of size 2 × 1, under the no-slip boundary condition, (2, 1) and (2,2) are the primary and secondary modes respectively. However, under the free-slip boundary condition, the corresponding modes are (2, 1) and (1, 1) respectively [14,99].
(iii) Verma et al. [99] have constructed a group-theoretic argument to decipher the reversing and non-reversing modes during a reversals. The structure of the groups is related to the Klein group.
(iv) Convection in a cylinder exhibit reversals that have similar behaviour as above.
Brown and Ahlers [15] termed such reversals as cessation-led reversals. Note however that during a cessation-led reversal, the secondary modes become significant, hence the kinetic energy does not vanish.
(v) Cylindrical convection exhibits another kind of flow reversals, called rotation-led reversals, in which the large-scale structure rotates azimuthally [15,64,104]. This rotation is due to the azimuthal rotation symmetry of the system. Such phenomena is also observed in a cylindrical annulus [68].
The magnetic field reversals in dynamo, and the velocity field reversals in Kolmogorov-flow also involve nonlinear interactions among the large-scale structures of the flow. Thus, these reversals share certain similarities with the flow reversals of RBC.

Summary
In this short review we describe the spectral and large-scale properties of buoyancydriven turbulence-stably-stratified flows and Rayleigh-Bénard convection. A summary of the results covered in this review is as follows: (i) The stably-stratified turbulence (SST) is nearly isotropic for Froude number Fr 1.
(ii) For Fr 1, SST exhibits Kolmogorov scaling, i.e. E u (k) ∼ k −5/3 , due to the dominance of the nonlinear term over the buoyancy.
(iii) For Fr 1, SST is quasi two-dimensional, and the kinetic-energy spectrum exhibits a combination of k −5/3 and k −3 .
(iv) Turbulence in Rayleigh-Bénard convection (RBC) has strong similarities with the hydrodynamic turbulence, e.g, it exhibits constant energy flux and k −5/3 energy spectrum in the inertial range.
(v) In RBC turbulence, the pressure gradient accelerates the flow, while the buoyancy is balanced by the viscous dissipation. This observation is consistent with the Kolmogorov-like phenomenology observed for RBC.
(vi) The aforementioned phenomenology of RBC turbulence is expected to work for other buoyancy-driven flows in which buoyancy feeds the kinetic energy. Some of the examples of such flows are bubbly turbulence, non-Boussinesq thermally-driven flows in stars, turbulent buoyancy-driven exchange flows in a vertical pipe [2], etc.
(vii) The scaling of the Reynolds and Nusselt numbers of RBC are well described by the models of Grossmann and Lohse [34][35][36][37][38], and by the recent formula of Pandey et al. [73].
We conclude this paper with a remark that several issues remain unresolved in buoyancy-driven turbulence, for example, the scaling of Nusselt number, the role of the boundary layer and bulk to turbulence, existence of the ultimate regime. Highresolution simulation, advanced experiments, and carefully modelling may resolve these outstanding questions in future.