# [

###### Abstract

Rayleigh–Bénard convection, i.e. the flow of a fluid between two parallel plates that is driven by a temperature gradient, is an idealised setup to study thermal convection. Of special interest are the statistics of the turbulent temperature field, which we are investigating and comparing for three different geometries, namely convection with periodic horizontal boundary conditions in three and two dimensions as well as convection in a cylindrical vessel, in order to work out similarities and differences. To this end, we derive an exact evolution equation for the temperature probability density function (PDF). Unclosed terms are expressed as conditional averages of velocities and heat diffusion, which are estimated from direct numerical simulations. This framework lets us identify the average behaviour of a fluid particle by revealing the mean evolution of fluid of different temperatures in different parts of the convection cell. We connect the statistics to the dynamics of Rayleigh–Bénard convection, giving deeper insights into the temperature statistics and transport mechanisms. We find that the average behaviour is described by closed cycles in phase space that reconstruct the typical Rayleigh–Bénard cycle of fluid heating up at the bottom, rising up to the top plate, cooling down and falling down again. The detailed behaviour shows subtle differences between the three cases.

Turbulent RB convection described by projected dynamics in phase space]Turbulent Rayleigh–Bénard convection described by projected dynamics in phase space
J. Lülff , M. Wilczek, R. J. A. M. Stevens, R. Friedrich, D. Lohse]Johannes Lülff^{†}^{†}thanks: Email address for correspondence: johannes.luelff@uni-muenster.de, Michael Wilczek, Richard J. A. M. Stevens, Rudolf Friedrich^{†}^{†}thanks: deceased, Detlef Lohse

## 1 Introduction

In Rayleigh–Bénard convection, a fluid enclosed between two horizontal plates is heated from below and cooled from above, which induces a flow and thereby enhanced heat transport between the plates. This simple setup is the benchmark system to study thermal convection, which is important in nature and technical applications. Prominent examples include convection in the oceans and the atmosphere or plate tectonics in the mantle of the earth. Depending on the control parameters, the Rayleigh–Bénard system displays a variety of different patterns and flow regimes, ranging from laminar to highly turbulent flows.

Apart from special cases like laminar convection, an analytical solution does not exist as turbulent flows are remarkably hard to handle analytically. Because of their importance, a deeper understanding of turbulent convective flows is still desired, despite the inability to solve the basic equations analytically. To achieve this, many different approaches have been pursued. The heat transport as function of the control parameters is of particular interest and is well described by the Grossmann–Lohse theory (Grossmann & Lohse, 2000, 2001; Ahlers et al., 2009; Grossmann & Lohse, 2011, 2012; Petschel et al., 2013; Stevens et al., 2013). There have also been studies on the turbulence properties of Rayleigh–Bénard convection, by e.g. characterising the statistics of temperature readings of thermal probes in the convection cell (Yakhot, 1989; Ching, 1993; Ching et al., 2004; Shang et al., 2008) or by examining the heat transport mechanisms and the large-scale circulation by Eulerian (Bailon-Cuba et al., 2010; Petschel et al., 2011; van der Poel et al., 2011; Ahlers et al., 2012) or Lagrangian (Gasteuil et al., 2007; Schumacher, 2009) approaches. An overview of recent progress on Rayleigh–Bénard convection can be found in Ahlers et al. (2009); Lohse & Xia (2010); Chillà & Schumacher (2012).

In this paper we will describe Rayleigh–Bénard convection with the full single-point temperature statistics using the temperature probability density function (PDF). This in turn gives us information about the dynamics of the convecting fluid. To this end, we will derive an evolution equation for the temperature PDF, feed in numerical data to complete our ansatz and obtain through the statistics a description of the mean dynamics of fluid particles that travel around in the convection cell. A similar method has been used to describe the statistics in homogeneous isotropic turbulence before, see Wilczek & Friedrich (2009); Wilczek et al. (2011); Friedrich et al. (2012). In this paper, we will generalise and extend the work presented by Lülff et al. (2011), where we first introduced the aforementioned method to Rayleigh–Bénard convection.

We start by deriving the framework in the most general form. Since Rayleigh–Bénard setups usually contain a number of symmetries that can be utilised to simplify the problem, we will apply our framework to three different showcases of three- and two-dimensional convection with homogeneous horizontal directions (i.e., periodic boundaries) and three-dimensional convection in a cylinder. All three cases have different statistical symmetries and show slightly different dynamics. The differences between two- and three-dimensional convection and also between fixed side walls and periodic horizontal boundaries are discussed, for example, in van der Poel et al. (2013, 2014). We will use the PDF methods presented here to further work out similarities and differences between these three cases and give a comprehensive description of the statistics and the dynamics of Rayleigh–Bénard convection.

Since the derivation of our framework utilises the basic equations of Rayleigh–Bénard convection, it can be considered as an ansatz from first principles. The basic equations that govern Rayleigh–Bénard convection are the Oberbeck–Boussinesq equations (Oberbeck, 1879; Boussinesq, 1903) for the velocity and temperature field :

(1a) | ||||

(1b) | ||||

(1c) |

Here, the equations have been non-dimensionalised by the heat diffusion time , the vertical height and the heat difference between the upper and lower plate. This introduces the Rayleigh number and the Prandtl number as the control parameters. The vertical coordinate lies in the range and the temperature takes values . Another control parameter that is often taken into account is the aspect ratio which indicates the lateral over the vertical extent of the system.

The remainder of this paper is structured as follows. In Sec. 2 we will briefly recount our method, i.e. derive an equation for the temperature PDF and connect it to the description of the dynamics of the system. This general theory is then applied to three different Rayleigh–Bénard geometries in Sec. 3, namely three- and two-dimensional convection with homogeneous horizontal directions, and three-dimensional convection in a closed cylindrical vessel with . Section 4 closes the article with an interpretation and discussion of the findings.

## 2 Statistical Description of Heat Transport

Our idea to describe Rayleigh–Bénard convection is to start from the temperature PDF.
Therefore we want to derive an equation that describes the temperature PDF and use it to gain insights into the dynamics of the system.
This ansatz is generally referred to as *PDF methods* (Pope, 1984, 2000) or the *Lundgren–Monin–Novikov hierarchy* (Lundgren, 1967; Monin, 1967; Novikov, 1968).
We now give a short overview of this derivation; a more detailed discussion of the framework can be found in Lülff et al. (2011); Wilczek & Friedrich (2009); Wilczek et al. (2011); Friedrich et al. (2012).
Similar equations have been derived for turbulent reactive flows by Pope (1985).
However, there the unclosed terms are modelled instead of estimated from the numerics as in our case.

### 2.1 PDF Methods

The starting point is the definition of the temperature PDF as an ensemble average,

(2) |

where the PDF describes the probability to find fluid of temperature at position and time . Accordingly, is the sample space variable, while is a realisation of the temperature field. The averaging process can be considered as an ensemble average; later on, ensemble averages are evaluated from the numerics by a suitable volume and time average.

Since the definition in Eq. (2) includes an actual realisation of the temperature field, it is now possible to calculate spatial and temporal derivatives of the PDF, i.e. and . These derivatives contain unclosed terms in the form of conditional averages , where, e.g., the appearing conditionally averaged velocity is a function of the sample space variables , and that tells us what the mean velocity is for fluid of given temperature, position and time.

Putting the aforementioned derivatives together and rearranging them gives the desired evolution equation that describes the temperature PDF:

(3a) | ||||

(3b) |

The left-hand side of (3a) can be seen as the convective derivative of the PDF , while the right-hand side of (3a) contains the conditional average of the convective derivative of the temperature field. Since is a realisation of the temperature field, it obeys the Oberbeck–Boussinesq equations, so in Eq. (3b) we replaced the convective derivative of the temperature field by the right-hand side of the Oberbeck–Boussinesq equation (1c).

Above we obtained an evolution equation that links the shape of the temperature PDF to the conditionally averaged velocity and heat diffusion , which have to be supplied externally; in our case, we estimate them from simulations later on.

### 2.2 Method of Characteristics

The above evolution equation (3) that determines the temperature PDF is a first-order partial differential equation. That means that we can apply the method of characteristics (Courant & Hilbert, 1962; Sarra, 2003) which lets us identify the average behaviour of fluid as it travels through phase space.

In a nutshell, by applying the method of characteristics to the evolution equation, one can identify trajectories (the so-called *characteristic curves* or just *characteristics*) in phase space, along which the partial differential equation for the temperature PDF transforms into a set ordinary differential equations for and .
The phase space is spanned by the variables that the temperature PDF depends upon, i.e. , and .
The characteristics are defined by the conditional averages,

(4) |

This states that the characteristics are solutions of Eq. (4) that follow the vector field on the right-hand side of the above equation; the vector field is regarded as the phase space velocity.
From the last line of Eq. (4), , it becomes clear that the parametrisation of the characteristics in phase space, i.e. the arc length, is the same as the time of the system – a *fast* movement in phase space therefore really has to be seen in the temporal sense.

It is now important to notice that, since the characteristics are governed by the conditionally averaged vector field, they show the average behaviour of a fluid parcel in phase space.
In other words, the characteristics can be seen as the mean evolution of an ensemble of Lagrangian particles that share the same coordinates in phase space.
This is what Pope (1985, Sec. 4.5) refers to as *conditional particles* – quasi-particles following the conditionally averaged vector field (4) that show the mean Lagrangian evolution and that have the same statistics as Lagrangian particles.
By examining the conditionally averaged vector field and the resulting characteristic curves, one can investigate the mean transport properties of fluid through phase space and gain insight into the mean heat transport properties of Rayleigh–Bénard convection.
Since the characteristics are trajectories in phase space, the framework can be seen as a quasi-Lagrangian description, but it has to be stressed that it is achieved by utilising the statistics of Eulerian fields alone.

Along the characteristics, the partial differential equation (3) becomes an ordinary differential equation which can be integrated. Thus, the temperature PDF along a certain characteristic evolves according to

(5) |

Here, the integral is a line integral along a characteristic from to . The integral kernel is the phase space divergence evaluated at the phase space position given by the characteristic for time , i.e. . This equation tells us that the temperature PDF along the characteristic that connects the initial point with the point in phase space changes according to the integrated phase space divergence. As an alternative interpretation, Eq. (2.2) determines how the temperature PDF for is traced back to the PDF for . While we will not further investigate Eq. (2.2) in the numerical results of the next section, we included it for the sake of completeness.

Up to now, we have kept the description as general as possible. But usually, a Rayleigh–Bénard setup has a number of statistical symmetries, which simplify the problem, i.e. the phase space dimension is reduced and the estimation of the unknown conditional averages from numerical simulation is simplified. In the next section, we will apply the framework that has been outlined in this section to three different Rayleigh–Bénard geometries with different symmetries and discuss the findings.

## 3 Result from DNS

In this section, we focus on three Rayleigh–Bénard geometries, i.e. three-dimensional convection with periodic horizontal boundaries (Sec. 3.1), two-dimensional convection with periodic horizontal boundaries (Sec. 3.2), and three-dimensional convection in a cylindrical vessel with (Sec. 3.3).

### 3.1 Three-dimensional Convection with Periodic Horizontal Boundaries

First, we consider three-dimensional convection with periodic horizontal boundaries in the statistically stationary state. A snapshot of the instantaneous temperature field taken from the numerics can be seen in Fig. 1. The parameters of the simulation are , , and the aspect ratio of the periodic box is . The two horizontal plates have a constant temperature and a no-slip velocity boundary condition. The numerical setup is a tri-periodic pseudospectral direct numerical simulation, where the boundary conditions are enforced by volume penalization methods (Lülff et al., 2011; Angot et al., 1999; Schneider, 2005; Keetels et al., 2007). The equidistant resolution in , and direction is grid points, and in the vertical direction grid points fall into the boundary layers, according to the criteria given by Shishkina et al. (2010). We calculated the statistical quantities from an ensemble consisting of snapshots, and the snapshots where taken free-fall time units apart.

Statistically stationary Rayleigh–Bénard convection in this geometry is homogeneous in horizontal directions. This means that the statistical quantities only depend on the temperature and the vertical coordinate and not on the horizontal coordinates and or time . Thus, the temperature PDF and the conditional averages read and , and the phase space becomes two-dimensional.

In Fig. 2, the height-resolved mean, standard deviation, skewness and kurtosis of the temperature field are shown. As is well known, the mean temperature is almost constant in the bulk and has a steep gradient towards the hot and cold boundaries at and . The standard deviation takes its highest values close to the boundaries and decreases towards the centre of the bulk, indicating a temperature PDF that is broadening towards the boundaries. The height-resolved skewness takes its highest absolute values near the boundaries and decreases linearly as function of height in the bulk. This can be interpreted as hot fluid that is beginning to cool down on its way from the lower to the upper plate (and vice versa). The kurtosis indicates that, apart from the boundaries, the temperature PDF is more peaked and shows stronger tails than the Gaussian distribution.

When the simplifications resulting from the statistical symmetries are incorporated into the general framework presented in Sec. 2, Eq. (3) that defines the PDF becomes

(6) |

while the vector field (4) of the characteristics reads

(7) |

The PDF and the characteristic curves are therefore defined by the conditional averages of vertical velocity and heat diffusion. The next step is to estimate the conditional averages from the numerics while taking the statistical symmetries into account. Subsequently, the characteristics are obtained by integrating Eq. (7) for arbitrary initial conditions in phase space. Obviously only initial conditions where the PDF and the conditional averages are defined, i.e. where there have been any events at all, can be considered.

When integrating the characteristics for many starting positions, we observed that they tend to converge to what at first seems to be similar to a limit cycle. This cycle is shown in Figs. 5 and 6; later on, we will come back to a description of the dynamics and behaviour of Rayleigh–Bénard convection that these figures offer. In contrast to a limit cycle, though, one would expect that the characteristics form concentric closed curves: To see this, let the phase space be densely seeded by the conditional particles described in Sec. 2.2. As the density of the conditional particles following the characteristics is proportional to the temperature PDF, and a limit cycle acts as an attractor for the conditional particles, the temperature PDF should converge towards a -function that is non-vanishing on the cycle and zero everywhere else. This in turn stands in contrast to the fact that we are considering statistically stationary systems and that the temperature PDF is clearly not a -function, and thus it follows that the characteristics cannot converge to a limit cycle but must form concentric closed curves.

We find that the observed convergence is caused by the flawed conditionally averaged vector field estimated from the numerics by a binning process. The noise inherent to the binning violates the solenoidality of the probability flux (phase space velocity times PDF) as demanded by Eq. (6), and thus the imperfect binned vector field contains many localised sinks where the characteristics converge.

By smoothing the binned data through a convolution with a Gaussian kernel and projecting the vector field onto the solenoidal part (i.e., enforcing the validity of Eq. (6)), we are indeed able to find the expected concentric closed curves, as exemplified in Fig. 3.

Here the horizontal axis corresponds to the temperature coordinate and the vertical axis to the vertical coordinate of the phase space; the background colour coding gives the temperature PDF . For every starting point located on the -axis, the characteristics perform a closed loop in counter-clockwise direction that shows how particles on average evolve through phase space. By tracing its course, one is able to reconstruct the typical Rayleigh–Bénard cycle a conditional particle undergoes: A fluid element near the lower plate first heats up and then starts to move up towards the cold plate. During its upward travel it slowly cools down and then becomes much colder when it is close to the top plate before it falls down again towards the lower plate while beginning to heat up and starting what we call the RB cycle all over again. The cycles for fluid starting at more moderate temperatures (i.e., near ) show a smaller amplitude in both - and -direction. We do not find closed circles located further outwards, because the characteristics would visit areas of the phase space were no events are recorded in the numerics, and thus the vector field is undefined there. This usually happens near the vertical boundaries of the phase space (i.e., near the plates) where the support of the vector field becomes very narrow. Furthermore, we remark that the precise appearance of the closed cycles slightly varies with the parameters of the post-processing described above (e.g., the amount of smoothing); nevertheless, the description of the qualitative behaviour we are aiming at is found to be robust.

Without the aforementioned projection to keep the vector field solenoidal in the presence of noise, the characteristics converge in some parts of the phase space. To study this behaviour, and also the near-wall regions, we seed characteristics on the -line and integrate them backwards in time, as shown in Fig. 4.

The comparison of the characteristics in the regions and immediately reveals the convergence for the non-solenoidal case, because the density of the characteristics on the right side (i.e., at later times) is higher than on the left side. Furthermore, the characteristics on the far right enter the boundary at when integrated backwards in time, or, to put it in other words, characteristics leave the boundaries very close to each other and then become less dense when following them forward in time. This corresponds to the fact that at , the temperature PDF becomes a -function located at due to the Dirichlet boundary condition. In fact, the characteristics should end precisely in this point, but this behaviour is not resolved here. Likewise, by symmetry considerations, the characteristics approaching the boundary from the far left side should also enter the point , which we do not observe. We speculate that a numerical resolution that goes beyond what is demanded by the criterion of Shishkina et al. (2010) may be needed to capture the correct behaviour of the characteristics in this singular region of phase space.

The cycle to which the characteristics converge in the presence of noise is shown in Fig. 5. This cycle is very similar to the closed curves shown in Fig. 3; in fact, it is almost completely embedded between two adjacent closed curves. While the solenoidal projection helps to find the required closed curves, the conditionally averaged vector field still containing the noise gives the same general picture of the RB cycle. Additionally, while removing divergences, the projection may introduce unpredictable systematic errors as a side effect, especially in regions of the phase space where the support of the conditional averages and the PDF changes rapidly. Therefore, as we want to focus on the qualitative features of Rayleigh–Bénard convection that the RB cycle as well as the vector field represent, we will from now on only consider the cases without the solenoidal projection and use the found cycle as one generic representative of the family of closed concentric curves. The same applies to the convection cases discussed in Secs. 3.2 and 3.3 where the characteristics also tend to converge due to imperfections induced by noise.

We now come back to the discussion of the qualitative features that the conditionally averaged phase space velocity and the RB cycle describe.
In Fig. 5, the RB cycle is shown together with the temperature PDF and in Fig. 6 together with the phase space velocity.
From the first-mentioned figure, it is interesting to see that hot fluid on the RB cycle has the highest phase space speed in the range , i.e. hot rising fluid is fastest in the lower half of the convection vessel, and likewise for cold fluid in the upper half due to the up-down symmetry of Rayleigh–Bénard convection.
As a side note, from now on, whenever we describe a fluid process, the *reversed* process – interchanging hotcold, bottomtop, updown etc. – is also implied.

The temperature PDF in Fig. 5 shows that the temperature distribution changes with the vertical coordinate and contracts to a -function at the fixed temperature boundaries. Furthermore, one can map the shape of the distribution to the higher moments in Fig. 2, i.e. the peaks of the standard deviation near the boundaries and the linear dependence of the skewness on the height in the bulk region. We note that the isocontours of the PDF do not lie tangent to the vector field or the RB cycle because the divergence of the phase space velocity, , is non-vanishing.

In Fig. 6, the black and white background colour corresponds to the phase space speed, and the temperature is given by the colour of the arrows.
This display lets us identify how fluid behaves in different parts of the phase space.
Near the boundaries, fluid of all temperatures displays high phase space speeds, while in the bulk only fluid of *intense* temperature (i.e., deviating strongly from the mean) has high speeds.
This supports the previous finding that hot fluid on the RB cycle has its highest speed in the lower half.
Fluid that has the mean temperature (cf. upper-left panel of Fig. 2) is found to be at rest because no buoyancy acts on it and it does not heat up or cool down.
As can be expected *a priori*, the vector field shows that the main movement in -direction, i.e. heating and cooling, takes place near the boundaries, while the main movement in vertical direction happens in the bulk.

### 3.2 Two-dimensional Convection with Periodic Horizontal Boundaries

The next case to investigate is two-dimensional Rayleigh–Bénard convection with periodic horizontal boundaries. The parameters are , and , and the numerical scheme that is used is identical to the one from Sec. 3.1. Again, the two horizontal plates are no-slip walls of fixed temperature. The numerical resolution is equidistant grid points with grid points in the vertical direction falling into the boundary layers (cf. Shishkina et al. (2010)), and the ensemble consists of snapshots separated by free-fall time units.

A snapshot of the temperature field is shown in Fig. 7, and one can see coherent structures in the form of four plume hot spots (two at the top, two at the bottom) and four convection rolls, even at this intermediate Rayleigh number. Also, localised round blobs of hot and cold fluid can be found. The statistical symmetries in this system are identical to the ones for the three-dimensional periodic case discussed before, which means the phase space becomes two-dimensional and the statistics depend on and only.

Figure 8 shows the first four height-resolved standardised moments of temperature. While the mean temperature profile has the same shape as the one from three-dimensional convection (cf. Fig. 2), the higher moments show subtle new features. For the three-dimensional case the moment profiles as function of the height are smoother than for the two-dimensional case. Especially the skewness shows transitions and is in the bulk not as linear as for the three-dimensional case. We link this to the coherent structures in the form of plume hot spots in the two-dimensional case because the position of the transitions in the moments corresponds to the vertical size of the plume hot spots (cf. Fig. 7). In the plume hot spots, there is a re-cycling of hot fluid, which means that fluid that is hotter than the mean temperature profile is trapped near the hot bottom plate for some time instead of being advected upwards directly, cf. Sugiyama et al. (2010). This hot trapped fluid causes a temperature distribution that is near the lower plate strongly skewed towards higher temperatures. Above the hot spot, only a sharp jet of hot fluid remains which results in a flatter profile of the skewness. In three-dimensional convection the trapping mechanism of plume hot spots is missing: Loosely speaking, in three-dimensional convection, there is another lateral dimension into which the fluid can escape and be advected away, forming the sheet-like plumes that can also be found in Fig. 1 (Schmalzl et al. (2004) and van der Poel et al. (2013) also discuss differences in flow structures between two- and three-dimensional convection). Therefore, the entrapment seen in two-dimensional convection does not occur in three dimensions, which means that a strong mechanism that in two dimensions traps hot fluid near the bottom is missing for the three-dimensional case.

When we estimate the conditional averages for two-dimensional convection from the simulation and then numerically calculate the characteristics, we again find that the dynamics in phase space resemble a closed cycle. This RB cycle is shown in Figs. 9 and 10. While the RB cycle displays the same basic cycle of fluid heating up at the bottom, moving upwards, cooling down at the top plate and falling down again, there are some differences as compared to the three-dimensional case (Figs. 5 and 6).

The most striking new feature is a kink in the RB cycle in the lower left and upper right corners. As can be seen from the colour coding of the cycle, this is also the region of the lowest phase space speed. Therefore, we also link this region to the re-cycling areas discussed above because fluid trapped in the plume hot spots undergoes almost no net vertical movement and needs more time to heat up in comparison to fluid that is in direct contact with the much hotter bottom plate. A similar argument is discussed by van der Poel et al. (2013).

Another difference with the three-dimensional case is that there is a bulge in the temperature PDF towards higher temperatures (around , ). This bulge is due to hotter-than-average fluid that gathers near the bottom plate and is therefore compatible with the interpretation of the re-cycling fluid from above; also, this bulge gives a direct impression of the high skewness values found in this region (cf. Fig. 8).

Figure 10 shows the vector field of the phase space velocities together with its norm (coded in black and white). The phase space velocities are more heterogeneously distributed as compared to the three-dimensional case, e.g. the high speeds in the bulk for intense temperatures are more pronounced (cf. Fig. 5). These strong vertical movements in the bulk lead to higher phase space speeds as compared to three dimensions (see the colour scale in Figs. 6 and 10). We think this can only in part be attributed to the difference in Rayleigh numbers ( vs. ), but is also due to the coherent structures found in two dimensions, i.e. plume hot spots as localised events of intense temperature. It is also found that in the kink region the RB cycle passes through the region of lowest phase space speed, which can be understood as the average dynamics being slowed down in the re-circulating plume hot spots.

### 3.3 Three-dimensional Convection in Cylindrical Vessel

The last Rayleigh–Bénard geometry under investigation is a closed cylindrical vessel. The control parameters are , and (diameter over height). All the walls are no-slip, and the horizontal plates are of constant temperature while the side walls are thermally insulating. The ensemble consists of snapshots that are obtained from direct numerical simulation using a second-order finite difference scheme on a staggered cylindrical grid (Verzicco & Camussi, 2003) with a resolution of grid points (with , and being the azimuthal, radial and vertical coordinate, respectively). The boundary layer contains grid points in the vertical direction. The snapshots are separated by free-fall time unit. Figure 11 shows a snapshot of the temperature field.

Rayleigh–Bénard convection in a cylinder has statistical symmetries that are different from the former two cases of horizontally homogeneous convection. In addition to the temperature and vertical position the statistics now also depend on the radial position in the cylinder. Here corresponds to the cylinder axis and the sidewall is at .

In Fig. 12 the --resolved first four standardised moments of the temperature distribution are shown. The horizontal and vertical axes correspond to the radial coordinate and the vertical coordinate , respectively. The mean temperature profile is almost constant in the bulk, and only near the hot and cold plate (and to a lesser extent near the sidewalls) a deviation can be seen. Like in the former cases, the standard deviation of the temperature takes its highest values near the horizontal plates and falls off towards the middle of the convection cell with a local minimum at . Also, it can be seen that this local minimum is less pronounced near the side walls, i.e. for high . After rising up to its maximum value at , the skewness varies monotonically with increasing . This has also been found for the former two cases. Regarding the radial dependence, the skewness falls off towards the sidewalls, indicating a less asymmetric temperature distribution there, while its highest values are found near the cylinder axis. Although the statistics are less converged for the kurtosis (especially near due to the cylindrical geometry), one can see that the highest values correspond roughly to the extrema of the skewness. These high values of skewness and kurtosis near the bottom wall can be attributed to hot localised plumes detaching from the hot bottom plate and piercing into the colder fluid of the bulk. We also note that the absolute values of skewness and kurtosis are higher than in the former two cases (Figs. 2 and 8), which can be understood on dimensional grounds: In horizontally periodic convection, we averaged over all horizontal directions and therefore averaged out the sharp maxima that can be seen in the cylindrical case (Fig. 12) where the statistics are resolved additionally in the horizontal coordinate .

When inserting the temperature PDF and the conditional averages into the general framework from Sec. 2, the PDF-defining equation (3) becomes

(8) |

while the characteristics (4) read

(9) |

In comparison to the former two cases, we now deal with a three-dimensional phase space where the additional dimension is related to the radial movement. From Eq. (9) one sees that the radial coordinate of the characteristics evolves according to the conditional average of radial velocity .

We now again turn to the integration of the characteristics, following Eq. (9). Although cylindrical convection is intrinsically different from the former two cases of horizontally periodic convection (due to three- vs. two-dimensional phase space), we still find that the average dynamics of fluid parcels are described by a closed, twisted loop in phase space that shares common features with the former two. The cycle is shown in Figs. 13 and 14 as the slender figure-eight-shaped curve.

In the left panel of Fig. 13, the background shows the mean temperature (cf. Fig. 12) and the figure-eight-shaped curve shows a projection of the RB cycle into the --plane. The third coordinate of the RB cycle, the temperature , is colour coded. The temperature scale corresponds to the minimal and maximal temperature ( and ) the RB cycle takes. When tracing the cycle, one can again identify the Rayleigh–Bénard cycle of the horizontally periodic convection cases, superposed with an additional inwards and outwards motion: Starting with fluid of mean temperature that is quickly heating up at the bottom, it then begins to rise up and move inwards into the bulk until and , where it goes outwards and starts to cool down. At the maximal , the fluid cools down quickly and then falls towards the lower plate while moving inwards, thus starting the RB cycle all over again. Additionally, one can see that the hot fluid rising from the lower plate steadily cools down when it crosses the bulk of almost uniform temperature; this is related to the monotonically decreasing skewness of temperature that can be seen in Fig. 12.

The difference of the temperature of the RB cycle and the background temperature (cf. colour coding of these two in Fig. 13 (left)) shows that the regions where the temperature of the RB cycle deviates most from the mean background temperature are the regions of high buoyancy and correspond to the regions of main vertical movement in the bulk. To elaborate on this, the right panel of Fig. 13 shows the standard deviation of the temperature field in the background, and the colour coding of the cycle shows the absolute deviation of its temperature coordinate from the mean temperature. The deviation of the temperature of a fluid particle on the cycle from the surrounding mean temperature determines its mean buoyancy, so the right panel tells us how strong the buoyancy acts; the highest values for hot rising fluid are found in the lower half of the convection cell. In comparison, the mean deviation of fluid from the mean temperature profile (shown in the background as the standard deviation of temperature) is much weaker. To summarise, from the left panel of Fig. 13 one can see in which direction the buoyancy acts, while the right panel shows its strength.

The vector field of the characteristics is shown in Fig. 14. Due to the difficulty to display a vector field in three-dimensional phase space, we show slices of the phase space velocity in the - plane at different . The panels I–VIII show one cycle of a fluid parcel travelling along the RB cycle, with its temperature colour coded as in Fig. 13 (left). Additionally to the RB cycle described above, Fig. 14 also reveals the average behaviour of fluid in different parts of the convection cell, conditioned on its temperature. The arrows show an - slice of the vector field of the characteristics (9) at the coordinate of the cycle. Also, the black and white background colour indicates the phase space speed, i.e. how fast a fluid parcel travels through phase space (with white being the fastest movement). The arrows indicate the mean movement of fluid of a given temperature in different regions of the convection cell.

Panel I shows that cold fluid (here, ) has the highest speeds in the bulk and near the sidewall. Near the cold top plate, cold fluid is mainly transported towards the outer wall. For , the direction of movement is slightly tilted towards the cylinder axis. Cold fluid that falls down along the sidewall is deflected towards the inner cylinder at around due to a corner flow that propels cold fluid upwards along the side wall; notice that at and , up- and downwelling cold fluid collides. This feature can be understood as cold plumes that are formed at the upper plate and are swept towards the sidewalls. The plumes then fall down along the sidewall until they hit the hot fluid at the bottom plate where they are directed inwards. These cold plumes that fall down along the sidewall can also be seen in Fig. 11.

The fluid from panel III is less cold () and has overall lower speeds, but still shows the same features as in panel I, e.g. cold fluid is swept along the upper plate outwards and falls down along the side wall until it hits the upwelling corner flow. The vector fields for fluid of mean temperature (, panel IV and VIII) are symmetric around and show an almost vanishing velocity in the bulk. Near the horizontal plates, fluid is swept towards the sidewalls and from there vertically towards the -line.

The rest of the panels complete one run of the RB cycle and due to the up-down symmetry contain the same information already described. Still, from all eight panels it can be seen that fluid for all temperatures has high phase space speed near the sidewalls, which can be attributed to plumes that are guided along the outer walls of the cylinder, and also has high speeds near the bottom and top plate which is due to the vigorous temperature contrast between fluid and horizontal plates that leads to high speeds regarding the coordinate. Also, we have to stress here that the corner flows show up due to the investigation being conditioned on the temperature and not due to some large-scale structure that may be present in the fields; this structure is lost in the azimuthal averaging process when obtaining the vector field from Eq. (9). Therefore, in our analysis the corner flows are statistical structures that are not necessarily related to structures in the flow.

## 4 Summary

In this paper, we analysed the turbulent flow in Rayleigh–Bénard convection on the basis of statistical quantities like the temperature PDF and conditionally averaged fields.
We derived that the mean path a fluid particle takes through phase space (spanned by temperature and spatial coordinates) is defined by so-called *characteristics*, i.e. trajectories in phase space that follow the conditionally averaged vector field composed of heat diffusion and fluid velocities.
Thereby, we could characterise the dynamics and flow features that occur in turbulent convection cells from a statistical point of view, i.e. from averaged quantities like the temperature distribution and its moments as well as statistics conditioned on temperature and spatial position.

By estimating the aforementioned vector fields for three different Rayleigh–Bénard geometries while utilising their symmetries and then integrating the characteristics we described the mean dynamics that fluid particles undergo, i.e. we could describe how fluid of different temperatures behaves in different regions of the convection volume. We also distinguished regions of high and low transport through phase space. For all geometries there are high phase space speeds for intense temperatures in the bulk (which we attribute to localised events of intense temperatures and high speeds, i.e. plumes) as well as high speeds near the horizontal plates for all temperatures, while for the case of cylindrical convection the phase space speed also takes high values near the wall of the cylinder. This we interpret as plumes that are directed along the insulating sidewalls. In the conditionally averaged vector field of the cylinder, we could furthermore identify corner flows near the sidewalls for fluid of different temperatures. Cold fluid experiences a corner flow near the bottom plate while showing no corner-flow near the upper plate and vice-versa. Additionally, we described the higher moments of the temperature distributions, where we could link features of the moments to coherent structures that appear in turbulent flows.

When we then obtained the characteristics by integrating trajectories through the conditionally averaged vector field, we found that for all different convection setups, the characteristics form closed cycles in phase space. These cycles reconstruct the typical Rayleigh–Bénard cycle a fluid particle undergoes on average, i.e. fluid is heated up at the bottom and rises upwards while slightly cooling down until it hits the upper plate, where it cools down fast and falls down to the lower plate while slightly heating up, thus starting the cycle all over again. In the cylindrical case, where there is another phase space dimension corresponding to horizontal movement, the fluid shows an additional inwards and outwards movement while following the RB cycle. The method thus allows to further pin-point and quantify the differences and similarities between Rayleigh–Bénard convection in two- and three-dimensional periodic boxes and in three-dimensional convection in a cylindrical cell.

For so-called homogeneous Rayleigh–Bénard convection (Lohse & Toschi, 2003; Calzavarini et al., 2005) – thermal convection with horizontally periodic boundary conditions as well as periodic boundary conditions in vertical direction together with an imposed temperature gradient driving the flow – we would expect quite different behaviour as such a system does not have boundary layers, but represents pure bulk turbulence. Furthermore, the statistics do not depend on the spatial coordinates and the phase space becomes one-dimensional, which is fundamentally different from the three cases considered in the present paper. The analysis of homogeneous Rayleigh–Bénard convection would actually be more in line with the work by Yakhot (1989) and Ching (1993), who analyse with the help of conditional averages the PDF of experimental time series obtained from temperature probe measurements. In these works, in a sense the phase space is also one-dimensional due to the lack of any spatial coordinate. Thus, while their ansatz can be used to e.g. describe the transition from soft to hard turbulence, no information about the dynamics and spatial structures can be obtained from the statistics, contrary to the method we proposed here.

The authors thank Prof. Stephen B. Pope for his insightful remarks regarding the interpretation of our results and the referees for their constructive input.
JL acknowledges fruitful discussions with and valuable advice by Oliver Kamps, Anton Daitche and Theodore Drivas as well as funding by the DFG (Deutsche Forschungsgemeinschaft) under grant FR 1003/10-1.
MW acknowledges support from the DFG under project WI 3544/2-1 and WI 3544/3-1.
RJAMS was supported by the *Fellowships for Young Energy Scientists* (YES!) of FOM.
Parts of the computations have been performed on PROVIDE and Huygens (SURFsara).

## References

- Ahlers et al. (2009) Ahlers, G., Grossmann, S. & Lohse, D. 2009 Heat transfer and large scale dynamics in turbulent Rayleigh–Bénard convection. Rev. Mod. Phys. 81 (2), 503–537.
- Ahlers et al. (2012) Ahlers, G., He, X., Funfschilling, D. & Bodenschatz, E. 2012 Heat transport by turbulent Rayleigh–Bénard convection for and : aspect ratio . New J. Phys. 14 (10), 103012.
- Angot et al. (1999) Angot, P., Bruneau, C.-H. & Fabrie, P. 1999 A penalization method to take into account obstacles in incompressible viscous flows. Numer. Math. 81 (4), 497–520.
- Bailon-Cuba et al. (2010) Bailon-Cuba, J., Emran, M. S. & Schumacher, J. 2010 Aspect ratio dependence of heat transfer and large-scale flow in turbulent convection. J. Fluid Mech. 655, 152–173.
- Boussinesq (1903) Boussinesq, J. 1903 Théorie Analytique de la Chaleur. Paris: Gauthier-Villars.
- Calzavarini et al. (2005) Calzavarini, E., Lohse, D., Toschi, F. & Tripiccione, R. 2005 Rayleigh and Prandtl number scaling in the bulk of Rayleigh–Bénard turbulence. Phys. Fluids 17 (5), 055107.
- Chillà & Schumacher (2012) Chillà, F. & Schumacher, J. 2012 New perspectives in turbulent Rayleigh–Bénard convection. Eur. J. Phys. E 35, 1–25.
- Ching (1993) Ching, E. S. C. 1993 Probability densities of turbulent temperature fluctuations. Phys. Rev. Lett. 70 (3), 283–286.
- Ching et al. (2004) Ching, E. S. C., Guo, H., Shang, X.-D., Tong, P. & Xia, K.-Q. 2004 Extraction of plumes in turbulent thermal convection. Phys. Rev. Lett. 93, 124501.
- Courant & Hilbert (1962) Courant, R. & Hilbert, D. 1962 Methods of Mathematical Physics Volume II. Wiley-Interscience.
- Friedrich et al. (2012) Friedrich, R., Daitche, A., Kamps, O., Lülff, J., Voßkuhle, M. & Wilczek, M. 2012 The Lundgren–Monin–Novikov hierarchy: Kinetic equations for turbulence. Compt. Rend. Phys. 13 (9–10), 929–953.
- Gasteuil et al. (2007) Gasteuil, Y., Shew, W. L., Gibert, M., Chillà, F., Castaing, B. & Pinton, J.-F. 2007 Lagrangian temperature, velocity, and local heat flux measurement in Rayleigh–Bénard convection. Phys. Rev. Lett. 99 (23), 234302.
- Grossmann & Lohse (2000) Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: a unifying theory. J. Fluid Mech. 407, 27–56.
- Grossmann & Lohse (2001) Grossmann, S. & Lohse, D. 2001 Thermal convection for large Prandtl numbers. Phys. Rev. Lett. 86 (15), 3316–3319.
- Grossmann & Lohse (2011) Grossmann, S. & Lohse, D. 2011 Multiple scaling in the ultimate regime of thermal convection. Phys. Fluids 23 (4), 045108.
- Grossmann & Lohse (2012) Grossmann, S. & Lohse, D. 2012 Logarithmic temperature profiles in the ultimate regime of thermal convection. Phys. Fluids 24 (12), 125103.
- Keetels et al. (2007) Keetels, G. H., D’Ortona, U., Kramer, W., Clercx, H. J. H., Schneider, K. & van Heijst, G. J. F. 2007 Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume-penalization: Convergence of a dipole-wall collision. J. Comput. Phys. 227 (2), 919–945.
- Lohse & Toschi (2003) Lohse, D. & Toschi, F. 2003 Ultimate state of thermal convection. Phys. Rev. Lett. 90, 034502.
- Lohse & Xia (2010) Lohse, D. & Xia, K.-Q. 2010 Small-scale properties of turbulent Rayleigh–Bénard convection. Annu. Rev. Fluid Mech. 42 (1), 335–364.
- Lülff et al. (2011) Lülff, J., Wilczek, M. & Friedrich, R. 2011 Temperature statistics in turbulent Rayleigh–Bénard convection. New J. Phys. 13 (1), 015002.
- Lundgren (1967) Lundgren, T. S. 1967 Distribution functions in the statistical theory of turbulence. Phys. Fluids 10 (5), 969–975.
- Monin (1967) Monin, A. S. 1967 Equations of turbulent motion. Prikl. Mat. Mekh. 31 (6), 1057–1068.
- Novikov (1968) Novikov, E. A. 1968 Kinetic equations for a vortex field. Sov. Phys.-Dokl. 12 (11), 1006–1008.
- Oberbeck (1879) Oberbeck, A. 1879 Über die Wärmeleitung der Flüssigkeiten bei Berücksichtigung der Strömungen infolge von Temperaturdifferenzen. Ann. Phys. Chem. 7, 271–292.
- Petschel et al. (2013) Petschel, K., Stellmach, S., Wilczek, M., Lülff, J. & Hansen, U. 2013 Dissipation layers in Rayleigh–Bénard convection: A unifying view. Phys. Rev. Lett. 110, 114502.
- Petschel et al. (2011) Petschel, K., Wilczek, M., Breuer, M., Friedrich, R. & Hansen, U. 2011 Statistical analysis of global wind dynamics in vigorous Rayleigh–Bénard convection. Phys. Rev. E 84 (2), 026309.
- van der Poel et al. (2014) van der Poel, E. P., Ostilla-Mónico, R., Verzicco, R. & Lohse, D. 2014 Effect of velocity boundary conditions on the heat transfer and flow topology in two-dimensional Rayleigh–Bénard convection. Phys. Rev. E 90, 013017.
- van der Poel et al. (2011) van der Poel, E. P., Stevens, R. J. A. M. & Lohse, D. 2011 Connecting flow structures and heat flux in turbulent Rayleigh–Bénard convection. Phys. Rev. E 84 (4), 045303.
- van der Poel et al. (2013) van der Poel, E. P., Stevens, R. J. A. M. & Lohse, D. 2013 Comparison between two- and three-dimensional Rayleigh–Bénard convection. J. Fluid Mech. 736, 177–194.
- Pope (1984) Pope, S. B. 1984 Calculations of a plane turbulent jet. AIAA Journal 22 (7), 896–904.
- Pope (1985) Pope, S. B. 1985 PDF methods for turbulent reactive flows. Prog. Energy Combust. Sci. 11, 119–192.
- Pope (2000) Pope, S. B. 2000 Turbulent Flows. Cambridge, England: Cambridge University Press.
- Sarra (2003) Sarra, S. A. 2003 The method of characteristics with applications to conservation laws. J. Onl. Math. Appl. 3, 1–6.
- Schmalzl et al. (2004) Schmalzl, J., Breuer, M., Wessling, S. & Hansen, U. 2004 On the validity of two-dimensional numerical approaches to time-dependent thermal convection. Europhys. Lett. 67 (3), 390–396.
- Schneider (2005) Schneider, K. 2005 Numerical simulation of the transient flow behaviour in chemical reactors using a penalisation method. Comput. Fluids 34 (10), 1223–1238.
- Schumacher (2009) Schumacher, J. 2009 Lagrangian studies in convective turbulence. Phys. Rev. E 79 (5), 056301.
- Shang et al. (2008) Shang, X.-D., Tong, P. & Xia, K.-Q. 2008 Scaling of the local convective heat flux in turbulent Rayleigh–Bénard convection. Phys. Rev. Lett. 100, 244503.
- Shishkina et al. (2010) Shishkina, O., Stevens, R. J. A. M., Grossmann, S. & Lohse, D. 2010 Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution. New J. Phys. 12 (7), 075022.
- Stevens et al. (2013) Stevens, R. J. A. M., van der Poel, E. P., Grossmann, S. & Lohse, D. 2013 The unifying theory of scaling in thermal convection: The updated prefactors. J. Fluid Mech. 730, 295–308.
- Sugiyama et al. (2010) Sugiyama, K., Ni, R., Stevens, R. J. A. M., Chan, T. S., Zhou, S.-Q., Xi, H.-D., Sun, C., Grossmann, S., Xia, K.-Q. & Lohse, D. 2010 Flow reversals in thermally driven turbulence. Phys. Rev. Lett. 105 (3), 034503.
- Verzicco & Camussi (2003) Verzicco, R. & Camussi, R. 2003 Numerical experiments on strongly turbulent thermal convection in a slender cylindrical cell. J. Fluid Mech. 477, 19–49.
- Wilczek et al. (2011) Wilczek, M., Daitche, A. & Friedrich, R. 2011 On the velocity distribution in homogeneous isotropic turbulence: correlations and deviations from Gaussianity. J. Fluid Mech. 676, 191–217.
- Wilczek & Friedrich (2009) Wilczek, M. & Friedrich, R. 2009 Dynamical origins for non-Gaussian vorticity distributions in turbulent flows. Phys. Rev. E 80 (1), 016316.
- Yakhot (1989) Yakhot, V. 1989 Probability distributions in high-Rayleigh number Bénard convection. Phys. Rev. Lett. 63 (18), 1965–1967.