Lagrangian studies of coherent sets and heat transport in constant heat flux-driven turbulent Rayleigh-B\'enard convection

We explore the mechanisms of heat transfer in a turbulent constant heat flux-driven Rayleigh-B\'enard convection flow, which exhibits a hierarchy of flow structures from granules to supergranules. Our computational framework makes use of time-dependent flow networks. These are based on trajectories of Lagrangian tracer particles that are advected in the flow. We identify coherent sets in the Lagrangian frame of reference as those sets of trajectories that stay closely together for an extended time span under the action of the turbulent flow. Depending on the choice of the measure of coherence, sets with different characteristics are detected. First, the application of a recently proposed evolutionary spectral clustering scheme allows us to extract granular coherent features that are shown to contribute significantly less to the global heat transfer than their spatial complements. Moreover, splits and mergers of these (leaking) coherent sets leave spectral footprints. Secondly, trajectories which exhibit a small node degree in the corresponding network represent objectively highly coherent flow structures and can be related to supergranules as the other stage of the present flow hierarchy. We demonstrate that the supergranular flow structures play a key role in the vertical heat transport and that they exhibit a greater spatial extension than the granular structures obtained from spectral clustering.


I. INTRODUCTION
Turbulent thermal convection, the essential mechanism by which heat is transported in many natural flows [1,2], has numerous geophysical [3], astrophysical [4] and technological [5] applications. Rayleigh-Bénard convection (RBC) is a paradigm of such natural thermal convection systems and consists basically of a fluid layer that is placed between two solid horizontal plates which uniformly heat and cool it from below and above, respectively [2]. Crucially, the formation of large coherent flow structures can be observed -similar to the natural convection flows -even in this simple system once a horizontally extended domain is provided [6][7][8][9][10][11]. In the Eulerian frame of reference, these long-living large-scale flow structures are termed (in the standard case where the plates possess spatially constant temperatures) turbulent superstructures as their characteristic horizontal scale exceeds the height of the convection layer. They consist of convection rolls and cells that are concealed in instantaneous velocity fields by turbulent fluctuations. However, they show up prominently after time-averaging of the velocity or temperature fields [7] and form the backbone of turbulent heat transport [12]. The situation is -ex-cept for the formation of an entire hierarchy of different long-living large-scale flow structures -roughly similar if the plates possess spatially constant temperature gradients or heat fluxes (more details follow further below) [10,11]. The description of convection in the Lagrangian frame of reference is connected to the material transport and thus opens the perspective of a more precise quantification of the complex spatio-temporal pathways that heat takes through the fluid layer, including a classification into spatial regions that contribute more or less to its transfer.
Lagrangian transport and mixing processes have been intensively studied by means of mathematical and computational tools from dynamical systems theory [13][14][15] and the frameworks are built around different, but related notions and definitions of coherent flow structures. The concept of a finite-time coherent set [16][17][18], a regularly shaped region in the fluid volume that only weakly mixes with its surrounding, is central to the present study. The boundaries of such regions can be identified within a geometric approach, where Lagrangian coherent structures represent minimal curves or surfaces that enclose coherent sets [14]. Coherent sets were originally introduced within a probabilistic approach based on arXiv:2304.02984v3 [physics.flu-dyn] 14 Sep 2023 (a) t = 100 (b) t = 300 (c) t = 700 (d) t = 1, 500 (e) t = 3, 500 1.0 0.5

2.0
T (x, y, z = 1 T /2, t) FIG. 1: Rayleigh-Bénard convection, which is driven by an applied constant heat flux, exhibits a gradual supergranule aggregation. This time series visualizes instantaneous temperature fields T close to the top plane (with δ T = 1/(2Nu) representing the thermal boundary layer thickness) for the simulation run that is at the focus of this present paper. The large and gradually growing patterns are termed supergranules, whereas the smaller superposed structures are termed granules and exhibit a time-independent size. Note that the domain's horizontal extent is 30 times its vertical extent and find more details on the configuration in section IV.
transfer operators [16,17] whereas recent methods make use of spatio-temporal clustering algorithms applied to Lagrangian trajectory data [19][20][21][22][23][24][25][26]. These aim at identifying coherent sets as groups of trajectories that remain close to each other during the time interval under investigation. Extensions have been proposed to study the long-term evolution -including the emergence and decay of coherent sets [27,28] -, while recent reviews discuss complex networks [29,30] and machine learning [31] approaches for fluid flows. Finite-time coherent sets in turbulent Rayleigh-Bénard convection, their relation to the Eulerian turbulent superstructures, and their role in heat transport have been studied using different Lagrangian approaches such as trajectory-based clustering [24,32,33], transfer operator methods [34], and evolutionary clustering [28]. In these works, aspect ratios (i.e., the ratio of the domain's horizontal to vertical extent) of Γ ∈ [2, 16] have been considered together with applied constant temperatures at the bottom and top plates (so-called Dirichlet boundary conditions). In contrast, when the convection is driven by a constant heat flux (representing Neumann boundary conditions), a gradual aggregation of smaller, but still large-scale convection cells -termed granules -to even larger structures -termed supergranules -is observed, see Figure 1. Note that the characteristic horizontal extension of the granular flow structures is O (H) (with H being the layer height), whereas the supergranular structures grow until eventually being limited by the horizontal extent of the domain [10]. This growth can only be interrupted once the fluid layer is subject to additional physical mechanisms such as rotation around the vertical axis [11]. However, this new mechanism has so far only been investigated in the Eulerian frame of reference [10,11]. The aim of the present paper is to extend the Lagrangian evolutionary framework [28] to this setting in order to address the following questions: 1. How do Lagrangian coherent sets or features relate to granules and supergranules in the heat flux-driven convection case?
2. What is the role of granules and supergranules as the building blocks of the structural hierarchy on the global turbulent heat transport across the convection layer?
3. Can we observe the gradual supergranule aggregation in the Lagrangian frame of reference despite the superposition of different long-living large-scale flow structures?
Particularly question 3 will require the evolutionary framework of the trajectory clustering analysis which we developed and outlined in [28] since the process will be affected by the slow transient aggregation.
The remainder of this paper is organized as follows. In Section II, we briefly review the concept of Lagrangian coherent features in flows and the corresponding computational methods, focusing on the constructions introduced in [28,35]. In Section III, we introduce the mathematical model for constant heat flux-driven Rayleigh-Bénard convection and briefly describe the results of our previous Lagrangian studies. Section IV provides objective Lagrangian characteristics of the flow and introduces subsequently methodical details of our following investigations. The results of the spectral clustering and the conditional heat transfer analysis are presented and discussed in Sections V and VI, respectively. We conclude with a summary of the findings and an outlook in Section VII.

II. LAGRANGIAN COHERENT FEATURES
From the nonlinear dynamics point of view, there are a number of different concepts that describe the notion of Lagrangian coherent behavior in flows. Neglecting molecular diffusion, the motion of passive fluid particles is described by the ordinary differential equatioṅ X = u (X, t) . (1) Here, u (X, t) is a sufficiently smooth time-dependent velocity field with state X ∈ X, time t ∈ R, and a compact subset X ⊂ R dim (dim ∈ {2, 3}) which is the physical domain of the flow. Lagrangian particle trajectories are obtained as solutions X (t) of this differential equation (1). Coherence can then be described in terms of the manner in which groups of Lagrangian trajectories behave. There are essentially three families of Lagrangian concepts for the identification of coherent features in a flow: geometric, probabilistic, and clustering approaches. Whereas the latter are data-based, the geometric and probabilistic methods have strong mathematical foundations in dynamical systems and ergodic theory. For discussions of the most frequently used frameworks, we refer to [13,15]. We briefly review the different concepts below before we describe an evolutionary trajectory network (clustering) approach [28] that will be the basis for our studies.

A. Methods for identifying coherent behavior
Geometric approaches aim at identifying Lagrangian coherent structures, i.e. material surfaces that extremize a certain stretching or shearing quantity such as measured by the finite-time Lyapunov exponent. These material surfaces act as barriers to advective mass transport [14]. The classical concept has been recently extended to deal with weakly diffusive and stochastic transport [36] and with transport of other quantities [37][38][39]. Lagrangian coherent structure methods require full knowledge of a smooth flow field, and so the extraction of evolving material surfaces in complex 3D flows remains challenging.
In contrast, probabilistic methods aim to identify fullvolume finite-time coherent sets that minimally mix with the surrounding phase space. Coherent sets can be efficiently identified from k leading singular vectors of sparse stochastic matrices that numerically approximate Perron-Frobenius (or transfer) operators, the latter of which describe the time evolution of densities in the phase or state space of the dynamical system. k is determined based on a spectral gap criterion. The structures of interest are then extracted from the corresponding k singular vectors via a hard clustering scheme such as kmeans or a sparse eigenbasis approximation [40].
Related work considers Fokker-Planck equations [41], a dynamic Laplacian approach [42] with a finite-elements approximation [43], and a geometric heat flow construction [18]. Bifurcations of coherent sets, such as their splits and mergers, have been studied both in transfer operator [44,45] and dynamic Laplacian settings [27]. In [46], the different time scales of coherent sets have been studied.
While clustering as a class of unsupervised machine learning algorithms allows to uncover underlying patterns in the data, different methods have already been proven to effectively extract coherent sets from Lagrangian trajectories. Given an appropriate similarity measure between trajectories (e.g. based on the average distance [19,20]), clustering results in groups of trajectories that are dynamically close to each other. First works in this context apply c-means clustering [19], spectral clustering [20,21,23], graph coloring [22], density-based clustering [24][25][26]47] and complex networks [23,48]. Interestingly, spectral clustering methods can be mathematically related to the transfer operator approaches [21]. Such methods allow us to extract and evaluate coherent sets for a given time window. To study, in contrast, the long-term evolution of coherent features -including the emergence and decay of coherent sets -, it is necessary to relax their material properties [27,28]. To this end, in [28] the trajectory-based network approach [23] has been extended to an evolutionary spectral clustering framework [49] which blends a time-dependent analysis over shorter time spans with a dynamic short-term memory.

B. Lagrangian trajectory networks and evolutionary spectral clustering
In the following, we will briefly review the concept of trajectory networks and spectral clustering approaches, focusing -due to their relevance to the present workon the constructions introduced in [23,28].
Suppose we are given N tracer trajectories X i (t) with i = 1 . . . N at discrete time instances t ∈ T = {0 . . . Θ}. Based on these trajectories, we construct an undirected network with Lagrangian particle trajectories serving as network nodes, where a link between two nodes is established if the respective trajectories come "close" to each other. There are different definitions of the distance between two trajectories and thus of "closeness". In [23], an unweighted network was considered with an adjacency matrix A ∈ {0, 1} N ×N , where A ij = 1 if min t∈T ∥X i (t) − X j (t)∥ < ε for i ̸ = j and A ij = 0 otherwise. ε > 0 is thereby some given threshold. In [28,50], the number of ε-close encounters was further included as link weights by first constructing instantaneous adjacency matrices and A ij,t = 0 otherwise, and then forming the symmetric network weight matrix W = t∈T A t , i.e., W ∈ N N ×N 0 . In other words, an entry W ij is large if X i and X j are close in the sense that they have many ε-close encounters during the time span T.
As the choice of ε is often tricky, especially when the given trajectory data is sparse, we propose and will use an alternative construction based on mutual K nearest neighbours. Here, the instantaneous adjacency matrix has the entry A ij,t = 1 if the particles X i (t) and X j (t) with i ̸ = j are K mutual nearest neighbours in the sense that X i (t) belongs to the K nearest neighbours of X j (t) FIG. 2: Illustration of the instantaneous network construction using mutual nearest neighbours. Here for each node its K = 3 neighbours are linked (bold and dashed lines), but only those links are kept where the K-neighbourhood is mutual (bold lines). For instance, the 3 nearest neighbours of the magenta node are blue, green and orange, however the magenta node does not belong to the 3 nearest neighbours of the orange node. So the link between orange and magenta (dashed) is not considered.
and simultaneously vice versa. Figure 2 provides an illustration of this idea. In this way, the link density as well as the connectivity can be controlled very conveniently by the choice of K. By including only mutual neighbourhoods, the adjacency matrices and thus the resulting weight matrix are symmetric just as in the εneighbourhood construction. We compare the different graph constructions based on K and ε in A. Node degrees are simple and useful network measures that evaluate the connections within the constructed graph. In our context, the unweighted node degree d i is defined by the number of different trajectories that are linked to node X i . That is, d i = N j=1 W ij where W ∈ {0, 1} N ×N is the adjacency matrix of the network, encoding only the presence of links and thus being related to W . This network measure can be related to the finitetime Lyapunov exponent [48] as well as to the trajectory encounter volume [51].
To identify coherent sets from the weighted trajectory network, we make use of spectral clustering [52,53] with the weight matrix W serving as a similarity matrix. In order to partition the trajectories into k groups, we aim to subdivide the graph such that the within-cluster similarity is maximized whereas the between-cluster similarity is simultaneously minimized. The solution to this is a k-way normalized cut. These optimization problems as formulated in [52,53] are NP-hard, but spectral relaxations have been proposed. We follow [53] and minimize the cost function over all orthogonal matrices Ξ ∈ R N ×k . D ∈ R N ×N represents a diagonal matrix with D ii = N j=1 W ij . This minimization problem is solved by the matrix Ξ if it consists of the eigenvectors corresponding to the k largest eigenvalues of the matrix D −1/2 W D −1/2 as columns.
Complementing the nodes' coordinates captured in these eigenvectors, gaps in the corresponding eigenvalue spectrum provide information on the intrinsic connectivity of the graph [21] and thus the potential numbers of coherent sets k [28] in a data-driven way. In our context, we use a simple spectral gap criterion, guided by the physical knowledge of the length scales of structures of interest. Different heuristics how to choose an appropriate k without a priori knowledge can be found e.g. in [40,54]. We note that the choice of k also effectively fixes the size of the structures to be identified.
To extract k coherent sets from these coordinates in eigenspace, we post-process the eigenvectors using the sparse eigenbasis approximation (SEBA) algorithm [40]. In a nutshell, this algorithm transforms the set of eigenvectors Ξ to a new set of vectors S ∈ R N ×k which spans approximately the same subspace but is significantly sparser. As different coherent features or sets are now separated into different sparse vectors (as the columns of S), the entry S ij can be interpreted as the likelihood of trajectory X i to belong to cluster j. In this work, we are only interested in the maximum likelihood of cluster affiliation of each trajectory (to belong to any cluster) and so we define subsequently the cluster membership indicator S max, i = max j S ij .
In [28], such a network-based approach was extended to study the evolution of large-scale coherent sets in turbulent flows over long time spans. We briefly review the evolving cluster approach as the one of the two proposed frameworks which will be used in the present study. For more details, we refer to [28].
The evolving cluster approach is based on the preserving cluster membership framework from [49]. Let W t denote the time-dependent weight matrices that are now computed over a shorter observation window of length ∆t ow ≪ Θ and centered around some time t ∈ Tmore precisely, W t = s∈Tt,∆tow A s with T t,∆tow = {t − ∆t ow /2, . . . , t + ∆t ow /2} for times t ∈ T. Now, for each t, we minimize (cf. [49]) which is (again) solved by Ξ t ∈ R N ×k if it contains the eigenvectors corresponding to the k largest eigenvalues of the matrix Note here that Ξ t−∆ts represent the eigenvectors obtained from the clustering at the previous time instance t − ∆t s , where ∆t s denotes the time shift. The obliviousness parameter µ regulates the importance of the current connectivity of the graph compared to the previous clustering, and using µ = 1 would lead to the usual clustering over the current time interval [t − ∆t ow /2, t + ∆t ow /2], similar to eq. (2). However, note that W t is (once µ < 1) no longer sparse due to Ξ t−∆ts Ξ T t−∆ts being a dense matrix. To overcome this issue, we use instead the sparse approximation S t−∆ts of the eigenspace spanned by Ξ t−∆ts (see also [28] for more details).
In the following, we denote by S max [X i (t)] the cluster membership indicator for the trajectory fragment X i (t) over the time interval T t,∆tow (i.e. centred at t), similarly

III. TURBULENT RAYLEIGH-BÉNARD CONVECTION
In this section, we introduce the present configuration of Rayleigh-Bénard convection and briefly review the results of previous Lagrangian studies.

A. Model definition and computational details
We study Lagrangian coherent sets in a threedimensional Cartesian domain using the model setting introduced in [10]. Therefore, the equations of motion are made dimensionless based on characteristic quantities of the system. In addition to the layer height H as unit of length, we adopt the free-fall time scale τ f = H/U f as the unit of time together with the corresponding freefall velocity U f = √ αgT char H. Here, α represents the isobaric expansion coefficient and g the acceleration due to gravity. T char is a characteristic temperature that depends on the thermal boundary conditions and thus on either the prescribed vertical temperature gradient β > 0 or the applied temperature difference ∆T > 0 between the plates. Hence, the governing equations follow in their non-dimensional description as where p, u, and T are the pressure, velocity, and temperature field, respectively. Two dimensionless parameters -the Prandtl number Pr and the Rayleigh number Ra -control the entire dynamics and summarize all the dimensional coefficients and parameters. These numbers are defined as with the kinematic viscosity and temperature diffusivity of the fluid, ν and κ, respectively. These equations are complemented by the numerical domain and its corresponding boundary conditions. With the aspect ratio as the ratio of the domain's horizontal to vertical extent, Γ = L/H, the domain extends (non-dimensionally in units of the layer height H) across x, y ∈ [−Γ/2, Γ/2] and z ∈ [0, 1]. The impermeable bottom and top planes at z ∈ {0, 1} obey free-slip mechanical boundary conditions whereas the lateral boundaries are closed and thermally insulating (walls), i.e., with n and t being the corresponding wall-normal and -tangential coordinates, respectively. In contrast to the classical scenario (as in [33]) where one applies constant temperatures at the top and bottom planes such that T (z = 0) = 1 and T (z = 1) = 0 with T char = ∆T , we prescribe here a constant vertical temperature gradient The characteristic temperature is thus T char = βH. Note that in this present case, the temperature at the top and bottom planes is allowed to vary locally -all that is fixed is the applied constant heat flux through the spatially constant temperature gradient. We solve the equations of motion using the spectral element solver nek5000 [55,56] and resolve all dynamically relevant scales down to the Kolmogorov scale η K . Any simulation starts with the fluid at rest and the linear conduction temperature profile, the latter of which is randomly and infinitesimally perturbed. The temporal advancement of the simulations exploits a second order backwards difference formula. More detailed information on numerical details can be found in our previous works [7, 10-12, 33, 56].
Finally, the Nusselt number Nu is the global measure to quantify the relevance of convective heat transfer across the fluid layer. For constant heat flux boundary conditions, it is given by [57] with A = Γ × Γ being the non-dimensional horizontal cross section. While the spectral element code works in the Eulerian frame of reference, the advection of N Lagrangian and thus massless tracer particles is an analysis in the Lagrangian frame of reference. These tracers follow the surrounding flow perfectly and so their trajectories X i (t) are given byẊ see again eq. (1). The velocity field components in this equation are therefore obtained from a spectral interpolation of the Eulerian fields to the tracer position X i (t). Hence, the spatial aspect of the advection of the Lagrangian particles is as accurate as the simultaneous temporal advancement of the Boussinesq equations (5) - (7) in the Eulerian frame of reference. Concerning the temporal aspect of this particle advection, a third-order Adams-Bashforth time stepping is exploited.
In addition to the global measure of turbulent heat transport (see above), we compute for each individual tracer particle a local Nusselt number [33] Nu in the Lagrangian frame of reference along its trajectory X i (t). To quantify eventually the overall heat transport of a trajectory fragment over the time set T t,∆tow centered around t, we compute the time-averaged local Nusselt number Note that the definition of the local Nusselt number is the same as in [33] even for the different thermal boundary conditions once the fields are re-scaled based on Nu.

B. Previous Lagrangian studies
Unlike in the present work where convection is driven by a constant heat flux, our previous Lagrangian analyses have been restricted to the setting of prescribed constant temperatures at the top and bottom boundaries. We summarize our major findings from these past studies in the following.
A first numerical analysis of turbulent superstructures in a large-aspect-ratio 3-dimensional convection cell by means of Lagrangian particle trajectories was carried out in [24] with Pr = 0.7, Ra = 10 5 , and Γ = 16. Coherent sets have been detected from spectral clustering of a weighted and undirected trajectory network from the trajectory points of Lagrangian particles while link weights have been determined by a mean dynamical distance between different particle trajectories. It was demonstrated that the resulting Lagrangian trajectory clusters, which were obtained by a subsequent k-means clustering, are related to bigger patches of the turbulent superstructures from the Eulerian frame of reference. Furthermore, the characteristic Lagrangian time and length scales of the superstructures were found to agree very well with their Eulerian counterparts. For longer time periods beyond the characteristic time, a density-based clustering (DBSCAN [58]) by means of time-averaged Lagrangian pseudo-trajectories was applied. A small coherent subset of the pseudo-trajectories has been obtained consisting of those Lagrangian particles that are trapped for long times in the core of the circulation rolls and are thus not subject to ongoing turbulent dispersion.
In [32], different trajectory-based Lagrangian methods to identify coherence in a simple 2-dimensional Rayleigh-Bénard benchmark problem with Pr = 10, Ra = 10 6 , and Γ = 4 have been compared and found to give consistent results. By introducing the concept of heat particles, we extended the notion of coherent sets to include heat transport. Our findings showed that convective rolls, which were identified by Lagrangian coherent sets, contribute least to the vertical heat transport. Furthermore, [34] demonstrated that also the applications of transfer operator and dynamic Laplacian approaches give consistent results in 2-and 3-dimensional RBC model systems and compare well with the structures identified in previous investigations [24,32].
Lagrangian heat transport in a complex large-aspect ratio 3-dimensional RBC system was first analysed and quantified in [33]. Here we considered a set of three different fluids with Pr ∈ {0.1, 0.7, 7} at Ra = 10 5 and Γ = 16 (one of these runs as in [24]). A diffusion distance, similar in construction to [21], was used to set up the network weight matrices and allowed us after applying a sparse eigenbasis approximation eventually to conclude that Lagrangian coherent sets contribute (independently of the fluid) significantly less -more precisely, up to one third less -to the vertical heat transport than their spatial counterparts.
While our approach in [33] accounted for the slow reorganisation of the turbulent superstructures over longer times only by repeating the analysis of relatively short time windows for several distinct times, [28] was eventually able to address the long-term dynamics of coherent sets in both 2-and 3-dimensional RBC simulations. In particular, we analysed a simple 2-dimensional setup at Pr = 7, Ra = 10 8 , and Γ = 8 as well as a complex 3dimensional setup at Pr = 0.7, Ra = 10 5 , and Γ = 16 (the latter as in [24]). The proposed evolutionary clustering framework eases the mathematical rigor of material transport and allows to monitor the long-term dynamics of Lagrangian coherent sets including their splits or mergers. These events were traced back to the changes of spectral properties of the corresponding network matrices. This study confirmed once more that coherent sets effectively suppress the turbulent heat transport across the domain, with their contribution to the turbulent heat transfer being reduced by about 30% in the twodimensional and 13% in the three-dimensional case, respectively. This evolutionary trajectory clustering framework from [28] represents -together with the physical insights obtained from [10] -the basis for our current study.

IV. LAGRANGIAN CHARACTERISTICS OF THE PRESENT CONVECTION FLOW AND THE INITIAL GRAPH CONSTRUCTION
At the centre of this present study is simulation run Nfs1 from [10] at Pr = 1 and Ra = 10, 430, which we repeat here in a laterally closed Γ = 30 cell (being in contrast to the periodic Γ = 60 domain in [10]). As introduced in sections I and III, this flow is driven by a constant heat flux and, consequently, exhibits a hierarchy of different long-living large-scale flow structures. In more detail, the granules -as the first stage in this hierarchy -offer a time-independent size of Λ G ∼ O (1), whereas the supergranules -as the second stage -grow gradually over time with a final statistically stationary size of Λ SG ≫ 1. In a horizontally periodic domain, Λ SG = Γ in the end. For the present setup, we can quantify the time-independent characteristic wavelength of the granular flow structures Λ G ≈ 5 . . . 6 via the corresponding spectral peak in the vertical velocity field close to the top plane. Although the simulation covers in total 3.5 × 10 3 τ f with a simultaneous advection of 512 2 × 2 Lagrangian particles, we restrict our subsequent analysis to 256 2 Lagrangian trajectories during the first 1 × 10 3 τ f at a temporal resolution of 1τ f . This will capture the major part of the aggregation process, see again Figure 1.
Lagrangian characteristics provide objective physical or data-driven insights into the actual flow at hand, and help eventually to set many of the hyper-parameters related to the clustering procedure. For this reason, we determine the characteristic turnover length and time scales of the flow, Λ to and τ to , respectively. To obtain these characteristics, we use the entire dataset of 512 2 × 2 particles across 3.5 × 10 3 τ f . On the one hand, the characteristic turnover wavelength Λ to is determined from the distribution of the individual turnover lengths λ ito , the latter of which represent four times the horizontal travel distance of every particle between two successive intersections of the midplane z = 0.5. Depending on the chosen typical measure (mean or median) of the resulting probability density function (PDF), see Figure  3 (a), this results in Λ to = 5.8 . . . 8.5 and agrees thus very well with Λ G . On the other hand, the characteristic turnover time τ to is similarly derived from the statistics of the individual turnover times t ito , the latter of which are probed by the particles passing the horizontal planes z = 0.2 and z = 0.8 to complete an entire turnover movement. Again, depending in the chosen typical measure, this yields τ to = 38 . . . 48, see Figure 3 (b).
We start our Lagrangian analysis at 100τ f ≈ 2 . . . 3τ to to allow for a proper distribution of the particles across the entire domain after the initial seeding. We set up our instantaneous adjacency matrices A t (see Section II B) by using K = 26 mutual nearest neighbours. A discusses and compares different choices of parameters for the network construction.
Based on τ to , we use here observation windows of width ∆t ow = 48 ≈ τ to and so the first window is centered around t = 124 - Figure 4 (a) visualizes the corresponding instantaneous temperature field close to the top plane.
As there is no historical information available yet, µ = 1 for this very first observation window and so we compute the leading 100 eigenpairs of only D In accordance with the estimated number of granular flow structures, which is based on the characteristic turnover wavelength Λ to , we pick the gap at k init = 66 (filled triangle) for the analysis in the main text. The evolution of structures corresponding to two further prominent gaps at k init = 15 (filled circle) and k init = 90 (filled square) are analysed and contrasted in B.
In order to decide for one of these potential numbers of coherent features, k init = {15, 66, 90}, we visualize all of them in Figure 4 (b -d). Our past studies [28,33] showed that up-and down-welling fluid separates coherent spatial regions, and so one should expect that the coherent features here are again influenced or separated by such regions -we thus superimpose the regions of strongly down-welling fluid as black contour lines. Together with our estimation of the number of present granular flow structures from above, this discards k init = 15 as a proper choice. In contrast, both k init = {66, 90} seem to extract features that are influenced by the incoherent down-welling fluid. Thus we select the first of these peaks, k init = 66, as the initial number of coherent sets for the upcoming evolutionary clustering. We will contrast the results for all of the k init = {15, 66, 90} in B.
Together with these mostly data-driven or physicsinformed decisions of several hyper-parameters, we move in the next section on to incorporate historical information in the evolutionary spectral clustering approach. First, we will evaluate the coherence of long-living largescale flow structures based on the maximum likelihood of cluster affiliation S max as an output of the evolutionary SEBA post-processing -this quantity is influenced by the incorporation of historical information. Second, we will contrast these results with the current structure of the network weight matrix in terms of the node degree d. Recall here that while this quantity does not take any historical information into account, it simultaneously does not require to solve the expensive eigenvalue problem. Finally, note in particular that high coherence implies where the respective reference values are deduced from ensemble means as outlined below.

V. EVOLUTIONARY SPECTRAL CLUSTERING ANALYSIS
In order to incorporate historical information from the previously extracted coherent sets in each subsequent clustering instance, one needs to shift the observation window by some small amount ∆t s and simultaneously incorporate historical information based on the obliviousness parameter µ -see again eq. (4). Based on a preliminary variation of µ ∈ {1.0, 0.95, 0.8} given a fixed ∆t s = 2, see C, we proceed here in the main text with µ = 0.95. While we continue to use a spectral gap criterion to detect the actually present number of coherent features, we allow this number only to vary by ±1 with each subsequent clustering instance. This allows to automate the entire remaining evolutionary approach. Given these choices of parameters and the features extracted from the initial graph at t = 124, we thus automatically construct W t for t = {126, 128, . . . } and obtain in this way 427 consecutive observation windows.
While the initial number of coherent features is motivated by the number of present granular flow structures as the first stage of the present hierarchy of granular and supergranular flow structures, these features are subsequently influenced by the previous clustering and allowed to evolve over time. Consequently, there may appear mergers, splits, or even establishments of completely new ones. In order to categorize the extracted evolving coherent sets or features with respect to the present hierarchy, we start by visualizing a time series of the instantaneous temperature field (which prominently indicates the gradual supergranule aggregation [10,11]) in Figure 6 (a -c) together with the maximum likelihood of cluster affiliation S max in panels (d -f). Note in case of the latter that regions of more intense color correspond to aspects of the flow with stronger coherence, whereas less intensely colored regions indicate some incoherent background. We superimpose the T = 0.5 contours of the temperature fields from panels (a -c) as solid lines to contrast the gradual growth of the supergranules with the evolving features. Although the locations of the coherent sets seem -to some extent -to be related to the steadily growing supergranules, there is a clear scale separation between them and the growing supergranular flow structures. Hence, our evolutionary algorithm keeps detecting the granular flow structures throughout the evolution of the flow. The implications on the spectrum and its timedependent gap can be found in B and C.
Intending to further investigate the interplay between the identified coherent sets (i.e., granules) and supergranules, we contrast regions of outstanding coherence with the instantaneous temperature field in Figure 7. In particular, we extract those trajectory fragments from  in panel (a). Although we find that the extracted coherent regions do not coincide with the large-scale supergranule pattern, finer temperature ridges -indicating up-and down-welling fluid in the form of detaching thermal plumes -can be found to be located between these coherent regions. This observation underlines once more the accordance of the extracted coherent sets or regions with the granular flow structures from the Eulerian frame of reference. Vice versa, one may similarly state that even the granule patterns leave a footprint in the temperature field. As such ridges are more intense when they correspond to the supergranular pattern, the evolution and location of the granular sets is consequently even influenced by the growing supergranules.
This above evolutionary clustering analysis is conditioned to a certain initial number of coherent features which we selected in the first observation window based on physical arguments about the size of the granular flow structures. In other words, it is fed with additional information on the (so far granular) flow structures of interest. An alternative option to extract coherent features without such a conditioning based on user-provided information is possible via the node degree d.
The node degree provides local information of the structure of the Lagrangian trajectory network weight matrix and represents thus an alternative measure of coherence. Hence, we repeat our analysis of the present coherence from above but in the following based on dthe result is contrasted in Figure 6 (g -i). This time, we find growing but less pronounced structures that seem to correlate with the thermal footprint of the supergranules -in particular, the growing hot and cold spots, i.e., the up-and down-welling parts of the supergranules, seem to be characterized by relatively small node degrees. This first impression is further supported by To quantify these visual impressions, we perform a statistical analysis of the temperature T , node degree d, and maximum likelihood of cluster affiliation S max in Figure  10. Note that this analysis exploits the entire dataset such that the result is independent of the selected time. While panels (a -c) show the individual distributions of the before-mentioned quantities together with their corresponding mean values, their joint probability density functions are provided in panels (d, e). Albeit all PDFs show wide distributions, the shapes of the peaks of the joint PDFs are quite similar. We highlight this pear-shaped peak via white dotted lines that surround the domain around the peak, the latter of which inte- These analyses based on two different measures of coherence (S max and d) show that it is possible to extract both stages of the present hierarchy of long-living largescale flow structures (granules and supergranules, respectively). Our previous studies have shown that supergranules contribute significantly to the heat transport [10], whereas coherent sets at the centres of convection rolls (in the case of turbulent superstructures) reduced the overall heat transport [28,33]. This apparent contradiction asks to study the heat transport for the present config- , which tend to be found more often in regions of extreme (hot or cold) temperatures. Both panels correspond to t = 750, see also Figure 6 (c, f, i).
uration once again depending on the chosen measure of coherence.

VI. CONDITIONAL HEAT TRANSFER ANALYSIS
A statistical analysis of the heat transfer properties along trajectory segments -coherent or incoherent, depending on the choice of the measure of coherence -allows to scrutinize our previous results. In Figure 9, we consider the time-dependent local Nusselt numbers averaged over different (coherent or incoherent) ensembles. The dark line represents the Nusselt number obtained from all particles (independent of their coherence) and offers ⟨Nu local ⟩ N,t = 4.02 ± 0.07 (indicating the mean value and standard deviation). Note that the standard deviation measures here the fluctuations in time of the time-dependent ensemble averages of the local Nusselt number which indicates that the overall Lagrangian heat transport does not change much over time. Based on the chosen measure of coherence and its time-dependent ref- , the evolution of the correspondingly averaged local Nusselt numbers changes.
On the one hand, a separation based on S max results in the two blue curves with ⟨Nu local ⟩ N S, coh ,t = 3.48 ± 0.09 and ⟨Nu local ⟩ N S, inc ,t = 4.53 ± 0.16. The local Nusselt number along trajectories in coherent regions is significantly smaller than those in mixing regions, implying that they participate only weakly in the turbulent heat transfer across the fluid layer. This is in agreement with our previous results in [28].
On the other hand, when we consider a separation based on d (yielding the red curves) we obtain ⟨Nu local ⟩ N d, coh ,t = 4.58 ± 0.12 and ⟨Nu local ⟩ N d, inc ,t = 3.35 ± 0.14. In contrast to our previous Lagrangian studies in [28], trajectories in coherent regions (in the sense of a node degree below average) contribute here significantly more to the vertical heat transport than those from the mixing regions. These results are confirmed when studying the time and ensemble averages which are summarized in Table I or the time-independent conditioned PDFs as shown in Figure 10.
Convective heat transfer across the fluid layer depends crucially on turnover movements of the individual particles. Thus, we finally determine the average turnover wavelengths depending on the chosen ensemble. If this is not conditioned to certain coherent or incoherent trajectories, this yields Λ to = 9.3 ± 6.9 -note that this result differs slightly from the results outlined in section IV as we restrict our clustering analysis to only a subset of the entire dataset. Additionally, we restrict the computation to coherent and incoherent regions based on both S max and d in the same manner as for the local Nusselt number above. The results are included in Table I and show that the turnover wavelengths of coherent sets are smaller (than those of incoherent sets) when coherence is measured with respect to S max , whereas they are larger once coherence is measured with respect to d.
This quantitative study underlines the relevance of the   however, the fraction of the heat transfer due to the supergranule steadily increased until the structure filled the entire cell. The above-mentioned apparent contradiction is thus indeed consistent with our previous Lagrangian and Eulerian analyses.

VII. CONCLUSION AND OUTLOOK
The present work investigates a Rayleigh-Bénard convection flow that is driven by a constant heat flux at the top and bottom boundaries and thus exhibits a hierarchy of flow structures ranging from granules to supergranules. The latter structures evolve out of the small-scale turbulence in a very slow aggregation process that lasts for several thousands of free-fall times and comes to an end only when the finite domain size prohibits any further aggregation. Such hierarchies of structures are found in several much more complex flows in nature, e.g., in convection in the interior of the Sun [4] where they are eventually constrained by rotation [4,11,59] (which is not considered here for simplicity).
In this present study, we obtained time-dependent coherent sets within a network-based framework based on Lagrangian trajectories, i.e., massless particles which follow the surrounding time-dependent velocity field perfectly. We analysed the flow based on two different measures of coherence. On the one hand, we demonstrated that granular coherent structures can be identified by an evolutionary spectral clustering approach and contribute only little to the overall heat transfer. This confirms our previous results [28,33] that have been obtained for the complementary thermal boundary conditions, i.e., for the case of applied constant temperatures at the top and bottom planes. On the other hand, an analysis of the node degree provided insights on the objectively most coherent flow structures present in the flow and extracted the supergranular flow structures instead. These structures were, in accordance with our previous Eulerian studies [10], found to contribute significantly to the overall heat transfer. tract here different stages of the hierarchy. The spectral clustering approach based on the maximum likelihood of cluster affiliation S max extracted the centres of convective granular flow structures (see also [33]), whereas the technically simpler approach based on the node degree d extracted the coherently up-and down-welling parts of the growing supergranular flow structures. Note that such spatial regions outside the centres of convection rolls correspond in the opposite classical case of turbulent superstructures to the incoherent background with high node degrees [28]. To summarize, this highlights that our evolutionary clustering algorithm is capable (i) to analyse the gradual evolution of granular Lagrangian coherent sets, and (ii) to detect the even larger transient supergranular flow structures by alternative measures of coherence. This result is promising since many convective flow phenomena are determined by such a transient evolution, see e.g. the diurnal cycle for atmospheric boundary layers.
In future work, we plan to extend the Lagrangian framework to experimental systems, especially from chemical process engineering, see e.g. [50] for a first study. In chemical reactors, the detailed knowledge of the hierarchy of coherent flow structures may help to identify zones of poor mixing and enable us to tailor the mixing and reaction processes. Mixing and entrainment processes in turbulent clouds are just another possible application of the present framework. Setting up the networks by means of K = 26 instantaneous mutual neighbours is inspired by a cubic grid consisting of 27 grid points, with the centre point having 26 neighbours in this sense. This results in a link density of about 0.003 for the respective network weight matrices, which is well in the range proposed by [60] in the related context of recurrence networks. This choice corresponds to a mean ε = 0.47 and compares well to our previous study in [28]. There, ε = 0.14 was chosen for the same number of trajectories but for a smaller domain of aspect ratio Γ = 16. The resulting difference in the density of the trajectory seeding by a factor 3.5 is directly reflected in the choice of the respective ε values. Figure 11 compares the initial eigenspectra of the spectral clustering for K = 26 and ε = 0.47 and confirms a very good correspondence. This figure includes also the spectra of the following investigations.
In contrast, the spectra for K = 6 (inspired by a 7point stencil) and ε = 0.28 (link density: 0.0007) do not display prominent gaps in the range of the smaller eigenvalues that are of interest in our study (i.e., 50 . . . 100 granular flow structures, see Figure 5) -any clustering would result in these cases in spurious clusters (which is not shown).
The spectra for K = 16 and ε = 0.4 (link density: 0.0018) are very similar to those of K = 26 and ε = 0.47. The corresponding clusters (not shown) compare also well, so that we can conclude that the network construction appears to be robust if sufficiently many neighbours are taken into account. Note that increasing K can be interpreted as increasing diffusion, resulting in smoother clusters and a stronger decay in the spectrum.
In [54], parameters of the spectral approach used in [20] are tuned with respect to optimizing a spectral gap. This is a useful ansatz when only a smaller number of distinct coherent sets are sought. In the setting of this present paper, where we are focusing on many sets with no clear spectral gap in the region of interest, such an ansatz is not feasible.

Appendix B: Evolutionary clustering for different k init
Here we present the results for different initial numbers of coherent features k init for the subsequent evolutionary clustering. In particular, we fix µ = 0.95 as well as the shift of the observation window ∆t s = 2 (equal to main text) and compare k init ∈ {15, 66, 90} corresponding to the major spectral gaps seen and highlighted in Figure  5 in the main text. The evolution of the eigenspectra is shown in Fig. 12. Mergers of coherent sets manifest as upward-moving spectral gaps in terms of the number of eigenvalues, whereas splittings correspond to downwardmoving gaps. The temporal evolution of the cluster indicator S max [X i (t)] is shown in Figure 13.
For k init = 90, the number and positions of the extracted coherent sets appear to be very similar to the case with k init = 66 (as was used in the main text). There are a number of merging events visible in the evolution of the spectrum, see Figure 12 (b, c), but the resulting clusters again represent the granules rather than the supergranules. The statistics of heat transport related to the coherent and incoherent structures is also very similar for k init = 66 and k init = 90, see Table II. For k init = 15, the numbers and positions of coherent sets are less dynamic. We observe spectral signatures of four merging events only, such that at the end of the monitoring 11 structures are remaining -the spectral gap can thus eventually be found after the 11th eigenvalue, see Figure 12 (a). The resulting incoherent regions tend to lie in the centers of the supergranules, where fluid is upand down-welling. However, the differences in heat transport contributions by the coherent and incoherent regions are not as pronounced as for k init = 66 and k init = 90, see again Table II. Consequently, even selecting k init = 15 in the initial setup of the evolutionary clustering algorithm does not allow to directly study the gradual supergranule aggregation based on the maximum likelihood of cluster affiliation S max given our choice of hyper-parameters. Here we study the influence of the obliviousness parameter µ that weights the historical and current costs in the evolutionary clustering framework (see eq. (4)). In particular, we fix k init = 66 as well as the shift of the observation window ∆t s = 2 (equal to main text) and compare µ ∈ {1, 0.95, 0.8}. The resulting eigenspectra obtained from the spectral clustering approach are shown in Figure 14. As µ < 1 introduces a spectral gap that becomes stronger for smaller µ, we observe less changes in the gaps of the spectrum when µ is decreased. This influences the stability of the clusters, as also demonstrated in Figure 15. While for the choice µ = 0.8 more structures persist than for larger µ, the qualitative picture of the evolving clusters appears to be relatively independent of the particular choice of µ. This is also supported by the quantitative heat transfer analysis in Table II which considers the different choices of k init and µ simultaneously.
The results of our parameter studies -both from B and C -indicate that the coherent sets as identified in the main text are robust in the sense that neither the qualitative pictures nor the heat transfer statistics change significantly. We have fixed the length of our observation window ∆t ow = 48 on physical grounds. Small changes in ∆t ow (not shown here) do not change the overall picture, but large changes do as the definition of coherence is crucially tied to the time scale under consideration. Our choice of ∆t s = 2 generates a considerable overlap of the subsequent time windows. We expect that decreasing the fixed shift of the observation window ∆t s will not have a strong impact on the results due to the observed robustness. In contrast, a moderate increases in ∆t s can be balanced by choosing a smaller µ. Further detailed investigations would be possible, but are beyond the scope of this paper.   II: Average local Nusselt numbers for the coherent and incoherent groups of particles measured by S max for different choices of k init and µ. We use ⟨S max ⟩ N,T = S max, ref ± σ S (mean and standard deviation) and identify the trajectory segments X i (t) for which S max [X i (t)] is smaller (incoherent) or larger (coherent) than the reference value S max, ref . The analysis shows that for k init = 15 the difference in the local Nusselt numbers for the coherent and incoherent ensembles is less pronounced as for the larger number of clusters. All parameter choices give very similar quantitative results, so that the method appears to be quite robust.  given a fixed shift of the observation window ∆t s = 2 and k init = 66 initial coherent sets. While we find that some coherent features disappear and others persist, the qualitative picture of evolving clusters is independent of the particular choice of µ. Note that times are chosen to coincide with Figure 6 and that the node degree from Figure 6 (g -i) is independent of µ.