# Staircase solutions and stability in vertically confined salt-finger
convection

Chang Liu1<EMAIL_ADDRESS>Keith Julien2  Edgar Knobloch1  1Department
of Physics, University of California, Berkeley, CA 94720, USA 2Department of
Applied Mathematics, University of Colorado, Boulder, CO 80309, USA

###### Abstract

Bifurcation analysis of confined salt-finger convection using single-mode
equations obtained from a severely truncated Fourier expansion in the
horizontal is performed. Strongly nonlinear staircase-like solutions having,
respectively, one (S1), two (S2) and three (S3) regions of well-mixed salinity
in the vertical direction are computed using numerical continuation, and their
stability properties are determined. Near onset, the one-layer S1 solution is
stable and corresponds to maximum salinity transport among the three
solutions. The S2 and S3 solutions are unstable but exert an influence on the
statistics observed in direct numerical simulations (DNS) in larger two-
dimensional (2D) domains. Secondary bifurcations of S1 lead either to tilted-
finger (TF1) or to traveling wave (TW1) solutions, both accompanied by the
spontaneous generation of large-scale shear, a process favored for lower
density ratios and Prandtl numbers ($Pr$). These states at low $Pr$ are
associated, respectively, with two-layer and three-layer staircase-like
salinity profiles in the mean. States breaking reflection symmetry in the
midplane are also computed. In 2D and for low $Pr$ the DNS results favor
direction-reversing tilted fingers (RTF) resembling the pulsating wave state
observed in other systems. Two-layer and three-layer mean salinity profiles
corresponding to RTF and TW1 are observed in 2D DNS averaged over time. The
single-mode solutions close to the high wavenumber onset are in an excellent
agreement with 2D DNS in small horizontal domains and compare well with 3D
DNS.

###### keywords:

Salt-finger convection, single-mode equations, double-diffusive convection

##  1 Introduction

Oceanographic measurements have widely reported staircase-like structures with
regions of nearly constant density in the vertical direction separated by
interfaces with a sharp density gradient (Tait & Howe, 1968, 1971; Schmitt et
al., 1987, 2005; Padman & Dillon, 1989; Muench et al., 1990; Zodiatis &
Gasparini, 1996; Morell et al., 2006; Timmermans et al., 2008; Fer et al.,
2010; Spear & Thomson, 2012). Such structures are typical of tropical and
subtropical regions where warm salty water often overlies cold fresh water and
have been ascribed to the presence of salt-finger convection. Indeed,
staircases have been observed in the western tropical Atlantic (Schmitt et
al., 1987, 2005), Tyrrhenian Sea (Zodiatis & Gasparini, 1996), and the
Mediterranean outflow in the Northeast Atlantic (Tait & Howe, 1968, 1971),
where conditions for salt-finger convection prevail. These staircases are
typically characterized by a large coherence length in the horizontal, much
larger than the step height in the vertical. For example, the Caribbean-Sheets
and Layers Transect (C-SALT) field programs show that the well-mixed layers
were 5–30 m thick and laterally coherent over scales of 200–400 km (Schmitt et
al., 1987). The presence of staircases also enhances the tracer mixing rate.
For example, the North Atlantic Tracer Release Experiments (NATRE) (Schmitt et
al., 2005) reveal a mixing rate in the western tropical Atlantic five times
larger than that in the eastern subtropical Atlantic. Spontaneous formation of
staircases has also been observed in idealized laboratory experiments (Linden,
1978; Krishnamurti, 2003, 2009) and direct numerical simulations (DNS)
(Piacsek & Toomre, 1980; Radko, 2003, 2005; Stellmach et al., 2011; Yang et
al., 2020).

Understanding the origin of spontaneous staircase formation is important for
parameterization and accurate modeling of oceanographic processes. Various
mechanisms for staircase formation have been proposed (Radko, 2013, §8). The
collective instability mechanism (Stern, 1969; Holyer, 1981) identifies
conditions for the formation of a staircase based on the instability of the
salt-finger field to internal gravity waves, but requires a closure model for
the Reynolds stress and the temperature and salinity fluxes based on
laboratory measurement and ocean observation. The onset of staircase formation
is also predicted by models with negative density diffusion (Phillips, 1972;
Posmentier, 1977). The related $\gamma$-instability of Radko (2003) requires a
parameterization of the Nusselt number and the flux ratio but shows good
agreement with 2D and 3D direct numerical simulation of salt-finger convection
(Radko, 2003; Stellmach et al., 2011). However, these predictions generally
only focus on the onset of a large-scale instability suggesting the appearance
of a staircase, but do not provide a detailed profile of the final staircase
or its parameter dependence. In this connection the model of Balmforth et al.
(1998) is of particular interest. The model parametrizes the coupling between
buoyancy flux and local turbulent kinetic energy but succeeds in generating
robust staircases. Nevertheless, the use of closure models employed in all
these predictions results in uncertainty in the applicable parameter regime.

This work aims to discuss an alternative mechanism for spontaneous staircase
formation from a bifurcation theory point of view, focusing on the computation
of strongly nonlinear staircase-like solutions and analyzing their stability.
Bifurcation analysis has been widely employed to provide insight into pattern
formation in fluid dynamics (Crawford & Knobloch, 1991). For example, a
secondary bifurcation of steady convection rolls in Rayleigh-Bénard convection
to tilted rolls was shown to be accompanied by the generation of large-scale
shear (Howard & Krishnamurti, 1986; Rucklidge & Matthews, 1996), resembling
both experimental observation (Krishnamurti & Howard, 1981) and direct
numerical simulations (Goluskin et al., 2014). A sequence of local and global
bifurcations of such tilted convection cells in magnetoconvection was shown to
lead to a pulsating wave characterized by periodic reversals in the direction
of the tilt and the accompanying large-scale shear (Matthews et al., 1993;
Proctor & Weiss, 1993; Proctor et al., 1994; Rucklidge & Matthews, 1996). In
the diffusive configuration in which cold fresh water overlies warm salty
water, a similar analysis found stable traveling waves near onset (Knobloch et
al., 1986) and provided insight into the transition to chaos (Knobloch et al.,
1992). Such chaotic or even fully developed turbulent states generally visit
neighborhoods of (unstable) steady, periodic, or traveling wave solutions, and
these visits leave an imprint on the flow statistics; see, e.g., Kawahara &
Kida (2001); van Veen et al. (2006) and the reviews by Kawahara et al. (2012);
Graham & Floryan (2021).

This work focuses on vertically confined salt-finger convection in order to
understand the interior between two well-mixed layers. While stress-free
velocity boundary conditions are suitable for understanding oceanographic
scenarios, no-slip boundary conditions are more relevant to laboratory
experiments (Hage & Tilgner, 2010). A wide range of bifurcation analyses of
related problems have been performed on vertically confined systems with
different boundary conditions including Rayleigh-Bénard convection. In the
salt-finger case the fluxes generated in a vertically confined system agree
quantitatively with those obtained in a vertically periodic domain when
normalized by the bulk conductive fluxes (Li & Yang, 2022) despite the
elimination of the elevator mode that is present in vertically periodic
domains; see, e.g., Stern (1969, §2.1); Holyer (1984, §2.1); Radko (2013,
§2.1).

In order to facilitate bifurcation analysis, we focus here on the single-mode
equations obtained from a severely truncated Fourier expansion in the
horizontal. Such single-mode equations reduce the two (or three) spatial
dimensions in the primitive equations into one vertical dimension, with the
dependence on the horizontal direction parameterized by a single assumed
horizontal wavenumber. Such a truncation may provide insight into salt-finger
convection as the underlying flow structures are dominated by well-organized
columnar structures; see, e.g., the visualization of confined salt-finger
convection by Yang et al. (2016a, figure 7). This single-mode formulation also
explicitly isolates and describes the interaction between horizontally
averaged modes (corresponding to the staircase or large-scale shear) and the
horizontal harmonics (corresponding to fingers), in the spirit of mean-field
theory; see, e.g, (Garaud, 2018, §3.2.1). Moreover, the single-mode equations
do not require any closure assumptions for the Reynolds stress and the
temperature and salinity fluxes to parametrize the feedback between the
fluctuations and horizontally averaged quantities, while preserving the
nonlinear interaction between them; see figure 1(b) below.

Single-mode equations have been widely employed to provide insight into
related problems. For example, single-mode solutions (also called
‘single-$\alpha$ mean-field theory’) for Rayleigh-Bénard convection with
stress-free boundary conditions have been shown to reproduce the mean
temperature profile expected at a high Rayleigh number (Herring, 1963). For
no-slip boundaries at high Prandtl numbers, the Nusselt number predicted from
this framework is within 20% of the experimental value, and the root-mean-
square velocity and temperature fluctuations also resemble the profiles seen
in experimental measurements (Herring, 1964). The single-mode theory has since
been extended to the time-dependent problem (Elder, 1969) and to more general
planforms (e.g., hexagonal) (Gough et al., 1975; Toomre et al., 1977), showing
qualitative agreement with experimental results. The vertical vorticity mode
was included within the single-mode equations with hexagon planform by (Lopez
& Murphy, 1983; Murphy & Lopez, 1984; Massaguer & Mercader, 1988; Massaguer et
al., 1990) and is excited beyond a secondary bifurcation of convection rolls.
Such single-mode equations have also proved to provide insight into many other
problems including rotating convection (Baker & Spiegel, 1975), plane
Poiseuille flow (Zahn et al., 1974) and double-diffusive convection (Gough &
Toomre, 1982; Paparella & Spiegel, 1999; Paparella et al., 2002). In
particular, DNS of single-mode solutions of salt-finger convection show that
fingers can tilt leading to the spontaneous formation of large-scale shear,
which may be steady or oscillate, together with a staircase-like profile in
the horizontally averaged salinity (Paparella & Spiegel, 1999; Paparella,
1997). Moreover, in certain asymptotic regimes single-mode solutions may be
exact; see, e.g., the high wavenumber (tall and thin) limit of Rayleigh-Bénard
convection (Blennerhassett & Bassom, 1994), convection in a porous medium
(Lewis et al., 1997) or salt-finger convection (Proctor & Holyer, 1986). Such
high wavenumber regimes are naturally achieved in convection with strong
restraints such as rapid rotation (Julien & Knobloch, 1997) or strong magnetic
field (Julien et al., 1999; Calkins et al., 2016; Plumley et al., 2018) and
the resulting exact single-mode equations have proved immensely useful for
understanding these systems; see the review by Julien & Knobloch (2007).

This work performs bifurcation analysis of vertically confined salt-finger
convection using single-mode equations. The resulting equations are solved for
the vertical structure of the solutions as a function of the density ratio,
the Prandtl number and the assumed horizontal wavenumber. We fix the
diffusivity ratio and thermal Rayleigh number and focus almost exclusively on
no-slip velocity boundary conditions. We found staircase-like solutions with
one (S1), two (S2) and three (S3) well-mixed mean salinity steps in the
vertical direction, all of which bifurcate from the trivial solution. In each
case salinity gradients are expelled from regions of closed streamlines,
leading to a well-mixed interior. Secondary bifurcations of S1 lead to either
tilted fingers (TF1) or traveling waves (TW1), both of which break the
horizontal reflection symmetry via the spontaneous formation of large-scale
shear. Secondary bifurcations of S2 and S3 lead to asymmetric solutions (A2
and A3) that spontaneously break midplane reflection symmetry.

The stability and accuracy of the single-mode solutions are further analyzed
with the assistance of 2D DNS. Near onset, the one-layer solution S1 is stable
and corresponds to maximum salinity transport among all known solutions, a
fact that is consistent with the prediction of the ‘relative stability’
criterion (Malkus & Veronis, 1958). The associated Sherwood number (or
salinity Nusselt number) near the high wavenumber onset is in excellent
agreement with DNS in small horizontal domains. In large domains the DNS
reveals a tendency to revert to the characteristic finger scale but the
wavenumber of the final state exhibits strong dependence on initial condition
(if the final state is steady) or the state may remain chaotic. The S2, S3,
A2, and A3 solutions are all unstable within the explored parameter regimes,
but may be imprinted on the statistics of the chaotic state. At lower density
ratios, the S1, S2 and S3 salinity profiles all sharpen and the large-scale
shear and tilt angle of the TF1 state increases. The TW1 solutions are also
present in this regime, while low Prandtl numbers favor direction-reversing
tilted fingers (RTF), a state in which the tilt direction and associated
large-scale shear reverse periodically as in the pulsating wave.

The remainder of this paper is organized as follows. §2 describes the problem
setup and the formulation of the single-mode equations. The equations are
solved using numerical continuation and their stability properties are
established for Prandtl number $Pr=7$ in §3. The tilted finger and traveling
wave states are analyzed for $Pr=0.05$ in §4. The paper concludes with a brief
summary and discussion in §5.

##  2 Single-mode equations for confined salt-finger convection

(a) (b)

Figure 1: (a) The salt-finger convection setup. The red color indicates a hot
plate while the blue color suggests a cold plate. The number of circles on top
and bottom plates suggests the salinity at the top is higher than at the
bottom. (b) Nonlinear interaction within the single-mode equations (7):
horizontal averages $\bar{(\cdot)}_{0}$ influence the horizontal harmonics
$\widehat{(\cdot)}$ through (7a)-(7d), while the horizontal harmonics
$\widehat{(\cdot)}$ contribute to the horizontal averages $\bar{(\cdot)}_{0}$
through (7e)-(7g).

We consider a fluid between two infinitely long parallel plates separated by a
distance $h$. The temperature and salinity at these two plates are maintained
at constant values with the top plate maintained at a higher temperature and
salinity, as illustrated in figure 1(a). The equation of state
$(\rho_{*}-\rho_{r*})/\rho_{r*}=-\alpha(T_{*}-T_{r*})+\beta(S_{*}-S_{r*})$ is
linear, with constant expansion/contraction coefficients $\alpha$, $\beta$ and
reference density, temperature, and salinity $\rho_{r*},T_{r*},S_{r*}$,
respectively. The subscript ∗ denotes a dimensional variable. In the following
we nondimensionalize the temperature $T_{*}$ by the temperature difference
between the top and bottom layer, $T=T_{*}/\Delta T$ ($\Delta T>0$), and
likewise for the salinity $S_{*}$, $S=S_{*}/\Delta S$ ($\Delta S>0$). Spatial
variables are nondimensionalized by the height $h$ of the layer while time and
velocity are nondimensionalized using the thermal diffusion time
$h^{2}/\kappa_{T}$ and the corresponding speed $\kappa_{T}/h$, respectively.
Here $\kappa_{T}$ is the thermal diffusivity. We decompose the temperature and
salinity into a linear base state and deviation,

$\displaystyle T=z+{\tilde{T}},\;\;S=z+{\tilde{S}},$ (1a-b)

and introduce the velocity field $\boldsymbol{u}:=(u,v,w)$ in Cartesian
coordinates $(x,y,z)$ with $z$ in the upward vertical direction. Dropping the
tildes and assuming the Boussinesq approximation the system is governed by

$\displaystyle\partial_{t}\boldsymbol{u}+\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u}=$
$\displaystyle
Pr\nabla^{2}\boldsymbol{u}-\boldsymbol{\nabla}p+PrRa_{T}\left(T-R_{\rho}^{-1}S\right)\boldsymbol{e}_{z},$
(2a) $\displaystyle\boldsymbol{\nabla}\cdot\boldsymbol{u}=$ $\displaystyle 0,$
(2b) $\displaystyle\partial_{t}T+\boldsymbol{u}\cdot\boldsymbol{\nabla}T+w=$
$\displaystyle\nabla^{2}T,$ (2c)
$\displaystyle\partial_{t}S+\boldsymbol{u}\cdot\boldsymbol{\nabla}S+w=$
$\displaystyle\tau\nabla^{2}S.$ (2d)

Here $p$ is the dimensionless pressure $p=p_{*}h^{2}/\kappa_{T}^{2}\rho_{r*}$
while $\boldsymbol{e}_{z}$ in (2a) is the unit vector in the vertical
direction associated with buoyancy. The governing parameters include the
Prandtl number, the diffusivity ratio, the density ratio, and the thermal
Rayleigh number defined by:

$\displaystyle
Pr:=\frac{\nu}{\kappa_{T}},\;\;\tau:=\frac{\kappa_{S}}{\kappa_{T}},\;\;R_{\rho}:=\frac{\alpha\Delta
T}{\beta\Delta S},\;\;Ra_{T}:=\frac{g\alpha\Delta Th^{3}}{\kappa_{T}\nu},$
(3a-d)

where $\nu$ is the viscosity and $\kappa_{S}$ is the salinity diffusivity.

We impose constant temperature and salinity as boundary conditions at the top
and the bottom plates:

$\displaystyle T(x,y,z=0,t)=T(x,y,z=1,t)=0,$ (4a) $\displaystyle
S(x,y,z=0,t)=S(x,y,z=1,t)=0,$ (4b)

while for the velocity we adopt no-slip boundary conditions:

$\displaystyle\boldsymbol{u}(x,y,z=0,t)=\boldsymbol{u}(x,y,z=1,t)=0.$ (5)

Periodic boundary conditions in the horizontal are imposed on all variables.

We next formulate single-mode equations following the procedure in Herring
(1963, 1964); Gough et al. (1975); Gough & Toomre (1982); Paparella & Spiegel
(1999); Paparella et al. (2002). The single-mode ansatz is

$\displaystyle S(x,y,z,t)=$
$\displaystyle\bar{S}_{0}(z,t)+\widehat{S}(z,t)\,e^{\text{i}(k_{x}x+k_{y}y)}+c.c.,$
(6a) $\displaystyle T(x,y,z,t)=$
$\displaystyle\bar{T}_{0}(z,t)+\widehat{T}(z,t)\,e^{\text{i}(k_{x}x+k_{y}y)}+c.c.,$
(6b) $\displaystyle\boldsymbol{u}(x,y,z,t)=$
$\displaystyle\bar{U}_{0}(z,t)\boldsymbol{e}_{x}+\widehat{\boldsymbol{u}}(z,t)\,e^{\text{i}(k_{x}x+k_{y}y)}+c.c.,$
(6c) $\displaystyle p(x,y,z,t)=$
$\displaystyle\bar{P}_{0}(z,t)+\widehat{p}(z,t)\,e^{\text{i}(k_{x}x+k_{y}y)}+c.c.,$
(6d)

where c.c. denotes the complex conjugate. Equation (6a) decomposes the
departure of the salinity from the linear profile into a horizontally averaged
quantity $\bar{S}_{0}(z,t)$ and a single harmonic in the horizontal direction
associated with the wavenumber pair $(k_{x},k_{y})$ and characterized by the
complex amplitude $\widehat{S}(z,t)$. The temperature is decomposed similarly.
To allow for mean flow in the horizontal we decompose the velocity into a
large-scale shear $\bar{U}_{0}(z,t)\boldsymbol{e}_{x}$ and a harmonic
associated with the same wavenumber pair $(k_{x},k_{y})$ as in equation (6c).
The horizontally averaged vertical velocity is zero based on the continuity
equation (2b) and the boundary conditions $w(z=0)=w(z=1)=0$. The horizontally
averaged quantities $\bar{(\cdot)}_{0}$ are real and the amplitudes of the
horizontal harmonics $\widehat{(\cdot)}$ are in general complex. Equations (6)
assume a horizontal planform in the form of square cells, an assumption that
also includes 2D rolls when $k_{y}=0$. Other planforms such as hexagons
generate additional self-interaction terms in the single-mode equations (Gough
et al., 1975) and we leave this extension to future work.

We now substitute (6) into the governing equations (2), dropping all harmonics
beyond the first, and balance separately the horizontally averaged components
and the harmonic components. We eliminate the horizontally averaged pressure
via
$-\partial_{z}\bar{P}_{0}+PrRa_{T}(\bar{T}_{0}-R_{\rho}^{-1}\bar{S}_{0})=0$
and eliminate the harmonic component of the pressure using the continuity
equation (2b). The resulting single-mode equations can then be expressed in
terms of the vertical velocity $w$ and vertical vorticity
$\zeta:=\partial_{y}u-\partial_{x}v$:

$\displaystyle\partial_{t}\widehat{\nabla}^{2}\widehat{w}+\text{i}k_{x}\bar{U}_{0}\widehat{\nabla}^{2}\widehat{w}-\text{i}k_{x}\bar{U}_{0}^{\prime\prime}\widehat{w}=$
$\displaystyle
Pr\widehat{\nabla}^{4}\widehat{w}+Pr\widehat{\nabla}_{\perp}^{2}Ra_{T}\left(\widehat{T}-R_{\rho}^{-1}\widehat{S}\right),$
(7a)
$\displaystyle\partial_{t}\widehat{\zeta}+\text{i}k_{x}\bar{U}_{0}\widehat{\zeta}+\text{i}k_{y}\bar{U}_{0}^{\prime}\widehat{w}=$
$\displaystyle Pr\widehat{\nabla}^{2}\widehat{\zeta},$ (7b)
$\displaystyle\partial_{t}\widehat{T}+\text{i}k_{x}\bar{U}_{0}\widehat{T}+\widehat{w}\partial_{z}\bar{T}_{0}+\widehat{w}=$
$\displaystyle\widehat{\nabla}^{2}\widehat{T},$ (7c)
$\displaystyle\partial_{t}\widehat{S}+\text{i}k_{x}\bar{U}_{0}\widehat{S}+\widehat{w}\partial_{z}\bar{S}_{0}+\widehat{w}=$
$\displaystyle\tau\widehat{\nabla}^{2}\widehat{S},$ (7d)
$\displaystyle\partial_{t}\bar{U}_{0}+\partial_{z}\left(\widehat{w}^{*}\widehat{u}+\widehat{w}\widehat{u}^{*}\right)=$
$\displaystyle Pr\partial_{z}^{2}\bar{U}_{0},$ (7e)
$\displaystyle\partial_{t}\bar{T}_{0}+\partial_{z}\left(\widehat{w}^{*}\widehat{T}+\widehat{w}\widehat{T}^{*}\right)=$
$\displaystyle\partial_{z}^{2}\bar{T}_{0},$ (7f)
$\displaystyle\partial_{t}\bar{S}_{0}+\partial_{z}\left(\widehat{w}^{*}\widehat{S}+\widehat{w}\widehat{S}^{*}\right)=$
$\displaystyle\tau\partial_{z}^{2}\bar{S}_{0},$ (7g)
$\displaystyle\widehat{u}=$
$\displaystyle\frac{\text{i}k_{x}\partial_{z}\widehat{w}}{k_{x}^{2}+k_{y}^{2}}-\frac{\text{i}k_{y}}{k_{x}^{2}+k_{y}^{2}}\widehat{\zeta},$
(7h)

where the superscript ∗ denotes a complex conjugate and
$\widehat{\nabla}^{2}:=\partial_{z}^{2}-k_{x}^{2}-k_{y}^{2}$,
$\widehat{\nabla}^{2}_{\perp}:=-k_{x}^{2}-k_{y}^{2}$,
$\widehat{\nabla}^{4}:=\partial_{z}^{4}-2(k_{x}^{2}+k_{y}^{2})\partial_{z}^{2}+(k_{x}^{2}+k_{y}^{2})^{2}$,
$\bar{U}_{0}^{\prime}:=\partial_{z}\bar{U}_{0}$ and
$\bar{U}_{0}^{\prime\prime}:=\partial_{z}^{2}\bar{U}_{0}$. The corresponding
boundary conditions for the salinity and temperature are:

$\displaystyle\widehat{S}(z=0,t)=\widehat{S}(z=1,t)=\widehat{T}(z=0,t)=\widehat{T}(z=1,t)$
(8a) $\displaystyle=\,$
$\displaystyle\bar{S}_{0}(z=0,t)=\bar{S}_{0}(z=1,t)=\bar{T}_{0}(z=0,t)=\bar{T}_{0}(z=1,t)$
(8b) $\displaystyle=\,$ $\displaystyle 0,$ (8c)

while the no-slip boundary conditions in (5) correspond to

$\displaystyle\widehat{w}(z=0,t)=\widehat{w}(z=1,t)=\partial_{z}\widehat{w}(z=0,t)=\partial_{z}\widehat{w}(z=1,t)$
(9a) $\displaystyle=\,$
$\displaystyle\widehat{\zeta}(z=0,t)=\widehat{\zeta}(z=1,t)=\bar{U}_{0}(z=0,t)=\bar{U}_{0}(z=1,t)$
(9b) $\displaystyle=\,$ $\displaystyle 0.$ (9c)

Symmetry name Formulation | Primitive equations in (2) | Single-mode in (7)
---|---|---
Midplane reflection |  $z\rightarrow 1-z$ $(w,T,S)\rightarrow-(w,T,S)$ |  $z\rightarrow 1-z$ $(\widehat{w},\widehat{T},\bar{T}_{0},\widehat{S},\bar{S}_{0})\rightarrow-(\widehat{w},\widehat{T},\bar{T}_{0},\widehat{S},\bar{S}_{0})$
Horizontal reflection |  $x\rightarrow-x$ $u\rightarrow-u$ |  $k_{x}\rightarrow-k_{x}$ $(\widehat{u},\bar{U}_{0},\widehat{\zeta})\rightarrow-(\widehat{u},\bar{U}_{0},\widehat{\zeta})$
Horizontal translation | $x\rightarrow x+\delta x$ | $\widehat{(\cdot)}\rightarrow\widehat{(\cdot)}e^{\text{i}k_{x}\delta x}$
Table 1: Symmetry properties of the primitive equations in (2) and the single-
mode equations in (7).

Equations (7a)-(7d) are the governing equation for the first harmonic in the
presence of the large-scale fields $\bar{U}_{0}$, $\bar{T}_{0}$ and
$\bar{S}_{0}$. In particular, equation (7a) is the Orr-Sommerfeld equation
modified by an additional buoyancy term, while (7b) is known as the Squire
equation for the vertical vorticity (Schmid & Henningson, 2012). Note that the
large-scale shear $\bar{U}_{0}$ is self-induced (i.e., a variable to be solved
for) instead of an imposed background shear. Thus the terms involving
$\bar{U}_{0}$ in equations (7a)-(7d) are also _nonlinear_. The remaining
equations (7e)-(7g) are the governing equations for the horizontally averaged
quantities, distorted by the harmonic fluctuations. The resulting nonlinear
interaction between horizontally averaged quantities and horizontal harmonics
is summarized schematically in figure 1(b). Equation (7h) is used to obtain
the horizontal velocity $\widehat{u}$ from the $\widehat{w}$ and
$\widehat{\zeta}$ required for the computation of the Reynolds stress in (7e).

Similar single-mode equations were previously used to study double-diffusive
convection focusing on the diffusive regime; see, e.g., Gough & Toomre (1982,
equations (3.7)-(3.11)), but did not include the coupling to the large large-
scale shear $\bar{U}_{0}$. A numerical simulation of the same single-mode
equation in a 2D configuration (i.e., with $k_{y}=0$) and stress-free boundary
conditions was employed to study double-diffusive convection in both the
fingering regime (Paparella & Spiegel, 1999) and the diffusive regime
(Paparella et al., 2002).

Here we focus on the bifurcation properties of the single-mode equations (7).
These reflect the symmetries of the primitive equations (2), including
midplane reflection, horizontal reflection, and horizontal translation, as
summarized in table 1. Note that we arbitrarily assume that the large-scale
shear $\bar{U}_{0}$ is generated in the $x$-direction within the single-mode
equations (7). This is appropriate in a 2D configuration, but in 3D the large-
scale shear can be oriented in principle in any horizontal direction, a
possibility that is also left for future work.

In the following we perform bifurcation analysis using pde2path (Uecker et
al., 2014; Uecker, 2021a) to obtain nonlinear solution profiles and analyze
their stability. The vertical direction is discretized using the Chebyshev
differentiation matrix (Weideman & Reddy, 2000) implemented following Uecker
(2021b). The number of grid points used, including the boundary, range from
$N_{z}=65$ for results at $R_{\rho}=40$ to $N_{z}=257$ for results at
$R_{\rho}=2$; our continuation in $R_{\rho}$ also uses $N_{z}=257$ grid
points. The tolerance of the maximal absolute value of the residue over each
vertical location ($L_{\infty}$ norm) is set as $10^{-6}$. The horizontal
translation symmetry within the single-mode equations requires an additional
phase condition to robustly fix the phase and obtain a unique solution, and is
implemented following Rademacher & Uecker (2017). The stability of each
solution is examined via the eigenvalues of the associated Jacobian matrix.
For validation we reproduced the results for single-mode equations for
Rayleigh-Bénard convection (Herring, 1963, 1964; Toomre et al., 1977), as well
as the high wavenumber asymptotic single-mode equations for Rayleigh-Bénard
convection (Bassom & Zhang, 1994, section 3) and porous medium convection
(Lewis et al., 1997, section 3). Selected solution profiles obtained from
pde2path were also validated against the results obtained from the nonlinear
boundary value problem (NLBVP) solver in Dedalus (Burns et al., 2020), and
selected eigenvalue computations from pde2path are also consistent with the
eigenvalue problem (EVP) solver in Dedalus (Burns et al., 2020). Finally,
selected predictions from the bifurcation diagram and stability properties
(figure 18(a)) were also validated against direct numerical simulations (IVP
solver) of the single-mode equations (7) using Dedalus (Burns et al., 2020).

We also performed 2D direct numerical simulations (DNS) of the primitive
equations (2) using Dedalus (Burns et al., 2020) to further analyze the
accuracy and stability of the solutions obtained from the bifurcation analysis
of the single-mode equations. We focus on 2D domains with periodic boundary
conditions in the horizontal and the no-slip boundary conditions (5) in the
vertical direction. We use a Chebyshev spectral method in the vertical
direction with $N_{z}=128$ grid points with a Fourier spectral method in the
horizontal direction with $N_{x}=128$ and a dealiasing scaling factor 3/2. To
check the accuracy we doubled the number of grid points in both $x$ and $z$
directions for selected results and confirmed that the resulting Sherwood
number $Sh$ in figure 16 (steady rolls resembling the S1 solution) does not
change up to eight decimal places. Time is advanced using a 3rd-order 4-stage
diagonal implicit Runge-Kutta (DIRK) scheme coupled with a 4-stage explicit
Runge-Kutta (ERK) scheme (RK443) (Ascher et al., 1997, section 2.8). We
mention that we do not expect the DNS results to coincide in all cases with
the single-mode results. This is because the former include, in principle, all
spatial harmonics of the fundamental wavenumber; it is precisely these that
are omitted from the single-mode theory. We are interested in identifying
parameter regimes in which the single-mode theory provides quantitatively
accurate results as well as regimes in which it fails. In all cases, the
horizontal domain is selected based on the wavenumber employed in constructing
the single-mode solutions, as described in detail in the next section.

##  3 Single-mode solutions and stability at $Pr=7$

Abbreviation | Description | Bifurcate from | Stability | $\bar{U}_{0}(z)$ | Color
---|---|---|---|---|---
S1 | Symmetric one-layer solutions | Trivial | S/U | $=0$ | Black ( )
S2 | Symmetric two-layer solutions | Trivial | U | $=0$ | Red ( )
S3 | Symmetric three-layer solutions | Trivial | U | $=0$ | Blue ( )
TF1 | Tilted fingers | S1 | S/U | $\neq 0$ | Magenta ( )
TW1 | Traveling waves | S1 | S/U | $\neq 0$ | Green ( )
A2 | Asymmetric two-layer solutions | S2 | U | $=0$ | Cyan ( )
A3 | Asymmetric three-layer solutions | S3 | U | $=0$ | Brown ( )
Table 2: Summary of solution features, stability, induced large-scale shear
$\bar{U}_{0}(z)$ and line color employed in the bifurcation diagrams (figures
2, 10, 13, 16 and 18). Symmetry is with respect to mid-plane reflection. S/U
indicates both stable and unstable solutions exist along this solution branch
and U indicates that only unstable solutions were found.

We start with the results for $Pr=7$, a Prandtl number value appropriate to
oceanographic applications, and fix $\tau=0.01$, corresponding to the
diffusivity ratio between salinity and temperature. The thermal Rayleigh
number is fixed at $Ra_{T}=10^{5}$ for a direct comparison with DNS results
(Yang et al., 2015, 2016a). Table 2 provides a summary of the solutions, all
of which are described in more detail below. All the solutions in the table
are steady-state solutions, including the traveling waves (TW1) which are
computed as steady states in a frame moving with the phase speed $c$ of the
wave. This speed solves a nonlinear eigenvalue problem. The bifurcation
diagrams presented in this work (figures 2, 10, 13, 16 and 18) all show the
time-averaged Sherwood number $Sh$ as a function of the horizontal domain size
$L_{x}$ or equivalently the fundamental wavenumber $k_{x}=2\pi/L_{x}$, where

$\displaystyle Sh:=\langle\partial_{z}S|_{z=0}\rangle_{h,t}+1$ (10)

with $\langle\cdot\rangle_{h,t}$ indicating average in both the horizontal
direction and over time. Without confusion, $\langle\cdot\rangle_{h,t}$ is
also referred to as a _mean_. In particular, for steady-state single-mode
solutions, we obtain

$\displaystyle
Sh=\langle\partial_{z}\bar{S}_{0}|_{z=0}\rangle_{t}+1=\partial_{z}\bar{S}_{0}|_{z=0}+1,$
(11)

where $\langle\cdot\rangle_{t}$ denotes time-averaging. We also use
$\langle\cdot\rangle_{h}$ to denote averaging in the horizontal. For time-
dependent solutions such as those obtained from DNS, these two types of
averaging yield distinct results.

In the bifurcation diagrams reported below, we use thick lines to indicate
stable solutions while thin lines represent unstable solutions. When the
maximum amplitude of the large-scale shear mode
$\underset{z}{\text{max}}\;|\bar{U}_{0}(z)|$ along a branch is of the order of
machine precision or smaller we report it as zero (see table 2). In all other
cases $\bar{U}_{0}(z)\neq 0$. When the large-scale shear is not generated,
$\bar{U}_{0}=0$, the harmonic quantities $\widehat{w}(z,t)$,
$\widetilde{u}(z,t):=-\text{i}\widehat{u}(z,t)$, $\widehat{S}(z,t)$ and
$\widehat{T}(z,t)$ are real for a suitably chosen phase. For solutions
associated with an induced large-scale shear $\bar{U}_{0}\neq 0$, we report
both real and imaginary parts of the harmonic profile denoted by
$\mathcal{R}e[\cdot]$ and $\mathcal{I}m[\cdot]$, respectively. We refer to the
base state in which all variables within the single-mode theory vanish as the
trivial state. This state corresponds to $Sh=1$. Except for the results with
stress-free velocity boundary conditions in subsection 3.2, all our results
are computed for no-slip boundary conditions.

###  3.1 Solution profiles and comparison with DNS at $R_{\rho}=40$

We start from the dynamics near onset when $R_{\rho}=40$, i.e., for
$R_{\rho}\tau\sim O(1)$. This is the parameter regime close to the onset of
salt-finger instability with $R_{\rho,{\rm crit}}=1/\tau$ in vertically
periodic domains (Radko, 2013, equation (2.4)) previously studied by Radko &
Stern (1999, 2000) and Xie et al. (2017, 2019). The main challenge of the
single-mode approach is the correct choice of the horizontal wavenumber
$k_{x}$, cf. Toomre et al. (1977). As a result, we start by analyzing the
bifurcation diagram as a function of $k_{x}$ and compare the results with DNS
in domains of length $L_{x}=2\pi/k_{x}$. This approach provides insight into
the effects of the domain size.

(a) 2D: $k_{y}=0$ (b) 3D: $k_{y}=k_{x}$

Figure 2: Bifurcation diagrams as a function of the fundamental wavenumber
$k_{x}$ from the single-mode equations (7) for (a) 2D: $k_{y}=0$, (b) 3D:
$k_{y}=k_{x}$, and the parameters $R_{\rho}=40$, $Pr=7$, $\tau=0.01$,
$Ra_{T}=10^{5}$. The black squares show the $Sh$ of the steady state reached
using 2D DNS in domains of size $L_{x}=2\pi/k_{x}$ (table 3); the black
pentagram corresponds to DNS with $L_{x}=2\pi$ displaying persistent chaotic
behavior.

Figure 2 shows the bifurcation diagram at $R_{\rho}=40$ with no-slip boundary
conditions. We consider both 2D results with $k_{y}=0$ in figure 2(a) and 3D
results with $k_{y}=k_{x}$ in figure 2(b) indicating the same aspect ratio in
two different horizontal directions ($x$ and $y$). Here, these two bifurcation
diagrams show almost the same behavior, and we can transform the 3D results in
figure 2(b) by defining an equivalent 2D horizontal wavenumber

$\displaystyle k_{x,2D}=\sqrt{k_{x}^{2}+k_{y}^{2}}$ (12)

to obtain the same 2D results as in figure 2(a). This suggests that in this
parameter regime, the 2D and 3D results are close to one another after a
rescaling of the horizontal domain and wavenumber. This can be readily seen
from the single-mode equation (7) where the horizontal wavenumbers always
appear as $k_{x}^{2}+k_{y}^{2}$ when a large-scale shear is not generated,
i.e., $\bar{U}_{0}=0$.

(a) (b) (c) (d)

(e) S1 (f) S2 (g) S3

(h) S1 (i) S2 (j) S3

Figure 3: Solution profiles from the single-mode equations (7) in the form of
S1, S2 and S3 when $k_{x}=8$, $k_{y}=0$, $R_{\rho}=40$, $Pr=7$, $\tau=0.01$,
$Ra_{T}=10^{5}$. The first row shows the profiles of (a) $z+\bar{S}_{0}$, (b)
$\widehat{S}$, (c) $\widetilde{u}$, and (d) $\widehat{w}$. The second row
shows the reconstructed total salinity using (6a) and (1) for (e) S1, (f) S2
and (g) S3 solutions. The third row shows isocontours of the streamfunction
for (h) S1, (i) S2 and (j) S3 solutions.

The solutions S1, S2 and S3 all bifurcate from the trivial solution. The S1,
S2 and S3 solution profiles with $k_{x}=8$, $k_{y}=0$ are shown in figure 3.
The horizontally averaged total salinity profiles $z+\bar{S}_{0}(z)$ in figure
3(a) show that these solutions are, respectively, associated with one, two and
three relatively well-mixed regions. These solution profiles thus resemble the
staircases observed in field measurements, e.g., (Schmitt et al., 1987).
Staircase-like solutions are also shown in some snapshots from direct
numerical simulations (Piacsek & Toomre, 1980, figure 2) and reproduced by
Zhang et al. (2018, figure 3). Recent DNS results show the coexistence of
multiple states with one, two, or three well-mixed regions using different
initial conditions (Yang et al., 2020, figure 2). In all three cases the well-
mixed regions correspond to large values of the salinity amplitude
$\widehat{S}$ as shown in figure 3(b), as well as a large vertical velocity
$\widehat{w}$ as shown in figure 3(d). In contrast, the horizontal velocity
peaks outside of the well-mixed region in each solution.

The second row of figure 3 reconstructs the total salinity using equations
(6a) and (1), while the third row of figure 3 shows the isocontours of the 2D
streamfunction computed as
$\psi(x,z)=\bar{\psi}_{0}(z)+\widehat{\psi}(z)e^{\text{i}k_{x}x}+c.c.$ with
$\bar{\psi}_{0}(z)=-\int_{0}^{z}\bar{U}_{0}(\xi)d\xi$ and
$\widehat{\psi}=\widehat{w}/(\text{i}k_{x})$. The isocontours are equispaced
between $\pm 0.9$ of the maximum value. The line ( ) is used for positive
(clockwise) streamlines while ( , blue) indicates negative (counterclockwise)
streamlines throughout this work. These results suggest that each well-mixed
region is associated with one downward and one upward moving plume. Since
salinity tends to be homogenized in regions of closed steady streamlines
(Rhines & Young, 1983) while salinity gradients are expelled from these
regions, the resulting mean salinity exhibits an overall staircase-like
profile.

Previous analysis of salt-finger convection in the spirit of single-mode
solutions (Radko & Stern, 2000; Proctor & Holyer, 1986) also found the S1 type
of solution, but S2 and S3 solutions were not found. This is likely because
these studies (Radko & Stern, 2000; Proctor & Holyer, 1986) focused on the
asymptotic behavior close to the onset of instability but did not go beyond
the first bifurcation point from the trivial solution.

TF1

(a) (b) (c) (d)

Figure 4: Solution profiles of the TF1 state from the single-mode equations
(7) with $k_{x}=4$, $k_{y}=0$, $R_{\rho}=40$, $Pr=7$, $\tau=0.01$,
$Ra_{T}=10^{5}$. Panels (a) and (b) show the profiles of $z+\bar{S}_{0}$ and
$\widehat{S}$, respectively, while (c) and (d) show the reconstructed total
salinity using (6a) and (1) and the isocontours of the streamfunction.

(a) (b) (c) (d)

(e)A2 (f)A3

Figure 5: Solution profiles of the A2 and A3 states from the single-mode
equations (7) with $k_{x}=8$, $k_{y}=0$, $R_{\rho}=40$, $Pr=7$, $\tau=0.01$,
$Ra_{T}=10^{5}$. The first row shows the profiles of (a) $z+\bar{S}_{0}$, (b)
$\widehat{S}$, (c) $\widetilde{u}$ and (d) $\widehat{w}$. The second row shows
the reconstructed total salinity using (6a) and (1) and the isocontours of the
streamfunction for (e) A2, (f) A3.

We also examined secondary bifurcations of the solutions S1, S2, and S3 in
figure 2(a), focusing on the bifurcation points closest to their high
wavenumber onset. The resulting secondary branches are of two types,
corresponding to tilted fingers (TF1) and asymmetric layer spacing (A2 and
A3). These states break the left-right reflection and the midplane reflection
symmetry, respectively, as revealed by the corresponding solution profiles
shown in figures 4 and 5. The mean salinity profile $z+\bar{S}_{0}$ for TF1 in
figure 4(a) is associated with two relatively well-mixed regions near
$z\in[0.2,0.4]$ and $z\in[0.6,0.8]$, while in the interior $z\in[0.4,0.6]$ the
profile is close to linear. The corresponding profile of $\widehat{S}$ in
figure 4(b) has both real and imaginary components with even and odd symmetry
with respect to the midplane, respectively, indicating that the harmonics are
no longer in phase in the vertical. The reconstructed total salinity shown in
figure 4(c) reveals that the finger is now tilted; see also the streamfunction
shown in figure 4(d). The tilted finger generates a nonzero large-scale shear
with $\underset{z}{\text{max}}\;\bar{U}_{0}(z)=3.16\times 10^{-4}$. The
profile of the large-scale shear is shown farther below (figures 15(a) and
19(b)) for parameter values for which it is much stronger. A similar tilted
finger state accompanied by large-scale shear was observed in earlier
simulations of the single-mode equations (Paparella & Spiegel, 1999) and
reported in observations from the North Atlantic Tracer Release Experiment
(St. Laurent & Schmitt, 1999, figure 3) as well as in laboratory experiments
on salt-finger convection; see, e.g., Taylor & Bucens (1989, figure 2) and
Krishnamurti (2009, figure 2). A related mechanism has also been widely
reported in experiments on Rayleigh-Bénard convection (Krishnamurti & Howard,
1981) as well as in direct numerical simulations of 2D Rayleigh-Bénard
convection (Goluskin et al., 2014); see the review by Siggia (1994).

The solution profiles in figures 5(a)-(d) indicate that the A2 and A3
solutions spontaneously break the midplane reflection symmetry; states
obtained via reflection in the midplane are therefore also solutions. In
particular, the profile $z+\bar{S}_{0}(z)$ no longer passes through
$(1/2,1/2)$. However, the reconstructed total salinity profile and the
isocontours of the streamfunction in figures 5(e)-(f) still resemble,
qualitatively, the S2 and S3 profiles that distinguish these two asymmetric
solutions. In particular, the A2 streamlines show two counter-rotating but
unequal rolls in the vertical while the A3 solution exhibits three counter-
rotating rolls, much as in the S2 and S3 states shown in figures 3(i)-(j).
Solutions that are asymmetric with respect to midplane reflection have been
seen in magnetoconvection with a depth-dependent magnetic diffusivity (Julien
et al., 2000, figure 17) but are a consequence of forced symmetry breaking.
Here, such asymmetric solutions originate from spontaneous symmetry breaking.

The stability of these solutions is indicated by thick (stable) and thin
(unstable) lines in figure 2. Near onset, i.e., in the large or small
wavenumber regime, S1 is stable and with the largest $Sh$ among all the
solutions shown in the figure, a finding that is broadly consistent with the
‘relative stability’ criterion of Malkus & Veronis (1958).

I.C. $L_{x}$ | $2\pi/18$ | $2\pi/16$ | $2\pi/14$ | $2\pi/12$ | $2\pi/10$ | $2\pi/8$ | $2\pi/6$ | $2\pi/4$ | $2\pi/2$ | $2\pi$ | $4\pi$
---|---|---|---|---|---|---|---|---|---|---|---
S1 | S1 | S1 | S1 | S1 | S1 | S1 (2) | S1 (2) | S1 (3) | S1 (7) | C | C
S2 | - | S1 | S1 | S1 | S1 | RTF | S1 (2) | S1 (3) | S1 (5) | - | -
S3 | - | - | S1 | S1 | S1 | S1 (2) | S1 (2) | S1 (3) | - | - | -
TF1 | - | - | - | - | - | - | S1 (2) | S1 (3) | S1 (6) | - | -
Table 3: The flow structures of 2D DNS simulations at $t=3000$ in domains of
size $L_{x}$ and initial condition (I.C.) constructed from S1, S2, S3 and TF1
solutions using the ansatz (6) with $k_{x}=2\pi/L_{x}$, $k_{y}=0$. RTF
indicates direction-reversing tilted fingers and C represents chaotic
behavior; ‘-’ indicates that a nonzero single-mode solution at
$k_{x}=2\pi/L_{x}$ is not present based on figure 2(a). The number
$n\in\mathbb{Z}$ inside a bracket indicates that the final horizontal
wavenumber reached by the solution is $k_{x}=2\pi n/L_{x}$, $n>1$, i.e. that
the evolution results in a changed wavenumber.

The dynamics of unstable solutions are typically not easy to isolate and
analyze without suitable initial conditions. Here, we use our single-mode
solutions as initial conditions for DNS to provide additional insight into
their stability. We set the horizontal domain size as $L_{x}\in[2\pi/18,4\pi]$
and then use the single-mode solution profile at $k_{x}=2\pi/L_{x}$,
$k_{y}=0$, to construct a 2D initial condition using the ansatz in (6). The
final state after $t=3000$ for different $L_{x}$ and initial conditions based
on S1, S2, S3, and TF1 solutions is summarized in table 3. For small
horizontal domains ($L_{x}\leq 2\pi/10$) the domain constrains the finger to
be tall and thin and the DNS results with S1, S2, and S3 initial conditions
all transition to a solution resembling the one-layer solution S1. This is
consistent with the stability observation in the bifurcation diagram in figure
2. The DNS also provides the final $Sh$ numbers and these values are plotted
using black squares in figure 2(a). The predicted $Sh$ from the single-mode
solution is close to the DNS results, especially close to the high wavenumber
onset of the S1 solution ($k_{x}=19.251$). The accuracy of the single-mode
equations close to the high wavenumber onset is a consequence of the strong
damping of modes with wavenumbers $k_{x}=2\pi n/L_{x}$ ($n\in\mathbb{Z}$ and
$n\geq 2$, e.g., $k_{x}=4\pi/L_{x}$, $k_{x}=6\pi/L_{x}$) that leaves only the
trivial solution at these wavenumbers (see figure 2)

(a) $\langle S\rangle_{h}(z,t)$ (b) $z+\langle S\rangle_{h,t}$

$z+S(x,z,t)$

(c) (d) (e)

Figure 6: 2D DNS showing (a) horizontally averaged salinity $\langle
S\rangle_{h}(z,t)$ and (b) time-averaged and horizontally averaged total
salinity $z+\langle S\rangle_{h,t}$. The second row shows isocontours of the
total salinity $z+S(x,z,t)$ at (c) $t=32$, (d) $t=35$, and (e) $t=38$. The
horizontal domain size is $L_{x}=2\pi/8$ with initial condition in the form of
a S2 single-mode solution with $k_{x}=2\pi/L_{x}=8$, $k_{y}=0$. The parameters
are $R_{\rho}=40$, $Pr=7$, $\tau=0.01$, $Ra_{T}=10^{5}$. See supplementary
video 1.

(a) $\langle S\rangle_{h}(z,t)$ (b) $z+S(x,z,t)$ at $z=0.1$

Figure 7: 2D DNS results showing (a) horizontally averaged salinity $\langle
S\rangle_{h}(z,t)$ and (b) total salinity $z+S(x,z,t)$ at $z=0.1$. The
horizontal domain size is based on $L_{x}=2\pi/8$ with S3 single-mode initial
condition with $k_{x}=2\pi/L_{x}=8$, $k_{y}=0$. The parameters are
$R_{\rho}=40$, $Pr=7$, $\tau=0.01$, $Ra_{T}=10^{5}$. See supplementary video
2.

$z+S(x,z,t)$

(a) I.C.: S1 (b) I.C.: S2 (c) I.C.: TF1

Figure 8: 2D DNS results for the total salinity $z+S(x,z,t)$ at $t=3000$ in a
horizontal domain of size $L_{x}=\pi$ initialized by (a) S1, (b) S2, and (c)
TF1 single-mode solutions with $k_{x}=2\pi/L_{x}=2$, $k_{y}=0$.

(a) $z+S(x,z,t)$ at $z=1/2$ (b)
$\langle|\mathcal{F}_{x}(S)|\rangle_{t}(z;k_{x})$ (c)
$\langle|\mathcal{F}_{x}(u)|\rangle_{t}(z;k_{x})$

(d) (e) $t=416$ (f) $t=437$ (g) $t=478$

Figure 9: 2D DNS results displaying (a) total salinity $z+S(x,z,t)$ at
$z=1/2$, (b) $\langle|\mathcal{F}_{x}(S)|\rangle_{t}(z;k_{x})$, (c)
$\langle|\mathcal{F}_{x}(u)|\rangle_{t}(z;k_{x})$, (d) $z+\langle
S\rangle_{h,t}(z)$, and snapshots of the horizontally averaged total salinity
$z+\langle S\rangle_{h}(z,t)$ at (e) $t=416$, (f) $t=437$ and (g) $t=478$. The
horizontal domain size is $L_{x}=2\pi$ initialized by a S1 solution with
$k_{x}=2\pi/L_{x}=1$, $k_{y}=0$. The governing parameters are $R_{\rho}=40$,
$Pr=7$, $\tau=0.01$, and $Ra_{T}=10^{5}$.

The final states in a domain of size $L_{x}=2\pi/8$ initialized by S2 and S3
show different flow structures. Figure 6 shows DNS results initialized with S2
and $L_{x}=2\pi/8$. The horizontally averaged salinity deviation $\langle
S\rangle_{h}(z,t)$ in figure 6(a) takes the requisite two-layer form for short
times but then transitions into a one-layer state with superposed oscillations
in time. After averaging over time and in the horizontal direction, the total
salinity profile $z+\langle S\rangle_{h,t}$ in figure 6(b) shows two well-
mixed regions, with a broad interior region that is close to a linear profile.
This behavior is similar to the TF1 single-mode solution as shown in figure
4(a). The mean salinity profile in figure 6(b) also resembles the mean
temperature profile in Rayleigh-Bénard convection with induced large-scale
shear (Goluskin et al., 2014, figure 5(b)), where the temperature profile
shows two well-mixed regions close to the boundary, and an interior that
approaches a linear profile with increasing Rayleigh number. The three
snapshots of the total salinity in figure 6(c)-(e) reveal that this behavior
is associated with a new state, a direction-reversing tilted finger (RTF).
Such behavior is also observed within simulations of the single-mode equations
(Paparella & Spiegel, 1999) (termed a layering instability there) and is
accompanied by an oscillating large-scale shear. The direction-reversing
tilted finger state also resembles (at least phenomenologically) the
‘pulsating wave’ found by Proctor et al. (1994, figure 7) and Matthews et al.
(1993, figure 5) in magnetoconvection. A comparison between convection with or
without a stabilizing magnetic field also highlights the stabilizing restoring
force that is important in sustaining this ‘pulsating wave’ state (Rucklidge &
Matthews, 1996), here the stabilizing temperature gradient. We have not
explored the origin of this state but note that similar ‘pulsating waves’ in
other systems originate from a sequence of local and global bifurcations of a
tilted convection roll (Matthews et al., 1993; Proctor et al., 1994; Rucklidge
& Matthews, 1996).

Figure 7 shows the DNS results with an S3 initial condition and $L_{x}=2\pi/8$
and displays the dynamics that are typically involved in transitioning to a
higher wavenumber S1 solution. The profile of $\langle S\rangle_{h}$ in figure
7(a) initially shows a three-layer structure but starts to oscillate prior to
a transition to the one-layer state S1 at $t\approx 30$. After $t\approx 160$,
this one-layer solution transitions to another one-layer solution but with a
smaller magnitude in $\langle S\rangle_{h}$. The total salinity near the lower
boundary, $z=0.1$, is shown in figure 7(b) and further explains this
transition scenario. The one-layer solution starts to tilt at $t\approx 110$
and begins to show oscillations at $t\approx 130$. Beyond $t\approx 160$, the
horizontal wavenumber doubles leading to the S1 (2) state in table 3. Such a
transition to a higher wavenumber S1 solution is also observed with horizontal
domains $L_{x}\in[2\pi/6,2\pi/2]$, albeit with different possibilities for the
final wavenumber. Note that such transitions between different horizontal
wavenumbers are not possible within the single-mode approach.

For $L_{x}=\pi$, the final state shows different horizontal wavenumbers
depending on the three different initial conditions S1, S2, and TF1 as
indicated in table 3. The total salinity at $t=3000$ in these three final
states is shown in figure 8. Although the horizontal domain size and the
governing parameters are the same, different initial conditions lead to
different final horizontal wavenumbers and the DNS results suggest that these
persist up to at least $t=3000$. This observation suggests that these three
different initial conditions S1, S2, and TF1 lie in the basin of attraction of
S1 solutions associated with different and larger horizontal wavenumbers. This
is of course a consequence of the preferred horizontal scale of the fingers.
Rolls with different horizontal wavenumber obtained from different initial
conditions are a familiar phenomenon in Rayleigh-Bénard convection (Wang et
al., 2020) and in spanwise rotating plane Couette flow (Xia et al., 2018; Yang
& Xia, 2021).

Solutions that do not reach a well-organized structure will be referred to as
chaotic. Such behavior is observed in large domains, e.g. $L_{x}=2\pi$ and
$L_{x}=4\pi$. The total salinity in the midplane for $L_{x}=2\pi$ is shown in
figure 9(a), and displays chaotic behavior with multiple downward and upward
plumes in the horizontal. In order to further characterize the horizontal
length scale, we computed $\langle|\mathcal{F}_{x}(S)|\rangle_{t}(z;k_{x})$
and $\langle|\mathcal{F}_{x}(u)|\rangle_{t}(z;k_{x})$, where
$\mathcal{F}_{x}(\cdot)$ is the Fourier transform in the $x$-direction and
$|\cdot|$ represents the modulus of the obtained complex Fourier coefficients.
The results are shown in figures 9(b)-(c). The $k_{x}=0$ component of
$\langle|\mathcal{F}_{x}(\cdot)|\rangle_{t}(z;k_{x})$ is related to the
amplitude of the horizontal mean $\bar{(\cdot)}_{0}$ while the $k_{x}\neq 0$
contributions to $\langle|\mathcal{F}_{x}(\cdot)|\rangle_{t}(z;k_{x})$ are
associated with the harmonic components $\widehat{(\cdot)}$ in the single-mode
ansatz (6). Figure 9(b) displays a peak at $k_{x}=0$ that corresponds to the
deviation of the mean salinity from its linear profile. A second peak around
$k_{x}\approx 9$ in figure 9(b) corresponds to a wavenumber that is close to
that providing the largest $Sh$ in the S1 solution shown in figure 2(a).
However, in the midplane, the peak wavenumber shifts to $k_{x}=13$. For the
horizontal velocity, the plot of $\langle|\mathcal{F}_{x}(u)|\rangle_{t}$ in
figure 9(c) also reveals a peak at $k_{x}=9$ but now at multiple locations
dominated by regions near the boundaries, with a lower amplitude near $z=1/3$
and $2/3$, the peak region of horizontal velocity $\widetilde{u}$ in the
three-layer solution S3 shown in figure 3(c). Furthermore,
$\langle|\mathcal{F}_{x}(u)|\rangle_{t}$ also shows a local peak at $k_{x}=1$
near the midplane, resembling the horizontal velocity $\widetilde{u}$ of the
two-layer solution S2 shown in figure 3(c). These observations suggest that
the unstable S2 and S3 solutions both play a role in the chaotic behavior and
that their properties may manifest themselves in the statistics of the chaotic
state, cf. Kawahara & Kida (2001) and the reviews by Kawahara et al. (2012)
and Graham & Floryan (2021).

In order to further quantify the role of unstable solutions, we plot in
figures 9(e)-(g) the horizontally averaged profiles at three different times
for comparison with the horizontally averaged time-averaged state in figure
9(d). Although the salinity profile in the latter state does not show a
significant staircase, the horizontally averaged salinity profiles at the
three times shown in figures 9(e)-(f) do in fact show a tendency towards three
and two steps, respectively, while the snapshot in figure 9(g) shows instead a
salinity profile that temporarily breaks the midplane reflection symmetry.
These results provide evidence that the chaotic state visits the neighborhood
of not only the unstable staircase-like solutions S2 and S3 but also of the
asymmetric solution. The time-averaged $Sh$ within $t\in[100,500]$ is
$Sh\approx 2.8693$, a value lower than the maximal $Sh$ for S1 but higher than
those for S2 and S3 at the suggested peak wavenumber $k_{x}=9$ as shown in
figure 2(a), and suggests that these states may play a role in determining the
time-averaged salinity transport in this chaotic state.

###  3.2 No-slip versus stress-free boundary conditions

(a) 2D: $k_{y}=0$ (b) 3D: $k_{y}=k_{x}$

Figure 10: Bifurcation diagram computed from the single-mode equations (7)
with the same parameter values as in figure 2 but stress-free velocity
boundary conditions at top and bottom.

(a) (b) (c) (d)

(e) S1 (f) S2 (g) S3

Figure 11: Solution profiles from the single-mode equations (7) in the same
parameter regime as figure 3 but with stress-free velocity boundary conditions
at top and bottom. The first row shows the profiles of (a) $z+\bar{S}_{0}$,
(b) $\widehat{S}$, (c) $\widetilde{u}$, and (d) $\widehat{w}$. The second row
shows the reconstructed total salinity using (6a) and (1) for (e) S1, (f) S2,
and (g) S3 solutions.

(a) (b) (c) (d)

Figure 12: Solution profiles on the secondary branches computed from the
single-mode equations (7) in the same parameter regime as figure 5 but with
stress-free velocity boundary conditions at the top and bottom. Panel (a)
shows the profiles of $z+\bar{S}_{0}$ of the A2 and A3 solutions, while panel
(b) displays the isocontours of the streamfunction of the A3 solution. Panels
(c) and (d) show $z+\bar{S}_{0}$ and the isocontours of the streamfunction for
a steady solution preserving midplane reflection symmetry generated in a
secondary bifurcation of the S3 solution.

Although no-slip boundary conditions are relevant for experiments (Hage &
Tilgner, 2010), stress-free boundary conditions are more appropriate for
oceanographic applications. Here we compare the bifurcation diagram in the no-
slip case with the corresponding diagram with stress-free boundary conditions.
For this purpose we replace the original no-slip boundary conditions in (5) by
the stress-free boundary conditions

$\displaystyle w(x,y,z=0,t)=w(x,y,z=1,t)=0,$ (13a)
$\displaystyle\partial_{z}u(x,y,z=0,t)=\partial_{z}u(x,y,z=1,t)=0,\;\;$ (13b)
$\displaystyle\partial_{z}v(x,y,z=0,t)=\partial_{z}v(x,y,z=1,t)=0.$ (13c)

Within the single-mode equations, the corresponding velocity boundary
conditions in (9) are thus replaced as:

$\displaystyle\widehat{w}(z=0,t)=\widehat{w}(z=1,t)=\partial_{z}^{2}\widehat{w}(z=0,t)=\partial_{z}^{2}\widehat{w}(z=1,t)$
(14a) $\displaystyle=\,$
$\displaystyle\partial_{z}\widehat{\zeta}(z=0,t)=\partial_{z}\widehat{\zeta}(z=1,t)=\partial_{z}\bar{U}_{0}(z=0,t)=\partial_{z}\bar{U}_{0}(z=1,t)$
(14b) $\displaystyle=\,$ $\displaystyle 0.$ (14c)

Note that the governing equations with stress-free boundary conditions admit
an additional Galilean symmetry: the large-scale shear can be shifted by an
arbitrary constant, $\bar{U}_{0}\rightarrow\bar{U}_{0}+C$. To obtain unique
solutions we therefore impose the additional constraint
$\int_{0}^{1}\bar{U}_{0}(z,t)dz=0$ following the procedure in Uecker (2021a,
§6.9).

Figure 10 shows the resulting bifurcation diagram for comparison with figure
2. The diagrams are very similar: all of the previously discussed solution
branches are still present. The main difference is that the stress-free
boundary conditions typically lead to a larger $Sh$ than the no-slip boundary
conditions at the same horizontal wavenumber. This behavior is also consistent
with the observation by Yang et al. (2016b), where 3D DNS results indicate
that the flow morphology is qualitatively similar but stress-free boundary
conditions display a larger $Sh$.

Figure 11 displays solution profiles and total salinity for the S1, S2 and S3
solutions at the same parameter values as in figure 3 but with stress-free
instead of no-slip boundary conditions. Here, the horizontally averaged total
salinity $z+\bar{S}_{0}$ still displays the one-layer, two-layer and three-
layer solution profiles in vertical. The horizontal structure of these states
is qualitatively similar to the no-slip results in figure 3. The main
difference is that the horizontal velocity $\widetilde{u}$ and the gradient of
vertical velocity $\partial_{z}\widehat{w}$ no longer vanish at the
boundaries.

Figure 12(a) displays the A2 and A3 solution profiles with stress-free
boundary conditions. The A2 solution exhibits a similar profile as in the no-
slip boundary condition case shown in figure 5(a), while the horizontally
averaged salinity in A3 differs from that found with the no-slip boundary
conditions. Figure 12(b) shows the isocontours of the streamfunction in A3,
and shows two large and one small recirculation cell in the vertical
direction, while figure 5(f) shows one large and two small cells in the
vertical. Here, the A3 solution in the stress-free boundary condition case
originates from a secondary bifurcation at $k_{x}=13.665$ on the S3 branch,
while there is another secondary bifurcation at $k_{x}=13.703$ (not shown in
figure 10) that is closer to the high wavenumber onset of the S3 state
$(k_{x}=15.573)$ and that preserves midplane reflection symmetry. The mean
salinity of the state shown in figure 12(c) resembles the mean temperature
profile in magnetoconvection (Julien et al., 2000, figure 7(b)), and the
isocontours of the streamfunction in figure 12(d) show one large and two small
closed streamlines that are similar to those in the A3 solution with no-slip
boundary conditions in figure 5(f). The fact that this additional secondary
bifurcation appears closer to the high wavenumber onset of S3 than the
bifurcation to the A3 branch likely contributes to the difference in the
solution profiles of the A3 state with stress-free and no-slip boundary
conditions.

###  3.3 Dependence on density ratio

Oceanographic conditions display a wide range of density ratios, and this
subsection therefore explores the dependence of our results on the density
ratio. In particular, the staircases observed in oceanographic observations,
laboratory measurements, and DNS are typically associated with the parameter
regime $R_{\rho}\sim O(1)$ (Schmitt et al., 1987; Krishnamurti, 2003, 2009;
Radko, 2003) that is far away from the onset of the salt-finger instability
$R_{\rho,{\rm crit}}=1/\tau$ with $\tau=0.01$.

2D: $k_{y}=0$ 3D: $k_{y}=k_{x}$

(a) (b)

(c) (d) (e) (f)

Figure 13: Bifurcation diagrams as a function of the wavenumber $k_{x}$ from
the single-mode equations (7) with (a) 2D: $k_{y}=0$ and (b) 3D:
$k_{y}=k_{x}$. The governing parameters are $R_{\rho}=2$, $Pr=7$, $\tau=0.01$,
and $Ra_{T}=10^{5}$. Panels (c)-(d) are a zoom of the 2D results in panel (a)
and panels (e)-(f) are a zoom of the 3D results in panel (b).

Figure 13 shows the bifurcation diagram with density ratio $R_{\rho}=2$ and
no-slip boundary conditions for (a) 2D results with $k_{y}=0$ and (b) 3D
results with $k_{y}=k_{x}$. Here, we focus on the S1, S2, S3, TF1 solution
branches, as well as a new traveling wave (TW1) solution branch that does not
exist when $R_{\rho}=40$. Here, the $Sh$ of the S1, S2, and S3 solutions in
the 2D and 3D configurations are still the same after converting the results
by rescaling the horizontal wavenumber according to (12). However, the
stability of the S1 solutions is now different and the secondary bifurcation
points of S1 that lead to TF1 and TW1 also differ. This can be seen in the
zoom of regions close to the appearance of TF1 and TW1 in figures 13(c)-(f).
Specifically, the bifurcation points leading to the TF1 and TW1 branches in 2D
are closer to the high wavenumber onset of the S1 solution than in 3D,
suggesting that the 2D configuration forms tilted fingers or traveling wave
solutions more readily. This is consistent with the suggestion that large-
scale shear is generated more easily in a 2D configuration than in 3D
(Goluskin et al., 2014).

(a) (b) (c) (d)

(e) S1 (f) S2 (g) S3

Figure 14: Solution profiles from the single-mode equations (7). The first row
shows the profiles of (a) $z+\bar{S}_{0}$, (b) $\widehat{S}$, (c)
$\widetilde{u}$ and (d) $\widehat{w}$. The second row shows the reconstructed
total salinity using (6a) and (1) for (e) S1, (f) S2 and (g) S3 solutions. The
parameter values are $k_{x}=17.5$, $k_{y}=0$, $R_{\rho}=2$, $Pr=7$,
$\tau=0.01$ and $Ra_{T}=10^{5}$.

The S1, S2 and S3 solution branches at $R_{\rho}=2$ in figure 13 exhibit much
larger $Sh$ than those at $R_{\rho}=40$ in figure 2. This is expected as a
lower density ratio indicates stronger destabilization by the salinity
gradient. Figure 14 shows profiles of the S1, S2 and S3 solutions at
$R_{\rho}=2$ and $k_{x}=17.5$, $k_{y}=0$. Here, the well-mixed region in the
horizontally averaged total salinity $z+\bar{S}_{0}(z)$ corresponds to almost
constant salinity compared with the solution profile at $R_{\rho}=40$ in
figure 3. Outside of the well-mixed region, the profile $z+\bar{S}_{0}(z)$ at
$R_{\rho}=2$ exhibits a sharper interface compared with that at $R_{\rho}=40$
in figure 3. This is consistent with the observation of well-defined
staircases in flow regimes with $R_{\rho}\sim O(1)$ in field measurements
(Schmitt et al., 1987) and laboratory experiments (Krishnamurti, 2003, 2009)
as well as in DNS (Radko, 2003). The vertical profiles of the salinity
harmonic $\widehat{S}$ in the S1, S2 and S3 states shown in figure 14(b) all
show a local peak at the region associated with sharp interfaces. Both the
horizontal velocity and the vertical velocity in figures 14(c)-(d) show a
larger magnitude compared with the $R_{\rho}=40$ results in figures 3(c)-(d).

(a) (b)TF1 (c)TW1

Figure 15: Solution profiles of (a) TF1 and (b) TW1 from the single-mode
equations (7) showing the associated large-scale shear $\bar{U}_{0}$ and the
isocontours of total salinity. The governing parameters are $k_{x}=17.5$,
$k_{y}=0$, $R_{\rho}=2$, $Pr=7$, $\tau=0.01$ and $Ra_{T}=10^{5}$.

We now focus on the solutions TF1 and TW1 originating from secondary
bifurcations of the S1 state. Figure 15(a) shows the associated large-scale
shear $\bar{U}_{0}(z)$ of TF1 and TW1. This shear is antisymmetric with
respect to the midplane for TF1 but symmetric for TW1. Figure 15(b) shows the
total salinity profile for a tilted finger with a larger tilt angle compared
with the tilt at $R_{\rho}=40$ shown in figure 4(e). Figure 14(c) shows the
total salinity of TW1 in the comoving frame, with structures that tilt in
opposite directions above and below the midplane. This alternating tilt
direction resembles the ‘wavy fingers’ tilted either left or right observed in
experiments (Krishnamurti, 2003, figure 7), although whether such ‘wavy
fingers’ travel depends on the boundary conditions in the horizontal.

The stability properties of these solutions are also indicated in figure 13
and are similar to those for high density ratios (figure 2). Here, the S1
solution is stable near the onset, but loses stability to tilted fingers. The
solutions S2 and S3 are still unstable within all of the current parameter
regime. The traveling wave appears to gain stability within a certain
wavenumber regime close to the low wavenumber onset, as indicated in the zoom
in figures 13(c) and 13(e). Since the single-mode solutions are expected to be
more accurate in the high wavenumber regime, cf. figure 2(a), we postpone a
study of the dynamics of TF1 and TW1 to the lower Prandtl number considered in
§4 where these states form at higher wavenumbers.

Figure 16:  (a) Bifurcation diagrams displaying $Sh$ as a function of the
density ratio $R_{\rho}$ from the single-mode equations (7) with wavenumbers
$k_{x}$ and $k_{y}$ given by the scaling law (15a). The other parameters are
$Pr=7$, $\tau=0.01$ and $Ra_{T}=10^{5}$. The black star $(*)$ is obtained from
3D DNS (Yang et al., 2015, table 1). The black square $(\square)$ is obtained
from 2D DNS results in a horizontal domain of size
$L_{x}=2\pi/k_{x,2D}=2\pi/(\sqrt{2}k_{x})$ with $k_{x}$ given by (15a). Panel
(b) plots the results in (a) in compensated form to exhibit the scaling
exponent.  Figure 17: Solution profiles from the single-mode equations (7)
displaying horizontally averaged total salinity (black solid lines) and
horizontally averaged total temperature (blue dashed lines). Solution profiles
are labeled by integers and correspond to $R_{\rho}=10,5,2,1,0.5,0.2$, and
$0.1$ from left to right. Other parameter values are $Pr=7$, $\tau=0.01$,
$Ra_{T}=10^{5}$, with $k_{x}=k_{y}$ given by the scaling law in (15a).

Figure 16 displays the dependence of the $Sh$ on the density ratio $R_{\rho}$
for $\tau=0.01$ for 3D results with $k_{y}=k_{x}$ based on the scaling law

$\displaystyle k_{x}=$ $\displaystyle
k_{y}=R_{\rho}^{0.125}\pi/(14.8211Ra_{S}^{-0.2428}),$ (15a) $\displaystyle
Ra_{S}:=$ $\displaystyle\frac{g\beta\Delta
Sh^{3}}{\kappa_{S}\nu}=\frac{Ra_{T}}{R_{\rho}\tau},$ (15b)

obtained through a least-squares fit of the finger width computed from 3D DNS
(Yang et al., 2016a, figure 13(b)). Here, we can see that the S1 solution
remains stable near onset, and then bifurcates to stable tilted fingers (TF1)
as $R_{\rho}$ decreases. The figure also shows the $Sh$ from 3D DNS data (Yang
et al., 2015, table 1) using black stars. We also perform 2D direct numerical
simulations in domains of size $L_{x}=2\pi/k_{x,2D}=2\pi/(\sqrt{2}k_{x})$ with
$k_{x}$ given by (15a). Note that our simulation results do not display a
smooth transition to the $Sh$ from the 3D results (Yang et al., 2015), which
is likely due to the intrinsic difference between 2D and 3D simulations and
possibly different domain sizes. However, the $Sh$ obtained from the S1
single-mode solutions agree well with the DNS results near the high density
ratio onset. At lower density ratios, the single-mode solution S1 and TF1 both
overpredict the DNS results. This may be because both are steady-state
solutions, while the DNS results exhibit time-dependent behavior which likely
reduces the time-averaged $Sh$ number.

We find that $Sh\sim R_{\rho}^{-0.30}$ for the single-mode solutions within
the range $R_{\rho}\in[0.1,2]$ indicating that $Sh\sim Ra_{S}^{0.30}$ as a
function of the salinity Rayleigh number defined in (15b) and fixed $Ra_{T}$
and $\tau$. This result agrees well the scaling law of Nusselt number $Nu\sim
Ra_{T}^{0.30}(\text{ln}Ra_{T})^{0.20}$ obtained from single-mode solutions for
Rayleigh-Bénard convection and wavenumber scaling as $k_{x}\sim Ra_{T}^{1/4}$
(Gough et al., 1975, p. 713). Figure 16(b) displays the compensated Sherwood
number $Sh\,R_{\rho}^{0.30}$, showing that the S1 solutions from both the
single-mode approach and DNS (Yang et al., 2015) follow the same $Sh\sim
R_{\rho}^{-0.3}$ scaling law. Note that the DNS results for salt-finger
convection in the asymptotic regime of high $Ra_{S}$ instead suggest $Sh\sim
Ra_{S}^{1/3}$ (Yang et al., 2015, figure 7(a)) while experimental results
suggest $Sh\sim Ra_{S}^{4/9}$ (Hage & Tilgner, 2010, figure 8). The parameter
regime considered here is within the regime $Ra_{S}\leq 10^{8}$, a value that
may be insufficient to reliably establish the asymptotic scaling; see e.g.,
(Yang et al., 2015, figure 7(a)) and (Yang et al., 2016b, figure 3(a)). The
TF1 solutions suggest the scaling law $Sh\sim R_{\rho}^{-0.27}$ or $Sh\sim
Ra_{S}^{0.27}$ for fixed $Ra_{T}$ and $\tau$. This is consistent with the
observation in Rayleigh-Bénard convection that the formation of the large-
scale shear decreases the scaling exponent $\eta$ of $Nu\sim Ra_{T}^{\eta}$;
see e.g., Goluskin et al. (2014, figure 4).

Figure 17 shows the corresponding mean salinity and temperature profiles at
the same parameter values as those used for the 3D DNS results in Yang et al.
(2015, figure 2(a)). The solution profiles from the single-mode equations
closely match the qualitative trend of the mean salinity and temperature
profiles as a function of the density ratio in the DNS results (Yang et al.,
2015, figure 2(a)). In particular, both the single-mode and DNS results show
that the mean salinity profile displays a well-mixed region in the layer
interior, and that both show a small overshoot (a thin stably stratified
region) near the top and bottom boundaries. The mean total temperature remains
close to a linear profile, however, and shows a visible deviation from a
linear profile only at $R_{\rho}=0.1$, a fact also consistent with DNS
observation (Yang et al., 2015, figure 2(a)).

##  4 Tilted salt fingers and traveling waves at $Pr=0.05$

The bifurcation diagram for a low density ratio, $R_{\rho}=2$, and $Pr=7$
displays a branch of steady tilted fingers (TF1) and a branch of traveling
waves (TW1), both of which bifurcate from the symmetric one-layer solution S1
at intermediate wavenumbers. Here, we study these states in the low Prandtl
number regime that pushes the bifurcations to TF1 and TW1 to higher
wavenumbers and hence closer to the high wavenumber onset, where we expect the
single-mode solutions to be accurate. Low Prandtl number salt-finger
convection is of interest in astrophysical applications, where molecular
weight gradients compete with thermal buoyancy but heat transport is dominated
by photon diffusion (Garaud, 2018). However, to keep the notation consistent
with the previous sections, we continue to use the symbol $S$ to represent the
higher molecular weight component and refer to it as salinity.

(a) 2D: $k_{y}=0$ (b) 3D: $k_{y}=k_{x}$

Figure 18: Bifurcation diagrams as a function of the wavenumber $k_{x}$ from
the single-mode equations (7) for (a) 2D: $k_{y}=0$ and (b) 3D: $k_{y}=k_{x}$.
Other parameters are $R_{\rho}=40$, $Pr=0.05$, $\tau=0.01$, and
$Ra_{T}=10^{5}$. The markers correspond to the $Sh$ associated with steady
states resembling S1 ($\square$), TF1 ( $\triangledown$) and TW1 ($\triangle$)
solutions from 2D DNS in domains of size $L_{x}=2\pi/k_{x}$ (table 4). The
black pentagrams corresponding to the $Sh$ at $L_{x}=2\pi/k_{x}$ ($k_{x}=1$,
2, 4) indicate chaotic behavior.

Figure 18 shows the bifurcation diagram for the parameter values used in
figure 2 but with $Pr=0.05$. Here, we focus on the solution branches S1, TF1
and TW1. Compared with the $Pr=7$ results in figure 2, we see that the $Sh$
associated with the S1 solution remains the same, but the bifurcation to TF1
($k_{x}=17.593$) now occurs much closer to the high wavenumber onset of the S1
solution ($k_{x}=19.251$). Moreover, the traveling wave branch TW1 that is
present here does not appear at $R_{\rho}=40$, $Pr=7$ in figure 2. Evidently
low Prandtl numbers favor spontaneous formation of large-scale shear, as found
in DNS by Radko (2010, figures 1-2) and Garaud & Brummell (2015) as well as in
a reduced model valid in the asymptotic limit of low $\tau$ and low $Pr$ (Xie
et al., 2019).

A comparison between 2D and 3D results at $Pr=0.05$ in figure 18 shows that
the bifurcation point of TF1 and TW1 is closer to the high wavenumber onset of
S1 in the 2D configuration. The difference between 2D and 3D at $Pr=0.05$ in
figure 18 is more evident than that at $Pr=7$ in figure 2, as is the case in
Rayleigh-Bénard convection (van der Poel et al., 2013). Recalling that both
TF1 and TW1 are associated with the formation of large-scale shear (see table
2, and figures 15(a) and 19(b)), this difference suggests that the 2D
configuration favors the formation of large-scale shear, cf. (Garaud &
Brummell, 2015).

(a) (b) (c)TF1 (d)TW1

(e)TF1 (f)TW1 (g)TF1 (h)TW1

Figure 19: Solution profiles from the single-mode equations (7) for the TF1 and TW1 solutions displaying (a) $z+\bar{S}_{0}$, (b) $\bar{U}_{0}$ and the horizontal velocity $\widetilde{u}$ for (c) TF1 and (d) TW1. The second row shows the vertical vorticity mode for (e) TF1 and (f) TW1 and the isocontours of the total salinity for (e) TF1 and (f) TW1 using (6a) and (1). Other parameter values are $k_{x}=k_{y}=6.873$, $R_{\rho}=40$, $Pr=0.05$, $\tau=0.01$ and $Ra_{T}=10^{5}$.  I.C. $L_{x}$ | $2\pi/18$ | $2\pi/16$ | $2\pi/14$ | $2\pi/12$ | $2\pi/10$ | $2\pi/8$ | $2\pi/6$ | $2\pi/4$ | $2\pi/2$ | $2\pi$ | $4\pi$
---|---|---|---|---|---|---|---|---|---|---|---
S1 | S1 | TF1 | RTF | RTF | RTF | C | RTF (2) | C | C | C | C
TF1 | - | TF1 | RTF | RTF | RTF | C | RTF (2) | C | C | - | -
TW1 | - | TF1 | RTF | RTF | RTF | C | RTF (2) | C | - | - | -
Table 4: The flow structures from 2D DNS at $t=3000$ in domains of different
sizes $L_{x}$ and initial conditions (I.C.) constructed from the single-mode
approximations to the S1, TF1 and TW1 states using the ansatz (6) with
$k_{x}=2\pi/L_{x}$, $k_{y}=0$. RTF indicates direction-reversing tilted
fingers and C represents chaotic behavior; ‘-’ indicates that a nonzero
single-mode solution at $k_{x}=2\pi/L_{x}$ is not present based on figure
18(a). The number $n\in\mathbb{Z}$ inside a bracket indicates the horizontal
wavenumber $k_{x}=2\pi n/L_{x}$ reached at $t=3000$ if different from the
initial wavenumber $n=1$.

Figure 19 shows the solution profiles in the TF1 and TW1 states. Here we
choose a 3D wavenumber $k_{y}=k_{x}=6.873$ close to the wavenumber that
maximizes the $Sh$ in the S1 solution. Similar to the earlier observation of
tilted fingers in figure 4, the mean total salinity profile $z+\bar{S}_{0}(z)$
of TF1 shows two relatively well-mixed regions with a linear profile in the
middle of the layer. In contrast, the TW1 shows instead a three-layer
structure in $z+\bar{S}_{0}(z)$. Panel (b) in figure 19 displays the
associated large-scale shear $\bar{U}_{0}(z)$, and exhibits a similar
structure to that observed at $R_{\rho}=2$, $Pr=7$ in figure 15(a). Figures
19(c)-(d) display the vertical profile of the first harmonic of the horizontal
velocity $\widetilde{u}$ in the TF1 and TW1 states showing local peaks close
to $z=1/3$, $1/2$, and $2/3$, in contrast to the corresponding S1 solution
that displays peaks near the boundaries. These properties are in turn
reflected in the staircase profile of the mean salinity shown in figure 19(a).

Furthermore, in 3D the TF1 and TW1 solutions are also associated with nonzero
vertical vorticity as shown in figures 19(e)-(f). Here, the generation of the
vertical vorticity originates from the large-scale horizontal shear
$\bar{U}_{0}$, which provides a source term for the vertical vorticity
equation (7b). Vertical vorticity is also generated in bifurcations of steady
convection rolls (resembling S1) in single-mode solutions of Rayleigh-Bénard
convection on a hexagonal lattice (Lopez & Murphy, 1983; Murphy & Lopez, 1984;
Massaguer & Mercader, 1988; Massaguer et al., 1990) with a source term
provided by self-interaction on this lattice (e.g. Massaguer et al. (1990,
equation (2.2b))). As the current single-mode formulation is limited to roll
or square planforms, this source term is not present. Instead, it is the
associated large-scale shear $\bar{U}_{0}$ that provides the source term for
vertical vorticity, a possibility not included in the previous work (Lopez &
Murphy, 1983; Murphy & Lopez, 1984; Massaguer & Mercader, 1988; Massaguer et
al., 1990).

The salinity profiles of TF1 and TW1 are shown in figures 19(g)-(h). These
panels show qualitative similarity with the TF1 and TW1 shown in figures
15(b)-(c), where the TF1 solution is tilted in one direction, while the TW1 is
tilted in opposite directions above and below the midplane. However, the TF1
and TW1 at $R_{\rho}=2$, $Pr=7$ in figure 15 show a relatively well-mixed
region in the interior, $z\approx 1/2$, compared with the TF1 and TW1 at the
higher density ratio $R_{\rho}=40$ but small Prandtl number, $Pr=0.05$, shown
in figure 19.

The stability of these solutions is indicated in figure 18 by thick (stable)
and thin (unstable) lines. In 2D TF1 is stable near the onset of this branch,
but becomes unstable at an intermediate wavenumber. The traveling wave is
instead unstable near the onset, but then gains stability over an interval of
intermediate wavenumbers before losing it again, cf. the TW1 branch at
$R_{\rho}=2$ and $Pr=7$ in figure 13. In 3D the TF1 branch is always stable,
as in the high Prandtl number regime, cf. figures 2 and 13, while the TW1
branch is always unstable, cf. figure 18(b).

We next compare these results with those from 2D DNS following our analysis of
the $Pr=7$ results in table 3. The domain length $L_{x}$ in the horizontal is
taken from $[2\pi/18,4\pi]$ and the initial conditions are, respectively,
constructed using the S1, TF1, and TW1 single-mode solutions with
$k_{x}=2\pi/L_{x}$ and $k_{y}=0$. In this $L_{x}$ interval single-mode theory
predicts that S1, TF1 or TW1 may be stable. In table 4 we record the state
reached in each case at $t=3000$. We find that the DNS returns S1 or TF1
depending on initial conditions but also supports direction-reversing tilted
fingers (RTF) and chaotic states (C). In particular, with horizontal domain
$L_{x}=2\pi/18$, the final state exhibits S1 behavior consistent with the
single-mode prediction and the $Sh$ from DNS overlaps with the prediction from
single-mode equations as shown in figure 18(a). For $L_{x}=2\pi/16$, the final
state in table 4 shows TF1, which is the only stable solution in single-mode
theory at $k_{x}=16$. The corresponding $Sh$ at $L_{x}=2\pi/16$ from DNS also
overlaps with the TF1 single-mode solution as shown in figure 18(a).

(a) $\langle u\rangle_{h}(z,t)$ (b) (c)

Figure 20: 2D DNS of a traveling wave showing (a) $\langle u\rangle_{h}(z,t)$,
(b) $z+\langle S\rangle_{h,t}(z)$ and (c) $\langle u\rangle_{h,t}(z)$ averaged
over $t\in[100,300]$ (black lines) for comparison with $z+\bar{S}_{0}$ and
$\bar{U}_{0}$ from the single-mode TW1 branch (red dashed lines). The
parameters are $R_{\rho}=40$, $Pr=0.05$, $\tau=0.01$, and $Ra_{T}=10^{5}$. The
initial condition is a TW1 solution with $L_{x}=2\pi/14$ for the DNS; the
single-mode solutions are computed for $k_{x}=2\pi/L_{x}=14$, $k_{y}=0$.

In order to demonstrate the fidelity of the single-mode TW1 solutions near the
high wavenumber onset, figure 20 compares one DNS run with $L_{x}=2\pi/14$ and
a TW1 single-mode initial condition. The solution takes the form of a
traveling wave up to $t\approx 400$, and consequently we compare the time
average of this state over $t\in[100,300]$ with the corresponding single-mode
results. Specifically, figure 20(a) shows the horizontally averaged horizontal
velocity $\langle u\rangle_{h}(z,t)$ corresponding to the large-scale shear
mode $\bar{U}_{0}(z,t)$ in the single-mode equations. Figures 20(b)-(c)
demonstrate an essentially perfect agreement between the $z+\langle
S\rangle_{h,t}$ and $\langle u\rangle_{h,t}$ profiles computed from the DNS in
this time interval and their counterparts $z+\bar{S}_{0}$ and $\bar{U}_{0}$
computed from the TW1 branch of the single-mode solutions at
$k_{x}=2\pi/L_{x}=14$ and the same parameters. The $Sh$ also agrees with the
single-mode prediction as indicated in figure 18(a). At $t\approx 400$, the
TW1 state undergoes an abrupt transition to an oscillatory state with an
antisymmetric mean shear profile suggesting that this state is a tilted finger
state whose tilt direction reverses periodically in time, i.e., it is
direction-reversing tilted finger (RTF). Similar behavior was also found at
$Pr=7$ with $L_{x}=2\pi/8$ initialized by an S2 solution; see table 3 and
figure 6. All RTF states in tables 3 and 4 were identified based on the
direction-reversing behavior of the associated large-scale shear $\langle
u\rangle_{h}(z,t)$.

(a) $z+S(x,z,t)$ at $z=1/2$ (b) $\langle u\rangle_{h}(z,t)$

(c) (d) (e) (f)

Figure 21: 2D DNS showing space-time plots of (a) $z+S(x,z,t)$ at $z=1/2$ and
(b) $\langle u\rangle_{h}(z,t)$. Panels (c) and (d) show the profiles of
$z+\langle S\rangle_{h,t}(z)$ averaged over $t\in[80,180]$ for the traveling
wave episode, and $t\in[280,380]$ for direction-reversing tilted fingers, and
compared with $z+\bar{S}_{0}$ from the single-mode TW1 and TF1 solutions.
Panel (e) shows $\langle u\rangle_{h,t}$ averaged over $t\in[80,180]$ and
compared with $\bar{U}_{0}$ for the single-mode TW1 branch, while (f) shows
$\langle u\rangle_{h}(z,t)$ at $t=332$, 333 and 334 for comparison with
$\bar{U}_{0}$ from the single-mode TF1 solution. The initial condition is a S1
solution with $L_{x}=2\pi/10$, and other parameters are $R_{\rho}=40$,
$Pr=0.05$, $\tau=0.01$, $Ra_{T}=10^{5}$. The single-mode solutions are
calculated for $k_{x}=2\pi/L_{x}=10$, $k_{y}=0$. See supplementary video 3.

(a) $z+S(x,z,t)$ at $z=1/2$ (b)
$\langle|\mathcal{F}_{x}(S)|\rangle_{t}(z;k_{x})$ (c)
$\langle|\mathcal{F}_{x}(u)|\rangle_{t}(z;k_{x})$

(d) $\langle u\rangle_{h}(z,t)$

Figure 22: 2D DNS results showing (a) total salinity at $z=1/2$, (b)
$\langle|\mathcal{F}_{x}(S)|\rangle_{t}(z;k_{x})$, (c)
$\langle|\mathcal{F}_{x}(u)|\rangle_{t}(z;k_{x})$ and (d) $\langle
u\rangle_{h}(z,t)$. The horizontal domain size is $L_{x}=2\pi$ initialized by
the corresponding S1 solution. The parameters are $R_{\rho}=40$, $Pr=0.05$,
$\tau=0.01$, and $Ra_{T}=10^{5}$.

Direction-reversing tilted fingers are also observed as the final state in
larger horizontal domains with $L_{x}=2\pi/12$, $2\pi/10$ and $2\pi/6$ as
shown in table 4. Figure 21 shows the DNS results with $L_{x}=2\pi/10$ and
initial condition S1. The total salinity $z+S$ at the midplane $z=1/2$ in
figure 21(a) shows that the initial state is sinusoidal in space and steady,
corresponding to the prescribed S1 solution. After $t\approx 80$, the S1 state
transitions into a traveling wave and subsequently, at $t\approx 220$, into an
RTF state. Figure 21(b) shows that $\langle u\rangle_{h}(z,t)=0$ at $t=0$ but
then this quantity becomes nonzero and midplane-symmetric at the transition to
the TW state with superposed periodic oscillations. Figures 21(c)-(d) present
the horizontally averaged total salinity profiles, respectively averaged over
$t\in[80,180]$ and $t\in[280,380]$, corresponding to the traveling wave and
the direction-reversing tilted fingers, and compared with $z+\bar{S}_{0}$ for
the single-mode TW1 and TF1 solutions with $k_{x}=2\pi/L_{x}$ and the same
parameter values. Here, both the DNS results and the single-mode solutions
show a modest three-layer staircase characteristic of the traveling wave
state, while both show a two-layer staircase with a linear profile in the
middle that is associated with direction-reversing tilted fingers in DNS and
the tilted finger state of the single-mode equations.

Figure 21(e) shows $\langle u\rangle_{h,t}$ from the DNS averaged over
$t\in[80,180]$ corresponding to the traveling wave, and compares it with
$\bar{U}_{0}$ from the corresponding single-mode TW1 solution. Here, the
large-scale shear from the single-mode theory overpredicts the amplitude
compared to that observed in DNS. This is expected as figure 21(b) indicates
that the state is in fact a modulated traveling wave with superposed
oscillations in $\langle u\rangle_{h}$ in $t\in[80,180]$. Moreover, as shown
previously in figures 2(a) and 18(a), the accuracy of the single-mode equation
is expected to decrease with decreasing wavenumber (here $k_{x}=10$), a fact
that may account for the difference in the solution profiles between DNS and
single-mode theory. This difference in the large-scale shear may also
contribute to the difference in the mean salinity profiles in figure 21(c) and
the associated $Sh$ in figure 18(a). As the TF1 solution in the single-mode
equation is steady, we compare the large-scale shear $\bar{U}_{0}$ with three
different snapshots of $\langle u\rangle_{h}(z,t)$ at $t=332$, 333, 334 in
figure 21(f). The shape and magnitude of $\langle u\rangle_{h}$ at $t=333$ is
close to a TF1 single-mode solution, but the maximum amplitude of the large-
scale shear of RTF observed in 2D DNS is larger than $\bar{U}_{0}$ from the
TF1 single-mode solution.

For larger $L_{x}$, the solution may transition to a double wavenumber RTF,
see $L_{x}=2\pi/6$ in table 4, and such a transition to a higher wavenumber
flow structure was already described for $Pr=7$ as summarized in table 3 and
figures 7 and 8. Since the final state does not exhibit well-organized flow
structures when $L_{x}\in[2\pi/4,4\pi]$ we refer to it as chaotic. Figure 22
presents (a) the total salinity at the midplane, (b) the Fourier-transformed
salinity deviation $\langle|\mathcal{F}_{x}(S)|\rangle_{t}$, and (c) the
Fourier-transformed horizontal velocity
$\langle|\mathcal{F}_{x}(u)|\rangle_{t}$. Here, we select $L_{x}=2\pi$ and the
initial condition is a S1 solution, so the only difference between the results
in figures 22 and 9 is the value of the Prandtl number. The total salinity at
the midplane $z=1/2$ in figure 22(a) at $Pr=0.05$ does not show a very clear
difference from figure 9(a) at $Pr=7$, while the Fourier spectrum of the
salinity at $P=0.05$ in figure 22(b) is broader than the same quantity at
$Pr=7$ in figure 9(b). The Fourier-transformed horizontal velocity in figure
22(c) shows a greater difference from the corresponding result at $Pr=7$ in
figure 9(c), where the peaks near $z=1/3$ and $z=2/3$ and the peak at $z=1/2$
at $k_{x}\approx 1$ have magnitudes similar to those near the boundary. This
suggests that both TF1 and TW1 are involved in the chaotic behavior. Moreover,
$\langle|\mathcal{F}_{x}(u)|\rangle_{t}$ shows a nonzero value at $k_{x}=0$
corresponding to the large-scale shear in figure 22(c), while the white region
in $\langle|\mathcal{F}_{x}(u)|\rangle_{t}$ at $k_{x}=0$ in figure 9(c) for
$Pr=7$ indicates a value close to zero. This is consistent with the fact that
low Prandtl numbers favor the formation of large-scale shear. Indeed, $\langle
u\rangle_{h}(z,t)$ in figure 22(d) reaches a maximum value of 0.065, while the
same quantity has the maximum amplitude $5.37\times 10^{-4}$ when $Pr=7$.

At other times, the large-scale shear $\langle u\rangle_{h}(z,t)$ in figure
22(d) resembles the large-scale shear $\bar{U}(z)$ associated with either TF1
or TW1 in figure 19(b) based on the parity of $\langle u\rangle_{h}(z,t)$ with
respect to $z=1/2$. The observed time series suggests that the solution visits
the neighborhoods of both the traveling waves and the tilted fingers in a
chaotic manner. The mean $Sh$ associated with this state in domains of size
$L_{x}=\pi/2$, $\pi$ and $2\pi$, and initialized using the S1 single-mode
solution are shown in figure 18(a), and fall within the range predicted by S1,
TF1, and TW1 single-mode states. Note that for these relatively large domain
sizes, transition to a higher wavenumber mode is also possible but cannot be
predicted within the single-mode approach. The comparison in figure 18(a)
supports the role of the (unstable) TF1 and TW1 solutions in determining the
salinity transport in this chaotic state.

##  5 Conclusions and future work

This work performed bifurcation analysis of vertically confined salt-finger
convection using single-mode equations obtained from a severely truncated
Fourier expansion in the horizontal. The resulting equations were solved for
the vertical structure of the solutions as a function of the density ratio,
the Prandtl number, and the assumed horizontal wavenumber. We fixed the
diffusivity ratio and thermal Rayleigh number and focused almost exclusively
on the case of no-slip velocity boundary conditions. We computed strongly
nonlinear staircase-like solutions having one (S1), two (S2) and three (S3)
well-mixed mean salinity regions in the vertical direction. These bifurcate in
successive bifurcations from the trivial solution. In each case, salinity
gradients are expelled from regions of closed streamlines resembling the
mechanism described by Rhines & Young (1983), leading to a well-mixed
interior. Secondary bifurcations of S1 lead to tilted fingers (TF1) or
traveling waves (TW1) both of which spontaneously break reflection symmetry in
the horizontal owing to the spontaneous generation of large-scale shear.
Secondary bifurcations from S2 and S3 lead to asymmetric solutions (A2 and A3)
that spontaneously break symmetry with respect to the horizontal midplane.

The stability and relevance of the single-mode solutions we obtained were
further analyzed with the assistance of 2D DNS. Near onset, the one-layer
solution S1 is stable and corresponds to maximum salinity transport among the
three solutions, an observation consistent with the prediction of ‘relative
stability’ criterion (Malkus & Veronis, 1958). The associated Sherwood number
near the high wavenumber onset is in excellent agreement with DNS in small
horizontal domains. In larger domains DNS reveals that the final state reverts
to a higher wavenumber S1 state closer to the natural finger scale or exhibits
chaotic behavior. In general the final state exhibits strong dependence on
initial conditions. The S2, S3, A2 and A3 solutions are all unstable within
the parameter regime explored, although the chaotic solutions we found appear
to visit neighborhoods of these unstable solutions at different instants in
time. With stress-free velocity boundary conditions the S1, S2, S3, TF1, A2
and A3 solutions persist whenever the density ratio and the Prandtl number are
high enough, but exhibit larger Sherwood numbers.

At lower density ratios, the S1, S2 and S3 solutions exhibit sharper staircase
structure, while the TF1 solution displays stronger large-scale shear and tilt
angle. In addition S1 also bifurcates to the TW1 state, a bifurcation that is
not present at high density ratio. The scaling of the S1 Sherwood number with
the density ratio closely matches the DNS results within the currently
explored parameter regime, while the corresponding change in the mean salinity
and temperature profiles resembles that seen in DNS.

The dynamics of the secondary TF1 and TW1 states in the low Prandtl number
regime also show excellent agreement with the 2D DNS results in small domains,
likely because they appear, in this regime, quite close to the high wavenumber
onset of S1. These states are typically associated with two- and three-layer
mean salinity profiles, respectively. The final state seen in low Prandtl
number 2D DNS is a prominent direction-reversing tilted finger (RTF) in which
the tilt and the associated large-scale shear reverse periodically with time
in a manner that resembles the ‘pulsating wave’ state identified in
magnetoconvection (Matthews et al., 1993; Proctor et al., 1994). In smaller
domains the time-averaged salinity and shear profiles computed from 2D DNS
resemble the two-layer and three-layer mean salinity and shear profiles
corresponding to RTF and TW1, respectively, while in larger domains such
profiles are still evident but only episodically.

The bifurcation diagrams of single-mode solutions shown here also shed light
on the difference between 2D and 3D. While the corresponding S1, S2, S3, A2
and A3 states in 3D can be obtained from the 2D results by a simple wavenumber
rescaling, this is no longer so for states associated with large-scale shear,
i.e. TF1 and TW1. The resulting difference is most prominent at low density
ratio or low Prandtl numbers. A comparison shows that 2D favors the generation
of large-scale shears, a fact consistent with DNS observations (Garaud &
Brummell, 2015).

The bifurcation diagram presented here is far from complete. Whether other
primary or secondary bifurcations possess the potential to provide stable
staircase solutions is an interesting question for future exploration. The
secondary Hopf and global bifurcations are also of interest since they provide
additional insight into the origin of the direction-reversing tilted finger
state and of chaotic salt-finger convection, respectively. Extension of the
present framework to 3D states with hexagonal coordination in the horizontal
is also of interest, as is a study of the doubly- or triply-periodic
configuration with periodic boundary conditions in the vertical as well as in
the horizontal, a formulation that is likely to have greater relevance to
oceanography. Finally, density staircases have also been widely observed in
the diffusive regime of double-diffusive convection (cold fresh water on top
of warm salty water) (Timmermans et al., 2008; Radko, 2016; Yang et al., 2022)
as well as in stratified shear flow (Oglethorpe et al., 2013; Taylor & Zhou,
2017; Lucas et al., 2017, 2019) and elsewhere. It is of interest to explore
the applicability of the present bifurcation-theoretic approach to the study
of staircase states in these flow regimes.

## Acknowledgment

This work was supported by the National Science Foundation under Grant Nos.
OCE 2023541 (C.L. and E.K.) and OCE-2023499 (K.J.).

## Declaration of Interests

The authors report no conflict of interest.

## References

  * Ascher et al. (1997) Ascher, U. M., Ruuth, S. J. & Spiteri, R. J. 1997 Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Appl. Numer. Math. 25, 151–167.
  * Baker & Spiegel (1975) Baker, L. & Spiegel, E. A. 1975 Modal analysis of convection in a rotating fluid. J. Atmos. Sci. 32, 1909–1920.
  * Balmforth et al. (1998) Balmforth, N. J., Llewellyn Smith, S. G. & Young, W. R. 1998 Dynamics of interfaces and layers in a stratified turbulent fluid. J. Fluid Mech. 355, 329–358.
  * Bassom & Zhang (1994) Bassom, A. P. & Zhang, K. 1994 Strongly nonlinear convection cells in a rapidly rotating fluid layer. Geophys. Astrophys. Fluid Dyn. 76, 223–238.
  * Blennerhassett & Bassom (1994) Blennerhassett, P. J. & Bassom, A. P. 1994 Nonlinear high-wavenumber Bénard convection. IMA J. Appl. Math. 52, 51–77.
  * Burns et al. (2020) Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D. & Brown, B. P. 2020 Dedalus: A flexible framework for numerical simulations with spectral methods. Phys. Rev. Res. 2, 023068.
  * Calkins et al. (2016) Calkins, M. A., Julien, K., Tobias, S. M., Aurnou, J. M. & Marti, P. 2016 Convection-driven kinematic dynamos at low Rossby and magnetic Prandtl numbers: single mode solutions. Phys. Rev. E 93, 023115.
  * Crawford & Knobloch (1991) Crawford, J. D. & Knobloch, E. 1991 Symmetry and symmetry-breaking bifurcations in fluid dynamics. Annu. Rev. Fluid Mech. 23, 341–387.
  * Elder (1969) Elder, J. W. 1969 The temporal development of a model of high Rayleigh number convection. J. Fluid Mech. 35, 417–437.
  * Fer et al. (2010) Fer, I., Nandi, P., Holbrook, W. S., Schmitt, R. W. & Páramo, P. 2010 Seismic imaging of a thermohaline staircase in the western tropical North Atlantic. Ocean Sci. 6, 621–631.
  * Garaud (2018) Garaud, P. 2018 Double-diffusive convection at low Prandtl number. Annu. Rev. Fluid Mech. 50, 275–298.
  * Garaud & Brummell (2015) Garaud, P. & Brummell, N. 2015 2D or not 2D: the effect of dimensionality on the dynamics of fingering convection at low Prandtl number. Astrophys. J. 815, 42.
  * Goluskin et al. (2014) Goluskin, D., Johnston, H., Flierl, G. R. & Spiegel, E. A. 2014 Convectively driven shear and decreased heat flux. J. Fluid Mech. 759, 360–385.
  * Gough et al. (1975) Gough, D. O., Spiegel, E. A. & Toomre, J. 1975 Modal equations for cellular convection. J. Fluid Mech. 68, 695–719.
  * Gough & Toomre (1982) Gough, D. O. & Toomre, J. 1982 Single-mode theory of diffusive layers in thermohaline convection. J. Fluid Mech. 125, 75–97.
  * Graham & Floryan (2021) Graham, M. D. & Floryan, D. 2021 Exact coherent states and the nonlinear dynamics of wall-bounded turbulent flows. Annu. Rev. Fluid Mech. 53, 227–253.
  * Hage & Tilgner (2010) Hage, E. & Tilgner, A. 2010 High Rayleigh number convection with double diffusive fingers. Phys. Fluids 22, 076603\.
  * Herring (1963) Herring, J. R. 1963 Investigation of problems in thermal convection. J. Atmos. Sci. 20, 325–338.
  * Herring (1964) Herring, J. R. 1964 Investigation of problems in thermal convection: rigid boundaries. J. Atmos. Sci. 21, 277–290.
  * Holyer (1981) Holyer, J. Y. 1981 On the collective instability of salt fingers. J. Fluid Mech. 110, 195–207.
  * Holyer (1984) Holyer, J. Y. 1984 The stability of long, steady, two-dimensional salt fingers. J. Fluid Mech. 147, 169–185.
  * Howard & Krishnamurti (1986) Howard, L. N. & Krishnamurti, R. 1986 Large-scale flow in turbulent convection: a mathematical model. J. Fluid Mech. 170, 385–410.
  * Julien & Knobloch (1997) Julien, K. & Knobloch, E. 1997 Fully nonlinear oscillatory convection in a rotating layer. Phys. Fluids 9, 1906–1913.
  * Julien & Knobloch (2007) Julien, K. & Knobloch, E. 2007 Reduced models for fluid flows with strong constraints. J. Math. Phys. 48, 065405\.
  * Julien et al. (1999) Julien, K., Knobloch, E. & Tobias, S. 1999 Strongly nonlinear magnetoconvection in three dimensions. Physica D 128, 105–129.
  * Julien et al. (2000) Julien, K., Knobloch, E. & Tobias, S. M. 2000 Nonlinear magnetoconvection in the presence of strong oblique fields. J. Fluid Mech. 410, 285–322.
  * Kawahara & Kida (2001) Kawahara, G. & Kida, S. 2001 Periodic motion embedded in plane Couette turbulence: regeneration cycle and burst. J. Fluid Mech. 449, 291–300.
  * Kawahara et al. (2012) Kawahara, G., Uhlmann, M. & van Veen, L. 2012 The significance of simple invariant solutions in turbulent flows. Annu. Rev. Fluid Mech. 44, 203–225.
  * Knobloch et al. (1986) Knobloch, E., Deane, A. E., Toomre, J. & Moore, D. R. 1986 Doubly diffusive waves. Contemp. Math. 56, 203–216.
  * Knobloch et al. (1992) Knobloch, E., Proctor, M. R. E. & Weiss, N. O. 1992 Heteroclinic bifurcations in a simple model of double-diffusive convection. J. Fluid Mech. 239, 273–292.
  * Krishnamurti (2003) Krishnamurti, R. 2003 Double-diffusive transport in laboratory thermohaline staircases. J. Fluid Mech. 483, 287–314.
  * Krishnamurti (2009) Krishnamurti, R. 2009 Heat, salt and momentum transport in a laboratory thermohaline staircase. J. Fluid Mech. 638, 491–506.
  * Krishnamurti & Howard (1981) Krishnamurti, R. & Howard, L. N. 1981 Large-scale flow generation in turbulent convection. Proc. Natl. Acad. Sci. 78 (4), 1981–1985.
  * Lewis et al. (1997) Lewis, S., Rees, D. A. S. & Bassom, A. P. 1997 High wavenumber convection in tall porous containers heated from below. Q. J. Mech. Appl. Math. 50, 545–563.
  * Li & Yang (2022) Li, J. & Yang, Y. 2022 On the wall-bounded model of fingering double diffusive convection. arXiv preprint p. arXiv:2204.03142.
  * Linden (1978) Linden, P. F. 1978 The formation of banded salt finger structure. J. Geophys. Res.: Oceans 83 (C6), 2902–2912.
  * Lopez & Murphy (1983) Lopez, J. M. & Murphy, J. O. 1983 Time-dependent thermal convection. Publ. Astron. Soc. Aust. 5, 173–175.
  * Lucas et al. (2017) Lucas, D., Caulfield, C. P. & Kerswell, R. R. 2017 Layer formation in horizontally forced stratified turbulence: connecting exact coherent structures to linear instabilities. J. Fluid Mech. 832, 409–437.
  * Lucas et al. (2019) Lucas, D., Caulfield, C. P. & Kerswell, R. R. 2019 Layer formation and relaminarisation in plane Couette flow with spanwise stratification. J. Fluid Mech. 868, 97–118.
  * Malkus & Veronis (1958) Malkus, W. V. R. & Veronis, G. 1958 Finite amplitude cellular convection. J. Fluid Mech. 4, 225–260.
  * Massaguer & Mercader (1988) Massaguer, J. M. & Mercader, I. 1988 Instability of swirl in low-Prandtl-number thermal convection. J. Fluid Mech. 189, 367–395.
  * Massaguer et al. (1990) Massaguer, J. M., Mercader, I. & Net, M. 1990 Nonlinear dynamics of vertical vorticity in low-Prandtl-number thermal convection. J. Fluid Mech. 214, 579–597.
  * Matthews et al. (1993) Matthews, P. C., Proctor, M. R. E., Rucklidge, A. M. & Weiss, N. O. 1993 Pulsating waves in nonlinear magnetoconvection. Phys. Lett. A 183, 69–75.
  * Morell et al. (2006) Morell, J. M., Corredor, J. E. & Merryfield, W. J. 2006 Thermohaline staircases in a Caribbean eddy and mechanisms for staircase formation. Deep Sea Res. Part II 53, 128–139.
  * Muench et al. (1990) Muench, R. D., Fernando, H. J. S. & Stegen, G. R. 1990 Temperature and salinity staircases in the northwestern Weddell Sea. J. Phys. Oceanogr. 20, 295–306.
  * Murphy & Lopez (1984) Murphy, J. O. & Lopez, J. M. 1984 The influence of vertical vorticity on thermal convection. Aust. J. Phys. 37, 179–196.
  * Oglethorpe et al. (2013) Oglethorpe, R. L. F., Caulfield, C. P. & Woods, A. W. 2013 Spontaneous layering in stratified turbulent Taylor–Couette flow. J. Fluid Mech. 721, R3.
  * Padman & Dillon (1989) Padman, L. & Dillon, T. M. 1989 Thermal microstructure and internal waves in the Canada Basin diffusive staircase. Deep Sea Res. Part A 36, 531–542.
  * Paparella (1997) Paparella, F. 1997 A few steps toward staircases. Double-Diffusive Processes, 1996 Summer Study Program in Geophysical Fluid Dynamics pp. 232–247.
  * Paparella & Spiegel (1999) Paparella, F. & Spiegel, E. A. 1999 Sheared salt fingers: Instability in a truncated system. Phys. Fluids 11, 1161–1168.
  * Paparella et al. (2002) Paparella, F., Spiegel, E. A. & Talon, S. 2002 Shear and mixing in oscillatory doubly diffusive convection. Geophys. Astrophys. Fluid Dyn. 96, 271–289.
  * Phillips (1972) Phillips, O. M. 1972 Turbulence in a strongly stratified fluid—Is it unstable? In Deep Sea Res. Oceanogr. Abstr., , vol. 19, pp. 79–81. Elsevier.
  * Piacsek & Toomre (1980) Piacsek, S. A. & Toomre, J. 1980 Nonlinear evolution and structure of salt fingers. Marine Turbulence, Proceedings of The 11th International Liège Colloquium on Ocean Hydrodynamics 28, 193–219.
  * Plumley et al. (2018) Plumley, M., Calkins, M. A., Julien, K. & Tobias, S. M. 2018 Self-consistent single mode investigations of the quasi-geostrophic convection-driven dynamo model. J. Plasma Phys. 84, 735840406.
  * van der Poel et al. (2013) van der Poel, E. P., Stevens, R. J. A. M. & Lohse, D. 2013 Comparison between two-and three-dimensional Rayleigh–Bénard convection. J. Fluid Mech. 736, 177–194.
  * Posmentier (1977) Posmentier, E. S. 1977 The generation of salinity finestructure by vertical diffusion. J. Phys. Oceanogr. 7, 298–300.
  * Proctor & Holyer (1986) Proctor, M. R. E. & Holyer, J. Y. 1986 Planform selection in salt fingers. J. Fluid Mech. 168, 241–253.
  * Proctor & Weiss (1993) Proctor, M. R. E. & Weiss, N. O. 1993 Symmetries of time-dependent magnetoconvection. Geophys. Astrophys. Fluid Dyn. 70, 137–160.
  * Proctor et al. (1994) Proctor, M. R. E., Weiss, N. O., Brownjohn, D. P. & Hurlburt, N. E. 1994 Nonlinear compressible magnetoconvection Part 2. Streaming instabilities in two dimensions. J. Fluid Mech. 280, 227–253.
  * Rademacher & Uecker (2017) Rademacher, J. D. M. & Uecker, H. 2017 Symmetries, freezing, and Hopf bifurcations of traveling waves in pde2path. https://www.staff.uni-oldenburg.de/hannes.uecker/pde2path/ .
  * Radko (2003) Radko, T. 2003 A mechanism for layer formation in a double-diffusive fluid. J. Fluid Mech. 497, 365–380.
  * Radko (2005) Radko, T. 2005 What determines the thickness of layers in a thermohaline staircase? J. Fluid Mech. 523, 79–98.
  * Radko (2010) Radko, T. 2010 Equilibration of weakly nonlinear salt fingers. J. Fluid Mech. 645, 121–143.
  * Radko (2013) Radko, T. 2013 Double-diffusive Convection. Cambridge University Press.
  * Radko (2016) Radko, T. 2016 Thermohaline layering in dynamically and diffusively stable shear flows. J. Fluid Mech. 805, 147–170.
  * Radko & Stern (1999) Radko, T. & Stern, M. E. 1999 Salt fingers in three dimensions. J. Mar. Res. 57, 471–502.
  * Radko & Stern (2000) Radko, T. & Stern, M. E. 2000 Finite-amplitude salt fingers in a vertically bounded layer. J. Fluid Mech. 425, 133–160.
  * Rhines & Young (1983) Rhines, P. B. & Young, W. R. 1983 How rapidly is a passive scalar mixed within closed streamlines? J. Fluid Mech. 133, 133–145.
  * Rucklidge & Matthews (1996) Rucklidge, A. M. & Matthews, P. C. 1996 Analysis of the shearing instability in nonlinear convection and magnetoconvection. Nonlinearity 9, 311–351.
  * Schmid & Henningson (2012) Schmid, P. J. & Henningson, D. S. 2012 Stability and Transition in Shear Flows. Springer Science & Business Media.
  * Schmitt et al. (2005) Schmitt, R. W., Ledwell, J. R., Montgomery, E. T., Polzin, K. L. & Toole, J. M. 2005 Enhanced diapycnal mixing by salt fingers in the thermocline of the tropical Atlantic. Science 308 (5722), 685–688.
  * Schmitt et al. (1987) Schmitt, R. W., Perkins, H., Boyd, J. D. & Stalcup, M. C. 1987 C-SALT: An investigation of the thermohaline staircase in the western tropical North Atlantic. Deep Sea Res. Part A 34, 1655–1665.
  * Siggia (1994) Siggia, E. D. 1994 High Rayleigh number convection. Annu. Rev. Fluid Mech. 26, 137–168.
  * Spear & Thomson (2012) Spear, D. J. & Thomson, R. E. 2012 Thermohaline staircases in a British Columbia fjord. Atmos. Ocean 50, 127–133.
  * St. Laurent & Schmitt (1999) St. Laurent, L. & Schmitt, R. W. 1999 The contribution of salt fingers to vertical mixing in the North Atlantic Tracer Release Experiment. J. Phys. Oceanogr. 29, 1404–1424.
  * Stellmach et al. (2011) Stellmach, S., Traxler, A., Garaud, P., Brummell, N. & Radko, T. 2011 Dynamics of fingering convection. Part 2 The formation of thermohaline staircases. J. Fluid Mech. 677, 554–571.
  * Stern (1969) Stern, M. E. 1969 Collective instability of salt fingers. J. Fluid Mech. 35, 209–218.
  * Tait & Howe (1968) Tait, R. I. & Howe, M. R. 1968 Some observations of thermo-haline stratification in the deep ocean. Deep Sea Res. Oceanogr. Abstr. 15, 275–280.
  * Tait & Howe (1971) Tait, R. I. & Howe, M. R. 1971 Thermohaline staircase. Nature 231 (5299), 178–179.
  * Taylor & Bucens (1989) Taylor, J. & Bucens, P. 1989 Laboratory experiments on the structure of salt fingers. Deep Sea Res. Part A 36, 1675–1704.
  * Taylor & Zhou (2017) Taylor, J. R. & Zhou, Q. 2017 A multi-parameter criterion for layer formation in a stratified shear flow using sorted buoyancy coordinates. J. Fluid Mech. 823, R5.
  * Timmermans et al. (2008) Timmermans, M.-L., Toole, J., Krishfield, R. & Winsor, P. 2008 Ice-Tethered Profiler observations of the double-diffusive staircase in the Canada Basin thermocline. J. Geophys. Res. Oceans 113, C00A02.
  * Toomre et al. (1977) Toomre, J., Gough, D. O. & Spiegel, E. A. 1977 Numerical solutions of single-mode convection equations. J. Fluid Mech. 79, 1–31.
  * Uecker (2021a) Uecker, H. 2021a Numerical Continuation and Bifurcation in Nonlinear PDEs. SIAM.
  * Uecker (2021b) Uecker, H. 2021b pde2path without finite elements. https://www.staff.uni-oldenburg.de/hannes.uecker/pde2path/ .
  * Uecker et al. (2014) Uecker, H., Wetzel, D. & Rademacher, J. D. M. 2014 pde2path: A Matlab package for continuation and bifurcation in 2D elliptic systems. Numer. Math. Theory Methods Appl. 7, 58–106.
  * van Veen et al. (2006) van Veen, L., Kida, S. & Kawahara, G. 2006 Periodic motion representing isotropic turbulence. Fluid Dyn. Res. 38, 19–46.
  * Wang et al. (2020) Wang, Q., Verzicco, R., Lohse, D. & Shishkina, O. 2020 Multiple states in turbulent large-aspect-ratio thermal convection: What determines the number of convection rolls? Phys. Rev. Lett. 125, 074501.
  * Weideman & Reddy (2000) Weideman, J. A. & Reddy, S. C. 2000 A MATLAB differentiation matrix suite. ACM Trans. Math. Software 26, 465–519.
  * Xia et al. (2018) Xia, Z., Shi, Y., Cai, Q., Wan, M. & Chen, S. 2018 Multiple states in turbulent plane Couette flow with spanwise rotation. J. Fluid Mech. 837, 477–490.
  * Xie et al. (2019) Xie, J.-H., Julien, K. & Knobloch, E. 2019 Jet formation in salt-finger convection: a modified Rayleigh–Bénard problem. J. Fluid Mech. 858, 228–263.
  * Xie et al. (2017) Xie, J.-H., Miquel, B., Julien, K. & Knobloch, E. 2017 A reduced model for salt-finger convection in the small diffusivity ratio limit. Fluids 2 (1), 6.
  * Yang & Xia (2021) Yang, X. & Xia, Z.-h. 2021 Bifurcation and multiple states in plane Couette flow with spanwise rotation. J. Fluid Mech. 913, A49.
  * Yang et al. (2020) Yang, Y., Chen, W., Verzicco, R. & Lohse, D. 2020 Multiple states and transport properties of double-diffusive convection turbulence. Proc. Natl. Acad. Sci. 117 (26), 14676–14681.
  * Yang et al. (2015) Yang, Y., van der Poel, E. P., Ostilla-Mónico, R., Sun, C., Verzicco, R., Grossmann, S. & Lohse, D. 2015 Salinity transfer in bounded double diffusive convection. J. Fluid Mech. 768, 476–491.
  * Yang et al. (2016a) Yang, Y., Verzicco, R. & Lohse, D. 2016a Scaling laws and flow structures of double diffusive convection in the finger regime. J. Fluid Mech. 802, 667–689.
  * Yang et al. (2016b) Yang, Y., Verzicco, R. & Lohse, D. 2016b Vertically bounded double diffusive convection in the finger regime: Comparing no-slip versus free-slip boundary conditions. Phys. Rev. Lett. 117, 184501.
  * Yang et al. (2022) Yang, Y., Verzicco, R., Lohse, D. & Caulfield, C. P. 2022 Layering and vertical transport in sheared double-diffusive convection in the diffusive regime. J. Fluid Mech. 933, A30.
  * Zahn et al. (1974) Zahn, J.-P., Toomre, J., Spiegel, E. A. & Gough, D. O. 1974 Nonlinear cellular motions in Poiseuille channel flow. J. Fluid Mech. 64, 319–346.
  * Zhang et al. (2018) Zhang, X., Wang, L.-l., Lin, C., Zhu, H. & Zeng, C. 2018 Numerical study on tilting salt finger in a laminar shear flow. Phys. Fluids 30, 022110.
  * Zodiatis & Gasparini (1996) Zodiatis, G. & Gasparini, G. P. 1996 Thermohaline staircase formations in the Tyrrhenian Sea. Deep Sea Res. Part I 43, 655–678.