# Shear flows and their suppression at large aspect ratio. Two-dimensional
simulations of a growing convection zone

J. R. Fuentes  A. Cumming  Department of Physics and McGill Space Institute,
McGill University, 3600 rue University, Montreal, QC H3A 2T8, Canada

###### Abstract

We investigate the onset and evolution of zonal flows in a growing convective
layer when a stably-stratified fluid with a composition gradient is cooled
from above. This configuration allows the study of zonal flows for a wide
range of values of the Rayleigh number, $Ra$, and aspect ratio of the
convection zone within a given simulation. We perform a series of 2D
simulations using the Boussinesq approximation, with  aspect ratio of the
computational domain between $1$ and $5$, and Prandtl number $Pr=0.1$, 0.5, 1,
and $7$. We find that for square domains zonal flows appear when the aspect
ratio of the convective layer is smaller than two, and the evolution of the
system depends on the Prandtl number. For $Pr\leq 1$, the fluid experiences
bursts of convective transport with negligible convective transport between
bursts. The magnitude and frequency of the bursts are smaller at low $Pr$,
which suggests that the bursting regime is stronger in a narrow range around
$Pr=1$, as observed in previous studies of thermal convection. For $Pr=7$, the
structure of the flow consists of tilted convective plumes, and the convective
transport is sustained at all times. In wider domains, the aspect ratio of the
convective zone is always much larger than two and zonal flows do not appear.
These results confirm and extend to fluids with stable composition gradients
previous findings on thermal convection. The fact that zonal flows can be
avoided by using computational domains with large aspect ratios opens up the
possibility of 2D studies of convective overshoot, layer formation and
transport properties across diffusive interfaces.

##  I Introduction

Large scale horizontal flows ( hereafter) arise in many geophysical and
astrophysical fluids, particularly in convective systems where the fluid
motions are approximately two-dimensional,  such as rotating flows in
planetary atmospheres, where the Taylor-Proudman effect suppresses fluid
motions parallel to the axis of rotation [1, 2]. These flows are of great
interest because of their chaotic non-linear behavior and also because of the
effects they have on mixing and transport properties.

Thermal convection in two-dimensions at large Rayleigh number ($Ra$111See
Sect. II for a definition.) has been used extensively as the canonical model
to study zonal flows  in convection. A remarkable result from previous studies
of zonal flows is that the nature of the fluid and the transport properties
depend strongly on the Prandtl number ($Pr$1) [3, 4, 5, 6]. In particular, it
has been shown that once the convective flow is affected by large-scale
horizontal motions at large $Ra$, for low $Pr$ ($\lesssim 1$) it undergoes
strong oscillations and heat transport occurs in chaotic bursts, whereas for
higher $Pr$ the flow does not burst and  vertically-sheared thermal plumes
dominate the structure of the flow at all times. In both regimes, once the
zonal flow sets in, its net effect is the decrease of the vertical heat
transport. Interestingly, bursts and sheared convective plumes have also been
found in two-dimensional simulations of fingering convection [e.g., 7, 8, 9],
and as with zonal flows in thermal convection, they decrease the vertical
transport of heat and solute.

Thompson [10] suggested that zonal flows in two-dimensional convection are
sustained against viscous diffusion by Reynold stresses. By perturbing
convective cells with a small shear flow, the shear advects convective plumes,
tilting them with respect to the vertical. The Reynolds stresses associated
with the tilted plumes transport momentum across the vertical, enhancing the
shear flow and therefore the amplitude of the shear perturbations. This
mechanism is called tilting instability and energy is supplied directly from
the convective scales to the zonal flow. Although a different mechanism has
been proposed to explain the zonal flow [see, e.g., 11, who proposed that
zonal flow generation can be due to an energy transfer from small to large
scale motions], the tilting mechanism has been shown to be operating in
numerical simulations and theoretical models [12, 3].

Zonal flows have also been observed in old laboratory experiments of Rayleigh-
Benard convection in water confined to a horizontal layer [13, 14]. However,
the properties and the different regimes described above have not been
observed in later experiments nor in three-dimensional numerical simulations
[15, 16], unless an external mechanism  induces them [17, 18]. Furthermore,
Hartlep  _et al._ [15] demonstrated that convective rolls of large horizontal
extent can arise at large $Ra$ in three-dimensional simulations, and suggested
them as a better interpretation of the observations, rather than the zonal
flows that appear in two-dimensional convection.

Besides the importance of characterizing the morphology and transport
properties of zonal flows in convection, the effects discussed above raised
the question of whether two-dimensional simulations are appropriate to study
convection problems in general. Since three-dimensional simulations are
computationally expensive, two-dimensional simulations are often adopted to
study convection in different systems. By comparing the outcomes of 2D and 3D
simulations of fingering convection at low $Pr$, Garaud and Brummell [7]
pointed out that 2D simulations are not appropriate since the shear effects of
artificial zonal flows significantly decrease the vertical transport, and
change the structure and morphology of the flow. However, recent work by
Fitzgerald and Farrell [5] and Wang  _et al._ [19] have shown that for the
case of pure thermal convection, zonal flows are not sustained in two-
dimensional simulations using domains of large aspect ratio.

Whether zonal flows appear in fluids with stable composition gradients has not
been explored yet. If they do, it is not clear if they can be suppressed at
large aspect ratio as for the case of thermal convection. The answers to these
questions are important for example for two-dimensional studies of turbulent
mixing and layer formation in double-diffusive convection.

In this work, we report two-dimensional simulations of convection driven by a
constant heat flux at the top boundary in a stable fluid with a solute
gradient. This configuration is particularly useful to study the onset and
evolution of zonal flows. The solute gradient stabilizes the fluid against
overturning convection in the whole fluid domain, leaving a convective layer
whose thickness increases with time. The aspect ratio of the convective layer
in a given simulation decreases in time as the thickness of the layer grows.
This allows the study of convection for a wide range of $Ra$ and aspect ratio
using a smaller grid of simulations. Further, it provides an opportunity to
study the transition to sheared convection once the zonal flow arises. We
perform simulations at $Pr$ ranging from 0.1 to 7 in order to explore the
bursting and non-bursting regimes. We also perform simulations with fluid
domains of different aspect ratio to confirm that zonal flows and their shear
effects do not appear at large aspect ratio, as observed in previous studies
of thermal convection.

The paper is organised as follows. Sect. II presents the model and the
numerical code used to perform the simulations. In Sect. III.1 we study the
onset of zonal flows and the bursting and non-bursting regimes in simulations
with aspect ratio of one, i.e., square domains ($L=H$). In Sect. III.2 we show
that strong horizontal flows and their effects vanish when the width of the
domain is increased ($L\geq 2H$). We conclude in Sect. IV.

##  II Model and numerical method

We study the onset and evolution of zonal flows in a two-dimensional
convective layer that grows inward by incorporating fluid from below. These
simulations are based on our previous study [20] that was focused on the rate
at which the convection zone grows inwards. However, in the previous work we
deliberately excluded those simulations that developed shearing convection. In
the simulations, we model the fluid under the Boussinesq approximation [21],
which is appropriate when the density variations are small respect to the
background density ($\rho/\rho_{0}\ll 1$). The fluid domain is a Cartesian box
in two-dimensions ($x$,$z$) of width $L$ and height $H$, with periodic
boundary conditions  in the horizontal direction. The top and bottom
boundaries are impermeable and stress-free, with no composition flux through
them, no heat flux at the bottom, and a constant heat flux $F_{0}$ at the top.
Initially, the fluid starts with a uniform temperature and a constant solute
gradient $d\overline{S}_{0}/dz=-\delta S_{0}/H<0$, such that the solute
concentration is two times larger at the bottom compared to the top. The
evolution of the system is driven by the imposed heat flux $F_{0}$ at the top
boundary, which cools the fluid from the top and quickly forms a convective
layer that grows inwards.

As we mentioned above, this numerical set-up is ideal to study the onset of
the zonal flow in the simulations. First, two dimensionality prevents
transverse perturbations that can reduce horizontal fluid motions. Second,
periodic boundary conditions on the side boundaries do not confine the fluid
in the horizontal direction. Third, the free-slip top and bottom boundaries
apply no shear stresses to slow down the horizontal flow. Further, since the
onset of the zonal flow depends on the aspect ratio of the convective layer,
the time-dependent nature of the convective layer thickness is ideal to study
the critical aspect ratio at which the zonal flow arises.

Following Sect. 2 in [20], we non-dimensionalize the variables using as
characteristic length and time scales the box height, $H$, and the thermal
diffusion time across the box, $H^{2}/\kappa_{T}$ (where $\kappa_{T}$ is the
thermal diffusivity). Solute concentration is measured in units of the initial
solute contrast across the box, $\delta S_{0}$, and temperature is in units of
the initial solute contrast, $\beta\delta S_{0}/\alpha$ (here $\beta$ and
$\alpha$ are the coefficients of solute and thermal contraction-expansion,
respectively). The resulting dimensionless equations are

$\displaystyle\nabla\cdot\bm{v}=0\,,$ (1) $\displaystyle\dfrac{\partial
T}{\partial t}=-(\bm{v}\cdot\nabla)\,T+\nabla^{2}T\,,$ (2)
$\displaystyle\dfrac{\partial S}{\partial
t}=-(\bm{v}\cdot\nabla)\,S+\tau\nabla^{2}S\,,$ (3)
$\displaystyle\dfrac{\partial\bm{v}}{\partial
t}=-(\bm{v}\cdot\nabla)\,\bm{v}-\nabla
P+\mathcal{R}_{T}Pr\left(\dfrac{F_{0}}{F_{\rm
crit}}\right)^{-1}\bm{\hat{z}}+Pr\nabla^{2}\bm{v}\,,$ (4)

with boundary conditions

$\displaystyle w\,\big{|}_{z=0,1}=0\,,\hskip 5.69046pt\dfrac{\partial
u}{\partial z}\,\bigg{|}_{z=0,1}=0\,,\hskip 5.69046pt\dfrac{\partial
S}{\partial z}\,\bigg{|}_{z=0,1}=0\,,$ (5) $\displaystyle\dfrac{\partial
T}{\partial z}\,\bigg{|}_{z=0}=0\,,\hskip 5.69046pt\dfrac{\partial T}{\partial
z}\,\bigg{|}_{z=1}=-\dfrac{F_{0}}{F_{\rm crit}}\,.$ (6)

Table 1: Dimensionless parameters used in the simulations. The first and second columns correspond to the diffusivity ratio and Prandtl number, respectively. The third column contains the input cooling flux at the top boundary, and the last two columns contain the modified Rayleigh numbers for temperature and solute, respectively. $\\#$ | $\tau$ | $Pr$ | $F_{0}/F_{\rm crit}$ | $\mathcal{R}_{T}$ | $\mathcal{R}_{S}$
---|---|---|---|---|---
1 | 0.1 | 0.1 | 5.4 | $4\times 10^{12}$ | $7.5\times 10^{11}$
2 | 0.1 | 0.1 | 10.8 | $8\times 10^{12}$ | $7.5\times 10^{11}$
3 | 0.1 | 0.5 | 5.4 | $8\times 10^{11}$ | $1.5\times 10^{11}$
4 | 0.1 | 0.5 | 10.8 | $1.6\times 10^{12}$ | $1.5\times 10^{11}$
5 | 0.1 | 1 | 5.4 | $4\times 10^{11}$ | $7.5\times 10^{10}$
6 | 0.1 | 1 | 10.8 | $8\times 10^{11}$ | $7.5\times 10^{10}$
7 | 0.1 | 7 | 5.4 | $5.76\times 10^{10}$ | $1.06\times 10^{10}$
8 | 0.1 | 7 | 10.8 | $1.15\times 10^{11}$ | $1.06\times 10^{10}$

The dimensionless parameters that control the simulations are $F_{0}/F_{\rm
crit}$, the Prandtl number ($Pr$), the diffusivity ratio ($\tau$) and a
modified Rayleigh number ($\mathcal{R}_{T}$), defined respectively as

$\displaystyle\dfrac{F_{0}}{F_{\rm
crit}}=F_{0}\left(k\dfrac{\beta}{\alpha}\dfrac{\delta
S_{0}}{H}\right)^{-1}\,,$ (7) $\displaystyle Pr=\dfrac{\nu}{\kappa_{T}}\,,$
(8) $\displaystyle\tau=\dfrac{\kappa_{S}}{\kappa_{T}}\,,$ (9)
$\displaystyle\mathcal{R}_{T}=\dfrac{\alpha
gH^{3}}{\kappa_{T}\nu}\left(\dfrac{F_{0}H}{k}\right),$ (10)

where $k$ is the thermal conductivity, $\nu$, and $\kappa_{S}$ are the viscous
and solute diffusivity, respectively, and $g$ is the magnitude of the
acceleration due to gravity. As we note in [20], the product
$\mathcal{R}_{T}(F_{0}/F_{\rm crit})^{-1}$ can be re-written as

$\displaystyle\mathcal{R}_{T}\left(\dfrac{F_{0}}{F_{\rm
crit}}\right)^{-1}=\mathcal{R}_{S}=\dfrac{\beta gH^{3}\delta
S_{0}}{\kappa_{T}\nu},$ (11)

which resembles the traditional Rayleigh-number

$Ra=\dfrac{\alpha gH^{3}\delta T_{0}}{\kappa_{T}\nu}\,,$ (12)

but for solute. The parameters used in our simulations are given in Table 1.

The governing differential equations are solved with the Dedalus spectral code
[22] on a Chebyshev (vertical) and Fourier (horizontally-periodic) domain in
which the physical grid dimensions are 3/2 the number of modes. Based on a
resolution study, we use 1024 modes in each direction.

##  III Analysis and results

###  III.1 Zonal flow and sheared convection in square domains

As soon as the heat flux at the top boundary turns on, a thermal boundary
layer develops and becomes convective. The recently-formed convective layer is
composed of many convective cells, with horizontal size approximately the
height of the layer. As the layer increases its thickness by incorporating
fluid from below, both the aspect ratio of the flow and the number of
convective cells within the layer decrease. Once the aspect ratio of the
growing convective layer reaches a critical size, convective plumes become
vertically-sheared by a large scale horizontal flow. We find the critical
aspect ratio lies between $1.2-1.7$, being smaller at $Pr=0.1$ (Fig. 1a).

Figure 1: Panel (a): Ratio between the width of the domain, $L$, and the
thickness of the convective layer, $h$ (i.e., the aspect ratio of the flow
within the convective layer) as a function of time. The dashed horizontal
lines denote the value of the aspect ratio when the zonal flow sets in. Panel
(b): Rayleigh number as a function of time. Once the zonal flow sets in
(vertical lines), the magnitude of $Ra$ varies slowly with time (flat region).
In both panels, the results are shown for simulations at
$F_{0}/F_{\mathrm{crit}}=10.8$. Colors and symbols distinguish between
simulations at different $Pr$, as shown in the legends. The rest of the
simulations behave in a similar way but on different time scales and
magnitudes.

In this problem convection is driven by the temperature difference across the
thermal boundary layer due to the imposed heat flux at the top, and the
convective layer grows in time. This means that the Rayleigh number $Ra$ grows
in time and its magnitude is set by the thickness of the convective layer. We
find that once the zonal flow sets in, the convective layer stops growing and
$Ra$ saturates to a roughly constant value, $Ra\sim 10^{8},\,5\times
10^{8},\,10^{9},$ and $5\times 10^{9}$ for $Pr=7,\,1,\,0.5$, and 0.1,
respectively. To illustrate this behaviour, Fig. 1(b) shows the temporal
evolution of $Ra$ for the cases using $F_{0}/F_{\rm crit}=10.8$.

The stalling of the convective layer growth is explained by the reduced
vertical transport of sheared convective plumes. In a normal situation absent
of strong horizontal flows, kinetic energy in the convective plumes or eddies
is used to lift and mix fluid from below, increasing the convection zone
thickness [20, 23]. However, in the situation considered here, the zonal flow
disperses the convective plumes and takes energy from them. As a result, the
flow near the convective boundary does not have enough energy to lift and mix
fluid from below, thereby the convection zone does not grow. The latter can be
inferred from the horizontally-averaged flux profiles considering just the
contribution from convective motions, which in our dimensionless form is given
by $\overline{F}^{\,\,\rm conv}_{H}=\overline{wT}$, where $\bar{\cdot}$
denotes the average over the horizontal direction. We observe in Fig. 2 that
after the zonal flow arises ($t\geq 0.01\,t_{\rm diff}$), the convective heat
flux near the interface decreases. Further, the position of the interface
slightly changes with time, meaning that the rate at which the outer
convection zone grows becomes smaller. We note how the zonal flow affects the
shape of the convective flux profiles at $t=0.012\,t_{\rm diff}$, changing
from being roughly linear within the outer convection zone (the expected
profile not affected by the zonal flow) to a profile with roughly two
different cooling rates (smaller near the bottom of the layer).

Figure 2: Horizontally-averaged profiles of the convective flux divided by the
imposed heat flux, $\overline{F}^{\,\,\rm conv}_{H}/F_{0}$. Results are shown
at different times for the case at $Pr=1$, $F_{0}/F_{\rm crit}=10.8$. The
vertical lines set the position of the convective boundary at the times when
profiles are shown. Colors distinguish between different times. We recall that
for this case, from $t\geq 0.01\,t_{\rm diff}$ the fluid is affected by the
zonal flow.

Figure 3 shows profiles of the $x$-component of the flow velocity normalized
to its maximum magnitude, $\overline{u}/\max|\overline{u}|$. The profiles are
computed after the zonal flow appears. As expected, since the boundaries do
not impose stress on the flow, we find the maximum of the flow at the top
boundary. The horizontal flow at the bottom of the convective layer, which is
oppositely directed to the flow at the top, reaches about 0.4-0.5 of the
maximum flow velocity, suggesting that there must be transport of momentum
across the convective boundary. We do not find any trend regarding the
direction of the zonal flow, in some of the simulations positive (negative)
momentum is transported upward (downward), and vice-versa. However, we find
for all the cases that once the zonal flow sets in, its direction is never
reversed during the simulation.

Figure 3: Vertical profiles of the horizontal flow velocity, $\overline{u}$,
normalized to its maximum magnitude, $\max|\overline{u}|$. Colors and line-
styles distinguish between different $Pr$. All the cases correspond to
simulations using $F_{0}/F_{\rm crit}=10.8$. The horizontal lines denote the
position of the convective boundary. Note that for $Pr=7$, there is also
contribution from the advection of secondary convective layers in the fluid
(see Figs. 7c and 8c-d), and therefore the horizontal flow velocity does not
decay to zero towards the bottom of the box.

Figure 4: Panels (a)-(d): Temporal evolution of the horizontal and vertical
rms velocity (dark and light colors) for different $Pr$ using $F_{0}/F_{\rm
crit}=5.4$. Note that panels do not share the same scale in the time axis. The
time series for the cases using $F_{0}/F_{\rm crit}=10.8$ are similar but the
instability appears earlier and the magnitude of $\langle u^{2}\rangle^{1/2}$
is higher.

The onset of the zonal flow can be observed in the time series of the
horizontal and vertical rms velocities, $\langle u^{2}\rangle^{1/2}$ and
$\langle w^{2}\rangle^{1/2}$, respectively, where $\langle\cdot\rangle$
denotes the average over the entire domain. As an example, we show in Fig. 4
the temporal evolution of $\langle u^{2}\rangle^{1/2}$ and $\langle
w^{2}\rangle^{1/2}$, for all cases using $F_{0}/F_{\rm crit}=5.4$. We observe
that during the early evolution, both rms velocities increase with time having
roughly the same magnitude. This is expected since in this stage the
convection zone grows and the fluid flow is dominated by cellular motions.
However, once the instability that triggers the zonal flow arises, there is a
large increase in $\langle u^{2}\rangle^{1/2}$, and a substantial decrease in
$\langle w^{2}\rangle^{1/2}$, meaning that a significant fraction of the work
done by buoyancy forces is transformed to kinetic energy but enhancing mainly
the horizontal fluid motions. Note that for $Pr\leq 1$, the time series of the
rms velocities exhibit oscillations or bumps, whereas for $Pr=7$ it does not.
The difference in the behaviour of the rms velocities distinguish between the
bursting and non-bursting regime of the system.

The bursting and non-bursting regimes have substantial differences in the
vertical transport. These differences are more clear when we look into the
time series of the heat flux averaged over the whole domain, which in
dimensionless form is given by

$\langle F_{H}\rangle=\langle wT\rangle-\langle dT/dz\rangle\,,$ (13)

where the first and second terms correspond to the convective and diffusive
heat flux, respectively.

Figure 5: Time series of the averaged heat flux (divided by the imposed flux
at the top), $\langle F_{H}\rangle/F_{0}$ for $Pr=0.1,\,0.5,\,1$, and 7
(panels a, b, c, and d, respectively). Colors distinguish between
$F_{0}/F_{\mathrm{crit}}=5.4$ (dark) and $F_{0}/F_{\mathrm{crit}}=10.8$
(light). Note that panels do not share the same scale in the time axis. In all
panels the dashed lines corresponds to $\langle F_{H}\rangle/F_{0}=0.5$, i.e.,
the expected averaged ratio if the whole fluid cools at a constant rate.

Figure 5 shows time series of the ratio $\langle F_{H}\rangle/F_{0}$ for all
our simulations. We observe for all the cases that as the convection zone
grows, the convective contribution to the heat flux increases with time and
dominates the magnitude of $\langle F_{H}\rangle/F_{0}$. Once the zonal flow
arises and is strong enough to disperse convective plumes and reduce the
vertical kinetic energy, $\langle F_{H}\rangle/F_{0}$ decays. The subsequent
evolution of the heat flux is different depending on $Pr$. For the cases
$Pr\leq 1$, the vertical transport occur through discrete bursts whose
intensity and frequency increase with $Pr$. Each burst is separated by a
quiescent phase in which $\langle F_{H}\rangle/F_{0}\approx 15-20\%$ of its
value before the onset of the zonal flow.

(a) (a) Quiescent phase

(b) (b) Pre-burst phase

(c) (c) Burst phase

(d) (d) Quiescent phase

Figure 6: Instantaneous snapshots of the temperature field during a single
burst, for the case $Pr=1$, $F_{0}/F_{\rm crit}=10.8$. Blue and red colors
represent cold and hot fluid, respectively. The color scale is not the same is
all panels.

Figure 6 shows the different stages of a burst. In the quiescent phase (a),
convective plumes are constantly dispersed by the zonal flow. This flow
suppresses convective instabilities until it has decayed sufficiently for
convective plumes to appear again (b) and the fluid suddenly overturns (c).
The fluid is energized by buoyancy and circulation motions, but once again,
kinetic energy is transferred to the zonal flow, convective plumes get
dispersed, and a new quiescent phase begins (d).

Figure 7: Horizontally-averaged profiles of the ratio between the heat flux
and the imposed flux, $\overline{F}_{H}/F_{0}$. Panel (a) shows the profiles
during the quiescent phase, whereas Panel (b) shows the profiles for a
particular burst during the active phases. Both panels show results for a
typical burst in the simulation $Pr=1$, $F_{0}/F_{\rm crit}=10.8$. Panel (c)
shows the profiles during the sheared convection phase observed in all
simulations at $Pr=7$. The results are shown for the case $Pr=7$,
$F_{0}/F_{\rm crit}=5.4$ at $t=0.065$. In all panels, the green, blue, and
orange lines correspond to the total, advective, and diffusive contribution to
the flux, respectively. The dotted line corresponds to the expected profile if
the whole fluid cools down at a constant rate, i.e.,
$\overline{F}_{H}/F_{0}=z$.

Note that although $\langle F_{H}\rangle/F_{0}$ decays during the quiescent
phase, its value is not negligible since the averaged heat flux across the box
is more than 10$\%$ of the imposed heat flux in all the cases. The reason for
this can be explained using the flux profiles during the quiescent phase (see
Fig. 7a). We observe that near the top boundary there are still convective
motions that contribute to heat transport. This is expected since the fluid
surrounding the top boundary is constantly cooling down by the imposed flux,
thereby it has a permanent energy source to undergo convection. We also
observe a smaller contribution to heat transport due to convective motions at
the bottom of the outer convective layer ($z\approx 0.4$) and below it due to
secondary convective layers at $z\approx 0.25$ (Fig. 8a-b). We recall that the
flux time series in Fig. 5 take into account the flux averaged over the whole
box, however, even if we just consider the average over the outer convective
layer, $\langle F_{H}\rangle/F_{0}$ would be still non-negligible. On the
contrary, during the bursting phase as soon as the zonal flow weakens, the
much colder fluid at the top sinks catastrophically in the way of a Rayleigh-
Taylor plume, increasing significantly the heat flux (Fig. 7b).

Finally, the cases $Pr=7$ are different, we do not observe bursts and rather
than cellular motions, convective plumes vertically sheared dominate the flow
within the convective layer. In those cases, the contribution of the
convective flux to the total flux is significant at all times (Fig. 7c).
Further, as for the cases $Pr\leq 1$, the diffusive and convective flux
profiles in Fig. 7c (and also the snapshots in Fig. 8a-b) show that secondary
convective layers form in the fluid, being responsible of the subsequent
increase in $\langle F_{H}\rangle/F_{0}$.

Since this work is focused on zonal flows and their properties, an analysis of
layer formation and its evolution is going to be presented in a future paper.
However, it is worth mentioning that the zonal flow is able to destabilize the
convective boundary to Kelvin-Helmholtz instabilities (clearly seen in Fig.
8c), causing strong mixing and eventually forming secondary convective layers
[this mechanism is known as thermohaline shear instability, see, e.g., 24,
25].

(a) (a) $Pr=1$, $F_{0}/F_{\rm crit}=10.8$, $t=0.0190\,t_{\rm diff}$

(b) (b) $Pr=1$, $F_{0}/F_{\rm crit}=10.8$, $t=0.0198\,t_{\rm diff}$

(c) (c) $Pr=7$, $F_{0}/F_{\rm crit}=5.4$, $t=0.06\,t_{\rm diff}$ (d) (d)
$Pr=7$, $F_{0}/F_{\rm crit}=5.4$, $t=0.07\,t_{\rm diff}$

Figure 8: Instantaneous snapshots of the temperature field for the cases
$Pr=1$, $F_{0}/F_{\rm crit}=10.8$ (panels and a and b), and $Pr=7$,
$F_{0}/F_{\rm crit}=5.4$ (panels c and d). Blue and red colors represent cold
and hot fluid, respectively. The color scale is not the same for all panels.
The times were chosen to show that secondary layers form and contribute the
the averaged heat transport (in the quiescent phase for the case of panels a
and b, and during the non-bursting regime for the case of panels c and d).

###  III.2 Suppression of zonal flows in wider domains

Figure 9: Panel (a): Temporal evolution of the horizontal and vertical rms
velocities ($\langle u^{2}\rangle^{1/2}$ and $\langle w^{2}\rangle^{1/2}$,
using solid and dotted lines, respectively), for different domain width $L$
(using different colors). Panel (b): Time series of the averaged heat flux
divided by the imposed flux at the top ($\langle F_{H}\rangle/F_{0}$), for
simulations using different domain widths ($L=H$ and $L=4H$, as shown in the
legends). In both panels, the results correspond to experiments with $Pr=1$,
$F_{0}/F_{\rm crit}=5.4$.

We have shown that for square domain simulations (i.e., $L/H=1$), the zonal
flow always appears when the aspect ratio of the convective layer, $L/h$, is
smaller than approximately two (Fig. 1a). For illustrative purposes, we show
in Fig. 9(a) the effect of increasing the width of the domain on the evolution
of the rms velocities for the case $Pr=1$, $F_{0}/F_{\rm{crit}}=5.4$. We find
that for $L/H\geq 2$, both the horizontal and vertical rms velocities increase
gradually in time, without the instabilities observed for $L=H$. Further, the
curves lie on top of each no matter the value of $L$. Fig. 9(b) shows the
effect of increasing the width of the domain on the time series of the
averaged heat transport. We find that the bursting regime disappears at large
aspect ratio and the system evolves toward a state in which the whole fluid
cools at a constant rate with $\langle F_{H}\rangle/F_{0}\approx 0.5$. We
perform additional simulations for $Pr=0.1$, 1, and 7 and find the same
behaviour. During the whole evolution of the simulations, the growing
convective layer has an aspect ratio that is always much larger than two and
thereby zonal flows are not expected to arise.

The suppression of the zonal flow has important consequences on the evolution
of the flow. First, the vertical transport is always significant and never
through quasi-periodic bursts. Second, the spatial structure of the flow is
different. Whereas the square domain simulations have flows dominated by
sheared convective plumes and bursts, the wider domain simulations exhibit
convective cells which persist in time. Third, the convective layer never
stops growing and reaches the bottom of the box mixing the whole fluid.

##  IV Summary and Discussion

We studied the onset and evolution of zonal flows when a convective layer
propagates into a fluid with a stable composition gradient. We considered
different values of the Prandtl number, $Pr=0.1$, 0.5, 1, and 7. Our goal was
to provide a novel way to study zonal flows and shear effects since the
growing convective layer allows exploration of a wide range of values of the
Rayleigh number and aspect ratios. Our results confirm and extend to
convection with stable composition gradients at low $Pr$ previous findings in
experiments of thermal and fingering convection.

In summary:

  1. 1.

In square domain simulations ($L/H=1$), zonal flows and their shear effects
always arise as long as the aspect ratio of the convective layer is smaller
than approximately two. The critical aspect ratio for the onset zonal flows
seems to depend on $Pr$, being smaller at low $Pr$ (Fig 1a).

  2. 2.

Once the zonal flow sets in, convective plumes get dispersed by the flow and
the convective layer stops growing.  The stalling of the convection zone
results in a saturation of the Rayleigh number toward a constant value in
time, $Ra\sim 10^{8}-5\times 10^{9}$, where the smallest and largest values
correspond to the cases using $Pr=7$, and 0.1, respectively (Fig. 1b).

  3. 3.

As found in previous numerical simulations of pure thermal convection, the
morphology and evolution of the flow depends on the Prandtl number. On the one
hand, for $Pr\leq 1$ the flow organizes into discrete bursts in which
convective plumes suddenly overturn quasi-periodically, with smaller transport
between bursts (Figs. 5a-c, 6, and 7a). On the other hand, for $Pr=7$ the flow
consists of sheared convective plumes instead of bursts, and the convective
transport is sustained at all times (Figs. 5d and 7c).

  4. 4.

The bursting regime is stronger at $Pr=1$ and weakens for $Pr<1$ (Fig. 7).

  5. 5.

We observe the formation of secondary convective layers in all the simulations
considered in this work (Fig. 8). The new layers contribute to the averaged
heat transport at all times for both the non-bursting and bursting regimes
(Fig 7).

  6. 6.

For wider domains ($L/H\geq 2$), the aspect ratio of the convective layer is
always much larger than two and zonal flows never develop during the evolution
of the simulations (Figs. 9a-b). The absence of large scale horizontal flow
means that the the growth of the convective layer is uninterrupted, and the
fluid fully mixes in all the cases considered.

We have shown that zonal flows also arise in the time-dependent problem of a
convective layer propagating into a stable fluid. Our work differs from
previous studies by considering a stable composition gradient and the fact
that convection is driven by a constant heat flux at the top boundary rather
than a constant temperature contrast across the fluid depth.

We find the onset of the zonal flow when the aspect ratio of the convection
zone is smaller than two and $Ra\sim 10^{8}$ \- $10^{10}$. The values of $Ra$
are narrower than in previous work. For example, Goluskin  _et al._ [6] found
zonal flows can arise when $Ra\sim 10^{4}$ \- $10^{10}$. A possible
explanation could be the fact that in our problem convection occurs just in a
portion of the fluid domain and the bottom of the convection zone is not
stress-free, meaning that horizontal fluid motions at the bottom of the layer
are decreased due to the interaction with the motion-less (stable) fluid
below. On the contrary, in Goluskin  _et al._ [6] convection occurs in the
whole box, and the stress-free boundaries enhance at all times the horizontal
fluid motions, which is favourable for instabilities that give rise to the
zonal flow.

Despite the differences mentioned above, the zonal flow and its effects are
similar to those reported in previous work. In particular, the transition from
the bursting to the non-bursting regime seems to occur at $Pr\simeq 1$ no
matter how convection is driven in the system of whether the fluid has
composition gradients. We did not explore in detail the range of $Pr$ in which
the system reaches the bursting regime. However, we find the intensity and
frequency of the bursts to decrease for $Pr<1$.

Recent work by Wang  _et al._ [19] explored in detail the influence of the
aspect ratio of the domain on the evolution of the zonal flow in thermal
convection. By imposing initial conditions consistent with a linear shear-
flow, and pure convective rolls, Wang  _et al._ [19] found that the zonal flow
only persists or arises when the aspect ratio of the domain is smaller than a
certain value depending on $Ra$ and $Pr$. For larger values, simulations
initialized with convective cells do not develop sheared-flows, and the ones
initialized with zonal flows experiment a transition to convective cells.
Those results support previous findings by Fitzgerald and Farrell [5], who
demonstrated that in fluid domains of large aspect ratio, the tilting
instability that enhances zonal flows can be suppressed due to the non-linear
interaction of horizontal modes of the velocity field. We did not perform a
spectral analysis to reveal the physical mechanism behind the zonal flow, but
our simulations including composition gradients seem to exhibit the same
behaviour. We find that zonal flows only appear when the aspect ratio of the
convective zone is smaller than two, and that it can be suppressed using wider
fluid domains (see Figs. 1 and 9).

The zonal flow and its effects are a problem for two-dimensional studies of
convective transport and mixing in fluids where strong shear flows are not
expected to appear. For example, layer formation and transport across
diffusive interfaces in double-diffusive convection [7, 19]. However, the fact
that zonal flows do not appear at large aspect ratio suggests that two-
dimensional simulations could still be useful to study convection and related
problems involving stable composition gradients (such as convective
overshoot). We did not explore different boundary conditions, however,
additional simulations with the numerical set-up of Garaud and Brummell [7]
but using a much larger aspect ratio would be of great interest in order to
see if zonal flows can also be avoided in that situation.

###### Acknowledgements.

This work was supported by an NSERC Discovery Grant. J.R.F. acknowledges
support from a McGill Space Institute (MSI) Fellowship. A. C. and J. R. F. are
members of the Centre de Recherche en Astrophysique du Québec (CRAQ) and the
Institut de recherche sur les exoplanètes (iREx). This research was enabled in
part by support provided by Calcul Québec (calculquebec.ca), and Compute
Canada (www.computecanada.ca). Computations were performed on Graham and
Béluga.

## References

  * Evonuk and Glatzmaier [2004] M. Evonuk and G. A. Glatzmaier, 2D Studies of Various Approximations Used for Modeling Convection in Giant Planets, Geophysical and Astrophysical Fluid Dynamics 98, 241 (2004).
  * Heimpel and Aurnou [2007] M. Heimpel and J. Aurnou, Turbulent convection in rapidly rotating spherical shells: A model for equatorial and high latitude jets on Jupiter and Saturn, Icarus 187, 540 (2007).
  * Finn [1993] J. M. Finn, Nonlinear interaction of Rayleigh-Taylor and shear instabilities, Physics of Fluids B 5, 415 (1993).
  * van der Poel  _et al._ [2014] E. P. van der Poel, R. Ostilla-Mónico, R. Verzicco, and D. Lohse, Effect of velocity boundary conditions on the heat transfer and flow topology in two-dimensional rayleigh-bénard convection, Phys. Rev. E 90, 013017 (2014).
  * Fitzgerald and Farrell [2014] J. G. Fitzgerald and B. F. Farrell, Mechanisms of mean flow formation and suppression in two-dimensional Rayleigh-Bénard convection, Phys. Fluids 26, 054104 (2014).
  * Goluskin  _et al._ [2014] D. Goluskin, H. Johnston, G. R. Flierl, and E. A. Spiegel, Convectively driven shear and decreased heat flux, J. Fluid Mech. 759, 360–385 (2014).
  * Garaud and Brummell [2015] P. Garaud and N. Brummell, 2D or Not 2D: The Effect of Dimensionality on the Dynamics of Fingering Convection at Low Prandtl Number, Astrophys. J. 815, 42 (2015).
  * Xie  _et al._ [2019] J.-H. Xie, K. Julien, and E. Knobloch, Jet formation in salt-finger convection: a modified Rayleigh-Bénard problem, J. Fluid Mech. 858, 228 (2019).
  * Garaud  _et al._ [2019] P. Garaud, A. Kumar, and J. Sridhar, The Interaction between Shear and Fingering (Thermohaline) Convection, Astrophys. J. 879, 60 (2019).
  * Thompson [1970] R. Thompson, Venus’s General Circulation is a Merry-Go-Round., Journal of Atmospheric Sciences 27, 1107 (1970).
  * Bian and Garcia [2003] N. H. Bian and O. E. Garcia, Confinement and dynamical regulation in two-dimensional convective turbulence, Physics of Plasmas 10, 4696 (2003).
  * Howard and Krishnamurti [1986] L. N. Howard and R. Krishnamurti, Large-scale flow in turbulent convection: a mathematical model, J. Fluid Mech. 170, 385 (1986).
  * Malkus [1954] W. V. R. Malkus, Discrete Transitions in Turbulent Convection, Proceedings of the Royal Society of London Series A 225, 185 (1954).
  * Krishnamurti and Howard [1981] R. Krishnamurti and L. N. Howard, Large-scale flow generation in turbulent convection, Proceedings of the National Academy of Sciences 78, 1981 (1981).
  * Hartlep  _et al._ [2003] T. Hartlep, A. Tilgner, and F. H. Busse, Large scale structures in rayleigh-bénard convection at high rayleigh numbers, Phys. Rev. Lett. 91, 064501 (2003).
  * Anders and Brown [2017] E. H. Anders and B. P. Brown, Convective heat transport in stratified atmospheres at low and high Mach number, Phys. Rev. Fluids 2, 083501 (2017).
  * Massaguer  _et al._ [1992] J. M. Massaguer, E. A. Spiegel, and J. P. Zahn, Convection-induced shears for general planforms, Phys. Fluids 4, 1333 (1992).
  * Matthews  _et al._ [1996] P. C. Matthews, A. M. Rucklidge, N. O. Weiss, and M. R. E. Proctor, The three-dimensional development of the shearing instability of convection, Phys. Fluids 8, 1350 (1996).
  * Wang  _et al._ [2020] Q. Wang, K. L. Chong, R. J. A. M. Stevens, R. Verzicco, and D. Lohse, From zonal flow to convection rolls in Rayleigh-Bénard convection with free-slip plates, arXiv e-prints , arXiv:2005.02084 (2020), arXiv:2005.02084 [physics.flu-dyn] .
  * Fuentes and Cumming [2020] J. R. Fuentes and A. Cumming, Penetration of a cooling convective layer into a stably-stratified composition gradient: entrainment at low Prandtl number, arXiv e-prints (2020).
  * Spiegel and Veronis [1960] E. A. Spiegel and G. Veronis, On the Boussinesq Approximation for a Compressible Fluid., Astrophys. J. 131, 442 (1960).
  * Burns  _et al._ [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Research 2, 023068 (2020).
  * Molemaker and Dijkstra [1997] M. J. Molemaker and H. A. Dijkstra, The formation and evolution of a diffusive interface, J. Fluid Mech 331, 199–229 (1997).
  * Radko [2016] T. Radko, Thermohaline layering in dynamically and diffusively stable shear flows, J. Fluid Mech. 805, 147–170 (2016).
  * Garaud [2017] P. Garaud, The formation of diffusive staircases, J. Fluid Mech. 812, 1–4 (2017).