# Flow states and heat transport in Rayleigh–Bénard convection with different
sidewall boundary conditions

Philipp Reiter1These authors contributed equally<EMAIL_ADDRESS>Xuan Zhang1${\ddagger}$ Olga Shishkina1<EMAIL_ADDRESS>1Max Planck
Institute for Dynamics and Self-Organization,
Am Fassberg 17, 37077 Göttingen, Germany

###### Abstract

This work addresses the effects of different thermal sidewall boundary
conditions on the formation of flow states and heat transport in two- and
three-dimensional Rayleigh–Bénard convection (RBC) by means of direct
numerical simulations and steady-state analysis for Rayleigh numbers $Ra$ up
to $4\times 10^{10}$ and Prandtl numbers $Pr=0.1,1$ and $10$. We show that a
linear temperature profile imposed at the conductive sidewall leads to a
premature collapse of the single-roll state, whereas a sidewall maintained at
a constant temperature enhances its stability. The collapse is caused by
accelerated growth of the corner rolls with two distinct growth rate regimes
determined by diffusion or convection for small or large $Ra$, respectively.
Above the collapse of the single-roll state, we find the emergence of a
double-roll state in two-dimensional RBC and a double-toroidal state in three-
dimensional cylindrical RBC. These states are most prominent in RBC with
conductive sidewalls. The different states are reflected in the global heat
transport, so that the different thermal conditions at the sidewall lead to
significant differences in the Nusselt number for small to moderate $Ra$.
However, for larger $Ra$, heat transport and flow dynamics become increasingly
alike for different sidewalls and are almost indistinguishable for
$Ra>10^{9}$. This suggests that the influence of imperfectly insulated
sidewalls in RBC experiments is insignificant at very high $Ra$ \- provided
that the mean sidewall temperature is controlled.

###### keywords:

Rayleigh–Bénard convection, Turbulent convection, Computational methods

##  1 Introduction

Understanding thermally induced convection as it arises in the earth’s
atmospheric/oceanic circulations and deducing its fundamental aspects from
laboratory experiments is an ongoing endeavour which motivated numerous
experimental and theoretical studies. In this realm, Rayleigh–Bénard
convection (RBC), i.e. a fluid held between two parallel plates heated from
below and cooled from above, is the most thoroughly investigated model system
to study the complex physics behind natural convection such as pattern
formation and the transition to turbulence (Bodenschatz et al., 2000; Ahlers
et al., 2009b; Lohse & Xia, 2010).

Most of the early theoretical advances were made by considering the system as
infinitely extended in the lateral direction. For instance, conventional
linear-stability analysis predicts the formation of two-dimensional rolls
(Chandrasekhar, 1961), while a weakly non-linear analysis reveals the
stability regimes of these rolls and their path to subsequent oscillatory or
stationary type bifurcations (Schlüter et al., 1965; Busse, 1967, 1978). In
laboratory experiments, however, we must resort to laterally confined systems
where our understanding is far less complete. In particular, when the lateral
size of the container is close to or less than the height of the cell, the
presence of sidewalls plays an important role (Roche, 2020; Shishkina, 2021).
Therefore, this study focuses on the effects of different thermal sidewall
boundary conditions on heat transfer and the emergence of different flow
states.

Different sidewalls are known to affect the critical Rayleigh number $Ra_{c}$
above which convection sets in (Buell & Catton, 1983; Hébert et al., 2010),
and perfectly conducting sidewalls have been found to delay the onset compared
to adiabatic sidewalls. In an attempt to better understand the flow regimes
above onset, bifurcation analyses were performed in a cubic domain for
adiabatic (Puigjaner et al., 2004) and perfectly conducting sidewalls
(Puigjaner et al., 2008). The bifurcation diagrams for the conducting
sidewalls are generally more complex, and double-toroidal states predominate
over the classical single-roll structure found for adiabatic sidewalls.
Sidewalls also have a strong influence on pattern formation (Cross &
Hohenberg, 1993; de Bruyn et al., 1996; Bodenschatz et al., 2000) and
different sidewall boundary conditions lead to differences in observable
patterns even in cells with large aspect ratio (Hu et al., 1993).

In RBC experiments, spurious sidewall heat fluxes are a major practical
difficulty that can substantially bias global heat transport measurements.
Ahlers (2000) reported that naive sidewall corrections can overstate Nusselt
number measurements by up to $20\%$ and underestimate the scaling of the
Nusselt number $Nu$ with respect to the Rayleigh number $Ra$ ($Nu\sim
Ra^{\lambda}$) reflected in the reduction of the scaling exponent $\lambda$ by
about $2\%$, underscoring the importance of more sophisticated sidewall
corrections. Roche et al. (2001) further emphasized this conclusion by showing
that the sidewall corrections can be considerably larger than assumed, leading
to scaling exponents closer to the turbulent scaling of $Nu\sim Ra^{1/3}$
(Grossmann & Lohse, 2000, 2001, 2004) than previously measured. Probably the
most important question in convection today is whether the ultimate regime in
confined geometries has the same scaling as predicted for unbounded domains,
i.e. $Nu\sim Ra^{1/2}$ (up to different logarithmic corrections), as proposed
by Kraichnan (1962) and Grossmann & Lohse (2011). Another important question
is when and how exactly the transition to the ultimate regime takes place in
confined geometries. Laboratory experiments (Chavanne et al., 1997; Niemela et
al., 2000; Chavanne et al., 2001; Ahlers et al., 2009a, 2012; He et al., 2012;
Urban et al., 2014; Roche, 2020) in this extremely high $Ra$ regime are
notoriously difficult to perform and potentially sensitive to several unknowns
of the system, one of which is the influence of imperfectly isolated/adiabatic
sidewalls.

Numerical simulations were performed incorporating thermal conduction in the
solid sidewall to clarify the differences between an ideal adiabatic setup and
a finite thermal conductivity sidewall (Verzicco, 2002; Stevens et al., 2014;
Wan et al., 2019). The results of these studies suggest that different thermal
properties of the sidewall alter the mean flow structure, leading to
significant differences in global heat transport in the low to mid $Ra$ range.
However, this effect vanishes for larger $Ra$, at least when the sidewall
temperature is constant and maintained at the arithmetic mean of upper and
lower plate temperatures. Conversely, if the sidewall temperature deviates
from the arithmetic mean, differences in heat transport persist even for large
$Ra$. This indicates that it is more important to keep the environment at the
correct temperature than to shield the interior of the cell from its
surroundings.

Despite extensive previous work, the spatial distribution of flow and heat
transport in confined geometries with different thermal boundary condition has
not been exhausted, especially the conditions related to real experimental
sidewall boundary conditions. In the present work, we investigate RBC with the
following thermal sidewall boundary conditions: adiabatic, constant
temperature (isothermal) and linear temperature. In the first part of the
results, we focus on a steady-state analysis based on an adjoint descent
algorithm (Farazmand, 2016) to identify different flow states, their
properties and their evolution over $Ra$. In the second part, the analysis is
complemented and extended to higher $Ra$ into the turbulent regime by a set of
DNS for a 2D box and 3D cylindrical setup, covering a range of
$10^{3}<Ra<10^{11}$ and $10^{3}<Ra<10^{9}$, respectively, aiming for a more
complete picture. We first present our numerical methods, discuss the results
and conclude with our main findings.

##  2 Numerical methods

###  2.1 Governing equations

The dimensionless control parameters in RBC are the Rayleigh number
$\mbox{{Ra}}\equiv\alpha g\Delta H^{3}/(\kappa\nu)$, the Prandtl number
$\Pran\equiv\nu/\kappa$, and the width-to-height aspect ratio of the box,
$\Gamma\equiv L/H$. Here, $\alpha$ denotes the isobaric thermal expansion
coefficient, $\nu$ the kinematic viscosity, $\kappa$ the thermal diffusivity
of the fluid, $g$ the acceleration due to gravity, $\Delta\equiv T_{+}-T_{-}$
the difference between the temperatures at the lower ($T_{+}$) and upper
($T_{-}$) plates, $H$ the distance between the parallel plates (the container
height), and $L$ the length of the container or the diameter in the case of a
cylindrical setup. In this study, we focus on variations with $Ra$, while
$Pr=1$ is fixed for most results in this paper except for a $Pr$-dependence
study in section 4.5, and $\Gamma=1$ is held constant throughout the study.

The governing equations in the Oberbeck–Boussinessq approximation for the
dimensionless, incompressible velocity ${\bf u}$, temperature $\theta$ and
kinematic pressure $p$ read as follows:

$\displaystyle\partial{\bf u}/\partial t+{\bf u}\cdot{\bm{\nabla}}{\bf
u}+{\bm{\nabla}}{p}$ $\displaystyle=$
$\displaystyle\sqrt{Pr/Ra}{\bm{\nabla}}^{2}{\bf u}+{\theta}{\bf e}_{z},$
$\displaystyle\partial{\theta}/\partial t+{\bf u}\cdot{\bm{\nabla}}{\theta}$
$\displaystyle=$ $\displaystyle
1/\sqrt{PrRa}{\bm{\nabla}}^{2}{\theta},\quad{\bm{\nabla}}\cdot{\bf u}=0.$ (1)

The equations were made dimensionless using the free-fall velocity
$u_{ff}\equiv(\alpha g\Delta H)^{1/2}$, the free-fall time $t_{ff}\equiv
H/u_{ff}$, the temperature difference $\Delta\equiv T_{+}-T_{-}$ between
bottom ($T_{+}$) and top ($T_{-}$) plates and $H$ the cell height. Here ${\bf
e}_{z}$ is the unit vector in the vertical $z$-direction. This set of
equations is solved with the direct numerical solver goldfish, which uses a
fourth-order finite volume discretization on a staggered grid and a third
order Runge–Kutta time scheme. The code has been widely used in previous
studies and validated against other direct numerical simulation codes (Kooij
et al., 2018; Reiter et al., 2021a).

###  2.2 Boundary conditions

$(a)$$(b)$$(c)$$(d)$ Figure 1:  2D Numerical setup of $(a)$ adiabatic, $(b)$
linear and $(c)$ constant sidewall temperature boundary conditions. $(d)$
Sketch of cylindrical domain. Profiles next to $(b)$ and $(c)$ show the
imposed sidewall temperature distribution.

We study 2D RBC in a square box and 3D RBC in a cylindrical domain. The setups
and profiles of the sidewall (SW) boundary conditions (BCs) used are shown in
figure 1. The adiabatic, linear and constant conditions for the sidewall
region $\delta V_{S}$ are defined by

adiabatic: $\displaystyle\quad\partial\theta/\partial{\chi}=0,$ (2) linear:
$\displaystyle\quad\theta=\theta_{+}+z\left(\theta_{-}-\theta_{+}\right),$ (3)
constant:
$\displaystyle\quad\theta=\begin{cases}\frac{-k(2z-1)}{k+2z}\left(\theta_{+}-\theta_{m}\right),&0\leq
z\leq 1/2,\\\
\frac{k(2z-1)}{k-2z+2}\left(\theta_{-}-\theta_{m}\right),&1/2<z\leq
1,\end{cases}$ (4)

with the temperature of the lower plate $\theta_{+}=1/2$, the temperature of
the upper plate $\theta_{-}=-1/2$, their arithmetic mean $\theta_{m}=0$,
$z\equiv z/H\in[0,1]$ and $\chi=x$ for box and $\chi=r$ for cylinder,
respectively. As for the constant temperature conditions, most of the sidewall
is kept at a nearly uniform temperature ($\theta_{m}$), except for the
transition regions in the vicinity of the top and bottom plates to ensure a
smooth temperature distribution. The parameter $0<k\ll 1$ in eq. (4) defines
the thickness of the transition layer. Here we used $k=0.01$, which gives a
fairly sharp albeit sufficiently smooth transition, as can be seen in figure 1
$(c)$. Moreover, the velocity no-slip conditions apply to all walls, i.e.
$\mathbf{u}\evaluated{}_{\text{missing}}{wall}=0$.

###  2.3 Adjoint descent method

A complementary analysis to direct numerical simulations is the study of the
Boussinesq equations by means of its invariant solutions. Hopf (1948)
conjectured that the solution of the Navier–Stokes equations can be understood
as a finite but possibly large number of invariant solutions, and turbulence
from this point of view is the migration from the neighbourhood of one
solution to another. While highly chaotic systems seem hopelessly complex to
understand, laminar or weakly chaotic flows can often be captured quite well
with this approach. In this work, we focus solely on solutions for steady-
states (equilibrium).

Determining steady-state solutions can be quite difficult, especially when the
number of dimensions is large as it is the case for most fluid mechanical
problems. The most commonly used numerical method for this task is Newton’s
method, which usually uses the generalized minimal residual (GMRES) algorithm
to solve the corresponding systems of linear equations (Saad & Schultz, 1986).
This method generally shows fast convergence rates when the initial estimate
is close to the equilibrium point. However, if the initial estimate is too far
from the equilibrium, Newton’s method often fails. In particular, for fluid
mechanics, the basin of attraction of Newton’s method can be quite small,
making the search for steady-states highly dependent on the initial guess.
Here we consider an alternative approach recently proposed by Farazmand (2016)
based on an adjoint method. Farazmand (2016) has shown that this adjoint-
descent method can significantly improve the chance of convergence compared to
the Newton–descent method, and thus more reliably capture equilibrium states
from a given initial state, but at the cost of a generally slower convergence
rate. A detailed derivation of the algorithm can be found in Farazmand (2016).
Below we sketch the idea of the method.

Suppose we want to find equilibrium solutions of a particular PDE (in our case
the Boussinessq equations)

$\partial_{t}{\bf u}=F({\bf u}),$ (5)

with ${\bf u}={\bf u}(\mathbf{x},t)$. The equilibrium’s of F(u) can be
generally unstable and therefore difficult to detect. The idea is to search a
new PDE, i.e.

$\partial_{\tau}{\bf u}=G({\bf u}),$ (6)

which solutions always converge to the equilibrium solutions of (5) when the
fictitious time $\tau$ goes to infinity

$\norm{F({\bf u})}_{\mathcal{A}}^{2}\rightarrow
0\quad\text{as}\quad\tau\rightarrow\infty,$ (7)

with the weighted energy norm
$\norm{\cdot}_{\mathcal{A}}\equiv\langle\cdot,\cdot\rangle_{\mathcal{A}}\equiv\langle\cdot,\mathcal{A}\cdot\rangle$
for a certain real self-adjoint and positive definite operator $\mathcal{A}$.
$F({\bf u})$ evolves along a trajectory ${\bf u}^{\prime}$ in accordance with

$\frac{1}{2}\partial_{\tau}\norm{F({\bf u})}_{\mathcal{A}}^{2}=\langle\delta
F({\bf u},{\bf u}^{\prime}),F({\bf u})\rangle_{\mathcal{A}},$ (8)

where $\delta F({\bf u},{\bf u}^{\prime})\equiv\lim\limits_{\varepsilon\to
0}\frac{F({\bf u}+\varepsilon{\bf u}^{\prime})-F({\bf u})}{\varepsilon}$ of
$F({\bf u})$ is the functional Gateaux derivative at ${\bf u}$ in the
direction ${\bf u}^{\prime}$. In the Newton-descent method, the search
direction ${\bf u}^{\prime}$ is approximated from $\delta F({\bf u},{\bf
u}^{\prime})=-F({\bf u})$ by using, for example, a GMRES iterative algorithm.
For the adjoint-descent method, on the other hand, we rewrite eq. (8) in the
form

$\frac{1}{2}\partial_{\tau}\norm{F({\bf u})}_{\mathcal{A}}^{2}=\langle{\bf
u}^{\prime},\delta F^{\dagger}({\bf u},F({\bf u}))\rangle_{\mathcal{A}},$ (9)

where $\delta F^{\dagger}$ is the adjoint operator of the functional
derivative $\delta F$. For ${\bf u}^{\prime}=-\delta F^{\dagger}({\bf
u},F({\bf u}))$ one guarantees that $\norm{F({\bf u})}_{\mathcal{A}}^{2}$
decays to zero along the trajectory ${\bf u}^{\prime}$, since then
$\frac{1}{2}\partial_{\tau}\norm{F({\bf u})}_{\mathcal{A}}^{2}=-\norm{\delta
F^{\dagger}({\bf u},F({\bf u}))}_{\mathcal{A}}^{2}$. Letting ${\bf u}$ evolve
along the adjoint search direction ensures the convergence to an equilibrium,
thus we find the desired PDE $G({\bf u})\equiv{\bf u}^{\prime}$, i.e.

$G({\bf u})=-\delta F^{\dagger}({\bf u},F({\bf u})).$ (10)

The choice of the norm $\norm{\cdot}_{\mathcal{A}}$ is important for the
algorithm to be numerically stable and is explained in more detail in the
appendix. As mentioned, the operator $\mathcal{A}$ should be real-valued,
positive-definite and self-adjoint. Following Farazmand (2016), we use an
operator $\mathcal{A}$ that is closely related to the inversed Laplacian, i.e.
$\mathcal{A}=(I-\alpha{\bm{\nabla}}^{2})^{-1}$ where $I$ is the identity
operator and $\alpha$ is a non-negative scalar parameter. For $\alpha=0$ this
norm converges to the $L^{2}$-norm and for $\alpha>0$ it effectively dampens
smaller scales and provides a better numerical stability.

The linear adjoint equations for the Boussinesq equations (1) read

$\displaystyle-\partial_{\tau}{\bf u}$
$\displaystyle=\left({\bm{\nabla}}\tilde{{\bf
u}}^{\prime\prime}+({\bm{\nabla}}\tilde{{\bf
u}}^{\prime\prime})^{\text{T}}\right){\bf
u}+\theta{\bm{\nabla}}\tilde{\theta}^{\prime\prime}-{\bm{\nabla}}p^{\prime\prime}+\sqrt{Pr/Ra}{\bm{\nabla}}^{2}\tilde{{\bf
u}}^{\prime\prime},$ $\displaystyle-\partial_{\tau}\theta$ $\displaystyle={\bf
u}\cdot{\bm{\nabla}}\tilde{\theta}^{\prime\prime}+1/\sqrt{PrRa}{\bm{\nabla}}^{2}\tilde{\theta}^{\prime\prime}+\mathbf{\mathbf{e}}_{z}\cdot\tilde{{\bf
u}}^{\prime\prime},$ $\displaystyle{\bm{\nabla}}\cdot{\bf u}^{\prime\prime}$
$\displaystyle=0,\quad{\bm{\nabla}}\cdot{\bf u}=0$ (11)

(see derivations in the appendix). Here the double prime fields ${\bf
u}^{\prime\prime}$ and $\theta^{\prime\prime}$ denote the residuals of the
Navier–Stokes eq. (1), i.e.

$\displaystyle{\bf u}^{\prime\prime}$ $\displaystyle\equiv-{\bf
u}\cdot{\bm{\nabla}}{\bf u}-{\bm{\nabla}}p+\sqrt{Pr/Ra}{\bm{\nabla}}^{2}{\bf
u}+\mathbf{\mathbf{e}}_{z}\theta,$ $\displaystyle\theta^{\prime\prime}$
$\displaystyle\equiv-{\bf
u}\cdot{\bm{\nabla}}\theta+1/\sqrt{PrRa}{\bm{\nabla}}^{2}\theta.$ (12)

and $\tilde{{\bf u}}^{\prime\prime}\equiv\mathcal{A}\mathbf{{\bf
u}}^{\prime\prime}$ as well as
$\tilde{\theta}^{\prime\prime}\equiv\mathcal{A}\mathbf{\theta}^{\prime\prime}$.
For simplicity, let $\mathbf{q}\equiv({\bf u},\theta)$, then the adjoint
descent method consists of three steps

  1. 1.

Find the residuals $\mathbf{q}^{\prime\prime}$ according to eq. (12).

  2. 2.

Solve $\tilde{\mathbf{q}}^{\prime\prime}=\mathcal{A}\mathbf{q}^{\prime\prime}$
for $\tilde{\mathbf{q}}^{\prime\prime}$.

  3. 3.

Update $\mathbf{q}$ according to eq. (11).

In step (i), we solve the time-stepping eq. (1), where we use a standard
pressure projection method and treat the diffusion term implicitly. The time
step size $\Delta t$ can be chosen independently of the artificial time step
size $\Delta\tau$ of the adjoint equations. For step (ii), using the energy
norm $\norm{\cdot}_{\mathcal{A}}$ with the operator
$\mathcal{A}=(I-\alpha{\bm{\nabla}}^{2})^{-1}$, we solve the Helmholtz-type
equation
$(I-\alpha{\bm{\nabla}}^{2})\tilde{\mathbf{q}}^{\prime\prime}=\mathbf{q}^{\prime\prime}$.
The integration of the adjoint equations in step (iii) is similar to step (i),
but all terms are treated explicitly. Through tests, we found that the
artificial time step $\Delta\tau$ can be chosen much larger than $\Delta t$ in
some cases, i.e. for large $Ra$.

The boundary conditions of $\tilde{{\bf u}}^{\prime\prime}$ and
$\tilde{\theta}^{\prime\prime}$ result from integration by parts in the
derivation of the adjoint equations. Evaluation of the adjoint operator of the
diffusion terms yields

$\displaystyle\int_{V}\tilde{{\bf u}}^{\prime\prime}{\bm{\nabla}}^{2}{\bf
u}^{\prime}dV=\int_{V}{\bf u}^{\prime}{\bm{\nabla}}^{2}\tilde{{\bf
u}}^{\prime\prime}dV+\int_{S}{\bf u}^{\prime}({\bm{\nabla}}\tilde{{\bf
u}}^{\prime\prime}\cdot\mathbf{n})dS-\int_{S}\tilde{{\bf
u}}^{\prime\prime}({\bm{\nabla}}{\bf u}^{\prime}\cdot\mathbf{n})dS,$ (13)

where we see the occurrence of two additional boundary terms (the last two
terms) evaluated on the boundary domain $S$. The first boundary term vanishes
since the search direction ${\bf u}^{\prime}$ is zero on the boundaries. The
second term can be eliminated if we also choose homogeneous Dirichlet boundary
conditions for the adjoint field $\tilde{{\bf u}}^{\prime\prime}$ on $S$. The
same logic applies to homogeneous Neumann conditions. For the pressure field
$p^{\prime\prime}$, we apply Neumann boundary conditions conditions on all
walls. In this study, all flow states showed good overall convergence
($\norm{F({\bf u})}_{\mathcal{A}}^{2}\leq 10^{-5}$) and the velocity fields
where almost divergence free ($\norm{\divergence{{\bf u}}}_{L^{2}}\leq
10^{-3}$). However, the rigorous verification of the chosen pressure BCs has
yet to be performed. Another interesting point, reserved for later
investigation, is whether a vorticity-streamfunction formulation might be
better suited to resolve issues with the boundary conditions.

Figure 2: Convergence of the adjoint-descent method for three different $Ra$,
starting from the same initial field. The time-step size for which the
algorithm is just stable increased with $Ra$, i.e., for these cases we used
$\Delta\tau=0.5$ ($Ra=10^{4}$), $\Delta\tau=2.0$ ($Ra=10^{5}$) and
$\Delta\tau=5.0$ ($Ra=10^{6}$). All three cases converged to large-scale
circulation flow states as described in section 3.2.

For the steady-state analysis, we use a Galerkin method with Chebyshev bases
in $x$ and $z$ directions and a quasi-inverse matrix diagonalization strategy
for better efficiency (Shen, 1995; Julien & Watson, 2009; Oh, 2019; Mortensen,
2018). The code is publicly available (Reiter, 2021). We use an implicit
backward Euler time discretization and alias the fields using the $2/3$ rule
by setting the last $1/3$ high-frequency spectral coefficients to zero after
evaluating the nonlinear terms. When used as a direct numerical solver, we
found excellent agreement with our finite-volume code goldfish. In addition,
the steady-states from the adjoint descent method showed excellent agreement
with those found by an alternative Newton–GMRES iteration. Figure 2 shows the
convergence rates for three different $Ra$, starting from the same initial
state. Overall, we find that the convergence chance is improved over the
Newton-descent method, although the convergence rate suffers and larger $Ra$
are either not feasible with the current approach as implemented in our code
or diverge after some time. Therefore, we restrict the steady-state analysis
to flows in the range $Ra\leq 10^{7}$ and investigate larger $Ra$ using direct
numerical simulations. One conceivable problem with the current approach is
that the currently used energy norm with the operator
$\mathcal{A}\equiv(I-\alpha{\bm{\nabla}}^{2})^{-1}$ dampens smaller scales in
order to increase the stability of the algorithm. But for larger $Ra$, smaller
scales become important to resolve the boundary layers sufficiently, so the
algorithm is likely to take longer to converge or the damping of the smaller
scales is too severe to reach convergence overall. Using smaller values of
$\alpha$ could lead to better results in that case, as it emphasizes smaller
scales more. Preliminary analysis suggests that $\alpha=10^{-3}$ leads to
better convergence to a steady-state than $\alpha=1$, but requires smaller
time steps $\delta\tau$, which currently makes it too costly to apply to a
wider range of parameters. In the future, the convergence rate might be
improved by employing a hybrid adjoint-descent and Newton-GMRES approach, as
proposed by Farazmand (2016). Alternative gradient optimization techniques are
also conceivable to boost convergence speed.

##  3 Steady-state analysis

$(a)$$(b)$$(c)$Adiabatic SWLinear SWConstant SW Figure 3: Growth rates
$\sigma$ as determined from linear stability analysis for the four most
unstable modes at the onset of convection in the 2D cell for $(a)$ adiabatic,
$(b)$ linear and $(c)$ constant sidewall boundary conditions. The most
unstable modes are schematically depicted above each graph with the
corresponding colour. The critical Rayleigh numbers for the current
convection, $Ra_{c}$, are marked with errors.

In this section, we study steady-states in 2D RBC for $Ra\leq 10^{7}$. In what
follows, we refer to flow states as single or multiple solutions connected by
inherent symmetries of the system. For example, the single-roll state (SRS) in
2D can exist in two forms, either circulating clockwise or counterclockwise,
but is considered as a single flow state that is invariant under reflection.
Steady-state solutions of the SRS state have been investigated in laterally
periodic flows with stress-free velocity boundary conditions on the horizontal
walls (Wen et al., 2015, 2020b) and with no-slip BCs (Waleffe et al., 2015;
Sondak et al., 2015; Wen et al., 2020a; Kooloth et al., 2021). Bifurcations
and different flow states have already been studied in laterally unbounded RBC
(Zienicke et al., 1998), in laterally bounded RBC for a cubic domain
(Puigjaner et al., 2008) and a 2D square domain (Venturi et al., 2010). Here
we focus on the onset of convection, the SRS and a vertically stacked double-
roll state (DRS) in two-dimensional RBC for three different sidewall BCs as
shown in figure 1.

###  3.1 Onset of Convection

In RBC, there is a critical Rayleigh number $Ra_{c}$ above which the system
bifurcates from the conduction state to coherent rolls. We calculate $Ra_{c}$
using a linear stability analysis described in more detail in Reiter et al.
(2021b). For adiabatic or linear (conductive) sidewall BCs, the conduction or
base state is characterized by a linear temperature profile in the vertical
direction with zero velocity field and independence from control parameters.
However, for a constant temperature sidewall distribution, a convective flow
is already present. In this case, we perform a steady-state search before
analyzing the local stability around this equilibrium point.

Figure 3 shows the linear growth rates of the four most unstable modes, which
resemble the first four Fourier modes as depicted in the same figure. All
three BCs initially bifurcate from the conduction state to a single roll
state. Adiabatic sidewalls lead to a lower critical Rayleigh number compared
to isothermal sidewalls, which is to be expected (Buell & Catton, 1983). The
onset for the adiabatic sidewall occurs at $Ra_{c}\approx 2.7\times 10^{3}$
which agrees well within our resolution limit with Venturi et al. (2010), who
reports a critical $Ra$ of about $2582$. The onset for the linear SW occurs at
$5.1\times 10^{3}$ and the onset for the constant SW occurs slightly later at
$5.6\times 10^{3}$. This indicates that the interaction of the convective
field - as present for the constant sidewall BC - with the unstable modes is
weak and its influence on the onset is small.

###  3.2 Single-roll (states $\mathcal{S}_{A}^{1}$, $\mathcal{S}_{L}^{1}$,
$\mathcal{S}_{C}^{1}$)

$(a)$

Adiabatic SW

$(b)$

Linear SW

$(c)$

Constant SW

Figure 4:  Single roll state for $(a)$ adiabatic ($Ra=10^{6}$), $(b)$ linear
($Ra=9\times 10^{4}$) and $(c)$ constant ($Ra=10^{6}$) sidewall temperature
boundary conditions. Contours (streamlines) represent the temperature
(velocity) field.

The single roll state (SRS) is arguably the most important state in RBC. It is
the first mode to appear above the conduction state, as we have just seen, and
prevails even up to largest $Ra$ in the form of large-scale circulation (LSC)
on turbulent superstructures (Zhu et al., 2018; Reiter et al., 2021a). The SRS
is stable and time-independent for small $Ra$ but oscillatory, chaotic, or
even completely vanishing for larger $Ra$, as we will show in section 4.3.
Here we analyze its properties before collapse and show that the growth of
secondary corner rolls plays an important role in its destabilization and that
this process can be both suppressed and enhanced by different sidewall
boundary conditions.

Figure 4 shows the temperature and velocity fields of the SRS for different
sidewall BCs. For all three BCs we can identify a large primary roll
circulating counter-clockwise and two secondary corner rolls. The corner rolls
are most pronounced for the linear sidewall BC and the primary roll is nearly
elliptical. The dimensionless heat-flux is expressed in form of the Nusselt
number $Nu\equiv\sqrt{RaPr}F_{f}H/\Delta$ with the heat-flux $F_{f}$ entering
the fluid and the imposed temperature difference $\Delta$. $F_{f}$ can be
defined in different ways, especially in the presence of sidewall heat-fluxes.
Averaging the temperature equation in eq. (1) over time, one obtains

$\displaystyle\divergence\mathbf{F}=0,\quad\mathbf{F}\equiv{\bf
u}\theta-1/\sqrt{RaPr}{\bm{\nabla}}{\theta},$ (14)

from which it follows that the total heat flux must vanish through the
boundaries $S=\delta V$, i.e. $\int_{S}(F\cdot\mathbf{n})dS=0$. For isothermal
sidewall BCs, asymmetric flow states with net nonzero sidewall heat-fluxes are
possible; in this case the heat fluxes through the bottom and top plates would
deviate from each other. However, in the present study, we found that all
sidewall heat fluxes are approximately equal to zero when integrated
vertically and the temperature gradient at the bottom plate is approximately
equal to the temperature gradient at the top plate. Therefore, we define $Nu$
based on the lower (hot) plate at $z=0$:

$Nu\equiv-\frac{1}{A_{+}}\int_{S_{+}}\frac{\partial\theta}{\partial
z}d{S_{+}},$ (15)

with the bottom plate domain $S_{+}$ and its surface area $A_{+}$. The
dimensionless momentum transport is given by the Reynolds number

$Re\equiv\sqrt{Ra/Pr}\sqrt{\langle\mathbf{U}^{2}\rangle_{V}}L,$ (16)

based on total kinetic energy of the mean field velocity $\mathbf{U}$. Here,
$\langle\cdot\rangle_{V}$ denotes a volume average.

$(a)$$(b)$$(c)$Adiabatic SWLinear SWConstant SW Figure 5:  Nusselt number $Nu$
for the single-roll states for $(a)$ adiabatic, $(b)$ linear and $(c)$
constant sidewall temperature boundary conditions.  $(a)$$(b)$$(c)$Adiabatic
SWLinear SWConstant SW Figure 6:  Reynolds number $Re$ for the single roll
states $\mathcal{S}_{A}^{1}$, $\mathcal{S}_{L}^{1}$ , $\mathcal{S}_{C}^{1}$.
$(a)$ adiabatic, $(b)$ linear and $(c)$ constant sidewall temperature boundary
conditions.

In the laminar regime, where the dissipation of velocity and temperature field
is determined by the contributions of the boundary layers, we expect the total
heat and momentum scaling $Nu\sim Ra^{1/4}$ and $Re\sim Ra^{1/2}$ (Grossmann &
Lohse, 2000), respectively. Figure 5 shows that the former scaling shows up
only for a very limited $Ra$ range and only for the adiabatic boundary
conditions. The SRS of the linear sidewall BCs is stable only up to $Ra\leq
10^{5}$, then the corner rolls become strong enough to lead to a collapse of
the SRS. The stability region where the steady-states converge is too small to
observe an unperturbed scaling. On the other hand, for the constant sidewall
boundary conditions, corner roll growth is less dominant. In this case, the
reason why $Nu$ scaling deviates from $1/4$, is that heat entering through the
bottom/top can immediately escape through the sidewalls in the form of a
”short-circuit”, which dominates the lower $Ra$ regime and is the reason why
$Nu$ is relatively large for small $Ra$. For the adiabatic sidewall BC, we
observe $Nu\sim Ra^{0.25}$ for $10^{4}\leq Ra\leq 3\times 10^{5}$, followed by
$Nu\sim Ra^{0.16}$ for $3\times 10^{5}\leq Ra\leq 10^{6}$. Similarly, the
growth of the corner rolls disturbs the convection wind, and $Nu$ deviates
from the ideal $1/4$ scaling. Looking at the $Re$ vs. $Ra$ scaling in figure
6, we find the theoretically predicted scaling of $1/2$ is better represented
in comparison and the different sidewall boundary conditions deviate less
among themselves. This suggests that momentum transport is less affected by
changing sidewall boundary conditions than heat transport.

####  3.2.1 Growth of corner rolls

The SRS is stable up to a certain $Ra$ limit. Above this limit, it may
fluctuate, reverse orientation, or even disappear altogether. This process
occurs at $Ra\approx 10^{6}$ for the adiabatic and constant temperature
sidewall BCs and at $Ra\approx 10^{5}$ for the linear sidewall BC. While up to
this event the dynamic behaviour of the three different sidewall BCs is
qualitatively very similar, from there on it differs. The constant sidewall BC
case shows a time dependence, but remains in the SRS state without changing
its orientation. The adiabatic and linear sidewall BCs, on the other hand,
enter a more chaotic regime of regular and chaotic flow reversals (Xi & Xia,
2007; Sugiyama et al., 2010), some of which are discussed in section 3.3. Of
greatest importance here appears to be the presence and magnification of
secondary corner rolls (CRs).

$(a)$$(b)$diffusion$+$$(c)$buoyancy$+$$=0$$(d)$convection Figure 7:  $(a)$
Steady-state vorticity field, velocity streamlines and corner roll size
$\delta_{CR}$ defined as a distance from the corner to the closest stagnation
point at the plate for $Ra=7\times 10^{5}$ and adiabatic sidewalls, and
vorticity balance contributions according to eq. (17) in the corner roll
domain, i.e., $(b)$ diffusion, $(c)$ buoyancy and $(d)$ convection. The same
contour levers were used for $(b-d)$.

Figure 7 $(a)$ shows the vorticity field and stream-function contour of two-
dimensional RBC with adiabatic sidewalls at $Ra=7\times 10^{5}$. The existence
of two corner vortices is apparent. Here we define their size $\delta_{CR}$
based on the zero crossing, or stagnation point, of the vorticity
$\omega\equiv\partial_{x}u_{z}-\partial_{z}u_{x}$ at the top plate, cf.
Shishkina et al. (2014). To understand the processes involved in the formation
of the corner rolls, we write down the evolution equation for vorticity

$\partial_{t}\omega=\underbrace{-{\bf
u}\cdot{\bm{\nabla}}\omega}_{\text{convection}}+\underbrace{\sqrt{Pr/Ra}{\bm{\nabla}}^{2}\omega}_{\text{diffusion}}+\underbrace{\partial_{x}\theta}_{\text{buoyancy}}.$
(17)

It is evident that for steady-states ($\partial_{t}\omega=0$) there must be an
equilibrium between convection, diffusion and buoyancy forces. The three
corresponding fields are shown in figure 7 $(b-d)$ zoomed in on the corner
roll region. For this particular $Ra$, all three contributions appear to be
significant. We evaluate the size of the corner rolls (figure 8) and analyse
contributions of diffusion, buoyancy, and convection for all $Ra$ (figure 7).
For this purpose, we evaluate the absolute values of the volume averages for
each term in the corner roll region, e.g.,
$\langle|\partial_{x}\theta|\rangle_{V_{CR}}$ represents the strength of the
buoyancy term in the corner roll volume $V_{CR}$, as shown in figure 7 $(c)$.
The constant BC yields a notable exception because multiple corner rolls can
exist. This can be sensed from figure 4 $(c)$. For small $Ra$, the corner roll
are dominant in the lower right and upper left corner, where the LSC detaches
(ejects). For the other two BCs, these rolls are not present. Looking at eq.
(17), we realize that the presence of a horizontal temperature gradient can
lead to the formation of vortex structures. This condition is present for the
constant BCs, e.g., in the lower right corner, where the hot LSC detaches
while the temperature is kept constant at zero, resulting in a (strong)
negative temperature gradient. The two more ”classical” corner rolls first
appear at larger $Ra$, but soon take over in size, as can be seen in figure 8.

$(a)$$(b)$$(c)$Adiabatic SWLinear SWConstant SW Figure 8:  Growth of the
corner roll size $\delta_{CR}$ for $(a)$ adiabatic, $(b)$ linear and $(c)$
constant sidewall temperature boundary conditions. Adiabatic BC show two
distinct regions, a buoyant dominated regime and a regime where convective
influx leads to a more rapid increase. For the constant BC, the corner rolls
appear first in the plume ejecting corner (bottom right and upper left in
figure 4) which is represented by the open symbols in $(c)$, and only for
larger $Ra$ do they appear in the plume impacting region (closed symbols).
$(a)$$(b)$$(c)$Adiabatic SWLinear SWConstant SW Figure 9:  Strength of the
vorticity balance contributions diffusion (black circles), buoyancy (orange
diamonds) and convection (purple pluses) in the corner roll region, according
to eq. (17). $(a)$ adiabatic, $(b)$ linear and $(c)$ constant sidewall
temperature boundary conditions. Adiabatic BC show two distinct regions, a
buoyancy dominated regime and a regime where convective influx leads to a more
rapid increase. For the constant BC, the corner rolls appear first in the
plume ejecting corner (main figure $c$) and only for larger $Ra$ do they
appear in the plume impacting region (inset $c$).

The adiabatic and linear sidewall BCs each yield only two corner rolls. These
are present from the onset of convection and grow until the collapse of the
SRS (figure 8). The main difference between the two is that for the adiabatic
sidewall, the corner rolls initially grow monotonically with respect to $Ra$,
whereas for the linear sidewall BCs, the corner rolls are already considerable
large as soon as the SRS is present. Moreover, they also grow faster with
respect to $Ra$ ($\delta_{CR}\sim Ra^{0.3}$) and soon cover almost $40\%$ of
the width of the cell. Their large initial size combined with faster growth is
the reason for premature SRS instability in linear sidewall BCs. Figure 9
$(b)$ shows that vorticity formation for the entire $Ra$ range is mainly
governed by buoyancy and balanced by diffusion. Assume the hot plumes carry
warm fluid to the upper plate where it meets a cold sidewall, generating
strong lateral gradients in the upper right corner and consequently vorticity,
according to eq. (17).

In the adiabatic case, on the other hand, the sidewall is warmer close to the
corner, which leads to less vorticity generation by lateral temperature
gradients and therefore smaller corner rolls. In the low $Ra$ regime, the
corner rolls of the adiabatic sidewall are also governed by buoyancy, with a
growth of the corner rolls of $\delta_{CR}\sim Ra^{0.21}$ (figure 8 $a$). This
can be understood by dimensional arguments. Assume convection can be neglected
in eq. (17), which is justified from the results in figure 9 $(a)$. Thus we
obtain $\sqrt{Pr/Ra}{\bm{\nabla}}^{2}\omega=\partial_{x}\theta$, or, in terms
of a characteristic temperature $\theta_{CR}$ and a characteristic vorticity
$\Omega_{CR}$, we have
$\nu\frac{\Omega_{CR}}{\delta_{CR}^{2}}\sim\frac{\theta_{CR}}{\delta_{CR}}$,
and thus

$\delta_{CR}\sim\sqrt{\frac{Pr}{Ra}}\frac{\Omega_{CR}}{\theta_{CR}}.$ (18)

The evaluation (not shown here) of the characteristic vorticity in the corner
roll regions by means of their root mean square value unveiled $\Omega\sim
Ra^{0.7}$. Assuming further that the temperature $\theta_{CR}$ is
approximately constant over $Ra$, we obtain $\delta_{CR}\sim Ra^{0.20}$, which
agrees remarkably well with $\delta_{CR}\sim Ra^{0.21}$. Figure 8 $(a)$
discloses a transition at $Ra\approx 3\times 10^{5}$ , above which the corner
roll growth accelerates exhibiting a scaling of $\delta_{CR}\sim Ra^{0.49}$.
Figure 9 $(a)$ indicates that convective processes begin to affect vorticity
generation. Figure 7 $(d)$ reveals a region with strong convective vorticity
current with the same sign as the buoyancy forces, which enhances the
vorticity generation in this region (figure 7 $c$). We interpret that above a
certain $Ra$ the primary roll of the SRS begins to feed the corner rolls until
they become strong enough, eventually leading to the collapse of the SRS
itself. We would like to note that the current analysis describes steady-
states up to $Ra\leq 10^{6}$. An opposite trend was observed for larger $Ra$
by Zhou & Chen (2018), who found a slow shrinkage of the corner rolls that
scales approximately with $\sim Ra^{-0.085}$. It would be interesting to
consolidate these results in future studies.

###  3.3 Double-roll ($\mathcal{S}_{A}^{2}$, $\mathcal{S}_{L}^{2}$)

$(a)$$(b)$

Adiabatic SW

Linear SW

Figure 10:  Double-roll state (DRS) for $(a)$ adiabatic and $(b)$ linear.
Contours (streamlines) represent the temperature (velocity) field.

Having discussed the properties of the SRS state, we proceed to the double-
roll state (DRS) as shown in figure 10. It consists of two vertically stacked
hot and cold circulation cells rotating in opposite directions with an almost
discrete temperature jump in the mid plane. The DRS was not identified as an
equilibrium for the constant sidewall BCs, so we will discuss it exclusively
for the adiabatic and linear sidewall setup. The DRS can coexist with the SRS,
but is generally found at larger $Ra$. Here we have tracked it in the range
$10^{5}\leq Ra<7\times 10^{6}$ for adiabatic and $10^{5}\leq Ra<4\times
10^{6}$ for linear sidewall BCs. This range is consistent with Goldhirsch et
al. (1989) who described a roll-upon-roll state in 2D RBC for $Pr=0.71$ at
$Ra\approx 10^{5}$, but interestingly it was not found for $Pr=6.8$.

$(a)$$(b)$Adiabatic SWLinear SW Figure 11:  Nusselt number $Nu$ for double-
roll states $\mathcal{S}_{A}^{2}$ and $\mathcal{S}_{L}^{2}$. $(a)$ adiabatic
and $(b)$ linear sidewall temperature boundary conditions.
$(a)$$(b)$Adiabatic SWLinear SW Figure 12:  Maximum peak frequency
$f_{\text{max}}$ and average frequency $\overline{f}$ determined from $Nu(t)$
for double-roll states $\mathcal{S}_{A}^{2}$ and $\mathcal{S}_{L}^{2}$ for
$(a)$ adiabatic and $(b)$ linear sidewall temperature boundary conditions.

From figure 11 we see that $Nu$ scales close to $Nu\sim Ra^{1/4}$, which
corresponds to laminar scaling for RBC flows governed by boundary layer
dissipation. Compared to the single-roll state, it is less effective in
transporting heat from wall to wall, as evidenced by an overall smaller $Nu$.
This is actually to be anticipated, since one roll of the DRS can be
conceptually viewed as a half-height, half-temperature gradient RBC system,
implying a $16$ times smaller effective $Ra$. However, this factor most likely
overestimates the difference, since the mid plane velocity is much closer to a
free-slip flow than a no-slip flow and the aspect ratio is two rather than
one. In reality, a DRS has about the same $Nu$ as a SRS with a $6$ times
smaller $Ra$.

The DRS is found to be time-independent (stable) only for the adiabatic
sidewall BCs for $Ra\leq 4\times 10^{5}$. For other $Ra$ it is either
periodically oscillating or chaotic. In figure 12 we show characteristic
frequencies of the DRS obtained by initializing DNS simulation with the
steady-state solutions and evaluating the frequency spectra of $Nu(t)$. The
frequency is presented in free-fall time units. The DRS oscillates with a
frequency of about $0.1$ for $Ra\leq 10^{6}$ for both the adiabatic and linear
setups, i.e., about one cycle every $10$ time units. This cycle corresponds to
about half the circulation time of a cell, i.e., the characteristic velocity
of the circulation is about $0.09\sim 0.11$ and its size is $\approx 2L$.
Thus, the DRS oscillation frequency seems to be initially tied to the
circulation time. When $Ra$ exceeds $10^{6}$, we see the emergence of a more
chaotic behavior. Despite increasing turbulence, the DRS state persists and
does not show transition to a SRS state for $Ra<10^{7}$. In section 4.3 we
will see that for larger $Ra$ the DRS state is eventually replaced by a single
roll LSC again.

The DRS state is not merely an equilibrium solution, but more fundamentally
there is a regime in $Ra$ where the DRS is the preferred flow state to which
all initial states tested in this work tend towards. Starting from random
perturbations, one usually first finds a SRS, which soon goes through a series
of flow reversals and restabilizations until it evolves to the DRS state. This
process is depicted in an SRS-DRS phase space picture in figure 13. The
horizontal axis represents the SRS, and the vertical axis represents the DRS.
This process is qualitatively the same for adiabatic and linear sidewall
boundary conditions. We do not address the flow reversal process, as it is
described in more detail in Xi & Xia (2007); Sugiyama et al. (2010); Castillo-
Castellanos et al. (2016); Zhao et al. (2019), but note that the intermediate
flow fields bear striking resemblance to the proper orthogonal decomposition
modes presented in Podvin & Sergent (2015, 2017). We want to stress that the
transition time is surprisingly long. It can take up to several thousand free-
fall time units for the flow to settle in the DRS state, so it can be missed
if the observation window is too small.

$(a)$$(b)$Adiabatic SWLinear SW Figure 13: Phase space trajectories from a
single-roll ($\mathcal{S}_{A}^{1}$/$\mathcal{S}_{L}^{1}$) to a double-roll
state ($\mathcal{S}_{A}^{2}$/$\mathcal{S}_{L}^{2}$) for $(a)$ adiabatic
sidewall at $Ra=2\times 10^{6}$ and $(b)$ linear sidewall BCs at $Ra=1.5\times
10^{5}$.

##  4 Direct numerical simulations

In addition to the steady-state analysis, we performed a series of DNS of RBC
for 2D in a square and 3D in a cylinder with $\Gamma=1$ and $Pr=1$, covering
$Ra$ from the onset of convection to $4.64\times 10^{10}$ and $10^{9}$,
respectively. The highest $Ra$ in 2D was simulated on a $1024^{2}$ grid with
at least $15$ grid points in the thermal boundary layer and performed for
several thousand free-fall time units, ensuring adequate spatial resolution
and temporal convergence. The largest simulation for the cylindrical setup was
performed on a $N_{r}\times N_{\varphi}\times N_{z}=128\times 256\times 320$
grid, with about $10$ points inside the thermal and viscous boundary layers
and the averaging statistics were collected for at least $600$ free-fall time
units.

###  4.1 Vertical temperature profiles

$(a)$$(b)$$(c)$

2D Box

Adiabatic SWLinear SWConstant SW$(d)$$(e)$$(f)$

3D Cylinder

Adiabatic SWLinear SWConstant SW Figure 14:  Mean temperature profile for
cases with $(a,d)$ adiabatic, $(b,e)$ linear and $(c,f)$ constant sidewall
boundary conditions for ($a$-$c$) 2D box and ($d$-$f$) cylinder.

Figure 14 shows the horizontally averaged temperature profiles
$\langle\theta\rangle_{A}$ for all conducted simulations. We first remark the
similarity between 2D and 3D. For example, both show the feature of a weakly
stabilizing positive temperature gradient in the mid plane for small $Ra$ and
adiabatic boundary conditions (figures 14 a,d). This phenomenon is often found
in the interior of the bulk (Tilgner et al., 1993; Brown & Ahlers, 2007; Wan
et al., 2019) and is caused by the thermal signature of the LSC. As the
thermal plume of the LSC climbs up along the sidewall, it penetrates deeper
into the bulk, thus hot (cold) plumes carry their signature into the top
(bottom) part of the cell, which can result in a slightly positive temperature
gradient in the center of the bulk.

Another important detail is the apparent non-monotonicity of the profiles in
the intermediate $Ra$ range, which is most pronounced for the linear sidewall
BCs (figure 14 b,e) and also occurs for the 2D adiabatic BCs. The temperature
profiles initially drop sharply and then level of at about a quarter of the
cell height before dropping sharply again in the cell center. This behaviour
was also observed in Stevens et al. (2014). These profiles are reminiscent of
the DRS state (see section 3.3) and indeed caused by transitions in the flow
structures, which we analyse in section 4.3 in more detail. Finally, all
simulations for larger $Ra$ show the classical RBC profile with steep
temperature gradients at the bottom and top plates and a well-mixed
homogeneous bulk.

###  4.2 Vertical sidewall heat flux profiles

$(a)$$(b)$

2D Box

Linear SWConstant SW$(c)$$(d)$

3D Cylinder

Linear SWConstant SW Figure 15:  Comparison of the lateral sidewall heat flux
$Nu_{sw}$ for cases ($a,c$) linear and ($b,d$) constant sidewall boundary
conditions in ($a,b$) 2D box and ($c,d$) cylinder.

Next we analyse the horizontal heat flux through the vertical sidewall
$Nu_{sw}$ which is more elaborately defined in the appendix A. This is shown
in figure 15 for the linear and constant BCs, while the sidewall heat flux of
the adiabatic BC is obviously zero. The linear and constant BCs show two
opposite trends. The constant setup has the largest temperature gradients for
small $Ra$ and almost vanishing gradients for large $Ra$. This can be
understood from the temperature profiles in figure 14 $(c,f)$. As $Ra$
increases, the bulk is more efficiently mixed and the temperature distribution
becomes nearly constant, hence the temperature in the cell becomes more
similar to the sidewall temperature imposed by the BCs. On the other hand, the
linear sidewall BC corresponds exactly to the temperature profile before the
onset of convection and from then on its contrast increases more and more,
which is reflected in the relatively strong vertical temperature gradients for
large $Ra$. However, all profiles are symmetrical around the center and
consequently, although heat flows in and out locally, there is no net heat
flux through the vertical sidewalls. This is supported by the fact that in our
simulations the temperature gradients at the top and bottom plates were nearly
equal, linked by the heat flux balance

$Nu_{c}-Nu_{h}+\zeta\langle Nu_{sw}\rangle_{z}=0$ (19)

with $\zeta=\frac{1}{\Gamma}$ for the 2D box and $\zeta=\frac{4}{\Gamma}$ for
the cylindrical setup (see appendix A). Lastly, we detect at least two
transitions in $Nu_{sw}$ for the linear sidewall BCs (figure 15 $a,c$). These
are consistent with the transitions in the temperature profiles discussed in
the previous section and are elucidated in more detail in the following.

###  4.3 Mode analysis

It is generally difficult to compare the dynamics of flows in different,
possibly even turbulent, states without restricting the underlying state
space. Therefore, in this section we analyze the DNS results by projecting
each snapshot onto four distinct modes and evaluate time averages and standard
deviations.

Starting with the 2D simulations, a common choice for the mode are the first
four Fourier modes, see e.g. Petschel et al. (2011) and (Wagner & Shishkina,
2013), i.e.

$\displaystyle u_{x}^{m,k}$ $\displaystyle=-\sin(\pi mx/L)\cos(\pi kz/H),$
$\displaystyle u_{z}^{m,k}$ $\displaystyle=\cos(\pi mx/L)\sin(\pi kz/H).$ (20)

For the cylinder, the choice of modes is less obvious. In this work, we follow
Shishkina (2021) and use a combination of Fourier modes in $z$ and $\varphi$
direction and Bessel functions of the first kind $J_{n}$ of order $n$ in $r$
for the radial velocity component $u_{r}$ and the vertical velocity component
$u_{z}$. The first two (non-axisymmetric) modes are

$\displaystyle u_{r}^{1,k}$ $\displaystyle=J_{0}(\alpha_{0}r/R)\cos(\pi
kz/H)e^{i\varphi},$ $\displaystyle u_{z}^{1,k}$
$\displaystyle=J_{1}(\alpha_{1}r/R)\sin(\pi kz/H)e^{i\varphi},$ (21)

and the axisymmetric modes are

$\displaystyle u_{r}^{2,k}$ $\displaystyle=J_{1}(\alpha_{1}r/R)\cos(\pi
kz/H),$ $\displaystyle u_{z}^{2,k}$
$\displaystyle=-J_{0}(\alpha_{0}r/R)\sin(\pi kz/H),$ (22)

where $\alpha_{n}$ is the first positive root of the Bessel function $J_{n}$
for Dirichlet boundary conditions on the sidewall ($u_{r}$) and the $k$-th
positive root of the derivative of the Bessel function $J_{n}^{\prime}$ for
Neumann boundary conditions $(u_{z})$. The non-axisymmetric modes are complex-
valued to account for different possible azimuthal orientations. Ultimately,
however, we are only interested in the energy content and not the orientation
of the modes, so we evaluate their magnitude. We note further, that a vertical
slice through the cylindrical modes is very similar to the first four 2D
Fourier modes, albeit with a slightly different dependence in the radial
direction. For this reason, we use the same notation for the cylindrical modes
as for the Fourier modes in 2D. More precisely, we have
$F_{1}\equiv(u_{r}^{1,1},u_{z}^{1,1})$,
$F_{2}^{=}\equiv(u_{r}^{1,2},u_{z}^{1,2})$,
$F_{2}^{\parallel}\equiv(u_{r}^{2,1},u_{z}^{2,1})$ and
$F_{4}\equiv(u_{r}^{2,2},u_{z}^{2,2})$. Having defined the modes, we project
the velocity field ${\bf u}$ of several snapshots onto a mode ${\bf u}^{m}$
and evaluate the energy content $\mathcal{P}$ of each mode according to

$\displaystyle\mathcal{P}\equiv\frac{\int_{V}{\bf u}{\bf
u}^{m}dV}{\int_{V}{\bf u}^{m}{\bf u}^{m}dV},$ (23)

and analyse the time average and standard deviation of $\mathcal{P}$.

$\mathcal{F}_{1}$$\mathcal{F}_{2}^{=}$$\mathcal{F}_{2}^{\parallel}$$\mathcal{F}_{4}$$(a)$$(b)$$(c)$Adiabatic
SWLinear SWConstant SW Figure 16: Energy and standard deviation of the
projection of flow field snapshots onto the modes defined by eq. (20) for the
2D box and $(a)$ adiabatic, $(b)$ linear and $(c)$ constant sidewall
temperature boundary condition for the 2D box. Below: Streamlines, coloured by
vertical velocity, of the modes $\mathcal{F}_{1}$, $\mathcal{F}_{2}^{=}$,
$\mathcal{F}_{2}^{\parallel}$ and $\mathcal{F}_{4}$.

The energy of the individual Fourier mode for the 2D box is shown in figure
16. Above the onset of convection, only the first Fourier mode (single-roll)
contains a considerable amount of energy. Because of its similarity to the
SRS, this mode will be referred to as the SRS-mode. Following the stable SRS,
we find for adiabatic and linear sidewall BCs a flow regime that changes from
the SRS to a roll-upon-roll second Fourier mode
($\mathcal{F}_{2}^{\parallel}$) state. This state embodies the DRS state,
which we discussed in section 3.3. The $F_{2}^{=}$ regime, or DRS regime, is
found in the range $10^{6}<Ra\leq 10^{7}$ for an adiabatic sidewall and
$10^{5}\leq Ra\leq 10^{7}$ for a linear sidewall BC. In contrast, the DRS
regime is absent for a constant sidewall BC. As a reminder, this state could
not be found as an equilibrium solution for the constant sidewall boundary
condition either, which is in line with its absence in DNS. The next regime
can be regarded as a weakly chaotic SRS regime, with the SRS mode again
dominating but being transient and a substantial amount of energy is contained
in the $F_{4}$ (4-roll) mode, indicative of dynamically active corner rolls.
Finally, above $Ra\approx 10^{9}$ there exists another surprisingly sharp
transition. This regime is different from the others as now all Fourier modes
contain a significant amount of energy and exhibit strong fluctuations. An
inspection of the flow fields revealed an abundance of small-scale plumes and
strong turbulent dynamics. Most remarkably, in this regime all three sidewall
BCs show a very similar mode signature, i.e., they become increasingly alike,
or in other words, RBC becomes insensitive to sidewall BCs for large $Ra$.

$(a)$$(b)$$(c)$Adiabatic SWLinear SWConstant
SW$\mathcal{F}_{1}$$\mathcal{F}_{2}^{=}$$\mathcal{F}_{2}^{\parallel}$$\mathcal{F}_{4}$
Figure 17: Energy and standard deviation of the projection of flow field
snapshots onto the modes defined by eq. (22) and (21) for $(a)$ adiabatic,
$(b)$ linear and $(c)$ constant sidewall temperature boundary condition for
the cylinder. Below: Streamlines, coloured by vertical velocity, of the modes
$\mathcal{F}_{1}$, $\mathcal{F}_{2}^{=}$, $\mathcal{F}_{2}^{\parallel}$ and
$\mathcal{F}_{4}$.

Moving on to the mode analysis for the cylindrical setup, shown in figure 17,
we see a very similar picture as for the 2D box with some noticeable
differences. First, for the constant BC setup we note that the onset of
convection is significantly later than in the 2D case, while the other two
setups show a closer similarity with the 2D case. The cylindrical setup might
be more sensitive to the BCs of the sidewalls in general, since the ratio of
sidewall area to cell volume ratio is larger than in the 2D box and therefore
the sidewall temperature likely has a larger impact on the interior.

Another difference between the cylindrical and 2D box setup is, that the
adiabatic setup does not show a transition to a regime with a vanishing SRS;
rather, the SRS mode is the most dominant mode over all $Ra$. In contrast, the
linear sidewall BC possess a striking similarity to the observations in 2D.
Above $Ra\approx 10^{5}$ it undertakes a transition from a SRS-dominated
regime to a $F_{4}$-dominated regime. The $F_{4}$-mode is axissymmetric and
has a double-donut, or double-toroidal shape. Similar flow states were found
in a bifurcation analysis by Puigjaner et al. (2008) in a cubic domain with
the same lateral boundary conditions. Here, its existence range extends over
$10^{5}\leq Ra\leq 10^{8}$. The double-donut state can be considered as the
counterpart of the DRS state in 2D RBC, although we see that it outlasts its
2D analog by about a decade in $Ra$. At the highest $Ra$ available, the SRS
again dominates for all BC configurations considered, although the amount of
energy and the strength of the fluctuations are somewhat different for the
different BCs. At this points, we can only conjecture from their trend and our
findings in 2D that their deviation will decrease for even larger $Ra$ in the
high-turbulence/high-$Ra$ regime.

We conclude that there exist at least five different flow regimes: conduction
state, stable SRS, DRS (or double-donut state in the cylindrical setup),
weakly chaotic SRS and highly turbulent state. We find the constant isothermal
sidewall generally enhances the SRS dominance, while a linear isothermal
sidewall BC suppresses the SRS in the mid $Ra$ regime and induces the DRS or
double-donut state. Moreover, although we find strong differences in the flow
dynamics in the small to medium $Ra$ range, but these differences eventually
disappear and the system becomes increasingly insensitive to the type of
sidewall BC at high $Ra$.

###  4.4 Heat transport

Lastly, the global heat transport is discussed. The results are shown in
figure 18. For the 2D setup, we include the results from the steady-state
analysis from the first part of this study. Here, we find a very good
agreement between $Nu$ of the DNS and steady-states for the SRS mode as well
as for the DRS state for adiabatic sidewalls. However, the DRS state for
linear sidewalls shows slightly larger $Nu$ in the DNS. This is because the
DRS state is an unstable equilibrium solution that can oscillate strongly,
which apparently enhances heat transport properties.

We find that $Nu$ degrades strongly when switching from a SRS- to a DRS-
dominated regime at $Ra\approx 10^{5}$ (linear) and $Ra\approx 10^{6}$
(adiabatic) for the 2D domains (figure 18$a$). In contrast, this does not
occur for the cylindrical setup as it transitions from the SRS to the double-
toroidal state (figure 18$b$). In fact, this flow transition is hardly
observed in the evolution of heat transport.

In the high $Ra$ regime, the heat transport in the the cylindrical setup is
found to be more efficient than in the 2D setup, with about $30\%$ larger
$Nu$. This agrees well with the observations of van der Poel et al. (2013).
Both setups show $Nu\sim Ra^{0.285}$ scaling at the largest studied $Ra$. We
also observe that $Nu$ becomes independent of the choice of sidewall BCs for
high $Ra$. This agrees with Stevens et al. (2014), at least when the sidewall
temperature is equal to the arithmetic mean of bottom and top plate
temperature. If this condition is violated, Stevens et al. (2014) has shown
that $Nu$ differences will exist even for high $Ra$. This indicates that the
effects of an imperfectly insulated sidewall tend to be small in experiments
when the mean temperature of the sidewall is well controlled.

$(a)$$(b)$2D Box3D Cylinder Figure 18:  Nusselt number Nu for cases with
different sidewall boundary conditions in $(a)$ 2D simulations, $(b)$ 3D
simulations. For comparison, open symbols shows heat transport in a periodic
2D domain with $\Gamma=2$ by Johnston & Doering (2009) $(a)$ and for
cylindrical setup with adiabatic sidewalls, $\Gamma=1$ and $Pr=0.7$ conducted
by Emran & Schumacher (2012) $(b)$. Dashed lines in $(a)$ show the results
from the steady-state analysis.

###  4.5 Prandtl number dependence

The previous analysis focused on fluids with $Pr=1$, but thermal convection is
relevant in nature in a wide variety of fluids and many experiments are
conducted in water ($Pr\approx 4$) or in liquid metals ($Pr\ll 1$) (Zwirner et
al., 2020). Therefore, we now explore the $Pr$ parameter space with $Pr=0.1,1$
and $10$ for $Ra$ up to $10^{9}$ in the 2D RBC setup.

The Nusselt number is shown in figure 19. We observe a collapse of all data
points for all studied boundary conditions at large $Ra$. However, the
collapse for large $Pr$ is achieved earlier, at $Ra\gtrapprox 10^{7}$, whereas
the differences between $Pr=1.0$ and $Pr=0.1$ are small. Both indicate heat
transport invariance for $Ra\gtrapprox 10^{8}$. This suggests that the size of
the thermal boundary layer $\lambda_{\theta}$ plays a crucial role. For small
$Pr$ we expect larger thermal boundary layers, which extend further into the
bulk and thus have a stronger influence on the system. As $\lambda_{\theta}$
gets smaller, the coupling between the sidewall and bulk disappears, and so do
the differences in heat transport. And although our results show a small
$Pr$-dependence, the main message remains. Experiments with very high $Ra$ are
not affected by different thermal sidewall BCs, regardless of whether they are
performed in a low $Pr$ or high $Pr$ medium.

$(a)$$(b)$$(c)$$Pr=0.1$$Pr=1$$Pr=10$ Figure 19:  Nusselt number $Nu$ for $(a)$
$Pr=0.1$, $(b)$ $Pr=1$ and $(c)$ $Pr=10$ in 2D RBC with different thermal
sidewall BCs.

##  5 Conclusion

We have investigated the influence of three different lateral thermal boundary
conditions, i.e., adiabatic, linearly distributed in the vertical direction
and constant (isothermal) ones, on heat transport and flow states in two- and
three-dimensional Rayleigh-Bénard convection (RBC) using direct numerical
simulation and steady-state analysis. The steady-state analysis is based on an
adjoint-descent method (Farazmand, 2016). We found superior convergence chance
in the laminar and weakly laminar regime compared to Newton’s method, but did
not achieve convergence at larger $Ra$. Further studies on the proper boundary
conditions, the choice of the energy norm and or a combination with Newton’s
method are needed to further explore the potential of the method in the study
of convective flows.

Investigation of the stability of the single-roll state (SRS) revealed that a
linear temperature distribution at the sidewall leads to a premature collapse
of the SRS compared to adiabatic BCs. In contrast, the stability of the SRS
was enhanced by the introduction of constant temperature sidewall BCs. We find
that in 2D and for linear and adiabatic sidewall BCs, the collapse of the SRS
is followed by a regime in which the preferred flow state is a double-roll
state (DRS), where one roll is located on top of the other. The DRS can be
found for adiabatic and linear BCs in the regime $10^{6}<Ra\leq 10^{7}$ and
$10^{5}\leq Ra\leq 10^{7}$, respectively, and is associated with suppressed
heat transport. The DRS can be stable, it can oscillate periodically with a
frequency of $\approx 0.1$ free-fall time unit, or it can be chaotic for
larger $Ra$. In 3D cylindrical simulations, a similar flow transition occurs.
Imposing linear sidewall BCs leads to the emergence of a double-toroidal
structure, that prevails over a wide range of $Ra$, i.e., $10^{5}\leq Ra\leq
10^{8}$. Unlike in 2D, the double-toroidal structure does not lead to a heat
transport recession.

We confirmed that the collapse of the SRS in 2D RBC is strongly related to the
enlarging of corner rolls. Examining the setup with adiabatic sidewalls, there
seem to be two regimes with distinct corner roll growth rates. For small $Ra$,
the vorticity balance is dominated purely by diffusion and buoyancy in the
form of lateral temperature gradients. In this regime, the size of the corner
roll $\delta_{CR}$ grows as $\delta_{CR}\sim Ra^{0.21}$, which is consistent
with dimensional analysis. For larger $Ra$, the convective flux starts to be
of significance and the growth of the corner roll accelerates to
$\delta_{CR}\sim Ra^{0.49}$ before the SRS finally collapses and slowly
transforms to the DRS state, undergoing several cycles of flow reversals and
restabilization.

Analysis of global heat transport and the flow dynamics have shown that for
$Ra\leq 10^{8}$ there are significant differences between the various sidewall
BCs. However, for larger $Ra$ and for various $Pr$ these differences disappear
and the different sidewall BCs become globally - in terms of their integral
quantities - and dynamically similar. In this context, Verzicco & Sreenivasan
(2008) and Johnston & Doering (2009) showed that regardless of imposition of
fixed temperature or fixed heat flux at the bottom/top plates, high $Ra$ show
similar heat transport. Thus, together with our results, we can conclude that
the effects of different boundary conditions, at the sidewalls or at the
top/bottom plates, are limited for experiments with high $Ra$. However, there
are exceptions. For example, when the sidewall temperature differs from the
mean fluid temperature, larger $Nu$ differences can occur (Stevens et al.,
2014). Thus, in experiments at high Rayleigh numbers, it appears to be more
important to control the mean sidewall temperature than to ensure perfectly
insulating conditions. However, close to the onset of convection, the sidewall
thermal boundary conditions significantly influence the flow organization and
heat transport in the system.

## Acknowledgement

This work was supported by the Deutsche Forschungsgemeinschaft, grants
Sh405/10, Sh405/8, Sh405/7 (SPP 1881 Turbulent “Superstructures”). The authors
also acknowledge Leibniz Supercomputing Centre (LRZ) for providing computing
time.

## Declaration of Interests

The authors report no conflict of interest.

##  Appendix A Heat flux

The temperature equation for an incompressible fluid in dimensional units is

$\displaystyle{\partial}{\theta}/{\partial}t+{\bm{\nabla}}\cdot{(\bf
u\theta)}$ $\displaystyle=\kappa{\bm{\nabla}}^{2}{\theta}.$ (24)

Averaging equation $\eqref{eq:T}$ over time yields the following relations for
the heat flux $\mathbf{F}$:

$\displaystyle\divergence\mathbf{F}=0,\quad\mathbf{F}\equiv{\bf
u}\theta-\kappa{\bm{\nabla}}{\theta}.$ (25)

Using the divergence theorem we obtain

$\displaystyle\int_{S}\mathbf{F}\cdot\mathbf{n}dS=0,$ (26)

which states that the net heat flux through the walls must be zero. Expressing
the heat fluxes by the Nusselt number and decomposing the contribution of the
surface integral into those for a lower plate heat flux $Nu_{h}$, for an upper
plate heat flux $Nu_{c}$ and for a side wall heat flux $Nu_{sw}$, we write

$\displaystyle Nu_{c}-Nu_{h}+\zeta\langle Nu_{sw}\rangle_{z}=0,$ (27)

where $\langle\cdot\rangle_{z}$ denotes a vertical mean and $\zeta$ a
geometric factor defining the ratio of the sidewall surface to the bottom/top
plate surface, which is $\zeta=1/\Gamma$ for the 2D box and $\zeta=4/\Gamma$
for the cylindrical setup. Note that the lateral heat flux $Nu_{sw}$ is
$z$-dependent as it was shown in section 4.2. For the 2D box this is

$\displaystyle Nu_{sw}=\frac{H}{\Delta}\left[\frac{\partial\theta}{\partial
x}\evaluated{}_{x=L}-\frac{\partial\theta}{\partial
x}\evaluated{}_{x=0}\right]$ (28)

and for the 3D cylinder setup it is

$\displaystyle
Nu_{sw}=\frac{H}{2\pi\Delta}\int_{0}^{2\pi}\frac{\partial\theta}{\partial
r}\evaluated{}_{r=R}d\varphi.$ (29)

##  Appendix B Thermal dissipation rate

Multiplying equation $\eqref{eq:T}$ with $\theta$ and averaging over time
yields

$\displaystyle\frac{1}{2}{\partial}_{t}\theta^{2}+\frac{1}{2}{\bm{\nabla}}\cdot{({\bf
u}\theta^{2})}$ $\displaystyle=\kappa\theta{\bm{\nabla}}^{2}{\theta}.$ (30)

Taking a time and volume average of $\eqref{eq:T2}$, the time derivative and
the convective part (for impenetrable walls) vanish and using the relation
$({\bm{\nabla}}\theta)^{2}=\divergence{(\theta{\bm{\nabla}}\theta)}-\theta{\bm{\nabla}}^{2}\theta$
we obtain

$\displaystyle\kappa\int_{V}\overline{({\bm{\nabla}}\theta)^{2}}dV=\kappa\int_{V}\divergence{(\overline{\theta{\bm{\nabla}}\theta})}dV,$
(31)

where an overbar denotes a time average and
$\varepsilon_{\theta}=\kappa({\bm{\nabla}}\theta)^{2}$ is known as the thermal
dissipation rate. Using the divergence theorem once more, we find the relation
between the total thermal dissipation rate and the wall heat fluxes

$\displaystyle\int_{V}\overline{\varepsilon_{\theta}}dV=\kappa\int_{S}(\overline{\theta{\bm{\nabla}}\theta})\cdot\mathbf{n}dS.$
(32)

For clarification, writing eq. $\eqref{eq:eth}$ more explicitly and only for
2D Cartesian coordinates, we get

$\displaystyle\langle\overline{\varepsilon_{\theta}}\rangle_{V}$
$\displaystyle=\frac{\kappa}{V}\left(L\left[\langle\overline{\theta\partial_{z}\theta}\rangle_{x}\right]_{z=0}^{z=H}+H\left[\langle\overline{\theta\partial_{x}\theta}\rangle_{z}\right]_{x=0}^{x=L}\right),$
(33)

with the horizontal and vertical average $\langle\cdot\rangle_{x}$ and
$\langle\cdot\rangle_{z}$, respectively. In RBC, the temperatures of the upper
and lower plates are spatially homogeneous, i.e. $\theta_{h}=\frac{\Delta}{2}$
and $\theta_{c}=-\frac{\Delta}{2}$, and assuming that the vertical wall fluxes
are equal (which is not necessarily the case for non-adiabatic sidewalls, but
has been shown to be true in all our simulations), i.e.,
$\partial_{z}\theta_{c}=\partial_{z}\theta_{h}$, then

$\displaystyle\langle\overline{\varepsilon_{\theta}}\rangle_{V}$
$\displaystyle=\frac{\kappa}{V}\left(-L\Delta\langle\partial_{z}\theta_{h}\rangle_{x}+H\left[\langle\overline{\theta\partial_{x}\theta}\rangle_{z}\right]_{x=0}^{x=L}\right),$
$\displaystyle\langle\overline{\varepsilon_{\theta}}\rangle_{V}$
$\displaystyle=\frac{\kappa\Delta^{2}}{H^{2}}Nu+\frac{\kappa}{L}\left[\langle\overline{\theta\partial_{x}\theta}\rangle_{z}\right]_{x=0}^{x=L}.$
(34)

This results in
$\langle\overline{\varepsilon_{\theta}}\rangle_{V}=\frac{\kappa\Delta^{2}}{H^{2}}Nu$
for adiabatic sidewalls or for zero temperature sidewalls, but adds an
additional term to the $\varepsilon_{\theta}-Nu$ relation otherwise. A
comparison of $Nu$ and $\varepsilon_{\theta}$ is shown in figure 20. The
virtual discontinuity of $\varepsilon_{\theta}$ for the linear sidewall
temperature reflects the reordering of the flow structures as explained in the
main part of this study, but surprisingly $Nu$ shows a rather smooth change in
this regime.

Figure 20:  Comparison of $Nu$ (closed symbols) and thermal dissipation rate
$\varepsilon_{\theta}$ (open symbols) in the 2D box. The connection between
thermal dissipation and $Nu$ is given in equation (34).

##  Appendix C Adjoint descent

###  C.1 Derivation

Following Farazmand (2016), we define the right-hand side of the Navier-Stokes
equations as the vector $\mathbf{F_{0}}$, i.e.

$\mathbf{F_{0}}({\bf q})=\begin{pmatrix}-{\bf u}\cdot{\bm{\nabla}}{\bf
u}-{\bm{\nabla}}p+\nu{\bm{\nabla}}^{2}{\bf
u}+\mathbf{\mathbf{e}}_{z}\theta\\\\[3.0pt] -{\bf
u}\cdot{\bm{\nabla}}\theta+\kappa{\bm{\nabla}}^{2}\theta\\\\[3.0pt]
{\bm{\nabla}}\cdot{\bf u}\end{pmatrix}.$ (35)

The functional Gateaux derivative $\delta F({\bf u},{\bf
u}^{\prime})\coloneqq\lim\limits_{\varepsilon\to 0}\frac{F({\bf
u}+\varepsilon{\bf u}^{\prime})-F({\bf u})}{\varepsilon}$ of equation (35) is

$\delta F({\bf q},{\bf q}^{\prime})=\begin{pmatrix}-{\bf
u}^{\prime}\cdot{\bm{\nabla}}{\bf u}-{\bf u}\cdot{\bm{\nabla}}{\bf
u}^{\prime}-{\bm{\nabla}}p^{\prime}+\nu{\bm{\nabla}}^{2}{\bf
u}^{\prime}+\mathbf{\mathbf{e}}_{z}\theta^{\prime}\\\\[3.0pt] -{\bf
u}^{\prime}\cdot{\bm{\nabla}}\theta-{\bf
u}\cdot{\bm{\nabla}}\theta^{\prime}+\kappa{\bm{\nabla}}^{2}\theta^{\prime}\\\\[3.0pt]
{\bm{\nabla}}\cdot{\bf u}^{\prime}\end{pmatrix}.$ (36)

We want to find the adjoint operator $\delta F^{\dagger}$ of equation (36)
with respect to the inner-product

$\langle{\bf q},{\bf
q}^{\prime}\rangle_{\mathcal{A}}=\int_{\mathcal{D}}\left({\bf
q}\cdot\mathcal{A}{\bf q}^{\prime}\right)\text{d}\bf x.$ (37)

The adjoint $\delta F$ of equation (36) with respect to the inner product
(37), with $\tilde{{\bf q}}\equiv\mathcal{A}{\bf q}$, is derived as follows

$\displaystyle\langle\delta F({\bf q},{\bf q}^{\prime}),\tilde{{\bf
q}}^{\prime\prime}\rangle_{\mathcal{A}}$ $\displaystyle=$
$\displaystyle=\int_{V}\begin{pmatrix}-{\bf u}^{\prime}\cdot{\bm{\nabla}}{\bf
u}-{\bf u}\cdot{\bm{\nabla}}{\bf
u}^{\prime}-{\bm{\nabla}}p^{\prime}+\nu{\bm{\nabla}}^{2}{\bf
u}^{\prime}+\mathbf{\mathbf{e}}_{z}\theta^{\prime}\\\\[3.0pt] -{\bf
u}^{\prime}\cdot{\bm{\nabla}}\theta-{\bf
u}\cdot{\bm{\nabla}}\theta^{\prime}+\kappa{\bm{\nabla}}^{2}\theta^{\prime}\\\\[3.0pt]
{\bm{\nabla}}\cdot{\bf u}^{\prime}\end{pmatrix}\begin{pmatrix}\tilde{{\bf
u}}^{\prime\prime}\\\\[3.0pt] \tilde{\theta}^{\prime\prime}\\\\[3.0pt]
\tilde{p}^{\prime\prime}\end{pmatrix}\text{d}\bf x$
$\displaystyle=\int_{V}\begin{pmatrix}\left({\bm{\nabla}}\tilde{{\bf
u}}^{\prime\prime}+{\bm{\nabla}}\tilde{{\bf
u}}^{\prime\prime\text{T}}\right){\bf
u}+\theta{\bm{\nabla}}\tilde{\theta}^{\prime\prime}-{\bm{\nabla}}\tilde{p}^{\prime\prime}+\nu{\bm{\nabla}}^{2}\tilde{{\bf
u}}^{\prime\prime}\\\\[3.0pt] {\bf
u}\cdot{\bm{\nabla}}\tilde{\theta}^{\prime\prime}+\nu{\bm{\nabla}}^{2}\tilde{\theta}^{\prime\prime}+\mathbf{\mathbf{e}}_{z}\cdot\tilde{{\bf
u}}^{\prime\prime}\\\\[3.0pt] {\bm{\nabla}}\cdot\tilde{{\bf
u}}^{\prime\prime}\end{pmatrix}\begin{pmatrix}{\bf u}^{\prime}\\\\[3.0pt]
\theta^{\prime}\\\\[3.0pt] p^{\prime}\end{pmatrix}\text{d}\bf x$
$\displaystyle=\langle{\bf q}^{\prime},\delta F^{\dagger}({\bf q},\tilde{{\bf
q}}^{\prime\prime})\rangle_{\mathcal{A}},$ (38)

where the second line follows from integration by parts. Here we have
refrained from writing the boundary terms that follow from the integration by
parts step, since they can be eliminated by choosing the boundary conditions
on $\tilde{{\bf q}}^{\prime\prime}$ as discussed in section 2.3.

###  C.2 Choice of the norm

As mentioned in Farazmand (2016), the most obvious choice for the norm is the
$\text{L}^{2}$ norm, i.e. $\mathcal{A}=I$, where $I$ is the identity operator.
However, this norm is rather stiff and leads to restrictive small time steps.
As an alternative, Farazmand (2016) uses a norm related to the Laplacian,
which effectively smooths the $\tilde{{\bf q}}^{\prime\prime}$ field. Here we
use a similar norm based on the inversed Laplacian, i.e.
$\mathcal{A}=(I-\alpha{\bm{\nabla}}^{2})^{-1}$,

$\langle{\bf q},{\bf
q}^{\prime}\rangle_{{\bm{\nabla}}^{-2}}=\int_{V}\left({\bf
q}\cdot\mathcal{A}{\bf q}^{\prime}\right)\text{d}\bf x=\int_{V}\left({\bf
q}\cdot\tilde{{\bf q}}^{\prime}\right)\text{d}\bf x$ (39)

where $a$ is a positive constant. Then, $\tilde{{\bf q}}^{\prime}$ is obtained
as the solution of the Helmholtz equation

$(I-\alpha{\bm{\nabla}}^{2})\tilde{{\bf q}}^{\prime}={\bf q}^{\prime},$ (40)

which points out the smoothing property of this norm. In practice, we choose
$\alpha=1$. The choice of the operator for the energy norm is somewhat
arbitrary, but this peculiar choice leads to improved numerical stability
properties. Note that the operator $\mathcal{A}$ should be positive definite
and should commute with the divergence operator, i.e.
$\mathcal{A}({\bm{\nabla}}\cdot{\bf u})={\bm{\nabla}}\cdot\mathcal{A}{\bf u}$.

## References

  * Ahlers (2000) Ahlers, G. 2000 Effect of sidewall conductance on heat-transport measurements for turbulent Rayleigh–Bénard convection. Phys. Rev. E 63, 015303.
  * Ahlers et al. (2009a) Ahlers, G., Funfschilling, D. & Bodenschatz, E. 2009a Transitions in heat transport by turbulent convection at Rayleigh numbers up to $10^{15}$. New J. Phys. 11, 123001.
  * Ahlers et al. (2009b) Ahlers, G., Grossmann, S. & Lohse, D. 2009b Heat transfer and large scale dynamics in turbulent Rayleigh–Bénard convection. Rev. Mod. Phys. 81, 503–537.
  * Ahlers et al. (2012) Ahlers, G., He, X., Funfschilling, D. & Bodenschatz, E. 2012 Heat transport by turbulent Rayleigh–Bénard convection for $Pr\sim 0.8$ and $3\times 10^{12}\lesssim~{}Ra\lesssim 10^{15}$: Aspect ratio $\Gamma=0.50$. New J. Phys. 14, 103012.
  * Bodenschatz et al. (2000) Bodenschatz, E., Pesch, W. & Ahlers, G. 2000 Recent developments in Rayleigh–Bénard convection. Annu. Rev. Fluid Mech. 32, 709–778.
  * Brown & Ahlers (2007) Brown, E. & Ahlers, G. 2007 Temperature gradients, and search for non-Boussinesq effects, in the interior of turbulent Rayleigh–Bénard convection. Eur. Phys. Lett. 80, 14001\.
  * de Bruyn et al. (1996) de Bruyn, J. R., Bodenschatz, E., Morris, S. W., Trainoff, S. P., Hu, Y., Cannell, D. S. & Ahlers, G. 1996 Apparatus for the study of Rayleigh–Bénard convection in gases under pressure. Rev. Sci. Instrum. 67 (6), 2043–2067.
  * Buell & Catton (1983) Buell, J. C. & Catton, I. 1983 The effect of wall conduction on the stability of a fluid in a right circular cylinder heated from below. J. Heat Transfer 105, 255–260.
  * Busse (1967) Busse, F. H. 1967 The stability of finite amplitude cellular convection and its relation to an extremum principle. J. Fluid Mech. 30, 625–649.
  * Busse (1978) Busse, F. H. 1978 Non-linear properties of thermal convection. Rep. Prog. Phys. 41, 1929–1967.
  * Castillo-Castellanos et al. (2016) Castillo-Castellanos, A., Sergent, A. & Rossi, M. 2016 Reversal cycle in square Rayleigh–Bénard cells in turbulent regime. J. Fluid Mech. 808, 614–640.
  * Chandrasekhar (1961) Chandrasekhar, S. 1961 Hydrodynamic and hydromagnetic stability. Clarendon.
  * Chavanne et al. (1997) Chavanne, X., Chilla, F., Castaing, B., Hebral, B., Chabaud, B. & Chaussy, J. 1997 Observation of the ultimate regime in Rayleigh–Bénard convection. Phys. Rev. Lett. 79, 3648–3651.
  * Chavanne et al. (2001) Chavanne, X., Chillà, F., Chabaud, B., Castaing, B. & Hébral, B. 2001 Turbulent Rayleigh–Bénard convection in gaseous and liquid He. Phys. Fluids 13, 1300–1320.
  * Cross & Hohenberg (1993) Cross, M. C. & Hohenberg, P. C. 1993 Pattern formation outside of equilibrium. Rev. Mod. Phys. 65, 851–1112.
  * Emran & Schumacher (2012) Emran, M. S. & Schumacher, J. 2012 Conditional statistics of thermal dissipation rate in turbulent Rayleigh–Bénard convection. Eur. Phys. J. E 108, 35–42.
  * Farazmand (2016) Farazmand, M. 2016 An adjoint-based approach for finding invariant solutions of Navier–Stokes equations. J. Fluid Mech. 795, 278–312.
  * Goldhirsch et al. (1989) Goldhirsch, I., Pelz, R. B. & Orszag, S. A. 1989 Numerical simulation of thermal convection in a two-dimensional finite box. J. Fluid Mech. 199, 1––28.
  * Grossmann & Lohse (2000) Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: A unifying theory. J. Fluid Mech. 407, 27–56.
  * Grossmann & Lohse (2001) Grossmann, S. & Lohse, D. 2001 Thermal convection for large Prandtl numbers. Phys. Rev. Lett. 86, 3316–3319.
  * Grossmann & Lohse (2004) Grossmann, S. & Lohse, D. 2004 Fluctuations in turbulent Rayleigh–Bénard convection: The role of plumes. Phys. Fluids 16, 4462–4472.
  * Grossmann & Lohse (2011) Grossmann, S. & Lohse, D. 2011 Multiple scaling in the ultimate regime of thermal convection. Phys. Fluids 23, 045108\.
  * He et al. (2012) He, X., Funfschilling, D., Nobach, H., Bodenschatz, E. & Ahlers, G. 2012 Transition to the ultimate state of turbulent Rayleigh–Bénard convection. Phys. Rev. Lett. 108, 024502.
  * Hébert et al. (2010) Hébert, F., Hufschmid, R., Scheel, J. & Ahlers, G. 2010 Onset of Rayleigh–Bénard convection in cylindrical containers. Phys. Rev. E 81, 046318.
  * Hopf (1948) Hopf, E. 1948 A mathematical example displaying features of turbulence. Commun. Appl. Maths 1, 303–322.
  * Hu et al. (1993) Hu, Y., Ecke, R. & Ahlers, G. 1993 Convection near threshold for Prandtl numbers near 1. Phys. Rev. E 48, 4399–4413.
  * Johnston & Doering (2009) Johnston, H. & Doering, C. R. 2009 Comparison of Turbulent Thermal Convection between Conditions of Constant Temperature and Constant Flux. Phys. Rev. Lett. 102, 064501.
  * Julien & Watson (2009) Julien, K. & Watson, M. 2009 Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods. J. Comput. Phys. 228, 1480–1503.
  * Kooij et al. (2018) Kooij, G. L., Botchev, M. A., Frederix, E. M.A., Geurts, B. J., Horn, S., Lohse, D., van der Poel, E. P., Shishkina, O., Stevens, R. J. A. M. & Verzicco, R. 2018 Comparison of computational codes for direct numerical simulations of turbulent Rayleigh–Bénard convection. Comp. Fluids 166, 1–8.
  * Kooloth et al. (2021) Kooloth, P., Sondak, D. & Smith, L. M. 2021 Coherent solutions and transition to turbulence in two-dimensional Rayleigh–Bénard convection. Phys. Rev. Fluids 6, 013501\.
  * Kraichnan (1962) Kraichnan, R. 1962 Turbulent thermal convection at arbitrary Prandtl number. Phys. Fluids 5, 1374–1389.
  * Lohse & Xia (2010) Lohse, D. & Xia, K.-Q. 2010 Small-scale properties of turbulent Rayleigh–Bénard convection. Annu. Rev. Fluid Mech. 42, 335–364.
  * Mortensen (2018) Mortensen, M. 2018 Shenfun: High performance spectral Galerkin computing platform. J. Open Source Softw. 3, 1071\.
  * Niemela et al. (2000) Niemela, J. J., Skrbek, L., Sreenivasan, K. R. & Donnely, R. J. 2000 Turbulent convection at very high Rayleigh numbers. Nature 404, 837–841.
  * Oh (2019) Oh, S. 2019 An Efficient Spectral Method to Solve Multi-Dimensional Linear Partial Different Equations Using Chebyshev Polynomials. Mathematics 7, 90.
  * Petschel et al. (2011) Petschel, K., Wilczek, M., Breuer, M., Friedrich, R. & Hansen, U. 2011 Statistical analysis of global wind dynamics in vigorous Rayleigh–Bénard convection. Phys. Rev. E 84, 026309.
  * Podvin & Sergent (2015) Podvin, B. & Sergent, A. 2015 A large-scale investigation of wind reversal in a square Rayleigh–Bénard cell. J. Fluid Mech. 766, 172–201.
  * Podvin & Sergent (2017) Podvin, B. & Sergent, A. 2017 Precursor for wind reversal in a square Rayleigh–Bénard cell. Phys. Rev. E 95, 013112.
  * van der Poel et al. (2013) van der Poel, E. P., Stevens, R. J. A. M. & Lohse, D. 2013 Comparison between two- and three-dimensional Rayleigh–Bénard convection. J. Fluid Mech. 736, 177–194.
  * Puigjaner et al. (2004) Puigjaner, D., Herrero, J., Giralt, F. & Simó, C. 2004 Stability analysis of the flow in a cubical cavity heated from below. Phys. Fluids 16, 3639–3655.
  * 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. J. Fluid Mech. 598, 393–427.
  * Reiter (2021) Reiter, P. 2021 https://github.com/preiter93/rustpde.
  * Reiter et al. (2021a) Reiter, P., Shishkina, O., Lohse, D. & Krug, D. 2021a Crossover of the relative heat transport contributions of plume ejecting and impacting zones in turbulent rayleigh-bénard convection (a). EPL 134, 34002.
  * Reiter et al. (2021b) Reiter, P., Zhang, X., Stepanov, R. & Shishkina, O. 2021b Generation of zonal flows in convective systems by travelling thermal waves. J. Fluid Mech. 913, A13.
  * Roche (2020) Roche, P. E. 2020 The ultimate state of convection: a unifying picture of very high Rayleigh numbers experiments. New J. of Phys. 22, 073056.
  * Roche et al. (2001) Roche, P.-E., Castaing, B., Chabaud, B., Hébral, B. & Sommeria, J. 2001 Side wall effects in Rayleigh–Bénard experiments. Eur. Phys. J. B 24, 405–408.
  * Saad & Schultz (1986) Saad, Y. & Schultz, M. H. 1986 GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM J. Sci. Comput. 7, 856–869.
  * Schlüter et al. (1965) Schlüter, A., Lortz, D. & Busse, F. 1965 On the stability of steady finite amplitude convection. J. Fluid Mech. 23, 129––144.
  * Shen (1995) Shen, J. 1995 Efficient Spectral-Galerkin Method II. Direct Solvers of Second- and Fourth-Order Equations Using Chebyshev Polynomials. SIAM J. Sci. Comput. 16, 74–87.
  * Shishkina (2021) Shishkina, O. 2021 Rayleigh–Bénard convection: The container shape matters. Phys. Rev. Fluids 6, 090502.
  * Shishkina et al. (2014) Shishkina, O., Wagner, S. & Horn, S. 2014 Influence of the angle between the wind and the isothermal surfaces on the boundary layer structures in turbulent thermal convection. Phys. Rev. E 89, 033014.
  * Sondak et al. (2015) Sondak, D., Smith, L. M. & Waleffe, F. 2015 Optimal heat transport solutions for Rayleigh–Bénard convection. J. Fluid Mech. 784, 565–595.
  * Stevens et al. (2014) Stevens, R., Lohse, D. & Verzicco, R. 2014 Sidewall effects in Rayleigh–Bénard convection. J. Fluid Mech. 741, 1–27.
  * Sugiyama et al. (2010) Sugiyama, K., Ni, R., Stevens, R. J. A. M., Chan, T. S., Zhou, S.-Q., Xi, H.-D., Sun, C., Grossmann, S., Xia, K.-Q. & Lohse, D. 2010 Flow reversals in thermally driven turbulence. Phys. Rev. Lett. 105, 034503.
  * Tilgner et al. (1993) Tilgner, A., Belmonte, A. & Libchaber, A. 1993 Temperature and velocity profiles of turbulent convection in water. Phys. Rev. E 47, 2253–2257.
  * Urban et al. (2014) Urban, P., Hanzelka, P., Musilova, V., Kralik, T., Mantia, M. L., Srnka, A. & Skrbek, L. 2014 Heat transfer in cryogenic helium gas by turbulent Rayleigh–Bénard convection in a cylindrical cell of aspect ratio 1. New J. Phys. 16, 053042\.
  * Venturi et al. (2010) Venturi, D., Wan, X. & Karniadakis, G. 2010 Stochastic bifurcation analysis of Rayleigh–Bénard convection. J. Fluid Mech. 650, 391–413.
  * Verzicco (2002) Verzicco, R. 2002 Sidewall finite-conductivity effects in confined turbulent thermal convection. J. Fluid Mech. 473, 201–210.
  * Verzicco & Sreenivasan (2008) Verzicco, R. & Sreenivasan, K. R. 2008 A comparison of turbulent thermal convection between conditions of constant temperature and constant heat flux. J. Fluid Mech. 595, 203–219.
  * Wagner & Shishkina (2013) Wagner, S. & Shishkina, O. 2013 Aspect ratio dependency of Rayleigh–Bénard convection in box-shaped containers. Phys. Fluids 25, 085110.
  * Waleffe et al. (2015) Waleffe, F., Boonkasame, A. & Smith, L. M. 2015 Heat transport by coherent Rayleigh–Bénard convection. Phys. Fluids 27, 051702.
  * Wan et al. (2019) Wan, Z., Wei, P., Verzicco, R., Lohse, D., Ahlers, G. & Stevens, R. 2019 Effect of sidewall on heat transfer and flow structure in Rayleigh-–Bénard convection. J. Fluid Mech. 881, 218––243.
  * Wen et al. (2015) Wen, B., Chini, G. P., Kerswell, R. R. & Doering, C. R. 2015 Time-stepping approach for solving upper-bound problems: Application to two-dimensional Rayleigh–Bénard convection. Phys. Rev. E 92, 043012.
  * Wen et al. (2020a) Wen, B., Goluskin, D. & Doering, C. R. 2020a Steady Rayleigh–Bénard convection between no-slip boundaries, arXiv: 2008.08752.
  * Wen et al. (2020b) Wen, B., Goluskin, D., LeDuc, M., Chini, G. P. & Doering, C. R. 2020b Steady Rayleigh-–Bénard convection between stress-free boundaries. J. Fluid Mech. 905, R4.
  * Xi & Xia (2007) Xi, H.-D. & Xia, K.-Q. 2007 Cessations and reversals of the large-scale circulation in turbulent thermal convection. Phys. Rev. E 75, 066307.
  * Zhao et al. (2019) Zhao, J., Cai, W. & Jiang, Y. 2019 Study on corner vortex enlarging process of 2D square Rayleigh–Bénard cells filled with air in transient states. Int. J. Heat Mass Transfer 129, 599–609.
  * Zhou & Chen (2018) Zhou, W.-F. & Chen, J. 2018 Letter: Similarity model for corner roll in turbulent Rayleigh–Bénard convection. Phys. Fluids 30, 111705.
  * Zhu et al. (2018) Zhu, X., Mathai, V., Stevens, R., Verzicco, R. & Lohse, D. 2018 Transition to the ultimate regime in two-dimensional Rayleigh–Bénard convection. Phys. Rev. Lett. 120, 144502.
  * Zienicke et al. (1998) Zienicke, E., Seehafer, N. & Feudel, F. 1998 Bifurcations in two-dimensional Rayleigh–Bénard convection. Phys. Rev. E 57, 428–435.
  * Zwirner et al. (2020) Zwirner, L., Khalilov, R., Kolesnichenko, I., Mamykin, A., Mandrykin, S., Pavlinov, A., Shestakov, A., Teimurazov, A., Frick, P. & Shishkina, O. 2020 The influence of the cell inclination on the heat transport and large-scale circulation in liquid metal convection. J. Fluid Mech. 884, A18.