#  On the characterization of the convective heat flux in turbulent Rayleigh-
Bénard convection

B. Podvin1<EMAIL_ADDRESS>L. Soucasse1  F. Yvon2  1EM2C,
Centralesupélec, CNRS, Université Paris-Saclay 2LISN, CNRS, Université Paris-
Saclay

###### Abstract

We propose to use a clustering method, Latent Dirichlet Allocation or LDA
(Frihat et al. JFM 2021), to characterize the convective heat flux in
turbulent Rayleigh-Bénard convection. LDA provides a probabilistic
decomposition of a collection of fields into latent factors that can be used
as a generative model. The technique is applied to vertical mid-sections of a
cubic cell for Rayleigh numbers in the range $[10^{6},10^{8}]$. We show that
LDA provides a relevant representation of the convective heat flux and
displays good generative properties. The latent factors identified by the
decomposition, called motifs, consist of localized regions of fluid. The
average size of the motifs decreases with the Rayleigh number and their
spatio-temporal characteristics vary within the cell: frequent, elongated
structures are found within the vertical boundary layers, while rounder,
larger structures generally aligned with the velocity field are observed in
the entrainment zone. At the lower Rayleigh numbers, large, infrequent
structures are also present in the bulk, but spatial coherence is
significantly reduced at the highest Rayleigh number $Ra=10^{8}$. The most
often detected motifs are centered at a distance of $0.25H$ from the
horizontal plate and consist of elongated structures within the vertical
boundary layer. The prevalence of these motifs depends on the large-scale
organization of the flow. Using a simple physical argument, we propose a model
based on the characteristics of these motifs to predict the average
reorientation rate of the large-scale circulation.

###### keywords:

MSC Codes  (Optional) Please enter your MSC Codes here

##  1 Introduction

Rayleigh-Bénard convection, in which a fluid is heated from below and cooled
from above, represents an idealized configuration to study thermal convection
phenomena. These characterize a variety of applications ranging from
industrial processes such as heat exchangers to geophysical flows in the
atmosphere or the ocean. A central question is to determine how the heat
transfer depends on nondimensional parameters such as the Prandtl number
$Pr=\nu/\kappa$ where $\nu$ is the kinematic viscosity and $\kappa$ the
thermal diffusitivity, and the Rayleigh number

$Ra=\frac{g\beta\Delta TH^{3}}{\nu a},$ (1)

where $g$ is the gravity, $\beta$ is the thermal expansion coefficient,
$\Delta T$ the temperature difference and $H$ the cell dimension. The
Grossmann & Lohse (2000)’s theory constitutes a unified approach to address
this question. It is based on separating the contributions from the bulk
averaged thermal and kinetic dissipation rate into two subsets, one
corresponding to the boundary layers, and one corresponding to the bulk. This
theory was further refined in Grossmann & Lohse (2004), where the thermal
dissipation rate was split into a contribution from the plumes and a
contribution from the turbulent background. Through the action of buoyancy,
the thermal boundary layers generate plumes which create a large-scale
circulation, as evidenced by Xi et al. (2004), also called ”wind” (Castaing et
al., 1989). The distribution of temperature fluctuations depends on plume
clustering effects (Wang et al., 2022), but it is also affected by interaction
with turbulent fluctuations in the bulk, resulting in fragmentation (Bosbach
et al., 2012).

Shang et al. (2003) showed that thermal plumes carry most of the convective
heat flux and that the plume-dominated regions were located near the sidewalls
and the conducting surfaces. The morphology of plumes and its effect on the
heat transfer have been given careful attention. The plumes have a sheet-like
structure near the boundary layer and progressively become mushroom-like as
they move into the bulk region (Zhou et al., 2007). Shishkina & Wagner (2008)
found that very high values of the local heat flux were observed in regions
where the sheet-like plumes merged, constituting ”stems” for the mushroom-like
plumes developing in the bulk. The relative contributions of the plumes and
turbulent background vary with the Rayleigh number: Emran & Schumacher (2012)
have shown that the fraction of plume-dominated regions decreases with the
Rayleigh number, while that of background-dominated regions increases.

The extraction of plumes from the temperature fluctuations is therefore an
essential step for the understanding of thermal convection flows. Several
definitions have been used: some of the first criteria were based on the
skewness of the temperature derivative (Belmonte & Libchaber, 1996) or the
temperature difference (Zhou & Xia, 2010). Ching et al. (2004) have proposed
to use simultaneous measurements of the temperature and the velocity to define
the velocity of the plumes using conditional averaging. Following Huang et al.
(2013), van der Poel et al. (2015) identified plumes from both a temperature
anomaly and an excess of convective heat flux. Zhou et al. (2016) relied on
cliff-ramp-like structures in the temperature signals to determine the spatial
characteristics of plumes. Emran & Schumacher (2012) and Vishnu et al. (2022)
separated the plume from the background regions based on a threshold on the
convective heat flux. Shevkar et al. (2022) have recently proposed a dynamic
criterion based on the 2-D velocity divergence to separate plumes from
boundary layers.

As pointed out by Chilla & Schumacher (2012), this multiplicity of criteria
illustrates the difficulty of identifying coherent structures in a consistent
and objective manner, which is a long-running question in turbulent flows. To
this end, Proper Orthogonal Decomposition (POD) (Lumley, 1967) has proven a
useful tool to analyze large-scale fluctuations in Rayleigh-Bénard convection.
It has been used in particular to study reorientations of the large-scale
circulation (Bailon-Cuba et al., 2010; Foroozani et al., 2017; Podvin &
Sergent, 2015, 2017; Soucasse et al., 2019). Through spectral decomposition of
the autocorrelation tensor, POD provides a basis of spatial modes, also called
empirical modes, since they originate from the data. The modes are
energetically optimal to reconstruct the fluctuations. The POD modes typically
have a global support, which is well suited to capture the large-scale
organization of the flow. However, this can make physical interpretation
difficult as there is no straightforward connection between a mode and a local
coherent structure as a local structure is represented with a superposition of
many POD modes, a situation also observed in Fourier analysis.

As an alternative, Frihat et al. (2021) have recently adapted a probabilistic
method that can extract localized latent factors in turbulent flow
measurements. This method, Latent Dirichlet Allocation or LDA (Griffiths &
Steyvers, 2002; Blei et al., 2003), was originally developed in the context of
natural language processing, where it aims to extract topics from a collection
of documents. In this framework, documents are represented by a non-ordered
set of words taken from a fixed vocabulary. A word count matrix can be built
for the collection, where each column corresponds to a document, each line
corresponds to a vocabulary word and the matrix entry represents the number of
times the word appears in the document. LDA provides a probabilistic
decomposition of the word count matrix, based on latent factors called topics.
Topics are defined by two distributions: the distribution of topics within
each document (each document is associated with a mixture of topics, the
coefficients of the mixture sum up to one) and the distribution of vocabulary
words with each topic (each topic is represented by a combination of words,
the coefficients of which also sum up to one).

The method has been adapted for turbulent flows as follows: we consider a
collection of snapshots of a scalar field taken over a 2D domain discretized
into cells. The equivalent of a document is therefore a snapshot, and the
cells (or snapshot pixels) constitute the vocabulary. The digitized values of
the scalar field over the cells in a snapshot are gathered into a vector which
is formally analogous to a column of the word count matrix. The ”topics”
produced by the decomposition, called motifs, correspond to fixed (in the
Eulerian sense), spatially coherent regions of the flow. The method was
applied to the analysis of the turbulent Reynolds stress in wall turbulence
(Frihat et al., 2021) and to pressure anomalies (Fery et al., 2022), where it
was found to be well suited for the representation of sparse (intermittent)
data.

In this paper we examine how LDA can help characterize the spatial
organization of the convective heat flux and its dependence on the Rayleigh
number. To this end, the technique is applied to 2D sections of a cubic
Rayleigh-Bénard cell in the range of Rayleigh number $[10^{6},10^{8}]$. The
numerical configuration and the data set are described in Section 2. Proper
Orthogonal Decomposition (POD) and Latent Dirichlet Allocation (LDA) are
respectively presented in Section 3 and 4. We examine in Section 5 how LDA
compares with POD and the extent to which it is able to capture the general
features of the heat flux. The spatial characteristics of the motifs are
detailed in Section 6. We then show in Section 7 how the signature of the
large-scale organization of the flow can be tracked in the dominant motifs.
Our results are summarized in Section 8.

##  2 Numerical setting

###  2.1 Set-up

Numerical setup and associated datasets are the same as used in Soucasse et
al. (2019, 2021). 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 various
values of the Rayleigh number. The Prandtl number is set to 0.707. All
physical quantities are made dimensionless using the cell size $H$, the
reference time $H^{2}/(a\sqrt{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.

Navier–Stokes equations under Boussinesq approximation are solved using a
Chebyshev collocation method (Xin & Le Quéré, 2002; Xin et al., 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. Numerical parameters are given in Table 1 for the four considered
Rayleigh numbers $Ra=\\{10^{6};3~{}10^{6};10^{7};10^{8}\\}$. 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). A number of 1000 snapshots have been extracted from the simulations
for each Rayleigh number at a sampling period of 10 (at
$Ra=\\{10^{6};3~{}10^{6};10^{7}\\}$) or 5 (at $Ra=10^{8}$), in dimensionless
time units. It is worth noting that the time separation between the snapshots
is sufficient to describe the evolution of the large-scale circulation but is
not suited for a fine description of the plume emission or of the
reorientation process. For each Rayleigh number, a dataset satisfying the
statistical symmetries of the flow was then constructed from these 1000
snapshots, as it will be described in the next section.

$Ra$ | $(N_{x},N_{y},N_{z})$ | $N_{S}$ | $\Delta t$ | $\delta_{BL}$
---|---|---|---|---
$10^{6}$ | (81,81,81) | 1000 | 10 | 0.056
$3$ $10^{6}$ | (81,81,81) | 1000 | 10 | 0.042
$10^{7}$ | (81,81,81) | 1000 | 10 | 0.0297
$10^{8}$ | (161,161,161) | 1000 | 5 | 0.0167
Table 1: Characteristics of the datasets at various Rayleigh numbers: spatial
resolution $N_{x}$, $N_{y}$, $N_{z}$ in each direction of space, number of
snapshot $N_{S}$, snapshot sampling period $\Delta t$ and boundary layer
thickness $\delta_{BL}$.

###  2.2 Construction of the data set

At each Rayleigh number, the data set consisted of a collection of
$N_{S}=1000$ heat flux snapshots
$w(\underline{x},t_{k})\theta(\underline{x},t_{k})$, $k=1,\ldots N_{S}$. We
note that due to the velocity reference scale, the non-dimensional heat flux
varies like $NuRa^{-1/2}$. As in Soucasse et al. (2019), the data set was
first 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 for this Rayleigh number range: the large-
scale circulation settles in one of the two diagonal planes of the cube with
clockwise or counterclockwise motion. The evolution of the large-scale
circulation can be tracked through that of the $x$ and $y$ components of the
angular momentum of the cell
$\underline{L}=\int(\underline{x}-\underline{x_{0}})\times\underline{u}d\underline{x}$
with respect to the cell center $\underline{x_{0}}$. As Figure 1 shows at
$Ra=10^{7}$, the angular momentum along each horizontal direction oscillates
near a quasi-steady position for long periods of times - several hundreds of
convective time scales, before experiencing a rapid switch ($\mathcal{O}(10)$
convective time scales) to the opposite value, which corresponds to a
reorientation. Reorientations from one state to another occur during the time
sequence but each state is not necessarily equally visited. In order to
counteract this bias, we have built enlarged snapshot sets, obtained by the
action of the symmetry group of the problem on the original snapshot sets. The
symmetries are based on four independent symmetries $S_{x}$, $S_{y}$, $S_{z}$
and $S_{d}$ with respect to the planes $x=0.5H$, $y=0.5H$, $z=0.5H$ and $x=y$.
This generates a group of 16 symmetries for the cube, which should lead to a
16-fold in the number of snapshots. However, since we will exclusively
consider the vertical mid-planes $x=0.5H$ and $y=0.5H$, which are invariant
planes for respectively $S_{x}$ and $S_{y}$, the increase is reduced. The data
set aggregates 1000 snapshots on each of the planes $x=0.5H$ and $y=0.5H$,
each of which undergoes a vertical flip, a horizontal flip and a combination
of the two, yielding a total of $N_{s}=8000$ snapshots.

The LDA technique requires transforming the data into a non-negative, integer
field. The signal defined on a grid of $\tilde{N_{C}}$ cells was digitized
using a rescaling factor $s$ and positive and negative values were split onto
two distinct grids, leading to a field defined on $N_{C}=2\tilde{N}_{C}$
cells. We thus define

$\displaystyle\Phi_{k}(\underline{x}_{j})=\Phi(\underline{x}_{j},t_{k})$
$\displaystyle=$
$\displaystyle\mbox{Max}\left[\mbox{Int}[s\,w(\underline{x}_{j},t_{k})\theta(\underline{x}_{j},t_{k}),0\right]$
(2)
$\displaystyle\Phi_{k}(\underline{x}_{j+\tilde{N}_{C}})=\Phi(\underline{x}_{j+\tilde{N}_{C}},t_{k})$
$\displaystyle=$
$\displaystyle-\mbox{Max}\left[-\mbox{Int}[s\,w(\underline{x}_{j},t_{k})\theta(\underline{x}_{j},t_{k}),0\right]$
(3)

where $s>0$, $k\in[1,N_{S}]$ and $j\in[1,\tilde{N}_{C}]$ and
$\underline{x}_{j}$ represents the $j^{\mathrm{th}}$ cell location on the
plane $x=0.5H$ or $y=0.5H$. We note that throughout the paper, the total field
will directly be represented on the physical grid of size $\tilde{N}$ from the
renormalized difference
$[\Phi(\underline{x}_{j},t_{k})-\Phi(\underline{x}_{j+\tilde{N}_{C}},t_{k})]/s$.

##  3 POD analysis

###  3.1 Method

Proper Orthogonal Decomposition (Holmes et al., 2002) makes it possible to
write a collection of $N_{S}$ spatial fields $\Phi(\underline{x}_{j},t_{k})$
defined on $N_{C}$ grid points, as a superposition of spatial modes
$\phi_{n}(\underline{x})$, the amplitude of which varies in time:

$\Phi(\underline{x}_{j},t_{k})=\sum_{n=1}^{N_{S}}a_{n}(t_{k})\phi_{n}(\underline{x}_{j}),$
(4)

with $k\in[1,N_{S}]$ and $j\in[1,N_{C}]$. The amplitudes $a_{n}(t_{k})$ are
solution of the eigenvalue problem

$C_{kp}a_{n}(t_{p})=\lambda_{n}a_{n}(t_{k}),$ (5)

where $C$ is the temporal autocorrelation matrix

$C_{kp}=\frac{1}{N_{S}}\sum_{j=1}^{N_{C}}\Phi(\underline{x}_{j},t_{k})\Phi(\underline{x}_{j},t_{p})$
(6)

The eigenvalues $\lambda_{n}$, such that
$\lambda_{1}>\lambda_{2}>\lambda_{3}>\ldots$, represent the respective
contribution of the modes to the total variance. If we consider the $p$ most
energetic modes, the reconstruction based on $p$ modes minimizes the
$L_{2}$-norm error between the set of snapshots and the projection of the set
of snapshots onto a basis of size $p$.

###  3.2 Application to the convective heat flux

POD is applied to the digitized signal $\Phi$ defined in equations (2) and
(3). The first three POD modes and POD coefficients are shown in Figure 2 for
$Ra=10^{7}$, where black vertical and horizontal lines indicate the thickness
of the boundary layers. We checked that the first mode corresponds to the mean
flow. The mode is most important in a region close to the wall, with a maximum
within the vertical boundary layer at a height of about $z\sim 0.1H$. The
second mode corresponds to a dissymetry between the vertical sides and is most
important at mid-height in the region outside the boundary layers. The third
mode is both antisymmetric in the vertical and in the horizontal direction. It
is maximum at the edge of the vertical boundary layers, at a vertical distance
of about $0.25H$ from the horizontal surfaces. The pattern it is associated
with corresponds to a more intense flux along a diagonal (bottom of one side
and top of the opposite side) and a less intense flux along the opposite
diagonal. As evidenced by application of a moving average performed over 200
convective time units (about 4 recirculation times $T_{c}$, as was determined
in Soucasse et al. (2019)), the evolution of the amplitude at large time
scales matches that of the horizontal angular momentum components $L_{x}$ and
$L_{y}$ (compare with Figure 1), unlike the two dominant modes. This mode
therefore appears to be the signature of the large-scale circulation, where
the flux is more intense in the lower corner of the cell as hot plumes rise on
one side and in the upper corner of the opposite side of the cell as cold
plumes go down.

##  4 Latent Dirichlet Allocation

###  4.1 Principles

We briefly review the principles of Latent Dirichlet Allocation and refer the
reader to Frihat et al. (2021) for more details. LDA is an inference approach
to identify latent factors in a collection of observed data, which relies on
Dirichlet distributions as priors. We first recall the definition of a
Dirichlet distribution, which is a multivariate probability distribution over
the space of multinomial distributions. It is parameterized by a vector of
positive-valued parameters
$\underline{\alpha}=(\alpha_{1},\alpha_{2},\ldots,\alpha_{N})$ as follows

$p(\theta_{1},...,\theta_{N};\alpha_{1},...,\alpha_{N})=\frac{1}{B(\underline{\alpha})}\prod_{n=1}^{N}\theta_{n}^{\alpha_{n}-1},$
(7)

where $B$ is a normalizing factor, which can be expressed in terms of the
Gamma function $\Gamma$:

$B(\underline{\alpha})=\frac{\prod_{n=1}^{N}\Gamma(\alpha_{n})}{\Gamma(\sum_{n=1}^{N}\alpha_{n})}.$
(8)

The components $\\{\alpha_{n},n=1\dots N\\}$ of $\underline{\alpha}$ control
the sparsity of the distribution: values of $\alpha_{n}$ larger than unity
correspond to evenly dense distributions, while values lower than unity
correspond to sparse distributions. Here $\theta$ will represent either the
motif-cell distribution or the snapshot-motif distribution.

As mentioned above, the data to which LDA is applied consists of a collection
of non-negative, integer fields that are defined in equations (2) and (3). For
each snapshot $k$, the integer value $\Phi_{k}(\underline{x}_{j})$ measured at
cell $j$ is interpreted as an integer count of the cell $j$. The key is to
interpret this integer count as the number of times cell $j$ appears in the
description of snapshot $k$. A snapshot $k$ is therefore defined as a list of
tuples of the form $(\underline{x}_{j},\Phi_{k}(\underline{x}_{j}))$.

The main assumptions of LDA are the following:

  1. 1.

Each snapshot consists of a mixture of $N_{T}$ latent factors called motifs.
$N_{T}$ is a user-defined parameter (analogous to a number of clusters).

  2. 2.

Each motif $n$ is associated with a multinomial distribution over the grid
cells $\psi_{n}$ so that the probability to observe the $j^{\mathrm{th}}$ grid
cell located at $\underline{x}_{j}$ given the motif $n$ is $\psi_{n}(x_{j})$.
The distribution $\psi_{n}$ is modelled with a Dirichlet prior parameterized
with an $N_{C}$-dimensional vector $\underline{\eta}$. Low values of
$\eta_{l}$ mean that the motif is distributed over a small number of cells.

  3. 3.

Each snapshot $k$ is associated with a distribution $b^{k}$ over motifs such
that the probability that motif $n$ is present in snapshot $k$ is $b_{n}^{k}$,
which we will rewrite as $b_{n}^{k}=b_{n}(t_{k})$. This distribution is
modelled with a $N_{T}$-dimensional Dirichlet distribution of parameter
$\underline{\alpha}$. The magnitude of $\alpha$ characterizes the sparsity of
the distribution. Low values of $\alpha_{l}$ mean that relatively few motifs
are observed in the snapshot.

###  4.2 Implementation

The snapshot–motif distribution $b^{k}$ and the motif-cell distribution
$\psi_{n}$ are determined from the observed snapshots
$\Phi_{k}(\underline{x})$ and constitute $N_{T}$ and $N_{C}$-dimensional
categorical distributions. Finding the distributions $b^{k}$ and $\psi_{n}$
that are most compatible with the observations constitutes an inference
problem. The problem can be solved either with a Markov chain Monte-Carlo
(known as MCMC) algorithm such as Gibbs sampling (Griffiths & Steyvers, 2002),
or by a variational approach (Blei et al., 2003), which aims to minimize the
Kullback–Leibler divergence between the true posterior and its variational
approximation. In both cases, the computational complexity of the problem of
the order of $(N_{C}N_{S}N_{T})$.

The solution a priori depends on the number of motifs $N_{T}$ as well as on
the values of the Dirichlet parameters $\underline{\alpha}$ and
$\underline{\eta}$. Special attention was therefore given to establish the
robustness of the results reported here. Non-informative default values were
used for the Dirichlet parameters i.e. the prior distributions were taken with
symmetric parameters equal to $\forall n,\alpha_{n}=1/N_{T}$ and $\forall
l,\eta_{l}=1/N_{C}$. Practical implementation was performed in Python using
gensim (Rehurek & Sojka, 2011). No significant change was observed in the
results when the value of the quantization $s$ was high enough (however it had
to be kept reasonably low in order to limit the computational time). Although
multiple tests were carried out for varying values of $s\in[40,600]$, all
results reported in this paper were obtained with $s=600$. Analyses were also
performed for varying numbers of motifs $N_{T}$, ranging from 50 to 400.

###  4.3 LDA as a generative process

The standard generative process performed by LDA with $N_{T}$ motifs is the
following.

  * $(i)$

For each motif $n$, a $N_{C}$-dimensional cell–motif distribution $\psi_{n}$
is drawn from the Dirichlet distribution of parameter $\underline{\eta}$.

  * $(ii)$

To generate snapshot $k$:

    * $(a)$

a $N_{T}$-dimensional snapshot–motif distribution $b^{k}$ is drawn according
to a Dirichlet distribution parameterized by $\underline{\alpha}$

    * $(b)$

a total integer count $\Phi_{T}(t_{k})=\sum_{j}\Phi(\underline{x}_{j},t_{k})$
is drawn This number corresponds to the total number of cell integer counts
associated with snapshot $k$. $\Phi_{T}(t_{k})$ is typically sampled from a
Poisson distribution that matches the statistics of the original database.

    * $(c)$

for each $i=1,\ldots,\Phi_{T}(t_{k})$:

      * *

a motif $n$ is selected from $b(t_{k})$ (since $b$ represents the probability
that motif $n$ is present in the $k^{\mathrm{th}}$ snapshot)

      * *

once this motif $n$ is chosen, a cell $j$ is selected from $\psi_{n}$ (since
$\psi_{n}(\underline{x}_{j})$ represents the probability that cell $j$ is
present in motif $n$)

The snapshot $k$ then represents the set of $\Phi_{T}$ cells $j$ that have
been drawn and can be reorganized as a list of $N_{C}$ cells with integer
counts $\Phi(\underline{x}_{j},t_{k})$.

In fluid mechanics applications ((Frihat et al., 2021; Fery et al., 2022)),
sampling from the motif-cell distribution (step c)) can be replaced with a
faster step, where the contribution of each motif $n$ to snapshot $k$ is
directly obtained from the motif-cell distribution $\psi_{n}$ and the
distribution $b(t_{k})$ as
$\Phi_{T}(t_{k})b_{n}(t_{k})\psi_{n}(\underline{x}_{j})$. Figure 3 illustrates
the LDA generative process on a $4\times 3$ grid for three topics.

###  4.4 Interpretation and evaluation criteria

By construction, the decomposition identifies fixed regions of space over
which the intensity of the scalar field is important at the same time. There
is therefore no direct connection between the motifs and plumes, which are
Lagrangian structures travelling and possibly changing in shape and
orientation through the shell. LDA motifs only aim to detect the Eulerian
signature of these structures.

Each motif $n$ can be characterized in space through the motif-cell
distribution $\psi_{n}$ (which integrates to 1 over the cells) and which will
sometimes referred to as the motif in the absence of ambiguity. Each
distribution has a maximum value $\psi^{max}$ and a maximum location
$\underline{x}_{n}^{max}$ such that
$\psi_{n}(\underline{x}^{max})=\psi^{max}$. One can also define a
characteristic area $\Sigma_{n}$ using

$\Sigma_{n}=\int_{\Omega}\mathbbm{1}_{\\{\psi_{n}\geq\psi_{n}^{max}/e\\}}d\Omega$
(9)

where $\Omega$ represents the plane of analysis and the factor $1/e\sim 0.606$
is an arbitrary factor chosen by analogy with a Gaussian distribution. Each
motif can also be characterized in time through the snapshot-motif
distribution $b_{n}^{k}$, which we will also call prevalence throughout the
paper. The motifs can be ordered by their time-averaged prevalence $\langle
b_{n}\rangle=\frac{1}{N_{S}}\sum_{k=1}^{N_{S}}b_{n}(t_{k})$ where
$\langle\cdot\rangle$ represents a time average.

A reconstruction of the field can be obtained by using the inferred motif-cell
distribution and snapshot-motif distribution to provide what we will call the
LDA-Reconstructed field, defined as

$\Phi_{R}(\underline{x}_{j},t_{k})=\Phi_{T}\sum_{n=1}^{N_{T}}b_{n}(t_{k})\psi_{n}(\underline{x}_{j})$
(10)

where $\Phi_{T}$ represents the sum of the field values over the cells. To
evaluate the relevance of the decomposition, one can compute for each snapshot
$k$ the instantaneous correlation coefficient $C_{k}$ between a given field
$\Phi$ and its reconstruction $\Phi_{R}$ defined as

$C_{k}(\Phi,\Phi_{R})=\frac{\int\Phi(\underline{x},t_{k})\Phi_{R}(\underline{x},t_{k})d\underline{x}}{\left(\int\Phi^{2}(\underline{x},t_{k})d\underline{x}\int\Phi_{R}^{2}(\underline{x},t_{k})d\underline{x}\right)^{1/2}}$
(11)

A global measure of the reconstruction is then given by $\langle
C\rangle=\frac{1}{N_{S}}\sum_{k=1}^{N_{S}}C_{k}$, the average value of $C$
over all snapshots.

##  5  Evaluation of LDA for reconstruction and generation

###  5.1 Reconstruction

We first evaluate to which extent the LDA decomposition provides an adequate
reconstruction of the heat flux. Figure 4 (left) shows how the instantaneous
value of the correlation coefficient $C_{k}(\Phi,\Phi_{R})$ depends on the
integral of the field $I_{k}=\int\Phi(\underline{x},t_{k})d\underline{x}$. The
Rayleigh considered is $Ra=10^{7}$ and the number of topics is $N_{T}=100$,
but the same trend was observed for all other Rayleigh numbers as well as all
other values of $N_{T}$. Lower values of the correlation were associated with
lower values of the total integrated heat flux, which illustrates that the LDA
representation is suited to capture extreme events.

Figure 4 (right) presents the mean correlation coefficient $\langle C\rangle$
for different number of motifs and different Rayleigh numbers on the vertical
planes. The minimum value for the lower number of topics and the highest
Rayleigh number was 0.8, which shows the relevance of the decomposition.
Unsurprisingly, the correlation increases with the number of topics. It also
decreases with the Rayleigh number, which is consistent with an increase in
the complexity of the flow.

Figure 5 compares an original snapshot at $Ra=10^{7}$ (based on the digitized
signal) with different reconstructions: i) the LDA-reconstruction based on
$N_{T}=100$ motifs, ii) the reconstruction limited to the 20 most prevalent
topics (for this particular snapshot), iii) the POD-based reconstruction based
on the first 20 modes. By construction, POD provides the best approximation of
the field for a given number of modes. Since the distribution of the heat flux
is intermittent in space and time, only a limited number of motifs is
necesssary to reconstruct the flow. We note that little difference was
observed between the full LDA reconstruction and the reconstruction limited to
the 20 most prevalent motifs, which highlights the intermittent nature of the
field. The relative error between the original and the reconstructed field is
29% for the full LDA reconstruction, 34 % when the 20 most prevalent modes are
retained in the reconstruction. In contrast, limiting the POD to 20 global
modes slightly lowers the quality of the reconstruction, with a global error
of 38%. It should be noted that the 20 dominant POD modes correspond to an
average over all snapshots, while the 20 most prevalent LDA modes are selected
for that specific snapshot. On average, the reconstructed field based on
keeping the 20 most prevalent motifs differed by less than 10% from the full
100-mode reconstruction and the average correlation coefficient $C$ decreased
from 0.89 to 0.83. This shows that LDA can provide a compact representation of
the local heat flux.

###  5.2 Generation

The ability to generate statistically relevant synthetic fields is of interest
for a number of applications, such as accelerating computations or developing
multi-physics models. As a generative model, LDA makes it possible to produce
such a set of fields, the statistics of which can be compared with those of
the original fields used to extract the motifs, as well as with those of the
corresponding LDA-Reconstructed fields. It would also be useful to compare the
generated LDA data set with one generated using POD. To this end, we generated
two sets of 4000 new fields using both LDA and POD. The same number
$N_{T}=100$ of POD modes and LDA motifs to generate the datasets. The plane in
which the data is generated is assumed to be the $(y,z)$ plane. The different
fields to be compared are therefore the following:

  1. 1.

the original (digitized) field $\Phi$ defined in section 3.2 with equations
(2) and (3)

  2. 2.

the LDA-Reconstructed field (LDA-R) as defined in equation (10)

  3. 3.

the LDA-Generated field (LDA-G): the field is constructed by sampling
prevalences $\tilde{b}_{n}(t_{k})$ from snapshot-motif distributions and then
reconstructing

$\Phi^{LDA-G}(\underline{x}_{j},t_{k})=\Phi_{T}\sum_{n=1}^{N_{T}}\tilde{b}_{n}(t_{k})\psi_{n}(\underline{x}_{j}),$
(12)

where $\Phi_{T}$ is a random variable obtained by sampling a Poisson
distribution with the same statistics as the original database.

  4. 4.

the POD-Generated field (POD-G): the field is constructed by independently
sampling $N_{T}$ POD mode amplitudes $\tilde{a}_{n}$ from the POD amplitudes
of the original database

$\Phi^{POD-G}(\underline{x}_{j},t_{k})=\sum_{n=1}^{N_{T}}\tilde{a}_{n}(t_{k})\phi_{n}(\underline{x}_{j})$
(13)

The time-averaged fields corresponding to the different databases are compared
in Figure 6. A good agreement is observed for all datasets, with global errors
of 4%, 8% and 3% for respectively the LDA-reconstructed, the LDA-generated and
the POD-generated datasets. Although it provides the lowest error (as could be
expected), the POD-generated data set overestimates negative values in the
core of the cell.

For a given location $(y_{0},z_{0})$, we defined spatial autocorrelation
functions in the horizontal and vertical directions as:

$R_{y}(y,y_{0},z_{0})=\frac{\langle\Phi(y,z_{0},t)\Phi(y_{0},z_{0},t)\rangle}{\langle\Phi(y_{0},z_{0},t)^{2}\rangle}$
(14)
$R_{z}(z,y_{0},z_{0})=\frac{\langle\Phi(y_{0},z,t)\Phi(y_{0},z_{0},t)\rangle}{\langle\Phi(y_{0},z_{0},t)^{2}\rangle}$
(15)

The autocorrelation functions are displayed in Figure 7 for the selected
locations indicated in Figure 6, which correspond to regions of high heat
flux. We can see that that in all cases, the flux remains correlated over much
longer vertical extents than in the horizontal direction. Both the LDA-
reconstructed and the POD-generated autocorrelations approximate the original
data well \- again, by construction, POD-based fields are optimal to
reconstruct second-order statistics. The LDA-generated autocorrelation is not
as close to the original one, but still manages to capture the characteristic
spatial scale over which the fields are correlated.

One-point pdfs of the convective flux are represented in Figure 8 for the same
selected locations (again, indicated in Figure 6). POD-generated fields tend
to overpredict lower values of the convective flux while underpredicting
higher values and therefore do not capture well the intermittent features of
the heat flux. The LDA-generated fields display a better agreement with the
original fields and are in particular able to capture the exponential tails of
the distributions.

##  6 Spatial organization of the motifs

We now describe the spatial organization of the motifs through the motif-cell
distribution $\psi_{n}$. The general trends reported below held for all values
of $N_{T}$ considered, which ranged from 50 to 400. For all Rayleigh numbers,
most LDA motifs were found to be associated with a positive flux (i.e were
associated with the first $\tilde{N}$ cells in the decomposition). Negative
(counter-gradient) motifs were also identified, but their average prevalence
was generally very small (at most 10 % of that of the dominant motif). We
therefore chose to focus only on the motifs with a positive contribution to
the heat flux. Figure 9 (left) displays these motifs for three different
Rayleigh numbers for $N_{T}=100$. The motif-cell distribution is materialized
by a black line corresponding to the iso-probability contour of
$0.606\psi_{n}^{max}$, which can be compared with the average value of the
heat flux at this location. For all Rayleigh numbers, the motifs are clustered
in the regions of high heat flux, close to the vertical walls. Within the
vertical boundary layers, motifs are elongated in shape. Outside the vertical
boundary layers, the motifs are more isotropic and tend to increase in shape
as one moves away from the walls. Outside the horizontal boundary layers, the
motif-cell distributions are elongated in the direction of the wind, with a
horizontal orientation in the center of the cell, and a gradual vertical shift
closer to the walls. Large motifs are found in the bulk at $Ra=10^{6}$ and
$Ra=10^{7}$ (it was also the case at $Ra=3$ $10^{6}$). In contrast, fewer,
smaller motifs are found in the bulk at $Ra=10^{8}$ in the central region
$y\in[0.2,0.8]H$, signalling a loss of spatial coherence in the bulk at this
Rayleigh number.

In general, the motif size seems to decrease with the Rayleigh number. This is
confirmed by Figure 10, which represents the average motif area as a function
of their distance from the vertical walls. In order to avoid the influence of
the horizontal plates, we only considered the motifs located at a vertical
distance larger than $0.07H$ from the horizontal walls (i.e. outside the
horizontal entrainment region). The size of the symbols shown in the picture
is proportional to the fraction of motifs over which the average was
performed. Results were relatively robust with respect to the number of topics
$N_{T}$, although some dependence on $N_{T}$ is observed in the center of the
cell. Within the boundary layer, the motif area grows quadratically, which
means that the characteristic size of the motif is proportional to the wall
distance. We note that a similar scaling was found for turbulent eddies in
pressure-gradient driven turbulence such as channel flow (Frihat et al.,
2021). Further away from the vertical wall, after a short plateau at the edge
of the boundary layer, a slower increase in the motif size was observed with a
rate that increased with the Rayleigh number, so that the motif area was about
the same (on the order of $0.02H^{2}$) for all Rayleigh numbers in the center
of the cell. This suggests the presence of a double scaling for the motifs:
one based on the boundary layer thickness, and one based on the cell size. The
decrease in size with the Rayleigh number appears consistent with a dependence
on the boundary layer thickness but also with an increase of the fragmentation
by the bulk turbulent fluctuations, in agreement with the literature (Bosbach
et al., 2012; van der Poel et al., 2015). The difference observed at the
highest Rayleigh number also signals that the flow is still evolving and has
not reached an asymptotic state.

##  7 Dominant motifs

###  7.1 Description

Owing to the symmetry of the database, the motifs in the vertical plane
$(x,z)$ (resp. $(y,z)$) should approximate the symmetry $x\rightarrow 1-x$
(resp. $y\rightarrow 1-y$), and $z\rightarrow 1-z$ (complete symmetry cannot
be expected owing to the stochastic nature of the decomposition). The four
dominant motifs at $Ra=10^{7}$ are represented in Figure 11 (left) for
$N_{T}=100$. They consist of elongated structures lying mostly in the boundary
layer, and located at a vertical distance of about 0.25 H from the horizontal
walls. They can be interpreted as the wall imprint of hot plumes rising in the
boundary layer (resp. cold plumes descending in the boundary layer). Although
the positions and sizes of the four motifs may slightly vary from one to the
other, their features are generally similar and a characteristic motif can be
obtained from taking the average over all four motifs. Figure 11 (right)
represents this characteristic motif for the various Rayleigh numbers. We can
see that the dominant motifs are always located mostly within the boundary
layer, with a maximum at a height of about $0.25$ $H$. Their characteristic
width $l_{y}$ was found to decrease as $Ra^{-0.23\pm 0.04}$, which matches the
scaling of the boundary layer thickness.

The evolution of the snapshot-motif distribution, or motif prevalence, is
represented in Figure 12 for $Ra=10^{7}$. We can see that the behavior of the
motif prevalence depends on the sign of the global momentum represented in
Figure 1. When a moving average of $T_{m}=150$ time units, corresponding to 3
recirculation times, was applied, two quasi-stationary states could be
identified in each plane (and are materialized by the dashed horizontal black
lines indicated in Figure 12). The two states appear to correspond to the two
possible directions of the angular momentum component. Streamlines of the flow
conditionally averaged on the higher prevalence value of $b_{1}$ are
represented in figure 13 left). They indicate that for the higher
characteristic value of the prevalence, which we will denote $b_{+}$, the
motif is associated with the large-scale circulation while it is associated
with the corner vortex on the opposite side for the lower prevalence value,
denoted as $b_{-}$, as summarized in figure 13 right).

This indicates that information about the large-scale reorientation can be
extracted from local measurements. Two states, $I_{+}$ and $I_{-}$,
respectively corresponding to the large-scale circulation and corner vortex
can be defined from the prevalence of the dominant motif $b_{1}$ using

$I_{+}=\\{k|\langle b_{1}(t_{k})\rangle_{T_{m}}>\langle b_{1}\rangle\\}\mbox{
and }I_{-}=\\{k|\langle b_{1}(t_{k})\rangle_{T_{m}}<\langle b_{1}\rangle\\}$
(16)

where $\langle\cdot\rangle_{T_{m}}$ represents the moving average over
$T_{m}$. The average prevalences conditioned on $I_{+}$ and $I_{-}$ are
respectively $b_{+}$ and $b_{-}$.

Figure 14 displays the histogram of the prevalence of the dominant motif
$b_{1}$ (motifs 2 to 4 displayed similar features). At all Rayleigh numbers,
the total distribution is characterized by two distinct lobes, which
correspond to the absence and the presence of the motif in the snapshot. The
relative importance of the lobes therefore provides an indirect measure of the
motif intermittency, which can be related to plume emission. The ratio of
motif presence to motif absence was about 0.5-0.6 in the range of Rayleigh
numbers considered \- no significant variation was observed with the Rayleigh
number.

However, further insights can be obtained by examining the respective
contributions of the $I_{+}$ and $I_{-}$ states to the distribution of
$b_{1}$, which are also represented in Figure 14. For all Rayleigh numbers,
$I_{+}$ states contribute more to the higher-value lobe than $I_{-}$ states,
while $I_{-}$ contributes more to the lower-value lobe. This suggests that
plumes are emitted at a lower frequency in the corner rolls than in the large-
scale circulation. Moreover, the relative contributions of the $I_{+}$ and the
$I_{-}$ states vary non-monotonically with the Rayleigh number. In the higher-
value lobe, the relative contribution of $I_{-}$ appears to increase
relatively to $I_{+}$ with more high values of $I_{-}$ at $Ra=3$ $10^{6}$,
while $I_{-}$ represents more low values at $Ra=10^{8}$. In the lower-value
lobe, the contribution of $I_{+}$ is least at $Ra=3$ $10^{6}$ and largest at
$Ra=10^{8}$. These observations suggest that both the intensity of the large-
scale circulation and that of the corner roll appear to change with the
Rayleigh number, in agreement with the findings of Vishnu et al. (2020).

###  7.2 A model for the reorientation time scale

A simple model can be made to link these observations with the dynamics of
reorientations. We can assume that the conditionally averaged prevalence of
the dominant motif in the region close to the wall $b_{\pm}$ is directly
proportional to the emission rate of plumes, which can be modelled as a
Poisson point process. This means that the time separating two plume ejections
$T_{\pm}$ follows an exponential distribution with mean $1/b_{\pm}$, where $+$
and $-$ respectively characterize the large-scale circulation ($I_{+}$) and
the corner vortex ($I_{-}$) states. $b_{\pm}$ therefore represents the
parameter of the exponential distribution. A reorientation can be associated
with the event where the corner vortex becomes stronger than the large-scale
circulation state, i.e. the time separating two emissions in the corner vortex
state becomes smaller than that separating two emissions in the large-scale
circulation state. This event can occur independently in either one of the two
horizontal directions $x$ or $y$.

One can show that the probability $p$ that this event occurs at any given time
is given by

$p=p(T_{-}>T_{+})=\frac{b_{-}}{b_{+}+b_{-}}$ (17)

Owing to the memoryless nature of the exponential distribution, this holds for
the time separating an arbitrary number of emissions, in particular over a
characteristic time $T_{s}$ sufficiently long to reverse the circulation in
that direction. $T_{s}$ should be on the order of the recirculation time
$T_{c}$ so that we have $T_{s}=\beta T_{c}$ with $\beta\lessapprox 1$. If
$f_{c}$ is the recirculation frequency, one would then expect the frequency
between reorientations $f_{r}$ to depend on $p$ and $f_{c}$ following

$f_{r}=2p\beta^{-1}f_{c},$ (18)

where the factor $2$ comes from the fact that a reorientation can occur in
each direction. Figure 15 (right) compares for different Rayleigh numbers the
probability $p$ with the ratio of the frequency between reorientations and the
recirculation frequency estimated in (Soucasse et al., 2021). We see that a
very good agreement is obtained between the variations of the average
reorientation rate and the measure of the relative intensity of the large-
scale circulation and corner vortices. We note that the largest discrepancy is
observed for the highest Rayleigh number, for which the reorientation rate is
very low and therefore cannot be determined with good precision from the DNS.
The value of $\beta$ used in the figure was determined empirically and was
found to be $1/5.6\sim 0.18$. $T_{s}=\beta T_{c}$ is therefore close to the
characteristic time it takes for a fluid particle to travel up the vertical
side of the cell, which can be estimated as $T_{c}/(2(1+\sqrt{2}))\sim
0.21T_{c}$ (see Soucasse et al. (2019)). This suggests that an estimate for
the reorientation rate can be obtained by comparing directly the average
prevalence of the motif associated with the large-scale circulation with that
of its counterpart in the corner structure. This could be of particular
interest in cases where the observation time is smaller than the expected
reorientation time, a situation that is often encountered in - but not limited
to - numerical simulations.

##  8 Conclusion

We have applied a new analysis technique, Latent Dirichlet Allocation, to
characterize the spatio-temporal organization of the local heat flux. As a
soft clustering method, LDA extracts probabilistic latent factors in the form
of spatially localized motifs, from a collection of instantaneous fields. The
technique was applied to the vertical mid-plane of a Rayleigh-Bénard cubic
cell in a range of Rayleigh numbers in $[10^{6},10^{8}]$. The method was found
to be relatively robust with respect to the user-defined parameters. It was
able to provide good reconstructions of the original database and to generate
new data sets that captured key statistics of the snapshots.

The motif size and prevalence were found to vary across the cell. Different
regions could be identified: the vertical boundary layers with frequent,
elongated motifs, an entrainment region, where motifs had a rounder shape, but
were generally oriented in the direction of the large-scale circulation, and
the bulk region with infrequent motifs. A lack of spatial coherence was
observed in the bulk at the highest Rayleigh number. The size of the motifs
generally decreased with the Rayleigh number.

For all Rayleigh numbers, the dominant motifs consisted of elongated vertical
structures located mostly within the vertical boundary layer, at a height of
$0.25H$. The width of these motifs scaled with the boundary layer thickness.
The prevalence of motifs was found to depend on the large-scale organization
of the flow: two states could be identified - one corresponding to the large-
scale circulation and one to a corner roll structure. The two states were
characterized by different average prevalences which varied non-monotonically
with the Rayleigh number. A simple model was able to relate the prevalences of
the dominant motif associated with the two states with the average
reorientation rate of the large-scale circulation in the cell. This suggests
that the model could be used as a predictor of this rate in cases where few or
even no reorientations are observed.

## References

  * Bailon-Cuba et al. (2010) Bailon-Cuba, J., Emran, M. S. & Schumacher, J. 2010 Aspect ratio dependence of heat transfer and large-scale flow in turbulent convection. Journal of Fluid Mechanics 655, 152–173.
  * Belmonte & Libchaber (1996) Belmonte, Andrew & Libchaber, Albert 1996 Thermal signature of plumes in turbulent convection: The skewness of the derivative. Phys. Rev. E 53, 4893–4898.
  * Blei et al. (2003) Blei, David, Ng, Andrew & Jordan, Michael I. 2003 Latent Dirichlet allocation. Journal of Machine Learning Research 3, 993–1022.
  * Bosbach et al. (2012) Bosbach, Johannes, Weiss, Stephan & Ahlers, Guenter 2012 Plume fragmentation by bulk interactions in turbulent Rayleigh-Bénard convection. Phys. Rev. Lett. 108, 054501\.
  * Castaing et al. (1989) Castaing, B., Gunaratne, G., Heslot, F., Kadanoff, L., Libchaber, A., Thomae, S., Wu, X.-Z., Zaleski, S. & Zanetti, G. 1989 Scaling of hard thermal turbulence in Rayleigh-Bénard convection. Journal of Fluid Mechanics 204, 1–30.
  * Chilla & Schumacher (2012) Chilla, F. & Schumacher, J. 2012 New perspectives in turbulent Rayleigh-Bénard convection. European Physical Journal E 35 (7).
  * Ching et al. (2004) Ching, E. S.C., Guo, H., Shang, X.D., Tong, P. & Xia, K.-Q. 2004 Extraction of plumes in turbulent thermal convection. Physical Review Letters 93.
  * Emran & Schumacher (2012) Emran, Mohammad S. & Schumacher, Joerg 2012 Conditional statistics of thermal dissipation rate in turbulent Rayleigh-Bénard convection. European Physical Journal E 35 (10).
  * Fery et al. (2022) Fery, Lucas, Dubrulle, Bérengere, Podvin, Bérengere, Pons, Flavio & Faranda, Davide 2022 Learning a weather dictionary of atmospheric patterns using latent Dirichlet allocation. Geophysical Research Letters 49 (9), e2021GL096184, e2021GL096184 2021GL096184, arXiv: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2021GL096184.
  * Foroozani et al. (2017) Foroozani, N., Niemela, J. J., Armenio, V. & Sreenivasan, K. R. 2017 Reorientations of the large-scale flow in turbulent convection in a cube. Physical Review E 95, 033107\.
  * Frihat et al. (2021) Frihat, Mohamed, Podvin, Bérengère, Mathelin, Lionel, Fraigneau, Yann & Yvon, François 2021 Coherent structure identification in turbulent channel flow using latent Dirichlet allocation. Journal of Fluid Mechanics 920, A27.
  * Griffiths & Steyvers (2002) Griffiths, Thomas L. & Steyvers, Mark 2002 A probabilistic approach to semantic representation. In Proceedings of the 24th Annual Conference of the Cognitive Science Society.
  * Grossmann & Lohse (2000) Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: a unifying theory. Journal of Fluid Mechanics 407, 27–56.
  * Grossmann & Lohse (2004) Grossmann, S & Lohse, D 2004 Fluctuations in turbulent Rayleigh-Bénard convection: The role of plumes. Physics of Fluids 16 (12), 4462–4472.
  * Holmes et al. (2002) Holmes, P., Lumley, J.L., Berkooz, Gal & Rowley, C.W. 2002 Turbulence, Coherent Structures, Dynamical Systems and Symmetry. Cambridge University Press.
  * Huang et al. (2013) Huang, Shi-Di, Kaczorowski, Matthias, Ni, Rui & Xia, Ke-Qing 2013 Confinement-induced heat-transport enhancement in turbulent thermal convection. Phys. Rev. Lett. 111, 104501\.
  * Lumley (1967) Lumley, J.L. 1967 The structure of inhomogeneous turbulent flows. In Atmospheric Turbulence and Radio Wave Propagation (ed. A.M Iaglom & V.I Tatarski), pp. 221–227. Nauka, Moscow.
  * Podvin & Sergent (2015) Podvin, B. & Sergent, A. 2015 A large-scale investigation of wind reversal in a square Rayleigh-Bénard cell. Journal of Fluid Mechanics 766, 172–201.
  * Podvin & Sergent (2017) Podvin, B. & Sergent, A. 2017 Precursor for wind reversal in a square Rayleigh-Bénard cell. Physical Review E 95 (013112).
  * Puigjaner et al. (2008) Puigjaner, D., Herrero, J., Simó, C. & Giralt, F. 2008 Bifurcation analysis of steady Rayleigh-Bénard convection in a cubical cavity with conducting sidewalls. Journal of Fluid Mechanics 598, 393–427.
  * Rehurek & Sojka (2011) Rehurek, Radim & Sojka, Petr 2011 Gensim–python framework for vector space modelling. NLP Centre, Faculty of Informatics, Masaryk University, Brno, Czech Republic 3 (2).
  * Shang et al. (2003) Shang, X.-D., Qiu, X.-L., Tong, P. & Xia, K.-Q. 2003 Measured local heat transport in turbulent Rayleigh-Bénard convection. Phys. Rev. Lett. 90, 074501.
  * Shevkar et al. (2022) Shevkar, Prafulla P., Vishnu, R., Mohanan, Sanal K., Koothur, Vipin, Mathur, Manikandan & Puthenveettil, Baburaj A. 2022 On separating plumes from boundary layers in turbulent convection. Journal of Fluid Mechanics 941, A5.
  * Shishkina et al. (2010) Shishkina, O., Stevens, R. J. A. M., Grossmann, S. & Lohse, D. 2010 Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution. New Journal of Physics 12, 075022.
  * Shishkina & Wagner (2008) Shishkina, Olga & Wagner, Claus 2008 Analysis of sheet-like thermal plumes in turbulent Rayleigh-Bénard convection. Journal of Fluid Mechanics 599, 383–404.
  * Soucasse et al. (2019) Soucasse, L., Podvin, B., Rivière, Ph. & Soufiani, A. 2019 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.
  * Soucasse et al. (2021) Soucasse, L., Podvin, B., Rivière, Ph. & Soufiani, A. 2021 Low-order models for predicting radiative transfer effects on Rayleigh-Bénard convection in a cubic cell at different Rayleigh numbers. submitted to Journal of Fluid Mechanics 917, A5.
  * van der Poel et al. (2015) van der Poel, E. P., Verzicco, R., Grossmann, S. & Lohse, D. 2015 Plume emission statistics in turbulent Rayleigh-Bénard convection. Journal of Fluid Mechanics 772, 5–15.
  * Vishnu et al. (2020) Vishnu, Venugopal T., De, Arnab Kumar & Mishra, Pankaj Kumar 2020 Dynamics of large-scale circulation and energy transfer mechanism in turbulent Rayleigh-Bénard convection in a cubic cell. Physics of Fluids 32, 095115.
  * Vishnu et al. (2022) Vishnu, Venugopal T., De, Arnab Kumar & Mishra, Pankaj Kumar 2022 Statistics of thermal plumes and dissipation rates in turbulent Rayleigh-Bénard convection in a cubic cell. International Journal of Heat and Mass Transfer 182.
  * Wang et al. (2022) Wang, Yin, Wei, Yongze, Tong, Penger & He, Xiaozhou 2022 Collective effect of thermal plumes on temperature fluctuations in a closed Rayleigh-Bénard convection cell. Journal of Fluid Mechanics 934, A13.
  * Xi et al. (2004) Xi, H.-D., Lam, S. & Xia, K. 2004 From laminar plumes to organized flows: the onset of large-scale circulation in turbulent thermal convection. Journal of Fluid Mechanics 503, 47–56.
  * Xin et al. (2008) Xin, S., Chergui, J. & Le Quéré, P. 2008 3D spectral parallel multi-domain computing for natural convection flows. In Parallel Computational Fluid Dynamics, Lecture Notes in Computational Science and Engineering book series (ed. Springer), , vol. 74, pp. 163–171.
  * Xin & Le Quéré (2002) Xin, S. & Le Quéré, P. 2002 An extended Chebyshev pseudo-spectral benchmark for the 8:1 differentially heated cavity. Numerical Metholds in Fluids 40, 981–998.
  * Zhou et al. (2007) Zhou, Quan, Sun, Chao & Xia, Ke-Qing 2007 Morphological evolution of thermal plumes in turbulent Rayleigh-Bénard convection. Physical Review Letters 98 (7).
  * Zhou & Xia (2010) Zhou, Quan & Xia, Ke-Qing 2010 Physical and geometrical properties of thermal plumes in turbulent Rayleigh-Bénard convection. New Journal of Physics 12 (7), 075006.
  * Zhou et al. (2016) Zhou, Sheng-Qi, Xie, Yi-Chao, Sun, Chao & Xia, Ke-Qing 2016 Statistical characterization of thermal plumes in turbulent thermal convection. Physical Review Fluids 1 (5).

Figure 1:  Evolution of the horizontal components of the angular momentum at $Ra=10^{7}$. The vertical black lines correspond to reorientations of the large-scale circulation. $n=1$ |
---|---
$n=2$ |
$n=3$ |
Figure 2: POD dominant modes and amplitudes in the vertical mid-plane at
$Ra=10^{7}$. Left : POD modes $\phi_{n}$ , Right : POD amplitudes $a_{n}$
associated with plane $x=0.5H$ (in blue) and plane $y=0.5H$ (in red). The
vertical black lines correspond to changes in the component of the angular
momentum. The darker line corresponds to a moving average over 200 convective
units (4 $T_{c}$).

Figure 3:  Schematics of the LDA generative model illustrated here for a field defined on 12 cells and generated with 3 motifs (corresponding to the blue, red and purple colours). A snapshot $k$ is represented as a set of integer values defined on an array of cells. Let $\Phi_{T}$ be the sum of the integer values over the cells, which can be represented as a stack of $\Phi_{T}$ tokens of unit value. Each token is assigned in turn to a cell as follows: a motif is selected by sampling the snapshot-motif distribution $b^{k}$ corresponding to this snapshot. Once the motif $n$ is selected, a cell is selected by sampling the motif-cell distribution $\psi_{n}$. At the end of the process, the number of tokens at cell $j$ yields the value of the field $\Phi_{k}(\underline{x}_{j})$. As an alternative, the field can be generated by summing directly the contributions of all distributions $\psi_{n}$ weighted by the corresponding component $b_{n}$ (see also text). |
---|---
Figure 4:  left) Instantaneous correlation coefficient between the projected and the true field as a function of the integral convective heat flux for $N_{T}=100$ and $Ra=10^{7}$; right) Average correlation coefficient $<C(\Phi,\Phi^{R})>$ as a function of the Rayleigh number and of the number of topics considered for both mid-planes.  instantaneous field $\phi$ | 100 motifs | 20 dominant motifs | 20-mode POD
---|---|---|---
|  |  |
Figure 5:  Example of an instantaneous snapshot and its reconstructions at $Ra=10^{7}$. From left to right: original field, LDA-reconstructed field using $N_{T}=100$ motifs, LDA-reconstructed field using the 20 (instantaneously) most prevalent motifs, POD-reconstructed field using the 20 (on average) most energetic modes.  original field | LDA-projected field | POD-generated field | LDA-generated field
---|---|---|---
|  |  |
Figure 6:  Time-averaged value of the convective heat flux for different databases at $Ra=10^{7}$. From left to right: original fields, LDA-reconstructed (LDA-R) fields using $N_{T}=100$ motifs, POD-generated (POD-G) fields using 100 modes, LDA-generated (LDA-G) fields using $N_{T}=100$ modes.  A | B | C | D
---|---|---|---
|  |  |
|  |  |
Figure 7:  Autocorrelation of the convective heat flux at selected locations (see Figure 6). Top row: $R_{y}$ (horizontal direction); Bottom row: $R_{z}$ (vertical direction). See text (equations (14) and (15)) for definition.  A | B | C | D
---|---|---|---
|  |  |
Figure 8:  Probability density function of the convective heat flux at the selected locations indicated in Figure 6. $Ra=10^{6}$ | $Ra=10^{7}$ | $Ra=10^{8}$
---|---|---
|  |
Figure 9:  Spatial distribution of the positive-flux motifs in the vertical mid-plane for $N_{T}=100$. The motifs are materialized by a black line corresponding to a probability contour of $0.606\,\psi_{n}^{max}$. The vertical lines correspond to the boundary layer thickness. The time-averaged convective heat flux is represented in the background.  $Ra=10^{6}$ |  $Ra=3$ $10^{6}$ | $Ra=10^{7}$ | $Ra=10^{8}$
---|---|---|---
|  |  |
Figure 10:  Distribution of motif areas (see definition in equation (9)) in the vertical mid-plane with the distance from the lateral walls at varying Rayleigh numbers. The size of the symbols shown in the picture is proportional to the fraction of motifs over which the average was performed. The black solid lines indicate the boundary layer thickness. The dashed lines have slope 2. |
---|---
Figure 11:  left) Dominant motifs at Rayleigh numbers $Ra=10^{7}$ for $N_{T}=100$. The contour lines correspond to 0.1 $\psi^{max}$ and 0.3 $\psi^{max}$. The motif labels correspond to those of figure 12; right) Characteristic dominant motif at different Rayleigh numbers $N_{T}=100$. Isocontours of $\psi_{1}$ at $0.1i\psi^{max},i=2,\ldots 9$. The black lines correspond to the boundary layer thickness.  |
---|---
Figure 12:  Evolution of the snapshot-motif distributions $b_{n}$ for the four dominant motifs (see figure 11 for labels) at $Ra=10^{7}$ and for $N_{T}=100$: left) plane $x=0.5H$ right) plane $y=0.5H$. The thick line corresponds to a moving average over 150 convective units (3 recirculation times $T_{c}$). The horizontal dashed lines correspond to the values $b_{-}=0.017$ and $b_{+}=0.035$. The vertical lines correspond to the changes in angular momentum. |
---|---
Figure 13:  Left) Streamlines of the flow conditionally averaged on the high prevalence value of $b_{1}$; Right) Schematics of the cell organization in the vertical mid-plane - the large-scale circulation (in red) corresponds to the $I_{+}$ state while the corner structure (in black) corresponds to the $I_{-}$ state. $Ra=10^{6}$ |  $Ra=3$ $10^{6}$
---|---
|
$Ra=10^{7}$ | $Ra=10^{8}$
|
Figure 14:  Distribution of the dominant motif prevalence $b_{1}$ for
different Rayleigh numbers and $N_{T}=100$.

Figure 15:  Probability $p(T_{+}>T_{-})$ (see text) and comparison with ratio
of reorientation to recirculation time scale at different Rayleigh numbers -
the rescaling factor is $\beta^{-1}=5.6$.

[Acknowledgements] This work was granted access to the HPC resources of IDRIS
under the allocation 2023- AD012A62062R1 made by GENCI. We thank Anouar
Soufiani and Philippe Rivière for helpful discussions about the manuscript. We
are also grateful to Jean-Michel Dupays, Rémy Dubois and Camille Parisel for
technical support and helpful discussions.

[Declaration of Interests]The authors report no conflict of interest.


We will prove this inequality by contradiction. If this is not the case, then
there exists some sequence $\\{\theta^{m}:\,m\geq 1\\}$ such that
$D(\theta^{m})\rightarrow 0$ as $m\to\infty$ and

$\displaystyle\lim_{m\rightarrow\infty}\frac{1}{D(\theta^{m})}\,\int|p(x|\theta^{m})-p^{*}(x)|dx=0.$
(46)

Since $D(\theta^{m})\rightarrow 0$ as $m\to\infty$ and the number of possible
permutations of $[K^{*}]$ is finite, there is a sub-sequence
$\\{\theta^{m_{\ell}}\\}$ and a permutation
$\sigma^{\prime}\in\mathcal{S}_{K^{*}}$ such that
$w^{m_{\ell}}_{\sigma^{\prime}(k)}\rightarrow w^{*}_{k}$ and
$\eta^{m_{\ell}}_{\sigma^{\prime}(k)}\rightarrow\eta^{*}_{k}$ for all
$k=1,...,K^{*}$ as $\ell\to\infty$. Without loss of generality, we may assume
that $\sigma^{\prime}$ is the identity permutation and substitute the sequence
with its sub-sequence, so that we can write
$D(\theta^{m})=\sum_{k=1}^{K^{*}}\big{(}|w_{k}-w_{k}^{*}|+w_{k}|\eta_{k}-\eta_{k}^{*}|\big{)}$
for all sufficiently large $m$. Now we can apply a Taylor expansion as follows
to obtain:

$\displaystyle p(x|\theta^{m})-p^{*}(x)$ $\displaystyle\
=\sum_{k=1}^{K^{*}}w_{k}^{m}\big{(}g(x;\eta^{m}_{k})-g(x;\eta_{k}^{*})\big{)}+\sum_{k=1}^{K^{*}}(w_{k}^{m}-w_{k}^{*})g(x;\eta_{k}^{*})$
$\displaystyle\
=\sum_{k=1}^{K^{*}}w_{k}^{m}\Big{[}\big{(}\eta^{m}_{k}-\eta_{k}^{*}\big{)}\nabla
g(x;\eta_{k}^{*})+O\big{(}\eta^{m}_{k}-\eta_{k}^{*}\big{)}^{2}\Big{]}+\sum_{k=1}^{K^{*}}(w_{k}^{m}-w_{k}^{*})g(x;\eta_{k}^{*})$
$\displaystyle\
=\sum_{k=1}^{K^{*}}w_{k}^{m}\big{(}\eta^{m}_{k}-\eta_{k}^{*}\big{)}\nabla
g(x;\eta_{k}^{*})+\sum_{k=1}^{K^{*}}(w_{k}^{m}-w_{k}^{*})g(x;\eta_{k}^{*})+R_{m},$
(47)

where the remainder term satisfies $R_{m}/D(\theta^{m})\rightarrow 0$ as
$m\rightarrow\infty$. Note that equation (46) will then imply

$\displaystyle\lim_{m\rightarrow\infty}\frac{p(x|\theta^{m})-p^{*}(x)}{D(\theta^{m})}=\lim_{m\rightarrow\infty}\bigg{\\{}\sum_{k=1}^{K^{*}}\frac{w_{k}^{m}\big{(}\eta^{m}_{k}-\eta_{k}^{*}\big{)}}{D(\theta^{m})}\,\nabla
g(x;\eta_{k}^{*})+\sum_{k=1}^{K^{*}}\frac{w_{k}^{m}-w_{k}^{*}}{D(\theta^{m})}\,g(x;\eta_{k}^{*})\bigg{\\}}=0$

for almost every $x$. Note that the term inside the second limit is a linear
combination of finite many $g(x;\eta_{k}^{*})$ and $\nabla g(x;\eta_{k}^{*})$
terms. Since the parameter space is compact under Assumption A1, each linear
combination coefficient has a convergent sub-sequence. Without loss of
generality, we may substitute the sequence with its sub-sequence, and denote
the respective limits as

$\displaystyle\alpha_{k}=\lim_{m\to\infty}\frac{w_{k}^{m}\big{(}\eta^{m}_{k}-\eta_{k}^{*}\big{)}}{D(\theta^{m})}\
\ \mbox{and}\ \
\beta_{k}=\lim_{m\to\infty}\frac{w_{k}^{m}-w_{k}^{*}}{D(\theta^{m})}.$

Note that at least one of these limiting coefficients
$\\{\alpha_{k}\\}_{k=1}^{K^{*}}$ and $\\{\beta_{k}\\}_{k=1}^{K^{*}}$ does not
vanish, since otherwise, we will have

$\displaystyle
1=\lim_{m\rightarrow\infty}\frac{D(\theta^{m})}{D(\theta^{m})}=\lim_{m\rightarrow\infty}\frac{\sum_{k=1}^{K^{*}}w_{k}^{m}|\eta^{m}_{k}-\eta_{k}^{*}|+\sum_{k=1}^{K^{*}}|w_{k}^{m}-w_{k}^{*}|}{D(\theta^{m})}=\sum_{k=1}^{K^{*}}\big{(}|\alpha_{k}|+|\beta_{k}|\big{)}=0,$

which is a contradiction. Therefore, we constructed some
$\\{\alpha_{k}\\}_{k=1}^{K^{*}}$ and $\beta=(\beta_{1},\ldots,\beta_{K^{*}})$
such that they are not all zero, and

$\displaystyle\sum_{k=1}^{K^{*}}\Big{[}\alpha_{k}g(x;\eta_{k}^{*})+\beta^{T}\nabla
g(x;\eta_{k}^{*})\Big{]}=\lim_{m\rightarrow\infty}\frac{p(x|\theta^{m})-p^{*}(x)}{D(\theta^{m})}=0$

for almost every $x$. Note that the preceding display and the weak
identifiability condition together would imply that $\alpha_{k}=0$,
$\beta_{k}=\vec{0}$ for all $k=1,...,K^{*}$, which contradicts the fact that
at least one of the coefficients is not zero. Therefore, (45) holds. Since if
$\frac{1}{2}\int|p(x|\theta)-p^{*}(x)|dx=0$, then $\sum_{k=1}^{K^{*}}w_{k}=1$
and by Assumption A3, $D(\theta)=0$; if $D(\theta)=0$, it also implies
$\sum_{k=1}^{K^{*}}w_{k}=1$ and from Assumption A3 again we have
$\frac{1}{2}\int|p(x|\theta)-p^{*}(x)|dx=0$. Therefore,
$\frac{1}{2}\int|p(x|\theta)-p^{*}(x)|dx=0$ if and only if $D(\theta)=0$.
Because both $D(\theta)$ and $\frac{1}{2}\int|p(x|\theta)-p^{*}(x)|dx$ are
continuous in $\theta$, we can further obtain (following the same argument
from which we derived (11)) from equation (45) that for each possible
$\theta$,

$\displaystyle\frac{1}{2}\int|p(x|\theta)-p^{*}(x)|dx\gtrsim D(\theta).$ (48)

Step 2: Recall that we have already shown in Theorem 2 and Corollary 5 that
when $\phi_{0}<(d+1)/2$, there are $K^{*}$ of the mixing weights
$\overline{w}_{j}$’s being bounded away from 0, and the rest $K-K^{*}$ of the
mixing weights $\overline{w}_{j}$’s will empty out at the rate $(\log
n)^{\rho_{1}}/n$. Without loss of generality, we may assume that the first
$K^{*}$ of the weights $\overline{w}_{j}$’s do not vanish and denote
$\bm{w}_{K^{*}}=(w_{j})_{j=1}^{K^{*}}$ and
$\bm{\eta}_{K^{*}}=(\eta_{j})_{j=1}^{K^{*}}$. Then for any fixed
$\phi_{0}<(d+1)/2$, we have

$\displaystyle d_{\rm
TV}(p(\cdot\,|\,\overline{\bm{w}},\bm{\eta}),p^{*}(\cdot))$ $\displaystyle\
=\frac{1}{2}\int\Big{|}\sum_{j=1}^{K}\overline{w}_{j}g(x;\eta_{j})-\sum_{k=1}^{K^{*}}w_{k}^{*}g(x;\eta_{k}^{*})\Big{|}dx$
$\displaystyle\
\geq\frac{1}{2}\int\Big{|}\sum_{j=1}^{K^{*}}\overline{w}_{j}g(x;\eta_{j})-\sum_{k=1}^{K^{*}}w_{k}^{*}g(x;\eta_{k}^{*})\Big{|}dx-\frac{1}{2}\sum_{j=K^{*}+1}^{K}\overline{w}_{j}$
$\displaystyle\
\geq\frac{1}{2}\int\Big{|}\sum_{j=1}^{K^{*}}\overline{w}_{j}g(x;\eta_{j})-\sum_{k=1}^{K^{*}}w_{k}^{*}g(x;\eta_{k}^{*})\Big{|}dx-\frac{C(\log
n)^{\rho_{1}}}{n}$ $\displaystyle\ \gtrsim
D(\overline{\bm{w}}_{K^{*}},\bm{\eta}_{K^{*}})-\frac{C(\log
n)^{\rho_{1}}}{n}.$

where the last step is due to inequality (48). Recall that we have
$\sup_{\eta,\eta^{\prime}\in\Omega}|\eta-\eta^{\prime}|\leq L$ from Assumption
A1, which leads to the bound
$D(\overline{\bm{w}}_{K^{*}},\bm{\eta}_{K^{*}})\leq\sum_{k=1}^{K^{*}}(w_{k}+w_{K}^{*}+w_{k}L)\leq
2+L$. Now since $\big{(}C(\log n)^{\rho_{1}}/n\big{)}^{2}>0$, we can then get

$d^{2}_{TV}(p(\cdot\,|\,\overline{\bm{w}},\bm{\eta}),p^{*}(\cdot))\gtrsim
D^{2}(\overline{\bm{w}}_{K^{*}},\bm{\eta}_{K^{*}})-\frac{2C(L+2)(\log
n)^{\rho_{1}}}{n}.$

Step 3: In this step, we show that with high probability, the variational
posterior distribution of each non-emptying out mixture component parameter
$\eta_{k}$ is concentrated within a neighborhood of $\overline{\eta}_{k}$ with
radius $Mn^{-1/2}$. Recall that from Corollary 5, for each non-emptying out
component $k$, we have $\widehat{n}_{k}\geq c_{k}n$ for some positive constant
$c_{k}$. In this step, we may again assume without loss of generality that the
first $K^{*}$ of the weights $\overline{w}_{j}$’s do not vanish. Recall that
from Theorem 1, we have that

$\displaystyle\widehat{q}_{\eta_{k}}(\eta_{k})\propto\pi(\eta_{k})\exp\left[\sum_{i=1}^{n}\widehat{p}_{ik}\left(\eta_{k}^{T}T(x_{i})-T_{0}(x_{i})-A(\eta_{k})\right)\right],$

where $\sum_{i=1}^{n}\widehat{p}_{ik}=\widehat{n}_{k}$. Recall that
$\widehat{\eta}_{k}$, defined in equation (27), denotes the $\eta_{k}$ that
maximizes the exponent in the preceding display. We will first show below that

$\displaystyle\min_{1\leq j\leq
K^{*}}|\widehat{\eta}_{j}-\eta^{*}_{k}|\lesssim\frac{\log
n}{\sqrt{\widehat{n}_{k}}}$ (49)

holds for each $1\leq k\leq K^{*}$.

From equations (28) and (29), and Assumption A2, we can get

$\displaystyle\pi(\eta_{k})\exp\left[\sum_{i=1}^{n}\widehat{p}_{ik}\left(\eta_{k}^{T}T(x_{i})-T_{0}(x_{i})-A(\eta_{k})\right)-\ell_{nk}(\widehat{\eta}_{k})\right]\geq\pi(\eta_{k})\exp\Big{[}-\frac{b\widehat{n}_{k}}{2}(\eta_{k}-\widehat{\eta}_{k})^{2}\Big{]},$

which implies that for any positive number $C$,

$\displaystyle\int_{\\{|\eta_{k}-\widehat{\eta}_{k}|_{\infty}\leq
C(b\widehat{n}_{k})^{-1/2}\\}}\widehat{q}_{\eta_{k}}(\eta_{k})d\eta_{k}\geq\inf_{\\{|\eta_{k}-\widehat{\eta}_{k}|_{\infty}\leq
C(b\widehat{n}_{k})^{-1/2}\\}}\frac{\pi(\eta_{k})}{\pi(\widehat{\eta}_{k})}\Bigg{(}\frac{a\big{(}2\Phi(C)-1\big{)}}{b}\Bigg{)}^{d/2},$

where $\Phi$ denotes the CDF of the standard normal distribution. In
particular, by choosing $C=\Phi^{-1}(3/4)$, the integral above can be bounded
from below by $c_{\pi}(a/2b)^{d/2}$, for some positive constant $c_{\pi}$ when
$n$ is sufficiently large. Therefore, for constant $C_{k}=C(d/bc_{k})^{1/2}$,
we have the following inequality,

$\displaystyle\int_{\\{|\eta_{k}-\widehat{\eta}_{k}|_{\infty}\leq
C_{k}n^{-1/2}\\}}\widehat{q}_{\eta_{k}}(\eta_{k})d\eta_{k}\geq
c_{\pi}(a/2b)^{d/2}.$

We can now prove (49) by contradiction. Suppose $\min_{1\leq j\leq
K^{*}}|\widehat{\eta}_{j}-\eta^{*}_{k^{\prime}}|\geq\max_{1\leq j\leq
K^{*}}C_{j}\log n/n^{1/2}$ holds for some $k^{\prime}\leq K^{*}$, then each
$\eta_{j}$ in the set
$\mathcal{B}:=\bigcap_{j=1}^{K^{*}}\\{\bm{\eta}:|\eta_{j}-\widehat{\eta}_{j}|\leq
C_{j}n^{-1/2}\\}$ cannot be close to $\eta_{k}^{*}$ by the triangle
inequality, i.e.,

$\displaystyle\min_{1\leq j\leq K^{*}}|\eta_{j}-\eta^{*}_{k}|\geq\min_{1\leq
j\leq
K^{*}}\big{(}|\widehat{\eta}_{j}-\eta^{*}_{k^{\prime}}|-|\eta_{j}-\widehat{\eta}^{*}_{k}|\big{)}\geq\min_{1\leq
j\leq K^{*}}\frac{C_{j}(\log n-1)}{n^{1/2}}$

holds for any $\eta=(\eta_{j}:\,j\in[K^{*}])\in\mathcal{B}$. Since from Step
1, we obtain $d^{2}_{TV}(p(x|\overline{\bm{w}},\bm{\eta}),p^{*}(x))\gtrsim
D^{2}(\overline{\bm{w}}_{K^{*}},\bm{\eta}_{K^{*}})-2C(L+2)(\log
n)^{\rho_{1}}/n$, which furhter implies

$\displaystyle\int_{\Omega^{K}}h^{2}\big{(}p(x\,|\,\overline{\bm{w}},\bm{\eta}),p^{*}(x)\big{)}\widehat{q}(\bm{\eta})\,d\bm{\eta}$
$\displaystyle\ \geq\int_{\mathcal{B}}d_{\rm
TV}^{2}\big{(}p(x\,|\,\overline{\bm{w}},\bm{\eta}),p^{*}(x)\big{)}\widehat{q}(\bm{\eta})\,d\bm{\eta}$
$\displaystyle\
\gtrsim\int_{\mathcal{B}}D^{2}(\overline{\bm{w}}_{K^{*}},\bm{\eta}_{K^{*}})\widehat{q}(\bm{\eta})\,d\bm{\eta}-2CL(\log
n)^{\rho_{1}}/n$ $\displaystyle\ \gtrsim(\log
n-1)n^{-1/2}\widehat{Q}(\mathcal{B})-2CL(\log n)^{\rho_{1}}/n$ $\displaystyle\
\geq(\log n-1)\,n^{-1/2}\big{(}c_{\pi}(a/2b)^{d/2}\big{)}^{K^{*}}-2CL(\log
n)^{\rho_{1}}/n.$

This is a contradiction with Lemma 9, which implies the left-hand side of the
preceding display to be at most $\mathcal{O}_{P}(n^{-1})$. This proves
inequality (49).

Lastly, we show that the variational posterior mode $\widehat{\eta}_{k}$ is
close to the variational posterior mean $\overline{\eta}_{k}$ for each
$k\in[K^{*}]$. In fact, since $\overline{\eta}_{k}$ is given by

$\displaystyle\overline{\eta}_{k}=\frac{\int\eta_{k}\pi(\eta_{k})\exp\\{\ell_{nk}(\eta_{k})-\ell_{nk}(\widehat{\eta}_{k})\\}d\eta_{k}}{\int\pi(\eta_{k})\exp\\{\ell_{nk}(\eta_{k})-\ell_{nk}(\widehat{\eta}_{k})\\}d\eta_{k}},$

we may use the Laplace formula again to obtain

$\displaystyle\left|\frac{\int\eta_{k}\pi(\eta_{k})\exp\\{\ell_{nk}(\eta_{k})-\ell_{nk}(\widehat{\eta}_{k})\\}d\eta_{k}}{(2\pi/n_{k})^{\frac{d}{2}}}-\frac{\widehat{\eta}_{k}\pi(\widehat{\eta}_{k})}{\sqrt{I(\widehat{\eta}_{k})}}\right|\leq\frac{C}{n_{k}}.$

This inequality combined with (29) leads to
$|\widehat{\eta}_{k}-\overline{\eta}_{k}|\leq C/n_{k}$; and the conclusion
follows by combining this with (49) and using the fact that
$\widehat{n}_{k}\asymp n$ for $k\in[K^{*}]$.