Resolving the fine-scale structure in turbulent Rayleigh-Benard convection

We present high-resolution direct numerical simulation studies of turbulent Rayleigh-Benard convection in a closed cylindrical cell with an aspect ratio of one. The focus of our analysis is on the finest scales of convective turbulence, in particular the statistics of the kinetic energy and thermal dissipation rates in the bulk and the whole cell. The fluctuations of the energy dissipation field can directly be translated into a fluctuating local dissipation scale which is found to develop ever finer fluctuations with increasing Rayleigh number. The range of these scales as well as the probability of high-amplitude dissipation events decreases with increasing Prandtl number. In addition, we examine the joint statistics of the two dissipation fields and the consequences of high-amplitude events. We also have investigated the convergence properties of our spectral element method and have found that both dissipation fields are very sensitive to insufficient resolution. We demonstrate that global transport properties, such as the Nusselt number, and the energy balances are partly insensitive to insufficient resolution and yield correct results even when the dissipation fields are under-resolved. Our present numerical framework is also compared with high-resolution simulations which use a finite difference method. For most of the compared quantities the agreement is found to be satisfactory.


Introduction
Turbulent fluid motion in nature and technology is frequently driven by sustained temperature differences [5]. Applications range from cooling devices of chips to convection in the Earth and the Sun. Turbulent Rayleigh-Bénard convection (RBC) is the paradigm for all these convective phenomena because it can be studied in a controlled manner, but it still has enough complexity to contain the key features of convective turbulence in the examples just mentioned. RBC in cylindrical cells has been studied intensely over the last few years in several laboratory experiments, mostly in slender cells of aspect ratio smaller than or equal to unity in order to reach the largest possible Rayleigh numbers or to resolve the detailed mechanisms of turbulent heat transport close to the walls [1,5]. Direct numerical simulations have also grown such that the detailed dynamical and statistical aspects of the involved turbulent fields and their characteristic structures can now be unraveled in detail.
The key question in RBC is the mechanism of turbulent transport of heat and momentum. Since the fluid motion is driven by a constant temperature difference between the top and bottom plates, thin boundary layers of temperature and velocity will form on these walls as well as on the side walls of the cell. A deeper understanding of the global transport mechanisms is possible only if we understand the dynamical coupling between the boundary layers and the rest of the flow in the bulk of the cell. While the boundary layers are strongly dominated by the presence of mean gradients of the temperature and velocity fields, the bulk of the convection layer is well mixed by the turbulence such that mean gradients of the involved turbulent fields remain subdominant compared to the local fluctuations. The flow at hand is thus strongly inhomogeneous, at least in the vertical direction, so it can be expected that the smallest dynamically relevant scales will differ when moving from the isothermal walls to the bulk.
Central to our understanding of the statistics of the turbulent transport is the role of the gradient fields of velocity and temperature which fluctuate extremely strongly at the small scales of the flow. This is a unique property of all turbulent flows. Dissipation rate fields -which measure the magnitude of these gradient fluctuations and are still inaccessible in experiments with respect to their three-dimensional structure [37] -are thus at the core of a deeper understanding of turbulence as a whole.
In the present work we want to make a further step forward with DNS of RBC by resolving fine scales never accessed before, both in the bulk and boundary layers, in order to study the statistics of the gradient fields, their joint extreme events, the statistical effect of rare high-amplitude events as well as Rayleigh and Prandtl number variation. A spectral element method is used to conduct the numerical studies [9,10]. It combines the flexibility in terms of mesh geometry that is inherent to every finite element method with the exponentially fast convergence of a spectral method. We will show that several tests which have been applied in DNS in the past are insensitive with respect to insufficient resolution. These tests, which are based on global averages of mean dissipation rates and the global mean heat fluxes, give correct results although the fine-scale structure of the turbulence is still under-resolved particularly in the bulk of the convection cell. In order to address these questions in detail we will present a comprehensive statistical analysis of the temperature and velocity gradient fields, in particular the related dissipation rates and dissipation scales.
It is crucial to resolve all the dynamically important scales to represent the flow faithfully when carrying out Direct Numerical Simulations (DNS) which involve no subgrid-scale parametrization. Several attempts have been made in order to derive resolution criteria starting with the pioneering work by Grötzbach [16], subsequent refinements of this criterion [32,2,28] and works with a focus to the fine resolution of the boundary layer dynamics [36,27,22]. Only recently the focus of DNS studies was shifted towards the bulk in a cubic convection cell [17] with a discussion of the scaling properties and statistics of the dissipation fields.
It is well-known that the gradients of the turbulent fields are most sensitive to insufficient resolution. Superfine resolution simulations in isothermal box turbulence [24,25,26] and in turbulent shear flows [14] have led to some enlightening results on the distribution of the finest scales in such flows and their relation to the smallscale intermittency. This intermittency is known to be coupled tightly to two highly fluctuating dissipation rates, one of the kinetic energy and the other of the thermal variance. The thermal dissipation rate is defined as where T (x, t) is the temperature field and κ the thermal diffusivity. The kinetic energy dissipation rate is defined as with the turbulent velocity field u(x, t) and the kinematic viscosity ν. The mean kinetic energy dissipation rate is related to the mean Kolmogorov scale η K which is the smallest mean scale when ν ≤ κ. The symbol · denotes an ensemble average which is calculated in numerical simulations as a volume-time average. In case of ν > κ, the smallest mean scale is determined by the (active) scalar field known as the Batchelor scale [3], η B . Both scales are defined as Here, P r is the Prandtl number and given by Both dissipation fields can be expected to fluctuate strongly exceeding their means by orders of magnitude [30,24]. Therefore it was suggested to generalize the classical dissipation and diffusion scales to local dissipation and diffusion scales [31,23] which are given by (see also a discussion in Hamlington et al. [14]) Both scales will pick up the highly intermittent fluctuations of the dissipation rates and can thus become smaller, but also larger than the mean scales which were defined in (3). Local dissipation scales have been studied in convection experiments by Zhou and Xia [42]. One main finding was that the distribution of the scales can be described by the same tools as in isothermal box turbulence. In the present work we will also access these scales and compare their distribution in different parts of the convection cell. The manuscript is organized as follows. In Sec. 2, we will discuss in brief the equations of motion and some central relations which will become necessary for our data analysis. Furthermore we briefly review the existing resolution criteria. First it is shown that global balance equation checks are insensitive to insufficient resolution. We also compare the results with a second-order finite difference method [34,35] which has been one workhorse in RBC over the past decade. In Sec. 3 we report our results. We study the statistics of the dissipation rate fields and calculate the local dissipation scales. Then we present a comparison of dissipation rate fields and scales as a function of Rayleigh number and Prandtl number using our very highest resolution. We summarize our findings and give a brief outlook at the end.

Boussinesq equations and further non-dimensional relations
We solve the three-dimensional Boussinesq equations numerically. The height of the cell H, the free-fall velocity U f = √ gα∆T H and the imposed temperature difference ∆T are used to rescale the equations of motion. The three control parameters of Rayleigh-Bénard convection are the Rayleigh number Ra, the Prandtl number P r and the aspect ratio Γ = D/H with the diameter D. This results in the following dimensionless form of the equations of motioñ ∂T ∂t where The variable g stands for the acceleration due to gravity and α is the thermal expansion coefficient. Throughout the study we set Γ = 1. Times are measured in free-fall time units, T f = H/(gα∆T ). At all walls of the simulation volume V no-slip boundary conditions for the fluid are applied,ũ = 0. The side walls are adiabatic, i.e., the normal derivative of the temperature field vanishes, ∂T /∂n = 0. The top and bottom plates are held at fixed temperaturesT = 0 and 1, respectively. In response to the input parameters Ra, P r and Γ, a turbulent heat flux from the bottom to the top plate is established. It is determined by the Nusselt number which is defined as Based on the volume average, we find N u V = 1 + √ RaP r ũ zT V,t which has to equal N u(z) for allz ∈ [0, 1]. The non-dimensional expressions for the two dissipation rate fields, (x, t) and T (x, t) are given by the following expressions: and thus˜ The kinetic energy dissipation rate is defined as and thus˜ Using equation (5) gives for the cases of P r ≤ 1. To simplify (15) we used the definition of the Grashof number Gr = (U f H) 2 /ν 2 = Ra/P r. The Batchelor scale follows as for P r > 1. For completeness, we also list two exact relations that can be derived from the balances of the turbulent kinetic energy and the scalar variance. They are given by [29] If we make use of (17), Eqns. (15), and (16) translate to Similar to the studies by Stevens et al. [32], we will use eqns. (17) to test different grid resolutions at a given set of parameters and define relative errors that measure the difference between the left and right hand sides of Eqns. (17) ‡ ) ‡ Note that the relative errors are fairly sensitive to the averaging time because of the large fluctuations in Nusselt number that can occur for these turbulent systems.
In the following, we will continue with the dimensionless quantities and omit the tildes for convenience.

Numerical methods
For the DNS studies in the present work two different numerical methods are used and compared, a second-order finite difference scheme and a spectral element method.

Finite difference method
The Boussinesq equations (6)-(8) are discretized on a staggered grid with a second-order finite difference scheme (FDM) which was developed by Verzicco and Orlandi [34,35]. The pressure field p is determined by a two-dimensional Poisson solver after applying a one-dimensional Fast Fourier Transform (FFT) in the azimuthal direction. The time advancement is done by a third-order Runge-Kutta scheme. The grid spacings are non-equidistant in the radial and vertical directions. In the vertical direction, the grid spacing is close to Tschebycheff collocation points. The simulation code is parallelized with MPI in combination with OpenMP.

Spectral element method
The equations are numerically solved by a spectral element method (SEM) [9] which has been adapted to our problem. The code employs second order time-stepping, using the backward difference formula BDF2 which results at time step n and for a step width δt in the following set of discrete equations (see also Eqns. (6)-(8))D P r RaÂ u n + 3 2δtB u n + ∇p n = f n u , with the corresponding boundary conditions. Here,D is the divergence operator and A the stiffness matrix which contains the Laplace terms. The quantityB is the mass matrix which will contain the Gauss-Lobatto-Legendre weights and the determinants of the Jacobian caused by the mapping to the deformed elements as diagonal entries. In order to arrive at (20)- (22) the Boussinesq equations are transformed into a weak formulation similar to other Galerkin methods. They are then discretized with the particular choice of spectral basis functions [6] which will be given further below. These basis functions allow for an exact evaluation of the integrals in the scalar products on the basis of the Gauss integration theorem. All flow fields are given in the Sobolev space H 1 (V ) in which the functions and their derivatives are square integrable. For this space it holds that C 1 (V ) ⊂ H 1 (V ) ⊂ C 0 (V ) [6]. The right hand sides, f n u and f n T of Eqns. (22), incorporate remaining terms from the BDF2 time derivative, the nonlinear convection which is obtained by second order extrapolation from steps n − 1 and n − 2, and the buoyancy.
The resulting linear, symmetric Stokes problem is solved implicitly. This system is split, decoupling the viscous and pressure steps into independent symmetric positive definite subproblems which are solved either by Jacobi (viscous) or multilevel Schwartz (pressure) preconditioned conjugate gradient iteration. Fast parallel solvers based on direct projection [33] or more scalable algebraic multigrid [12] are used for the coarsegrid solve that is part of the pressure preconditioner. For stabilization of the SEM, we perform de-aliasing by the use of over-integration of the convective term by a factor of either N +5 or 3(N +1)/2, where N is the polynomial order. We also filter out 5% of the energy in the N th mode for additional stabilization (see [11] for further information).
The basis functions ψ k (x) are Lagrangian interpolation polynomials of order N and composed of Legendre polynomials P k for the present study. They are given by with the Gauss-Lobatto-Legendre points ξ k . In contrast to a classical spectral or pseudospectral method the evaluation of spatial derivatives translates into matrix multiplications which have to be highly optimized (see the appendix for further details). The expansion in the three-dimensional case with a reference element Ω = [−1, 1] 3 is based on the tensor product formulation of the basis functions In the simulation the elements that sum up to the volume V are deformed. Hence an additional mapping (Jacobian) from the reference element to all elements needs to be incorporated. Clearly, the mapping of the coordinates and the matching of the velocity and temperature fields between elements enhances the numerical effort in comparison to the second order FDM. We estimated that production runs on the same number of cores for the same system size would be approximately 10 times slower. In turn, gradient fields are calculated on each element separately with an exponentially fast convergence.

Existing resolution criteria for direct numerical simulations
The first estimation of spatial resolution requirements for direct numerical simulations of Rayleigh-Bénard convection were made by Grötzbach [16]. His criteria for confined convection cells consisted of (I) resolving the steep gradients in the velocity and temperature near the walls with a sufficient vertical grid width distribution and (II) resolving the smallest relevant turbulence elements with a sufficiently small mean grid width. Based on tests of Nusselt number with a spectral code, Grötzbach's first criterion requires at least 3 nodes within the thermal boundary layer thickness for Prandtl numbers on the order of one or larger. For much smaller Prandtl numbers, more nodes may be necessary as the viscous boundary layer becomes much thinner than the thermal boundary layer. The second criterion translates to a relation between the mean grid width∆ and the mean dissipation or diffusion scale. For P r ≤ 1 this relation Figure 1. Left: Display of the horizontal primary node structure as used for runs SEM1 to SEM4 displayed in Table 1. Right: Display of the vertical primary node mesh for runs SEM6 to SEM9 (see Table 2). The stretching factors r are r = 0.91 for N e,z = 32, r = 0.95 for N e,z = 64, and r = 0.97 for N e,z = 96.
is∆ ≤ π η K and for P r ≥ 1 it is∆ ≤ π η B . Grötzbach then assumes that η K and η B can be approximated by the mean kinetic energy dissipation rate as in Eq. (3). By using an argument similar to Eqns. (17) this leads to the following global criteria on the grid widths §: The criteria of Grötzbach were revised by Stevens et al. [32] based on DNS results using the second order finite difference method also used for comparison in this paper [34,35]. They systematically found the Nusselt number to be overestimated in poorly resolved simulations, especially when the plume dynamics were not properly resolved. They suggested changing the mean grid width criteria to one that instead holds for the largest grid width in any spatial dimension, since the Kolmogorov length needs to always be resolved in order to properly characterize the flow. A similar perspective was developed in Bailon-Cuba et al. (2010). Although Stevens et al. (2010) did not determine any exact resolution criteria, they did compute the volume averaged dissipation rates and T and compared these values to the globally computed § The factor of π can be rationalized to our view by the resolution criteria as formulated for pseudospectral box turbulence simulations (see e.g. [23]). There, k max η K ≥ 1 should be satisfied with the maximum resolved wave number (after de-aliasing) k max = √ 2N x /3 and N x = N y = N z equals the number of grid points in each direction. The standard box length is then the periodicity length of the Fourier modes L x = 2π and thus N x = 2π/∆.  Table 1). We compare the mean temperature profile T (z) A,t , the mean convective flux profile √ RaP r u z T A,t , vertical profiles of mean thermal dissipation T (z) A,t and mean kinetic energy dissipation (z) A,t , a z-dependent Kolmogorov scale η K (z) A,t and how well the Grötzbach criterion is satisfied plane by plane ∆z e (z)/ η K (z) A,t . The dashed line in the lower right panel marks ∆z e (z)/ η K (z) A,t = π. The dashed line in the upper mid panel is the Nusselt number from run FDM3 (see Table 2 or Ref. [2]). The inset in the top right panel magnifies the vertical profile of the mean thermal dissipation rate.
Nusselt number as we have done in equation (17). They found that for high enough Rayleigh number (≥ 10 9 ), even though the Grötzbach criteria was technically followed, and equation (17) was well-satisfied for the viscous dissipation rate, equation (17) was not as well-satisfied for the thermal dissipation rate.
A further revision was conducted by Shishkina et al. [28], who used the theoretical Prandtl-Blasius (PB) theory to derive a lower bound on the number of nodes required to be placed in both the thermal and the viscous boundary layers such that the estimated Kolmogorov lengths in the boundary layers are adequately resolved. For higher Rayleigh number, this minimum bound is much larger than that suggested by Grötzbach. For example, for our parameter range (Pr = 0.7), Shishkina et al. suggest a minimum of 5 nodes for Ra = 2 × 10 7 but increasing to 9 nodes for Ra = 2 × 10 9 .
We will also discuss our own results in light of these criteria, including Grötzbach and the revisions by Stevens and Shishkina. However, we will take this analysis one step further by investigating the implications of resolving not only the global but also the local dissipation scales.  Table 1. Parameters of the different spectral element simulations SEM1 to SEM4. The runs have an identical primary node mesh, but different polynomial order on each element. We display the order N of the Legendre polynomials, the total number of spectral elements, N e , the number of spectral elements with respect to z direction, N e,z , the number of grid cells resulting from primary and secondary nodes with respect to z direction, N z = N e,z N , and the Nusselt numbers N u(z = 0), N u(z = 1) and N u V . Furthermore we list the relative errors Λ T and Λ v (see Eqns. (19)). All four runs are at Ra = 10 8 , Γ = 1 and P r = 0.7.

Statistical properties for resolutions with different polynomial orders
In correspondence with the so-called h-type and p-type spectral element methods (SEM), two routes of modification of the resolution exist. In the h-type SEM the number of primary elements, N e , is varied, in the p-type SEM one changes the polynomial degree N of the basis functions on each element and keeps the number of elements fixed. In the following, we summarize efforts in both directions in order to study resolution effects for the gradient fields. Since the grid is non-uniform in all three directions the side lengths of an element are functions of the three coordinates, i.e. ∆x e (x, y, z), ∆y e (x, y, z) and ∆z e (x, y, z). Figure 1 (left) shows a view of the horizontal primary element mesh. The coarsest elements are always found at the cell center line. Particular emphasis was given here to the vertical resolution since this is the important direction for the correct resolution of the BLs. The formula that has been chosen to determine the element boundaries in the vertical direction is given by the following geometric scaling for the upper half of the cell with the scaling factor r In correspondence with the up-down-symmetry this relation has to be applied for the lower half as well. Equidistant vertical meshing corresponds to r = 1. Figure 1 (right) demonstrates the resulting vertical meshing for different numbers of primary element nodes, N e,z . In the appendix, we describe one way to obtain an optimal non-equidistant grid with respect to z, in other words, an optimal scaling factor r in (27). Figure 2 shows important statistical quantities for a variation in correspondence with a p-type refinement. Results are obtained for different polynomial orders but the same primary element mesh. The results are summarized in Table 1. In all cases shown, derivative-based quantities are evaluated spectrally on each element and no derivatives are taken across boundaries. All runs are conducted at a Rayleigh number Ra = 10 8 with  Table 1. Data are normalized with respect to mean Kolmogorov length. The definitions are given in (30) and (31). Right: Mean dissipation rates of thermal variance, T V,t , and kinetic energy, V,t , as functions of the polynomial order. a primary element mesh as displayed in Figure 1. On average, we ran our simulations for at least 30 free-fall times T f to ensure that the system had settled into its relaxed state, and then we continued the evolution for at least 75 free-fall times (in case of the biggest DNS), outputting on average at least 80 statistically independent snapshots.
Both the table and the figure indicate that insufficient spectral resolution is manifested in multiple ways, but is not necessarily obvious when looking at standard quantities, e.g. the ingredients for the turbulent heat flux. The graphs for the mean temperature profile T (z) A,t , the convective heat flux √ RaP r u z T A,t and even the Nusselt numbers which are obtained in different ways do not indicate a resolution effect at first glance. However the large magnitude of the relative errors of run SEM1 which has been used to test the dimensionless energy balances (19) is definitely caused by the insufficient resolution. The run SEM1 was one of our longest, taking about 1200 T f . The statistical analysis is based on 192 turbulence statistically independent threedimensional snapshots separated by 6 free-fall times each.
An increase to N = 5 in run SEM2 improves the convergence of the energy balances drastically. Nevertheless, the plane-averaged mean profiles of the thermal and the kinetic energy dissipation rates still display an insufficient resolution which is present for the next higher polynomial order, N = 7, as well. This becomes visible by the discontinuities at the element boundaries, especially near the center of the cell. The lower right panel of Fig. 2 relates a refined Kolmogorov-type scale We can refine the classical Grötzbach criterion (26) to where ∆z is the vertical grid spacing, recognizing now element mesh and collocation grid on each element. This is exactly what produces the characteristic shape of all the curves in the lower right panel. Such a criterion has been suggested already by [2]. It shows clearly that all orders N ≤ 7 result in grid spacings that are too coarse in the bulk region of the convection cell. In the appendix, we explain how insufficient resolution can cause the spike structures in the vertical profiles of the dissipation rates and related quantities by means of a convergence test for a simple analytical profile. The artifacts at the element boundaries which we see for the SEM are due to insufficient resolution and hence the failure of the derivatives to match at the boundaries, since this SEM method enforces only the continuity of the functions at the boundaries. This gives rise to a clear criterion for resolution: when the system is sufficiently resolved, all spikes in both dissipation profiles completely disappear. In classical finite difference methods numerical diffusion and dispersion will suppress such spikes. In addition it is shown in the appendix that similar observations as in Fig. 2 follow when the h-type route of grid refinement is followed, i.e., refining the primary mesh at fixed system parameters (Ra, P r, Γ) and a given polynomial order. To conclude this part, while some standard indicators for sufficient resolution which have been discussed in previous works [32,28] are all well-satisfied, a closer look at the dissipation fields indicates clearly that the spatial resolution is not sufficient, in particular in the bulk of the cell. The artifacts in the mean vertical profiles of the gradient fields do not completely disappear even when the order is increased to N = 11 as demonstrated in the inset in the top right panel of Fig. 2.
Compared to previous resolution studies of fluid turbulence in periodic boxes [24] and shear flow turbulence channels [14], the situation in the present RB case is more complex. On the one hand, the turbulent flow is inhomogeneous in all space dimensions. This causes space-dependent statistical properties of the turbulent fields and their derivatives. On the other hand, the computational grid is non-uniform in all three directions as described already above. Although we refine the grid towards all walls, the regions where one expects the largest amplitude of the derivatives, it is not necessarily assured that both, steepest gradient and finest grid cells, coincide. In this situation one can however define the coarsest and finest grid spacing in the whole cell or a subvolume to get a global indication of the quality of resolution. This is done by the following geometric means ∆ max = 3 max x∈I ∆x e (x, y, z) max y∈I ∆y e (x, y, z) max z∈Iz ∆z e (x, y, z) with I = [−0.5, 0.5] or a subinterval and I z = [0, 1] or a subinterval, such as the bulk volume V b which is given further below in the text. In Fig. 3 (left) we display the minimum and maximum grid spacing for the whole cell obtained by (30) and (31). Furthermore, we show the minimum resolution in the bulk region where we defined a subvolume V b = {x = (r, φ, z) | 0 ≤ r ≤ 0.3; 0.2 ≤ z ≤ 0.8}. The right panel confirms what we have discussed already above, that the mean dissipation rates level off for N ≥ 5 although vertical profiles are still not sufficiently well resolved. In Figure 4 we display the probability density functions (PDFs) of the fields T and . The upper row shows data which have been obtained in the whole cell, the lower row those for the subvolume V b in the bulk. It can be seen for all four panels that with increasing polynomial order more very-high-amplitude events are resolved and that the tail is further stretched out. The better resolution manifests in significantly less scatter at the largest amplitudes. Even more pronounced are the resolution effects in the bulk (lower row). We observe now for both dissipation rates the same systematic trend. The tail of the stretched exponential distribution is fatter for higher polynomial order. This latter finding is also in agreement with previously reported spectral resolution studies for homogeneous isotropic box turbulence as reported in [24]. Note that the tails of the exponents in our figure do not always increase in an even manner with resolution. One sees for example, a jump in the lower left panel of Figure 4 in going from N=3 to N=5 and then again from N=7 to N=11. This is understandable in light of Figure 12, where we see that high-amplitude events can increase the tails significantly. Since the system is chaotic as well as turbulent, simulations done for the same Rayleigh number Figure 5. Comparison of spectral element runs SEM9 and SEM10 with finite difference run FDM1 (see Table 2 for specifications). The same quantities are plotted as in Figure 2. The dashed line in the lower right panel marks ∆z(z)/ η K (z) A,t = π. The dashed line in the upper mid panel is the Nusselt number of the FDM run from [2]) with N u = 64.3. The inset in the top right panel magnifies the thermal dissipation rate profiles for FDM1 and SEM10.

Run
Ra  Table 2. Parameters of the different spectral element simulations. Runs SEM7, SEM8 and SEM10 have an order of the Legendre polynomials N = 11, runs SEM5, SEM6 and SEM9 use N = 7. We show the Rayleigh number Ra, the Prandtl number P r, the total number of spectral elements, N e , the number of spectral elements with respect to z direction, N e,z , the number of grid planes (primary and secondary nodes) with respect to z direction, N z , the total number of grid cells N e N 3 , and the Nusselt number N u V . Furthermore we list the relative errors Λ T and Λ v (see Eqns. (19)). All simulations are conducted for Γ = 1 and P r = 0.7. The finite difference run FDM1 is conducted for N φ × N r × N z = 721 × 361 × 621, respectively. but different resolution are statistically different, so some of them could have more high amplitude events than others. Longer simulation times would help smooth this out.

Comparison with the finite difference method
The comparison of DNS runs at Ra = 10 9 and P r = 0.7, i.e. runs SEM9, SEM10 and FDM1 from Table 2, is displayed in Figure 5. The resolution of FDM1 has been chosen twice as fine as in the run from Bailon-Cuba et al. [2] in order to get a comparable number of grid points with respect to SEM9. The thermal boundary layer for FDM1 is resolved with 21 grid planes. We see that the agreement of the mean vertical profiles is very good. The value of the Nusselt numbers for the higher resolution runs match now to within three significant figures, also with the data from Ref. [36]. In Table 2, we list also the corresponding relative errors Λ T and Λ v which are larger for FDM1 than those for SEM9 and SEM10, particularly for the kinetic energy dissipation rate in our study. The latter has been evaluated in correspondence with (13) which has to be applied for inhomogeneous flows. The difference is supported by the deviation in the vertical profiles of (z) A,t in the lower left panel of Fig. 5. We also note that Stevens et al. (see their Table 1 in [32]) reported similar errors which were in their case however mostly detected for the thermal energy dissipation rate and not only for the lowest resolutions. The distribution of the local amplitudes of both dissipation fields is compared in Fig. 6. Both panels show that the deviations arise mostly for the outer tails where the extreme fluctuations are captured. In case of the thermal dissipation rate both PDFs remain closely together for nearly the whole range. The kinetic energy dissipation rate data start to differ for roughly 15-20% of the maximum amplitude which might be one reason for the larger values of Λ v . The agreement in the low-amplitude part of the PDFs is good as shown in both insets.
Furthermore, we find excellent agreement between the two codes when comparing global transport properties. For example, the Nusselt number and the globally averaged thermal dissipation rates for both the whole cell and the bulk volume V b as shown in Fig.  7 agree quite well. We varied our Rayleigh number between 10 6 and 10 9 and compared with fits of the FDM data from [2], [8] and [7], respectively. We did need to use a different prefactor for the bulk-averaged thermal dissipation rate since our subvolume V b was chosen differently.
To estimate the effect of the size of our subvolume on the thermal dissipation rates, we show two addtional data sets in the right panel of Fig. 7: one for a smaller subvolume 8V b /27 and the other for an even smaller subvolume of V b /27 (all centered about the middle of the cell). The general trend is for the thermal dissipation rates to slightly decrease as the subvolume decreases. Also, our uncertainty becomes larger as the subvolume decreases. We estimated the uncertainty in the mean T V,t values by computing the difference between the mean taken over the entire time series and the mean taken over only the latter half of the time series.
The fits to the data sets corresponding to the smallest subvolumes are T 8V b /27,t = (0.21 ± 0.07)Ra −0.40±0.02 and T V b /27,t = (0.25 ± 0.12)Ra −0.42±0.03 . Kaczorowski and Xia [17] also studied the scaling of subvolume-averaged thermal dissipation rates in a similar range of Rayleigh numbers using a small subvolume (V /64) but for a Prandtl number of 4.38. Our exponent disagrees with theirs of T V,t = 43.9Ra −0.78 . We do see a trend towards a larger exponent as our subvolume decreases, but our largest exponent still disagrees with [17] even when including our estimates of numerical uncertainty.

Very-high-resolution runs at different Rayleigh numbers
In the following section, we want to discuss a series of very-high-resolution runs in more detail. All the runs with their resolution are displayed in Table 2. We first compare runs at P r = 0.7 spanning a Rayleigh number range from 10 7 to 10 9 .
Snapshots of high-amplitude regions of both dissipation fields are shown in Figure 8. The data are given in logarithmic units. Both dissipation rates form smooth sheet-like structures in the bulk, in particular the thermal dissipation rate. The very fine resolution is clearly obvious from the absence of ripples at the isosurfaces of both dissipation fields. In Figure 9 we show horizontal slices of both dissipation fields at fixed height z. Again the data is given in logarithmic units to highlight the variation. The top row is for Figure 7. Comparison of global transport properties between SEM and FDM. Left: Nusselt versus Rayleigh number for 10 6 < Ra < 10 9 . The data are compensated by the power law N u = 0.145 × Ra 0.294 which was a fit to the data as reported in Bailon-Cuba et al. [2]. Right: mean thermal dissipation as a function of the Rayleigh number for different subvolumes. The first two data sets are compared with fits to the FDM data (shown as dashed lines). For the whole cell V we take former results from Ref. [8], for the bulk volume V b we compare with data from Ref. [7]. In this case the prefactor is different since the subvolume V b was chosen differently. The last two series are obtained in smaller subvolumes and are fitted by power laws as given in the legend and shown as solid lines.

detachment.
The distribution of the locally fluctuating dissipation scales η K (x) as defined in (5) is shown Fig. 10 for runs SEM7, SEM8 and SEM10. The scales have been analyzed in the whole cell with volume V as well as in a bulk region which is defined by V b . The definition (5) has been chosen for this analysis which can be straightforwardly applied to the non-uniform grids that have been used for all DNS. An alternative definition of local dissipation scales which is based on velocity increments was suggested in [39,40]. In Ref. [14], it was shown how both distributions can be related to each other. It can be observed first that the scales in the whole cell cover a wider range, both, to the large-and small-scale end (see top left panel) which is centered around the most probable value which is always close to mean dissipation scale η K V,t which is calculated following (18). This finding is also in agreement with previous DNS results [7,8] which show that dissipation rates have significantly higher amplitudes in the boundary layers. We also see that the right tail ends of the distributions in the whole cell decrease with Rayleigh number. It demonstrates that the scales in turbulent RBC become finer as the Rayleigh number increases. This argument is also supported by the fact that the differences between the distributions in V and V b become smaller. In the top right  Table 2.
figure, we zoom into the left tail end for all six data sets. The smallest local dissipation scales are associated with the largest dissipation events which arise for very steep local gradients. With increasing Rayleigh number these contributions become larger, i.e. the left tail becomes fatter. A similar behavior, however much less pronounced, can be observed if one restricts the analysis to the bulk volume V b . In Figure 11, we display the PDFs of both dissipation rates in the whole cell and in the bulk. For both rates, it can be clearly seen that the major contribution to the high-amplitude events comes from the boundary layer regions. This has been studied already in [7]. Figure 12 illustrates the sensitivity of the statistics with respect to a single extreme event that was monitored in the course of the simulation run and can be identified as a large scale plume sweeping through the bulk volume. It causes a large instantaneous thermal dissipation which is not easily detectable in the mean dissipation T (t) V b . Only in the fourth moment of the thermal dissipation 4 T which is taken in the bulk volume does this strong event become clearly visible as seen in the right panel of figure 12. Note also that the fourth moment in the whole cell is also fairly insensitive to this particular high-amplitude bulk event. The impact of this single event on the statistics is shown in the left panel of Figure 12. As expected the tail is stretched significantly.
In Fig. 13 we display a time sequence of the dynamics which is associated with this single extreme event. The upper panel of the figure replots the moments of the thermal dissipation rate obtained in V b with respect to time but on a finer time scale. It can be observed that the entire event lasts less than two free-fall times. We find that the fourth-order moment increases by about three orders of magnitude within this short period. Figures 12 and 13 clearly show that the statistical fingerprint of this strong event is best detected in the higher-order moments. The bottom panels of Figure 13 show vertical slice images of the thermal dissipation rate and the temperature corresponding to four different times in the evolution of this event. Clearly visible is the pronounced hot plume rising and then detaching from the bottom plate which generates steep temperature gradients and thus a large amplitude of the thermal dissipation rate.

Very-high-resolution run at higher Prandtl number
Lastly, we compare the gradient statistics at a given Rayleigh number for two different Prandtl numbers. Runs SEM7 and SEM7a are conducted at Ra = 10 7 and P r = 0.7 (air) and 6 (water). The data in Table 2 indicates already that the resolution requirements remain the same for the enhancement of the Prandtl number by a factor of nearly 10. For P r > 1, the mean diffusion scale of the temperature field, η B is smaller than mean Kolmogorov scale, η K , since a viscous-convective range on scales smaller than the Kolmogorov scale builds up. Figure 14 displays the distributions of η k (x, t). We observe again that the range of varying scales is larger when the data are taken in the whole cell in comparison to the bulk volume. Figure 11. Probability density function (PDF) of the dissipation rates for runs SEM7, SEM8 and SEM10 as given in the legend and Table 2. We compare the PDFs obtained in the whole convection cell (upper panels) and those obtained in the bulk which is defined as the subvolume V b (lower panels). Thermal dissipation rates are displayed in the left column, energy dissipation rates in the right one.  We have not rescaled the distributions by the corresponding mean dissipation scale since we want to point to the shift of both PDFs for P r = 6. This means that the local dissipation scales are larger as a whole than for the case of P r = 0.7. This behavior looks counter-intuitive at first glance, particularly from the perspective of passive scalar mixing at increasing Prandtl (or Schmidt) number [23]. There one detects increasingly finer diffusion scales for the passive scalar leaving however the local dissipation scales unchanged. For turbulent convection, we estimated already in the introductory part that the relatively slow falloff of η B ∼ P r −1/8 as we progress from P r ≈ 1 to P r 1. It is the active nature of the temperature field which causes this different behavior in the convection case as compared to passive mixing. A temperature field at a higher Prandtl number exists on finer scales than a velocity field, obeying narrower plume structures which causes a weaker driving of the fluid motion resulting in less steep velocity gradients and consequently larger local dissipation scales. This argumentation is also supported by a comparison of the PDFs of both dissipation fields as shown in Fig. 15. High-amplitude events and tails are shifted to smaller magnitudes in all cases.

Joint Statistics
The joint statistics of both dissipation rates is displayed in Fig. 16. We show the joint and normalized probability density function which is given by The contour levels are plotted in logarithmic values as indicated by the color bar. Similar to [17] for RB convection or to [15] for a channel flow, the support of Π has an ellipsoidal form with a tip at the joint high-amplitude events. The joint PDF P ( T , ) is here normalized by the corresponding single quantity PDFs, P ( ) and P ( T ). In this way we highlight the correlations between both fields. If Π( T , ) is larger than unity then the correlation is larger than if the two dissipation rate fields were statistically independent. It can be observed that the support of the joint PDF for P r = 6 is shifted to smaller amplitudes in comparison to P r = 0.7, which is in agreement with Fig. 15. In both cases, the high-amplitude events are correlated strongest, exceeding the corresponding mean amplitudes by at least two orders of magnitude.

Summary and discussion
We have computed both global and local measures of dissipation and heat transport from high resolution direct numerical simulations of turbulent Rayleigh-Bénard convection using a spectral element method. We find that the global measures of heat transport, such as Nusselt number, time-averaged temperature profiles, and volume-averaged dissipation rates, are fairly insensitive to insufficient resolution, as long as the mean Kolmogorov length is resolved. However, if one computes instead plane-averaged or even more local dissipation rates, one finds that the Grötzbach criteria (or something even more stringent as in (33)) needs to be satisfied for every grid point in order to have the system properly resolved. The main effects of a poorly resolved simulation are that some of the largest dissipation (both thermal and viscous) scales in the system are not resolved, especially in the bulk where the computational grid is coarsest. Our investigations suggest that the refined SEM analysis which we conducted to study the statistics of dissipation fields require at least This follows e.g. from the data displayed in Figure 4 for the largest Rayleigh number. It is clear that such a criterion can be applied a posteriori only. Recall also that the horizontal spacing was always finer in the present cases such that a geometric mean remains smaller than ∆z(z). We have also compared our SEM results with a FDM code and find excellent agreement for global quantities such as Nusselt number and temperature profiles, and even fair agreement with globally averaged dissipation rates. The only discrepancy is with the vertical profiles of the mean dissipation rates, which disagree by about 9%.
Once we determined our resolution criteria, we then compared local dissipation rates ( , T ) and the local dissipation scale η k as a function of Rayleigh number for our sufficiently-resolved simulations. Local dissipation scales can be considered as a generalization of the mean Kolmogorov dissipation scale which incorporate the spatially intermittent nature of the energy dissipation field. Local scales below the Kolmogorov scales are related to strong local gradients or high-amplitude dissipation events. We find that the local dissipation scales in the entire cell have a wider range of values than the dissipation scales in the bulk of the cell. But in all cases, there is a fairly wide range of dissipation scales both above and below the mean Kolmorgov dissipation scale. The range of these local scales is a manifestation of the intermediate dissipation range (IDR) which exists in the crossover region between the inertial and viscous range. The IDR was developed in the multifractal formalism [20,19,13,4]. Similar to previous studies in box turbulence and channel flow turbulence, this range increases as the Rayleigh number grows. We have shown here that the dissipation scales on the left end of the PDFs become smaller as Rayleigh number increases, and correspondingly the probability of largest dissipation scales decreases. We also found, by looking at the fourth moment of the thermal dissipation rate, that high-amplitude but rare dissipation events can dominate the tails of the PDFs of the thermal dissipation rates. This highlights the sensitivity of turbulent RBC to such rare, but extreme events and calls for caution when generalizing statistical quantities in turbulent RBC. We also computed the joint statistics of the kinetic energy and thermal dissipation rates and find that the high amplitude events are the most strongly correlated.
Finally we compared results at two different Prandtl numbers. The range of local dissipation scales becomes smaller when P r > 1 which is in line with smaller amplitudes of both dissipation rates. Our estimates (18) indicate that the resolution demands grow significantly when the Prandtl number is decreased starting from P r ≈ 1. Equation (18) suggests a stronger Prandtl number dependence on the dissipation scales, namely η K ∼ P r 3/8 , for cases decreasing from P r ≈ 1 to P r 1 than for those which increase from P r ≈ 1 to P r 1 (where η B ∼ P r −1/8 ). On the numerical side, a second challenge appears that is related to the high diffusivity of temperature field and which has been discussed recently for the case of passive scalar mixing at very low Schmidt number [41]. An explicit time advancement becomes increasingly demanding since the scalar relaxes increasingly faster. Preliminary studies suggest e.g. that for the Prandtl number of mercury (P r = 0.021) a mesh is necessary that equals the one which we used for P r = 0.7 for a Rayleigh number larger by a factor of one hundred.

Appendix. Additional resolution studies
Sensitivity with respect to vertical element spacing In this section, we describe in brief one way to obtain an optimal scaling factor r given by (27). As r becomes smaller, the elements become more clustered towards the boundary plates as shown in Figure 17. Table 3 summarizes ten different test runs at fixed Ra, P r and Γ. In all cases the horizontal mesh (see again Figure 1) and the total number of elements N e remain unchanged. We varied the polynomial order N and r only. Figure 18 shows time series of the Nusselt numbers for the different values of r obtained at z = 0, i.e., N u(t) = −∂ T A /∂z| z=0 . The values N u(t) fluctuate about their temporal means. These fluctuations do not decrease when N is increased, i.e. when the resolution is improved (not shown). A systematic effect for an increase of r is clearly visible in the insets of both figures, where we report the time averages of N u(t) at both plates with the error bars corresponding to the standard deviation σ. Neither an equidistant nor a strongly non-uniform grid are preferable since they give the largest discrepancy in Nusselt number. There is a trade-off between resolving the boundary layers (non-uniform grid) and the bulk (equidistant grid). Based on our analysis here, scaling factors of about r ≈ 0.9 seem to be the optimum and were kept for the rest of the studies. For the present studies, we have chosen r = 0.91 and have also matched finer primary element meshes correspondingly.  [2] and the magenta dashed line from [36]. The number of elements is the same in both series. Top: polynomial order N = 3 with runs T1, T3, T5, T7 and T9. Bottom: N = 5 with runs T2, T4, T6, T8 and T10. The insets in both figures display N u as obtained by a time average at z = 0 (blue circles) and 1(red stars) as well as the corresponding error bars. For comparison we add again N u ± σ from [2,36] in the same color style as in the main figures. All data are for Ra = 10 9 and Γ = 1 (see Table 3).  Table 3. Parameters of the different spectral element test runs T1 to T10 with different vertical spacing. We display the order N of the Legendre polynomials, the geometric stretching factor r in the vertical direction, the total number of spectral elements, N e , the number of spectral elements with respect to z direction, N e,z , and the number of grid cells resulting from primary and secondary nodes with respect to z direction, N z = N e,z N . In all cases, Ra = 10 9 , Γ = 1 and Pr = 0.7.

Run N
(N e , N e,z )  Table 4. Parameters of the different spectral element simulations T11 to T13. The three runs have different primary node meshes, but the same polynomial order N = 7. We display the order N of the Legendre polynomials, the total number of spectral elements, N e , the number of spectral elements with respect to z direction, N e,z , the number of grid cells resulting from primary and secondary nodes with respect to z direction, N z = N e,z N , and the Nusselt numbers N u(z = 0), N u(z = 1) and N u V . Furthermore we list the relative errors Λ T and Λ v (see Eqns. (19)). All runs are at Ra = 10 9 , Γ = 1 and P r = 0.7. Note that T13 equals SEM9.

Complementary series of resolution tests
A complementary series of resolution tests in comparison to those reported in Sec. 3.1 is presented in Table 4 and Figure 19. In the test runs T11 to T13 we varied the primary element meshes and left the polynomial order of each element unchanged. The outcome from this series is similar to what was already demonstrated in the main text. While the vertical profiles for the mean temperature or the mean convective heat flux are practically equal, differences manifest for the gradient fields (see Figure 19 for details). . Resolution tests for Ra = 10 9 , Γ = 1, and N = 7 using different primary meshes (runs T11 to T13 in Table 3). We compare the same quantities as in Figure 2. The dashed line in the lower right panel marks ∆z e (z)/η K (z) = π. The dashed line in the upper mid panel is the Nusselt number from run FDM1 as displayed in Table 2.

Spatial derivatives in the spectral element method
In order to illustrate how we take very accurate derivatives in the SEM, we use as an example a one-dimensional case with the reference element Ω = [−1, 1]. A spectral approximation of a function u e ∈ L 2 w (Ω) (with w(x) being a positive weight function) can be written as follows where ψ k (x) is the kth order basis function and the (N+1) points are the nodes of the Gauss-Lobatto quadrature. They are determined by the Gauss-Lobatto integration theorem [6]. For the approximation, one has to take a set of polynomials which form an orthogonal system of the underlying Hilbert space of square-integrable functions L 2 w (Ω). The first derivative of the function u e (x) at the GLL points is Starting from Eq. (23) together with the relation (1 − x 2 )L (x) = 0 for x = ξ k and substituting the Legendre differential equation ((1 − x 2 )L (x)) = −N (N + 1)L N (x) the derivative becomes We also recall that L N (−1) = (−1) N and L N (1) = 1 for all N . Thus, the derivative at the boundary x = ξ 0 is given by at x = ξ k for 0 < k < N by and at x = ξ N by In Figure 20 we summarize the results for a simple function u(x) = cos(πx/2). With a view to dissipation rates we are interested in the accuracy for quantities that contain (du/dx) 2 . In the left panel of Figure 20, we compare the derivative as obtained from (37)- (39). We see that the errors at the boundary result in strong overshoots at the element boundary which are amplified by the second power of the derivatives as in the dissipation rates. The mid panel of Figure 20 repeats the analysis for a primary element node mesh of half the size obtained here byx → (x − 1)/2. An increase in resolution of the primary element node mesh reduces the errors significantly. The error can be quantified by which is shown in the right panel of Figure 20. The exponential convergence with respect to the polynomial order N is clearly demonstrated in the right panel of the figure. Thus a combination of both increasing the polynomial order and the node mesh leads to an accurate calculation of spatial moments.