Temperature Statistics in Turbulent Rayleigh-B\'enard Convection

Rayleigh-B\'enard convection in the turbulent regime is studied using statistical methods. Exact evolution equations for the probability density function of temperature and velocity are derived from first principles within the framework of the Lundgren-Monin-Novikov hierarchy known from homogeneous isotropic turbulence. The unclosed terms arising in the form of conditional averages are estimated from direct numerical simulations. Focusing on the statistics of temperature, the theoretical framework allows to interpret the statistical results in an illustrative manner, giving deeper insight into the connection between dynamics and statistics of Rayleigh-B\'enard convection. The results are discussed in terms of typical flow features and the relation to the heat transfer.


Introduction
Rayleigh-Bénard convection is a paradigm of a pattern forming system far from equilibrium. Convective fluid motion in a vessel is induced by a vertical temperature gradient between the bottom and top boundaries due to buoyancy forces. In dependence on this temperature gradient, the geometry of the experiment and the fluid properties, a whole zoo of instabilities has been observed ranging from laminar, spatially coherent convective motion over spatially ordered but temporally chaotic up to highly turbulent fluid motion. We refer the reader to reviews available on the topic [1,2].
Recently, much efforts have been devoted to the analysis of turbulent Rayleigh-Bénard (RB) convection both by experimental as well as theoretical means [3,4]. Direct numerical simulations allow one to consider the dynamical and statistical properties of RB turbulence and the transitions between different types of flows in fine detail.
It is obvious that the analysis of turbulent convective fluid motion has to be based on a combination of tools from dynamical systems theory, statistical physics, and the theory of stochastic processes. A necessary step is the statistical formulation of the underlying basic fluid dynamic equations, which for the most simple case are the Oberbeck-Boussinesq equations for the velocity field u(r, t), the temperature field T (r, t), and the pressure field p(r, t): ∂ ∂t + u · ∇ T (r, t) = ∆T (r, t) ∂ ∂t + u · ∇ u(r, t) = − ∇p(r, t) + Pr [∆u(r, t) + Ra T (r, t)e z ] (1.1) ∇ · u(r, t) = 0 The equations have been nondimensionalized using the Rayleigh number Ra = αg∆T h 3 νκ , which is a dimensionless measure of the temperature gradient across the fluid layer (with thermal expansion coefficient α, gravitational acceleration g, outer temperature difference ∆T , and distance of top and bottom plate h), as well as the Prandtl number Pr = ν κ as the ratio of kinematic viscosity ν to heat conductivity κ of the fluid. Thus, the vertical spatial coordinate obeys z ∈ [0, 1], and the boundary conditions of the temperature at the bottom and top plate are T (z = 0) = 1 2 and T (z = 1) = − 1 2 . For the velocity, no-slip boundary conditions u(z = 0) = u(z = 1) = 0 are assumed. These equations are solved numerically by a suitably designed penalization approach, described in section 5. A snapshot of the temperature field is exhibited in figure 1. The statistical analysis is based on joint probability density functions (PDFs) for the temperature and the velocity at a single point in space and time. The basic fluid dynamic equations require the validity of certain relations among these PDFs. From these relations, corresponding expressions relating the various moments of the fields can be derived. For the case of incompressible turbulence, these relations have been formulated by Lundgren, Monin and Novikov [6,7,8], and Ulinich and Lyubimov [9] and are sometimes known as the Lundgren-Monin-Novikov (LMN) hierarchy. They are directly related to Hopf's functional equation, which can be viewed as the basic statistical formulation of the Navier-Stokes equation in the Eulerian framework [10]. Similar relations can be derived for the corresponding Lagrangian quantities [11]. It is evident that an analogous treatment is feasible for the Oberbeck-Boussinesq equations.
In the present article we will use this approach in order to analyze the single-point temperature probability density function for stationary turbulent RB convection. Our analysis combines direct numerical simulations with the relation of the LMN hierarchy for the single-point PDF. The result is a partial differential equation for the temperature PDF. The derivation of this relation is outlined in section section 2.
The aim is to formulate an equation which characterizes this PDF. Starting from the Oberbeck-Boussinesq equations, we derive an evolution equation for the singlepoint joint probability density function of velocity and temperature along the lines of Lundgren, Monin and Novikov [6,7,8] for the case of incompressible turbulence. This equation is unclosed due to the fact that it contains unclosed expressions which can be related to fluid pressure, viscous dissipation and heat diffusion. However, these expressions can be treated by introducing conditional averages, which can be extracted from direct numerical simulations. This leads to a partial differential equation for the joint temperature-velocity PDF. A similar approach has been performed by Novikov [12,13], and more recently by Wilczek et al. for the PDFs of vorticity [14] and velocity [15] for stationary, isotropic turbulence. On the other hand, modeling of unclosed terms is also a possible method, as performed by e.g. Pope [16,17]. As we shall indicate, the analysis of the evolution equation for the temperature PDF yields a comprehensive description of the dynamical processes in RB convection.
The article is structured as follows: In section 2 we will derive the evolution equation for the temperature-velocity joint PDF. In section 3, we reduce the joint PDF to the temperature PDF and make use of statistical symmetries to cut down the complexity of the evolution equation. Then we present a descriptive way to deal with this equation involving the method of characteristics. Section 4 briefly discusses connections to the Nusselt number, relevant for the heat transport. These theoretical results are complemented by results from direct numerical simulations, which will be discussed in section 5, followed by a summary in section 6.

Single-Point PDF
We are interested in the joint temperature-velocity probability distribution f (τ, v; r, t) and want to derive the corresponding evolution equation. Formally, the probability density function is obtained as a suitable average over the so-called fine-grained probability distribution f (τ, v; r, t) = δ(τ − T (r, t)) δ(v − u(r, t)). (2.1) It is important to distinguish between sample space variables τ , v and the corresponding realizations of the temperature and velocity fields T (r, t), u(r, t). Therefore, one could think of f as the PDF of one particular realization of the fields. Also, the notation of the arguments in f (τ, v; r, t) emphasizes the difference between the sample space variables τ , v and the coordinates r, t -the PDF is normalized with respect to the sample space variables, the coordinates are just parameters. The full probability density function is now obtained as an ensemble average over all possible realizations of the temperature and velocity fields: The brackets · denote the ensemble average, in contrast to spatial averages · V and · A over the whole fluid volume, or a horizontal plane at height z, respectively. The definition of the fine-grained PDF (2.1) can be differentiated with respect to the space and time variables, giving as the spatial gradient, and an analogous equation for the temporal derivative ∂ ∂t f . Note that the operators ∂ ∂τ and ∇ v act on f . Now, multiplying (2.3) by u i and adding the temporal derivative allows us to make use of the basic Oberbeck-Boussinesq equations (1.1). This results in the desired evolution equation for the fine-grained PDF: Performing the ensemble average of this equation in order to arrive at an equation for the full PDF (2.2), one encounters the closure problem of turbulence, since the unclosed averages ∆T f , −∇p f and ∆u f show up. The LMN hierarchy ansatz would mean to treat these terms via a coupling to the two-point PDF; the evolution equation of this two-point PDF would in turn introduce a coupling to the three-point PDF, and so on [6].
Instead of introducing this hierarchy of coupled evolution equations for the multipoint PDFs, our strategy [14,15] is to express the unclosed terms as conditional averages, since these are accessible to direct numerical simulations. The result is a partial differential equation governing the joint temperature-velocity PDF and relating its shape as a function of space point r to the conditional averages. The functional form of these conditional averages is a signature of the underlying dynamical processes of RB convection.
Introducing the conditional averages we arrive at the following relation for the single-point probability distribution f (τ, v; r, t): Since at the boundaries of the RB cell the velocity and temperature fields are statistically sharp quantities, the probability distribution has to obey the conditions for arbitrary x, y. We note in passing that a different version of the evolution equation can be obtained by introducing the Laplacian of the PDF. Specializing in the case Pr = 1, it is possible to re-express the conditional averages ∆T |τ, v, r, t and ∆u|τ, v, r, t via the relation Here and in the following, Einstein's summation convention over repeated indices is used. Again, we can introduce conditional expectations, where the arguments of the conditional averages have been abbreviated as = τ, v, r, t. The resulting evolution equation reads A somewhat more complicated equation holds for Pr = 1. This relationship shows that the single-point joint temperature-velocity PDF f (τ, v; r, t) is essentially determined by the conditionally averaged dissipation-like terms (∇T ) 2 | , ∇T · ∇u j | and ∇u i · ∇u j | as well as the conditional pressure gradient −∇p| .

Single-Point Temperature PDF and Implications of Statistical Symmetries
In the following we shall restrict our attention to the reduced temperature probability distribution and its evolution equation. As it will turn out, already this equation gives insight into the connection between RB dynamics and temperature statistics, besides obviously describing the temperature statistics itself. Also, because the final evolution equation of temperature PDF involves scalar valued functions only, a numerical approach is easily feasible.
The reduced temperature PDF is obtained by integrating out the velocity part: Starting from equation (2.6) we obtain the simple equation where we have performed the integration with respect to the velocity v. Thereby, we had to introduce the conditional averages u|τ, r, t and ∆T |τ, r, t . Alternatively, one could re-enact the derivation performed in section 2 for the fine-grained PDF of the temperature h(τ ; r, t) = δ(τ − T (r, t)). Analogous to (2.8), we can derive a further relation using the identity This relation also follows from (2.7). With this equation, we can now summarize the equation for the temperature PDF h(τ ; r, t) in the form Here, the conditional average of the term (∇T ) 2 comes up, which is related to the Nusselt number. Details will be discussed in section 4. The evaluation of the conditional averages appearing in (3.2) is greatly simplified by considering convection that is statistically stationary in time and has periodic horizontal boundaries (i.e. is homogeneous with respect to the horizontal coordinates). Under the assumption of a statistically stationary flow, the PDFs and therefore the conditional averages cannot depend on the time variable; also the dependency on the horizontal coordinates drops out. So instead of dealing with statistical quantities that depend on r and t, we simply have to retain the z-dependence.
Let us first consider the determining equation for the temperature PDF h(τ ; z) in the form obtained from (3.2). This equation in principle has to be solved together with the boundary conditions h(τ ; 0) = δ τ − 1 2 and h(τ ; 1) = δ τ + 1 2 and with appropriately modeled expressions for the conditional averages u z |τ, z and ∆T |τ, z . This direct approach will not be conducted in the present paper though, but can be taken as a starting point for future modeling.
Instead, this first order partial differential equation (PDE) can be analyzed with the help of the method of characteristics [18]. Applying this method, one can find curves τ (s), z(s) in the τ -z-phase space parameterized by s along which the PDE (3.5) transforms into an ordinary differential equation which can be integrated. This approach will be sketched in the following.
Writing h(s) = h τ (s); z(s) and calculating the derivative of it gives The PDE (3.5) is re-expressed in the form Comparing these two equations identifies the characteristic curves as solutions of Along these curves, the PDE (3.5) becomes which can be integrated to This equation describes the evolution of the PDF along a trajectory τ (s), z(s) starting at point τ (s 0 ), z(s 0 ) in phase space. A particularly appealing property of this formalism is that it allows to interpret the statistical results in an illustrative manner, because the characteristics, i.e. trajectories in τ -z-phase space, show the evolution of the "averaged" physical process.
It is tempting to interpret the characteristics as a kind of Lagrangian dynamics of a tracer particle inside the RB cell. However, the dynamic of a tracer particle is stochastic, whereas the characteristics defined by (3.8) describe purely deterministic trajectories and, thus, take the stochastic properties into account only in an averaged way. In a sense, the characteristics describe the averaged evolution of an ensemble of fluid particles that are defined by their initial condition in the τ -z-plane.
Thinking of turbulent RB convection with some physical intuition, one can expect certain features from the statistical quantities introduced in this section. The conditionally averaged vertical velocity u z |τ, z should show positive correlation with the temperature, i.e. it should mirror the well-known fact that hot fluid rises up and cold fluid sinks down. Also the no-slip boundaries should be recognizable for z ≈ 0 and z ≈ 1, respectively. The absolute value of the heat diffusive term ∆T |τ, z should be highest near the boundaries because of the sharp change of the temperature profile. As the characteristics in a way describe the average path a fluid particle takes through τ -zphase space, the typical Rayleigh-Bénard cycle of fluid heating up at the bottom, rising up, cooling down at the top and sinking down again should find its correspondence in the statistical quantities describing the evolution of the PDF. Actual numerical studies of these quantities will be discussed in detail in section 5.
In the 1990s, an approach similar to the presented one was undertaken by Yakhot and followers, where PDFs of various quantities in a stationary flow were expressed as integrals over conditionally averaged variables. Considered quantities included passive scalars by Sinai and Yakhot [19], active scalars, such as temperature fluctuations by Yakhot [20] and temperature increments by Ching [21], and even general functions of an arbitrary quantity measured in the flow by Pope and Ching [22,23]. The conditional averages included spatial and temporal derivatives of these quantities.
In contrast to these works we do not assume homogeneity in the spatial coordinates, i.e. we still have the z-dependency present in our PDF equations. This allows us to discuss the PDFs with respect to the z-coordinate and observe qualitatively different statistics in different regions of the flow. As a result, we are able to see the differing behaviour of bulk and boundary parts of the convection cell and how these are connected to the conditional averages and their dependency on the vertical position. This has not been addressed in the literature.

Connection to the Heat Transport
The Nusselt number Nu as the ratio of convective to conductive heat transfer plays a key role in the analysis of RB convection. It serves as a measure of how efficient heat can be transported through the convection cell.
Though we did not consider the Nusselt number so far, we note that there is an interesting connection between Nu and the conditional averages that appear in our calculations. Because the Nusselt number may be defined as the volume average Nu = (∇T ) 2 V in non-dimensional units, the following expression comes up: Therefore, (∇T ) 2 |τ, r can be viewed as a conditional Nusselt number density. We point out that this quantity should be of considerable interest for the evaluation of theories concerning the Rayleigh number dependency of the Nusselt number based on a decomposition of the heat transport into bulk and boundary contributions, which underlies the Grossmann-Lohse theory [24] outlined in the review of Ahlers et al. [4]. Also, this term is linked with the temperature dissipation rate, which is discussed in [25] with respect to temperature PDFs. In a similar manner, we can employ the temperature-velocity joint PDF f (τ, v; r) to derive an equation relating Rayleigh and Nusselt number. From the relation Ra(Nu −1) = (∇u) 2 V it is straightforward to see that These two exact relations underline the importance of the conditional averages of (∇T ) 2 and (∇u) 2 that naturally come up in our derivations.

Numerical Results
The benefit of our theoretical approach is that we can easily provide it with measurements and data in the form of numerical results. To this end, we solve the basic Oberbeck-Boussinesq equations (1.1) with a standard dealiased pseudospectral code on a three-dimensional equidistant Cartesian grid with periodic boundary conditions. For an introduction into this topic the reader is referred to [26,27]. Periodic boundaries are required in horizontal direction, but in vertical direction Dirichlet conditions for velocity and temperature (i.e. no-slip boundaries of constant temperature) are needed. They are enforced by a volume penalization ansatz [28,29,30]: 1+d]. Inside the fluid domain Ω ⊂ Ω c the unaltered Oberbeck-Boussinesq equations are solved, while in the appended extra regions Ω c \Ω a strong exponential damping (− 1 η u and − 1 η θ, respectively, with η 1) is added to the evolution equations (1.1) of velocity and temperature that damps the fields to zero. By simulating the deviation from the linear temperature profile, θ(r, t) := T (r, t) + z − 1 2 , instead of the temperature itself, the desired boundary conditions read u = 0 and θ = 0 for z = 0 and z = 1. This change of variables T → θ allows us to make use of the volume penalization approach in a straightforward manner.
The reason we choose this numerical scheme instead of often used Chebyshev-based codes, described for example in [26] and used in e.g. [31,32], is because it allows for almost arbitrary shaped boundaries and sidewalls. Although this feature is not used in the present paper due to the required horizontal homogeneity, it even allows to simulate cylindrical vessels on a (numerically cheap) Cartesian grid. A more detailed report on this will be published in the future.
Our theoretical derivation relies on the concept of ensemble averages. Of course, through our numerics we can only access a finite subset of all possible ensemble members. So due to the statistical symmetries and by assuming ergodicity, the ensemble average is substituted for a combined volume-time average in the numerics. Likewise, the volumeaverages · V and · A introduced earlier are actually evaluated as combined volume and time averages. Here, time averaging means averaging over 1250 statistically independent snapshots of the fields. The line plots below only show the parts of the statistical quantities where the statistics converged, i.e. where a significant number of events was obtained.
The simulation was conducted for the parameters Ra = 4.33 × 10 7 and Pr = 1. The computational domain Ω c is resolved with N x × N y × N z = 256 × 256 × 192 gridpoints on an equidistant Cartesian grid, where the fluid domain Ω ⊂ Ω c is represented by 256 × 256 × 128 gridpoints. Thus, the aspect ratio is Γ = 2, with the two horizontal dimensions being identical. The Nusselt number is estimated to Nu = 24.2. We mention that we repeated the analysis presented below with a simulation of aspect ratio Γ = 4 and found basically the same features, due to the fact that for sufficiently large Γ the periodic boundaries do not have a significant influence on the flow.  deviation. Figure 2(b) shows slices in τ -direction indicated by the dashed lines in figure 2(a). In the color plot, one clearly observes the sharp change of the temperature PDF from a δ-function at the boundaries across the boundary layer to a shape exhibiting larger tails in the bulk. In addition to these tails, another feature of the PDF is the hump close to the τ = 0-line. This hump corresponds to the most probable value of the temperature. One expects two different dynamical features to be responsible for this special shape; a tempting explanation would be to attribute the hump to the background temperature field of mean temperature, and account the wings for large |τ | for plumes that carry fluid that is much colder or hotter than the surrounding fluid. An evidence for this is the shape of the PDF close to the bottom or top boundaries (but still outside the boundary layers): At z = 4δ T , though the most probable temperature value is moved slightly towards lower temperatures, the PDF exhibits a large tail at high temperatures. The interpretation is that mostly cold fluid gathers in the lower regions of the bulk, being almost at rest (compare the region of u z |τ, z = 4δ T in figure 3(b) corresponding to the hump), while very hot fluid is a more rare event, because hot fluid is convected away quickly due to plume dynamics. The reason why the hot fluid takes greater temperature values than the cold fluid (in terms of absolute value) is that very cold fluid detaching from the top plate already heats up on its way down. PDFs of the same shape for Rayleigh numbers of the same order are reported in [25], where also the dependence on the vertical coordinate is taken into account. The experimental data in [33,21] shows a more pronounced exponential shape of the temperature PDF, which can be attributed to the difference in the Rayleigh numbers which are several orders of magnitude above ours; the numerical data in [25] suggests that the PDFs become more exponential with increasing Rayleigh number. Figure 3 and 4 exhibit the conditional averages introduced in section 3. One can clearly observe the features that were suggested in the aforementioned section. The conditional vertical velocity u z |τ, z is high (low) for hot (cold) fluid respectively, and the no-slip boundary conditions manifest in the fact that u z |τ, z is close to zero for z ≈ 0 and z ≈ 1. Additionally, one observes a stripe close to the τ = 0-line of almost vanishing vertical velocity which coincides with the reddish core (the hump, i.e. the most probable value) of the temperature PDF in figure 2(a). The interpretation is that fluid that is as hot as the mean temperature is neutrally buoyant and neither moves up nor down. Another striking feature is the sudden increase of the vertical velocity for high τ near the boundary layer, i.e. for z = δ T , which we attribute to rising plumes that detach from the hot bottom plate. Again, it must be stressed that these interpretations hold in an averaged sense. Figure 4 shows that the conditional heat diffusion term ∆T |τ, z is (in terms of absolute value) highest at the boundaries, with the term being positive (negative) at the hot bottom (cold top) plate. On the contrary, in the bulk the absolute value is high (low) for very cold (hot) fluid, i.e. in the wings of the temperature PDF. Additionally, the τ -slice near the boundary in figure 4(b) shows an under-and overshoot. The connection of these unique features to the RB dynamics has yet to be understood. By combining the two aforementioned conditional averages to the vector field (3.8) that defines the characteristics as suggested in section 3, one arrives at the vector field depicted in figure 5 -one of our central results. It is easy to interpret this graph by tracing the vector field; one can qualitatively reconstruct the typical RB cycle of fluid heating up at the bottom, rising up while starting to cool down, cooling down drastically at the top plate, falling down towards the bottom plate while warming up a bit and heating up again at the bottom. It is especially illustrative to see that the main contribution of cooling and heating (i.e. biggest movement in τ -direction of phase space) takes place near the boundaries, highlighting the importance of the boundary layers, while obviously the biggest movement in z-direction occurs in the bulk.
Yet one has to consider, for example, that although hot fluid rises up very quickly (referring to the vectors pointing upwards at the right side in figure 5), this does not contribute much to heat transport because these events occur rarely, as indicated by the temperature PDF shown along with the vector field governing the characteristics.
In figure 6, the conditional heat dissipation rate (∇T ) 2 |τ, z is shown, which can be interpreted as a Nusselt number density according to (4.1). The noteable features are a pronounced minimum near the most probable temperature, again coinciding with the reddish core of the PDF. Also, in the boundary layer this quantity attains huge values (note the logarithmic scaling in the color plot), which highlights the fact that the boundary layer contributes much to the heat transport. A similar shape of a related quantity is reported in [25] (using the deviation from the mean temperature profile instead of the temperature itself), though there the conditional average is taken over the whole fluid volume and is hence lacking the z-dependency.

Summary
In the present work, we have analyzed the single-point temperature PDF on the basis of the Lundgren-Monin-Novikov hierarchy by truncating the hierarchy on the first level via the introduction of conditional averages. We have first derived the evolution equation of the full joint PDF of temperature and velocity. Then we focused on the temperature PDF only, which is the central point of our paper, and obtained an evolution equation for it by reducing the joint PDF equation. We assumed rather weak symmetry conditions of statistical stationarity in time and homogeneity in lateral spatial directions; these conditions should be fulfilled at reasonably high aspect ratios even for closed vessels, i.e. are a good approximation of experimental setups in the bulk of the flow. Under these symmetry considerations, the evolution equation of the temperature PDF becomes fairly simple. The arising conditional averages of temperature diffusion ∆T |τ, z and vertical velocity u z |τ, z are estimated by direct numerical simulations using a suitably designed penalization approach, and features of them are discussed. It shows that expected features such as properties of the temperature and velocity boundary layers, correlation of temperature and velocity and so on are related to the form of these conditional averages that naturally come up in our derivations.
The evolution equation of the temperature PDF is readily treated by the method of characteristics. Due to the applied symmetry conditions, the phase space which describes our system becomes two-dimensional, spanned by temperature τ and vertical coordinate z. Because of this reduced dimensionality of the system, the method of characteristics yields a descriptive view of the RB dynamics, resulting in the vector field describing the evolution in τ -z-phase space. The characteristics, i.e. trajectories in τ -z-phase space, are found to reproduce the typical cycle of a fluid parcel. The regions of the main transport in τ -direction have been identified as the boundary layers, while the major movement in z-direction takes place in the bulk. This highlights the importance of the boundary layers to the heat transport. The relation of heat transport in terms of the Nusselt number to the conditional averages introduced in our derivation is briefly discussed, leading us to the definition of a Nusselt number density. It would be very interesting to obtain the statistical quantities describing the evolution of the PDF directly from experiments, e.g. from measurements of instrumented particles as described in [34]. Future efforts will be to not only use the characteristics as an illustrative way to describe the mean movement in phase space, but actually calculate the PDF of temperature from the integral representation (3.10). Also, modeling of the conditional averages, which are up to now estimated from direct numerical simulations, might be feasible; an intermediate step would be to discuss the quantities not in the turbulent case, but close above the bifurcation from heat conduction to convection, where analytical solutions of temperature and velocity fields are available. Though an easy illustration in the form of trajectories in two-dimensional space will not be achievable in the case of the joint PDF, this approach is nevertheless promising and planned for the future, because we hope that already the form of the conditional averages will give insight into the connection of the statistics to the RB dynamics. An intermediate step would be to concentrate on the joint PDF of temperature and vertical velocity, which should among others relate to the dynamics of plumes.