Supergranule aggregation for constant heat flux-driven turbulent convection

Turbulent convection processes in nature are often found to be organized in a hierarchy of plume structures and flow patterns. The gradual aggregation of convection cells or granules to a supergranule which eventually fills the whole horizontal layer is reported and analysed in spectral element direct numerical simulations of three-dimensional turbulent Rayleigh-B\'{e}nard convection at an aspect ratio of 60. The formation proceeds over a time span of more than $10^4$ convective time units for the largest accessible Rayleigh number and occurs only when the turbulence is driven by a constant heat flux which is imposed at the bottom and top planes enclosing the convection layer. The resulting gradual inverse cascade process is observed for both temperature variance and turbulent kinetic energy. An additional analysis of the leading Lyapunov vector field for the full turbulent flow trajectory in its high-dimensional phase space demonstrates that turbulent flow modes at a certain scale continue to give rise locally to modes with longer wavelength in the turbulent case. As a consequence successively larger convection patterns grow until the horizontal extension of the layer is reached. This instability mechanism, which is known to exist near the onset of constant heat flux-driven convection, is shown here to persist into the fully developed turbulent flow regime thus connecting weakly nonlinear pattern formation with the one in fully developed turbulence. We discuss possible implications of our study for observed, but not yet consistently numerically reproducible, solar supergranulation which could lead to improved simulation models of surface convection in the Sun.

Turbulent convection processes in nature are often found to be organized in a hierarchy of plume structures and flow patterns. The gradual aggregation of convection cells or granules to a supergranule which eventually fills the whole horizontal layer is reported and analysed in spectral element direct numerical simulations of three-dimensional turbulent Rayleigh-Bénard convection at an aspect ratio of 60. The formation proceeds over a time span of more than 10 4 convective time units for the largest accessible Rayleigh number and occurs only when the turbulence is driven by a constant heat flux which is imposed at the bottom and top planes enclosing the convection layer. The resulting gradual inverse cascade process is observed for both temperature variance and turbulent kinetic energy. An additional analysis of the leading Lyapunov vector field for the full turbulent flow trajectory in its high-dimensional phase space demonstrates that turbulent flow modes at a certain scale continue to give rise locally to modes with longer wavelength in the turbulent case. As a consequence successively larger convection patterns grow until the horizontal extension of the layer is reached. This instability mechanism, which is known to exist near the onset of constant heat flux-driven convection, is shown here to persist into the fully developed turbulent flow regime thus connecting weakly nonlinear pattern formation with the one in fully developed turbulence. We discuss possible implications of our study for observed, but not yet consistently numerically reproducible, solar supergranulation which could lead to improved simulation models of surface convection in the Sun.

I. INTRODUCTION
Turbulent convection, the essential mechanism by which heat is transported in natural flows, manifests often in a hierarchy of structures and flow patterns. Clusters of clouds over the warm oceans in the tropics on Earth [1] or giant storm systems in the atmospheres of the big gas planets Jupiter [2] and Saturn [3] illustrate this phenomenon. One of the most prominent astrophysical examples is the convection zone in the outer 30% of the Sun [4]. Convection cells are termed granules if they have an extension of G ∼ 10 3 km and a lifetime of about 10 minutes. These granules form the basic pattern at the solar surface where a heat flux drives convection [5,6]. Spectral observations reveal supergranules with extensions of SG ∼ 30 G and a lifetime of a day as the next larger building block in this hierarchy, detected either by line shifts in optical observations [7], helioseismology [8] or granule tracking [9]. Giant cells that extend across the whole convection zone could be a third stage in this hierarchy [10][11][12], but this is still an open question. Different physical effects have been proposed for the formation of supergranules, such as helium recombination in the upper convection zone [13], self-organisation of granules [9], or dynamical constraints by deeper convection at scales SG that are affected by the slow rotation of the Sun [14]. Numerical simulations that model convection and try to predict the spectral measurements have still been unable to develop supergranules self-consistently [15,16].
Here, we demonstrate the aggregation of granules to a large-scale supergranule in the simplest setting of convection without additional physical processes such as radiation, rotation or magnetic fields involved in heat and momentum transfer. This turbulent Rayleigh-Bénard convection (RBC) case in the Boussinesq limit is often considered as the paradigm for convective turbulence with its many facets [17,18]. Our three-dimensional direct numerical simulations (DNS) differ in three important ways from the majority of numerical studies in RBC: (1) they are subject to constant heat flux boundary conditions at the top and bottom; (2) they require simulations that are run on the order of 10 4 convective free-fall time units and even more; and (3) they are conducted in sufficiently extended layers. Layers with a fixed aspect ratio Γ = L/H = 60 with the horizontal length L and the layer height H are considered in this paper. We observe the gradual supergranule formation for all accessible Rayleigh numbers up to Ra ≈ 7.7×10 7 , a dimensionless measure for the vigor of convective turbulence. Their formation proceeds despite the fact that the flow becomes fully time-dependent and turbulent. This is in a regime for which one would not expect a pattern coherence across the whole domain, given that Ra is far beyond the critical Rayleigh number Ra c for the onset of the primary linear flow instability [19,20] or subsequent secondary instabilities of the onset pattern [21,22].
We confirm the continued gradual aggregation trend into the fully turbulent regime by the Lyapunov vector field of the largest Lyapunov exponent [23] of the tur-bulent states (see refs. [24,25,27] for similar analyses to characterize pattern defects in weakly nonlinear convection). The Lyapunov vector analysis probes here the growth of linear instabilities and of the corresponding scales of our high-dimensional nonlinear dynamical system. In particular, we show that the leading Lyapunov vector field becomes coarser as time proceeds which suggests that the turbulent flow remains unstable at a given scale with respect to longer-wavelength instabilities until the domain size has been reached. Indeed such an ongoing inverse cascade of energy and thermal variance can be clearly shown by the power spectra. The supergranule becomes visible in the velocity and temperature once a timewindowed averaging is applied that suppresses the turbulent fluctuations and the faster converging and diverging flows in the small-scale convective granules. We find supergranules independent of the boundary conditions of the velocity field.
Furthermore, we show that the supergranule is absent in DNS with constant temperature boundary conditions at the top and bottom planes of the RBC layer. For these cases, the formation of the recently comprehensively investigated turbulent superstructures [28][29][30][31][32][33][34][35][36][37][38] -well-ordered patterns of temperature and velocity with characteristic convection roll widths up to Λ/2 ∼ 3 − 4H [36-38] -takes place. Λ is the characteristic scale or wavelength of the pattern. Big velocity field condensates have been studied in two-dimensional [40] and quasi-two-dimensional fluid turbulence [41][42][43] to analyse the dependence of the inverse cascade on the energy injection scale. These settings are different to RBC where the driving of the fluid motion proceeds by thermal plumes that have a typical width of the order of the thermal boundary layer thickness δ T H L. A slowly progressing clustering of thermal plumes in RBC has been studied in von Hardenberg et al. [31] for L ≤ 12πH reaching roll widths of Λ/2 ∼ 3H for similar Rayleigh and Prandtl numbers. The generation of a large-scale anisotropy in turbulent convection for free-slip velocity conditions at the walls requires additional rotation about a horizontal axis in the threedimensional case as shown in ref. [44].
Hurle et al. [20] studied the linear stability of a two-dimensional thermal convection layer at rest for the constant flux case analytically. They detected a critical wavenumber k c = 0 and a critical Rayleigh number Ra c = 6! = 720 for no-slip velocity conditions and Ra c = 5! = 120 for free-slip conditions. This implies that the pair of counterrotating convection rolls at the onset of convection will always extend to the largest possible wavelength Λ = L in a finite cell. Instabilities of finiteamplitude convection rolls for Rayleigh numbers slightly above Ra c showed that each mode is unstable to one longer wavelength [22]. Interestingly, this gradual aggregation process has not been observed in previous turbulent RBC simulations with constant flux boundary conditions [45][46][47], most probably because they were conducted in smaller aspect ratio domains and for shorter total integration times.
Our study suggests that the mechanisms of supergranule formation in a simple convection flow are related to linear instabilities in the turbulent flow that give rise to longer-wavelength structures. Even though the RBC flow operates at Rayleigh numbers that are up to nearly 5 orders of magnitude above the critical Rayleigh number for the onset of the primary instability, a cell with the longest wavelength is still formed without any additional physical mechanism at work in the unstably stratified layer. Our investigation can thus shed a new light on the fundamentals of solar granulation processes.

II. NUMERICAL ANALYSIS
We consider here the simplest turbulent convection configuration, the three-dimensional Boussinesq case which couples the temperature field T (x, t) and the velocity vector field u(x, t) in an incompressible fluid [17,18] with u = (u x , u y , u z ) and x = (x, y, z). In this case the mass density is a linear function of the temperature deviation from the reference value. The Cartesian domain V = L × L × H applies periodic boundary conditions in both lateral directions, x and y for all fields. Regarding the vertical direction, the following boundary conditions are applied at the bottom and top plates at z = 0, H. For the velocity field, these are either no-slip (ns) or free-slip or stress-free (fs) conditions. They are given by (fs) ∂u x ∂z = ∂u y ∂z = 0 and u z = 0 .
with β > 0. We adopt as units of length and time the layer height H and the free-fall time t f = H/U f with the freefall velocity U f to rescale the equations in a dimensionless form. The latter is defined by The quantity α is the isobaric expansion coefficient and g is the acceleration due to gravity. In case of Neumann (N) boundary conditions U f = αgβH 2 . The dimensionless equations of motion follow as ∂T ∂t with the dimensionless pressure field p. The Prandtl and Rayleigh numbers are given by with the kinematic viscosity of the fluid ν and the temperature diffusivity κ. The presentation of the results is done in dimensionless units, e.g. 0 ≤ z ≤ 1.
The equations of motion are solved numerically with the spectral element method nek5000 [48,49]. The polynomial order N on each element and the total spectral element number N e are chosen properly such that the steep gradients near the top and bottom walls and the Kolmogorov scale η K can be resolved sufficiently (see [49] for more details). We varied the order N at fixed Rayleigh and Prandtl numbers to verify that the supergranule formation, the mean profiles of temperature, and the temperature variance spectra are unaffected (see appendix A for further details on the resolution tests). The boundary conditions (ns), (fs) for the velocity field and (D), (N) for the temperature field can be combined to four different sets of runs at several Rayleigh numbers. All four groups of boundary conditions were investigated. The main focus of our presentation will be on the series Nfs1 to Nfs4 which is listed in Table 1. This combination of boundary conditions comes closest to the solar convection case, that motivates our study. Note also that the Rayleigh numbers for cases Dfs and Nfs are related to each other by Since the runs with Dirichlet conditions served always as the starting point, the values of Ra N follow as given in Table 1.
The turbulent heat transfer across the convection layer is determined by the dimensionless Nusselt number which is given for constant flux boundary conditions (N) by and for the constant temperature conditions (D) by The symbol · A,t denotes a combined average over the horizontal cross section A = L 2 and time t. In comparison the global turbulent momentum transport is quantified by the Reynolds number, which is for numerical studies defined as III. SUPERGRANULE FORMATION Figure 1 illustrates the final stages of the simulations Nfs1 to Nfs4 at four different Rayleigh numbers between 10 4 Ra N 10 8 (see table 1). The top row shows the temperature contour snapshots taken in the final stages of the simulations close to the upper surface of the layer. A supergranule, in our case a pair of counterrotating rolls that fit into the domain, continues to exist into the fully developed turbulent regime and becomes clearly visible in all cases as a hotter background structure of the temperature field. Superposed is a fine-scale granule pattern that would also be observable for other RBC cases. It is related to the instability of fragments of the thermal boundary layer and the related thermal plume formation. The bottom row displays the corresponding surface velocity streamlines, again viewed from the top. They are obtained after a time-averaging over 500 convective time units centered around the temperature plots. Figure 2 underlines the different character of the largescale patterns of the Dirichlet and Neumann cases at the related Rayleigh numbers Ra ∼ 10 5 − 10 6 . We also de- where · t stands for a time average over an interval of 500 free-fall times t f . The horizontal cuts which are always taken at the height z 0 = 1 − δ T /2 decompose the data of the Dirichlet case into a time-averaged superstructure pattern with a characteristic wavelength that agrees with those in refs. [36,37] and a fine skeleton of plumeridges. This is in contrast to the Neumann case, where the supergranule is now clearly revealed. Interestingly, the instantaneous temperature fluctuation patterns for the Neumann case are more similar to those of the corresponding Dirichlet case. In appendix B we demonstrate furthermore that our observation is not altered by a substitution of free-slip with no-slip conditions which is consistent with the behaviour at the onset of convection [20].  Table 1. All fields are displayed for the whole cross-section of size L × L = 60 × 60. Figure 3 quantifies the gradual formation of the supergranule by quantities in the physical and Fourier spectral space for the case Nfs3 at Ra N = 3.93 × 10 6 . We note here that all simulations started with a perturbation of the linear diffusive equilibrium profile with random noise where the fluid is at rest. The very initial evolution where the flow relaxes into a statistically stationary regime for the fluctuations of the turbulent fields is discarded. In the course of the nonlinear evolution, a growth of the scale of the convection patterns, here the temperature field close to top plate, is observed in panels (a-d) which shows similarities to a phase separation process that has been analysed in binary-fluid mixtures [50,51]. The slow aggregation proceeds over ∼ 10 4 t f ; we find that the time this large-scale structure formation takes grows with respect to Rayleigh number (see also Fig. 1). This time span is longer than a vertical diffusion time scale τ v = √ Pr Ra t f , but significantly shorter than a horizontal diffusion scale Panels (e-h) plot the squared magnitude of the twodimensional Fourier transforms with respect to the horizontal coordinates x, y of the temperature. We display |T (k x , k y , z 0 , t 0 )| 2 , in logarithmically increasing contour levels at four different times t 0 for z 0 = 1 − δ T /2. The data correspond to those in the top row of the figure. We observe the slow transformation from the ring-like maximum, which is also observed in the Dirichlet case [28,37,38], to a condensate in the four (next neighbor) discrete planar wavevectors k 1,2 = (±k min , 0) and k 3,4 = (0, ±k min ) around the horizontally homogeneous mode with k = 0 (see panel h). They correspond to the largest wavelength the convection cell can take in a domain, namely Λ = Γ and thus k min = 2π/Γ ≈ 0.1. We have verified that a similar, but faster, aggregation takes place in boxes at smaller aspect ratio and the behavior is the same when the analysis is repeated in the midplane.
The panels in the bottom row of Fig. 3 show the azimuthally averaged spectra for different quantities: (i) temperature E T T (k, z = z 0 ) as in top and middle row, (ii) vertical velocity component E uzuz (k, z = z 0 ), and (iii) convective heat flux E uzT (k, z = z 0 ). For all temperature fields, that enter the spectral analysis, the area mean T (z 0 ) A,t is subtracted. Spectra are given by for uv = {T T, u z u z , u z T }. All quantities display clearly an accumulation of spectral density at the lowest wavenumber suggesting an inverse cascade process that leads to the formation of the supergranule.
Despite this strong aggregation in physical and Fourier space, the global heat transfer and thus the Nusselt number Nu N remain on average constant over the whole time period of the simulation in all runs as shown in the inset of Fig. 4 for one example. Figure 4 demonstrates, however, that the slow formation of the supergranule shifts the relative fraction of convective heat flux in the course of the evolution. We apply a filter in Fourier space and and the rest 1 − J kmin (z = 0.5, t) with k min = 2π/Γ as stated already above. It can be seen in all cases how the contribution of the supergranule structure gains importance for later times and reaches a statistically steady transport regime which is indicated in all runs. The share of the structure to the global transport drops from nearly 40% for the lowest to about 25% for the highest Ra which underlines its relevance for the turbulent heat transfer across the convection layer.

IV. LEADING LYAPUNOV VECTOR ANALYSIS
A better understanding of the physical mechanisms of the aggregation can be obtained by applying a technique that is well-established in dynamical systems -the Lyapunov analysis [23]. In this framework, the evolution of the turbulent convection flow corresponds to a trajectory in a very high-dimensional phase or state space (strictly speaking this state space is infinite-dimensional). The state of the fluid flow at time t is given by a column vector y(t) = (u(x k , t), T (x k , t)) which has k = 1, . . . , N e N 3 entries and thus N dof = 4N e N 3 1 is the number of degrees of freedom in the present numerical model (see also The sensitivity of the trajectory with respect to infinitesimal perturbations or in other words the tendency to develop new instabilities out of the present (fully turbulent) state can be probed by the strength of exponential separation of two initially very close trajectories, y(t) and y(t) + δy(t), of the present turbulent flow. The corresponding linearized equations to (15) are given by δẏ(t) = J (y(t))δy(t) with the Jacobian J = ∂F /∂y. Here, δy(t) = (δu(x k , t), δT (x k , t)) is the infinitesimal perturbation field to the original trajectory of the flow. In detail, this gives the following set of linearized Boussinesq equations for our study ∂δT ∂t where the pressure perturbation field δp is determined by δu with a Poisson equation similar to the original incompressible case (5)- (7). The determination of the spectrum of the first n of the total of N dof Lyapunov exponents, λ 1 ≥ λ 2 ≥ . . . λ n and their corresponding Lyapunov vector fields requires the simultaneous numerical solution of n versions of (16)-(18) with different initial perturbations δy 1 (0), δy 2 (0), . . . , δy n (0) in combination with the original equations (5)- (7). The computational complexity of this task in an extended turbulent convection flow limits us here to the leading Lyapunov exponent λ 1 (t) and the corresponding leading Lyapunov vector field which encodes the locations in the fluid volume and associated scales of the flow patterns that become unstable first. As pointed out by Levanger et al. [27], the local magnitude of the components of the Lyapunov vector indicate the sensitivity of these local regions with respect to perturbations. It thus contains the essential information that we need to explain the aggregation. The leading Lyapunov exponent λ 1 (t) is given by with a norm that is determined by The resulting spatial ridge patterns of the components of the leading Lyapunov vector reflect a critical wavelength (or scale) across which the instability is triggered. We have tested our algorithm for a RBC flow in the weakly nonlinear regime at Ra = 5 × 10 3 Ra c where this technique has been established [24][25][26][27], see appendix C.  ). The wavenumber range on both axes magnifies the smallest value to demonstrate the variance accumulation in the largest possible modes. The color bar is again the same for all panels in this row. (i-l) Additionally azimuthally averaged Fourier spectra with respect to k = (k 2 x + k 2 y ) 1/2 . The Fourier transform of three different quantities is displayed as indicated in the legend in panel (j): temperature variance spectrum, co-spectrum of the turbulent convective heat flux, and the kinetic energy spectrum with respect to the vertical velocity component. Figure 5 compares the mean vertical profiles of the temperature field, T (z, t) A , and the absolute value of temperature component of the (renormalized) leading Lyapunov vector field δỹ 1 which is denoted as |δT (z, t)| A for cases with Dirichlet and Neumann boundary conditions. For simplicity, we denoted the fourth component of δỹ 1 again by δT . We take the absolute value as we are interested in the magnitude only. Here · A denotes an average over the cross section plane A = L 2 . While the mean temperature profile of Dfs2 has a zero slope in the bulk, the one for Nfs2 is slightly stably stratified. Such a subadiabaticity of the temperature is considered as a possible origin of the emergence of supergranules in the Sun [11,16]. The mean profiles of the Lyapunov vector component are found to differ qualitatively as reported in the same figure. In case of Dirichlet boundary conditions, the most sensitive region is located at the top of the thermal boundary layer where plume mixing starts. This is in contrast to Neumann boundary condition case, where the bottom and top planes are found to be by far the most susceptible with respect to small perturbations.
We select the plane that corresponds to one of the two local maxima in Fig. 5(a) at z 0 ≈ 0.94, i.e., where the instability has the largest magnitude, and show in Fig. 6 (a-h) that local instabilities and thus the maxima of |δT (x, y, z 0 , t)| are found close to (but not at) the downflow regions in Dfs2. They also remain connected to local creation or annihilation of defects of the flow patterns for the turbulent and fully time-dependent Dirichlet boundary case. Despite operating in the turbulent regime of the flow, these instabilities are thus very similar to what is found for the weakly nonlinear regime. We display therefore a magnification of a short dynamical sequence in panels (b-d, f-h). In the corresponding turbulent Neumann boundary case Nfs2 in Fig. 6 (i-p) at Ra N = 2.04×10 5 , the structure of the temperature perturbation |δT (x, y, z 0 , t)| is different. The plane that was now selected is taken close to the top wall at z 0 = 1 − δ T /2 in correspondence with Fig. 5(b). One observes a connected pattern of highamplitude ridges of δT with a coarser spacing indicating a larger scale of instability. It is also observed now that the locally most unstable regions coincide with the downflow regions thus stabilizing the bulk regions (see also [11,52] for similar mechanisms in solar case). Figure 7 provides additional results on the instabilities as well as on the corresponding scales. First, we show time series of λ 1 (t) and Nu(t) for Dfs2 (a) and Nfs2 (e).
Both values vary with respect to time about the mean values which is typical for the turbulent flow case. Furthermore, Fourier spectra of the temperature, E T T (k), and the temperature component of the leading Lyapunov vector, E δT δT (k), are shown. Data are taken in planes close to the top where maximum magnitudes of δT are found. The Dirichlet case displays a local maximum at k = k m that coincides for both spectra (see panels (b,c)). This local spectral peak corresponds to a characteristic wavelength which is given byλ which remains unchanged in time at a value ofλ δT δT ≈ 4 and thus corresponds to the characteristic extension of the turbulent superstructures which have been discussed in Pandey et al. [37] and Fonda et al. [34]. The results differ for the Neumann case as displayed in panels (f-h) of the figure. The spectrum E δT δT (k) has a large-wavenumber bump at k m ≈ 2 which is indicated by the dashed line in Fig. 7(h). These instabilities corresponds to the fine granule patterns which are for example seen in Fig. 2(b). While this local maximum remains unchanged in time, a second one moves gradually towards larger wavelengths, see again Fig. 7(h). This suggests that the turbulent flow develops instabilities at an increasingly larger wavelength. The process is ceased when the system size is reached and the nonlinear process of supergranule formation is completed. We stress once more that this behavior is fundamentally different to the case with Dirichlet boundary conditions for the temperature field.

V. DISCUSSION AND PERSPECTIVE
Our main motivation was to demonstrate the gradual formation of a salient large-scale convection pattern on a time scale larger than the vertical diffusion time Pr Ra t f and that eventually fills the whole convection domain in a Rayleigh-Bénard convection setup. Following solar convection, this structure is termed supergranule. We showed that this formation proceeds only in case of constant heat flux boundary conditions (also known as Neumann conditions) at the top and bottom planes of the layer, independently of no-slip or free-slip boundary conditions for the velocity field. Surprisingly the supergranule is still observed when the flow is in the state of fully developed turbulence as being the case, at least for runs Nfs3 and Nfs4. We mention here simulations of compressible photospheric convection (with similar boundary conditions) by Rincon et al. [53] at Γ = 42. The authors report the formation of a dominant convective mode and conclude that the simulations could not be run long enough to study a further aggregation. Our results confirm that long simulation times are necessary and demonstrate that this dominant convection mode is eventually a supergranule which can be seen even in the simpler (incompressible) RBC setup.
As discussed, the critical mode at onset of Rayleigh-Bénard convection in an infinitely extended layer with FIG. 6. Lyapunov vector fields for Dirichlet and Neumann boundary conditions. The time evolution of patterns in the temperature field T (x, y, z0, t) (a-d,i-l) and the leading Lyapunov vector temperature perturbation field |δT (x, y, z0, t)| (e-h,m-p) is visualised. We show data for run Dfs2 at z0 ≈ 1 − 2/3δT ≈ 0.94 (a-h) and run Nfs2 at z0 = 1 − δT /2 (i-p). The value of z0 ≈ 0.94 in Dfs2 corresponds with one of the two maxima in Fig. 5(a). Panels (a,i) show T and panels (e,m) δT over the whole cross-section of size L × L = 60 × 60. The remaining subplots enlarge a highlighted region of interest of size 10 × 10. The times of the corresponding snapshots are indicated above. See also appendix C for the same analysis of Dfs2 in the midplane.
Neumann conditions at the top and bottom is k c = 0. This implies that a single pair of counter-rotating convection rolls fills a domain with a finite periodicity length L at onset, a behavior which is found in this work to persist to Ra N Ra N,c . We confirm the behavior by detecting the accumulation of kinetic energy and thermal variance in the four next neighboring Fourier modes to the (critical) zero mode with wavelength L. The determination of the leading Lyapunov vector field and the subsequent spectral analysis of its temperature component, demonstrates clearly that the flow structures at a given scale give rise to an instability at a next bigger wavelength and thus to a spatially larger new flow structure. This inverse cascade continues until the horizontal periodicity length is reached for the present setup. In the solar case a further physical process will limit the supergranule size.
We thus demonstrate that the structure formation mechanism, which was described in Chapman and Proctor [22] above the onset of convection in the weakly nonlinear regime, persists far into the turbulent range. A possible next step would be to derive effective amplitude equations, now for the perturbations about a fully turbulent state. This will include turbulent closures and certainly requires simplifying assumptions, but could be done along lines of a very recent work by Ibbeken et al. [54]. Our Lyapunov vector analysis answers furthermore a question left open in [37]: the generation of turbulent superstructures in the Dirichlet case is a local pattern instability with a scale of the size of a pair of counter-rotating mean circulation rolls, hereλ δT δT ≈ 4. In contrast, the Neumann case proceeds slowly to a global wavelength instability by a cascading process withλ δT δT → L.
Finally, we return to the initial example of solar convection where the fixed heat flux at the top is connected with the well-known solar luminosity L . This flux is the main driver of convection and thus the formation of granules and supergranules in the upper convection zone. Our study showed that already these boundary conditions alone generate a large-scale convection roll pair, i.e., without additional magnetic fields, changes in chemical composition and the strong compressibility effects. As the typical scale ratio SG / G ≈ 30 of the solar convection case is equal to the diameter ratio of our supergranule to granule roll for the prescribed layer extension, namely (Λ/2)/H ≈ 30, we want to compare now characteristic velocities and evolution times of Nfs4 with the solar data given in the introduction. We thus decompose u(x, t) = U (x) + u (x, t) with U (x) = u(x, t) t . The ratio of the corresponding root mean square velocities u rms /U rms ≈ 5.8 comes close to the velocity ratio of v G /v SG ≈ 6 [10]. When the lifetime of a granule is estimated by the mean turnover time of a Lagrangian tracer across the layer with a value of t to ≈ t u ∼ 20 [55], one arrives at t U /t u ∼ 10 4 /20 ∼ 500 which is at least of the same order of magnitude as τ SG /τ G ≈ 144 [10].
Clearly, this approximate agreement should not be overestimated as the solar convection zone contains a much more complicated physics at a much larger Rayleigh num-ber and an extremely small Prandtl number, Pr 10 −6 [4]. Nevertheless, our simple convection model might still turn out to be fruitful to better interpret the solar observations as we were able to reveal a basic instability mechanism in this class of turbulent flows that leads to a large-scale flow structure. It is thus also a good starting point for a step-by-step increase of complexity towards the solar case that can test how the supergranule formation is affected by an inclusion of further physical processes. RaDPr uzT A,t (a), the kinetic energy dissipation rate A,t (b), and the thermal dissipation rate T A,t (c).

APPENDIX A -SPECTRAL RESOLUTION TESTS
In spectral element methods, the resolution depends on two factors: the total number of spectral elements (N e ) and the polynomial order of the spectral expansion on each element and in each spatial direction (N ). The method belongs to the bigger class of exponentially fast converging hp-FEM (FEM=finite element method) where equations are solved on elements that fill the volume by a piecewise polynomial approximation of the solutions. One can vary the size of the element h (here N e ) or the polynomial degree p (here N ). In our flow, the polynomial order and spectral element number has to be chosen properly such that the steep gradients near the top and bottom walls and the Kolmogorov scale η K can be resolved sufficiently. Sufficient resolution is established once, dz(z)/ η K (z) A,t π/2 for Pr 1. Here dz(z) is the vertical spectral element extension. This criterion was suggested and tested in ref. [49]. To show the convergence of our results, we perform first a resolution test with respect to two different polynomial orders N , as shown in Fig. 8. The production run setup is at N e = 160, 000 and N = 11 (case Dns2). It can be seen that the same results can be achieved already with lower polynomial order of N = 7 and that the curves for both N collapse. We conclude that the spectral resolution with N = 11 is sufficient.
Furthermore, we demonstrate that the gradual supergranule formation is not a resolution effect. This is done in a smaller cell of aspect ratio Γ = 15 to accelerate the formation process. We use the parameters of simulation run Nfs2 with Ra N = 2.04×10 5 and the corresponding run Dfs2 with Ra D = 3.85 × 10 4 . We apply the same spectral element resolution as in the main text, which translates to N e = 10, 000. In the corresponding comparison runs, we double the number of elements in each space direction which leads to N e = 80, 000. The supergranule evolves in the long-term dynamics in case of Neumann boundary conditions, while there is no such effect for Dirichlet boundary conditions. The results are summarized in Fig.  9. In panel (a), it is seen that the convective heat flux profiles collapse on each other for both pairs of runs. The Fourier spectra in panels (b,c), which are taken for the Dirichlet run from 750 to 950 t f and for the Neumann run from 2,500 to 2,700 t f , display the aggregation in the latter case which agrees very well with the results for Γ = 60. In the Neumann boundary case, the supergranule is already fully developed in both runs. Note that in most panels the curves collapse onto each other.

APPENDIX B -VELOCITY BOUNDARY CONDITIONS
We compare in Fig. 10 four different combinations of temperature and velocity field boundary conditions for snapshots of the temperature field at z = 1 − δ T /2 and at a late stage of the dynamical evolution. The panels in the leftmost column coincide with those in Fig. 1. It can be seen that the instantaneous temperature patterns have very different characteristics for the Neumann and Dirichlet cases. It is also seen that the change of temperature boundary conditions is the essential one (and not the change of the velocity boundary conditions) that leads to the supergranule. All fields are visualized for the whole cross-section of size L × L = 60 × 60. In order to have the same distance from the onset of convection, we take Ra D = 3.85 × 10 3 , 3.85 × 10 4 , and 3.85 × 10 5 for runs Dfs1, Dfs2, and Dfs3, respectively. The corresponding two series with Neumann boundary conditions follow by eq. (9). Thus Ra N = 2.23 × 10 4 , 4.34 × 10 5 , and 8.31 × 10 6 for runs Nns1, Nns2, and Nns3, respectively. The corresponding Rayleigh numbers for Nfs1, Nfs2, and Nfs3 are listed in Table 1.

APPENDIX C -LYAPUNOV VECTOR DETERMINATION
We provide in Fig. 11 two series of snapshots that show the evolution of the temperature field T (x, y, z, t) together with the temperature component of the corresponding leading Lyapunov vector, δT (x, y, z, t), for a combination of Dirichlet and no-slip boundary conditions. Panels (a)-(h) are taken in the weakly nonlinear regime at Ra D = 5.0 × 10 3 . This run is a test of our routine as it can be compared with results in the listed references, e.g. Egolf et al. [24] or Scheel et al. [25]. As in those refer- FIG. 9. Spectral resolution study with respect to the number of spectral elements Ne at fixed polynomial order N = 11. A comparison of Dirichlet and Neumann boundary conditions is shown as indicated in the legend at the top. The aspect ratio is Γ = 15, the Prandtl number is fixed at Pr = 1, and free-slip boundary conditions are applied at the plates for the velocity field in all simulations. The lateral boundary conditions are periodic. The Rayleigh number for the case of constant flux RaN is agrees with case Nfs2, whereas RaD is equal to that of case Dfs2 (see Table 1). Profiles in panel (a) show convective heat transfer uzT A,t. The data for the Neumann case have been multiplied by a factor of 10 for better visibility. Plots (b,c) display the time-and azimuthally averaged Fourier spectra of temperature and vertical velocity component in the midplane. The time average is performed over a time span of 200 convective free-fall times in all cases. ences, the Lyapunov vector field highlights the regions of instability, where the defect formation is observable as a bright spot. The time-averaged leading Lyapunov exponent λ 1 = 252±1. Panels (i)-(p) are for the turbulent flow case Dfs2 which is also discussed in the main text. Again, local maxima of the Lyapunov vector field correspond to a defect generation. The time-averaged leading Lyapunov exponent λ 1 = 2703 ± 7. All data which are shown in this figure are taken in the midplane z = 0.5. The appearance of the localised defect is clearly detectable in the leading Lyapunov vector field for both cases, see panels (c,g,k,o) fo the figure.