Lagrangian coherent sets in turbulent Rayleigh-B\'enard convection

TL;DR

This paper investigates Lagrangian coherent sets in turbulent Rayleigh-Bénard convection using graph Laplacian-based methods, comparing them with Eulerian analysis, and explores their role in heat transport at high Rayleigh numbers.

Contribution

It introduces and compares three Lagrangian analysis techniques for identifying coherent structures in turbulent convection and relates these to Eulerian structures and heat transfer.

Findings

Lagrangian and Eulerian coherent sets are largely disjoint.

The three methods produce consistent coherent structures.

Coherent sets influence heat transport in turbulent convection.

Abstract

Coherent circulation rolls and their relevance for the turbulent heat transfer in a two-dimensional Rayleigh--B\'{e}nard convection model are analyzed. The flow is in a closed cell of aspect ratio four at a Rayleigh number ${\rm Ra}=10^6$ and at a Prandtl number ${\rm Pr}=10$. Three different Lagrangian analysis techniques based on graph Laplacians -- distance spectral trajectory clustering, time-averaged diffusion maps and finite-element based dynamic Laplacian discretization -- are used to monitor the turbulent fields along trajectories of massless Lagrangian particles in the evolving turbulent convection flow. The three methods are compared to each other and the obtained coherent sets are related to results from an analysis in the Eulerian frame of reference. We show that the results of these methods agree with each other and that Lagrangian and Eulerian coherent sets form basically a disjoint union of the flow domain. Additionally, a windowed time-averaging of variable interval length is performed to study the degree of coherence as a function of this additional coarse graining which removes small-scale fluctuations that cause trajectories to disperse quickly. Finally, the coherent set framework is extended to study heat transport.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}$$

$$\frac{\partial \theta }{\partial t}+\mathbf{u}\cdot \nabla \theta$$

$$\nabla \cdot \mathbf{u}$$

$$T(x,z,t)=T_{\mathrm{bottom}}-\frac{\Delta T}{H}z+\theta (x,z,t),$$

$$\mathrm{Ra}$$

$$\mathrm{Pr}$$

$$u_x(x,z,t)$$

$$u_z(x,z,t)$$

$$\theta (x,z,t)$$

$$\frac{d{\mathbf{x}}_i}{dt}=\mathbf{u}({\mathbf{x}}_i,t).$$

$${\mathbf{\Phi }}_{t_0}^t:\mathbf{x}(t_0)\mapsto \mathbf{x}(t)={\mathbf{\Phi }}_{t_0}^t(\mathbf{x}(t_0),t_0)$$

$$\rho (\mathbb{A},\mathbb{B})=\frac{\mu (\mathbb{A}\cap {\mathbf{\Phi }}^{-1}(\mathbb{B}))}{\mu (\mathbb{A})}+\frac{\mu ({\mathbb{A}}^C\cap {\mathbf{\Phi }}^{-1}({\mathbb{B}}^C))}{\mu ({\mathbb{A}}^C)},$$

$$\left\{{\mathbf{x}}_i(t_k)\in {\mathbb{R}}^d|i\in \{1,\dots ,N\},k\in \{0,\dots ,\mathcal{T}\}\right\}\subset {\mathbb{R}}^d,$$

$$Ly=\lambda Dy$$

$$Q^{\mathrm{nw}}y=\lambda y,\text{where }{Q}^{\mathrm{nw}}=D^{-1}{L}^{\prime }.$$

$$A_{ij}=\left\{\begin{array}[]{ll}1,&\|{\bm{x}}_{i}(t_{k})-{\bm{x}}_{j}(t_{k})\|\leq\varepsilon\\ &\text{ for at least one }k\in\{0,\ldots,\mathcal{T}\}\text{ and }i\neq j,\\ 0,&\text{otherwise}.\end{array}\right.$$

$$D_{ij}=\left\{\begin{array}[]{ll}\sum_{j^{\prime}}A_{ij^{\prime}},&i=j\\ 0,&\text{otherwise}.\end{array}\right.$$

$$\overline{d}({\mathbf{x}}_i,{\mathbf{x}}_j)=\frac{1}{\mathcal{T}+1}\sum _{k=0}^{\mathcal{T}}\Vert {\mathbf{x}}_i(t_k)-{\mathbf{x}}_j(t_k)\Vert ,$$

$$k_{\delta }(\mathbf{x},\mathbf{y})=\exp \left(-\frac{\Vert \mathbf{x}-\mathbf{y}{\Vert }^2}{\delta }\right)\text{mathds}{1}_{[0,r]},$$

$$P^{\mathrm{dm}}=\frac{1}{\mathcal{T}+1}\sum _{k=0}^{\mathcal{T}}{P}_{t_k},$$

$$Q^{\mathrm{dm}}=\frac{1}{\delta}\big{(}P^{\mathrm{dm}}-I\big{)},$$

$$K_{t_{k},ij}=k_{\delta}\big{(}{\bm{x}}_{i}(t_{k}),{\bm{x}}_{j}(t_{k})\big{)}.$$

$$P_{t_k}=D_{t_k}^{-1}{K}_{t_k},$$

$$P^{\mathrm{dm}}=\frac{1}{\mathcal{T}+1}\sum _{k=0}^{\mathcal{T}}{P}_{t_k},$$

$$\frac{\partial S}{\partial t}={\tilde{\kappa }}_L{\mathbf{\Delta }}^{t}S,$$

$$\mathcal{D}(\mathbf{x},t)=D{\mathbf{\Phi }}_{t_0}^t(\mathbf{x}{)}^{-1}D{\mathbf{\Phi }}_{t_0}^t(\mathbf{x}{)}^{-\top },$$

$$\bar{\mathbf{\Delta }}=\frac{1}{\mathcal{T}+1}\sum _{k=0}^{\mathcal{T}}{\mathbf{\Delta }}^{t_k},$$

$$K_{ij}^k=-{\int }_{\Omega }⟨\nabla {\phi }_i^k,\nabla {\phi }_j^k⟩dx.$$

$$M_{ij}={\int }_{\Omega }{\phi }_i^0\cdot {\phi }_j^{0}dx$$

$${\delta }_T=\frac{1}{2\mathrm{Nu}}<\epsilon <\frac{r_{roll}}{2}.$$

$$\varepsilon\in\left\{\begin{array}[]{ll}\left[0.03,0.07\right]&:\;\tau=2020t_{f}\,,\\ \left[0.085,0.25\right]&:\;\tau=2200t_{f}\,.\end{array}\right.$$

$${\hat{\mathbf{x}}}_i(t_{\hat{k}})=\frac{1}{{\mathcal{T}}_W}\sum _{j=\hat{k}}^{\hat{k}+{\mathcal{T}}_W-1}{\mathbf{x}}_i(t_j),$$

$$\frac{\partial T}{\partial t}=\frac{1}{\sqrt{\mathrm{PrRa}}}{\nabla }^{2}T-\nabla \cdot (\mathbf{u}T)=:-\nabla \cdot \mathbf{\Psi },$$

Full text

I Introduction

Thermal turbulent convection acts as one essential driving mechanism in many turbulent flows in nature spanning a wide range of examples starting from stellar interiors Miesch2005 via planetary cores Christensen2006 to atmospheric motion Stevens2005 and transport dynamics in lakes and ponds Bouffard2019 . An idealized model of thermal convection is Rayleigh–Bénard convection (RBC), in which a fluid layer placed between two solid horizontal plates is uniformly heated from below and cooled from above Chilla:EPJE2012 . This particular setting contains already many of the properties which can be observed in natural flows. One is the formation of large-scale coherent patterns when RBC is investigated in horizontally extended domains Hartlep2003 ; Hardenberg2008 ; Bailon2010 ; Emran2015 ; Stevens2018 ; Pandey:NatCommun2018 . These coherent sets, which have been detected in the Eulerian frame of reference, are termed turbulent superstructures as the characteristic horizontal scale extends the height of the convection layer. In thermal convection flows, they consist of convection rolls and cells that may, however be concealed in instantaneous velocity fields by turbulent fluctuations. Large-scale circulations in Rayleigh–Bénard convection exist also in smaller domains or cells and have been analyzed for example by proper orthogonal decomposition (POD) Bailon2011 ; Podvin2012 ; Verma2017 (see also the Appendix of Verma Verma2018 for POD in RBC). In the Eulerian frame, large-scale patterns show up prominently either after time averaging or as the primary POD modes, for both, temperature and velocity fields. This is illustrated in Figure 1 where two coherent circulation rolls are present with narrow regions of upwelling hot and downwelling cold fluid.

At the core of the data-driven analysis in the Lagrangian frame of reference is the concept of a coherent set FrLlSa10 ; Froyland2013 ; Allshouse2015 ; Karrasch2017 , a region in the fluid volume that only weakly mixes with its surrounding and which often stays regularly shaped (non-filamented) under the evolution by the flow. Such regions can be determined in two ways, by set-oriented DellnitzFroylandJunge2001 ; DellnitzJunge2002 or manifold-based methods (see Allshouse2015 ; Hadjighasem2017 for recent reviews of both concepts). The manifold-based approach comprises Lagrangian coherent structures (LCS), i.e., minimal curves in two dimensions and surfaces in the three-dimensional case that enclose coherent sets Haller2015 . This framework was extended recently to include weak diffusion across the manifold Haller2018 .

Coherent sets were originally introduced based on transfer operators FrLlSa10 ; Froyland2013 . These are linear operators that evolve densities under the action of the flow. Coherent sets can be identified from the leading singular functions of this operator. More recently, in Ref. Froyland2015 they have been characterized as sets which possess a minimal boundary-to-volume ratio for the entire flow duration. Different approaches have been introduced recently that make use of spatio-temporal clustering algorithms applied to Lagrangian trajectory data Froyland_Padberg_2015 ; Hadjighasem2016 ; Banisch2017 ; Schlueter2017 ; Padberg2017 ; Schneide2018 . These algorithms aim at identifying coherent sets as groups of trajectories that remain close to each other in the time interval under investigation.

In this work, we will focus on the latter of these Lagrangian approaches. The three methods that we are going to apply characterize coherence via regularized linear operators that are directly approximated on the basis of the Lagrangian trajectory data in the convection flow. These are the (i) transfer operator which is regularized by a diffusion kernel Banisch2017 , the (ii) graph Laplacian operator that characterizes a network of Lagrangian trajectories in the flow Padberg2017 and the (iii) dynamic Laplacian which characterizes sets with minimal averaged boundary-to-volume ratio Froyland2015 ; FroylandJunge2018 or equivalently almost invariant sets of a time-dependent heat flow Karrasch2017 . In all cases, gaps in the discrete eigenvalue spectrum of the operator under consideration give an indication of the number of coherent sets. As will be seen, all three Lagrangian methods detect the same core regions of the large-scale circulation rolls as the coherent sets in which fluid particles remain together for the longest time. With progressing time these coherent sets get increasingly smaller in their spatial extent as expected for turbulent flows.

A second aspect of this work is therefore to extend these Lagrangian concepts and to perform the analysis on data which are averaged in time over a window of variable length. Similar to the Eulerian studies which were mentioned at the beginning of this introduction, we want to investigate coherence as a function of this additional coarse graining which removes small-scale fluctuations in the flow that typically cause a fast separation of Lagrangian trajectories that are initially close together.

A final aspect of this work is to adapt the presented analysis directly to the transfer of heat. The presented Lagrangian methods can then be used to investigate heat coherence in RBC.

Here, we study RBC in a two-dimensional closed box of aspect ratio four. Note that the large- and small-scale quantities show similar scalings in two- and three-dimensional RBC Schmalzl:EPL2004 ; Poel:JFM2013 ; Pandey:Pramana2016 for large Prandtl numbers. Therefore, very-long-time temporal evolution of the convective flow configurations has been studied in two-dimensional settings Petschel:PRE2011 ; Pandey:PRE2018 . The objective of this work is to take such a simple two-dimensional turbulent flow at a moderate Rayleigh number and to demonstrate and compare the Lagrangian concepts and ideas.

In Sec. II we introduce the numerical model and the data set. Section III gives a further motivation for coherence. The Lagrangian methods which we will compare are introduced in Sec. IV and applied to the convection data in Sec V. Heat coherence is discussed in Sec VI and we conclude in Sec. VII.

II Rayleigh–Bénard convection flow

Conservation of mass, momentum, and internal energy lead to equations which govern the dynamics of RBC. In the Boussinesq approximation Chilla:EPJE2012 they are given in a non-dimensional form by

[TABLE]

where ${\bm{u}}=(u_{x},u_{z})$, $\theta$, and $p$ are the velocity, temperature deviation, and pressure fluctuation fields, respectively. Note that the temperature fluctuation from the linear conductive profile is related to the total temperature field $T$ as

[TABLE]

where $T_{\rm bottom}$ is the temperature at the bottom plate. Equation (4) is given here in physical units. Equations (1–3) were nondimensionalized using the height of the simulation domain $H$ of as the characteristic length scale, the free-fall velocity $u_{f}=\sqrt{\alpha g\Delta TH}$ as the characteristic velocity, and the temperature difference $\Delta T$ between the top and bottom plates as the characteristic temperature. The main governing parameters of RBC are the Rayleigh number ${\rm Ra}$ and the Prandtl number ${\rm Pr}$. The Rayleigh number signifies the strength of thermal driving force compared to dissipative forces, and the Prandtl number is the ratio of the kinematic viscosity and thermal diffusivity of the fluid. They are defined as

[TABLE]

where $\alpha,\nu,\kappa$ are the thermal expansion coefficient, the kinematic viscosity, and the thermal diffusivity of the fluid, respectively. The acceleration due to gravity $g$ points downwards.

We assume that the fields $u_{x},u_{z},\theta\in H^{1}(\Omega\times[0,\tau])$ with $\Omega=[0,L_{x}]\times[0,H]$ and total integration time $\tau$. Here, $H^{1}$ is the Sobolev space of square integrable functions with square integrable derivatives. Equations (1–3) are solved using a pseudospectral solver Tarang Verma:Pramana2013 in a two-dimensional box of aspect ratio $\Gamma=L_{x}/H=4$. Stress-free (or free-slip) boundary conditions for the velocity field are employed at all the walls. This implies that the corresponding normal velocity component and the normal derivative of the tangential velocity component vanish to zero, respectively. For the temperature field, isothermal (adiabatic) boundary conditions are applied in the vertical (horizontal) direction. To satisfy these boundary conditions, the temperature fluctuation and velocity components are expanded in sine and cosine basis functions. This results to

[TABLE]

where ${\bm{k}}=(k_{x},k_{z})$ is the wave vector Pandey:PRE2018 with $k_{x}=m\pi/L_{x}$ and $k_{z}=n\pi/H$; $m,n$ being integers. We perform direct numerical simulation for $\mathrm{Pr}=10,$ $\mathrm{Ra}=10^{6}$, and $\Gamma=4$ using $513\times 129$ uniformly spaced grid points. The time advancement is done using fourth-order Runge–Kutta method (RK4), and the fields are de-aliased using the 2/3 rule. We refer to Pandey:PRE2018 ; Verma:POF2015 ; Verma:Pramana2013 for more details on the numerical simulations.

The presented analyses require Lagrangian particle tracks which are evaluated together with the turbulent flow. Each individual tracer particle $i$ is advected in the velocity field corresponding to

[TABLE]

We simulate $i=1\dots N$ particle trajectories with $N=5000$. Time integration is done again by the RK4 method. The interpolation of the velocity field on the particle position applies cubic splines.

We start our simulation with random noise for velocity and temperature fields, ${\bm{u}}$ and $\theta$, as the initial condition and continue until a statistically steady state after 2000 free-fall times $t_{f}$ is reached. Here, $t_{f}=H/u_{f}$. The time-averaged flow structure exhibits a pair of counter-rotating circulation rolls as shown in the bottom panel of Figure 1. Hot fluid rises in the central region whereas cold fluid falls near the sidewalls.

The velocity and temperature fields at all the grid points are written out every $0.1\,t_{f}$.

III Lagrangian coherence in a turbulent flow

In the dynamical systems perspective, we consider the turbulent convection flow as a mapping in the state or phase space $\Omega\subset\mathbb{R}^{2}$. Let

[TABLE]

denote the flow map, which takes fluid particles from their initial location ${\bm{x}}_{0}={\bm{x}}(t_{0})$ at time $t_{0}$ to their spatial location ${\bm{x}}(t)$ at time $t$ in correspondence with the velocity field ${\bm{u}}$ of the RBC flow. This mapping is given by the differential equations (10) such that ${\bm{\Phi}}_{t_{0}}^{t}({\bm{x}}(t_{0}),t_{0})$ solves the corresponding initial value problem. The advecting flow is simultaneously determined by solving the Boussinesq equations for the RBC flow. The phase space $\Omega$ is equipped with a reference measure $\mu$ and a sequence of non-singular time-dependent flow maps ${\bm{\Phi}}^{t_{0}+\Delta t}_{t_{0}},\ldots,{\bm{\Phi}}^{t_{0}+n\Delta t}_{t_{0}+(n-1)\Delta t}$ for $n$ time steps $\Delta t$. We construct a single flow map via a successive application of flow maps over smaller time steps $\Delta t$, namely as ${\bm{\Phi}}={\bm{\Phi}}^{t_{0}+\Delta t}_{t_{0}}\circ\ldots\circ{\bm{\Phi}}^{t_{0}+n\Delta t}_{t_{0}+(n-1)\Delta t}$.

An overarching goal is to detect and locate slow mixing dynamical structures. These structures should be macroscopic in size and by ”slow mixing” we have in mind a geometric mixing rate that is slower than $1/\Lambda$ where $\Lambda$ is the largest positive Lyapunov exponent that measures the exponential separation of initially infinitesimally distant neighboring particles. Thus, such slow mixing cannot be explained by local stretching, but is instead due to the way in which the dynamics acts globally. Following Froyland2015 , we wish to partition the state space $\Omega=\mathbb{A}\cup\mathbb{A}^{C}$ into a disjoint union of $\mathbb{A}$ and its complement $\mathbb{A}^{C}$ at initial time and $\Omega=\mathbb{B}\cup\mathbb{B}^{C}$ at final time such that we optimize the coherence ratio which can be thought of as

[TABLE]

and that quantifies how much of $\mathbb{A}$ ends up in $\mathbb{B}$ after applying the flow map ${\bm{\Phi}}$ (augmented by a small random perturbation, to be precise) with respect to the reference volume of $\mathbb{A}$ and similarly for $\mathbb{A}^{C}$ and $\mathbb{B}^{C}$. In a nutshell, Eq. (12) quantifies in our example the area content of the subset that stays connected from a Lagrangian point of view.

We seek for an optimal 2-partition into slow mixing (or coherent) sets and the remainder which becomes rapidly well-mixed under the action of convective turbulence. Coherent means now that these sets are almost-invariant under the combined forward–backward dynamics of the flow map ${\bm{\Phi}}$ (including a small amount of diffusion for regularization). The exact version of (12) and the connection to the transfer operator $\mathcal{P}$ and its singular values or functions can be found in Froyland2013 .

Central to manifold-based concepts of coherence are LCS Haller2015 . These are material surfaces that extremise a certain stretching or shearing quantity, such as measured by finite-time Lyapunov exponents (FTLEs) Haller2001 . Ridges in the FTLE field from a forward time computation highlight repelling LCS. In Figure 2 we show the forward time FTLE fields for the two different time windows considered in section V.3. While for the short time interval (top) isolated ridges can be observed, this is no longer the case for the long time interval (bottom). This makes it difficult to unambiguously define coherent sets from the FTLE ridges (see also Schneide2018 for related observations), and therefore, in the present work, we focus on set-based approaches. In Figure 2 we also show the results of a set-based computation (in this case using the minimum distance spectral trajectory clustering method introduced in section IV.1), which identifies coherent sets as the core regions of the large-scale circulation rolls that appear to be characterized by low FTLE values. We come back to this point in Sec. V.3.

IV Graph Laplacian–based coherent structure detection

All algorithms to be introduced work with the following set of Lagrangian data. Let the positions of $N$ particles ${\bm{x}}_{i}$ at $\mathcal{T}+1$ time instances $t_{k}$ be given, i.e., the trajectory data consists of

[TABLE]

where ${\bm{x}}_{i}(t_{k+1})={\bm{\Phi}}_{t_{k}}^{t_{k+1}}{\bm{x}}_{i}(t_{k})$, for $i\in\{1,\ldots,N\}$ and $k=0,\ldots,\mathcal{T}-1$. Trajectories in a coherent set will stay close to each other over a long time in contrast to trajectories not belonging to this specific coherent set. Diverging trajectories indicate filamentation of a set which has a large diffusive outward transport with respect to the dynamics. Based on this notion the algorithms evaluate the $\varepsilon$-neighborhood of the trajectories using distinct models of the diffusion process. If possible, this is followed by the construction of a rate matrix $Q$, which will be introduced in the following sections. The solution of an eigenvalue problem yields eigenvalues $\lambda_{\ell}$ which satisfy $0=\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}$. The eigenvectors corresponding to the dominant eigenvalues (i.e. close to zero) are used to cluster the trajectory data into coherent sets.

IV.1 Distance spectral trajectory clustering

The idea of the network–based analysis, published in Padberg2017 , is to interpret each Lagrangian trajectory $\{{\bm{x}}_{i}(t_{k})\}_{k}$ as node of a network consisting of $N$ nodes. A link between two nodes $\{{\bm{x}}_{i}(t_{k})\}_{k}$ and $\{{\bm{x}}_{j}(t_{k})\}_{k}$ is created if and only if the minimum distance of the two trajectories for at least one time instance $t_{k}$ is smaller than a pre-specified cutoff radius $\varepsilon$. Thereby, network measures as, e.g., the node degree or local clustering coefficient can be used to distinguish coherent sets from incoherent flow Padberg2009 ; Donner2009 ; Banisch2019 . In order to partition the network into independent sets we use a version of the normalized cut approach ShiMalik2000 , which is based on spectral graph theory. This approach aims at maximizing the intra-cluster connectivity and simultaneously minimizing the inter-cluster connectivity. An approximate solution of the normalized cut problem can be achieved solving the generalized eigenvalue problem ShiMalik2000

[TABLE]

where $L=D-A$ is the non-normalized graph Laplacian, $D$ is the degree matrix and $A$ is the adjacency matrix. Here, we vary the approach in the way that we use $L^{\prime}=A-D$ and, for the sake of direct comparability with the second method, we state (14) in the equivalent form 111Assuming $D$ has no zeros on the diagonal, i.e., every trajectory “meets” another at least once. If not, that is there is some $i$ with $D_{ii}=0$, we can either delete the trajectory from the set completely, or if we want to keep it as an independent entity not connected to any other, we can set the $i$-th row of $Q^{\mathrm{nw}}$ equal to the $i$-th row of the identity matrix. The latter step produces a cluster containing only the $i$-th trajectory.,

[TABLE]

Note that the eigenvalues of (14) and (15) differ only in their signs. The multiplicity $m$ of the first eigenvalue $\lambda_{1}=0$ equals the number of connected components in the network. We construct the network such that the zero eigenvalue is simple. Then the number of eigenvalues $\lambda_{2},\ldots\lambda_{m}$ close to [math] determines the number of weakly linked subgraphs. The eigenvectors corresponding to these eigenvalues are used to extract $m$ clusters. In this paper, this post-processing is done by the k-means algorithm Lloyd1982 . The corresponding pseudo-code is given in Algorithm 1.

We note that the adjacency matrix in Algorithm 1 could be constructed as a weighted matrix instead of the current binary one. One possibility is to measure the average distance between two trajectories

[TABLE]

and set $A_{ij}=\overline{d}({\bm{x}}_{i},{\bm{x}}_{j})^{-1}$, $i\neq j$ as was suggested in Hadjighasem2016 . In general, relaxing the binary structure of $A$ by introducing weights is a refinement of the dynamical information contained in $A$, and is used in the method discussed next.

IV.2 Time-averaged diffusion maps

The theory of diffusion maps introduced by Coifman and Lafon in Coifman2006 has been successfully applied to a variety of nonlinear dimensionality reduction problems. In Ref. Banisch2017 the framework of diffusion maps is used to analyze transport in dynamical systems and to find coherent sets solely based on possibly sparse or incomplete Lagrangian trajectory data. The idea is to introduce a diffusion process on the data points and detect points that can be reached more easily from one another by the diffusion process. This is done via the eigenvectors of this diffusion process (or operator) which provide intrinsic coordinates of the data set.

The time-averaged diffusion map algorithm from Banisch2017 (called “space-time diffusion map” therein) proceeds as follows. Given the trajectory data as in (13), the algorithm looks for tight bundles of trajectories. To achieve this, a diffusion map matrix $P_{t_{k}}$ is constructed at every time instance $t_{k}$ by row-normalizing a similarity matrix $K_{t_{k}}$. This matrix is constructed using a rotation invariant kernel which is given by

[TABLE]

with the characteristic function $\mathds{1}_{[0,r]}$ introducing a cutoff at some radius $r$. Hence, the similarity matrix $(K_{t_{k}})_{ij}=k_{\delta}({\bm{x}}_{i}(t_{k}),{\bm{x}}_{j}(t_{k}))$ is only dependent on the Euclidean distances $\|{\bm{x}}_{i}(t_{k})-{\bm{x}}_{j}(t_{k})\|$ in $\mathbb{R}^{2}$ between the point pairs ${\bm{x}}_{i}(t_{k}),{\bm{x}}_{j}(t_{k})$. The parameter $\delta$ can be seen as the strength (or duration) of the diffusion, and characterizes what will be considered close; hence its square root is a length scale of spatial resolution. In practice $\delta$ together with the kernel function will determine a cutoff radius $r$, beyond which the similarity is negligibly small (here $r=\sqrt{2\delta}$). This is similar to $\varepsilon$ in Sec. IV.1, and allows for efficient numerical computation. Now, $P_{t_{k}}$ corresponds to a random walk (Markov chain) constructed on all existing data points (or particle positions) at $t_{k}$ with two points getting a higher jump probability the closer they are. Using time averaging,

[TABLE]

we obtain the space-time diffusion matrix $P^{\mathrm{dm}}$ which describes a Markov chain on all trajectories. We note, that averaging the matrices $P_{t_{k}}$ yields different results to using diffusion maps with the averaged distances $\overline{d}({\bm{x}}_{i},{\bm{x}}_{j})$ from above; and has a different interpretation (cf. Banisch2017 ).

Again, for the reason of comparability, we will work with the rate matrix

[TABLE]

where $I$ denotes the identity matrix. By construction, the eigenvalues $\lambda_{\ell}$ of $Q^{\mathrm{dm}}$ satisfy $0=\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}$. Just as in Sec. IV.1, the dominant eigenvectors of $Q^{\mathrm{dm}}$ are now used to cluster the trajectory data into coherent sets.

We remark that the algorithm, as introduced here in Algorithm 2, corresponds to a particular choice of scaling in diffusion maps, the so-called $\alpha=0$ case Coifman2006 ; Banisch2017 . Originally in ref. Banisch2017 a different scaling was used. We refer to that work for the details, and note only that for uniformly distributed data there is no practical difference. With the correct scaling, it is shown in Ref. (Banisch2017, , Thm. 3 and eq. (21)) that if the number of trajectories goes to infinity, $Q^{\mathrm{dm}}$ converges to the so-called dynamic Laplacian Froyland2015 ; FroylandKwok2017 ; Karrasch2017 whose eigenvectors characterize Lagrangian coherent sets. This makes the time-averaged diffusion-map method a consistent data-based approach for finding coherent sets.

IV.3 Discrete dynamic Laplacian

A third approach to the detection of Lagrangian coherent sets from sparse and possibly incomplete Lagrangian trajectory data has been developed in FroylandJunge2018 . It is based on the geometric idea Froyland2015 that Lagrangian coherent sets can be characterized by a small boundary-to-volume ratio: Whenever the length of the boundary of some advecting set $\bm{\Phi}_{t_{0}}^{t}(\mathbb{A})$, $\mathbb{A}\subset\Omega$, is small in relation to its area consistently for all times $t$, diffusive transport of some passive scalar over its boundary (induced by small random perturbations to $\bm{\Phi}_{t_{0}}^{t}$) will be small. Consequently, the coherence ratio (12) of the pair $\mathbb{A},\bm{\Phi}_{t_{0}}^{t}({\mathbb{A}})$ will be large even in the presence of small perturbations to $\bm{\Phi}_{t_{0}}^{t}$.

Equivalently, cf. Karrasch2017 , we can characterize a Lagrangian coherent set as a material set which is almost-invariant under the flow of the Lagrangian diffusion equation Thiffeault2004 for some diffusive scalar quantity $S$. In Lagrangian coordinates, this equation is given by Karrasch2017

[TABLE]

with ${\bm{\Delta}}^{t}S={\bm{\nabla}}\cdot\left({\cal D}({\bm{x}},t){\bm{\nabla}}S\right)$ and the dimensionless $\tilde{\kappa}_{L}\ll~{}1$. Here, the advection by the flow map that deforms a material set has been encoded in a Lagrangian eddy diffusivity (i.e. the inverse Cauchy-Green strain tensor) which is given by

[TABLE]

where $D{\bm{\Phi}}_{t_{0}}^{t}$ is the Jacobian of the flow map.

In order to compute coherent sets via this approach, we first need to remove the time-dependence from the diffusion operator ${\bm{\Delta}}^{t}$. This can be achieved by time-averaging, i.e. by considering the operator

[TABLE]

called the dynamic Laplacian Froyland2015 .

In a second step, one considers the eigenproblem $\bar{\bm{\Delta}}v~{}=~{}\lambda v$ with appropriate boundary conditions (here, we have used Neumann boundary conditions). Lagrangian coherent sets are then given by sublevel sets of the eigenvectors at the leading eigenvalues of $\bar{\bm{\Delta}}$, i.e. those closest to [math] (cf. Dellnitz1999 ; DellnitzFroylandJunge2001 ; Froyland2015 for more details).

The dynamic Laplacian $\bar{\bm{\Delta}}$ in this eigenproblem can be discretized by standard finite element methods (FEM), leading to the generalized eigenproblem $\bar{K}v=\lambda Mv$ with eigenvalues $0=\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}$. The stiffness matrix $\bar{K}$ is the average $\bar{K}=\frac{1}{\mathcal{T}+1}\sum_{k=0}^{\mathcal{T}}K^{k}$ of the stiffness matrices at each time $t_{k}$,

[TABLE]

Here, the functions $\varphi^{k}_{i}$ are the finite element basis functions on a triangulation of the data points ${\bm{x}}_{i}(t_{k})$. For the mass matrix $M$, it suffices to compute

[TABLE]

only at the initial time $t_{0}$.

A typical finite element in step 1 of Algorithm 3 is the linear triangular Lagrange element, i.e., the mesh in step 3 consists of triangles and the basis functions are the piecewise linear nodal basis functions. This computational mesh can be constructed as the Delaunay triangulation of the data points $\{{\bm{x}}_{i}(t_{k}),i=1,\ldots,N\}$ at each time step. Whenever the data set has a very irregular hull, an $\alpha$-shape will be more appropriate. Typically, $\alpha$ is then chosen minimal such that the triangulation is connected (cf. alphashape in Matlab). We refer to Ref. FroylandJunge2018 and to the documentation of the packages FEMDL (Matlab) and CoherentStructures.jl (Julia), which are available from github, for more details and examples.

The triangulation induces a network of trajectories for each time step. There is an interpretation of this discretization of the dynamic Laplacian from a graph-Laplacian perspective, detailed in Appendix A.1.

Finally, we want to state that there is no direct link to some rate matrix $Q$ as in the other two approaches since the inverse of $M$ will be in general not sparse.

V Comparison of different Lagrangian methods

In order to compare the methods introduced above, we choose different perspectives. We compare the approaches first on a theoretical level (V.1) and regarding their structure (V.2). Finally we apply all these three methods to the same trajectory data set introduced in Sec. II and compare their results (V.3).

V.1 Theoretical level

We start by noting that $Q^{\mathrm{nw}}$ and $Q^{\mathrm{dm}}$ are so-called rate matrices. A rate matrix $Q$ has the properties $\sum_{j}Q_{ij}=0$ for all $i$ and $Q_{ij}\geq 0$ for $i\neq j$. It defines a time-continuous Markov chain in the following sense. The process being in state $i$ at current time, (i) will stay in $i$ for a random amount of time $\tau$, where $\tau$ is an exponentially distributed random variable with rate $-Q_{ii}$, 222A random variable is exponentially distributed with rate $\rho$, if the probability that $\tau>t$ for any $t>0$ is $e^{-\rho t}$. Equivalently, the probability density of $\tau$ is given by $\smash{f_{\tau}(t)=\rho^{-1}e^{-\rho t}}$. and (ii) then will jump to the next state $j$ randomly (and independently of $\tau$) with probability $Q_{ij}/|Q_{ii}|$. Thus, the larger the absolute value of $Q_{ii}$, the faster the jump occurs on average. The quantity $1/|Q_{ii}|$ is called the (average/expected) holding time of state $i$.

We refer to Appendix A.2 for further technical details, and only note here briefly that the jump processes $Q^{\mathrm{nw}}$ and $Q^{\mathrm{dm}}$ introduce jumps of mean length $\mathcal{O}(\varepsilon)$ and $\mathcal{O}(\sqrt{\delta})$, respectively, and this distance governs the finest scales they can resolve. Coherent sets below these scales can not be detected. The FEM-based discretization of the dynamic Laplacian has no explicit scale parameter, however the average diameter of the triangulation can implicitly be viewed as such.

V.2 Structure

The descriptions in Sec. IV suggest that the algorithms construct similar objects and have partially even similar steps. All three algorithms are by their very nature frame independent (as they only use mutual distances of trajectories), fast, and construct sparse matrices that encode the space-time behavior. They require the following three parameters as user input:

(1) A subset $\Theta\subset\{0,\ldots,\mathcal{T}\}$ of all time instances at disposal, Not every coherent set might be present for the entire observation range of a fully non-autonomous flow, thus a natural choice is to restrict the time interval of consideration to $\Theta=\{k_{0},\ldots,k_{1}\}$, where $0\leq k_{0}<k_{1}\leq\mathcal{T}$. Also, if the sampling time step $\max_{k}|t_{k}-t_{k-1}|$ is small compared to typical dynamic time scales, the data can be subsampled with respect to time without essentially altering the results.

(2) Proximity parameters $\delta$ and $\varepsilon$ (not for the FEM approach). These govern the minimal spatial scales on which the methods can detect coherence. In the FEM method this scale will implicitly be defined by the size of the elements resulting from the triangulation. For uniformly distributed Lagrangian particles the diameter of the resulting triangulation can be interpreted as the implicit closeness parameter.

(3) The number of clusters. Gaps in the spectrum of $Q^{\mathrm{dm}}$, $Q^{\mathrm{nw}}$ or $(\bar{K},M)$ can indicate natural choices for the number of clusters. Some arguments for this approach are collected in Denner2017 .

All methods discussed here are highly suitable for sparse and incomplete data (see Banisch2017 , in particular section V therein, Padberg2017 and FroylandJunge2018 ).

V.3 Numerical

In the following, we discuss how the three methods compare to each other when applied to the turbulent convection flow data introduced in Sec. II for times $t>2000t_{f}$. In view to the Lagrangian coherence, we will look at a short time interval of ${\cal T}=200$ steps and a long time interval of ${\cal T}=2000$ steps with $0.1t_{f}$ per step. We thus start at $t=2000.0t_{f}$ and end at $\tau=2020.0t_{f}$ and $\tau=2200.0t_{f}$, respectively. These time intervals are chosen since the average circulation (or eddy turnover) time of one roll was found to be $t\approx 20t_{f}$, calculated using the maximum circumference of the roll covering half of the box and the root mean square velocity (see Pandey:NatCommun2018 for more details). Thus, the shorter (longer) time interval approximately corresponds to $1$ circulation ($10$ circulations). We initialize $N=5000$ uniformly distributed particles in order to analyze mass transport. The particles are advected using the snapshot files of the velocity field and for comparability all three algorithms are applied to the same trajectory data set.

At this point we would like to seize the opportunity to explain some options for the choice of the similarity parameters $\delta$ and $\varepsilon$. Some general criteria for the choice of $\delta$ and $\varepsilon$ are (i) sparsity of $Q^{\mathrm{dm}}$ and $Q^{\mathrm{nw}}$ to achieve for example $5\%$ non zero entries, (ii) a stable spectrum which implies that for a series of $\delta$, $\varepsilon$ the dominant spectrum of $Q^{\mathrm{dm}}$ or $Q^{\mathrm{nw}}$ does not vary qualitatively thus conserving relative distances and gaps of eigenvalues, (iii) edge density for the graph- or network-based approaches Donner2010 .

Aside from those technical constraints, there are physical reasons that account for the turbulent convection flow under consideration. For Rayleigh–Bénard convection that is considered here, we use the Nusselt number ${\rm Nu}$, a dimensionless number that quantifies the global turbulent heat transfer across the plane, and the expected nondimensional radius of a convection roll $r_{roll}$ to obtain bounds for both cutoff scales:

[TABLE]

We recall here that $\varepsilon=\sqrt{2\delta}$. The lower bound in (25) is determined by the mean thickness of the thermal boundary layer, $\delta_{T}=0.028$ for the present Nusselt number of ${\rm Nu}=17.7$. Recall that all length scales are expressed in units of height $H$. The viscous boundary layer thickness can be determined for stress-free boundary conditions by a method that has been suggested by Petschel et al. Petschel2016 . They determined a so-called dissipation layer thickness from the intersection points (close to the bottom and top walls) of the line- and time-averaged profile of the kinetic energy dissipation rate with its plane mean. The analysis for the present flow gives a viscous boundary layer thickness of $\delta_{v}\approx 0.039$ which is slightly larger than the thermal boundary layer thickness.

The dynamics further restricts the choice of both closeness parameters. In the short time interval, the trajectories cover only one circulation and very few Lagrangian particles will switch from the left roll to the right roll or vice versa (see Fig. 1). Thus, we have to take a small cut-off parameter in order to capture relevant dynamical structures. For the long time interval, the spectra of eigenvalues are typically not stable in case of a small closeness parameter. To satisfy these restrictions for both methods, we get bounds on $\varepsilon$ for the short and long time intervals of

[TABLE]

As we will see further below, our choices satisfy $\delta_{T}<\delta_{v}<\varepsilon$.

In the following we want to visualize the differences and similarities of the results of the three algorithms. In general, the coherent sets are expected to be larger for the network-based analysis compared to the time-averaged diffusion maps method. The reason for this behavior is that the network-based approach does not take into account how close trajectories pass by (as long as it is closer than the threshold $\varepsilon$) and how long they reside in the vicinity of each other, thus having a less pronounced and hence larger “dynamic neighborhood”. As the strength of the diffusion decreases exponentially with increasing distance between particles diffusion the transition from coherent structure to incoherent flow is more easily detected for diffusion maps. Consequently, the definition of coherence is more strict in the diffusion maps method compared to the network-based method and this will lead to smaller coherent sets.

V.3.1 Short time interval up to $\tau=2020t_{f}$

For the particular parameter choices $\delta=0.002$ and $\varepsilon=0.0632$, we now show the results for the short time span in Figures 3, 4 and 5. For visualization purposes we plot the eigenvalues of $2\varepsilon^{-1}Q^{\mathrm{nw}}$.

The second eigenvector gives a separation of the left and right circulation rolls. The fifth eigenvector (resp. third for the FEM approach) gives a separation of the gyre cores and the background. We omit the plots for the third and fourth eigenvectors of $Q^{\mathrm{nw}}$ and $Q^{\mathrm{dm}}$ since they correspond to sub-partitions of the left and the right rolls only. They basically split each roll into halves, i.e., bottom and top or left an right for the third or fourth eigenvector, respectively. This occurrence can be explained in different ways. These (third and fourth) eigenvalues can be considered as higher multiplicities of the second eigenvalue interacting with the numerics. Furthermore these sub-partitions are still valid coherent sets which can be explained as follows. Within the selected time span the main motion is the rotation. However, even though most particles complete one rotation and the rotational speed varies with radius, this difference is not large enough to effectively separate particles from the inner (core) regions and the outer (background) regions. Therefore, the sub-partitions of the circulation rolls are more coherent than the gyre cores. As we will see later for the longer time span these sub-partitions do not occur, implying that they are then less coherent than the gyre cores.

V.3.2 Long time interval up to $\tau=2200t_{f}$

In the following we visualize the results for the particular parameter choices of $\delta=0.005$ and $\varepsilon=0.1$ in Figures 6, 7 and 8. Again, for visualization purposes we plot the eigenvalues of $2\varepsilon^{-1}Q^{\mathrm{nw}}$. We still get the separation of the left and the right side from the second eigenvector and the separation of the gyre cores and the background from the third eigenvector. We also note that the spectral gap is more prominent in the eigenvalue spectrum of $Q^{\mathrm{dm}}$ compared to $Q^{\mathrm{nw}}$. This could be advantageous in case of an unknown number of coherent sets.

V.3.3 Comparison of the three methods

Figure 9 compares the results by highlighting the trajectories that have been assigned differently by the network based and diffusion map method. This is done with a symmetrical difference of coherent sets $A^{\rm nw}$ and $A^{\rm dm}$, given by $(A^{\rm nw}\setminus A^{\rm dm})\cup(A^{\rm dm}\setminus A^{\rm nw})$. For simplicity we only plot the symmetric differences for the $Q^{\mathrm{nw}}$ and $Q^{\mathrm{dm}}$ results.

Furthermore, we observe a correlation between the distances between consecutive eigenvalues (see Fig. 6).

In the short time setting the coherent sets resulting from the FEM method are comparable in size and shape to the ones detected by the other two methods. For the long time interval, however, they are significantly smaller. This might be explained by the original construction of the dynamic Laplacian via dynamic isoperimetry in Ref. Froyland2015 : the method detects sets which keep a small boundary-to-volume ratio over the entire time interval of the flow evolution. It is thus a stricter criterion than the ones in the other two approaches.

Figure 10 displays the magnitude of the time-averaged velocity field. Small values can be identified in the elliptic centers of the counter-rotating rolls, the hyperbolic regions in the center of the top and bottom plates where Lagrangian particle pairs separate, and the four corners. The elliptic centers are the regions where the particles stay closely together for the longest time.

V.3.4 Lagrangian particle advection in time-averaged flow

In the Eulerian frame of reference, temporal averaging has to be performed in order to reveal the coherent large-scale patterns in the flow clearly Pandey:NatCommun2018 . In the following, we adapt these ideas to the present Lagrangian analysis. We therefore carry out a time-averaging in the following sense

[TABLE]

where the number of time steps $\mathcal{T}_{W}$ depends on the chosen window size $\Delta t$ in free-fall time units. The time-averaged trajectories $\{\hat{{\bm{x}}}_{i}(t_{\hat{k}})\}_{\hat{k}}$ are assembled such that they represent the same time interval of length $200t_{f}$, i.e., $\hat{k}$ depends on $\Delta t$. Here, we apply the minimum distance spectral trajectory clustering method (IV.1) to the time-averaged trajectories. Equivalent results can be achieved by applying the time-averaged diffusion maps method. The degree of coherence of the independent sets, interpreted as the size of the spectral gap, is improved by the time-averaging. Figure 11 shows the eigenvalues for the original setting, i.e., with no time-averaging, and for time-averaged trajectories with window sizes $\Delta t=5t_{f}$ and $14t_{f}$ (which corresponds with ${\cal T}_{W}=50$ and 140 respectively). For $\Delta t=5t_{f}$ a prominent gap is visible between the third and fourth eigenvalues. This corresponds to three almost decoupled sets, i.e., dynamically independent flow regions, which are the two cores of the convection rolls and the background. With increasing window size transitions from one side of the domain to the other side, which occur rarely for individual trajectories in the original setting, are removed and the Lagrangian time-averaged trajectories $\{\hat{{\bm{x}}}_{i}(t_{\hat{k}})\}_{\hat{k}}$ mostly remain in the initial flow region. Thereby, a segmentation into inner and outer cores is formed on both sides of the domain, corresponding to four almost decoupled sets. The spectral gap is shifted to appear between the fourth and fifth eigenvalues. This segmentation of the time-averaged trajectories is visualized in Figure 12. The transition from three to four almost decoupled sets occurs around a window size of $\Delta t=10t_{f}$. We can thus conclude that an additional averaging enhances the long-term coherence of the sets significantly even for the present two-dimensional case for which the small-scale dispersion by turbulence remains moderate.

VI Heat coherence in convection flow

VI.1 Temperature from a passive scalar perspective

The dynamical structure of heat transport can be analyzed from different perspectives. Although temperature is not a passive scalar, but a prognostic variable, once the evolution of the full system is known (i.e., the velocity $\bm{u}$ and total temperature $T$ fields are computed), we can view heat as a passive scalar with evolution governed by the advection-diffusion equation (see also eq. (2))

[TABLE]

with boundary conditions that we have defined in Sec. II. Here, ${\bm{\Psi}}$ denotes the total dimensionless heat flux vector and the average $\langle\Psi_{z}\rangle_{x,t}={\rm Nu}/\sqrt{\rm RaPr}$. The profile $T_{0}=1-z$ is the dimensionless linear equilibrium temperature profile (see also eq. (4)). Furthermore, $\Gamma_{D}$ and $\Gamma_{N}$ are parts of the domain boundary, where prescribed constant temperature (i.e., Dirichlet conditon on $\Gamma_{D}$) and insulating wall (i.e., Neumann condition on $\Gamma_{N}$) boundary conditions are applied. To recall, an insulated side wall implies that the normal derivative vanishes, $\smash{\frac{\partial T}{\partial n}}=0$.

This suggests that coherence with respect to heat could be analyzed by the approach from Sec. III. The situation is more delicate, however, as the identification of “heat packages” with particles has to be done properly. One issue is, for instance, that heat is not a conserved quantity and can enter and exit at the top and bottom boundaries, $\Gamma_{D}$. Thus the boundary conditions in (28) have to be incorporated in the analysis. We start by briefly discussing a popular line of approaches that is however inappropriate for heat-coherence investigations.

In the so-called “heatline” approaches Kimura1983 ; Costa2006 ; Speetjens2012 ; Pratt2016 ; Balasuriya2018 , the temperature field is considered as a concentration field driven by some velocity field ${\bm{u^{\mathrm{H}}}}$. In order to get the correct flux, ${\bm{\Psi}}={\bm{u^{\mathrm{H}}}}T$ has to hold, giving

[TABLE]

Note that the system (28) is translation-invariant in the following sense: For any constant $T_{c}\in\mathbb{R}$ and given solution $T$ of (28), the translated field $T^{\prime}=T+T_{c}$ solves (28) if the initial and Dirichlet boundary conditions are translated accordingly, i.e., $T^{\prime}_{\Gamma}=T_{\Gamma}+T_{c}$ and $T^{\prime}_{0}=T_{0}+T_{c}$. This means intuitively, that the evolution does not change even if we express temperature with respect to an arbitrary reference value.

In contrast, the “heat velocity” field ${\bm{u^{\mathrm{H}}}}$ changes under translation of the temperature field, giving completely different “heat trajectories”. Nevertheless, by construction, the global temperature field $T$ is advected correctly by ${\bm{u^{\mathrm{H}}}}$. This means, that the usefulness of the “heat velocity field” is restricted to considerations involving the global temperature distribution, but internal fluctuations of heat and the internal transport of heat are biased by the choice of reference temperature.

It should be remarked at this point that in Speetjens2012a ; Speetjens2012 this problem is alleviated by splitting convective and conductive contributions of heat transport, where the translational invariance is taken up completely by the conductive part. To the remaining (convective) flux a velocity field ${\bm{u}}^{\mathrm{conv}}={\bm{\Phi}}^{\mathrm{conv}}/T^{\mathrm{conv}}$ can be assigned, which then describes the convective heat transfer in a Lagrangian manner. However, as we would like to describe structures governing the entire heat transport (and not just its convective part), we take a different route in the following.

VI.2 Randomly evolving heat packages and induced transport

The microscopic evolution of heat can be described by advection and stochastic fluctuations. Let us consider the stochastic differential equation

[TABLE]

where ${\bm{x}}(t)$ denotes the random position of a particle driven by the drift ${\bm{u}}$ and by white noise ${\bm{\eta}}$ (i.e., $\mathrm{Cov}[{\bm{\eta}}(t),{\bm{\eta}}(s)]=\delta(t-s)$, where $\delta$ denotes the Dirac delta distribution). The probability distribution $p({\bm{x}},t)$ of {{\color[rgb]{0,0,1}\bm{x}}}(t) satisfies the Fokker–Planck equation

[TABLE]

If we replace $p$ by $T$, this is identical with the evolution equation in (28). Thus, scaling $T$ to have integral one (setting the total heat to be unity), we can express the evolution of heat as the evolution of the probability distribution of an ensemble of particles, evolving mutually independently with their motions governed by (30). This ensemble has initial distribution $T_{0}$. Single particles can thus be viewed as “heat packages”, each carrying one unit of heat.

We will illustrate some adaptations necessary according to the boundary conditions in the Rayleigh–Bénard convection problem. Neumann zero boundary conditions naturally translate into reflecting particles at the side walls $\Gamma_{N}$. Furthermore, we have $T=0$ at the top and $T=1$ for the temperature in dimensionless units. These conditions simply translate into absorbing every trajectory hitting the top lid, and re-injecting a new one at a random position at the bottom lid. Trajectories are reflected once they try to exit at the bottom lid. Note that the dynamical equation (30) is independent of the choice of a reference temperature, which is well aligned with the intuition that the dynamics of single heat packages should not depend on the reference temperature. The total heat transport depends on the reference temperature through the initial distribution $T_{0}$. The absorption and re-injection of trajectories at the boundaries results in some trajectories with a shorter life span than the considered total integration time since they leave before the end or are seeded later.

Combining (Banisch2017, , Theorem 3) with (Froyland2014, , Sec. 4.1) shows that eigenvectors of the newly obtained matrix $P^{\mathrm{dm}}$ are relaxed solutions of the problem (cf. (12))

[TABLE]

where $\mathbb{A}\subset\{1,\ldots,N\}$ is a subset of all trajectories (including trajectories seeded later), $\mathbb{A}^{C}$ is its complement, $T(\mathbb{A})$ is the total heat content of the set $\mathbb{A}$ (as every particle is a heat package with one unit of heat, this is the cardinality of $\mathbb{A}$), and $T(\mathbb{A}\to\mathbb{B})$ is the heat moving from set $\mathbb{A}$ to set $\mathbb{B}$ under the heat dynamics. Thus, $T(\mathbb{A}\to\mathbb{A})$ describes the heat remaining in $\mathbb{A}$. Solutions of (32) attempt to partition the domain into two subdomains between which there is as little heat exchange as possible, while the normalization by $T(\mathbb{A})$ and $T(\mathbb{A}^{C})$ avoids highly unbalanced partitions where one of the sets contains almost all of the heat.

While the problem (32) is clearly not invariant under translation with respect to a reference temperature, in certain cases it can be argued that it would give highly consistent solutions for any reference value. Let us assume that both $T(\mathbb{A})\approx T(\mathbb{A}^{C})$, and that the trajectories in $\mathbb{A}$ and $\mathbb{A}^{C}$ cover approximately the same area in physical space. Then problem (32) and the problem $\min_{\mathbb{A}}\{T(\mathbb{A}\to\mathbb{A}^{C})+T(\mathbb{A}^{C}\to\mathbb{A})\}$ have approximately the same solutions, while the latter describes total heat exchange between the two subdomains. Now, the physically relevant quantity is the net heat flow, which is the difference of $T(\mathbb{A}\to\mathbb{A}^{C})$ and $T(\mathbb{A}^{C}\to\mathbb{A})$, and this is almost invariant under changing the reference measure, if the spatial domains occupied by $\mathbb{A}$ and $\mathbb{A}^{C}$ are almost of the same size, because the contributions due to translation by a reference value cancels out.

VI.3 Heat coherence in two-dimensional Rayleigh–Bénard convection

In the following, this methodology is applied to the present two-dimensional Rayleigh–Bénard system. We conduct the analysis for the diffusion maps method (IV.2) only since the other two methods gave similar results as we saw in Sec. V. We first look again at a small time interval which means here $\mathcal{T}=100$ steps which corresponds to a time interval of $10t_{f}$.

If we look at the initial and the final time slices (see Figs. 13 and 14, respectively) we can see a separation of the right and left sides of the domain. Furthermore, the color of the particles that leave the domain fast and particles that enter the domain late is green, which implies that they have a zero value in the eigenvector. These Lagrangian particles correspond to thermal plumes as seen in the temperature field snapshot in Figure 15.

The algorithm identifies them as less relevant for the overall dynamics due to their short life span, which in turn implies that they are highly relevant for the effective heat transport from the bottom to the top.

We now look at a longer sequence of $\mathcal{T}=2000$ steps which corresponds to a time span of $200t_{f}$. The second eigenvector (see Figure 16, top panel) reveals now that the cores of the convection rolls are the most coherent features for the temperature evolution. This implies that there is almost no effective heat transport through the cores aside from diffusion which will become increasingly subdominant as the Rayleigh number grows.

It is important to note that even though the most heat-coherent sets do not contribute to vertical heat transport in this example, other cases may arise as well. If there would be big bubbles of hot fluid moving (slowly, compared to the considered time span) from the bottom to the top, the algorithm would identify them as heat coherent.

VII Conclusion

The main objectives of this work were to discuss coherence in a simple two-dimensional turbulent convection flow from a Lagrangian point of view and to relate it to the more frequently used perspective of the Eulerian frame of reference. We therefore compared three different Lagrangian approaches, (i) minimum distance spectral-clustering based, (ii) diffusion-map based, and (iii) dynamic Laplacian based analysis of Lagrangian data. We find that all three methods identify similar coherent sets, the core regions of the two convection rolls in our example flow. These are the areas in which neighboring Lagrangian particles would remain close to each other for the longest time before dispersion by small-scale turbulence would tear them apart from each other. Our analysis shows that these regions are the complementary to those which would be highlighted in a time-averaging procedure in the Eulerian frame of reference, namely the ridges of hot upwelling and cold downwelling fluid between these two circulation rolls. However, the notion of coherence is less strict in the minimum distance spectral-clustering based analysis compared to the diffusion-map based analysis and it varies with interval length for the dynamic Laplacian based analysis.

Furthermore, we introduced the concept of time averaging in the Lagrangian frame and demonstrated that – similar to the Eulerian case, coherence of structures is improved.

Finally, we discussed the concept of heat coherence in the present setting. We therefore suggested an approach to analyze the transport of non-passive scalar quantities including boundary sources and sinks utilizing the theory of coherent sets. We find that effective heat transport only occurs outside the cores.

How do the present Lagrangian and the (standard) Eulerian description compare to each other? In RBC flows, the prominent structures are flow circulations between the top and bottom plates in connection with rising hot or falling cold thermal plumes. The Eulerian picture in Figure 1 reveals a space-filling coherent pattern with a characteristic length $\lambda$ (which is here simply the horizontal extension of a pair of counter-rotating rolls) and a characteristic time scale $\tau$ which is proportional to the typical turnover time inside a roll Pandey:NatCommun2018 . The Eulerian analysis highlights the regions of the flow that contribute most to the convective heat transfer, namely, the hot upwellings and the cold downwellings between the counter-rotating rolls. Lagrangian methods are connected to the material transport by the evolving flow. As shown in this work, all Lagrangian methods detect the core regions of the large-scale circulation rolls as the coherent sets in which fluid particles remain together for the longest time. Thus, they reveal regions that form the spatial complement to the Eulerian ones, those that contribute least to the turbulent heat transfer. One major advantage of the discussed Lagrangian analysis methods is their objectivity, i.e., coordinate-frame independence (see, e.g., Haller2015 ) that helps to identify barriers to the turbulent transport in the flow. Lagrangian methods provide complementary powerful tools to reduce the dynamics to that of a few relevant degrees of freedom. They may not only be applicable in numerical simulations, but also for the growing number of experimental techniques that provide Lagrangian particle tracks, e.g. Monaro2019 ; Cierpka2013 ; Schanz2016 .

The present example was a two-dimensional flow at a moderate Rayleigh number in a working fluid with a large Prandtl number. Such a setup is an appropriate starting point for Lagrangian studies as the temperature field obeys a small diffusivity and the magnitude of turbulent velocity fluctuations remains small. The Reynolds number which quantifies the turbulent momentum transfer is ${\rm Re}\approx 96$ in our case. As a part of the future work, we will extend the mathematical foundations of the present Lagrangian framework to temporally averaged turbulence fields and apply these techniques in three-dimensional settings for extended flows. This is necessary to compensate for the enhanced turbulent dispersion which is always probed by Lagrangian methods and which increases as Rayleigh numbers are increased or Prandtl numbers are decreased. These efforts are partly under way and will be reported elsewhere.

Acknowledgments

This work of C.S., M.S. and A.P. is supported by the Priority Programme DFG-SPP 1881 ”Turbulent Superstructures” of the Deutsche Forschungsgemeinschaft (DFG). We acknowledge computing resources by the Large Scale Project pr62se of the Gauss Centre for Supercomputing. The authors thank Christian Cierpka, Daniel Karrasch and Nathanael Schilling for helpful comments.

Appendix A Additional Material

A.1 Interpretation in terms of graph Laplacians.

The stiffness matrices $K^{k}$ are symmetric and have zero row and column sums FroylandJunge2018 . Its diagonal entries are positive and the off-diagonal entries non-positive in certain important cases (and also in all numerical experiments). Because of this, for the graph $G^{k}$ with nodes $X^{k}$ and edges defined by the edges of the triangulation, we can write $K^{k}=\Pi^{k}-W^{k}$, where

[TABLE]

are the weights assigned to edges $(i,j)$, and $\Pi^{k}_{ii}=K^{k}_{ii}$ is a diagonal matrix with the $i^{\rm th}$ diagonal element equal to the degree of node $i$ in $G^{k}$, namely $\Pi^{k}_{ii}=-\sum_{j}K^{k}_{ij}$. Thus, we can view $K^{k}$ as an (unnormalized) graph Laplacian.

Note that the entry $K^{k}_{ij}$ will decrease in magnitude as the distance between the associated data points ${\bm{x}}_{i}(t_{k}),{\bm{x}}_{j}(t_{k})$ increases until the Delaunay edge $(i,j)$ no longer exists in $G^{k}$.

When we solve the eigenproblem $\bar{K}v=\lambda Mv$, we normalize by the symmetric, nonnegative mass matrix $M$, i.e., based on local area or volume elements that neighboring data points enclose. This normalization is different to the standard graph Laplacian normalization (which is based on node degree only): $M$ is not diagonal and the small number of off-diagonal entries of $M$ coincide with the off-diagonal entries of $K^{0}$ and correspond to arcs $(i,j)$ that are in the graph $G^{0}$. In fact, normalizing by the mass matrix automatically handles nonuniformly distributed data, because if initial points ${\bm{x}}^{0}_{i},{\bm{x}}^{0}_{j}$ are far apart the value of $M_{ij}$ will be commensurately larger.

Note further that since a triangulation is used here, there is no free parameter (like the cutoff radius) to choose and that the method can always yield a decomposition of the entire domain $\Omega$ (more precisely, the convex hull of the data points) into coherent sets.

A.2 Comparison of methods based on the rate-matrix interpretation

From the point of view of rate matrices in Sec. V.1, on the one hand, $Q^{\mathrm{nw}}$ defines a process where every state has holding time $1$ (because $Q^{\mathrm{nw}}_{ii}=1$), and the process has equal probabilities $1/D_{ii}$ to jump to any of its “neighbors”, where $i$ and $j$ are neighbors if $A_{ij}=1$. We observe that the proximity parameter $\varepsilon$ can also be interpreted as diffusion strength (diffusion coefficient), as the holding times of the Markov process are all one, and the “jumps” cover an $\varepsilon$-neighborhood in space. This can be seen from noting that a Brownian diffusion $c\,B_{t}$ of strength $c$ has standard deviation $\sqrt{\mathrm{Var}[c\,B_{t}]}=c\sqrt{t}$, i.e., if it has strength $c=\varepsilon$, then it produces $\varepsilon$ mean deviation in unit time. Thus, in a first order approximation, we could interpret $Q^{\mathrm{nw}}$ as a diffusion on trajectories with strength $\varepsilon$.

On the other hand, $Q^{\mathrm{dm}}$ defines a random walk where the $i$-th state has holding time

[TABLE]

which follows from the construction noting that $K_{t_{k},ii}=1$. We observe immediately that the scale parameter $\delta$ features directly as a timescale. As the “range” of the kernel 333The cutoff radius $r$ is defined by $\smash{k_{\delta}({\bm{x}},{\bm{y}})=h(\rho^{2}/\delta)<\vartheta}$ for all $\|{\bm{x}}-{\bm{y}}\|=\rho>r$, with some small threshold parameter $\vartheta$. This directly yields $r=\mathrm{const}\cdot\sqrt{\delta}$. and thus the mean jump distance is $\mathcal{O}(\sqrt{\delta})$, a similar consideration as above shows that $Q^{\mathrm{dm}}$ can be interpreted as a diffusion of strength $1$. Returning to (33), we see that the holding time grows very large if the trajectory $i$ has only few and distant neighbors (i.e., it is unlikely to jump over to trajectories that are not alike), and approaches $\delta$ from above if the trajectory has many close neighbors. The jump probabilities are readily encoded in the entries of $P^{\mathrm{dm}}$, and due to the construction involving the similarity matrix, it is more likely that the process jumps to a neighbor which is closer on average (in time).

To summarize this theoretical comparison of the methods from Secs. IV.1 and IV.2, both use spectral clustering with matrices that are interpretable as rate matrices of certain Markov jump processes on the trajectories. By this, those trajectories belong to the same coherent set which are likely to be reached from one another; opposed to unlikely transitions from other trajectories. This behavior is often referred to as metastability or almost invariance, and its connection to the dominant eigenmodes of the process’ jump matrix is well analyzed Davies1982a ; Davies1982b ; Bovier2002 ; Froyland2005 ; Huisinga2006 ; Schuette2013 . The difference in the methods considered here relies in characteristics of the associated Markov processes: The network-based process jumps to its neighbors with equal probabilities (thus utilizing only a binary information about distances; whether they are smaller or larger than $\varepsilon$); while the diffusion-maps based process prefers to jump to closer neighbors (thus utilizing a more refined distance information). While the proximity parameter $\varepsilon$ can be viewed as strength of the Markov process (diffusion) generated on the trajectories by $Q^{\mathrm{nw}}$, the strength of the analogous diffusion generated by $Q^{\mathrm{dm}}$ is always one. As these diffusion processes are discrete in space (because of jumping between a finite set of trajectories), their mean jump distance governs the finest scales they can resolve: $\mathcal{O}(\varepsilon)$ and $\mathcal{O}(\sqrt{\delta})$. Coherent sets below these scales can not be detected.