# Dissipation-based proper orthogonal decomposition of turbulent Rayleigh-
Bénard convection flow

P.J. Olesen<EMAIL_ADDRESS>Department of Civil and Mechanical Engineering,
Technical University of Denmark, Koppels Allé 403, 2800 Kgs. Lyngby, Denmark
L. Soucasse  Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C,
91190 Gif-sur-Yvette, France.  Netherlands eScience Center, Science Park 402
(Matrix III), 1098 XH Amsterdam, Netherlands  B. Podvin  Université Paris-
Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, 91190 Gif-sur-Yvette, France.
C.M. Velte  Department of Civil and Mechanical Engineering, Technical
University of Denmark, Koppels Allé 403, 2800 Kgs. Lyngby, Denmark

###### Abstract

We present a formulation of proper orthogonal decomposition (POD) producing a
velocity-temperature basis optimized with respect to an $H^{1}$ dissipation
norm. This decomposition is applied, along with a conventional POD optimized
with respect to an $L^{2}$ energy norm, to a data set generated from a direct
numerical simulation of Rayleigh-Bénard convection in a cubic cell
($\mathrm{Ra}=10^{7}$, $\mathrm{Pr}=0.707$). The data set is enriched using
symmetries of the cell, and we formally link symmetrization to degeneracies
and to the separation of the POD bases into subspaces with distinct
symmetries. We compare the two decompositions, demonstrating that each of the
20 lowest dissipation modes is analogous to one of the 20 lowest energy modes.
Reordering of modes between the decompositions is limited, although a corner
mode known to be crucial for reorientations of the large-scale circulation is
promoted in the dissipation decomposition, indicating suitability of the
dissipation decomposition for capturing dynamically important structures.
Dissipation modes are shown to exhibit enhanced activity in boundary layers.
Reconstructing kinetic and thermal energy, viscous and thermal dissipation,
and convective heat flux, we show that the dissipation decomposition improves
overall convergence of each quantity in the boundary layer. Asymptotic
convergence rates are nearly constant among the quantities reconstructed
globally using the dissipation decomposition, indicating that a range of
dynamically relevant scales are efficiently captured. We discuss the
implications of the findings for using the dissipation decomposition in
modeling, and argue that the $H^{1}$ norm allows for a better modal
representation of the flow dynamics.

††preprint: AIP/123-QED

##  I Introduction

Flows driven by thermal buoyancy are both of fundamental interest and of
relevance to a wide range of applications within e.g. engineering, geophysics,
and astrophysics. Such flows are often studied through simplified
configurations such as the Rayleigh-Bénard cell, in which a fluid layer
between two horizontal plates is heated from the bottom and cooled from the
top. For a given geometry the relevant control parameters are the Rayleigh
number, $\mathrm{Ra}$, and the Prandtl number, $\mathrm{Pr}$,

$\displaystyle\mathrm{Ra}$ $\displaystyle=\frac{\beta g\Delta
TH^{3}}{\nu\kappa}\,,$ $\displaystyle\mathrm{Pr}$
$\displaystyle=\frac{\nu}{\kappa}\,,$ (1)

where $\beta$ is the thermal expansion coefficient, $g$ the acceleration due
to gravity, $\Delta T$ the temperature difference between the plates, $H$ the
fluid layer height, $\nu$ the kinematic viscosity, and $\kappa$ the thermal
diffusivity. The Rayleigh number measures the ratio of buoyancy and
dissipative effects, while the Prandtl number depends only on the fluid,
measuring the ratio of momentum diffusivity to thermal diffusivity. The
resulting behavior of the flow is described through the Nusselt number,
$\mathrm{Nu}$, and the Reynolds number, $\mathrm{Re}$:

$\displaystyle\mathrm{Nu}$ $\displaystyle=\frac{\left\langle
wT\right\rangle_{A,t}-\kappa\left\langle\partial_{z}T\right\rangle_{A,t}}{\kappa\Delta
TH^{-1}}\,,$ $\displaystyle\mathrm{Re}$ $\displaystyle=\frac{UH}{\nu}\,,$ (2)

where $\left\langle\cdot\right\rangle_{A,t}$ denotes the average over a
horizontal plane and time, $w$ is the vertical velocity component, $T$ the
temperature, $\partial_{z}$ the vertical derivaive, and $U$ a characteristic
mean velocity.

The Nusselt number is the non-dimensionalized heat flux, including both
convective and conductive contributions. A central question is to determine
the dependence of the heat transfer, characterized by a Nusselt number, on the
control parameters $\mathrm{Ra}$ and $\mathrm{Pr}$. Shraiman and Siggia (1990)
have established strong relationships between these parameters and the thermal
and viscous dissipation, valid in the case of adiabatic side walls, showing in
particular that $\mathrm{Nu}=\left\langle\epsilon_{\theta}\right\rangle$,
where $\left\langle\epsilon_{\theta}\right\rangle$ is the non-dimensional
volume-averaged total thermal dissipation, and that
$\left(\mathrm{Nu}-1\right)\mathrm{Ra}=\left\langle\epsilon_{u}\right\rangle$,
where $\left\langle\epsilon_{u}\right\rangle$ is the non-dimensional volume-
averaged total viscous dissipation.

Grossmann and Lohse (2000) have developed a successful theory for predicting
$\mathrm{Nu}$ from $\mathrm{Ra}$ and $\mathrm{Pr}$. It is based on a scaling
analysis of the relative contributions of the bulk and the boundary layer to
the global kinetic and thermal dissipation, taking into account the relative
thickness of the thermal and viscous boundary layers. Crucial dynamical
information is therefore contained in these boundary layers, as pointed out by
Siggia (1994), Chillà and Schumacher (2012), and more recently by Scheel and
Schumacher (2014), who identified strong inhomogeneities in the boundary
layers using local analysis. At the same time, boundary layers present a
challenge in the modeling of buoyancy-driven flows, as pointed out by Hanjalić
(2002), due to rapidly varying properties and the lack of universal scalings
for these regions. The theory was further refined by Grossmann and Lohse
(2004) to decompose the thermal bulk contribution into a background part and a
plume-dominated one. Plumes are coherent structures, akin to detached thermal
boundary layers, that carry a large part of the heat flux (Shang  _et al._ ,
2003) and that are associated with large thermal dissipation events (Emran and
Schumacher, 2012). Plume-based analysis was recently used by Vishnu, De, and
Mishra (2022) to investigate the statistics of viscous and thermal dissipation
rates in a cubic cell. Shishkina and Wagner (2008) found that plumes collect
into a large-scale circulation (LSC) at a sufficiently high Rayleigh number.
The connection between small-scale plumes and the onset of the LSC has also
been evidenced in experiments(Xi, Lam, and Xia, 2004). The presence of the LSC
is responsible for the failure of standard one-point turbulence closures
Hanjalić and Launder (2022).

Coherent pattern extraction from velocity and thermal fluctuations can be
carried out using proper orthogonal decomposition (POD)(Berkooz, Holmes, and
Lumley, 1993; Berkooz, 1996), a statistical technique introduced to turbulence
studies by Lumley (1967). The first such study of a convective flow was
carried out by Sirovich and Park (1990). Several configurations have been
studied since: Bailon-Cuba, Emran, and Schumacher (2010) applied POD to a
cylindrical configuration; Podvin and Sergent (2012); Verdoold, Tummers, and
Hanjalić (2009) to rectangular cavities; Podvin and Sergent (2015) to two-
dimensional cells; and Soucasse  _et al._ (2019) to cubic cells. In all these
studies, the decomposition extracts the velocity and temperature fluctuations
with the largest energy. However, as described above, evidence suggests that
the viscous and the thermal dissipation rates are the fundamental quantities
to characterize the flow. This relates to a more general issue affecting POD-
based models, namely that convergence in energy does not itself guarantee that
the model reproduces temporal dynamics accurately. In particular, the focus on
large-scale structures resulting from energy optimization contrasts with the
fundamentally multi-scale nature of turbulence dynamics(Bergmann, Bruneau, and
Iollo, 2009). Convergence of the gradients ($H^{1}$ norm) is a stronger
requirement than $L^{2}$ convergence, and may be necessary for ensuring that
the dynamics of the reduced-order system reproduce those of the full flow.
Aubry, Lian, and Titi (1993) showed that a truncation containing more than 99%
of the energy was not enough to capture the dynamics of the Kuramoto-
Sivashinsky equation. It is therefore of interest to investigate the ability
of POD-based methods to capture structures optimized with respect to $H^{1}$
norms rather than the conventional $L^{2}$ norms.

The formalism underlying POD can be modified to produce decompositions with
respect to any quantity that can be expressed formally as a norm on some
tensor-, vector-, or scalar-valued field derived from the decomposed data.
Using POD this way produces a basis spanning the space corresponding to the
field in question optimally with respect to the chosen norm. Examples of such
norms are enstrophy and dissipation, which are $H^{1}$ norms computed from the
velocity gradient field and strain rate tensor field, respectively. Sengupta
and Dey (2004) studied the bypass transition using an enstrophy-optimized
vorticity decomposition, and Lee and Dowell (2020) investigated a gradient
basis used in combination with a conventional energy-optimized velocity basis
as a method for stabilizing a reduced order model (ROM) of a two-dimensional
lid driven cavity flow. Another example of POD applied to non-conventional
data sets was presented by Schiødt  _et al._ (2022), who decomposed Lagrangian
velocity data of particles suspended in a turbulent flow.

The extended POD introduced by Borée (2003) makes it possible to select the
decomposed quantity and the optimization quantity independently, enabling the
computation of modes spanning any desired quantity optimally with respect to a
norm computed from a different quantity, subject to the constraints discussed
above. For example, it allows the computation of velocity modes optimized with
respect to temperature norm and vice versa, as done by Podvin and Sergent
(2015), or of velocity modes optimized with respect to the dissipation norm,
as demonstrated by Olesen  _et al._ (2023). This provides a method for educing
dynamically important structures that may be missed in a conventional POD
analysis. In the present work this approach is used for obtaining velocity-
temperature modes that are optimized with respect to either total energy or
total dissipation. In the light of the role of dissipation in Rayleigh-Bénard
convection discussed above, the dissipation optimization can be expected to
produce structures representing important aspects of the flow dynamics that
are not well captured in energy optimized decompositions.

In the present work we consider a cubic Rayleigh-Bénard cell, which has
previously been investigated with standard POD by Soucasse  _et al._ (2019).
We apply the dissipation optimized POD to the data set, and produce velocity-
temperature modes than can be compared directly to their energy-optimized
counterparts, in terms of both the per-mode large scale organization and the
convergence of energy, dissipation and heat flux globally and in the boundary
layers.

The remainder of this paper is laid out as follows. In Section II we present
the generic formalism for POD and the particular formulations for the energy
and dissipation optimized versions, as well as the data set upon which the
subsequent analysis is based. The data set is symmetrized based on the
geometrical symmetries of the cell, and we analyze the formal consequences
that this has on the decomposition in Section III. In Section IV we analyze
the results of the energy and dissipation decomposition, including spectra,
large-scale organization of POD modes, and resolved boundary layer structures.
The convergence of reconstructed quantities in considered in Section V. We
discuss the implications of the results for the applicability of the formalism
in a modeling context in Section VI, and summarize the conclusions in Section
VII.

##  II POD Formalism

POD was introduced to turbulence studies by Lumley (1967) as a method for
educing large-scale coherent structures in turbulent flows. It has found use
both as a valuable analytical tool for understanding mechanisms in flow
dynamics (see Bakewell Jr and Lumley (1967) and Herzog (1986) for early
examples), as a modeling component in the context of reduced-order models
(ROMs) (Aubry  _et al._ , 1988; Deane  _et al._ , 1991), and for flow
control(Ravindran, 2000; Ly and Tran, 2001). In general terms, POD produces an
orthogonal basis for a given data set such that the convergence of an
expansion of the data set in this basis is optimal in the mean with respect to
the norm on the data set. Formally, we let the data set
$\mathcal{Q}=\\{q_{m}\\}_{m=1}^{M}$ of size $M$ be a subset of a Hilbert space
$\mathcal{H}$ with the associated inner product
$\left\langle\cdot,\cdot\right\rangle$ and norm
$\left\lVert\cdot\right\rVert=\sqrt{\left\langle\cdot,\cdot\right\rangle}$.
The data set is equipped with an averaging operation,
$\left\langle\left\\{\cdot_{m}\right\\}_{m=1}^{M}\right\rangle$. Each mode in
the resulting orthogonal modal basis $\\{\varphi_{n}\\}_{n=1}^{N}$ maximizes
the mean projection of the data set,

$\displaystyle\varphi_{n}$
$\displaystyle=\operatorname*{arg\,max}_{\varphi\in\mathcal{H}_{n}}\frac{\left\langle\left\\{\left|\left\langle\varphi,q_{m}\right\rangle\right|\right\\}_{m=1}^{M}\right\rangle}{\left\lVert\varphi\right\rVert^{2}}\,,\quad
n=1,2,\ldots,N\,,$ (3)

where
$\mathcal{H}_{n}=\mathcal{H}\setminus\operatorname{span}{\left(\\{\varphi_{n^{\prime}}\\}_{n^{\prime}=1}^{n-1}\right)}$,
and $N$ is the number of modes needed to span the data set, i.e., its
effective dimensionality, which satisfies $N\leq M$.

The POD modes are eigenmodes of the POD operator $R$,

$\displaystyle R\varphi_{n}$ $\displaystyle=\lambda_{n}\varphi_{n}\,,\quad
R\varphi=\left\langle\left\\{\left\langle\varphi,q_{m}\right\rangle
q_{m}\right\\}_{m=1}^{M}\right\rangle\,.$ (4)

This operator is Hermitian and positive semi-definite by construction,
ensuring that the eigenmodes resulting from (4) form a complete orthogonal
basis for the data set and that the eigenvalues are real and non-negative. The
eigenpairs are conventionally ordered by decreasing eigenvalues, corresponding
to modes being ordered by decreasing mean amplitude. The expansion of elements
of the data set,

$\displaystyle q_{m}$
$\displaystyle=\sum_{n=1}^{N}a_{nm}\sigma_{n}\varphi_{n}\,,$ (5)

where $\sigma_{n}=\sqrt{\lambda_{n}}$ are the singular values, and expansion
coefficients are given by

$\displaystyle a_{nm}$
$\displaystyle=\frac{1}{\sigma_{n}}\left\langle\varphi_{n},q_{m}\right\rangle\,,$
(6)

is thus optimally robust against truncation, preserving as much of the mean
norm of the truncated expansion as is possible for any basis. The eigenvalue
$\lambda_{n}$ is the mean energy resolved by the mode $\varphi_{n}$, which can
be interpreted as a measure of the importance of the mode in the
reconstruction of the data set. It should be noted, however, that eigenvalues
measure the importance of modes only in the narrow sense defined by the inner
product; it is thus entirely possible for features of central importance for
the dynamics of the flow to be resolved only by modes with low eigenvalues.
One consequence of this is the difficulty of constraining _a priori_ the order
of POD-based ROMs needed to achieve the desired model accuracy.

The POD coefficients $a_{nm}$ in (5) are uncorrelated,

$\displaystyle\left\langle\left\\{a_{nm}a_{n^{\prime}m}\right\\}_{m=1}^{M}\right\rangle$
$\displaystyle=\delta_{nn^{\prime}}\,,$ (7)

where $\delta_{nn^{\prime}}$ denotes the Kronecker delta. Different
normalization conventions exist regarding the definition of POD coefficients
and the expansion in (5), with singular values often absorbed in coefficients.
In that case $\sigma_{n}$ disappears from (5) and (6), and the right-hand side
in (7) becomes $\lambda_{n}\delta_{nn^{\prime}}$.

The normalization convention chosen in the present work emphasizes the analogy
between the expansion in (5) and the singular value decomposition (SVD). The
scaled data matrix $Q\in\mathbb{R}^{N_{f}\times M}$ is formed with entries
given by $Q_{ij}=\frac{1}{\sqrt{M}}q_{j}(x_{i})$, where $N_{f}$ is the number
of degrees of freedom in the flow field. The SVD of $Q$ consists of the
decomposition

$\displaystyle Q$ $\displaystyle=\Phi\Sigma A^{\mathsf{T}}\,,$ (8)

where $\Phi\in\mathbb{R}^{N_{f}\times N}$ and $A\in\mathbb{R}^{M\times N}$ are
orthonormal matrices with POD modes and modal coefficients as their respective
columns, and $\Sigma\in\mathbb{R}^{N\times N}$ is a diagonal matrix with the
singular values along its diagonal. The eigenvalue problem in (4) can then be
compactly written as

$\displaystyle QQ^{\mathsf{T}}\Phi$ $\displaystyle=\Phi\Lambda\,,$ (9)

where $\Lambda=\Sigma\Sigma^{\mathsf{T}}\in\mathbb{R}^{N\times N}$ is a
diagonal matrix with eigenvalues along its diagonal.

Commonly, $N_{f}\gg M$, and the $N_{f}\times N_{f}$ matrix $QQ^{\mathsf{T}}$
becomes impractically large. The method of snapshots introduced by Sirovich
(1987) replaces the matrix eigenvalue problem in (9) with the snapshot
eigenvalue problem,

$\displaystyle Q^{\mathsf{T}}QA$ $\displaystyle=A\Lambda\,,$ (10)

where the matrix $Q^{\mathsf{T}}Q\in\mathbb{R}^{M\times M}$ presents a much
more manageable eigenvalue problem. After solving (10) for the coefficient
matrix $A$ and eigenvalue matrix $\Lambda$ the mode matrix can be recovered as
$\Phi=QA\Sigma^{-1}$.

In the original formulation by Lumley (1967) the data set was an ensemble of
independent flow realizations defined on a spatio-temporal domain. This method
is in principle capable of capturing the full spatial and temporal dynamics
governing the data set, although the size of the necessary data sets may make
it prohibitively expensive for fully three-dimensional time-dependent flows. A
number of more computationally tractable approaches have been suggested.
Commonly, the ensemble is formed from a set of measurement or simulation
snapshots, each representing a single instant in time. Only purely spatial
correlations can be captured from such a data set, as all information relating
to temporal structure is discarded in treating the time series as a
statistical ensemble. Zhang  _et al._ (2023) demonstrated a domain combining
spatial dimensions and phase for a two-dimensional periodic lid driven cavity
flow, forming the ensemble from different periods of the flow. The combined
spatial and temporal dynamics within a period could thus be described. Lumley
(1967) showed how dimensions exhibiting homogeneity and periodicity could be
eliminated from the domain by replacing $R$ with its Fourier transform along
these dimensions, yielding a separate eigenvalue problem for each wave number
or frequency. This approach has become known as spectral POD (Towne, Schmidt,
and Colonius, 2018) (not to be confused with the identically named method
proposed by Sieber, Paschereit, and Oberleithner (2016)). The method has also
been used to approximate POD modes along dimensions that were not manifestly
periodic, although this introduces certain deviations from the optimality
offered by true POD, cf. Hodžić, Olesen, and Velte (2022).

###  II.1 Energy POD

In this work we consider the decompositions based on two different ensembles:
one formed from a set of snapshots of velocity and temperature data, leading
to a decomposition optimized with respect to total energy (kinetic $+$
thermal); and one formed from snapshots of the corresponding strain rate
tensor and temperature gradient data, optimized with respect to total (viscous
$+$ thermal) dissipation. The velocity-temperature data set is described in
more detail in Section II.4. In either case the POD is fully specified by
defining the ensemble $\mathcal{Q}=\\{q_{m}\\}_{m=1}^{M}$ and the inner
product $\left\langle\cdot,\cdot\right\rangle$ used in (3), (4), and (6).

For the energy POD the ensemble
$\mathcal{Q}_{E}=\\{q_{m}\\}_{m=1}^{M}\subset\mathcal{H}_{E}$ consists of
combined velocity and temperature snapshots,

$\displaystyle q_{m}$
$\displaystyle=\left(u_{m},v_{m},w_{m},\theta_{m}\right)\,,$ (11)

which are elements in the Hilbert space

$\displaystyle\mathcal{H}_{E}$
$\displaystyle=\left\\{\varphi:\Omega\rightarrow\mathbb{R}^{4}\left|\left\langle\varphi,\varphi\right\rangle_{E}<\infty\right\\}\right.\,.$
(12)

This space is equipped with the inner product
$\left\langle\cdot,\cdot\right\rangle_{E}$ and norm
$\left\lVert\cdot\right\rVert_{E}$ corresponding to total (kinetic and
thermal) energy, defined as

$\displaystyle\left\langle\varphi,\psi\right\rangle_{E}$
$\displaystyle=\sum_{i=1}^{3}\int_{\Omega}\varphi^{i}\psi^{i}\,\mathrm{d}x+\gamma_{E}^{2}\int_{\Omega}\varphi^{\theta}\psi^{\theta}\,\mathrm{d}x\,,$
(13a) $\displaystyle\left\lVert\varphi\right\rVert_{E}$
$\displaystyle=\sqrt{\left\langle\varphi,\varphi\right\rangle_{E}}\,,$ (13b)

where $\varphi,\psi\in\mathcal{H}_{E}$. The first of the two terms in (13a)
accounts for kinetic energy, and the second term accounts for thermal energy.
The weighing factor $\gamma_{E}$ in (13a) is defined so as to ensure that
kinetic and thermal energy contribute equally to the norm when averaged over
the snapshot ensemble,

$\displaystyle\gamma_{E}$
$\displaystyle=\left(\frac{\left\langle\left\\{\sum_{i=1}^{3}\int_{\Omega}\left(q_{m}^{i}\right)^{2}\,\mathrm{d}x\right\\}_{m=1}^{M}\right\rangle}{\left\langle\left\\{\int_{\Omega}\left(\theta_{m^{\prime}}\right)^{2}\,\mathrm{d}x\right\\}_{m^{\prime}=1}^{M}\right\rangle}\right)^{\frac{1}{2}}\,.$
(14)

The total energy POD operator is then

$\displaystyle R_{E}\varphi$
$\displaystyle=\left\langle\left\\{\left\langle\varphi,q_{m}\right\rangle_{E}q_{m}\right\\}_{m=1}^{M}\right\rangle\,.$
(15)

Solutions to the eigenvalue problem resulting from setting $R=R_{E}$ in (4)
are temperature-velocity modes optimized with respect to the total energy.
This was done by Soucasse  _et al._ (2019) for the data set also considered in
this work, in which the results serve as a point of comparison for results of
the dissipation POD described in the following.

###  II.2 Dissipation POD

The dissipation POD is based on the idea of letting mean dissipation take on
the role played by mean energy in the energy POD formulated in the previous
section. The mean dissipation in question is formed by combining mean viscous
dissipation, $\left\langle\varepsilon_{u}\right\rangle$, and mean thermal
dissipation, $\left\langle\varepsilon_{\theta}\right\rangle$. Our task is
therefore to formulate the relevant ensemble, along with an inner product with
an associated norm that produces the total dissipation when applied to the
ensemble. For this purpose, we consider the strain rate tensor (SRT) and
thermal gradient,

$\displaystyle q_{m,\mathrm{srt}}^{ij}$
$\displaystyle=\frac{1}{2}\left(\nabla^{i}q_{m}^{j}+\nabla^{j}q_{m}^{i}\right)\,,$
(16a) $\displaystyle q_{m,\mathrm{tg}}^{i}$
$\displaystyle=\nabla^{i}\theta_{m}\,,$ (16b)

with $i,j\in\\{1,2,3\\}$. From the SRT and thermal gradient we obtain the mean
viscous and thermal dissipation, and, by combining these two, the total
dissipation,

$\displaystyle\left\langle\varepsilon_{u}\right\rangle$
$\displaystyle=\left\langle\left\\{\sum_{i,j=1}^{3}\int_{\Omega}\left|q_{m,\mathrm{srt}}^{ij}\right|^{2}\,\mathrm{d}x\right\\}_{m=1}^{M}\right\rangle\,,$
(17a) $\displaystyle\left\langle\varepsilon_{\theta}\right\rangle$
$\displaystyle=\left\langle\left\\{\sum_{i=1}^{3}\int_{\Omega}\left|q_{\mathrm{tg},m}^{i}\right|^{2}\,\mathrm{d}x\right\\}_{m=1}^{M}\right\rangle\,,$
(17b) $\displaystyle\left\langle\varepsilon_{\mathrm{tot}}\right\rangle$
$\displaystyle=\left\langle\varepsilon_{u}\right\rangle+\gamma_{D}^{2}\left\langle\varepsilon_{\theta}\right\rangle\,,$
(17c)

where, similarly to $\gamma_{E}$ in (14), $\gamma_{D}$ in (17c) ensures equal
mean contribution from viscous and thermal dissipation,

$\displaystyle\gamma_{D}$
$\displaystyle=\left(\frac{\left\langle\varepsilon_{u}\right\rangle}{\left\langle\varepsilon_{\theta}\right\rangle}\right)^{\frac{1}{2}}\,.$
(18)

The dissipation POD ensemble
$\mathcal{Q}_{D}=\\{q_{m}^{\prime}\\}_{m=1}^{M}\subset\mathcal{H}_{D}$ is
formed by joining the objects defined in (16),

$\displaystyle q_{m}^{\prime}$
$\displaystyle=\mathfrak{D}q_{m}=\left(q_{m,\mathrm{srt}},q_{m,\mathrm{tg}}\right)\,,$
(19)

to form elements in the Hilbert space
$\mathcal{H}_{D}=\mathcal{H}_{\mathrm{srt}}\times\mathcal{H}_{\mathrm{tg}}$.
Here, $\mathcal{H}_{\mathrm{srt}}$ and $\mathcal{H}_{\mathrm{tg}}$ are Hilbert
spaces containing $q_{m,\mathrm{srt}}$ and $q_{m,\mathrm{tg}}$, respectively,

$\displaystyle\mathcal{H}_{\mathrm{srt}}$
$\displaystyle=\left.\left\\{\varphi:\Omega\rightarrow\mathbb{R}^{3\times
3}\right|\left\langle\varphi,\varphi\right\rangle_{\mathrm{srt}}<\infty\right\\}\,,$
(20a) $\displaystyle\mathcal{H}_{\mathrm{tg}}$
$\displaystyle=\left.\left\\{\varphi:\Omega\rightarrow\mathbb{R}^{3}\right|\left\langle\varphi,\varphi\right\rangle_{\mathrm{tg}}<\infty\right\\}\,,$
(20b)

where the inner products $\left\langle\cdot,\cdot\right\rangle_{\mathrm{srt}}$
and $\left\langle\cdot,\cdot\right\rangle_{\mathrm{tg}}$ and their associated
norms are

$\displaystyle\left\langle\varphi_{\mathrm{srt}},\psi_{\mathrm{srt}}\right\rangle_{\mathrm{srt}}$
$\displaystyle=\sum_{i,j=1}^{3}\int_{\Omega}\varphi^{ij}_{\mathrm{srt}}\psi^{ij}_{\mathrm{srt}}\,\mathrm{d}x\,,$
(21a)
$\displaystyle\left\lVert\varphi_{\mathrm{srt}}\right\rVert_{\mathrm{srt}}$
$\displaystyle=\sqrt{\left\langle\varphi_{\mathrm{srt}},\varphi_{\mathrm{srt}}\right\rangle_{\mathrm{srt}}}\,;$
(21b)
$\displaystyle\left\langle\varphi_{\mathrm{tg}},\psi_{\mathrm{tg}}\right\rangle_{\mathrm{tg}}$
$\displaystyle=\sum_{i=1}^{3}\int_{\Omega}\varphi^{i}_{\mathrm{tg}}\psi^{j}_{\mathrm{tg}}\,\mathrm{d}x\,,$
(21c)
$\displaystyle\left\lVert\varphi_{\mathrm{tg}}\right\rVert_{\mathrm{tg}}$
$\displaystyle=\sqrt{\left\langle\varphi_{\mathrm{tg}},\varphi_{\mathrm{tg}}\right\rangle_{\mathrm{tg}}}\,,$
(21d)

where
$\varphi_{\mathrm{srt}},\psi_{\mathrm{srt}}\in\mathcal{H}_{\mathrm{srt}}$ and
$\varphi_{\mathrm{tg}},\psi_{\mathrm{tg}}\in\mathcal{H}_{\mathrm{tg}}$. The
operator $\mathfrak{D}:\mathcal{H}_{E}\rightarrow\mathcal{H}_{D}$ in (19) maps
velocity-temperature snapshots to the corresponding joint SRT and thermal
gradients.

The total dissipation inner product and norm can now be written as

$\displaystyle\left\langle\varphi^{\prime},\psi^{\prime}\right\rangle_{D}$
$\displaystyle=\left\langle\varphi_{\mathrm{srt}},\psi_{\mathrm{srt}}\right\rangle_{\mathrm{srt}}+\gamma_{D}^{2}\left\langle\varphi_{\mathrm{tg}},\psi_{\mathrm{tg}}\right\rangle_{\mathrm{tg}}\,,$
(22a) $\displaystyle\left\lVert\varphi^{\prime}\right\rVert_{D}$
$\displaystyle=\sqrt{\left\langle\varphi^{\prime},\varphi^{\prime}\right\rangle_{D}}\,,$
(22b)

with $\varphi^{\prime},\psi^{\prime}\in\mathcal{H}_{D}$. The total dissipation
POD operator is

$\displaystyle R_{D}\varphi^{\prime}$
$\displaystyle=\left\langle\left\\{\left\langle\varphi^{\prime},q_{m}^{\prime}\right\rangle_{D}q_{m}^{\prime}\right\\}_{m=1}^{M}\right\rangle\,.$
(23)

The decomposition of $\mathcal{Q}_{D}$ which is obtained using this operator
is optimized with respect to the norm $\left\lVert\cdot\right\rVert_{D}$
corresponding to total dissipation,

$\displaystyle\left\langle\varepsilon_{\mathrm{tot}}\right\rangle$
$\displaystyle=\left\langle\left\\{\left\lVert
q^{\prime}_{m}\right\rVert_{D}\right\\}_{m=1}^{M}\right\rangle\,.$ (24)

###  II.3 The extended snapshot method

The POD eigenvalue problems derived from the operators defined in (15) and
(23) are

$\displaystyle R_{E}\varphi_{E,n}$
$\displaystyle=\lambda_{E,n}\varphi_{E,n}\,,$ (25a) $\displaystyle
R_{D}\varphi_{D,n}^{\prime}$
$\displaystyle=\lambda_{D,n}\varphi_{D,n}^{\prime}\,.$ (25b)

The solutions of (25a) are energy eigenmodes which span $\mathcal{Q}_{E}$
optimally with respect to the inner product
$\left\langle\cdot,\cdot\right\rangle_{E}$, whereas solutions of (25b) are
dissipation eigenmodes spanning $\mathcal{Q}_{D}$ optimally with respect to
$\left\langle\cdot,\cdot\right\rangle_{D}$. The eigenmodes in the energy basis
$\\{\varphi_{E,n}\\}_{n=1}^{N}$ and those in the dissipation basis
$\\{\varphi_{D,n}^{\prime}\\}_{n=1}^{N}$ are elements of different Hilbert
spaces, and cannot be compared directly.

We generate velocity-temperature modes $\\{\varphi_{D,n}\\}_{n=1}^{N}$
corresponding to $\\{\varphi_{D,n}^{\prime}\\}_{n=1}^{N}$ using the extended
POD method introduced by Borée (2003), combined with the snapshot POD
formalism of Sirovich (1987). The resulting procedure is similar to what was
done by Podvin and Sergent (2015), as well as formally equivalent to the
procedure presented by Olesen  _et al._ (2023). In the snapshot formalism we
construct the snapshot matrices corresponding to $Q^{\mathsf{T}}Q$,

$\displaystyle\left[S_{E}\right]_{m,m^{\prime}}$
$\displaystyle=\frac{1}{M}\left\langle
q_{m},q_{m^{\prime}}\right\rangle_{E}\,,$ (26a)
$\displaystyle\left[S_{D}\right]_{m,m^{\prime}}$
$\displaystyle=\frac{1}{M}\left\langle
q_{m}^{\prime},q_{m^{\prime}}^{\prime}\right\rangle_{D}\,,$ (26b)

and solve the corresponding eigenvalue problems for the coefficient vectors
$\\{a_{p,n}\\}_{n=1}^{N}$,

$\displaystyle S_{p}a_{p,n}$ $\displaystyle=\lambda_{p,n}a_{p,n}\,,\quad
p\in\\{E,D\\}\,.$ (27)

Velocity-temperature modes are then constructed using the coefficients,

$\displaystyle\varphi_{p,n}$
$\displaystyle=\frac{1}{\sigma_{p,n}}\left\langle\left\\{a_{p,n,m}q_{m}\right\\}_{m=1}^{M}\right\rangle\,,\quad
p\in\\{E,D\\}\,,$ (28)

where $\sigma_{p,n}=\sqrt{\lambda_{p,n}}$.

The basis $\\{\varphi_{E,n}\\}_{n=1}^{N}$ formed this way is identical to the
eigenbasis obtained from (25a). In contrast, $\\{\varphi_{D,n}\\}_{n=1}^{N}$
is _not_ the eigenbasis obtained from (25b); the modes satisfy
$\varphi_{D,n}^{\prime}=\mathfrak{D}\varphi_{D,n}$, and form a complete (but
non-orthogonal) basis for the data set. The data set can therefore be
reconstructed using either basis as, cf. (5),

$\displaystyle q_{m}$
$\displaystyle=\sum_{n=1}^{N}a_{p,n,m}\sigma_{p,n}\varphi_{p,n}\,,\quad
p\in\\{E,D\\}\,.$ (29)

Both coefficient sets are uncorrelated, cf. (7),

$\displaystyle\left\langle\left\\{a_{p,n,m}a_{p,n^{\prime},m}\right\\}_{m=1}^{M}\right\rangle$
$\displaystyle=\delta_{nn^{\prime}}\,,\quad p\in\\{E,D\\}\,,$ (30)

allowing reconstruction of any second order mean quantity (such as energy,
dissipation, heat flux, and stresses) without the use of cross-modal terms.
This lets us decompose mean quantities into a sum of modal contributions,

$\displaystyle\Pi^{p}_{n}$
$\displaystyle=\frac{1}{M}\sum_{n^{\prime}=1}^{n}\lambda_{p,n^{\prime}}\Pi\left[\varphi_{p,n^{\prime}}\right]\,,$
(31)

where $\Pi\in\\{E_{u},E_{\theta},\varepsilon_{u},\varepsilon_{\theta},\Phi\\}$
identifies the reconstructed quantity (kinetic and thermal energy, viscous and
thermal dissipation, and convective heat flux), $p\in\\{E,D\\}$ designates the
decomposition, $n$ the order of the reconstruction (number of modes included),
and $\Pi[\varphi]$ denotes the reconstructed quantity computed from mode
$\varphi$. Mean quantities are recovered as
$\left\langle\Pi\right\rangle=\Pi^{p}_{N}$, which due to the completeness of
either basis is independent of $p$. This makes it possible to identify the
contribution from each mode to each mean profile, and, as we shall see, to
link this to the structures described by the mode.

###  II.4 Velocity and temperature fields

The numerical setup and associated data set are the same as those used in
Soucasse  _et al._ (2019). The configuration studied is a cubic Rayleigh-
Bénard cell filled with air, with isothermal horizontal walls and adiabatic
side walls. The air is assumed to be transparent and thermal radiation effects
are disregarded. Direct numerical simulations have been performed at a
Rayleigh number $\mathrm{Ra}=10^{7}$. The Prandtl number
$\mathrm{Pr}=\nu/\kappa$ is set to 0.707. All physical quantities are made
dimensionless using the cell size $H$, the reference time
$H^{2}/(\kappa\sqrt{\mathrm{Ra}})$ and the reduced temperature
$\theta=(T-T_{0})/\Delta T$, $T_{0}$ being the mean temperature between hot
and cold walls. Spatial coordinates are denoted $x$, $y$, $z$ ($z$ being the
vertical direction) and the origin is placed at a bottom corner of the cube.

The Navier–Stokes equations under the Boussinesq approximation are solved
using a Chebyshev collocation method (Xin and Le Quéré, 2002; Xin, Chergui,
and Le Quéré, 2008). Computations are made parallel using domain decomposition
along the vertical direction. Time integration is performed through a second-
order semi-implicit scheme. The velocity divergence-free condition is enforced
using a projection method. The spatial mesh is made of 81 Chebyshev
collocation points in each direction of space. We have checked that the number
of collocation points is sufficient to accurately discretize the boundary
layers according to the criterion proposed by Shishkina  _et al._ (2010). The
resulting Nusselt number averaged over time and over the horizontal plane is
$\mathrm{Nu}=16.24$, and the resulting Reynolds number based on the volume and
time average of the velocity magnitude is

$\displaystyle\mathrm{Re}$ $\displaystyle=\frac{\sqrt{\left\langle
u_{i}u_{i}\right\rangle_{v,t}\mathrm{Ra}}}{\mathrm{Pr}}=651\,.$ (32)

These values are in good agreement with reference literature results (Delort-
Laval  _et al._ , 2022).

A number of 1000 snapshots has been extracted from the simulation at a
sampling period of 10 in dimensionless time units. Isosurfaces of convective
heat flux and stream lines are shown in Figure 1 for the first snapshot
(Figure 1$a$), for the conditional mean of realizations with a particular
orientation of integrated angular momentum ($L_{x}<0$, $L_{y}>0$; Figure
1$b$), and for the full mean (Figure 1$c$) of the symmetry-enriched data set,
cf. Section III. The isosurfaces highlight cold ($\theta<0$) and warm
($\theta>0$) convective regions, where appreciable net upward convective heat
transfer ($w\theta>0.05$) occur due to warm fluid moving upwards in the case
of warm convective regions, and cold fluid moving downwards in the case of
cold convective regions. Figure 1$a,b$ shows a diagonally aligned roll
producing oppositely oriented convective regions on each side. This represents
the typical overall organization of the flow, with the orientation of the roll
shifting with irregular intervals. This contrasts with the mean flow shown in
Figure 1$c$, which does not represent a typical flow state.

Figure 1: The snapshot $m=1$ $(a)$, the conditional mean field given
$L_{x}<0$, $L_{y}>0$ $(b)$, and the full mean field $(c)$, shown as
isosurfaces of scaled convective heat flux corresponding to
$\Phi_{\mathrm{iso}}=0.05\cdot\max_{\Omega}\left(w\theta\right)$ and colored
by temperature $\theta$. In this and the similar plots in Figures 4 and 5 the
origin is located in the lower left corner of the cell, and the axes are
labeled as shown in $(a)$.

##  III Symmetries

As in Soucasse  _et al._ (2019) the data set is enriched by making use of the
statistical symmetries of the flow (Puigjaner  _et al._ , 2008). In the cubic
Rayleigh-Bénard cell, four quasi-stable states are available for the flow at
this Rayleigh number: the large-scale circulation (LSC) settles in one of the
two diagonal planes of the cube with clockwise or counterclockwise motion.
Reorientations from one state to another occur during the analyzed time
sequence, but each state is not necessarily equally visited. In order to
counteract this bias we construct an enlarged snapshot ensemble, obtained by
the action of the symmetry group of the problem on the original ensemble.

The symmetry group is formed from the four basic symmetry operations
$S_{k\in\\{x,y,z,d\\}}$(Puigjaner  _et al._ , 2008),

$\displaystyle S_{x}$
$\displaystyle:\left\\{\begin{array}[]{rcl}\left(x,y,z\right)&\mapsto&\left(1-x,y,z\right)\\\
\left(u,v,w,\theta\right)&\mapsto&\left(-u,v,w,\theta\right)\end{array}\right.\,,$
(33c) $\displaystyle S_{y}$
$\displaystyle:\left\\{\begin{array}[]{rcl}\left(x,y,z\right)&\mapsto&\left(x,1-y,z\right)\\\
\left(u,v,w,\theta\right)&\mapsto&\left(u,-v,w,\theta\right)\end{array}\right.\,,$
(33f) $\displaystyle S_{z}$
$\displaystyle:\left\\{\begin{array}[]{rcl}\left(x,y,z\right)&\mapsto&\left(x,y,1-z\right)\\\
\left(u,v,w,\theta\right)&\mapsto&\left(u,v,-w,-\theta\right)\end{array}\right.\,,$
(33i) $\displaystyle S_{d}$
$\displaystyle:\left\\{\begin{array}[]{rcl}\left(x,y,z\right)&\mapsto&\left(y,x,z\right)\\\
\left(u,v,w,\theta\right)&\mapsto&\left(v,u,w,\theta\right)\end{array}\right.\,.$
(33l)

The operations $S_{x}$, $S_{y}$, and $S_{z}$ describe mirroring in the planes
$x=0.5$, $y=0.5$, and $z=0.5$, respectively. This involves flipping the sign
of the corresponding velocity component, and, for $S_{z}$, also the sign of
the temperature fluctuation, due to the opposite thermal boundary conditions
imposed at the top and bottom by the cold and the hot plates. The operation
$S_{d}$ describes reflection in the diagonal plane $x=y$. The four symmetry
operations form a 16-element symmetry group $\mathcal{G}$, and the enriched
data set of $M=$16\,000$$ realizations is formed from the images of each of
the $\widetilde{M}=$1000$$ DNS realizations
$\\{\tilde{q}_{m}\\}_{m=1}^{\widetilde{M}}$ under each of the 16 symmetry
group elements,

$\displaystyle\left\\{q_{m}\right\\}_{m=1}^{M}$
$\displaystyle=\bigcup_{S\in\mathcal{G}}\left\\{S\tilde{q}_{m^{\prime}}\right\\}_{m^{\prime}=1}^{\widetilde{M}}\,.$
(34)

The symmetries in the data set lead to degeneracies in the resulting POD
bases. Here we show the formal origin and nature of these degeneracies using
the commutation relations involving the POD and symmetry operators. Since the
enriched data set is by construction invariant under each element in the
symmetry group, the POD operators $R_{p\in\\{E,D\\}}$ commute with the
symmetry operations. This implies that each POD operator shares an orthonormal
eigenbasis with each symmetry operation, such that for any symmetry $S_{k}$ a
POD eigenbasis $\\{\varphi_{p,n}\\}_{n=1}^{N}$ can be chosen satisfying
$S_{k}\varphi_{p,n}=\mu_{k,p,n}\varphi_{p,n}$. The symmetry operations are
unitary such that $\mu_{k,p,n}=\pm 1$, with the eigenvalue $+1$ indicating
symmetry and $-1$ indicating antisymmetry with respect to $S_{k}$. By
combining (33c) and (33l) we find

$\displaystyle S_{x}S_{d}$
$\displaystyle:\left\\{\begin{array}[]{rcl}\left(x,y,z\right)&\mapsto&\left(1-y,x,z\right)\\\
\left(u,v,w,\theta\right)&\mapsto&\left(-v,u,w,\theta\right)\end{array}\right.\,,$
(35c) $\displaystyle S_{d}S_{x}$
$\displaystyle:\left\\{\begin{array}[]{rcl}\left(x,y,z\right)&\mapsto&\left(y,1-x,z\right)\\\
\left(u,v,w,\theta\right)&\mapsto&\left(v,-u,w,\theta\right)\end{array}\right.\,,$
(35f)

implying that $\left[S_{x},S_{d}\right]\neq 0$, and by a similar argument,
$\left[S_{y},S_{d}\right]\neq 0$. From this it follows that $S_{x}$ and
$S_{y}$ do not share orthogonal eigenbases with $S_{d}$, which in turn implies
that the POD eigenbases are not unique, and that some form of degeneracy must
exist.

We can show that $S_{x}S_{d}=S_{d}S_{y}$ and $S_{y}S_{d}=S_{d}S_{x}$, from
which it follows that $[S_{x}+S_{y},S_{d}]=0$. Simultaneous eigenmodes of
$S_{x}$ and $S_{y}$ satisfying $\mu_{x,p,n}=\mu_{y,p,n}$ are therefore also
eigenmodes of $S_{d}$. Taken together, this means that degenerate POD
eigenspaces of dimension two exist and are spanned by modes with
$\mu_{x,p,n}=-\mu_{y,p,n}$. These degenerate pairs appear as pairs of repeated
eigenvalues in the POD spectra, as shown in Figure 2 in Section IV.

The pair of degenerate POD eigenmodes spanning a given eigenspace can be
chosen arbitrarily within that eigenspace, under the constraint of
orthonormality. Each degenerate pair can thus be formed from one mode
symmetric under $S_{x}$ ($\mu_{x,p,n}=+1$) and antisymmetric under $S_{y}$
($\mu_{y,p,n}=-1$), and one with the opposite configuration,
$\mu_{x,p,n^{\prime}}=-1$ and $\mu_{y,p,n^{\prime}}=+1$. These configurations
correspond to modes which are oriented along the $x$ and $y$ axis,
respectively, and which are not eigenmodes of $S_{d}$. The degeneracy is
responsible for resolving the different roll orientations, and modes in
degenerate eigenspaces are capable of supporting non-vanishing angular momenta
in the $xy$ plane, though it is not necessary that they do so. In this work we
use the eigenmodes of $S_{x}$ and $S_{y}$ to span degenerate eigenspaces,
resulting in modal angular momenta aligned with the $x$ and $y$ axes. One
could equally well have spanned the eigenspaces using pairs of $S_{d}$
eigenmodes with $\mu_{d,p,n}=+1$ and $\mu_{d,p,n^{\prime}}=-1$, aligning the
modal angular momenta with the diagonals. These modes would not be eigenmodes
of $S_{x}$ and $S_{y}$. While this choice would arguably have been in closer
agreement with the physical orientations attained by the flow (cf. Figure
1$a,b$, showing a typical diagonal-aligned state of the flow), the choice of
aligning degenerate pairs along the principal axes allows for direct
comparison with modes identified by Soucasse  _et al._ (2021). Since
$\mu_{d,p,n}$ is undefined whenever $\mu_{x,p,n}\neq\mu_{y,p,n}$ we are left
with 12 possible combinations of symmetry eigenvalues, or 12 distinct
isosymmetric families. Within the lowest 20 POD modes we find each of those
families represented by at least one mode, and in some cases by several modes.

The symmetry operations commute with the mapping $\mathfrak{D}$ from velocity-
temperature snapshots to SRT-thermal gradient snapshots in (19), and the
preceding discussion therefore applies to both the POD eigenmodes and to the
dissipation modes formed using (28).

##  IV POD spectra and modes

Figure 2 shows the POD spectra corresponding to the two decompositions. It
should be emphasized that while the energy eigenvalue $\lambda_{E,n}$ is the
mean _energy_ contribution of the energy mode $\varphi_{E,n}$, the dissipation
eigenvalue $\lambda_{D,n}$ is the mean _dissipation_ contribution of the
dissipation mode $\varphi_{D,n}$. Consequently, as the eigenvalues of the
different decompositions correspond to different norms, they are not directly
comparable. Both spectra are dominated by the lowest three modes, representing
circa 60% of the total energy for the energy decomposition, and 55% of the
total dissipation for the dissipation decomposition. The remaining parts of
the spectra decay slowly, with the energy POD spectrum exhibiting a slightly
faster rate of decay than the dissipation POD spectrum. With $n=100$ modes we
thus capture about 80% of the total energy and 65% of the total dissipation.
The asymptotic decay rates are determined by fitting the spectra to power
laws, $\lambda_{p,n}\sim n^{\alpha_{p}}$, in the interval $10\leq n\leq
10^{3}$. This yields the exponents $\alpha_{E}=-0.91$ and $\alpha_{D}=-0.79$
for the energy and dissipation spectra, respectively, confirming the more
rapid convergence of the energy POD spectrum.

Figure 2: POD spectra (solid lines, left axis) and cumulative spectra (dotted
lines, right axis) for PODs optimized with respect to total kinetic and
thermal energy ($\lambda_{E,n}$) and total viscous and thermal dissipation
($\lambda_{D,n}$). The spectra are normalized such that
$\sum_{n}\lambda_{E/D,n}=1$. Power laws fitted to the spectra in the interval
$10\leq n\leq 10^{3}$ are shown for the energy (black dashed line) and
dissipation spectra (black dotted line).

###  IV.1 Identification of modes

In order to facilitate a subsequent comparison between the modes produced by
the energy and dissipation decomposition we assign labels to the lowest 20
modes of each decomposition, adopting and extending the modal nomenclature
introduced for the first 12 energy optimized modes by Soucasse  _et al._
(2019, 2021). Each of the energy modes considered here is assigned a label
$G[\varphi_{E,n_{E}}]$, based on similarities to modes considered by Soucasse
_et al._ (2019, 2021), and the same label is then assigned to the dissipation
mode $\varphi_{D,n_{D}}$ selected among the lowest 20 modes that maximizes the
overlap defined using the energy inner product and norm in (13):

$\displaystyle\Gamma(n_{E},n_{D})$
$\displaystyle=\frac{\left|\left\langle\varphi_{E,n_{E}},\varphi_{D,n_{D}}\right\rangle_{E}\right|}{\left\lVert\varphi_{D,n_{D}}\right\rVert_{E}}\,.$
(36)

The overlap $\Gamma(n_{E},n_{D})\in[0;1]$ measures the degree of similarity
between modes of the different decompositions, with $\Gamma(n_{E},n_{D})=1$ in
the case when modes are identical up to normalization of $\varphi_{D,n_{D}}$,
while the overlap vanishes in case of orthogonal modes.

Table 1 summarizes all the modes investigated in this way in the present work,
their assigned label, their angular momentum and symmetry properties, and the
overlap of modes with $G[\varphi_{E,n_{E}}]=G[\varphi_{D,n_{D}}]$. The
assignment of shared labels between decompositions is found to be largely
consistent with shared large scale organization, as discussed in Section IV.3.

| I | II | III | IV | V
---|---|---|---|---|---
| $G[\varphi]$ | $n_{E}$ | $n_{D}$ | $L_{i}\neq 0$ | $\mu_{x}$ | $\mu_{y}$ | $\mu_{z}$ | $\mu_{d}$ | $\Gamma(n_{E},n_{D})$
| $M$ | 1 |  | $+$ | $+$ | $+$ | $+$ | 1.000
$\biggl{\\{}$ | $L_{x}$ | 2 | $x$ | $+$ | $-$ | $-$ |  | 0.9996
$L_{y}$ | 3 | $y$ | $-$ | $+$ | $-$ |
| $D$ | 4 |  | $-$ | $-$ | $+$ | $+$ | 0.9981
$\biggl{\\{}$ | $BL_{x}$ | 5 | 7 | $x$ | $+$ | $-$ | $-$ |  | 0.9103
$BL_{y}$ | 6 | 8 | $y$ | $-$ | $+$ | $-$ |
| $C$ | 7 | 5 |  | $+$ | $+$ | $+$ | $-$ | 0.9870
| $C^{\ast}$ | 8 | 6 |  | $+$ | $+$ | $-$ | $-$ | 0.9642
| $D^{\ast}$ | 9 |  | $-$ | $-$ | $-$ | $+$ | 0.9129
$\biggl{\\{}$ | $BL_{x}^{\ast}$ | 10 |  | $+$ | $-$ | $+$ |  | 0.9741
$BL_{y}^{\ast}$ | 11 |  | $-$ | $+$ | $+$ |
| $M^{\ast}$ | 12 |  | $+$ | $+$ | $-$ | $+$ | 0.9297
$\biggl{\\{}$ | $L_{x}^{\dagger}$ | 13 | 19 | $x$ | $+$ | $-$ | $-$ |  | 0.6037
$L_{y}^{\dagger}$ | 14 | 20 | $y$ | $-$ | $+$ | $-$ |
| $M^{\dagger}$ | 15 | 16 |  | $+$ | $+$ | $+$ | $+$ | 0.8949
| $L_{z}$ | 16 | 17 | $z$ | $-$ | $-$ | $+$ | $-$ | 0.6461
| $C^{\dagger}$ | 17 | 13 |  | $+$ | $+$ | $+$ | $-$ | 0.9446
$\biggl{\\{}$ | $K_{x}^{\dagger}$ | 18 | 14 |  | $+$ | $-$ | $+$ |  | 0.9092
$K_{y}^{\dagger}$ | 19 | 15 |  | $-$ | $+$ | $+$ |
| $D^{\dagger}$ | 20 | 18 |  | $-$ | $-$ | $-$ | $-$ | 0.8413
Table 1: Characteristics of the lowest 20 POD modes of either decomposition,
$\varphi_{E,n_{E}}$ and $\varphi_{D,n_{D}}$ with $n_{E/D}\leq 20$. Brackets on
the left designate degenerate pairs. Mode labels in group I are consistent
with those previously introduced by Soucasse  _et al._ (2019, 2020, 2021),
with additional modes designated using "$\dagger$". Group II shows mode
indices in the energy ($n_{E}$) and dissipation ($n_{D}$) optimized
decompositions; a single number is given when $n_{E}=n_{D}$. Group III lists
non-vanishing angular momentum components. Group IV shows symmetries ($+$) and
anti-symmetries ($-$) with respect to each of the four basic symmetry
operations $S_{x}$, $S_{y}$, $S_{z}$, and $S_{d}$, corresponding to
eigenvalues $\mu_{k}$ of the symmetry operation, cf.
$S_{k}\varphi=\mu_{k}\varphi$. Group V shows the overlap of modes from the two
decompositions as defined in (36) and shown in Figure 3. Since the overlap is
the same for each pair of modes in a degenerate set only a single value is
given for those. Boldface text denotes outlying overlap values. Rows are
grouped to aid reading.

The overlaps $\Gamma(n_{E},n_{D})$ for $n_{E},n_{D}\leq 20$ are shown in
Figure 3, where we have also highlighted pairs of modes for which
$G[\varphi_{E,n_{E}}]=G[\varphi_{D,n_{D}}]$. Of such pairs, several modes
($n_{E/D}\leq 4$ and $9\leq n_{E/D}\leq 12$) preserve their indices between
the decompositions, and exhibit large overlaps; these appear in Figure 3 as
framed red squares positioned on the diagonal. A number of modes are indexed
differently, but still exhibit a large overlap (framed red squares off the
diagonal); the significance of the reordering of energy modes 5–8 is discussed
below. Lastly, a few of the modes considered show only moderate overlap
(framed non-red squares). Non-framed squares denote overlaps between modes
designated with different labels. The overlap matrix is sparse, as the overlap
vanishes for mode pairs with non-identical symmetry signatures. Taken
together, this points to a large degree of similarity between the
decompositions, which will be investigated further in Section IV.3. It should
be emphasized, however, that the similarity between the decompositions tends
to diminish when considering modes higher than the ones investigated here.

Figure 3: Non-vanishing overlap magnitudes $\Gamma(n_{E},n_{D})$ between
energy and dissipation modes with $n_{E},n_{D}\leq 20$, as defined in (36).
Pairs with matching labels $G[\varphi_{E,n_{E}}]=G[\varphi_{D,n_{D}}]$ are
marked with black frames, and the diagonal $n_{E}=n_{D}$ is shown as a dashed
line.

A notable reordering of modes is the promotion of the $C$ and $C^{\ast}$ modes
from indices 7 and 8 in the energy decomposition to indices 5 and 6 in the
dissipation decomposition, ahead of the degenerate boundary roll mode pair
$BL_{x/y}$. The promotion is moderately significant, in that the eigenvalue of
the $BL_{x/y}$ modes exceed that of the $C$ modes by about 16% in the energy
decomposition, whereas the eigenvalue of the $C$ mode exceeds that of the
$BL_{x/y}$ mode by $BL_{x/y}$ modes by 37% in the dissipation decomposition.
The $C$ mode was found by Soucasse  _et al._ (2019) to have a destabilizing
effect on the LSC, and the promotion of this mode could thus be taken as an
indication that certain dynamically important structures are captured more
efficiently by the dissipation decomposition than by the energy decomposition.

For the modes $n_{E/D}\leq 20$ the ordering is to a large extent preserved
between the decompositions, although this trend breaks down when considering
higher modes. Lee and Dowell (2020) found that ordering energy-optimized modes
for a two-dimensional lid driven cavity flow by their modal contribution to
mean enstrophy led to only a small amount of reordering compared to the native
ordering of modes by energy, with the least amount of reordering found for
low-index modes. Reordering was found to increase when increasing
$\mathrm{Re}$. To the same effect, Olesen  _et al._ (2023) showed for a
channel flow that both the mean dissipation contribution of energy modes and
the mean energy contribution of dissipation modes tended to decrease with mode
number. In both cases this was taken as an indication that the flow in
question was as a whole dominated by a few structures spanning all relevant
scales, and that the lowest POD modes were capturing structures on this range
of scales. Such a lack of separation of scales could be a consequence of the
moderate Rayleigh number of the flow investigated in this study.

###  IV.2 Large-scale organization of energy modes

The energy decomposition of the data set used in this work has been considered
by Soucasse  _et al._ (2019, 2020, 2021), who described the structures
contained in the lowest several modes along with the dynamics arising from a
ROMs based on these modes. In this section we summarize and expand on the
structural analysis, so as to gain a point of comparison for subsequent
investigation of the corresponding dissipation modes.

The large-scale organization of the lowest 20 modes of each decomposition is
shown in Figures 4–5 as temperature-colored heat flux isosurfaces along with
stream lines. The lowest 20 energy modes are shown in order of increasing mode
number, each juxtaposed with the dissipation mode sharing the label of the
energy mode. Modes marked by "⋆" (not to be confused with "∗") represent one
of a degenerate pair, for which the complimentary mode (not shown) can be
obtained by rotating the mode shown by $\frac{\pi}{2}$ around the central
vertical axis.

Figure 4: Energy-optimized POD modes $\varphi_{E,1}$–$\varphi_{E,10}$
juxtaposed with their similar dissipation-optimized POD modes
$\varphi_{D,1}$–$\varphi_{D,10}$, shown as isocontours of scaled convective
heat flux corresponding to
$\Phi_{\mathrm{iso}}=0.05\times\max_{\Omega}\left(\varphi^{w}\varphi^{\theta}\right)$
and colored by scaled temperature
$\theta_{\mathrm{scale}}=\varphi^{\theta}/\max_{\Omega}\left(\left|\varphi^{\theta}\right|\right)$.
For modes marked with $\star$ ($c,d,g,h,o,p$) only one mode of a degenerate
pair is shown; the other mode in the pair is obtained by rotating the mode
shown by $\frac{\pi}{2}$ around the central vertical axis. Figure 5: Energy-
optimized POD modes $\varphi_{E,12}$–$\varphi_{E,20}$ juxtaposed with their
similar dissipation-optimized POD modes $\varphi_{D,12}$–$\varphi_{D,20}$,
shown as isocontours of scaled convective heat flux corresponding to
$\Phi_{\mathrm{iso}}=0.05\times\max_{\Omega}\left(\varphi^{w}\varphi^{\theta}\right)$
and colored by scaled temperature
$\theta_{\mathrm{scale}}=\varphi^{\theta}/\max_{\Omega}\left(\left|\varphi^{\theta}\right|\right)$.
For modes marked with $\star$ ($c,d,k,l$) only one mode of a degenerate pair
is shown; the other mode in the pair is obtained by rotating the mode shown by
$\frac{\pi}{2}$ around the central vertical axis.

The $M$, $L_{x/y}$, $D$, $BL_{x/y}$ and $C$ modes
($\varphi_{E,1}$–$\varphi_{E,7}$) shown in Figure 4$a,c,e,g,i$, are identical
to those obtained by Soucasse  _et al._ (2019), who used the same data set as
in the present work. The $M$ mode is approximately equal to the mean of the
field taken over the full data set, covering all possible orientations of the
LSC. It consists of a pair of counter-rotating toroidal structures which
captures most of the thermal stratification. The $L_{x/y}$ modes form a
degenerate pair, each mode describing a single roll along a principal axis;
they combine to form the diagonal rolls characterizing the LSC. The $D$ mode
consists of eight rolls, and transports fluid between the corners while
stabilizing circulation along the diagonal. Together, the $M$, $L_{x/y}$, and
$D$ modes form the LSC. The $BL_{x/y}$ and $C$ modes represent the dominant
fluctuations superposed on the LSC. Each of the $BL_{x/y}$ modes corresponds
to a co-rotating pair of longitudinal rolls that connect the core to the
boundary layers and modulate the $L_{x/y}$ modes. The $C$ mode corresponds to
a corner structure that brings in fluid from one roll to another, and plays an
important role in the LSC reorientations. The reader is referred to the paper
by Soucasse  _et al._ (2019) for a more detailed discussion of these modes.

Additionally, modes similar to the $C^{\ast}$, $D^{\ast}$, $BL_{x/y}^{\ast}$,
$M^{\ast}$, and $L_{z}$ modes ($\varphi_{E,8}$–$\varphi_{E,12}$ and
$\varphi_{E,16}$), shown in Figures 4$k,m,o$ and 5$a,g$, were identified and
discussed for $\mathrm{Ra}=10^{6}$ and $10^{8}$ by Soucasse  _et al._ (2021).
The $C^{\ast}$ and $D^{\ast}$ modes are somewhat similar to the $C$ and $D$
modes, respectively, except that each mode has rolls extending along the full
height of the cell, where the rolls in their non-asterisked counterparts are
confined to the upper and lower half-cells. Correspondingly, as seen in Table
1 the $C^{\ast}$ and $D^{\ast}$ modes are antisymmetric with respect to
$S_{z}$, while the $C$ and $D$ modes are symmetric with respect to $S_{z}$;
symmetries are otherwise similar between the respective modes. The
$BL_{x/y}^{\ast}$ modes form a pair of counter-rotating rolls (compare with
the _co_ -rotating rolls in the $BL_{x/y}$ modes). The angular momenta of the
rolls cancel out, meaning that the $BL_{x/y}$ modes as a whole carry no
angular momentum. As before the $BL_{x/y}$ and $BL_{x/y}^{\ast}$ modes are
linked by similar symmetries, except for symmetry vs antisymmetry with respect
to $S_{z}$. Likewise, the $M^{\ast}$ mode is an $S_{z}$-antisymmetric
counterpart to the $M$ mode, consisting of a single thermally unstratified
toroidal roll extending the full height of the cell. Lastly, the $L_{z}$ mode
consists of a vertical roll surrounded by horizontal rolls. This mode is the
only one out of those considered capable of supporting a non-zero vertical
angular momentum component, which is needed to resolve fluctuations in this
component away from zero.

The last set of energy modes to be considered has not been discussed in prior
literature, and consists of the modes $L_{x/y}^{\dagger}$, $M^{\dagger}$,
$C^{\dagger}$, $K_{x/y}^{\dagger}$, and $D^{\dagger}$
($\varphi_{E,13}$–$\varphi_{E,15}$ and $\varphi_{E,17}$–$\varphi_{E,20}$),
shown in Figure 5$c,e,i,k,m$. The $L_{x/y}^{\dagger}$ modes form a degenerate
pair which shares symmetry signatures with the $L_{x/y}$ and $BL_{x/y}$ pairs.
The $L_{x}^{\dagger}$ mode features a central roll along the $x$ axis, which
buckles around near the corners to form eight rolls perpendicular to the main
roll, one near each corner. Together the central roll and the cross rolls give
rise to temperate convection zones on the flanks of the central roll, while
the cross rolls drive cold and hot convection zones along the vertical edges
of the cell. The $M^{\dagger}$ mode shares its symmetry signature with the
mean mode $M$. It features a pair of counter-rotating toroidal structures,
localized in the lower and upper half of the cell. These produce vertically
split convection zones along the vertical center line, which are surrounded by
convection regions along the top and bottom edges as well as opposite
convection zones towards the middle of the vertical edges. Of all modes
considered here only the $M$ and $M^{\dagger}$ modes have non-vanishing
thermal energy at the top and bottom walls, implying that $M^{\dagger}$ serves
as a first correction to the $M$ mode in resolving the thermal stratification.

The $C^{\dagger}$ mode shares its symmetry signature with the $C$ mode, and
consists of four vertical rolls along the vertical edges, separated by pairs
of short horizontal rolls. It moves fluid between the vertical walls and the
center of the cell. The $K^{\dagger}_{x/y}$ modes form another degenerate
pair, with symmetry signatures identical to those of the $BL_{x/y}^{\ast}$
modes. The $K_{x}^{\dagger}$ mode consists of a horizontal roll near each of
the eight corners, aligned along the $y$ axis and with alternating rotations.
These are connected by an additional set of horizontal rolls lodged between
the first and orthogonal to these. The angular momentum contribution of all
the rolls cancel in each direction, and like the $BL_{x/y}^{\ast}$ modes the
$K^{\dagger}_{x/y}$ modes therefore carry no angular momentum. Lastly, the
$D^{\dagger}$ mode is unique among the modes considered in being antisymmetric
with respect to all of the basic symmetry operations. It consists of eight
vertical rolls aligned along the vertical sides and edges. The rolls produce
eight alternatingly hot and cold convective regions, each extending the full
height of the domain.

###  IV.3 Large scale organization of dissipation modes

The dissipation decomposition is motivated by the desire to better resolve
dynamically important structures which are not well captured by the energy
decomposition. As demonstrated in Section IV.1 the results of the two
decompositions differ both in the structures contained in the modes and in the
ordering of those modes, although as evidenced by the near-unity overlaps the
structural differences are mostly modest for the modes considered. In this
section we take a closer look at the differences in large-scale organization
between the modes, while details relating to boundary layer structures are
considered in Section IV.4.

Apart from reordering, the large-scale organization of the modes corresponding
to $n_{E/D}\leq 12$ show little change between the two decompositions, as seen
from Figures 4 and 5. Several modes ($M$, $L_{x/y}$, $D$, $C$,
$BL_{x/y}^{\ast}$, and $C^{\dagger}$) are found with virtually unchanged large
scale organization; these modes also generally exhibit overlaps very close to
unity, cf. Table 1. Some changes that are notable from Figures 4 and 5 include
a deformation of rolls in the $BL_{x/y}$ modes leading to pinching of the hot
and cold convection zones at mid height and elongation of the temperate
convection zones towards the top and bottom; a narrowing of rolls in the
$C^{\ast}$ and $D^{\ast}$ modes also around mid height; and a splitting of the
rolls at the vertical edges along the mid height plane in $M^{\ast}$.

The changes in modal structure when moving from the energy decomposition to
the dissipation decomposition impact the modal contributions to reconstructed
flow quantities. Vertical modal profiles of kinetic and thermal energy,
viscous and thermal dissipation, and convective heat flux corresponding to the
contribution of individual modes to the mean profiles are obtained by
averaging individual terms in (31) along the $x$ and $y$ directions. In Figure
6 we compare such profiles for some selected modes of each decomposition. The
profiles for the $M$ and $L_{x/y}$ modes (not shown) are virtually identical
between decompositions. For the $D$, $BL_{x/y}$, and $C$ modes the profiles
are very similar among the decompositions, showing only minor variations; the
greatest deviation for those modes is found for the $BL_{x}$ mode profiles,
which are shown in Figure 6$a$-$e$. For higher modes the profiles differ
between decompositions in a mostly consistent manner, with notable exceptions
highlighted below. While the shape of profiles varies somewhat between modes,
for each of the quantities considered the profiles obtained from dissipation
modes generally contribute more than the energy modes in the region near the
top and bottom. The contributions are reversed or similar in the core of the
cell.

Figure 6: Profiles of kinetic and thermal energy, viscous and thermal
dissipation, and convective heat flux associated with $BL_{x}$ modes
$\varphi_{E,5}$ and $\varphi_{D,7}$ ($a$–$e$), $L_{x}^{\dagger}$ modes
$\varphi_{E,13}$ and $\varphi_{D,19}$ ($f$–$j$), and $L_{z}$ modes
$\varphi_{E,16}$ and $\varphi_{D,17}$ ($k$–$o$). The profiles are symmetric
around the mid plane, and only $z$ values corresponding to the lower half of
the cell are shown. A dotted line marks the extent of the boundary layer,
$z=\delta=0.03$.

In the $L_{x}^{\dagger}$ dissipation mode in Figure 5$d$ the central roll is
widened compared to the energy mode in Figure 5$c$, and the cross rolls align
with the horizontal diagonals instead of being orthogonal to the roll axis.
This also causes a wider temperate convection zone, splitting the edge
convection zones vertically. The diagonally aligned cross rolls produce
additional convection zones above and below the central roll. The splitting of
the edge convection zone diminishes the thermal fluctuation and the heat flux
near mid-height of the cell, resulting in a reduced contribution from the
dissipation mode to the vertical mean thermal energy and convective heat flux
profiles shown in Figure 6$g,j$. The substantially different structures and
heat flux profiles found for the energy and dissipation $L_{x/y}^{\dagger}$
modes ($\varphi_{E,13}$ and $\varphi_{D,19}$) are reflected in the rather
modest overlap between the modes, at $\Gamma(13,19)=0.6037$.

In the $M^{\dagger}$ dissipation mode in Figure 5$f$ the opposite convection
zones near the corners are strengthened compared to the energy mode, and they
extend further towards the horizontal walls and mid plane of the domain. The
convection zone along the top and bottom edges largely disappear.

The $L_{z}$ modes shown in Figure 5$g,h$ exhibit substantial dissimilarities
between the decompositions. In the dissipation $L_{z}$ mode the vertical roll
is much narrower than in the energy mode, and the mode is instead dominated by
the surrounding horizontal rolls. The dissimilarity is also supported by the
mode overlaps; while the greatest overlap for the energy $L_{z}$ mode
($\varphi_{E,16}$) found among the modes investigated here was with the
$L_{z}$ dissipation mode ($\varphi_{D,17}$) and vice versa, at
$\Gamma(16,17)=0.6461$, both exhibit slightly greater overlaps with different
modes of higher indices; in particular, $\Gamma(30,17)=0.7053$, and
$\Gamma(16,72)=0.7066$. For all other modes considered here the maximum
overlapping mode was found within the set $n_{E/D}\leq 20$. The dissimilarity
between the $L_{z}$ modes is also found in the vertical profiles in Figure
6$k$–$o$, which differ rather substantially between the two modes. The kinetic
energy profile of the $L_{z}$ energy mode has a marked peak at mid-height,
which is entirely absent in the $L_{z}$ dissipation mode. Conversely, the
dissipation mode thermal energy profile peaks within the boundary layer, while
the corresponding contribution from the energy mode is much smaller throughout
the domain. This suggests that the energy mode is mainly kinetic in nature,
while the dissipation mode is mainly thermal. The dissipation mode
contribution to both viscous and thermal dissipation near the wall is notably
stronger than that of the energy mode, which agrees qualitatively with the
trends displayed by other modes. Finally, due to the diminished thermal
fluctuations in the energy mode, the dissipation mode contributes
significantly more to convective heat flux throughout the cell (except at mid
height) than does the energy mode.

While the primary rolls in the $K_{x/y}^{\dagger}$ energy mode (Figure 5$k$)
are aligned along the $y/x$ axis and joined by orthogonal horizontal rolls,
for the corresponding dissipation mode (Figure 5$l$) the first set of rolls
are instead aligned along the horizontal diagonals, and cross-rolls are
notably wider. In the $D^{\dagger}$ dissipation mode (Figure 5$n$) the rolls
along the vertical edges are split into pairs of co-rotating rolls.

To summarize, there is generally a large degree of similarity in the large-
scale organization between analogous modes of different decompositions. This
is especially striking for the lower modes in Figure 4, whereas the higher
modes in Figure 5 tend be more distinct, in agreement with the overlaps
summarized in Table 1. Part of the explanation for these similarities is the
symmetry-imposed constraints on mode mixing. The two bases produced by the
decompositions span the same Hilbert space, indicating that each mode of one
basis can be expanded as a linear combination of modes from the other basis.
As was mentioned in Section IV.1, the overlap vanishes for pairs of modes
which do not share their configuration of symmetry eigenvalues, implying that
only isosymmetric modes enter in this expansion. This prohibits mixing between
different families of isosymmetric modes, constraining the dissimilarity
between modes produced by different decompositions.

Even within such families we find mixing to be limited, however. Among all
modes considered the $L_{x/y}^{\dagger}$ modes have the smallest overlap
between analogous modes, with $\Gamma(13,19)=0.6037$ for the $L_{x}^{\dagger}$
modes $\varphi_{E,13}$ and $\varphi_{D,19}$. The modes exhibit clear
differences in large-scale structures between the decompositions, as shown in
Figure 5$c,d$. This dissimilarity is due to mixing within the isosymmetric
family $\\{L_{x/y},BL_{x/y},L_{x/y}^{\dagger},\ldots\\}$, including a
substantial overlap of $\Gamma(5,19)=0.5828$ between the $BL_{x}$ mode
$\varphi_{E,5}$ and the $L_{x}^{\dagger}$ mode $\varphi_{D,19}$. On the other
hand, little mixing happens between the $L_{x/y}$ and $BL_{x/y}$ modes, with
$\Gamma(2,7)=0.1255$ and $\Gamma(5,2)=0.0211$. Other isosymmetric families
include $\\{M,M^{\dagger},\ldots\\}$, $\\{C,C^{\dagger},\ldots\\}$, and
$\\{BL_{x/y}^{\ast},K_{x/y}^{\dagger},\ldots\\}$; naturally, the remaining
modes also belong to isosymmetric families, although those modes are each the
sole representative of their respective families within the set of modes
considered here. Mixing within the families considered here is generally
rather limited, which results in large-scale organization being mostly stable
across decompositions. Again, however, we remark that while this is true for
modes of low indices, the stability diminishes when considering higher index
modes.

The stability up to reordering of large-scale organization of modes under
varying decompositions was also found by Podvin and Sergent (2015) when
comparing kinetic and thermal energy modes for a two-dimensional square
Rayleigh-Bénard cell. They attributed the similarity between the
decompositions to the close coupling between velocity and temperature, as the
motion is entirely driven by temperature fluctuations. In the present work,
the observation that the first several modes of both energy and dissipation
decompositions describe structures much more similar than what is required by
symmetry constraints supports the hypothesis of a coupling between energetic
and dissipative structures discussed in Section IV.1. Due to the modest
Rayleigh and Reynolds numbers $\mathrm{Ra}=10^{7}$ and $\mathrm{Re}=651$, as
defined in (1) and (32), we expect only a moderate separation between
energetic and dissipative scales, leading to a tighter coupling between the
associated structures. We would thus expect a greater amount of mode
dissimilarity and reordering given higher values of $\mathrm{Ra}$ and
$\mathrm{Re}$.

We note that the remarkable stability of the LSC modes $M$, $L_{x/y}$, and $D$
across decompositions may further be attributed to the degree to which the
structures they encode dominate the flow.

The stability of large-scale structures between decompositions
notwithstanding, certain subtle but consistent changes are found. Dissipation
modes tend to contribute more to vertical mean profiles of kinetic and thermal
energy, viscous and thermal dissipation, and heat flux in the boundary layer
region than do energy modes. This trend is seen in Figure 6, and it is also
found in the corresponding profiles of other modes (not shown). Conversely, in
the core of the domain the contributions are either similar between
decompositions or diminished for dissipation modes. A similar observation was
made by Olesen  _et al._ (2023) for profiles of kinetic energy and of viscous
dissipation reconstructed using decompositions optimized with respect to these
quantities. This was ascribed to the different spatial distributions of the
quantities with respect to which the decompositions were optimized. Kinetic
energy is located mainly in the core of the cell, viscous dissipation is
concentrated along the walls (both horizontal and vertical), and thermal
energy and thermal dissipation peak near the top and bottom. The energy and
dissipation decompositions each give equal weight to the two terms of energy
and dissipation, respectively. This causes the energy decomposition to give
roughly equal emphasis to the core and boundary regions, and the dissipation
decomposition to emphasize structures in the boundary layer region.

Finally, we note that Figure 6 demonstrates the failure of the assumption used
in isotropic eddy-diffusivity models, as pointed out by Hanjalić (2002), that
heat flux locally scales as $\Phi\sim\sqrt{\left\langle
E_{u}E_{\theta}\right\rangle}$. For example, in the energy $BL_{x/y}$ mode
both $E_{u}$ and $E_{\theta}$ show local minima at $z=0.1$ (Figure 6$a,b$),
while $\Phi$ increases monotonically with $z$ up to around $z=0.2$ (Figure
6$e$). This is attributed to the unequal distribution of kinetic energy among
the available degrees of freedom.

###  IV.4 Boundary layer structures

As discussed in Section I, boundary layers are of special interest in
Rayleigh-Bénard convection. In this section we consider the boundary layer
structures resolved by energy and dissipation modes. Since viscous and thermal
dissipation are both concentrated in the boundary layer region we expect this
region to be particularly sensitive to the choice of decomposition. The
analysis presented here highlights some subtle but potentially important
differences between the decompositions.

In Figure 7 we show horizontal cross sections of the bottom boundary layer
($z=0.01$) of selected modes. The cross sections are colored by heat flux, and
also show the direction of the horizontal components of the flow field. The
LSC modes $M$, $L_{x/y}$, and $D$ were found in Section IV.3 to be largely
stable under the change of decomposition, and the cross sections obtained from
these modes (not shown) exhibit no discernible differences between the
decompositions. For the $BL_{x}$ modes shown in Figure 7$a,b$ the dissipation
mode exhibits an enhanced heat flux along the edges at $y=0$ and $y=1$, as
well as in separated pairs of patches along $x=0$ and $x=1$. These correspond
to the convective regions seen in Figure 4$g,h$, reflecting the extension of
the temperate convective regions discussed in Section IV.3. Figure 7$a,b$ thus
demonstrates that the changes in these features extend well into the boundary
layer. This extension and strengthening of features into the boundary layer in
dissipation modes compared to analogous energy modes represents a general
trend that is also seen for the $L_{x}^{\dagger}$ and $L_{z}$ modes shown in
Figure 7$e$–$h$, as well as for most of the modes $n_{E,D}\leq 20$ not shown
in Figure 7.

Figure 7: Horizontal cross sections taken near the bottom ($z=0.01$) of the
$BL_{x}$ modes $\varphi_{E,5}$ ($a$) and $\varphi_{D,7}$ ($b$), the $C$ modes
$\varphi_{E,7}$ ($c$) and $\varphi_{D,5}$ ($d$), the $L_{x}^{\dagger}$ modes
$\varphi_{E,13}$ ($e$) and $\varphi_{D,19}$ ($f$), and the $L_{z}$ modes
$\varphi_{E,16}$ and $\varphi_{D,17}$. Slices are colored by scaled convective
heat flux,
$\Phi_{\mathrm{scale}}=\varphi^{w}\varphi^{\theta}/\max_{\Omega}\left({\varphi^{w}\varphi^{\theta}}\right)$.
Arrows show the direction of the in-plane velocity.

A deviation from this trend is seen in the horizontal cross sections of the
bottom boundary layer of the $C$ modes shown in Figure 7$c,d$. We note little
difference between the two decompositions, although a slight narrowing of the
wedge-shaped regions of positive (upward) heat flux can be discerned. The
cross sections show a high level of activity in the boundary layer for both
modes, consisting of said wedge-shaped zones marginally separated from narrow
stretches of slight negative (downward) heat flux along the edges. Comparing
with Figure 4$i,j$ we see that the wedge-shaped zones correspond to downwards
movement of cold fluid near the $x=0$ and $x=1$ walls, and upwards movement of
warm fluid near the $y=0$ and $y=1$ walls. The zones of negative heat flux are
found to correspond to inversions of the vertical velocity, implying large
velocity gradients in the margins separating the regions. It was noted in
Section IV.1 that the $C$ mode is important in promoting LSC reorientations,
and that its promotion in the dissipation decomposition could reflect an
improved representation of the flow dynamics in this decomposition; in this
light it is interesting to note that the mode itself appears to change
relatively little in its large-scale organization as well as in its boundary
layer layout between the decompositions.

The $L_{x/y}^{\dagger}$ and $L_{z}$ modes were found in Section IV.3 to differ
substantially in their large scale organization between the decompositions.
Horizontal cross sections of the $L_{x}^{\dagger}$ modes are shown in Figure
7$e,f$. Moving from the energy to dissipation decomposition brings about
changes in the cross section reminiscent of those seen for $BL_{x}$ in Figure
7$a,b$, with regions of positive heat flux towards the corners being stronger
and more well-defined in the dissipation mode, as well as additional regions
appearing around $x=0.5$. This corresponds to the additional regions of
positive convection appearing in Figure 5$d$ compared to Figure 5$c$ as
discussed in Section IV.3.

We show horizontal cross sections of the $L_{z}$ mode in Figure 7$g,h$. While
the qualitative pattern observed changes little between the decompositions,
the convective regions are substantially enlarged and strengthened in the
dissipation mode. Comparing with Figure 5$g,h$ we see that the stronger
vertical movement in the dissipation mode compared to the energy mode causes
the convective regions to reach further towards the top and bottom.

To summarize, there is an general trend towards increased boundary layer
activity in the dissipation modes compared to energy modes, expressed through
more extensive and well-defined zones of substantial convective heat flux in
these modes.

The theory proposed by Grossmann and Lohse (2000, 2004) referred to in Section
I suggests that important modeling parameters include the relative
contributions of viscous and thermal dissipation in the boundary layers and in
the whole cell. This motivates analyzing the modes in terms of these
parameters. In Figure 8 we compare on a per-mode basis the ratio of boundary
layer contributions to the full-cell viscous (Figure 8$a$) and thermal
dissipation (Figure 8$b$) for energy and dissipation modes, as well as the
ratio of thermal to viscous dissipation (Figure 8$c$), for the LSC modes $M$,
$L_{x/y}$, and $D$; the $BL_{x/y}$ and $C$ modes describing the leading
fluctuation terms; and the $L_{z}$ mode which was found to deviate from
several of the trends established in Sections IV.3 and IV.4. The values
averaged over the modes (excluding the $M$ modes) are shown by dashed lines,
and eigenvalues are shown to indicate the relative importance of the modes in
the reconstructed flow.

Figure 8$a$ shows the fraction of total viscous dissipation located in the top
and bottom boundary layers. These contributions amount to 0.2–0.5, except for
the $M$ mode where it is around 0.1. Averaged over modes $1<n_{E/D}\leq 20$
(i.e., excluding $M$ modes) the mean fraction is 0.32 for energy modes and
0.41 for dissipation modes. The boundary layer contributions are enhanced for
most dissipative modes compared to energy modes, except for the $M$ and
$L_{x/y}$ modes where ratios are nearly identical; a more substantial
difference is observed for the $L_{z}$ modes. The trend of enhanced boundary
layer contributions for dissipation modes is maintained by the remaining modes
$n_{E/D}\leq 20$ (not shown), with the difference between the fractions for
dissipation modes and for energy modes tending to increase for higher modes.
The enhanced boundary layer contributions from dissipation modes reflect the
increased boundary layer activity that was found for most of the modes
considered in Figure 7.

Figure 8: Modal viscous $(a)$ and thermal $(b)$ dissipation fraction in
boundary layer, and modal ratio of thermal to viscous dissipation $(c)$ for
the modes up to $C$ (left axes). Dashed lines show the mean values for each
decomposition, excluding the $M$ modes. The relative importance of modes is
given by eigenvalues, which are plotted with lines (right axes).

The fraction of thermal dissipation contributed from the boundary layer is
shown in Figure 8$b$, where we also find approximately constant ratios around
0.6–0.8 among the modes considered, with the exception of the $L_{z}$ energy
mode for which the fraction is around 0.4. The average fractions for modes
$1<n_{E/D}\leq 20$ are 0.69 for energy modes and 0.76 for dissipation modes.
The thermal dissipation boundary layer fraction is thus generally on the order
of twice the viscous dissipation boundary layer fraction. This difference is
likely caused by the exclusion of the vertical boundary layers, which
contribute appreciable viscous dissipation and little thermal dissipation. The
difference in thermal dissipation boundary layer fraction between energy and
dissipation modes is generally small compared to the difference in the viscous
dissipation boundary layer fraction, with the most substantial difference
again found for the $L_{z}$ modes. The increased boundary layer activity for
dissipation modes compared to energy modes thus has a smaller relative impact
on the thermal dissipation boundary layer fraction than on the viscous
dissipation boundary layer fraction.

The ratios of thermal to viscous dissipation over the whole cell are shown in
Figure 8$c$. The $M$ modes stand out with a ratio of 5.1, their dissipation
being strongly dominated by the thermal contribution. This supports the
finding by Soucasse  _et al._ (2019) that the $M$ modes capture the thermal
stratification of the flow. The remaining modes (including the higher modes
not shown here) are dominated by viscous dissipation, with average ratios for
modes $1<n_{E/D}\leq 20$ of 0.26 for energy modes and 0.34 for dissipation
modes. For the modes up to and including $C$ ratios are mostly consistent
between energy and dissipation modes. Variability between the decompositions
increases for higher modes, in agreement with the trend of decreasing
similarity found in Sections IV.1 and IV.3. The $C$ modes, along with a few
higher modes, have substantially higher ratios than those found for $L_{x/y}$,
$D$, and $BL_{x/y}$ modes. It is possible that the destabilizing effect of the
$C$ modes on the LSC is linked to the relative importance of thermal
dissipation in this mode. The $L_{z}$ energy mode has a remarkably small
ratio, whereas the ratio of the corresponding dissipation mode matches the
average dissipation mode ratio. This matches the observation from Figure 6 in
Section IV.3 that the energy mode was mainly kinetic in nature, whereas the
dissipation mode was mainly thermal.

We see from the preceding analysis that dissipation modes in general tend to
have a greater fraction of both viscous and thermal dissipation localized in
the boundary layers than do energy modes. This is a consequence of the
increased activity in the boundary layer for dissipation modes. Also, while
modes of either decomposition are as a rule dominated by viscous dissipation,
dissipation modes tend toward a higher ratio of thermal to viscous
dissipation.

##  V Modal convergence

The energy and dissipation POD bases each allow reconstruction of second order
mean quantities using (31). Here we investigate the convergence of such
reconstructions of kinetic and thermal energy, viscous and thermal
dissipation, and convective heat flux. We consider the convergence of
quantities integrated over the full domain ($\Omega$) as well as over the
region consisting of the horizontal top and bottom boundary layers,
$BL=[0;1]\times[0;1]\times\left([0;\delta]\cup[1-\delta;1]\right)$,

$\displaystyle\Pi^{p,\Omega^{\prime}}(n)$
$\displaystyle=\int_{\Omega^{\prime}}\Pi^{p}_{n}\,dx\,,$ (37)

where $\Omega^{\prime}\in\left\\{\Omega,BL\right\\}$. Here we use
$\delta=0.03$ as an estimate of the boundary thickness, based on the kinetic
and thermal boundary layer thicknesses $\delta_{u}=0.027$ and
$\delta_{\theta}=0.031$(Shishkina  _et al._ , 2010).

We compare the convergence achieved for each integrated quantity using the
energy decomposition to that achieved with the dissipation decomposition. The
convergence of each reconstruction is shown in Figure 9. Considering the
cumulative values (dotted lines) we see that for all reconstructed quantities
the three lowest energy modes contribute the same as the corresponding
dissipation modes, as these modes are nearly identical (cf. Section IV). The
$M$ mode captures much of the thermal stratification in the flow, and
accordingly accounts for most of the thermal energy and thermal dissipation in
both the full domain and in the boundary layers.

Figure 9: Modal convergence of kinetic ($a,b$) and thermal energy ($c,d$),
viscous ($e,f$) and thermal dissipation ($g,h$), and heat flux ($i,j$) in the
full domain (top row) and in the boundary layer
($z\in[0;\delta]\cup[1-\delta;1]$ with $\delta=0.03$, bottom row). Fully drawn
lines show the per-mode contribution (left axes), and dotted lines the
cumulative value (right axes). Also shown are power laws fitted to per-mode
contributions in the interval $10\leq n\leq 10^{3}$. All values are normalized
with respect to their full reconstruction.

For the reconstruction of kinetic energy $E_{u}$ in the full domain, shown in
Figure 9$a$, the energy decomposition converges slightly faster than the
dissipation decomposition, as would be expected because of the energy
optimization. However, the converse is true for the reconstruction of kinetic
energy in the boundary layer (Figure 9$b$). The convergence of thermal energy
(Figure 9$c,d$) is very rapid, and the two decompositions show little
difference in convergence rates, in the full domain as well as in the boundary
layer. A very small lead can be ascertained for the energy decomposition in
the full domain.

For the reconstructions of viscous ($\varepsilon_{u}$) and thermal dissipation
($\varepsilon_{\theta}$) in the full domain (Figure 9$e,g$) the two
decompositions are again very similar, with a marginal advantage to the
dissipation decomposition; the same applies to the reconstruction of the
thermal dissipation in the boundary layer (Figure 9$h$). However, in the
reconstruction of viscous dissipation in the boundary layer (Figure 9$f$) the
dissipation decomposition exhibits a substantially more rapid convergence than
the energy decomposition.

The reconstruction of convective heat flux $\Phi$ in the full domain and in
the boundary layer (Figure 9$i,j$) repeat the pattern for kinetic energy, with
faster convergence using the energy decomposition in the full domain, and
using the dissipation decomposition in the boundary layer. However, given the
very limited amount of convective heat flux in the boundary layer the
convergence rate of convective heat flux in the boundary layer is of limited
importance.

The overall picture emerging from this is that in the boundary layer the
dissipation decomposition reconstructs all quantities investigated at least as
efficiently as the energy decomposition, while results are more mixed when
considering the full domain. This suggests that the dissipation decomposition
not only gives the most efficient representation of total dissipation, for
which it is optimized, but provides an overall improvement in the spectral
resolution of the boundary layer region.

The per-mode contributions for each quantity in Figure 9 approximately follow
power law decays in the range $10\leq n\leq 10^{3}$. The exponents have been
determined by fitting power functions to per-mode contributions in this
interval, and are summarized in Table 2. In agreement with the trends
identified above, the energy decomposition for the full cell exhibits faster
convergence for kinetic energy and convective heat flux compared to the
dissipation decomposition, while the convergence rates are nearly identical
for thermal energy. Conversely, the dissipation decomposition gives the faster
convergence rate for both viscous and thermal dissipation. In the boundary
layer the dissipation decomposition yields faster convergence for all
quantities than does the energy decomposition. This also includes thermal
energy, for which no substantial difference on convergence rate was apparent
from Figure 9$d$ due to the very fast initial convergence caused by the
thermally stratified $M$ mode. We observe that all quantities converge faster
in the boundary layer than in the full domain, with the exception of the
energy reconstruction of the convective heat flux, although as mentioned the
convective heat flux contribution from the boundary layer is negligible.
Finally, we note that the dissipation-based reconstructions in the full cell
converge with approximately the same rate ($\alpha\approx-0.8$) for all
quantities. This might be taken as an indication that the dissipation
decomposition is able to capture a wide range of scales more uniformly than is
the case for the energy decomposition.

$\Omega^{\prime}$ | $p$ | $E_{u}$ | $E_{\theta}$ | $\varepsilon_{u}$ | $\varepsilon_{\theta}$ | $\Phi$
---|---|---|---|---|---|---
$\Omega$ | $E$ | -0.98 | -0.79 | -0.75 | -0.59 | -0.95
$D$ | -0.84 | -0.78 | -0.81 | -0.79 | -0.79
$BL$ | $E$ | -1.08 | -0.83 | -0.99 | -0.79 | -0.71
$D$ | -1.27 | -1.04 | -1.15 | -1.02 | -0.91
Table 2: Asymptotic decay exponents $\alpha$ for the convergence plots shown
in Figure 9. The exponents are obtained by fitting power functions
$\Delta\Pi^{p,\Omega^{\prime}}\sim n_{p}^{\alpha}$ in the range $10\leq n\leq
10^{3}$ to per-mode contributions to each reconstruction.

We may also reconstruct full profiles using (31). Figure 10 shows the mode-by-
mode cumulative reconstruction of vertical profiles (with the mean field
profile subtracted) of kinetic and thermal energy, viscous and thermal
dissipation, and convective heat flux. In the figure, the color of each
reconstructed profile corresponds to the number of modes included for that
profile, and differences in the shape of each isochrome between the respective
decompositions signify local differences in convergence of the profile. Note
that each full reconstruction ($n=N$) exactly matches the corresponding mean
profile, illustrating the completeness of both bases.

Figure 10: Cumulative reconstructions of vertical profiles, with the mean
field profile $\Pi\left[\left\langle q\right\rangle\right]$ subtracted, of
kinetic energy ($a,b$), thermal energy ($c,d$), viscous dissipation ($e,f$),
thermal dissipation ($g,h$), and convective heat flux ($i,j$) reconstructed
using energy (top row) and dissipation decomposition (bottom row). Shown are
each profile using the first $n$ modes for each $n\leq 40$, every 10 profiles
for $40<n\leq 800$, and every 100 profiles for $800<n<$16\,000$$. Colors show
the number of modes included. The mean profile is shown as a dashed line, and
a dotted line at $z=\delta=0.03$ marks the extent of the boundary layer.
Profiles corresponding to $n=20$ are emphasized with darker lines, cf. the
legend on the upper right. The profiles for $n=20$ from the dissipation
decomposition (lower row plots) are reproduced in the upper row plots for
comparison. All profiles are vertically symmetric and are shown only for the
lower half of the domain.

As discussed above the lowest four modes are nearly identical, and the
profiles corresponding to modes 2, 3, and 4 in Figure 10 differ only
imperceptibly between the decompositions. The profiles corresponding to mode 1
are nearly identical to the subtracted mean field profile, and are not shown.
For higher modes differences can be ascertained by studying the details in
each plot. When a given isochrome extends further to the right in the energy
decomposition profile than in the corresponding dissipation decomposition
profile for a given vertical position it implies a faster convergence of the
energy decomposition at that position, and vice versa.

However, extracting meaningful differences in the reconstruction by inspecting
isochromes in Figure 10 can be difficult. Therefore, profiles corresponding to
$n=20$ have been marked in Figure 10, and the corresponding profiles from the
two decompositions are shown together for comparison in the top row plots in
Figure 10. Comparing the emphasized reconstructions for $n=20$ allows us to
see the variations in convergence with the vertical coordinate more clearly
for these particular reconstructions. The energy modes reconstruct more of the
emphasized energy and heat flux profiles near the center of the cell, while
little difference can be seen for the remaining quantities. The converse is
true near the lower wall, with more of each profile being reconstructed using
the dissipation decomposition than using the energy decomposition.

In agreement with what was found from the global convergence discussed above,
we find a trend of greater reconstruction efficiency near the boundary of the
domain using dissipation decomposition, and near the core of the domain using
the energy decomposition.

Figure 11: Ratio of reconstructed profiles using $n=20$ modes to the full mean
profile, $\frac{\Pi^{p}_{n}}{\left\langle\Pi\right\rangle}$ for $\Pi=E_{u}$
($a$), $E_{\theta}$ ($b$), $\varepsilon_{u}$ ($c$), $\varepsilon_{\theta}$
($d$), and $\Phi$ ($d$), and $p\in{E,D}$. A dotted line at $z=\delta=0.03$
marks the extent of the boundary layer. Gray regions mark the approximate
extent of the entrainment zone, $0.02\lesssim z\lesssim 0.1$. All profiles are
symmetric and shown only for the lower half of the cell.

##  VI Perspectives for modeling

An important motivation for considering the dissipation decomposition is its
potential for use as a modeling tool. Energy-based POD modes frequently serve
as the basis in Galerkin projections, in which governing equations are
projected on the chosen basis to produce a set of equations for the modal
coefficients(Rempfer, 2000). The number of coefficients, and thus the number
of equations to be solved, can be reduced by truncating the basis, projecting
the dynamics into the subspace spanned by the truncated basis. The suitability
of such ROMs relies on how well the projected dynamics approximate the full
dynamics of the flow, and whether the truncation error can be compensated for
by use of an additional modeling layer. One important source of truncation
error lies in the multi-scale nature of turbulence in combination with the
emphasis on large-scale structures assumed for energy-based POD. This is
particularly an issue in the case of high-$\mathrm{Re}$ flows, where the range
of scales can be considerable, and the coupling between energetic large-scale
structures and small-scale dissipative structures diminishes. This generally
causes an under-modeling of dissipation, for which a compensating closure
model will be required.

The dissipation POD is intended as an approach to tackle this issue. As
mentioned in the introduction, the energy-based POD corresponds to $L^{2}$
convergence of the reconstructed fields, while the dissipation-based POD
enforces $H^{1}$ convergence, which should ensure a better convergence of the
dynamics. One of the motivations for the dissipation-based POD is therefore to
capture significant dynamics more efficiently. While this is only guaranteed
asymptotically, we have seen this is already true at the lowest order. The
dissipation-based POD identifies the C mode as the most significant
modification of the LSC, which was not the case for the energy POD. Soucasse
_et al._ (2020) found the $C$ mode to be necessary for the energy POD-based
ROM in order to capture reorientations of the LSC. Linear stability analysis
of the ROM confirmed the destabilizing influence of this mode on the LSC.

More generally, we have established that dissipation-based POD provides a more
accurate reconstruction of all energy and dissipation quantities in the
boundary layers, the dynamics of which are crucial for the Rayleigh-Bénard
problem. As expected, the dissipation-based POD reconstruction of energy does
not converge as fast as the energy POD in the full cell. However, Table 2
shows that the asymptotic convergence rate is about the same for all energy
and dissipation quantities reconstructed in the full cell using the
dissipation based POD, suggesting optimality of the decomposition across the
full range of scales.

An advantage to conventional POD based Galerkin methods stems from the
orthogonality of POD modes, which leads to a decoupling of the coefficient
equations. The dissipation modes are not orthogonal with respect to the usual
inner product, creating a potential issue when applying the Galerkin
projection. However, a similar issue is faced e.g. when working in the
vorticity formulation of the Navier-Stokes equation using energy modes, as
their vorticity fields are not orthogonal. Rempfer and Fasel (1994) argued
that the vorticity fields of energy eigenmodes were sufficiently close to
vorticity eigenmodes that approximate orthogonality could be assumed, and
showed that this assumption lead to acceptable results. In Section IV.2 we
observed a general trend of similarity between energy and dissipation modes,
which gives some hope that at least an approximate orthogonality can be
achieved for isosymmetric modes, while non-isosymmetric modes are known to be
exactly orthogonal. We test this by computing the overlap matrix for the
lowest 30 dissipation modes, the entries of which are shown in Figure 12. The
modes are normalized such that the diagonal entries are unity. The off-
diagonal part of the matrix is sparse due to the symmetry-induced constraints
on mixing, and the non-zero off-diagonal entries are generally small compared
to the diagonal entries ($\lesssim 0.4$). These observations give reason to
believe that modal cross terms could be ignored in the present case as well,
lending credibility to the feasibility of using dissipation-based velocity (or
velocity-temperature) modes for ROMs. We note that Petrov-Galerkin models
allow the use of non-orthogonal bases (Taira  _et al._ , 2020), and that data-
driven learning methods such as SINDY (Brunton and Kutz, 2022) could be used
to determine the ROM in the case where these approaches fail.

Figure 12: Non-vanishing overlap magnitudes between normalized dissipation
velocity-temperature modes with $n_{D},n_{D}^{\prime}\leq 30$, computed with
the energy inner product $\left\langle\cdot,\cdot\right\rangle_{E}$.

As seen from (31), reconstruction of any second-order quantity can be obtained
by summing the individual contributions of each mode, weighted by its
eigenvalue. Figure 11 suggests that an estimate for the full quantities could
be obtained from a 20-mode reconstruction, using global rescaling factors
corresponding to the boundary layer and the core, with a possible matching law
corresponding to the intermediate entrainment zone $0.02\lesssim z\lesssim
0.1$. As a first approximation, almost all the thermal dissipation and thermal
energy, but only $\sim$60% of the viscous dissipation and kinetic energy were
captured by the 20-mode reconstruction in the boundary layer. Meanwhile,
$\sim$50% of both thermal and kinetic energy were reconstructed in the core,
against about respectively 15% and 30% for the thermal and viscous
dissipation. A logarithmic interpolation between the core and the boundary
layer values can be used in the entrainment region for all quantities, except
for the kinetic energy reconstruction, which increased up to 75% in the
entrainment zone. We note that all results reported here have been obtained at
$\mathrm{Ra}=10^{7}$. Soucasse  _et al._ (2021) showed that the energy-
optimized basis remained stable (with possible reordering among the modes) at
least over the range $10^{6}\leq\mathrm{Ra}\leq 10^{8}$, and proposed a
$\mathrm{Ra}$-based scaling for the evolution of the dominant eigenvalues.
This paves the way to an alternative POD-based approach to modeling based on
determining the evolution of POD eigenvalues with the Rayleigh number. Any
second-order quantity could be reconstructed over a range of Rayleigh numbers,
based on a set of dissipation-optimized modes and evolution laws for the
corresponding eigenvalues as well as for the velocity and temperature global
scaling factors.

##  VII Conclusions

We have presented a formulation of POD optimized with respect to viscous and
thermal dissipation, based on replacing the usual $L^{2}$ norm in the POD
optimization problem with an $H^{1}$ norm. This was motivated by the desire to
improve the ability of modal decompositions to capture dynamically important
structures, regardless of their energy content. The dissipation decomposition
produces an orthogonal basis for the combined strain rate tensor and thermal
gradient. A corresponding non-orthogonal velocity-temperature basis was
produced using the extended method of snapshots, along with an energy-
optimized orthogonal velocity-temperature basis using conventional snapshot
POD. This allowed direct comparison between the dissipation-optimized modes
and the energy-optimized POD basis. The two decomposition procedures were
applied to an enriched data set based on a DNS of a Rayleigh-Bénard flow in a
cubic cell, with Rayleigh number $\mathrm{Ra}=10^{7}$ and Prandtl number
$\mathrm{Pr}=0.707$. The underlying DNS data set of $1000$ realizations was
enriched based on the geometric symmetries of the convection cell, resulting
in a data set of $16\,000$ realizations. Analyzing commutation relations for
the POD and symmetry operators used for the enrichment led to a
characterization of the degeneracy in the POD spectra, and the enrichment was
formally shown to lead to a separation of each POD basis into distinct
subspaces characterized by different sets of symmetries and anti-symmetries,
constraining overlaps between modes of the different decompositions as well as
the non-orthognality of dissipation modes.

The decompositions were compared in detail for the lowest 20 modes of each,
for which several features were found to be shared between the two
decompositions. The lowest 20 dissipation modes were found to each have a
direct analog among the lowest 20 energy modes, with which both symmetry
configuration and large-scale organization were shared. Much of the ordering
between modes was found to be preserved between the decompositions. Both
decompositions were dominated by four modes contributing to the large-scale
circulation (LSC), labeled $M$, $L_{x/y}$, and $D$, which were nearly
identical between the decompositions. The corner roll mode $C$, which in
previous work was found to be associated with LSC reorientations, was found to
be the largest contributor to the fluctuations of the dissipation, but not
those of the energy. This lends support to the hypothesis that a dissipation-
based decomposition is an improvement over the energy-based decomposition in
identifying dynamically significant modes.

With some exceptions dissipation modes were generally found to exhibit greater
boundary layer activity and enhanced contribution to mean kinetic and thermal
energy, viscous and thermal dissipation, and convective heat flux in the
boundary layer compared to analogous energy modes, while showing similar or
decreased contributions to mean quantities in the core of the cell. Boundary
layer structures are known to play a crucial role in the dynamics of Rayleigh-
Bénard flows, and the enhanced modal representation of such structures
therefore suggests that the overall dynamics are more efficiently captured
using the dissipation decomposition.

Differences between dissipation modes and their energy analogs tended to
increase for higher modes, and particularly the roll modes $L_{x/y}^{\dagger}$
and the vertical mode $L_{z}$ exhibited substantial dissimilarities between
the decompositions. The fact that large-scale organization was generally
preserved between the decompositions was ascribed to a coupling between
energetic and dissipative structures, caused by the modest Rayleigh and
Reynolds numbers of the flow, at $\mathrm{Ra}=10^{7}$ and $\mathrm{Re}=651$.
The conservation of large-scale structures between the decompositions shows
that any improvement in the resolution of boundary layers and other small-
scale structure does not necessarily impact the ability of the decomposition
to capture large-scale structures, suggesting instead that a wide range of
scales is captured.

Modal reconstructions of global and boundary layer kinetic and thermal energy,
viscous and thermal dissipation, and convective heat flux were carried out
using each decomposition, and the convergence of the reconstructions was
analyzed. In all cases the asymptotic convergence rates of the boundary layer
reconstructions were enhanced using the dissipation decomposition compared to
the energy decomposition. The asymptotic convergence rate of global
reconstructions was similar for all quantities using the dissipation
decomposition, while it varied considerably between energy-based
reconstructions. The similar convergence rates of different quantities using
dissipation modes support the hypothesis that a wider range of relevant scales
are captured by dissipation modes, while the improved convergence in boundary
layers indicates that dynamically important boundary layer structures are more
efficiently captured by these modes.

Some possible drawbacks of applying dissipation modes as a basis for Galerkin
methods were also discussed. While it was argued that the non-orthogonality of
dissipation-optimized velocity-temperature modes present a potential
complication, the departure from orthogonality was found to be modest in the
presently studied case, suggesting that the impact of this should be
manageable. Based on reconstructed quantities we suggested a modeling approach
based on regionally defined scaling of the flow reconstructed using the
dominant eigenmodes, combined with $\mathrm{Ra}$-based scaling of dominant
eigenvalues.

We have thus shown results indicating that the dissipation-based decomposition
is able to capture dynamically important structures in a Rayleigh-Bénard
convection flow more efficiently than conventional energy-based
decompositions. The study was limited to a single configuration characterized
by a moderate Rayleigh number, and it would be interesting to determine the
effects of increasing the Rayleigh number. In particular, insight as to
whether and how the discrepancy between the decompositions increases with the
Rayleigh number would be valuable. Even at the moderate Rayleigh number
considered here results suggest that the method poses a promising alternative
to energy-based decompositions in the context of ROMs. Naturally, any
convincing assessment of the potential of this method must revolve around the
eventual implementation and evaluation of such models, which we hope will be
attempted in future studies.

###### Acknowledgements.

CMV acknowledges the support from the European Research Council: This project
has received funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 803419). CMV and PJO acknowledge the support from the Poul Due
Jensen Foundation: Financial support from the Poul Due Jensen Foundation
(Grundfos Foundation) for this research is gratefully acknowledged. The
authors gratefully acknowledge the computational and data resources provided
on the Sophia HPC Cluster at the Technical University of Denmark, DOI:
10.57940/FAFC-6M81. This work was granted access to the HPC resources of IDRIS
under the allocation 2022-A0102B00209 attributed by GENCI (Grand Équipement
National de Calcul Intensif). This work was also performed using computational
resources from the "Mésocentre" computing center of Université Paris-Saclay,
CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS
and Région Île-de-France (http://mesocentre.centralesupelec.fr/).

## Author declarations

### Conflict of interest

The authors have no conflicts to disclose.

### Author contributions

Peder J. Olesen: Conceptualization (equal), data curation (equal), formal
analysis (equal), methodology (equal), project administration (equal),
software (equal), visualization (lead), writing – original draft preparation
(lead), writing – review and editing (lead). Laurent Soucasse:
Conceptualization (equal), data curation (equal), formal analysis (equal),
methodology (equal), project administration (equal), software (equal), writing
– original draft preparation (support), writing – review and editing
(support). Bérengère Podvin: Conceptualization (equal), formal analysis
(equal), methodology (equal), project administration (equal), supervision
(equal), writing – original draft preparation (support), writing – review and
editing (support). Clara M. Velte: Conceptualization (equal), funding
acquisition (lead), project administration (equal), supervision (equal),
writing – review and editing (support).

### Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * Shraiman and Siggia (1990) B. I. Shraiman and E. D. Siggia, “Heat transport in high-Rayleigh-number convection,” Physical Review A 42, 3650 (1990).
  * Grossmann and Lohse (2000) S. Grossmann and D. Lohse, “Scaling in thermal convection: a unifying theory,” Journal of Fluid Mechanics 407, 27–56 (2000).
  * Siggia (1994) E. D. Siggia, “High Rayleigh number convection,” Annual review of fluid mechanics 26, 137–168 (1994).
  * Chillà and Schumacher (2012) F. Chillà and J. Schumacher, “New perspectives in turbulent Rayleigh-Bénard convection,” The European Physical Journal E 35, 1–25 (2012).
  * Scheel and Schumacher (2014) J. D. Scheel and J. Schumacher, “Local boundary layer scales in turbulent Rayleigh–Bénard convection,” Journal of fluid mechanics 758, 344–373 (2014).
  * Hanjalić (2002) K. Hanjalić, “One-point closure models for buoyancy-driven turbulent flows,” Annual review of fluid mechanics 34, 321–347 (2002).
  * Grossmann and Lohse (2004) S. Grossmann and D. Lohse, “Fluctuations in turbulent Rayleigh–Bénard convection: the role of plumes,” Physics of fluids 16, 4462–4472 (2004).
  * Shang  _et al._ (2003) X.-D. Shang, X.-L. Qiu, P. Tong, and K.-Q. Xia, “Measured local heat transport in turbulent Rayleigh-Bénard convection,” Physical review letters 90, 074501 (2003).
  * Emran and Schumacher (2012) M. S. Emran and J. Schumacher, “Conditional statistics of thermal dissipation rate in turbulent Rayleigh-Bénard convection,” The European Physical Journal E 35, 1–8 (2012).
  * Vishnu, De, and Mishra (2022) V. T. Vishnu, A. K. De, and P. K. Mishra, “Statistics of thermal plumes and dissipation rates in turbulent rayleigh-benard convection in a cubic cell,” International Journal of Heat and Mass Transfer 182, 121995 (2022).
  * Shishkina and Wagner (2008) O. Shishkina and C. Wagner, “Analysis of sheet-like thermal plumes in turbulent Rayleigh–Bénard convection,” Journal of Fluid Mechanics 599, 383–404 (2008).
  * Xi, Lam, and Xia (2004) H.-D. Xi, S. Lam, and K.-Q. Xia, “From laminar plumes to organized flows: the onset of large-scale circulation in turbulent thermal convection,” Journal of Fluid Mechanics 503, 47–56 (2004).
  * Hanjalić and Launder (2022) K. Hanjalić and B. Launder,  _Modelling Turbulence in Engineering and the Environment: Rational Alternative Routes to Closure_ (Cambridge University Press, 2022).
  * Berkooz, Holmes, and Lumley (1993) G. Berkooz, P. Holmes, and J. L. Lumley, “The proper orthogonal decomposition in the analysis of turbulent flows,” Annual review of fluid mechanics 25, 539–575 (1993).
  * Berkooz (1996) G. Berkooz,  _Turbulence, coherent structures, dynamical systems and symmetry_ (Cambridge University Press, 1996).
  * Lumley (1967) J. L. Lumley, “The structure of inhomogeneous turbulent flows,” Atmospheric turbulence and radio wave propagation , 166–178 (1967).
  * Sirovich and Park (1990) L. Sirovich and H. Park, “Turbulent thermal convection in a finite domain: Part I. Theory,” Physics of Fluids A: Fluid Dynamics 2, 1649–1658 (1990).
  * Bailon-Cuba, Emran, and Schumacher (2010) J. Bailon-Cuba, M. S. Emran, and J. Schumacher, “Aspect ratio dependence of heat transfer and large-scale flow in turbulent convection,” Journal of Fluid Mechanics 655, 152–173 (2010).
  * Podvin and Sergent (2012) B. Podvin and A. Sergent, “Proper orthogonal decomposition investigation of turbulent Rayleigh-Bénard convection in a rectangular cavity,” Physics of Fluids 24 (2012).
  * Verdoold, Tummers, and Hanjalić (2009) J. Verdoold, M. J. Tummers, and K. Hanjalić, “Prime modes of fluid circulation in large-aspect-ratio turbulent Rayleigh-Bénard convection,” Physical Review E 80, 037301 (2009).
  * Podvin and Sergent (2015) B. Podvin and A. Sergent, “A large-scale investigation of wind reversal in a square Rayleigh–Bénard cell,” Journal of Fluid Mechanics 766, 172–201 (2015).
  * Soucasse  _et al._ (2019) L. Soucasse, B. Podvin, P. Rivière, and A. Soufiani, “Proper orthogonal decomposition analysis and modelling of large-scale flow reorientations in a cubic Rayleigh–Bénard cell,” Journal of Fluid Mechanics 881, 23–50 (2019).
  * Bergmann, Bruneau, and Iollo (2009) M. Bergmann, C.-H. Bruneau, and A. Iollo, “Enablers for robust POD models,” Journal of Computational Physics 228, 516–538 (2009).
  * Aubry, Lian, and Titi (1993) N. Aubry, W.-Y. Lian, and E. S. Titi, “Preserving symmetries in the proper orthogonal decomposition,” SIAM Journal on Scientific Computing 14, 483–505 (1993).
  * Sengupta and Dey (2004) T. K. Sengupta and S. Dey, “Proper orthogonal decomposition of direct numerical simulation data of by-pass transition,” Computers & structures 82, 2693–2703 (2004).
  * Lee and Dowell (2020) M. W. Lee and E. H. Dowell, “Improving the predictable accuracy of fluid galerkin reduced-order models using two POD bases,” Nonlinear Dynamics 101, 1457–1471 (2020).
  * Schiødt  _et al._ (2022) M. Schiødt, A. Hodžić, F. Evrard, M. Hausmann, B. Van Wachem, and C. M. Velte, “Characterizing Lagrangian particle dynamics in decaying homogeneous isotropic turbulence using proper orthogonal decomposition,” Physics of Fluids 34 (2022).
  * Borée (2003) J. Borée, “Extended proper orthogonal decomposition: a tool to analyse correlated events in turbulent flows,” Experiments in fluids 35, 188–192 (2003).
  * Olesen  _et al._ (2023) P. J. Olesen, A. Hodžić, S. J. Andersen, N. N. Sørensen, and C. M. Velte, “Dissipation-optimized proper orthogonal decomposition,” Physics of Fluids 35 (2023).
  * Bakewell Jr and Lumley (1967) H. P. Bakewell Jr and J. L. Lumley, “Viscous sublayer and adjacent wall region in turbulent pipe flow,” The Physics of Fluids 10, 1880–1889 (1967).
  * Herzog (1986) S. Herzog,  _The large scale structure in the near-wall region of turbulent pipe flow_ (Cornell University, 1986).
  * Aubry  _et al._ (1988) N. Aubry, P. Holmes, J. L. Lumley, and E. Stone, “The dynamics of coherent structures in the wall region of a turbulent boundary layer,” Journal of fluid Mechanics 192, 115–173 (1988).
  * Deane  _et al._ (1991) A. Deane, I. Kevrekidis, G. E. Karniadakis, and S. Orszag, “Low-dimensional models for complex geometry flows: Application to grooved channels and circular cylinders,” Physics of Fluids A: Fluid Dynamics 3, 2337–2354 (1991).
  * Ravindran (2000) S. S. Ravindran, “A reduced-order approach for optimal control of fluids using proper orthogonal decomposition,” International journal for numerical methods in fluids 34, 425–448 (2000).
  * Ly and Tran (2001) H. V. Ly and H. T. Tran, “Modeling and control of physical processes using proper orthogonal decomposition,” Mathematical and computer modelling 33, 223–236 (2001).
  * Sirovich (1987) L. Sirovich, “Turbulence and the dynamics of coherent structures. I. Coherent structures,” Quarterly of applied mathematics 45, 561–571 (1987).
  * Zhang  _et al._ (2023) Y. Zhang, A. Hodžić, F. Evrard, B. Van Wachem, and C. M. Velte, “Phase proper orthogonal decomposition of non-stationary turbulent flow,” Physics of Fluids 35, 045109 (2023).
  * Towne, Schmidt, and Colonius (2018) A. Towne, O. T. Schmidt, and T. Colonius, “Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis,” Journal of Fluid Mechanics 847, 821–867 (2018).
  * Sieber, Paschereit, and Oberleithner (2016) M. Sieber, C. O. Paschereit, and K. Oberleithner, “Spectral proper orthogonal decomposition,” Journal of Fluid Mechanics 792, 798–828 (2016).
  * Hodžić, Olesen, and Velte (2022) A. Hodžić, P. J. Olesen, and C. M. Velte, “On the discrepancies between POD and Fourier modes on aperiodic domains,” arXiv preprint arXiv:2207.02550 (2022).
  * Xin and Le Quéré (2002) S. Xin and P. Le Quéré, “An extended Chebyshev pseudo-spectral benchmark for the 8:1 differentially heated cavity,” Numerical Metholds in Fluids 40, 981–998 (2002).
  * Xin, Chergui, and Le Quéré (2008) S. Xin, J. Chergui, and P. Le Quéré, “3D spectral parallel multi-domain computing for natural convection flows,” in  _Parallel Computational Fluid Dynamics, Lecture Notes in Computational Science and Engineering book series_ , Vol. 74, edited by Springer (2008) pp. 163–171.
  * Shishkina  _et al._ (2010) O. Shishkina, R. J. Stevens, S. Grossmann, and D. Lohse, “Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution,” New Journal of Physics 12, 075022 (2010).
  * Delort-Laval  _et al._ (2022) M. Delort-Laval, L. Soucasse, P. Rivière, and A. Soufiani, “Rayleigh–Bénard convection in a cubic cell under the effects of gas radiation up to $\mathrm{Ra}=10^{9}$,” International Journal of Heat and Mass Transfer 187, 122453 (2022).
  * Puigjaner  _et al._ (2008) D. Puigjaner, J. Herrero, C. Simo, and F. Giralt, “Bifurcation analysis of steady Rayleigh–Bénard convection in a cubical cavity with conducting sidewalls,” Journal of Fluid Mechanics 598, 393–427 (2008).
  * Soucasse  _et al._ (2021) L. Soucasse, B. Podvin, P. Rivière, and A. Soufiani, “Low-order models for predicting radiative transfer effects on Rayleigh–Bénard convection in a cubic cell at different Rayleigh numbers,” Journal of Fluid Mechanics 917, A5 (2021).
  * Soucasse  _et al._ (2020) L. Soucasse, B. Podvin, P. Rivière, and A. Soufiani, “Reduced-order modelling of radiative transfer effects on Rayleigh–Bénard convection in a cubic cell,” Journal of Fluid Mechanics 898, A2 (2020).
  * Rempfer (2000) D. Rempfer, “On low-dimensional Galerkin models for fluid flow,” Theoretical and Computational Fluid Dynamics 14, 75–88 (2000).
  * Rempfer and Fasel (1994) D. Rempfer and H. F. Fasel, “Dynamics of three-dimensional coherent structures in a flat-plate boundary layer,” Journal of Fluid Mechanics 275, 257–283 (1994).
  * Taira  _et al._ (2020) K. Taira, M. S. Hemati, S. L. Brunton, Y. Sun, K. Duraisamy, S. Bagheri, S. T. Dawson, and C.-A. Yeh, “Modal analysis of fluid flows: Applications and outlook,” AIAA journal 58, 998–1022 (2020).
  * Brunton and Kutz (2022) S. L. Brunton and J. N. Kutz,  _Data-driven science and engineering: Machine learning, dynamical systems, and control_ (Cambridge University Press, 2022).