Combined particle image velocimetry and thermometry of turbulent superstructures in thermal convection

Turbulent superstructures in horizontally extended three-dimensional Rayleigh-B\'enard convection flows are investigated in controlled laboratory experiments in water at Prandtl number $Pr = 7$. A Rayleigh-B\'enard cell with square cross-section, aspect ratio $\Gamma = l/h = 25$, side length $l$ and height $h$ is used. Three different Rayleigh numbers in the range $10^5<Ra<10^6$ are considered. The cell is accessible optically, such that thermochromic liquid crystals can be seeded as tracer particles to monitor simultaneously temperature and velocity fields in a large section of the horizontal mid-plane for long time periods of up to 6 h, corresponding to approximately $10^4$ convective free-fall time units. The joint application of stereoscopic particle image velocimetry and thermometry opens the possibility to assess the local convective heat flux fields in the bulk of the convection cell and thus to analyse the characteristic large-scale transport patterns in the flow. A direct comparison with existing direct numerical simulation data in the same parameter range of $Pr, Ra$ and $\Gamma$ reveals the same superstructure patterns and global turbulent heat transfer scaling $Nu(Ra)$. Slight quantitative differences can be traced back to violations of the isothermal boundary condition at the extended water-cooled glass plate at the top. The characteristic scales of the patterns fall into the same size range, but are systematically larger. It is confirmed experimentally that the superstructure patterns are an important backbone of the heat transfer. The present experiments enable, furthermore, the study of the gradual evolution of the large-scale patterns in time, which is challenging in simulations of large-aspect-ratio turbulent convection.


Introduction
The occurrence of temperature differences and the resulting buoyancy forces induce complex fluid motions in many geo-and astrophysical settings and technical systems (Kadanoff 2001;Verma 2018). This is one reason why thermally driven fluid flows have been investigated comprehensively in the past decades in theory, simulation, and laboratory experiment. Of central interest is a better understanding of the close link between structure formation and the resulting turbulent heat transport properties. In this context, the Rayleigh-Bénard model has established as one of the most frequently considered basic configurations (Bodenschatz et al. 2000;Ahlers et al. 2009;Chillà & Schumacher 2012). This simple convection model consists of a fluid layer between two parallel plates, which is uniformly heated from below and cooled from above, enclosed by adiabatic side walls. If the driving due to the vertical temperature gradient is large enough, turbulent Rayleigh-Bénard convection (RBC) develops in the layer that carries heat from the bottom to the top. Even though the model cannot represent the full complexity of most natural and technical systems, see e.g. Schumacher & Sreenivasan (2020) for the case of convection in the solar interior, it turned out to be very useful for studies of fundamental turbulence characteristics in thermal convection, which also motivates our present study.
In this work, turbulent superstructures of convection -a gradually evolving largescale order in horizontally extended layers with a coherence length larger than the layer height -are investigated in controlled laboratory experiments. The characteristic scale and dynamical evolution of superstructures is studied by a combination of stereoscopic particle image velocimetry and thermometry enabling a joint experimental determination of the local structure of the three velocity components u x , u y , and u z and the temperature T as fields in a larger horizontal cross sectional plane. This opens the way to an experimental investigation of the spatial structure of the local convective heat flux, its fluctuations, and thus also the global heat transfer across the layer. Our long-term experiments are conducted in water at a Prandtl number Pr = 7 in a cell with a square cross section and an aspect ratio Γ = l/h = 25, with l being the horizontal length and h the height of the fluid layer. Furthermore, we compare the experimental results with existing direct numerical simulation data for the same parameter range which will be detailed later. Characteristic length scales and times are found to be in the same range as the numerical results. We also compare the statistics of the local convective heat flux, the velocity and temperature fields, and discuss differences that can be traced back to resolution effects in the experiments.
The "large-scale order" in the form of superstructures in fully developed, horizontally extended convective turbulence has been analyzed with regard to the characteristic length scales for different Rayleigh and Prandtl numbers (Hartlep et al. 2003;Parodi et al. 2004;Pandey et al. 2018Pandey et al. , 2021Lenzi et al. 2021), the associated dependence of the turbulent heat transfer on the aspect ratio of the convection layer (Bailon-Cuba et al. 2010), and the effect of the aspect ratio on the characteristic pattern scale (von Hardenberg et al. 2008;Stevens et al. 2018). In the studies of Emran & Schumacher (2015) and Sakievich et al. (2016), it has been demonstrated that the large-scale structures slowly rearrange over time. The characteristic lengths and time scales of the superstructures, as determined by Pandey et al. (2018) on the basis of numerical simulations in the Eulerian frame of reference, have also been confirmed in a complementary analysis of Lagrangian trajectory clusters (Schneide et al. 2018;Vieweg et al. 2021b). Deep learning algorithms have been applied to parametrize the turbulent convective heat flux in global models by a reduction to a planar dynamical network assessing the amount of heat transported by superstructures in Fonda et al. (2019). Moreover, the coherence of the superstructures in the temperature and velocity field has recently been investigated in detail by Krug et al. (2020) and Blass et al. (2021), while the interplay between small-scale turbulent fluctuations and the large-scale turbulent superstructures has been studied by Green et al. (2020), Berghout et al. (2021), and Valori & Schumacher (2021).
All studies of turbulent superstructures cited above are based on direct numerical simulations of the Boussinesq equations. The number of experimental investigations of turbulent convection in large aspect ratio cells with Γ 1 is much smaller. Fitzjarrald (1976) measured co-spectra of the temperature and velocity in air at a Prandtl number of Pr = 0.7 with local sensor elements and determined a Rayleigh-number-dependence of the characteristic convection pattern scale. The large-scale structures in RBC have also been visualized in experiments with different working fluids at higher Prandtl and moderate Rayleigh numbers via the shadowgraph technique (Busse & Whitehead 1971, 1974Busse 1994) or photography (Krishnamurti & Howard 1981). Fluctuation profiles of the velocity and temperature across the layer were investigated by Adrian (1996). Turbulent superstructure patterns have been obtained in convection experiments in air by Kästner et al. (2018) and Cierpka et al. (2019) by time averaged velocity fields which were measured with particle image velocimetry, in short PIV (Raffel et al. 2018).
Our experiment has been set up to perform simultaneous measurements of the temperature and velocity field over long time intervals in a large section of the horizontal midplane of a Rayleigh-Bénard cell. We therefore suspend thermochromic liquid crystals (TLCs) in the flow that change their color with respect to the local temperature (Dabiri 2008;Moller et al. 2021). TLCs have also been used in RBC experiments by Zhou et al. (2007) to analyze the morphology of thermal plumes in a cylindrical cell with Γ = 1. However, these particles serve in the current case not only as temperature sensors, but also as tracer particles for stereoscopic PIV measurements and thus enable the simultaneous velocity and temperature measurements, see also Schmeling et al. (2014). On the basis of these measurements, the strong influence of the turbulent superstructures on the local convective heat flux is demonstrated. In addition, their horizontal extent as well as their reorganization over longer time intervals is analyzed. In order to assess both the experimental and numerical investigations of turbulent superstructures in RBC, the results of the measurements are compared with those of direct numerical simulations (DNS) performed previously by Pandey et al. (2018) and Fonda et al. (2019).
The outline of the article is as follows. In § 2, we present experimental details and give a short compact overview of the numerical simulations, which are used for comparison. § 3 follows with a characterization of the turbulent superstructures. § 4 presents the analysis of the local convective heat flux, which is followed by § 5 and § 6 on the characteristic length scales and long-term evolution of the large-scale patterns, respectively. The manuscript ends with a conclusion and an outlook in § 7.

Methods
The turbulent Rayleigh-Bénard flow is determined by the geometry of the flow domain, the working fluid and the strength of the thermal driving. These dependencies are quantified by dimensionless numbers, the already mentioned aspect ratio Γ = l/h, the Prandtl number Pr = ν/κ and the Rayleigh number Ra = αg∆T h 3 /(νκ). Here, ∆T is the temperature difference between the hot and cold plates at the bottom and top, respectively, g is the acceleration due to gravity, and α, ν and κ represent the thermal expansion coefficient, the kinematic viscosity and the thermal diffusivity of the working fluid, respectively. Since the aspect ratio, the Prandtl number and the Rayleigh number determine the flow, these dimensionless numbers will be similar to those of the numerical simulations for comparison.

Experimental investigation
For the experimental investigation of the turbulent superstructures in RBC a cuboid cell with dimensions of l × w × h = 700 mm × 700 mm × 28 mm, thus having an aspect ratio of Γ = 25, has been built. The experimental setup is explained in detail in Moller et al. (2021). A sketch of the experiment can be seen in figure 1.
While the heating plate at the bottom of the cell is made of aluminum, the sidewall and the cooling plate at the top are made of glass. As can be seen in figure 1, the temperature of the cooling plate is adjusted by an external cooling water circuit, which is covered with another glass plate at the top side. Since the working fluid in the cell is water for all the experiments in the present work, the entire flow domain is optically accessible, such that optical measuring techniques can be used to study the Rayleigh-Bénard flow. In our case, thermochromic liquid crystals (TLCs) are inserted into the flow as tracer particles to measure quantitatively the temperature distribution in the fluid. When illuminated with white light, these particles change their color appearance in dependence on the temperature and can therefore be used to measure the local temperature in the flow. Physical details of this mechanism can be found in Dabiri (2008). The TLCs are illuminated in the central horizontal bulk region through the transparent sidewall with a thin sheet of white light, which was collimated using Fresnel lenses with a minimal thickness between 3 mm and 4 mm over the entire field of view. More details on the design of the white light source and the light sheet optics can be found in Moller et al. (2020). For additional background information, the reader is also referred to Schmeling et al. (2014) or Anders et al. (2020).
The color appearance of the TLCs is recorded from the top with a color camera through the glass plates and the cooling water. Based on the hue value H, which represents the pure color shade in the HSV -color space (H -Hue, S -Saturation, V -Value), the temperature is determined from local hue value calibration curves via linear interpolation (Moller et al. 2021). Here, TLCs of type R20C20W (LCR Hallcrest Ltd) are used. According to the nominal specifications, which are only achieved for the same direction of observation and illumination, the TLCs start to get red at 20 • C and their color shade varies continuously across the visible wavelength spectrum until they appear in blue at 40 • C. However, when analyzing the color of the TLCs dispersed in a fluid and illuminated with a thin light sheet, the observation angle between the color camera and the light sheet ϕ cc must be varied towards an almost perpendicular arrangement. Since the applicable temperature range, in which the color shade changes from red to blue, shifts and decreases with increasing observation angles (Moller et al. 2019(Moller et al. , 2020, the angle must be adapted carefully. In this case, the angle has been adjusted to ϕ cc = 65 • as a final tradeoff, which allows to measure the temperature in the temperature range 18.9 • C T 20.8 • C with a mean uncertainty of less than 0.1 K (Moller et al. 2021).
Besides the color camera for the determination of the temperature field, two monochrome cameras in a stereoscopic arrangement are applied to record images of the TLCs in the horizontal mid-plane, such that not only the horizontal velocity components, but also the vertical velocity component can be measured by means of the stereoscopic particle image velocimetry (Prasad 2000;Raffel et al. 2018). The main motivation for using the stereoscopic setup of the cameras is that the simultaneous measurement of the vertical velocity component and the temperature allows to estimate the local convective heat flux in RBC, as it will be shown in the following sections. In order to assess the effects of the Rayleigh number on the turbulent superstructures, three distinct temperature differences between the heating and cooling plates are prescribed, yielding the Rayleigh numbers which are listed in table 1 together with further experimental details. This experimental setup is appropriate to conduct longterm investigations of RBC -a big advantage of experimental studies in comparison to direct numerical simulations, since it is not necessary to compute the current state of the flow with a high temporal and spatial resolution. As discussed in Bailon-Cuba et al. (2010), the numerical effort grows in addition with Γ 2 . In the present case this allows to characterize the reorganization of the slowly evolving turbulent superstructures without determining the temperature and velocity field continuously, meaning that temporarily observing the flow with the cameras is sufficient to capture this gradual dynamics.
For each of the experiments, a strict procedure is followed. Once the temperatures of the heating and cooling plates reach the stationary state, the TLCs are inserted into the working fluid. As the flow in the cell is disturbed by the seeding process, we wait for 45 min, so that the perturbation has decayed prior to the actual measurements. Afterwards, image series with a duration of 5 min are recorded with a frequency of f = 5 Hz every t rec = 20 min. The last image series is recorded after six hours of continuous operation, resulting in 19 series with 1500 images for each Rayleigh number. In order to keep the computing time of the evaluation in limits, the temperature and velocity fields are determined for each fifth image. However, the images in between are necessary as well, because for the measurements at the two highest Rayleigh numbers the temporal delay of ∆t = 1 s between each fifth image is too large to properly compute the velocity field via PIV. On the contrary, the displacement of the TLCs on the cameras' images during the time gap of ∆t = 1 s has turned out to be well suited for the determination of the velocity field at the smallest Rayleigh number. In addition, the temperature of both, the heating and the cooling plate was monitored by four Pt-100 temperature sensors in each plate during the whole experimental run to allow for a precise temperature control by the thermostats. The sensors that measure the heating plate temperature are mounted inside the aluminum heating plate close to the surface. Due to the heating plate's high thermal conductivity, this arrangement allows for a precise measurement of the heating plate temperature while not disturbing the flow. The temperature sensors at the cooling plate are glued on the glass surface that confines the fluid from the top. Thereby, the temperature of the cooling plate is measured directly at the solid-liquid interface.
The small temperature range of the TLCs and the intended large aspect ratio of the cell limit the range of accessible Rayleigh number Ra when turbulent convection under Boussinesq conditions is studied, see table 1.  Table 1: Parameter sets for the measurements in the Rayeigh-Bénard cell with Γ = 25 and water as the working fluid. The table lists the Rayleigh number Ra, the temperature of the heating plate T h and cooling plate T c , the non-dimensionalized size of the intersection of the field of view (fov) of the three camerasx fov ×ỹ fov = x fov /h × y fov /h, the number of interrogation windows N x × N y with an overlap of 50%, the resulting grid resolution ∆(x,ỹ), the free-fall time unit t f , the non-dimensionalized total measurement timeτ t = t total /t f after the initial transient period, and the number of the evaluated snapshots N s . It is noted that the Rayleigh numbers in the first column are always given as approximate values in the discussion of the results.

Direct numerical simulations for comparison
In fluid dynamics, customarily non-dimensional or dimensionless equations are studied. Therefore, we simulate the following equations, which represent the conservation laws of mass, momentum and energy in dimensionless form where all fields are indicated by an additional tilde (Chillà & Schumacher 2012;Verma 2018), In these equationsũ = (ũ x ,ũ y ,ũ z ),p andT denote the dimensionless velocity, pressure and temperature fields, respectively. The equations are made dimensionless using the following scales: the height of the flow domain h as the length scale, the temperature difference between the horizontal plates ∆T as the temperature scale, the free-fall velocity u f = √ αg∆T h as the velocity scale and the free-fall time t f = h/u f as the time scale. For our incompressible flow, the pressure field is computed using the velocity field by solving a Poisson equation. Therefore, the characteristic unit of the pressure field is related to that of the velocity and reads ρ 0 u 2 f , with ρ 0 being the constant mass density of the fluid. All scales used for non-dimensionalization are given in table 2.
We integrate the equations (2.1) -(2.3) by applying a spectral element solver Nek5000 (Fischer 1997) in a cuboid box with square cross-section of size 25h × 25h. The velocity field satisfies the no-slip condition at all boundaries, while the temperature field satisfies the isothermal condition at the horizontal walls and the adiabatic (or thermally insulated) boundary condition at the vertical side walls. We decompose the flow domain into N e spectral elements which are again small cuboids. They form a non-uniform spectral element mesh with smaller element spacing towards all faces in order to resolve the strong variations of the velocity and temperature in the near-wall region correctly. The turbulence fields within each element are discretized using Lagrangian interpolation polynomials of order N in each space direction using Gauss-Lobatto-Legendre collocation

Quantity
Characteristic unit Length points. Thus, the flow domain consists of N e N 3 mesh cells of non-uniform size, see also Scheel et al. (2013).
With a view to the comparability of the numerical and experimental study of turbulent superstructures we use direct numerical simulation data of RBC in a flow domain with Γ = 25 and Pr = 7. Furthermore, the Rayleigh number is adjusted to Ra = 10 5 and Ra = 10 6 in the simulations, because these Rayleigh numbers span the range of the experimental investigation. The simulations are started from the motionless conduction state with random perturbations and we wait until the initial transients have decayed and the flow has reached a statistically steady state. After that, we continue our simulations for a total of τ total = 2670 free-fall time units for Ra = 10 5 and τ total = 3033 free-fall time units for Ra = 10 6 . Finally, 268 and 560 equidistant snapshots were stored during this time interval for Ra = 10 5 and Ra = 10 6 , respectively. Ensuring the adequacy of the spatial and temporal resolutions is important in order to faithfully study a turbulent flow using DNS. In a high-Pr convective flow, the smallest length scale in the temperature field is finer compared to that in the velocity field (Silano et al. 2010;Pandey et al. 2014). Therefore, we ensure that the average size of the mesh cells remains smaller or comparable to the Batchelor length scale within each horizontal plane (Scheel et al. 2013).

Characterization of turbulent superstructures
Exemplary instantaneous temperature and vertical velocity component fields are shown in figure 2. The instantaneous fields in figures 2a and 2b are obtained from measurements in the horizontal midplane at Ra = 7 × 10 5 , while those in figures 2c and 2d result from the DNS at Ra = 10 6 at the same position. With respect to the measurement uncertainty the field of view of the three cameras applied for the measurements has been adjusted, yielding the intersection of the experimentally determined temperature and velocity fields in figure 2. For a better comparability the numerical results are restricted to the same section of the horizontal midplane.
As it can be seen in figures 2b and 2d, the Rayleigh-Bénard flow is characterized by strong spatial variations ofũ z with up-and downwellings due to rising and falling thermal plumes. The fingerprints of the flow can to some extent also be identified in the 10 15 x Figure 2: Visualization of the instantaneous (a-d) and the corresponding time-averaged fields (e-h) of the temperatureT and of the vertical velocity componentũ z in the mid-plane, obtained from the laboratory experiments (exp) for the Rayleigh number Ra = 7 × 10 5 (a,b,e,f) and from the direct numerical simulations (sim) for the Rayleigh number Ra = 10 6 (c,d,g,h). The circles within the figures 2e -2h indicate local events with a strongly pronounced vertical velocity component despite a small deviation from the average temperature as described in the text.
corresponding instantaneous temperature fields. However, one observes that larger areas with similar temperature ranges appear. A closer look reveals that these areas can also be detected in the instantaneous field of the vertical velocity component. This indicates that the flow does not only exhibit structures at smaller length scales, but also tends to develop patterns at somewhat larger scales. As the small-scale flow structures rearrange quickly and display stronger temporal fluctuations, they can be effectively removed in both fields by time-averaging. As a result, the large-scale patterns are revealed. The latter are denoted as turbulent superstructures and can be seen clearly in the time-averaged fields of the temperature and of the vertical velocity component in the lower row of figure 2.
In the evaluation of the experimental data, the full recording time of 5 min for each series of images has been taken as the time-averaging interval which translates into 125 free-fall units in the considered case with Ra = 7 × 10 5 . The numerical data have been evaluated according to Fonda et al. (2019), in which averaging times of 207 free-fall units and 179 free-fall units are given for Ra = 10 5 and Ra = 10 6 , respectively. Hence, the experimental data are averaged over a somewhat shorter time interval. However, the quantitative analysis on the basis of the time-averaged fields in § 5 and § 6 will not be affected by this slight difference, since the averaging time is sufficiently long to clearly reveal the turbulent superstructures, see also the supplementary information of Pandey et al. (2018).
The comparison of the time-averaged temperature and velocity fields in figure 2 shows that the turbulent superstructures yield sharp contours in the vertical velocity component, while giving smoother contours of the temperature field. In addition, it can be anticipated without further analysis that the characteristic length scale in the temperature field is somewhat larger than in the velocity field. This behaviour was also found in the DNS by Pandey et al. (2018). The effect can be explained by the fact that at the onset of convection, the temperature field and the vertical velocity field are perfectly synchronized (hot fluid dwells in upward direction whereas cold fluid sinks down) and thus show the same characteristic wavelengths. The synchronization decays with increasing Rayleigh number since the temperature field is not only advected by the vertical velocity component across the midplane, but also stirred by increasing velocity fluctuations in the horizontal directions. As a consequence the contours of the superstructures in the temperature field might appear broader given the fact that the thermal plumes are dispersed by the velocity fluctuations when rising or falling into the bulk of the layer. Locally this can lead to differences in the pattern scales of temperature and velocity fields. For example, this is visible around (x,ỹ) = (6, 7) in the experiment and around (x,ỹ) = (15, 16) in the simulation, where the temperature is close toT = 0.5 despite a strong vertical fluid motion in upward or downward direction (see also circles in figure 2). However, the turbulent superstructures in the temperature and velocity fields can be found in good qualitative agreement. The visual inspection of the results indicates that the experimentally determined turbulent superstructures seem to be somewhat larger, which will be detailed in § 5.
Comparing the temperature fields in figure 2, it can also be seen that the occurring temperatures span a wider range in the experiment. This is confirmed by the standard deviation of the temperature σ T in table 3. We suspect that this is an effect of the turbulent superstructures on the temperature distribution at the heating and cooling plate in the experimental setup, as the impinging superstructures yield some slightly warmer and colder regions on both plates. In contrast, the range of vertical velocity amplitudes is larger in the simulations as obvious from figure 2. A possible explanation for this observation is the higher spatial resolution of the DNS data. The in-plane resolution for the PIV measurements can be largely increased by advanced algorithms. For the current case we obtained an in-plane vector spacing of about ∆(x,ỹ) = 0.084, which corresponds to 2.35 mm. In the vertical direction, the resolution is given by the light sheet thickness, amounts to 3 to 4 mm and cannot further be improved. It is this vertical resolution that limits a better resolution in all spatial directions. Even though the experiment allows to resolve most of the small-scale structures, the peak velocities embedded in the fine structures are not fully captured due to the spatial averaging of the velocity in the particle image velocimetry measurements, such that the largest magnitudes are slightly smoothed out (Kähler et al. 2012). However, the difference between the numerical and experimental results regarding the velocity range is smaller as seen in table 3 and in figure 6, where the root-mean square velocity U rms = ( u 2 x + u 2 y + u 2 z ) 1/2 and the resulting Reynolds number Re = (Ra/Pr) 1/2 U rms are also listed for each simulation and experiment.
The root-mean-square velocities and Reynolds numbers of the DNS have also been computed on a coarser grid to investigate the effect of the smaller spatial resolution of the measurements. For this, the temperature and velocity data have been averaged in small interrogation volumes around the horizontal mid-plane. Hence, the numerical data have not only been averaged in the horizontal direction, but also to some extent in the vertical   .1), the mean characteristic wavelength of the turbulent superstructures determined from the temperature field λ T and from the field of the vertical velocity component λ uz . In case of the DNS, the quantities are determined from the horizontal mid-plane for a better comparability with the experimental data. In addition, the temperature and velocity data of the DNS are recomputed in small interrogation volumes around the horizontal midplane to estimate the effect of the smaller spatial resolution of the experimental data. The quantities resulting from the coarse-grained grid are denoted as σ Tco , U co, rms , Re co and Nu co .
direction, thereby taking into account that the light sheet for the illumination of the measurement area has a certain thickness. The size of the interrogation volumes has been chosen such that the spatial resolution of the numerical data for Ra = 10 5 and Ra = 10 6 roughly matches that of the experimental data for Ra = 2 × 10 5 and Ra = 7 × 10 5 , respectively. The resulting root-mean-square velocities U co, rms and Reynolds numbers Re co are also given in table 3. As expected, the magnitude of these quantities decreases due to the smaller spatial resolution. The numerical results are therefore fairly well in line with the experimental results.

Analysis of the local convective heat flux
The simultaneous measurement of the temperature and the vertical velocity component allows to determine the local convective heat flux in RBC. Measurements of the local convective heat flux have been performed in the past, e.g., by tracking a mobile Lagrangian temperature sensor (Gasteuil et al. 2007). Shang et al. (2004) combined laser-Doppler velocimetry (LDV) and temperature probe measurements to determine time series of the components of the local convective heat flux vector as a function of the cell radius. It was found that the most significant flux contribution stems from the large-scale circulation that rises or sinks along the sidewalls. The novelty of the present study is that the local heat transfer is investigated as a field in the mid-plane in the Eulerian frame of reference.
The Nusselt number Nu as the global measure of turbulent heat transfer is defined as Recall that Nu(z) takes the same value for any 0 z 1, even though the contributions to the heat transfer due to convective fluid motion (conv) and diffusion (diff) differ near the boundaries compared to those in the bulk, see e.g. Scheel & Schumacher (2014). At the horizontal mid-plane, the second term in eq. (4.1) is zero in the Boussinesq regime with its top-down symmetry. In dimensionless units, the two corresponding local currents are given by Since 0 T 1, we take a symmetric form which is obtained by decomposing the original temperature field into the linear diffusive equilibrium profile and the rest, (4.4) Thus, we define eventually a local convective Nusselt number atz = 1/2 by which is directly accessible as a field in the experiments. The uncertainty for the temperature measurements of 0.1 K results into an uncertainty for the individual values of the Nusselt number of about 14% for Ra = 2 × 10 5 , 7% for Ra = 4 × 10 5 and 4% for Ra = 7 × 10 5 . However, the mean values (see table 3) show a much lower uncertainty as ∼ 38000 data points for the two lower Rayleigh numbers and ∼ 85000 data points for the highest Rayleigh number are averaged. An exemplary instantaneous field of the local convective Nusselt number, which has been determined from the figures 3a and 3b, is depicted in figure 3c. These results of the measurement at Ra = 2 × 10 5 already indicate the effect of the superstructures on the local convective heat flux which appears even clearer in the corresponding time-averaged fields in the figures 3d, 3e, and 3f. It can be seen that the maximum values of the local convective Nusselt number occur in the regions, where the time-averaged temperature reaches its extreme values, corresponding to updrafts and downdrafts of the turbulent superstructures. It is therefore confirmed that the superstructures strongly enhance the heat flux from the bottom to the top side of the Rayleigh-Bénard cell. It can also be seen that multiple regions of strong heat transfer develop that are distributed homogeneously over the whole field of view. These regions are not constrained by the side walls which is different to Rayleigh-Bénard convection experiments in small aspect ratio cells (Shang et al. 2004).
The quantitative analysis of the local heat flux is done by the probability density functions (PDFs) of the local convective Nusselt number as shown in figure 4a for each of the three measurements and in figure 4b for the two DNS. In this case, the PDFs incorporate all the temporally resolved values obtained from the measurements and the DNS at each Rayleigh number. As expected, the PDFs are more extended for positive values, since the total heat transfer increases due to convective motion. It can also be seen that the negative tails increase with the Rayleigh number, which shows that in turbulent thermal convection the upward motion of cold fluid and the downward motion of warm fluid become more probable. The data also allows to evaluate the probability density functions of the normalized vertical convective heat current component j z /j z,rms as shown in figure 5a for each of the three measurements and in figure 5b for the two DNS. Shang et al. (2004) performed the same analysis and showed for a Γ = 1 cell that the PDFs fall on top of each other for a Rayleigh number range from Ra = 1.8 × 10 9 to 7.6 × 10 9 . The present PDFs of the normalized vertical convective heat current are supported on a similar interval (see their figure 9). The PDFs are also skewed to positive values since a mean heat flux goes from the bottom to the top. It can also be seen now that the positive tails on both panels collapse fairly well onto each other, both for experiments and simulations. However, the negative tails still differ even though being closer compared to figure 4. This effect is again attributed to the non-ideal boundary conditions that hinder the fast heat transport by conduction in the top plate. This effect becomes more dominant with increasing Rayleigh numbers as will be outlined in the following.
In order to further compare the experimental and numerical results regarding the heat flux, the spatial and temporal average of the local convective Nusselt number Nu loc are plotted in figure 6. Even though the contribution of the diffusive heat flux cannot be obtained from the experimental data, the resulting average of the convective flux with respect to time and area is representative for the total heat flux, since the observation area in the experiment is large enough. By applying the power law fit of Nu = 0.133 × Ra 0.298   An underestimated value in the experiment is possible, because the temperature and the velocity are measured with a lower spatial resolution and can also be considered as an average across the thickness of the light sheet used for the illumination of the tracer particles. The effect of these two aspects can be shown easily by means of the data of the numerical simulations. For this, the temperature and the vertical velocity component are averaged in small interrogation volumes around the horizontal midplane as it has been described in § 3 for the estimation of the root-mean-square velocity U co, rms and the Reynolds number Re co . If the coarse-grained temperature and the vertical velocity fields are applied to compute the global Nusselt number Nu co from the DNS data, the values converge to the experimentally determined Nusselt number Nu according to the results in table 3.
However, as the Nusselt numbers show, the convective heat flux in the experiment is still smaller in comparison to the simulations. This effect may be attributed to the boundary conditions. In the numerical simulations, perfect isothermal boundary conditions can be set, while temperature inhomogeneities at the boundaries cannot be fully impeded in an experiment. A dimensionless measure of the uniformity of the prescribed temperature at the boundaries is the Biot number Schindler et al. 2022), which is given by with k f , h and k w , h w being the heat conduction coefficient and the thickness of the fluid layer and of the wall material, respectively. The Biot number can also be interpreted as the ratio of the thermal resistance due to thermal conduction in the wall and the thermal resistance due to convective heat transfer at the wall. If Bi 1, then the limiting heat transfer mechanism is the convective heat transfer and not the conduction within the wall material, such that the isothermal boundary condition can be assumed. For the heating plate made of aluminum the Biot number is of the order of Bi ∼ 10 −3 and thus sufficiently small. However, the cooling plate is made of glass for the optical access. The glass plate has a thickness of h w = 8 mm and a thermal conductivity of k w = 0.78 W/mK. Assuming k f = 0.6 W/mK for the thermal conductivity of water at the mean fluid temperature of about T = 20 • C, the Biot number increases with the Rayleigh numbers adjusted in the experiment from Bi = 0.86 to Bi = 1.22. Therefore, the thermal resistance due to conduction in the wall is not negligible in this case and limits the heat transfer. A transparent material with a much higher thermal conductivity would be sapphire. However, such a sapphire plate of 70 × 70 cm 2 is not available off the shelf and would have to be produced individually as a growing aluminum-oxide monocrystal. Therefore, we had to stick with the non-ideal boundary conditions. Hot plumes from the bottom cannot transfer their thermal energy to the top plate in the same way as for the case of an isothermal boundary condition in the simulation. Due to continuity, the flow is forced to move downward with higher temperatures. This results in more events with a negative local convective Nusselt number as obvious from figure 4a, where the correspondence between experiment and simulation is good for positive values of the local convective Nusselt Number, but differs for negative ones. The occurrence of events transferring heat towards the heating plate can clearly be seen in figure 3c in the form of blue regions, which correspond to hot fluid streaming downward, for example at (x,ỹ) ≈ (14, 16), and cold fluid that moves upward, for example at (x,ỹ) ≈ (4, 9), indicated by the circles in the figures 3a -3c.
Due to the limitation of the convective heat transfer, which result from the material properties of the top plate, the local convective Nusselt number in the experiment is on average smaller than in the simulation, which is also visible by the larger positive tails for the DNS data in figure 4. The difference for the mean value increases with the Rayleigh number, as the resistance for heat conduction in the top plate remains constant and the turbulent convective heat transfer increases and thus the Biot number increases as well. In the present case, the ratio of the Nusselt number obtained from the DNS and from the experiment is in the same order as the Biot number, i.e. Nu/Nu exp ∼ Bi. This effect would also explain the smaller velocities in the experiment, as the driving of the flow is hampered. Moreover, Vieweg et al. (2021a) found that the variance of the temperature in the mid-plane is two orders of magnitude larger in the case of constant heat flux compared to the case of constant temperature at the horizontal boundaries †. This finding is in line with the larger standard deviation of the temperature σ T for the experimental data according to table 3.

Characteristic length scale of turbulent superstructures
The visual comparison of the temperature and velocity fields in figure 2 already indicates that turbulent superstructures in the experiments have a somewhat larger horizontal length scale than in the numerical simulations. For the quantitative estimation of the size of the superstructures in the experiment, the temperature field is considered in the following. As described in Pandey et al. (2018), the size of the superstructures can be determined via Fourier spectral analysis on the basis of the fieldΘ, which represents the temperature difference to the linear diffusive equilibrium profile at the height of the measurement plane. For simplicity, this field is denoted as the temperature field from here on. The characteristic horizontal length scale of the superstructures can be derived from the so-called power spectrum, which quantitatively represents the match between the temperature field and sinusoidal plane waves with different wavelengths and orientations.
A power spectrum of an exemplary time-averaged temperature field depicted in figure 7a can be seen in figure 7b. Here, the power spectrum PΘ is arranged in the way, such that the dimensionless wavenumberk = (k x ,k y ) increases from the center towards each of the edges. Hence, the entries of the power spectrum on the same circumference around the center represent features with the same wavenumber but a varying orientation. It will also be noted that the power spectrum is symmetric around the center, such that PΘ (k x , k y ) = PΘ (−k x , −k y ). Due to the time-averaging, the turbulent superstructures represent the prominent pattern in the temperature field, which is also obvious from the two distinctive maxima of the power spectrum. These two symmetric maxima occur, since the superstructures are in this case preferentially arranged along a certain direction, namely along the line connecting the two maxima. However, not for all data records the superstructures exhibit a preferential orientation along a specific direction. In some cases, they are also represented by a ring-like structure of enhanced spectral power rather than two prominent maxima. Therefore, the azimuthal average of the power spectrum can be taken to determine the circumference with the largest average, in order to identify the characteristic wavelength of the superstructures. For a better visualization of the length scales the average intensities I λ normalized with the maximum are shown in figure 7c in dependence on the dimensionless wavelengthλ = λ/h. The latter is related to the wavenumber viaλ = 2 π/k. As the gray curve shows, there are two distinctive local maxima. While the peak at λ ≈ 8.7 represents the superstructures, the peak atλ ≈ 29 results from the experimental boundary conditions. The second peak corresponds to the vertically aligned maxima close to the center in figure 7b and can be explained by the increasing temperature of the cooling water in its flow direction from the bottom to the top edge in figure 7a. A close look to this figure reveals that the measured temperature fields also exhibit the trend of larger temperatures along the cooling water flow in y -direction due to the slight increase of the cooling water temperature, which amounts to ∆T cw,y 0.15 K. In order to ensure that the characteristic wavelength of the superstructures is not determined based on the wrong maximum in any case, the large-scale trend is removed from the time-averaged temperature fields by computing the inverse Fourier Transform without considering the smallest wavenumbers close to the center of the power spectrum, meaning that the corresponding Fourier coefficients are set to zero. The limit of the wavenumber of k = (k 2 x +k 2 y ) (1/2) = 0.4 has been chosen, such that all the wavelengthsλ > 2π/0.4 = 15.7 are not taken into account. The latter span much wider horizontal dimensions compared to the turbulent superstructures. On the basis of the temperature field with the removed large-scale trend, which is depicted in figure 8a, the procedure for the determination of the average intensities I λ is then performed again. The resulting wavelength intensities with a clear peak at the characteristic wavelength of the superstructures can be also seen in figure 7c. This confirms the applicability of the filtering approach.
In the evaluation of the numerical data such a filtering is not required, because there are no deviations from the isothermal boundary condition yielding a large-scale trend of the temperature. Therefore, the power spectrum obtained from the time-averaged temperature field can directly be applied to derive the characteristic wavelength of the superstructures from the position of the maximum azimuthal average. The mean characteristic wavelength of the superstructures λ T obtained from the experimental and numerical investigations are given in table 3. At each Rayleigh number the 19 timeaveraged temperature fields resulting from the image series covering 5 min have been considered to determineλ T for the experiment, while the sliding time averages in the total runtime of the simulations have been used to calculateλ T for the DNS. The total runtimes and the number of equidistantly stored sliding time averages are given in Fonda et al. (2019).
For the estimation of the characteristic wavelength of the superstructures the timeaveraged field of the vertical velocity component has been considered, too. In principle the resulting wavelengths λ uz according to table 3 are obtained in the same way as for the temperature field, but the filtering approach of the experimentally determined velocity field has been changed slightly. It is obvious from figure 3e that some regular variations of the vertical velocity component remain on small length scales after the time-averaging. Thus the filtering has been modified, such that variations of the velocity on either much larger or smaller length scales compared to the turbulent superstructures are removed. The wavelengths in the range 4.8 λ uz 15.7 are not affected by the filtering. However, despite the additional filtering of small-scale variations of the velocity, the characteristic wavelengths determined from the time-averaged velocity field are somewhat smaller than those of the temperature field according to table 3. This is also the case for the fields obtained from the DNS, which might be caused due to the time averaging.
The results in table 3 show that the characteristic wavelength of the superstructures λ T increases with the Rayleigh number according to the experiments, while the wavelength nearly remains constant according to the DNS. Moreover, the characteristic wavelengths determined from the DNS are typically one third smaller. Considering the boundary conditions of the experiment this is also reasonable. Due to a thick insulation around the sidewall of the cell and the adaption of the temperature in the lab to the mean temperature of the working fluid the adiabatic boundary condition at the sidewall is well satisfied, but deviations from the isothermal boundary condition at the horizontal plates cannot be fully impeded in an experimental setup. Despite extensive efforts, in the experimental setup at hand especially the top plate made of glass for the optical accessibility cannot compensate the effect of the impinging flow structures without a temporal delay.
Based on the above discussion of the Biot number and on theoretical considerations (Hurle et al. 1967) it can be shown that the limited thermal conductivity of the plates enclosing the fluid layer strongly alters the flow in comparison to the classical RBC with isothermal boundary conditions. Essential is the ratio of the effective thermal conductivity of the working fluid k f to the thermal conductivity of the solid plates k w . The two extreme cases k f /k w → 0 and k f /k w → ∞, which correspond to the Dirichlet boundary condition with constant temperature at the plates and the Neumann boundary condition with constant heat flux at the plates, have been compared in the study from Vieweg et al. (2021a) on the basis of DNS for Γ = 60. While the first case yields the turbulent superstructures with a horizontal extent larger than the height of the cell, the second case results in structures, which gradually grow and aggregate to one large structure spanning even the entire flow domain. None of both extreme cases is matched perfectly in any experiment, however, here the temperature at the bottom plate is almost uniform due to the high thermal conductivity of aluminum and k f /k w 1. At the top plate which is made of glass, this ratio is k f /k w ≈ 1 when taking into account that the effective thermal conductivity of the fluid increases due to the convective motion. Therefore, the increase of the horizontal characteristic scale of the superstructures in comparison to the DNS might be caused by non-ideal boundary conditions at the cooling plate which are a mixture of constant temperature and constant heat flux.

Long-term evolution of the turbulent superstructures
Compared to the data of the DNS, which are often limited to total simulation times of several hundred free-fall units, the experimental data span approximately 5 × 10 3 up to 10 4 t f free-fall units for the different Rayleigh numbers. Hence, the ability to run the experiments over long time intervals opens the possibility to evaluate the reorganisation of the turbulent superstructures. The reorganization of the turbulent superstructures can be seen in figure 8, which shows two different time-averaged temperature fields in the horizontal mid-plane for the different Rayleigh numbers adjusted in the experiments. The time difference between these two temperature fields is 200 min, which corresponds to 2663 t f at Ra = 2 × 10 5 , to 3846 t f at Ra = 4 × 10 5 , and to 4988 t f at Ra = 7 × 10 5 . The differences in the distribution of warm and cold regions in the temperature field demonstrate the relocation of the superstructures.
On the basis of the numerical studies from Pandey et al. (2018) and Fonda et al. (2019) it can be estimated that the turbulent superstructures in water reorganize every 200 freefall units on average, if the Rayleigh number is in the range 10 5 Ra 10 6 . In order to assess the gradual long-term evolution of the superstructures, for each time-averaged field obtained from one of the 19 series of 1500 snapshots the central coefficient of the normalized two-dimensional cross-correlation is calculated according to Here, the running index i ∈ [0, 18] andt rec = t rec /t f is the dimensionless time interval between each of the 19 image series with t rec = 20 min as given in subsection 2.1. The Greek symbol α represents the quantities α = {T ,ũ z , Nu loc },X α the respective time-averaged, dimensionless fields, µX α and σX α the corresponding mean value and the standard deviation, respectively. Furthermore, · denotes the average. This analysis shows that for all Rayleigh numbers the temporal variation of the time-averaged temperature fields is not as pronounced as the variation of the time-averaged fields of the vertical velocity component and the local Nusselt number. If two successive time-averaged temperature fields are considered, respectively, the correlation coefficient remains on a high level of RΘΘ ≈ {0.5, 0.7, 0.8} for Ra = {2 × 10 5 , 4 × 10 5 , 7 × 10 5 }. This behaviour can again be attributed to the boundary conditions, as the deviation from the isothermal boundary condition at the top wall becomes stronger with increasing Rayleigh number. If the temperature value at the top plate is at a certain point larger than in the neighbourhood due to an impinging superstructure, this causes a stabilization of the superstructure at that position, such that a random reorganisation of the structures will take longer or is inhibited. In order to analyze the gradual reorganisation, the time-averaged fields of the temperature, vertical velocity component and local Nusselt number for Ra = 2 × 10 5 are evaluated, as the dimensionless temporal delay between each of the fields is the smallest for this Rayleigh number and amounts to aboutt rec = 266. The correlation coefficients have been normalized for this purpose and are shown in figure 9. Based on the normalized correlation coefficients, the typical time scales for the rearrangement of the superstructures are estimated to be τT ≈ 920, τũ z ≈ 320 and τ Nu loc ≈ 200 for the fields of the temperature, the vertical velocity component and the local Nusselt number, respectively. This confirms the aforementioned finding that Figure 9: Time correlation coefficient R αα for the experimental data obtained from the measurement at Ra = 2 × 10 5 . The typical time scales were estimated at a normalized correlation coefficient of R αα = 1/e, which is indicated with the horizontal dashed line.
in the present case the temperature patterns persist much longer in comparison to the characteristic patterns of the velocity field and of the corresponding local Nusselt number field. In order to further investigate the reorganisation of the turbulent superstructures and evaluate their role in heat transfer, a spatial and temporal median filter with a kernel size of nx × nỹ × nt = 9 × 9 × 3 was applied to the total data. In the smoothed fields of the temperature, the vertical velocity component and the local Nusselt number isosurfaces can be extracted that display the reorganisation in time effectively.
In figure 10, space-time isosurfaces of the time-averaged temperature field are shown forX Θ = ±0.13 at the top left and of the time-averaged vertical velocity component field forX uz = ±0.25 at the top right. Both fields are replotted in the lower row of the figure as a view from the top. In addition, isosurfaces ofX Nu loc = 0 are superposed in green colour within the panels showing the velocity isosurfaces. It is seen that the temperature field contours change gradually due to a certain reorganization of the superstructures. The isosurfaces partly bifurcate or end in the course of the time evolution. Furthermore, the isosurfaces of the vertical velocity component on the right side of the figure vary much stronger in space and time as it can be expected on the basis of the correlation coefficients in figure 9. In the view from the top, it can now be seen clearly that cells of strongly downwelling fluid in the middle (blue) and upwelling fluid in a kind of rim (red) do form. Even though these patterns remain relatively stable, especially for the case of the upwelling fluid a variation of the orientation of the ridges becomes apparent. These patterns do not necessarily coincide with the region of a larger temperature magnitude, which explains the decreased average of the local Nusselt number and the larger portion of negative local Nusselt numbers in the experiments, as obvious from table 3 and figure 4. It is interesting to note, that regions with negative local Nusselt numbers also remain for a certain timet ∼ O(10 3 ).
A more detailed look into one event with heat transfer towards the bottom plate can be seen in figure 11. In the top left panel a region with the temperature below the mean temperature is displayed. It is visible fort ∈ [1000, 3500]. If at the same time the vertical velocity points downwards, the local Nusselt number is positive. However, at the right hand side of the figure it can be seen in the isosurfaces of the vertical velocity component, Figure 10: Isosurfaces ofΘ (left) as well as ofũ z and Nu loc (right) for the experimental data obtained from the measurement at Ra = 2×10 5 . The lower row shows the isosurfaces from the top view.
that at a certain time instant the velocity points upwards and thus colder fluid will be transported towards the cooling plate. In this case, the local Nusselt number becomes negative, indicated by the green isosurface of the local Nusselt number fort ∈ [2500, 3500]. Later the velocity still points upwards but the structures in the temperature field have reorganized and the local temperature is not anymore colder than the mean temperature, such that the Nusselt number becomes positive again. From the current data it seems that this is a process that can be seen more often in the vicinity of the upwelling ridges. To summarize, a gradual evolution of the turbulent superstructures is clearly observable as our space-time analysis suggests.

Conclusion and outlook
In this work, the structure and evolution of turbulent superstructures in a horizontally extended Rayleigh-Bénard convection layer have been investigated in a controlled laboratory experiment. Particle image velocimetry and thermometry were combined for the measurements in a water-filled Rayleigh-Bénard cell with dimensions of l × w × h = 700 mm×700 mm×28 mm. The aspect ratio is Γ = l/h = 25, the Prandtl number Pr = 7, and the range of Rayleigh numbers 2 × 10 5 Ra 7 × 10 5 ; thus the experimental results are directly comparable with existing results of direct numerical simulations. Here, we accessed the field structure of the three velocity components and the temperature in their joint dynamical evolution in a wider observation area of about 451 mm×468 mm centered in the bulk of the cell to study the spatio-temporal organization of the turbulent heat Figure 11: Isosurfaces ofΘ (left) as well as of u z and Nu loc (right) in sections covering ∆x × ∆ỹ = 4 × 4 (top) and ∆x × ∆ỹ = 2 × 2 (bottom) for the experimental data obtained from the measurement at Ra = 2 × 10 5 .
transfer. Due to its transparent sidewalls and the transparent cooling plate made of glass the experimental setup provides full optical access to the turbulent flow. Simultaneous measurements of the temperature and velocity field in a large section of the horizontal midplane of the cell are thus feasible using thermochromic liquid crystals as tracer particles. Based on their temperature-dependent color shade the temperature field has been determined, while their temporal displacement has been evaluated to calculate both the horizontal and the vertical components of the velocity field in the observation plane.
The time-averaged fields of the temperature and the vertical velocity component, obtained from the present experiments and the DNS, clearly display turbulent superstructures of convection and agree qualitatively well. However, quantitative differences between the experimental and numerical results, in particular with respect to the characteristic wavelength of the large-scale patterns and the magnitude of the global Nusselt number, have been found. While the characteristic length scales of the turbulent superstructures are larger in the experiment, the values of the Nusselt number Nu are smaller compared to the results of the DNS. Both can be traced back to the deviations from the ideal isothermal boundary conditions at the top plate of the convection cell, which we discussed in detail in the text by an estimation of the corresponding Biot numbers. Despite the slight differences of the magnitude of the Nusselt number, the power law exponent of the transport law Nu(Ra) is in the same range as those of the DNS and other studies of thermal convection in water (Ahlers et al. 2009;Chillà & Schumacher 2012).
The joint measurement of the temperature and of the vertical velocity component has been applied to study the local convective heat flux in the experiment, thereby clearly demonstrating that the skeleton of turbulent superstructures contributes significantly to the turbulent heat transfer from the bottom to the top of the cell. The contribution of the turbulent superstructures to the heat flux is therefore quantified by means of the local convective Nusselt number, which reaches its maxima within the localized up-and downdrafts that compose to the superstructure patterns in the flow.
Since the measurements of the temperature and velocity fields cover up to approximately 10 4 convective free-fall times, our investigations also allow to study the reorganization of the superstructures for the different Rayleigh numbers. Their re-arrangement is analyzed here by means of the correlation coefficient of the time-averaged temperature fields. It is seen that the superstructures in the experiment have a preferential arrangement, most probably due to the slightly inhomogeneous temperature at the horizontal top plate. Nevertheless, the superstructures exhibit a clear gradual reorganization, which is in good agreement with the findings obtained in the previous numerical studies by Pandey et al. (2018) and Fonda et al. (2019).
Despite the non-ideal boundary conditions at the top plate, the presented results prove the concept of an experimental study of the long-term evolution of turbulent superstructures in Rayleigh-Bénard convection. A further improvement of the thermal boundary conditions at the top plate of the experimental facility is one point of our future work. An alternative approach to the determination of the structures, which are connected to the local convective heat flux, would be by three-dimensional, threecomponent volumetric (3D3C) measurements in a small volume. Again they would have to be combined with TLC measurements and require a larger number of simultaneously operating cameras. These studies are currently underway and will be reported elsewhere.