# Wall-attached convection under strong inclined magnetic fields

Shashwat Bhattacharya1<EMAIL_ADDRESS>Thomas Boeck1
Dmitry Krasnov1  Jörg Schumacher1,2  1Institute of Thermodynamics and Fluid
Mechanics, Technische Universität Ilmenau, Postfach 100565, D-98684 Ilmenau,
Germany 2Tandon School of Engineering, New York University, New York 11021,
USA

###### Abstract

We employ a linear stability analysis and direct numerical simulations to
study the characteristics of wall-modes in thermal convection in a rectangular
box under strong and inclined magnetic fields. The walls of the convection
cell are electrically insulated. The stability analysis assumes periodicity in
the spanwise direction perpendicular to the plane of the homogeneous magnetic
field. Our study shows that for a fixed vertical magnetic field, the
imposition of horizontal magnetic fields results in an increase of the
critical Rayleigh number along with a decrease in the wavelength of the wall
modes. The wall modes become tilted along the direction of the resulting
magnetic fields and therefore extend further into the bulk as the horizontal
magnetic field is increased. Once the modes localized on the opposite walls
interact, the critical Rayleigh number decreases again and eventually drops
below the value for onset with a purely vertical field. We find that for
sufficiently strong horizontal magnetic fields, the steady wall modes occupy
the entire bulk and therefore convection is no longer restricted to the
sidewalls. The above results are confirmed by direct numerical simulations of
the nonlinear evolution of magnetoconvection.

##  1 Introduction

Buoyancy-driven flows of electrically conducting fluids under the influence of
magnetic fields are a common occurrence in geophysical, astrophysical, as well
as in several technological applications. Such flows are called
magnetoconvection and their driving mechanism is the temperature dependence of
the fluid density, which results in spatial density variations leading to
buoyancy forces acting on the fluid. When such a fluid moves under the
influence of magnetic fields, electric currents are induced in the fluid due
to Faraday’s law, which, in turn, induce magnetic fields by the virtue of
Ampere’s law. These electric currents interact with the applied and induced
magnetic fields to generate a Lorentz force distribution that acts on the
fluid. Therefore, magnetoconvective flows are governed by equations of
conservation of mass, momentum, and thermal energy, along with Maxwell’s
equations for electromagnetism and Ohm’s law (Weiss & Proctor, 2014).
Magnetoconvection is encountered in the Sun, stars, and planetary dynamos. In
industries and technological applications, magnetoconvection is typically
encountered in liquid-metal batteries (Kelley & Sadoway, 2014; Shen & Zikanov,
2016; Kelley & Weier, 2018), cooling liquid-metal blankets in fusion reactors
(Mistrangelo et al., 2020, 2021), and magnetic stirring and braking of liquid
metal melts (Davidson, 1999, 2017; Lyubimov et al., 2010).

A simplified paradigm for magnetoconvection consists of a fluid layer that is
heated from below and cooled from above (Rayleigh-Bénard convection or RBC)
with imposed magnetic fields in different configurations. Typically, the
Boussinesq approximation is employed for modelling the above flows; this
approximation assumes that the flow is incompressible and the density
variations are negligible except in the buoyancy term in the momentum equation
(Chandrasekhar, 1981; Lohse & Xia, 2010; Chillà & Schumacher, 2012; Verma,
2018). Magnetoconvection is governed by the following nondimensional
parameters: i) Rayleigh number $\Ray$ – the ratio of buoyancy to dissipative
forces, ii) Prandtl number $\Pran$ – the ratio of kinematic viscosity to
thermal diffusivity, iii) Hartmann number $\Ha$ – the ratio of Lorentz to
viscous forces, and iv) the magnetic Prandtl number $\Pm$ – the ratio of
kinematic viscosity to magnetic diffusivity. The important nondimensional
output parameters of magnetoconvection are i) the Nusselt number $\Nu$ – the
ratio of the total heat transport to the diffusive heat transport, ii) the
Reynolds number $\Rey$ – the ratio of inertial to viscous forces, and iii) the
magnetic Reynolds number $\Rm$ – the ratio of induction to diffusion of the
magnetic field. In liquid-metal convection typically encountered in laboratory
experiments and most industrial applications, the magnetic Reynolds number is
sufficiently small such that the induced magnetic field is negligible compared
to the applied magnetic field and is thus neglected in the expressions of the
Lorentz force and Ohm’s law (Roberts, 1967; Davidson, 2017; Verma, 2019). Such
cases are referred to as quasi-static magnetoconvection where the induced
magnetic field adjusts instantaneously to the changes in velocity. In the
quasi-static approximation, there exists a one-way influence of the magnetic
field on the flow only.

Magnetoconvection has been studied theoretically in the past (for example,
Chandrasekhar, 1981; Houchens et al., 2002; Busse, 2008) as well as with the
help of experiments (for example, Nakagawa, 1957; Fauve et al., 1981; Cioni et
al., 2000; Aurnou & Olson, 2001; Burr & Müller, 2001; King & Aurnou, 2015;
Vogt et al., 2018, 2021; Zürner et al., 2020; Grannan et al., 2022) and
numerical simulations (for example, Liu et al., 2018; Yan et al., 2019;
Akhmedagaev et al., 2020a, b; Nicoski et al., 2022; Bhattacharya et al.,
2023). An application of horizontal magnetic fields causes the large-scale
rolls to become quasi two-dimensional and align in the direction of the field
(Fauve et al., 1981; Busse & Clever, 1983; Burr & Müller, 2002; Yanagisawa et
al., 2013; Tasaka et al., 2016; Vogt et al., 2018, 2021). These self-organized
flow structures reach an optimal state wherein the heat transport and
convective velocities increase significantly compared to convection without
magnetic fields (Vogt et al., 2021). In contrast, strong vertical magnetic
fields suppress convection (Chandrasekhar, 1981; Cioni et al., 2000; Zürner et
al., 2020; Akhmedagaev et al., 2020a, b). It must be noted that in a Rayleigh-
Bénard system, convection commences only above a certain critical Rayleigh
number, which is $\Ray_{c}\approx 1708$ for the case with infinite no-slip
horizontal walls (Chandrasekhar, 1981). For $\Ray<\Ray_{c}$, the heat transfer
occurs purely by diffusion. The critical Rayleigh number increases when a
vertical magnetic field is imposed and scales as $\Ray_{c}\sim\Ha^{2}$ in the
asymptotic limit of large Hartmann numbers.

The dynamics of convection under strong vertical magnetic fields become more
intricate in the presence of sidewalls. Houchens et al. (2002) and Busse
(2008) analytically showed that magnetoconvection near the sidewalls ceases at
Rayleigh numbers below the ones required to completely suppress convection in
the bulk. Several numerical and experimental studies on magnetoconvection with
sidewalls have also revealed the presence of residual wall-attached convection
at $\Ray<\Ray_{c}$ (Houchens et al., 2002; Liu et al., 2018; Akhmedagaev et
al., 2020a, b; Zürner et al., 2020; Teimurazov et al., 2023; McCormack et al.,
2023). These so-called wall modes were shown to exhibit a two-layered
structure and become more closely attached to the sidewalls with the increase
of Hartmann number (Liu et al., 2018).

There are a few studies only on convection with inclined magnetic fields which
motivates the present work (Hurlburt et al., 1996; Nicoski et al., 2022). The
results of Hurlburt et al. (1996) indicate that the mean flows tend to travel
in the direction of the tilt. Nicoski et al. (2022) observed qualitative
similarities between convection with inclined magnetic field and that with
vertical magnetic field in terms of the behavior of convection patterns, heat
transport, and flow speed. However, to the best of our knowledge, there are no
studies for the case with inclined magnetic fields where the Rayleigh number
is close to but less than the critical Rayleigh number. Therefore, in the
present work, we study thermal magnetoconvection in the wall-attached
convection regime and explore the effects of additional horizontal magnetic
fields on the wall modes. We use a combination of linear stability analysis
and direct numerical simulations to study the dependence of the horizontal
magnetic field strength, relative to the vertical magnetic field, on the wall-
mode structures and their impact on large-scale heat and momentum transport.

The outline of the paper is as follows. In § 2, we discuss the problem setup,
linear stability model, and the schemes for direct numerical simulations. The
linear stability analysis and the results of direct numerical simulations are
described in § 3. We conclude in § 4.

##  2 Numerical model

In this section, we discuss the mathematical model of our problem and the
numerics employed for the stability analysis and direct numerical simulations.
The study will be conducted under the quasi-static approximation, in which the
induced magnetic field is neglected as it is very small compared to the
applied magnetic field. This approximation is fairly accurate for
magnetoconvection in liquid metals (Davidson, 2017). Further, we employ
Boussinesq approximation, in which the variations in the density of the fluid
are ignored except in the buoyancy term in the momentum equation. Hence, the
flow is essentially treated as incompressible. The governing equations of
magnetoconvection under the above approximations are as follows:

$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=$ $\displaystyle 0$ (1)
$\displaystyle\frac{\partial\bm{u}}{\partial t}+\bm{u}\cdot\nabla\bm{u}$
$\displaystyle=$ $\displaystyle-\frac{\nabla p}{\rho}+\alpha
gT\hat{z}+\nu\nabla^{2}\bm{u}+\frac{1}{\rho}(\bm{j}\times{\bm{B}}),$ (2)
$\displaystyle\frac{\partial T}{\partial t}+\bm{u}\cdot\nabla T$
$\displaystyle=$ $\displaystyle\kappa\nabla^{2}T,$ (3) $\displaystyle\bm{j}$
$\displaystyle=$
$\displaystyle\sigma\left(-\nabla\phi+\bm{u}\times{\bm{B}}\right),$ (4)
$\displaystyle\nabla^{2}\phi$ $\displaystyle=$
$\displaystyle\nabla\cdot(\bm{u}\times{\bm{B}}),$ (5)

where $\bm{u}$, $\bm{j}$, $p$, $T$, and $\phi$ are the fields of velocity,
current density, pressure, temperature, and electrical potential respectively,
and ${\bm{B}}=(B_{x},B_{y},B_{z})$ is the applied magnetic field. In our work,
the magnetic fields are inclined along $y$-direction only, hence $B_{x}=0$.
Further, $\nu$ is the kinematic viscosity, $\kappa$ is the thermal
diffusivity, $\rho$ is the density, and $\sigma$ is the electrical
conductivity of the fluid. The last term in the momentum equation (2) is the
Lorentz force density. Equation (4) is Ohm’s law. The Poisson equation (5) for
the electric potential is a consequence of the charge conservation condition
$\nabla\cdot\bm{j}=0$.

In the following, we will discuss how the above equations have been employed
for our linear stability analysis and direct numerical simulations.

###  2.1 Linear stability model

We first discuss the derivation of the perturbation equations for our linear
stability analysis. The equations (1) to (5) are non-dimensionalized using the
cell height $H$ as the length scale, $\kappa/H$ as the velocity scale, the
temperature difference $\Delta$ between the two horizontal plates as the
temperature scale, and $B_{z}$, the vertical component of the applied magnetic
field. For this part of the analysis, we take the units that are typically
chosen for a linear stability analysis to end with a Prandtl number-
independent set of equations at the marginal stability threshold. The
characteristic units in the subsequent simulation part will differ. The non-
dimensionalized governing equations are as follows.

$\displaystyle\frac{1}{Pr}\left(\frac{\partial\bm{u}}{\partial
t}+(\bm{u}{\cdot\nabla})\bm{u}\right)=-{\nabla
p}+{\nabla}^{2}\bm{u}+\Ray\,\theta\bm{e}_{z}+\Ha_{z}^{2}\,\bm{j}\times\left(\bm{e}_{z}+R\bm{e}_{y}\right),$
(6) $\displaystyle\frac{\partial\theta}{\partial
t}+(\bm{u}{\cdot\nabla})\theta={\nabla}^{2}\theta+u_{z},$ (7)
$\displaystyle\nabla\cdot\bm{u}=0,$ (8)
$\displaystyle\bm{j}=-\nabla\phi+\bm{u}\times\left(\bm{e}_{z}+R\bm{e}_{y}\right),$
(9) $\displaystyle\nabla\cdot\bm{j}=0.$ (10)

In the above system of equations, $\Ray$ is the Rayleigh number, $\Pran$ is
the Prandtl number, $\Ha_{z}$ is the Hartmann number based on the vertical
component of the magnetic field, and $R$ is the ratio of the horizontal to the
vertical magnetic field strength. These quantities are the governing
parameters for our setup and are given by

$\Ray=\frac{\alpha g\Delta
H^{3}}{\nu\kappa},\quad\Pran=\frac{\nu}{\kappa},\quad\Ha_{z}=B_{z}H\sqrt{\frac{\sigma}{\rho\nu}},\quad
R=\frac{B_{y}}{B_{z}}.$ (11)

The quantity $\theta$ is the difference between the temperature and the linear
conduction profile, i.e.,

$T({\bm{x}},t)=\theta({\bm{x}},t)-z\,.$ (12)

Apart from $\Ray$, $\Pran$, $\Ha_{z}$, and $R$, the dynamics are also governed
by the aspect ratio $\Gamma$, which is the ratio of the length to the height
of the convection cell.

We examine a total of 12 cases of magnetoconvection in a bounded,
horizontally-extended domain of dimension $\Gamma\times 1$ in the $y$-$z$
plane. Two aspect ratios are considered: $\Gamma=2$ and $\Gamma=4$. For
$\Gamma=2$, we consider the cases of $\Ha_{z}=50$ and $\Ha_{z}=100$, whereas
for $\Gamma=4$, we consider only the case of $\Ha_{z}=50$. For each $\Ha_{z}$
analysed in our study, we vary $R$ from 0 (corresponding to a purely vertical
magnetic field) to 3 in steps of 1. The convection cell is periodic in
$x$-direction and consists of no-slip horizontal walls at $z=\pm 1/2$ and two
no-slip sidewalls at $y=\pm\Gamma/2$. All the walls are electrically
insulated. Each horizontal wall is at a constant temperature, with the bottom
wall at $T=0.5$ and the top wall at $T=-0.5$. The sidewalls are thermally
insulated with $\partial T/\partial\eta=0$, where $\eta$ is the direction
normal to the wall. A sketch of our setup is shown in figure 1.

Figure 1: A sketch of the Rayleigh-Bénard convection setup with inclined
magnetic field employed for our linear stability analysis.

In the present work we assume that the instability is of stationary type. The
nonlinear terms as well as the time derivatives in equations (6) and (7) are
therefore neglected. The momentum and continuity equations reduce to the
Stokes problem with additional buoyancy and Lorentz force. To avoid
complications stemming from the coupling between pressure and velocity we
choose a representation for the velocity that satisfies the continuity
equation automatically and eliminate the pressure term. The velocity field is
written as the curl of a vector streamfunction $\bm{\psi}$

$\bm{u}=\nabla\times\bm{\psi},$ (13)

and the gauge condition $\nabla\cdot\bm{\psi}=0$ is imposed to determine
$\bm{\psi}$ uniquely as in Priede et al. (2010). The dependence on $x$ is
represented by the normal mode ansatz with wavenumber $\beta$ for all fields,
e.g. $\theta(y,z)\exp(i\beta x)$ for the temperature perturbation. The gauge
condition allows one to express the $x$-component of $\bm{\psi}$ by

$i\beta\psi_{x}(y,z)=-\partial_{y}\psi_{y}(y,z)-\partial_{z}\psi_{z}(y,z).$
(14)

The velocity components then read

$\displaystyle u_{x}$
$\displaystyle=\partial_{y}\psi_{z}-\partial_{z}\psi_{y}\,,$ (15)
$\displaystyle u_{y}$
$\displaystyle=-i\beta\psi_{z}+\frac{i}{\beta}\left(\partial_{y}\partial_{z}\psi_{y}+\partial_{z}^{2}\psi_{z}\right)\,,$
(16) $\displaystyle u_{z}$
$\displaystyle=i\beta\psi_{y}-\frac{i}{\beta}\left(\partial_{y}^{2}\psi_{y}+\partial_{y}\partial_{z}\psi_{z}\right)\,.$
(17)

Equations for $\psi_{y}$ and $\psi_{z}$ are obtained by taking the curl of the
definition (13) and the momentum equation (6). They are

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\nabla^{2}\psi_{y}+\omega_{y},$ (18) $\displaystyle 0$
$\displaystyle=$ $\displaystyle\nabla^{2}\psi_{z}+\omega_{z},$ (19)
$\displaystyle 0$ $\displaystyle=$
$\displaystyle\nabla^{2}\omega_{y}-i\beta\,\Ray\,\theta+\Ha_{z}^{2}\left(-\partial_{y}\partial_{z}\phi-\partial_{z}u_{x}+R\left(-\partial_{y}^{2}\phi+\omega_{z}-i\beta
u_{y}\right)\right),$ (20) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\nabla^{2}\omega_{z}+\Ha_{z}^{2}\left(-\partial_{z}^{2}\phi+R\left(-\partial_{y}\partial_{z}\phi+\omega_{y}+i\beta
u_{z}+R\partial_{y}u_{x}\right)\right).$ (21)

The quantities $\omega_{y}$ and $\omega_{z}$ are the $y$\- and $z$-components
of the vorticity field $\nabla\times\bm{u}$. Equations for the remaining
quantities are

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\nabla^{2}\theta+u_{z},$ (22)
$\displaystyle 0$ $\displaystyle=$
$\displaystyle\nabla^{2}\phi-\omega_{z}-R\omega_{y}.$ (23)

The last equation (23) is obtained by substitution of Ohm’s law (9) into (10).
Combined with boundary conditions on the top wall, bottom wall and side walls
specified below, equations (18)–(23) represent a linear eigenvalue problem for
the Rayleigh number $\Ray$ that must be solved numerically. A suitable
discretization of this problem is obtained by a spectral collocation method
with Chebyshev polynomials $T_{n}(z)=\cos\\{n\arccos(z)\\}$. The scalar fields
such as $\theta$ are expanded as

$\theta(y,z)=\sum_{i}\sum_{k}\theta_{ik}\,T_{i}(2y/\Gamma)\,T_{k}(2z),$ (24)

where $-\Gamma/2\leq y\leq\Gamma/2$ and $-1/2\leq z\leq 1/2$. The Poisson
equations (20–23) and boundary conditions are imposed pointwise at the Gauss-
Lobatto collocation points

$y_{j}=\Gamma\cos(j\pi/N_{y})/2\quad(0\leq j\leq N_{y}),\qquad
z_{k}=\cos(k\pi/N_{z})/2\quad(0\leq k\leq N_{z}),$ (25)

where $N_{y}+1$ and $N_{z}+1$ are the number of expansion terms with respect
to $y$ and $z$.

The boundary conditions for the vector stream function and vorticity
components are determined with the help of equations (15)-(17). Zero normal
velocity on the horizontal walls requires $\psi_{y}=0$ and
$\partial_{z}\psi_{z}=0$. On the $y$=$\pm\Gamma/2$ sidewalls, the
corresponding conditions are $\psi_{z}=0$ and $\partial_{y}\psi_{y}=0$. The
tangential velocity vanishes on the sidewalls if $\omega_{y}=0$ and
$\partial_{y}\psi_{z}=\partial_{z}\psi_{y}$. On the top and bottom walls these
conditions are $\omega_{z}=0$ and $\partial_{y}\psi_{z}=\partial_{z}\psi_{y}$.
The remaining boundary conditions for (22) are $\theta=0$ on the top and
bottom walls and $\partial_{y}\theta=0$ on the $y$=$\pm\Gamma/2$ sidewalls.
The boundary condition for the electric potential supplementing equation (23)
is the homogeneous Neumann condition.

Since the boundary conditions (zero normal velocity) for (18) and (19) only
involve $\psi_{y}$ and $\psi_{z}$, respectively, one can represent the
expansion coefficients of $\psi_{y}$ and $\psi_{z}$ by linear invertible maps
through those of $\omega_{y}$ and $\omega_{z}$ (assuming the latter are
augmented by the zero boundary values to be imposed on either $\psi_{y}$,
$\psi_{z}$ or its normal derivatives). The expansions for $\psi_{y}$ and
$\psi_{z}$ therefore contain $N_{y}+3$ and $N_{z}+3$ terms, respectively. The
values of $\psi_{y}$, $\psi_{z}$ or its derivatives in equations (20)–(23)
(and associated boundary conditions) at the collocation points are represented
through expansion coefficients of $\omega_{y}$ and $\omega_{z}$ via these
linear invertible maps. The same can be done for the electric potential, which
is the sum of a contribution from $\omega_{z}$ and $\omega_{y}$. As a result
of the collocation approximation one obtains a vector $\bm{Y}$ of unknowns
containing the expansion coefficients of $\omega_{y}$, $\omega_{z}$, $\theta$
with a size of $3(N_{y}+1)(N_{z}+1)$ and a generalized linear eigenvalue
problem

$\bm{\mathcal{A}}\bm{Y}=\Ray\,\bm{\mathcal{B}}\bm{Y}.$ (26)

The method was implemented in Matlab (The MathWorks Inc., 2022) using the
default double precision. Notice that $\omega_{y}$, $\omega_{z}$ and $\phi$
are real variables. According to equations (15-17) and (20-22), $u_{y}$,
$u_{z}$ and $\theta$ would be purely imaginary quantities. They are considered
as real variables in the code and below. Problem (26) was solved with Matlab’s
eig routine to find all eigenvalues and eigenvectors. The routine also works
with a matrix $\bm{\mathcal{B}}$ whose rank is smaller than the rank of
$\bm{\mathcal{A}}$ (as it is the case for (26)). It associates the spurious
solutions that stem from equations not containing the eigenvalue $Ra$ with
infinite eigenvalues.

For the cases studied in this paper, the numerical resolution was typically
$N_{y}=70$, $N_{z}=60$ for $Ha_{z}=50$ and $N_{y}=90$, $N_{z}=70$ for
$Ha_{z}=100$ with aspect ratio $\Gamma=2$. A decrease of $N_{y}$ by 10
typically only resulted in a relative change of the first eigenvalue below
$10^{-5}$. The generation of the matrices and the solution of problem (26) for
a given wavenumber took about 20 hours for $N_{y}=90$, $N_{z}=70$ and about 6
hours for $N_{y}=70$, $N_{z}=60$ on an Intel Xeon E5 CPU.

###  2.2 Direct numerical simulations

We conduct direct numerical simulations (DNS) of our magnetoconvection setup
using a second-order finite difference code developed by Krasnov et al.
(2011). The governing equations are made dimensionless by using the cell
height $H$, the imposed temperature difference $\Delta$, and the free-fall
velocity $U=\sqrt{\alpha g\Delta H}$ (where $g$ and $\alpha$ are respectively
the gravitational acceleration and the volumetric coefficient of thermal
expansion of the fluid). The following non-dimensional equations are employed
for our DNS:

$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=$ $\displaystyle 0$ (27)
$\displaystyle\frac{\partial\bm{u}}{\partial t}+\bm{u}\cdot\nabla\bm{u}$
$\displaystyle=$ $\displaystyle-\nabla
p+T\hat{z}+\sqrt{\frac{\Pran}{Ra}}\nabla^{2}\bm{u}+\Ha_{z}^{2}\sqrt{\frac{\Pran}{Ra}}(\bm{j}\times{\bm{B}}),$
(28) $\displaystyle\frac{\partial T}{\partial t}+\bm{u}\cdot\nabla T$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{Ra\Pran}}\nabla^{2}T,$ (29)
$\displaystyle\bm{j}$ $\displaystyle=$
$\displaystyle-\nabla\phi+(\bm{u}\times{\bm{B}}),$ (30)
$\displaystyle\nabla^{2}\phi$ $\displaystyle=$
$\displaystyle\nabla\cdot(\bm{u}\times{\bm{B}}),$ (31)

We numerically solve equations (27) to (31). The Prandtl number $\Pran$ is
chosen to be 0.025, which is the same as that of mercury. We choose two
Hartmann numbers based on vertical magnetic field: $\Ha_{z}=100$ and
$\Ha_{z}=200$. For these Hartmann numbers, the corresponding critical Rayleigh
numbers ($\Ray_{c}$) for the case without horizontal walls obtained using the
linear stability analysis of Chandrasekhar (1981) is:

$\Ray_{c,\infty}=\begin{cases}1.245\times 10^{5},&\Ha_{z}=100,\\\ 4.48\times
10^{5},&\Ha_{z}=200.\end{cases}$ (32)

Since we are interested in the wall-attached convection regime, we choose the
Rayleigh number to be slightly below the critical Rayleigh numbers given by
(32). Thus, the chosen Rayleigh numbers are $\Ray=10^{5}$ for the case of
$\Ha_{z}=100$ and $\Ray=4\times 10^{5}$ for the case of $\Ha_{z}=400$. For
each $\Ha_{z}$, we vary $R$ from 0 to 3.

We choose the domain-size to be $\Gamma\times\Gamma\times H=4\times 4\times 1$
and employ a grid-resolution ranging from $1200\times 1200\times 300$ points
to $1800\times 1800\times 400$. The horizontal walls are at $z=\pm 1/2$ and
the sidewalls are at $x=\pm\Gamma/2$ and $y=\pm\Gamma/2$. The mesh is non-
uniform with stronger clustering of the grid points near the top and bottom
boundaries. The elliptic equations for pressure, electric potential, and the
temperature are solved based on applying cosine transforms in $x$\- and
$y$-directions and using a tridiagonal solver in the $z$-direction. The
diffusive term in the temperature transport equation is treated implicitly.
The time discretization of the momentum equation uses the fully explicit
Adams-Bashforth/Backward-Differentiation method of second order (Peyret,
2002). A constant time step size ranging from $5\times 10^{-5}$ to $1\times
10^{-4}$ free fall time unit was chosen for our simulations, which satisfied
the Courant–Friedrichs–Lewy (CFL) condition for our runs.

All the walls are rigid and electrically insulated such that the electric
current density $\bm{j}$ forms closed field lines inside the cell. The top and
bottom walls are held fixed at $T=-0.5$ and $T=0.5$ respectively, and the
sidewalls are adiabatic with $\partial T/\partial\eta=0$ (where $\eta$ is the
component normal to sidewall). All the simulations are initialized with the
linear conduction profile for temperature (which is a function of the
$z$-coordinate only) and a random noise of amplitude $A=0.001$ along the
$z$-direction for velocity. We run the simulations initially on a coarse grid
of $120\times 120\times 30$ points for 100 free-fall time units in which they
converge to a statistically steady state. Following this, we successively
refine the mesh to the required resolutions and allow the simulations to
converge after each refinement. Once the simulations reach the statistically
steady state at the highest resolution, they are run for another 20 to 21
free-fall time units and a snapshot of the flow field is saved after every
free-fall time unit.

Since all the walls are no-slip, thin velocity boundary layers are formed
adjacent to the walls. For our simulations to be well-resolved, an adequate
number of gridpoints need to be present in these boundary layers. It must be
noted that the boundary layer profiles are strongly influenced by the magnetic
fields. For a purely vertical magnetic field, these boundary layers are
categorized into Hartmann layers adjacent to the top and bottom walls and
Shercliff layers adjacent to the sidewalls. The thickness of Hartmann layers
is given by $\delta_{H}=1/\Ha_{z}$ and that for the Shercliff layers is given
by $\delta_{S}=1/\sqrt{\Ha_{z}}$. However, in our case, the magnetic field is
inclined with respect to the vertical direction; therefore, both the
horizontal walls and $y$=$\pm\Gamma/2$ sidewalls will have a mix of Hartmann
and Shercliff layers. On the other hand, the $x$=$\pm\Gamma/2$ sidewalls will
have purely Shercliff layers. Thus, we use $\delta$ only for representing the
boundary-layer thickness. For a conservative analysis, we estimate the
thicknesses of these boundary layers as follows.

$\delta=\begin{cases}\min(\frac{1}{\Ha_{z}},\frac{1}{\sqrt{\R\Ha_{z}}}),&\mathrm{Horizontal~{}walls},\\\
\min(\frac{1}{\sqrt{\Ha_{z}}},\frac{1}{\R\Ha_{z}}),&y\mbox{=}\pm\Gamma/2~{}\mathrm{sidewalls},\\\
\dfrac{1}{(\Ha_{z}^{2}(1+R^{2}))^{1/4}},&x\mbox{=}\pm\Gamma/2~{}\mathrm{sidewalls}.\end{cases}$
(33)

Table 1 lists the important parameters of our simulation runs. In this table,
we also report the number of points in the different velocity boundary layers.
We ensure that a minimum of 10 points is present in the boundary layers so as
to adequately resolve our simulation runs.

In the next section, we will discuss the results obtained from our stability
analysis and numerical simulations.

Runs | $\Ha_{z}$ | $\Ray$ | $R$ | Grid size | $n_{z}$ | $n_{y}$ | $n_{x}$
---|---|---|---|---|---|---|---
1 | 100 | $1\times 10^{5}$ | 0 | $1200\times 1200\times 300$ | 22 | 89 | 89
2 | 100 | $1\times 10^{5}$ | 0.3 | $1200\times 1200\times 300$ | 22 | 39 | 88
3 | 100 | $1\times 10^{5}$ | 1 | $1200\times 1200\times 300$ | 22 | 15 | 79
4 | 100 | $1\times 10^{5}$ | 2 | $1200\times 1200\times 300$ | 22 | 15 | 120
5 | 100 | $1\times 10^{5}$ | 3 | $1200\times 1200\times 300$ | 22 | 10 | 106
6 | 200 | $4\times 10^{5}$ | 0 | $1200\times 1200\times 300$ | 12 | 21 | 21
7 | 200 | $4\times 10^{5}$ | 0.3 | $1200\times 1200\times 300$ | 12 | 21 | 21
8 | 200 | $4\times 10^{5}$ | 1 | $1200\times 1200\times 300$ | 12 | 13 | 18
9 | 200 | $4\times 10^{5}$ | 2 | $1800\times 1800\times 400$ | 16 | 11 | 22
10 | 200 | $4\times 10^{5}$ | 3 | $1800\times 1800\times 400$ | 16 | 12 | 18
Table 1: Parameters of the simulations: the Hartmann number ($\Ha_{z}$) based
on the vertical magnetic field, the Rayleigh number ($\Ray$), the ratio $R$ of
the magnetic field strength along the horizontal to the vertical directions,
the grid-size, the number of points in the velocity boundary layers along the
horizontal walls ($n_{z}$), the $y$=$\pm\Gamma/2$ sidewalls ($n_{y}$), and the
$x$=$\pm\Gamma/2$ sidewalls ($n_{x}$). The Prandtl number is constant in all
cases at a value of $\Pran=0.025$.

##  3 Results

In this section, we make a detailed analysis of the structure of the wall-
modes and their impact on the heat and momentum transport. We will first
present the linear stability analysis which is followed by a discussion of the
results obtained from our direct numerical simulations.

###  3.1 Linear stability analysis

In this subsection, we discuss the results from our linear stability model. We
assume that the set of finite positive eigenvalues of problem (26) is sorted
in ascending order, i.e. the smallest one will be referred to as the first
eigenvalue etc.

In figures 2(a-c), we plot the first eigenvalue, denoted as Rayleigh number
$\Ray$, at which the instability sets in due to disturbances at wavenumber
$\beta$ for different values of $R$. The minimum value of $\Ray$ is the
critical Rayleigh number $\Ray_{c}$, and the wavenumber corresponding to
$\Ray_{c}$ is the most unstable wavenumber. We also plot $\Ray$ for the
infinite plane layer for the corresponding $\Ha_{z}$, where the flow is
assumed to be uniform in the $y$-direction. The figures show that $\Ray_{c}$
is smaller than that corresponding to the plane infinite layer
($\Ray_{c,\infty}$) for all $R$, $\Ha_{z}$, and aspect ratios considered in
this study. This clearly implies that sidewalls destabilize the
magnetoconvection system. Figure 3(a) exhibits the dependence of $\Ray_{c}$ on
$R$ for different aspect ratios and strengths of the vertical magnetic field.
It is evident from the above figure that for $R<1$, the critical Rayleigh
number increases for all aspect ratios. For the larger aspect ratio box
($\Gamma=4$), $\Ray_{c}$ continues to increase with $R$ beyond $R=1$ and
saturates at $R=2$. On the other hand, for the smaller aspect ratio box, the
critical Rayleigh number starts to decrease with $R$ for $R>1$, a trend that
does not seem to depend on $\Ha_{z}$.

In figure 3(b), we exhibit the plots of the most unstable wavelength $\lambda$
versus $R$. The most unstable wavelength is calculated as

$\lambda=\frac{2\pi}{\beta_{c}},$

where $\beta_{c}$ is the most unstable wavenumber. The figure shows that
$\lambda$, like $\Ray_{c}$, exhibits a non-monotonic variation with $R$ and
lies between 1.1 and 1.6 for the entire range of parameters considered in our
study. For small values of $R$, $\lambda$ decreases with $R$ and for the
larger aspect-ratio box, it saturates at $R=2$. For the smaller aspect-ratio
box, $\lambda$ increases with increasing $R$ beyond $R=2$.

Figure 2: Neutral stability curves for magnetoconvection with different values
of $R$ for (a) $\Gamma=2$, $\Ha_{z}=50$; (b) $\Gamma=4$, $\Ha_{z}=50$; and (c)
$\Gamma=2$, $\Ha_{z}=100$.

Figure 3: For $\Ha_{z}=50$, $\Gamma=2$ (red squares); $\Ha_{z}=50$, $\Gamma=4$
(blue circles); and $\Ha_{z}=100$, $\Gamma=2$ (black triangles): (a) critical
Rayleigh number $\Ray_{c}$, normalized with the same for infinite horizontal
layer, and (b) the wavelength of the most unstable mode at $\Ray_{c}$.

We now examine the behavior of the first three eigenvalues ($\Ray$) for our
cases. In figures 4(a)–(l), we plot the variations of these eigenvalues with
$\beta$. The neutral stability curve for the infinite plane convection layer
for the corresponding $\Ha_{z}$ is also shown in dashed curves. For $R=0$, it
can be seen in figures 4(a,e,i) that the minima of the third eigenvalue
(represented by green triangles) overlaps with that of the plane convection
layer. This indicates that the third eigenvalue corresponds to the onset of
bulk convection for different wavenumbers. The figures also indicate that
minima of the third eigenvalue, which corresponds to the critical Rayleigh
number for convection in the bulk, increases with $R$. This indicates that the
bulk convection is further suppressed as the horizontal magnetic field
increases.

The first and second eigenvalues (denoted by black squares and red circles
respectively) correspond to wall-attached convection, and their minima are
less than that for the third eigenvalue (corresponding to bulk convection).
These eigenvalues nearly overlap for small horizontal magnetic fields but
begin to diverge at large horizontal magnetic fields. These variations are
shown more explicitly in figures 5(a,b) in which we exhibit the variations of
the first two eigenvalues with $R$ for $\beta=4$ (near the most unstable
wavenumber). It is evident from these figures that the eigenvalues diverge
visibly for $R>R_{c}\approx 1.5$ for $\Gamma=2$ and for $R>R_{c}\approx 2.5$
for $\Gamma=4$. It is interesting to note that the value of $R$ for which the
eigenvalues begin to diverge does not seem to depend on $\Ha_{z}$.

Figure 4:  Dependence of the first three eigenvalues on the wavenumber $\beta$
for $\Gamma=2$ and $\Ha_{z}=50$ with (a) $R=0$, (b) $R=1$, (c) $R=2$, and (d)
$R=3$; $\Gamma=4$ and $\Ha_{z}=50$ with (e) $R=0$, (f) $R=1$, (g) $R=2$, and
(h) $R=3$; and $\Gamma=2$ and $\Ha_{z}=100$ with (i) $R=0$, (j) $R=1$, (k)
$R=2$, and (l) $R=3$. Also shown is the neutral stability curve (dashed blue
lines) for the corresponding infinite horizontal convection layer
($\Gamma=\infty$) with a purely vertical magnetic field ($R=0$). For the
infinite plane layer, the flow is assumed to be uniform in the $y$-direction.

Figure 5: Influence of (a) aspect ratio and (b) vertical Hartmann number on
the stability curves for 1st and 2nd eigenvalues as function of the ratio of
the horizontal to vertical magnetic fields at fixed wavenumber $\beta=4$.

Figure 6: Eigenvector field from the linear stability analysis. Contours of
the temperature deviation $\theta$ are shown. For $\beta=4$, $\Gamma=2$, and
$\Ha_{z}=50$: isocontours of $\theta$ corresponding to the first solution for
(a) $R=0$, (b) $R=1$, (c) $R=2$, and (d) $R=3$. For $\beta=4$, $\Gamma=2$, and
$\Ha_{z}=50$: isocontours of $\theta$ corresponding to the second solution for
(e) $R=0$, (f) $R=1$, (g) $R=2$, and (h) $R=3$. For $\beta=4$, $\Gamma=2$, and
$\Ha_{z}=100$: isocontours of $\theta$ corresponding to the first solution for
(i) $R=0$, (j) $R=1$, (k) $R=2$, and (l) $R=3$. For $\beta=4$, $\Gamma=2$, and
$\Ha_{z}=100$: isocontours of $\theta$ corresponding to the second solution
for (m) $R=0$, (n) $R=1$, (o) $R=2$, and (p) $R=3$. Also shown are the vector
plots of the corresponding magnetic fields.

Figure 7: Ratio of the difference between the critical Rayleigh numbers for
antisymmetric solution ($\Ray_{a}$) and the symmetric solution ($\Ray_{s}$) to
the critical Rayleigh number of the symmetric solution. The inset exhibits a
magnified version of the above plot with $R$ ranging from 0 to 0.8. Parameters
are $Ha_{z}=50$, $\Gamma=2$ and $\beta=4$.

Figure 8: Eigenvector field from the linear stability analysis. Contours of
the temperature deviation $\theta$ are shown. For $\Ha_{z}=50$, $\beta=4$, and
$\Gamma=4$: isocontours of $\theta$ corresponding to the first eigensolutions
for (a) $R=0$, (c) $R=1$, (e) $R=2$, and (g) $R=3$, and corresponding to
second eigensolutions for (b) $R=0$, (d) $R=1$, (f) $R=2$, and (h) $R=3$. Also
shown are the vector plots of the corresponding magnetic fields.

We will now examine the spatial structure of the wall-modes, i.e. the
eigenfunctions that correspond to the first and second eigenvalues. The
eigenfunctions are normalized such that the maximum of the vertical velocity
becomes equal to unity.

In figures 6(a)–(p), we exhibit the contour plots of the temperature
perturbation $\theta$ on the vertical $y$-$z$ plane corresponding to the first
and second eigenvalues for $\beta=4$, $\Gamma=2$, and different $\Ha_{z}$ and
$R$. The first and third rows of the figure correspond to the first
eigensolutions for $\Ha_{z}=50$ and $\Ha_{z}=100$, respectively, whereas the
second and fourth rows correspond to the second eigensolutions. The figures
show that as the horizontal magnetic field strength is increased, the wall
modes get elongated and tilt along the direction of the resultant magnetic
field. In fact, for $R\geq 2$, the wall modes are no longer confined adjacent
to the sidewalls; instead, they occupy almost the entire bulk and interact
with each other as explained later in this section.

Figures 6(a)–(p) also show that two there are two types of spatial structures
corresponding to the first two eigensolutions. The first type consists of a
hot plume ($\theta>0$) adjacent to the $y=-\Gamma/2$ sidewall and a cold plume
($\theta<0$) adjacent to $y=\Gamma/2$ sidewall; the corresponding
eigensolution will be referred to as antisymmetric solution, see panel (a).
The second type consists of cold (or hot) plumes adjacent to both the
sidewalls; the corresponding eigensolution will be referred to as symmetric
solution, see panel (e). For $R<1.5$, there is only a marginal difference
between the symmetric ($\Ray_{s}$) and antisymmetric ($\Ray_{a}$) eigenvalues.
As exhibited in figure 7, this difference is less than 0.5% for $R<0.8$ and is
approximately $10\%$ at $R=1.5$ for the case of $\Ha_{z}=50$ and $\Gamma=2$.
Thus, for $R<1.5$, there is nearly an equal preference for symmetric and
antisymmetric structures to develop at the onset of convection. However, these
eigenvalues deviate significantly once R exceeds the threshold $R_{c}\approx
1.5$, above which the eigenvalue corresponding to the symmetric solution
becomes significantly smaller. This implies that there is a stronger
preference for the symmetric structures to develop at the onset of convection
for $R>R_{c}$. In this regime of $R$, the symmetric eigensolution comprises of
a large plume developed by the merging of two cold (or hot) plumes adjacent to
the opposite walls as visible in figures 6(c), (d), (k), and (l). On the other
hand, as seen in figures 6(g), (h), (o), and (p), the antisymmetric
eigensolutions for $R>R_{c}$ comprise of a cold plume extending on the top of
a hot plume (or vice-versa) extended from the opposite wall, resulting in two
convection rolls on top of each other. Our observations imply that a clear
separation of the symmetric and antisymmetric solutions occurs when the wall
modes adjacent to the opposite walls for symmetric solutions start interacting
with each other. The merging of the plumes and the resultant formation of a
merged roll results in an increased heat and momentum transport and hence in
the decrease of the critical Rayleigh number.

Figures 8(a–h) exhibit the contour plots of $\theta$ on vertical $y$-$z$
midplane corresponding to the first and second eigensolutions for $\Gamma=4$,
$\beta=4$, and $\Ha_{z}=50$. These figures again show that the wall modes tend
to elongate along the direction of the resultant magnetic field. Similar to
the $\Gamma=2$ cases, the two eigenvalues correspond to symmetric and
antisymmetric solutions respectively. The wall modes are symmetric for the
first eigensolution and antisymmetric for the second. It can be recalled from
figure 5(a) that the first and second eigenvalues for $\Gamma=4$ box begin to
diverge at a higher value of $R$ compared to the $\Gamma=2$ box. A clear
separation occurs near $R_{c}\approx 2.5$; this is because owing to the larger
aspect ratio of the box, the plumes adjacent to opposite walls are able to
interact and merge only at higher tilts, and hence at larger $R$. The merged
plume is exhibited in figure 8(g) which displays the contours of $\theta$
corresponding to the symmetric solution for $R=3$.

Our analysis suggests that the nonmonotonic behaviour of $\Ray_{c}$ and
wavelength $\lambda$ with respect to $R$ is due to the increasing interaction
of the plumes on opposite sidewalls. As long as the opposite plumes do not
merge, an increase in the horizontal magnetic field component further
stabilizes the magnetoconvection system with $\Ray_{c}$ increasing and the
wavelength $\lambda$ decreasing with $R$. The system gets destabilized due to
the merging of the opposite plumes and the variations of $\Ray_{c}$ and
$\lambda$ with $R$ get reversed. In the next subsection, we will discuss the
results of direct numerical simulations of the nonlinear evolution of
magnetoconvection.

###  3.2 Results of direct numerical simulations

In this subsection, we analyze our steady-state DNS results and examine the
structures of the wall modes, their role in the heat and momentum transport,
and the effects of initial conditions on the formation of the wall modes.

####  3.2.1 Structure of the wall modes

We use our numerical data to study the spatial convection structures for
$\Ha_{z}=100$ and $\Ha_{z}=200$ with $R$ ranging from 0 to 3.

In figures 9(a)–(j), we exhibit the isosurfaces of $u_{z}=\pm 0.01$ for all
our runs. In figures 10(a)–(h), we display the contours of $u_{z}\geq 0.01$
(red) and $u_{z}\leq 0.01$ (blue); these contours give a visualization of the
upwelling and downwelling plumes respectively. These figures show the presence
of wall-attached convection with suppressed fluid flow in the bulk. Consistent
with the results of the stability analysis in the previous subsection, these
wall modes become elongated and align themselves along the direction of the
resultant magnetic field as the horizontal magnetic field is increased.
Figures 9 and 10 also show that with the exception of the case with
$\Ha_{z}=100$ and $R=3$, the wall modes are antisymmetric, that is, the plumes
are upwelling (or downwelling) adjacent to $y=-\Gamma/2$ wall and downwelling
(or upwelling) adjacent to $y=\Gamma/2$ wall. The wall modes occupy the entire
bulk for $R=3$, consistent with the stability analysis of the $\Gamma=4$ box.
For the above field ratio, the wall modes for the $\Ha_{z}=100$ case are
symmetric and the plumes adjacent to the opposite walls merge to form a large
plume as displayed in figures 9(e) and 10(d).

Figure 9: Results of direct numerical simulations for the vertical velocity
component. Isosurfaces of $u_{z}=0.01$ (red) and $u_{z}=-0.01$ (blue) for
$Ha_{z}=100$ with (a) $R=0$, (b) $R=0.3$, (c) $R=1$, (d) $R=2$, and (e) $R=3$,
and for $\Ha_{z}=200$ with (f) $R=0$, (g) $R=0.3$, (h) $R=1$, (i) $R=2$, and
(j) $R=3$.

Figure 10: Results of direct numerical simulations for the wall modes.
Contours of $u_{z}\geq 0.01$ (red) and $u_{z}\leq-0.01$ (blue) along with the
vector plots of $\bm{B}$ on $x=0$ midplane for $Ha_{z}=100$ with (a) $R=0$,
(b) $R=1$, (c) $R=2$, (d) $R=3$, and for $\Ha_{z}=200$ with (e) $R=0$, (f)
$R=1$, (g) $R=2$, (h) $R=3$. The wall modes align themselves along the
direction of $\bm{B}$.

It can also be observed in figures 9(a–j) that there are slight irregularities
in the wall-mode structures. These irregularities seem to be associated with
small temporal fluctuations ($\sim 1\%$) in integral quantities that remain
after our simulations reach an apparently stationary state. These fluctuations
indicate that the wall modes are still evolving, albeit slowly. However, it is
important to note that we could not observe any noticeable change in the
structures of the wall modes over the time-span of our simulations. Any change
in the spatial arrangement of the modes is likely to be visible only after the
simulations are run for several times the diffusion time scale, which is
around 150 to 200 free-fall time units. Hence our solutions can be considered
to be effectively steady.

We estimate the wavelength of the wall modes by visual inspection of the
vertical velocity isosurfaces in figures 9(a)–(j). We consider only those wall
modes that are adjacent to $y=\pm\Gamma/2$ sidewalls; these walls are not
parallel to the resultant magnetic field. The estimation of the wavelength is
done as follows. We count the number of upwelling (or downwelling) plumes
adjacent to both $y=-\Gamma/2$ and $y=\Gamma/2$ sidewalls and calculate their
average as $N_{p}$. The wavelength $\lambda_{wm}$ of the wall modes is given
by

$\lambda_{wm}=\frac{\Gamma}{N_{p}}.$ (34)

Figure 11: Variations of the wavelength $\lambda_{wm}$ of the wall modes with
$R$. The wavelength tends to decrease with an increase in the horizontal
magnetic field.

We plot the estimated $\lambda_{wm}$ versus $R$ in figure 11. The figure shows
that $\lambda_{wm}$ tends to decrease as the horizontal magnetic field
increases and saturates at large $R$. Further, the wavelength corresponding to
$Ha_{z}=200$ is smaller than that for $\Ha_{z}=100$ except for $R=3$. This is
in line with the fact that the threshold wavelength decreases with increasing
vertical magnetic field strength (Chandrasekhar, 1981). It can be recalled
that a similar trend was observed in the variation of the most unstable
wavelength with $R$ in our stability analysis in § 3.1. However, it must be
kept in mind that since the Rayleigh numbers considered in our DNS are higher
than the threshold Rayleigh number above which the wall modes appear, there is
always a chance for secondary modes with smaller wavenumbers to develop.

####  3.2.2 Heat and momentum transport

We will now explore the influence of the wall modes on the heat and momentum
transport. We compute the Nusselt ($\Nu$) and the Reynolds numbers ($\Rey$)
using our numerical data as follows:

$\displaystyle\Nu$ $\displaystyle=$ $\displaystyle 1+\sqrt{\Ray\Pran}\langle
u_{z}T\rangle_{V},$ (35) $\displaystyle\Rey$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\Ray}{\Pran}}U_{rms},$ (36)

where $U_{rms}=\sqrt{\langle u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\rangle_{V}}$ is the
root mean square velocity and $\langle\cdot\rangle_{V}$ denotes volume
averaging. The second term of the right hand side of (35) is the normalized
convective heat flux. We plot the normalized heat flux and the Reynolds number
versus $R$ for all our runs in figures 12(a) and (b) respectively. The figures
show that for $R<1$, both $\Nu$ and $\Rey$ decrease with $R$. This is
consistent with our conclusion from the stability analysis in § 3.1 that for
small values of $R$, an increase in the horizontal magnetic field results in
the stabilization of the magnetoconvection system. There is, however, no clear
trend in the variations of $\Nu$ and $\Rey$ for $R>1$.

Figure 12: Plots of (a) the nondimensional convective heat flux $\Nu-1$ and
(b) the Reynolds number $\Rey$ versus $R$. For $R\leq 1$, the heat and
momentum transport decreases as $R$ is increased.

Having studied the trends of the global heat and momentum transport with $R$,
we will now explore the spatial variation of heat and momentum transport for
different regimes of $R$. In figures 13(a)–(e), we plot twice the kinetic
energy $E=0.5(u_{x}^{2}+u_{y}^{2}+u_{z}^{2})$, averaged over the $x$-$z$
plane, versus $y$, which is the direction parallel to the horizontal magnetic
field. For $R=0$, there are two sharp peaks close to $y=-2$ and $y=2$ for both
$\Ha_{z}=100$ and $\Ha_{z}=200$; these peaks correspond to the wall modes. The
bulk region lies in between these two peaks. The bulk consists of several
smaller peaks and the kinetic energy in this region is much less compared to
the near-wall regions. As $R$ is increased, the height of the near-wall peaks
decreases and of those in the bulk increases. Thus, the distribution of
kinetic energy along $y$ becomes more uniform as $R$ is increased. It is a
known fact that in magnetohydrodynamic flows, the Lorentz force tends to
modify the flow-field so as to minimize the gradients of velocity along the
direction of the magnetic field (Davidson, 2017). Thus, in our case, as $R$ is
increased, the Lorentz force generated due to the $y$-component of the
magnetic field becomes strong and suppresses the gradients of velocity along
$y$.

In figures 13(f)–(j), we plot the local convective heat flux $u_{z}T$,
averaged over the $x$-$z$ plane, along $y$. Although the volume-averaged heat
flux is positive in thermal convection, the local heat flux for $R=0$
fluctuates between positive and negative values as one proceeds along $y$.
These fluctuations even out as $R$ is increased, and for $R=3$, the local heat
flux remains positive throughout with very small gradients along $y$. Again,
this is due to the strong Lorentz forces generated by the horizontal component
of the magnetic field which suppresses the gradients of velocity along the
magnetic field’s direction.

Figure 13: Spatial distribution of kinetic energy $E$ and the convective heat
flux $u_{z}T$. Plots of 2$\langle E\rangle_{x,z}=\langle
u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\rangle_{x,z}$ for (a) $R=0$, (b) $R=0.3$, (c)
$R=1$, (d) $R=2$, and (e) $R=3$. Plots of $\langle u_{z}T\rangle_{x,z}$ for
(f) $R=0$, (g) $R=0.3$, (h) $R=1$, (i) $R=2$, and (j) $R=3$. In the above,
$\langle\cdot\rangle_{x,z}$ represents averaging over $x$-$z$ plane.

Figure 14: Spatial distribution of kinetic energy $E$ and the convective heat
flux $u_{z}T$. Plots of 2$\langle E\rangle_{y,z}=\langle
u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\rangle_{y,z}$ for (a) $R=0$, (b) $R=0.3$, (c)
$R=1$, (d) $R=2$, and (e) $R=3$. Plots of $\langle u_{z}T\rangle_{y,z}$ for
(f) $R=0$, (g) $R=0.3$, (h) $R=1$, (i) $R=2$, and (j) $R=3$. In the above,
$\langle\cdot\rangle_{y,z}$ represents averaging over $y$-$z$ plane.

In figures 14(a)–(e), we plot twice the kinetic energy, averaged over the
$y$-$z$ plane, versus $x$, which is the direction perpendicular to the
horizontal magnetic field. In the absence of horizontal magnetic field
($R=0$), the variation of $2\langle E\rangle_{y,z}$ versus $x$ is similar to
that of $2\langle E\rangle_{x,z}$ versus $y$ due to the symmetry of the
problem. Again, there are two sharp peaks corresponding to the wall modes
close to $x=-2$ and $x=2$ for both $\Ha_{z}=100$ and $\Ha_{z}=200$. However,
unlike $\langle E\rangle_{x,z}$, the values of $\langle E\rangle_{y,z}$ at the
peaks recede only marginally as $R$ is increased from 0 to 2. As the wall
modes extend fully into the bulk at $R=3$, the above peaks recede sharply with
the value of $2\langle E\rangle_{y,z}$ at these peaks being close to its
values at the peaks in the bulk. It must be noted that $\langle
E\rangle_{y,z}$ continues to fluctuate as it is varied with $x$ at $R=3$ and
does not smoothen out unlike $\langle E\rangle_{x,z}$.

Figures 14(f–j) exhibit the variations of the local convective heat flux
$u_{z}T$, averaged over the $y$-$z$ plane, along $x$. The figures show that
$\langle u_{z}T\rangle_{y,z}$ fluctuates between positive and negative values.
It is clear from the figures that for the $\Ha_{z}=100$ case, the spatial
fluctuations increase with increasing $R$. This is because as discussed
before, the gradients along the $y$-direction are reduced as $R$ is increased.
Therefore, for large $R$, if the gradients along $z$ are also small, any
quantity averaged over $y$-$z$ plane will fluctuate between the quantity’s two
extremities. The case of $R=3$ for $\Ha_{z}=100$ consists of plumes from the
opposite $y$=$\pm\Gamma$ walls merged with each other; therefore, the
gradients of the velocity, temperature, and heat flux along the $z$-direction
are small. Thus, $\langle u_{z}T\rangle_{y,z}$ exhibits strong fluctuations
along the $x$-direction at $R=3$.

For $\Ha_{z}=200$ and $R\leq 2$, $\langle u_{z}T\rangle_{y,z}$ follows a
similar trend in that its fluctuations increase with increasing $R$. However,
for $R=3$, its fluctuations get suppressed and $\langle u_{z}T\rangle_{y,z}$
is mostly positive. Now, let us recall that unlike in the case of
$\Ha_{z}=100$, the structures for the case of $\Ha_{z}=200$ are antisymmetric.
Thus, although the gradients along the $y$-direction are small, there are
fluctuations along the $z$-direction for $\Ha_{z}=200$ due to the presence of
upwelling and downwelling plumes on top of each other. An averaging over the
$y$-$z$ plane cancels out the opposing effects of the upwelling and
downwelling plumes, thus resulting in the suppression of fluctuations of
$\langle u_{z}T\rangle_{y,z}$.

Finally, we analyze the contributions by the bulk and near-wall regions to the
total kinetic energy and heat flux of the system. Towards this objective, we
define the near-wall region as follows. For $R=0$ and for both the Hartmann
numbers, we determine $\Nu_{S}$, which is the Nusselt number averaged over
successively smaller concentric volumes $S=\Gamma\times[r_{y},\Gamma-
r_{y}]\times 1$. In this definition, $r_{y}$ is the normal distance from the
$y$=$\pm\Gamma/2$ sidewalls. We plot $\Nu_{S}$ versus $r_{y}$ in figure 15.
The figure shows that $\Nu_{S}$ initially decreases with $r_{y}$ upto a point
of local minima and then begins to increase with $r_{y}$. The distance between
the point of the local minima and the sidewall is taken as the width
$\delta_{w}$ of the near-wall region. These widths are computed to be
$\delta_{w}=0.37$ for $\Ha_{z}=100$ and $\delta_{w}=0.33$ for $\Ha_{z}=200$.
Although the above values of $\delta_{w}$ were computed for $R=0$, these
values will be assumed to hold for all $R$.

The total kinetic energy $\mathcal{E}$ and the convective heat flux
$\mathcal{H}$ can be expressed as the sum of their bulk and near-wall
contributions. Therefore,

$\displaystyle\mathcal{E}$ $\displaystyle=$
$\displaystyle\mathcal{E}_{bulk}+\mathcal{E}_{nw},$ (37)
$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\mathcal{H}_{bulk}+\mathcal{H}_{nw},$ (38)

where $\mathcal{E}_{bulk}$ and $\mathcal{E}_{nw}$ are respectively the bulk
and near-wall contributions to the total kinetic energy, and
$\mathcal{H}_{bulk}$ and $\mathcal{H}_{nw}$ are respectively the bulk and
near-wall contributions to the total heat flux. The total and bulk
contributions to these quantities are defined as

$\displaystyle\mathcal{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{-\Gamma/2}^{\Gamma/2}\int_{-\Gamma/2}^{\Gamma/2}\int_{-1/2}^{1/2}(u_{x}^{2}+u_{y}^{2}+u_{z}^{2})~{}dzdydx,$
(39) $\displaystyle\mathcal{E}_{bulk}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{-\Gamma/2+\delta_{w}}^{\Gamma/2-\delta_{w}}\int_{-\Gamma/2+\delta_{w}}^{\Gamma/2-\delta_{w}}\int_{-1/2}^{1/2}(u_{x}^{2}+u_{y}^{2}+u_{z}^{2})~{}dzdydx,$
(40) $\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\int_{-\Gamma/2}^{\Gamma/2}\int_{-\Gamma/2}^{\Gamma/2}\int_{-1/2}^{1/2}u_{z}T~{}dzdydx,$
(41) $\displaystyle\mathcal{H}_{bulk}$ $\displaystyle=$
$\displaystyle\int_{-\Gamma/2+\delta_{w}}^{\Gamma/2-\delta_{w}}\int_{-\Gamma/2+\delta_{w}}^{\Gamma/2-\delta_{w}}\int_{-1/2}^{1/2}u_{z}T~{}dzdydx.$
(42)

We compute the relative strengths of the bulk and near-wall kinetic energies
and heat fluxes using our numerical data and plot them versus $R$ in figures
15(a)–(d). For small values of $R$, the bulk contribution to the total kinetic
energy and heat flux is very small (less than 10%). This is expected because
convection is completely suppressed in the bulk at small $R$. It is clear from
the figures that as $R$ increases, the bulk contributions to the total kinetic
energy and heat flux increase. This is due to the fact that the wall modes
extend further into bulk as $R$ increases. In fact, for $R>2$, the bulk and
sidewall contributions become comparable to each other because the wall modes
at such high values of $R$ extend fully into the bulk. It is also interesting
to note that for $\Ha_{z}=100$, the bulk contribution to the heat flux
decreases as $R$ is increased from 2 to 3. The reason for this anomalous
behaviour is yet to be understood.

Figure 15: Determination of the near-wall region for our analysis from the
results of our DNS for $R=0$. Variation of the Nusselt number $\Nu_{S}$
averaged over successively smaller concentric volumes with the sidewall-normal
distance $r_{y}$. The width of the near-wall region is given by the point of
the first minima (represented as black dashed vertical line for $\Ha_{z}=100$
and blue dashed vertical line for $\Ha_{z}=200$) in the $\Nu_{S}$ profile

Figure 16: Results of direct numerical simulations: Relative strengths of the
bulk and boundary layer contributions to the total kinetic energy
$\mathcal{E}$ for (a) $\Ha_{z}=100$ and (b) $\Ha_{z}=200$, and to the total
heat flux $\mathcal{F}$ for (c) $\Ha_{z}=100$ and (d) $\Ha_{z}=200$. The bulk
contributions to kinetic energy and heat flux increase with increasing $R$.

####  3.2.3 Effects of initial conditions

The direct numerical simulations of magnetoconvection for $R=3$ resulted in
different structural arrangements of the convection plumes for $\Ha_{z}=100$
and $\Ha_{z}=200$. On one hand, we obtained a symmetric arrangement of the
plumes adjacent to the opposite walls merging with each other for
$\Ha_{z}=100$. On the other hand, an antisymmetric arrangement of the plumes
was obtained for $\Ha_{z}=200$ with the upwelling and downwelling plumes on
top of each other. In this subsection, we will explore the sensitivity of the
above results to initial conditions.

We conduct two more direct numerical simulations of magnetoconvection for
$R=3$; one with $\Ha_{z}=100$, $\Ray=10^{5}$, and the other with
$\Ha_{z}=200$, $\Ray=4\times 10^{5}$. However, we change the initial
conditions in this set of simulations as follows. The results of the old
simulation of $\Ha_{z}=200$ are taken as the initial condition for the new
simulation of $\Ha_{z}=100$. In the same way, the results of the old
simulation of $\Ha_{z}=100$ are taken as the initial condition for the new
simulation of $\Ha_{z}=200$. Both these simulations were allowed to run for 20
free-fall time units after reaching the steady state. Henceforth, we refer to
the old set of simulations as IC1 and the new set of simulations of $R=3$ as
IC2.

We examine the structure of the wall-modes for the new set of simulations by
plotting the vertical velocity isosurfaces for $u_{z}=\pm 0.01$ in figures
17(a) and (c). The figures show that the simulation of $\Ha_{z}=100$ with the
new initial conditions yields antisymmetric arrangement structures with
upwelling and downwelling plumes on top of each other. This is unlike the
result of the simulation with old initial conditions which resulted in a
symmetric arrangement of structures with merged plumes. Similarly, the
simulation of $\Ha_{z}=200$ with the new initial conditions yield symmetric
structures with merged plumes, unlike the case with old initial conditions
which resulted in antisymmetric arrangement of structures. Our results
indicate that that the solution of magnetoconvection with $R=3$ for
$\Ha_{z}=100$ and 200 is non-unique, and the arrangement of the wall modes
depends on the initial conditions.

Figure 17: Results of direct numerical simulations with new sets of initial
conditions (IC2): Isosurfaces of $u_{z}=0.01$ (red) and $u_{z}=-0.01$ (blue)
for $R=3$ and (a) $Ha_{z}=100$ and (c) $\Ha_{z}=200$. For comparison, the
vertical velocity isosurfaces using results of the simulations of $R=3$ with
the old initial conditions are shown for (b) $\Ha_{z}=100$ and (d)
$\Ha_{z}=200$

Figure 18: Results of direct numerical simulations including those using the
new initial conditions (filled markers) along with the old initial conditions
(unfilled markers): Plots of (a) the nondimensional convective heat flux
$\Nu-1$ and (b) the Reynolds number $\Rey$ versus $R$. The heat and momentum
transport for the case of $\Ha_{z}=200$, $R=3$ is significantly more for the
new initial condition (IC2) compared to the old initial condition (IC1)

We will now explore the impact of the arrangement of the plumes on the global
heat and momentum transport. In figures 18(a,b), we replot the computed values
of $\Nu-1$ and $\Rey$ for our runs with the old initial conditions and also
add the plots of these quantities computed using the results of our
simulations with new initial conditions. The figures show that for
$\Ha_{z}=200$, the Reynolds and Nusselt numbers computed using the solutions
of simulations with new initial conditions are significantly higher than those
corresponding to the old initial conditions. On the other hand, for
$\Ha_{z}=100$, the Reynolds and Nusselt numbers computed using the solutions
of simulations with new initial conditions are lower, albeit marginally, than
those corresponding to the old initial conditions. The above observations
clearly indicate that the solutions that are symmetric consisting of merged
plumes have higher heat and momentum transport. This is because merged plumes
give rise to larger convection rolls which, in turn, result in more efficient
transport of heat and momentum. It is worth noting that for $\Ha_{z}=200$, the
increase of Reynolds number at $R=3$ due to the merged plumes is so
significant that its value is larger than for $R=0$.

Our results indicate that the non-unique nature of the solutions for $R=3$ has
a profound effect on the heat and momentum transport, especially for
$\Ha_{z}=200$. We conclude in the next section.

##  4 Summary and conclusions

In this paper, we systematically examined the effects of strong inclined
magnetic fields on wall-attached magnetoconvection using a combination of
linear stability analysis and direct numerical simulations. The linear
stability analysis was conducted for lower Hartmann numbers and aspect-ratio
box whereas the DNS were conducted for higher Hartmann numbers and aspect
ratio. The ratio $R=B_{y}/B_{z}$ of the outer imposed horizontal to vertical
magnetic field was varied from 0 to 3.

The linear stability analysis revealed that the critical Rayleigh number
varies non-monotonically with the relative strength of the horizontal magnetic
field. The plumes at the onset of instability get elongated along the
direction the resultant magnetic field and extend more into the bulk as $R$ is
increased. At sufficiently high $R$, the plumes extend fully into the bulk.
The linear stability analysis further revealed the existence of non-unique
solutions at the onset at low $R$. One eigensolution corresponds to a
symmetric arrangement of hot or cold plumes adjacent to the opposite
sidewalls, whereas the other eigensolution corresponds to an antisymmetric
arrangement of hot or cold plumes. For the symmetric solution, when the
opposite plumes extend sufficiently into the bulk at a large $R$, they merge
into a single large plume. Below this critical value of $R=R_{c}$, the
symmetric and antisymmetric eigensolutions overlap; however, at $R>R_{c}$,
these solutions start to diverge with the symmetric eigensolution being more
unstable than the antisymmetric solution. The critical Rayleigh number
increases with $R$ for $R<R_{c}$ and decreases with $R$ for $R>R_{c}$.

The direct numerical simulations of the fully nonlinear regime reinforced the
observation from the linear stability analysis that the wall modes get
elongated along the direction of the resultant magnetic field and extend
further into the bulk as $R$ is increased. The heat and momentum transport was
observed to decrease with $R$ for $R<1$ but not show any visible trend for
$R>1$. The fluctuations of the local heat flux and kinetic energy along the
direction of the horizontal magnetic field get suppressed as $R$ is increased.
Since the wall modes extend more into the bulk as $R$ is increased, the
relative contribution to the total heat transport and kinetic energy by the
near-wall regions decrease with an increase of $R$.

The analysis from direct numerical simulations further revealed that at least
for $R=3$, the solutions are dependent on the initial conditions. The
solutions for $R=3$ can either comprise of an antisymmetric arrangement of
upwelling and downwelling plumes on top of each other, or a symmetric
arrangement of merged upwelling or downwelling plumes. The solution with
merged plumes corresponds to higher heat and momentum transport because the
merged plumes give rise to larger convection rolls and hence more efficient
heat transfer.

Our present work provides important insights into the dynamics of wall-
attached convection under inclined magnetic fields which will be a typical
situation in view to most applications. Our work may find applications in
several industrial flows such as cooling blankets in fusion reactors. Although
we worked on a small set of parameters, we expect our results to hold for
higher Rayleigh numbers as well. In the future, we plan to conduct a similar
analysis for fluids at different Prandtl numbers.

## Acknowledgements

The authors thank M. K. Verma and M. Brynjell-Rahkola for useful discussions.
The authors acknowledge the computing time provided by Leibniz Supercomputing
Center, Garching, Germany, under the project ‘pn49ma’. The authors further
acknowledge the Computing Center of Technische Universität Ilmenau for the
resources provided to them for the linear stability computations as well as
for the postprocessing and visualization of their simulation data.

## Funding

The work of S.B. is sponsored by a Postdoctoral Fellowship of the Alexander
von Humboldt Foundation, Germany.

## Declaration of interests

The authors report no conflict of interest.

## Author ORCIDs

Shashwat Bhattacharya https://orcid.org/0000-0001-7462-7680

Thomas Boeck https://orcid.org/0000-0002-0814-7432

Dmitry Krasnov https://orcid.org/0000-0002-8339-7749

Jörg Schumacher https://orcid.org/0000-0002-1359-4536

## References

  * Akhmedagaev et al. (2020a) Akhmedagaev, R., Zikanov, O., Krasnov, D. & Schumacher, J. 2020a Rayleigh-Bénard convection in strong vertical magnetic field: flow structure and verification of numerical method. Magnetohydrodynamics 56, 157–166.
  * Akhmedagaev et al. (2020b) Akhmedagaev, R., Zikanov, O., Krasnov, D. & Schumacher, J. 2020b Turbulent Rayleigh-Bénard convection in a strong vertical magnetic field. J. Fluid Mech. 895, R4.
  * Aurnou & Olson (2001) Aurnou, J. M. & Olson, P. L. 2001 Experiments on Rayleigh–Bénard convection, magnetoconvection and rotating magnetoconvection in liquid gallium. J. Fluid Mech. 430, 283–307.
  * Bhattacharya et al. (2023) Bhattacharya, S., Boeck, T., Krasnov, D. & Schumacher, J. 2023 Effects of strong fringing magnetic fields on turbulent thermal convection. J. Fluid Mech. 964, A31.
  * Burr & Müller (2001) Burr, U. & Müller, U. 2001 Rayleigh–Bënard convection in liquid metal layers under the influence of a vertical magnetic field. Phys. Fluids 13, 3247–3257.
  * Burr & Müller (2002) Burr, U. & Müller, U. 2002 Rayleigh–Bénard convection in liquid metal layers under the influence of a horizontal magnetic field. J. Fluid Mech. 453, 345–369.
  * Busse (2008) Busse, F. H. 2008 Asymptotic theory of wall-attached convection in a horizontal fluid layer with a vertical magnetic field. Phys. Fluids 20 (2), 024102.
  * Busse & Clever (1983) Busse, F. H. & Clever, R. M. 1983 Stability of convection rolls in the presence of a horizontal magnetic field. J. Theor. Appl. Mech. 2, 495–502.
  * Chandrasekhar (1981) Chandrasekhar, S. 1981 Hydrodynamic and Hydromagnetic Stability. Oxford: Dover publications.
  * Chillà & Schumacher (2012) Chillà, F. & Schumacher, J. 2012 New perspectives in turbulent Rayleigh-Bénard convection. Eur. Phys. J. E 35 (7), 58.
  * Cioni et al. (2000) Cioni, S, Chaumat, S & Sommeria, J 2000 Effect of a vertical magnetic field on turbulent Rayleigh-Bénard convection. Phys. Rev. Lett. 62 (4), R4520–R4523.
  * Davidson (1999) Davidson, Peter A. 1999 Magnetohydrodynamics in materials processing. Annu. Rev. Fluid Mech. 31 (1), 273–300.
  * Davidson (2017) Davidson, Peter A. 2017 An Introduction to Magnetohydrodynamics, 2nd edn. Cambridge: Cambridge University Press.
  * Fauve et al. (1981) Fauve, S., Laroche, C. & Libchaber, A. 1981 Effect of a horizontal magnetic field on convective instabilities in mercury. J. Physique Lett. 42 (21), L455.
  * Grannan et al. (2022) Grannan, A. M., Cheng, J. S., Aggarwal, A., Hawkins, E. K., Xu, Y., Horn, S., Sánchez-Álvarez, J. & Aurnou, J. M. 2022 Experimental pub crawl from Rayleigh-Bénard to magnetostrophic convection. J. Fluid Mech. 939, R1.
  * Houchens et al. (2002) Houchens, B. C., Witkowski, L. Martin & Walker, J. S. 2002 Rayleigh–bénard instability in a vertical cylinder with a vertical magnetic field. J. Fluid Mech. 469, 189–207.
  * Hurlburt et al. (1996) Hurlburt, N., Matthews, P. & Proctor, M. 1996 Nonlinear compressible convection in oblique magnetic fields. Ap. J. 457, 933.
  * Kelley & Sadoway (2014) Kelley, D. & Sadoway, D. R. 2014 Mixing in a liquid metal electrode. Phys. Fluids 26, 057102.
  * Kelley & Weier (2018) Kelley, D. & Weier, T. 2018 Fluid Mechanics of Liquid Metal Batteries. Appl. Mech. Rev. 70, 020801.
  * King & Aurnou (2015) King, E. M. & Aurnou, J. M. 2015 Magnetostrophic balance as the optimal state for turbulent magnetoconvection. Proc. Natl. Acad. Sci. U.S.A. 112, 990–994.
  * Krasnov et al. (2011) Krasnov, D., Zikanov, O. & Boeck, T. 2011 Comparative study of finite difference approaches in simulation of magnetohydrodynamic turbulence at low magnetic reynolds number. Comput. Fluids 50 (1), 46 – 59.
  * Liu et al. (2018) Liu, W., Krasnov, D. & Schumacher, J. 2018 Wall modes in magnetoconvection at high Hartmann numbers. J. Fluid Mech. 849, R2.
  * Lohse & Xia (2010) Lohse, Detlef & Xia, Ke-Qing 2010 Small-scale properties of turbulent Rayleigh–Bénard convection. Annu. Rev. Fluid Mech. 42 (1), 335–364.
  * Lyubimov et al. (2010) Lyubimov, D. V., Burnysheva, A. V., Benhadid, H., Lyubimova, T. P. & Henry, D. 2010 Rotating magnetic field effect on convection and its stability in a horizontal cylinder subjected to a longitudinal temperature gradient. J. Fluid Mech. 664, 108–137.
  * McCormack et al. (2023) McCormack, Matthew, Teimurazov, Andrei, Shishkina, Olga & Linkmann, Moritz 2023 Wall mode dynamics and transition to chaos in magnetoconvection with a vertical magnetic field, arXiv: 2308.15165.
  * Mistrangelo et al. (2020) Mistrangelo, C., Bühler, L. & Klüber, V. 2020 Three-dimensional magneto convective flows in geometries relevant for DCLL blankets. Fusion Eng. Des. 159, 111686\.
  * Mistrangelo et al. (2021) Mistrangelo, C., Bühler, L., Smolentsev, S., Klüber, V., Maione, I. & Aubert, J. 2021 MHD flow in liquid metal blankets: Major design issues, MHD guidelines and numerical analysis. Fusion Eng. Des. 173, 112795.
  * Nakagawa (1957) Nakagawa, Y 1957 Experiments on the inhibition of thermal convection by a magnetic field. Proc. R. Soc. Lond. 240, 108–113.
  * Nicoski et al. (2022) Nicoski, J. A., Yan, M. & Calkins, M. A. 2022 Quasistatic magnetoconvection with a tilted magnetic field. Phys. Rev. Fluids 7, 043504.
  * Peyret (2002) Peyret, R 2002 Spectral Methods for Incompressible Viscous Flows. New York: Springer.
  * Priede et al. (2010) Priede, Jānis, Aleksandrova, Svetlana & Molokov, Sergei 2010 Linear stability of Hunt’s flow. J. Fluid Mech. 649, 115–134.
  * Roberts (1967) Roberts, P. H. 1967 An introduction to magnetohydrodynamics. London: Longmans.
  * Shen & Zikanov (2016) Shen, Y. & Zikanov, O. 2016 Thermal convection in a liquid metal battery. Theor. Comput. Fluid Dyn. 30, 275–294.
  * Tasaka et al. (2016) Tasaka, Y., Igaki, K., Yanagisawa, T., Vogt, T., Zürner, T. & Eckert, S. 2016 Regular flow reversals in Rayleigh-Bénard convection in a horizontal magnetic field. Phys. Rev. E 93, 043109.
  * Teimurazov et al. (2023) Teimurazov, Andrei, McCormack, Matthew, Linkmann, Moritz & Shishkina, Olga 2023 Unifying heat transport model for the transition between buoyancy-dominated and lorentz-force-dominated regimes in quasistatic magnetoconvection, arXiv: 2308.01748.
  * The MathWorks Inc. (2022) The MathWorks Inc. 2022 Matlab version: 9.13.0 (r2022b).
  * Verma (2018) Verma, Mahendra K. 2018 Physics of Buoyant Flows: From Instabilities to Turbulence. Singapore: World Scientific.
  * Verma (2019) Verma, M. K. 2019 Energy trasnfers in Fluid Flows: Multiscale and Spectral Perspectives. Cambridge: Cambridge University Press.
  * Vogt et al. (2018) Vogt, T., Ishimi, W., Yanagisawa, T., Tasaka, Y., Sakuraba, A. & Eckert, S. 2018 Transition between quasi-two-dimensional and three-dimensional Rayleigh-Bénard convection in a horizontal magnetic field. Phys. Rev. Fluids 3, 013503\.
  * Vogt et al. (2021) Vogt, T., Yang, J.-C., Schindler, F. & Eckert, S. 2021 Free-fall velocities and heat transport enhancement in liquid metal magneto-convection. J. Fluid Mech. 915, A68.
  * Weiss & Proctor (2014) Weiss, N. O. & Proctor, M. R. E. 2014 Magnetoconvection. Cambridge: Cambridge University Press.
  * Yan et al. (2019) Yan, M., Calkins, M., Maffei, A., Julien, K., Tobias, S. & Marti, P. 2019 Heat transfer and flow regimes in quasi-static magnetoconvection with a vertical magnetic field. J. Fluid Mech. 877, 1186–1206.
  * Yanagisawa et al. (2013) Yanagisawa, T., Hamano, Y., Miyagoshi, T., Yamagishi, Y., Tasaka, Y. & Takeda, Y. 2013 Convection patterns in a liquid metal under an imposed horizontal magnetic field. Phys. Rev. E 88 (6), 063020.
  * Zürner et al. (2020) Zürner, T., Schindler, F., Vogt, T., Eckert, S. & Schumacher, J. 2020 Flow regimes of Rayleigh-Bénard convection in a vertical magnetic field. J. Fluid Mech. 894, A21.