Transition to ultimate Rayleigh-B\'{e}nard turbulence revealed through extended self similarity scaling analysis of the temperature structure functions

In turbulent Rayleigh-B\'{e}nard (RB) convection, a transition to the so-called ultimate regime, in which the boundary layers (BL) are of turbulent type, has been postulated. Indeed, at very large Rayleigh number $Ra \approx 10^{13}-10^{14}$ a transition in the scaling of the global Nusselt number $Nu$ (the dimensionless heat transfer) and the Reynolds number with $Ra$ has been observed in experiments and very recently in direct numerical simulations (DNS) of two-dimensional (2D) RB. In this paper we analyse the local scaling properties of the lateral temperature structure functions in the BLs of this simulation of 2D RB, employing extended self-similarity (ESS) (i.e., plotting the structure functions against each other, rather than only against the scale) in the spirit of the attached eddy hypothesis, as we have recently introduced for velocity structure functions in wall turbulence (Krug \textit{et al.}, \textit{J. Fluid Mech.}, vol. 830, 2017, pp. 797-819). We find no ESS scaling below the transition and in the near wall region. However, beyond the transition and for large enough wall distance $z^+>100$, we find clear ESS behaviour, as expected for a scalar in a turbulent boundary layer. In striking correspondence to the $Nu$ scaling, the ESS scaling region is negligible at $Ra = 10^{11}$ and well developed at $Ra = 10^{14}$, thus providing strong evidence that the observed transition in the global Nusselt number at $Ra \approx 10^{13}$ indeed is the transition from a laminar type BL to a turbulent type BL. Our results further show that the relative slopes for scalar structure functions in the ESS scaling regime are the same as for their velocity counterparts extending their previously established universality. The findings are confirmed by comparing to scalar structure functions in 3D turbulent channel flow.


Introduction
Thermal convection is relevant to a wide range of applications across various fields such as building ventilation (e.g. Linden 1999), atmospheric (e.g. Hartmann et al. 2001) or oceanic (e.g. Rahmstorf 2000) flows. The phenomenon is widely studied in terms of the paradigmatic case of Rayleigh-Bénard (RB) convection (see the reviews of Ahlers et al. 2009;Lohse & Xia 2010;Chillà & Schumacher 2012), in which two horizontal plates at a distance L are cooled from above and heated from below. In the presence of a gravitational acceleration g, a flow is driven with properties (for given Prandtl number P r ≡ ν/κ and aspect ratio Γ = D/L) that are controlled by the Rayleigh number Ra ≡ αg∆Θ/(νκ). Here, ν and κ are kinematic viscosity and the thermal expansion coefficient, respectively, D is some lateral length scale and ∆Θ ≡ Θ b − Θ t denotes the temperature difference between the bottom (held at Θ b ) and top (Θ t ) plates. The primary response parameter of the system is the resulting heat flux, which in its non-dimensional form is given by the Nusselt number N u relating the actual heat transfer to the one in the purely conductive case. The behaviour of the system at very high Ra is of interest in many applications, and theoretical work (Kraichnan 1962;Spiegel 1971;Grossmann & Lohse 2000, 2011 predicts the existence of a so called 'ultimate regime', in which the scaling N u ∝ Ra β switches from the classical β 1/3 (Malkus 1954) to β > 1/3. This transition is related to a change of the boundary layer (BL) structure, of both velocity and temperature, from a laminar to turbulent type. A transition in the N u scaling -as well as in the Reynolds number scaling-has indeed been observed experimentally, most convincingly by He et al. (2012bHe et al. ( ,a, 2015 and very recently also in a direct numerical simulation (DNS) of two-dimensional (2D) RB by Zhu et al. (2018). Results of the latter are reproduced in figure 1a, where the transition is evident from the change in slope of the (compensated) N u at Ra ≈ 10 13 . As expected from the theory and indicative of the turbulent nature of the BLs, Zhu et al. (2018) found logarithmic dependencies with respect to the distance from the wall for the mean temperature and velocity BLs. In the following, we will characterize the thermal BLs further by studying their lateral structure functions.
Such an analysis originated from Davidson et al. (2006), who pointed out that the k −1 spectral scaling (k being the streamwise wavenumber) of the streamwise velocity power spectrum predicted by the attached-eddy hypothesis (Townsend 1976;Perry & Chong 1982), is equivalent to -and more readily observed as-a ln(r/z) scaling of the second order longitudinal structure function. Here and in the following, r denotes the streamwise separation distance and z is the distance off the wall. de Silva et al. (2015) found that the ln(r/z) scaling also applies to velocity structure functions of arbitrary even order 2p such that for the so-called energy-containing range z < r δ (δ being the boundary layer thickness) (1.1) Here ∆u(r, z) is the velocity increment between two points at a distance z off the wall separated by r along the streamwise direction, U τ is the mean friction velocity and E u p , D u p are constants. We use superscript u to denote quantities relating to the velocity and θ when referring to the scalar later on. While the direct scaling according to (1.1) is only observed at relatively large Reynolds numbers Re τ = δU τ /ν ∼ O(10 4 ), de Silva et al. (2017) and Krug et al. (2017) have shown that a relative scaling is evident at much smaller Re τ if the so called extended self-similarity (ESS) framework is employed. In this case, the scaling is not analysed as a function of r but -in the spirit of the original ESS hypothesis by Benzi et al. (1993Benzi et al. ( , 1995 -relative to a structure function of different order. For an arbitrary reference order 2m this results in the 'ESS form' of (1.1), namely de Silva et al. (2017) and Krug et al. (2017) found that the linear scaling of (1.2) could not only be observed at relatively low Re τ well within the capabilities of current DNS, but also that the relative slopes D u p /D u m exhibit non-trivial, i.e. non-Gaussian, universal behaviour across various flow geometries such as flat plate boundary layers, pipe and channel flow and even Taylor-Couette flow. From our experience in Taylor-Couette flow  we learned that ESS-scaling according to (1.2) is not observed if large scale structures in the bulk (such as the Taylor rolls) contribute to the velocity component under investigation. This is also the case for the wall parallel velocity component in RB, as can be seen from the snapshots in figures 1(a,b). Since no equivalent scaling exists for the wall-normal component (Perry & Chong 1982), the 2D velocity field in RB is ruled out as a suitable candidate for this scaling analysis. Therefore the current study instead focuses on structure functions of the temperature field Here, ∆θ and Θ τ = −κ∂ z Θ| z=0 /U τ are analogous to ∆u and U τ , respectively. ESS-scaling of scalar structure functions in the energy containing range has not been considered yet. It is however conceivable that the situation will be similar to velocity structure functions: The theory underlying (1.1) and (1.2) is based on an inertial assumption, which implies that momentum transport is dominated by turbulent eddies that are larger than and therefore constrained by z η (where η is the Kolmogorov length scale) but smaller than δ. In the spirit of the Reynolds analogy and for P r ≈ 1 these eddies can be assumed to affect momentum and scalar similarly. Consequently, the relevant arguments based on the attached eddy hypothesis (see de Silva et al. 2015;Krug et al. 2017) can be expected to transfer to the scalar field as well. It is therefore our objective here to investigate whether and at what Ra the thermal boundary layers in RB convection exhibit an ESS scaling according to The expectation is that such a scaling regime should originate coinciding with the N uscaling transition, which has been link to the emergence of logarithmic boundary layers (Grossmann & Lohse 2011). Before moving ahead, we would like to mention that the attached eddy framework, from which (1.1) and (1.2) can be derived, is not new in the RB context. It has already been referred to by Ahlers et al. (2014), who investigate logarithmic dependencies of temperature profiles, and by He et al. (2014) in a study of f −1 temperature power spectra scalings the spatial equivalent to the k −1 scaling) in order to interpret the observations made.
Since the data set of Zhu et al. (2018) that will be employed for this endeavour is 2D, we will also check our results by using a three dimensional channel simulation as a reference and for comparison. A brief overview over both datasets will given in §2, before presenting our results in §3 and conclusions in §4.

Datasets
The simulations of Zhu et al. (2018) where performed using the second-order finitedifference code AFiD (van der Poel et al. 2015) on a two-dimensional domain with periodic boundary conditions in the lateral direction with Γ = 2 and P r = 1. While the original dataset spans six decades Ra ∈ [10 8 , 10 14 ], only four data points for Ra 10 11 are used here (see figure 1a). We refer to the original publication for further details on the numerical setup and validation of the results. One of the particularities of the 2D setup is that the large scale structures, which can be observed in the vector maps of figures 1(a,b), remain almost fixed in place allowing for simple temporal averaging. Since the temporal mean velocity gradient changes sign along the plate, we take the spatial mean of its absolute value when computing U τ in RB convection. For the present simulations we obtain L + /2 = LU τ /(2ν) ≈ [2400; 5700; 12300; 34400] at Ra = [10 11 ; 10 12 ; 10 13 ; 10 14 ], respectively.
A DNS of channel flow was performed using the fourth-order code described in Chung et al. (2014) at Re τ = hU τ /ν = 590, where h is half the channel height. The periodic (in streamwise and spanwise directions) box of size 12h × 4h × 2h was discretised by 640 × 320 × 240 grid points in the streamwise, spanwise, and wall-normal directions, respectively. A passive scalar with P r = 1 was added with values fixed to −0.5 at the bottom and 0.5 at the top wall as boundary conditions.
Convergence of the statistics computed from both datasets was checked by plotting the premultiplied probability density function of ∆θ at various r and found acceptable up to tenth order.

Direct analysis of the scaling in the energy-containing range in RB convection
We begin by plotting structure functions of second, fourth and tenth order from 2D RB convection as a function of r/z in figure 2. For all cases, we present results at three distances from the wall, namely z + = 30 (figure 2a-c), z + = 100 (figure 2d-f) and z + = 200 (figure 2g-i), where as usual the superscript + indicates normalization by ν/U τ . Clearly, there is no scaling according to a scalar equivalent of (1.1) for either Ra at the position closest to the wall ( figure 2a-c). However, in the other cases an approximately linear region appears for r/z > 1 making it tempting to fit the slopes directly. A fitting range 5z r 0.5L (indicated by crosses in the figure) captures this region quite well for all Ra and the resulting fits are shown as red dashed lines. No fits are computed at the lowest Ra where the fitting range becomes prohibitively small. applies to all panels. Red dashed lines represent fits in the range 5z r 0.5L, symbols mark r/z = 1 (circles) and the extent of the fitting range (crosses), respectively, for later reference.
A detailed comparison of the values of D θ p obtained in this way (figure 3a) reveals that the results depend on both Ra and z + in the investigated range. Generally, values of D θ p are higher at z + = 100 and this difference becomes larger with increasing Ra. In most cases, D θ p is very low and it is only at z + = 100 and Ra = 10 14 (and small p) that the values of D θ p are of comparable magnitude to results for D u p in high-Re turbulent boundary layers (TBLs). We point out that the decrease of D θ p with increasing p observed in some cases at higher p is unphysical and likely related to the insufficient fitting range at higher orders. We further emphasize that even at high Re a direct match between D θ p and D u p is not necessarily expected. It is well established (see e.g. Warhaft 2000, for passive scalars) that intermittency in the inertial range η r z is higher for scalars as compared to the velocity itself. This translates to lower scaling exponents ξ 2p in the corresponding scaling relation S θ 1 ∼ (r/z) ξ θ 2p /p as compared to the velocity counterpart ξ u 2p . From matching the inertial scaling with (1.1) at r = z de Silva et al. (2015) semi-empirically derived a linear relationship between D u p and ξ u 2p , suggesting that (disregarding other dependencies) these differences may persist in the energy-containing range investigated here. However, a definitive answer to this question will have to be based on high-Re data of the scalar field in wall-bounded turbulent flows.
While the direct analysis of the slopes D θ p remains inconclusive, additional insight can be gained from studying the second order structure function S θ 1 . In this case, the difference ∆S θ 1 = S θ 1 (r 2 ) − S θ 1 (r 1 ) between the structure function at two different separation distances r 1 and r 2 can be interpreted as the contribution of eddies with sizes in the range between r 1 and r 2 to the overall (scalar) energy. Consequently, when using the bounds of the log-linear scaling regime for r 1 , r 2 as indicated in figure 2d, this increment ∆S θ 1 , characterizes the contribution of the energy-containing range to the overall energy. Note that here we adopted the bounds of the fitting region used above for for simplicity, but other reasonable choices give qualitatively similar results. Figure  3(b), where ∆S θ 1 is normalized by the result obtained at Ra = 10 14 , shows that this quantity only increases mildly at low Ra. However, it rises steeply between Ra = 10 13 and Ra = 10 14 consistent with transitional behaviour in this range. No significant differences arise between z + = 100 and z + = 200 in this case.

Scalar ESS scaling in RB convection
We now focus on the relative scalings D θ p /D θ 1 according to (1.4). The ESS framework has been demonstrated to extend the scaling regime not only to low Re but also to a wider range of wall-normal distances. In particular, Krug et al. (2017) found convincing scaling for D u p /D u 1 as low as z + = 30. From figure 4a-c, it is clear that the same does not hold for the scalar structure functions in RB convection. Even at the highest Ra, there is no linear relationship between S θ 1 and S θ p for all orders considered at this location. The situation improves at z + = 100 ( figure 4d-f), where particularly at the highest Ra and with the y-axis cutoff fit at S 1 (r = 5z), symbols mark the locations of r/z = 1 (circles) and the fitting region indicated in figure 2 (crosses) for reference. low orders the curve begin to exhibit an approximately linear range. However, deviations become more apparent with increasing p and results at the higher orders in figures 4(e,f) demonstrate that the ESS scaling is not yet fully attained at this position. It is only at z + = 200 (figures 4g-i) that a convincing ESS scaling according to (1.4) is obtained up to tenth order. The scaling range is well established at Ra = 10 14 , already significantly decreases in size at Ra = 10 13 and is basically non-existent at Ra = 10 11 . This behaviour is in very good correspondence to the changes in the N u scaling in figure 1a, corroborating that the change in scaling observed there is indeed due to a transition in the BL structure from laminar (no scaling) to turbulent (with ESS scaling). It should be noted that, just as for velocity structure functions, ESS scaling at z + = 200 is observed for r/z 1 for all Ra, i.e. the spatial scaling range is the same. However, consistent with the results in figures 2 and 3(b), there is hardly any energy in this range at Ra below transition to the ultimate regime. Remarkably, the relative slopes attained for the temperature structure functions in 2D RB at z + = 200 appear to be the same as those measured for their velocity counterparts and reported in de . This is remarkable since as pointed out in §3.1, the directly measured slopes D θ p and D u p need not be the same. So this feature seems to be one of scalar fields.

Scalar ESS scaling in channel flow
The question remains why ESS scaling is only observed at a larger distance from the wall in the present case as compared to previous findings for velocity structure functions. At this point, possible explanations to be considered are that this might be either a consequence of the 2D setup, a property of RB convection or a feature of the scalar field. To address this, we present ESS results of scalar structure functions in turbulent channel flow in figure 5(a-c) at different orders. Indeed, the results are very similar to what was observed for 2D RB before in that there is no ESS scaling for z + = 30 and z + = 100. And again ESS scaling is recovered at z + = 200 with relative slopes matching those measured for velocity structure functions mirroring the observations made for RB convection.

Discussion and Conclusions
We have analysed the temperature boundary layers in the energy-containing range of 2D RB convection by means of temperature structure functions. Even though the RB structure functions exhibit log-linear scaling for r/z 1 when plotted against separation distance, the slopes remain small at low Ra and they were found to vary significantly with both Ra and z + , rendering the analysis inconclusive in this point. While a dependence of the slopes on Re τ and wall-normal position is also observed for velocity structure functions (see Krug et al. 2017), typically the log-linear scaling is much less evident in these cases compared to figure 2 at z + 100. Also for the scalar in channel flow investigated here (plots not shown) a direct scaling regime is not discernible such that it appears likely that the more prominent log-linear regimes in figure 2 are a consequence of the 2D setup. An important point remains however, that the contribution of the energycontaining range to the total energy increases significantly beyond Ra = 10 13 , which coincides with the transition in N u vs. Ra scaling.
The main finding of the paper is the clear evidence that the temperature structure function in the BLs of turbulent RB flows exhibit ESS scaling in the energy-containing range for large enough wall distances z + 100. The extent of the scaling range thereby reflects the behaviour of N u very well in increasing from non-existent at Ra = 10 11 (well in the classical regime) to considerable beyond the N u scaling transition at Ra = 10 14 . This provides further evidence that said transition is indeed related to a switch from laminar type to turbulent type thermal BLs.
Moreover, we establish that the relative slopes for scalar structure functions in the ESS-form are the same as those previously obtained when analysing velocity BLs. This is confirmed by comparing the RB results to those obtained in a planar channel geometry. In both cases, ESS scaling is only established at z + = 200 which is different from velocity structure functions and appears to be a feature of the scalar field. Interestingly enough, already Perry & Chong (1982) point out differences between the scalar and the velocity field. They argue that at the end of the 'lifetime' of an eddy, the vorticity contributions in the two rods of the assumed hairpin cancel such that no induced velocity field remains. However, such a cancellation does not apply for the scalar transported by the eddy, such that the 'debris', as they call it, of past eddies sets up a scalar background profile. If and how exactly this is related to the observations made here remains however unclear.
At least for the present cases, there also appears to be no difference between active (in RB flow) and passive (in the channel flow) scalars. In other regards, the deviations close to the wall underline in addition to the fact the observed values of D p /D m are sub-Gaussian ) that the ESS scaling and the values of the relative slopes are indeed non-trivial.
The work was financially supported by NWO-I and the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), both sponsored by the Netherlands Organization for Scientific Research (NWO). Part of the simulations were carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. We also acknowledge PRACE for awarding us access to Marconi based in Italy at CINECA under PRACE Project No. 2016143351 and the DECI resource Archer based in the United Kingdom at Edinburgh with support from the PRACE aisbl under Project No. 13DECI0246. We also acknowledge the support of the Australian Research Council and the McKenzie Fellowship program at the University of Melbourne.