A Reduced Model for Salt-Finger Convection in the Small Diffusivity Ratio Limit

Abstract

A simple model of nonlinear salt-finger convection in two dimensions is derived and studied. The model is valid in the limit of a small solute to heat diffusivity ratio and a large density ratio, which is relevant to both oceanographic and astrophysical applications. Two limits distinguished by the magnitude of the Schmidt number are found. For order one Schmidt numbers, appropriate for astrophysical applications, a modified Rayleigh–Bénard system with large-scale damping due to a stabilizing temperature is obtained. For large Schmidt numbers, appropriate for the oceanic setting, the model combines a prognostic equation for the solute field and a diagnostic equation for inertia-free momentum dynamics. Two distinct saturation regimes are identified for the second model: the weakly driven regime is characterized by a large-scale flow associated with a balance between advection and linear instability, while the strongly-driven regime produces multiscale structures, resulting in a balance between energy input through linear instability and energy transfer between scales. For both regimes, we analytically predict and numerically confirm the dependence of the kinetic energy and salinity fluxes on the ratio between solutal and thermal Rayleigh numbers. The spectra and probability density functions are also computed.

1. Introduction

Doubly-diffusive systems in which two components with different diffusivities contribute to buoyancy in opposite ways arise frequently in geophysics and astrophysics [1,2,3]: for example, heat and salt contribute to the stratification in the oceanic setting, while heat and chemical composition contribute in stellar astrophysics and sugar and salt in laboratory experiments. In the following, motivated by oceanic dynamics, we refer to the slowly-diffusing component as salinity and to the rapidly-diffusing component as temperature. Under appropriate conditions, the different diffusivities of these two components may destabilize an otherwise stably-stratified configuration and, hence, lead to mixing. Two scenarios are possible. In the first, temperature stratification is stabilizing, but the salinity stratification is destabilizing, a configuration that leads to a fingering instability. In the second scenario, where the salinity stratification is stabilizing while the temperature is destabilizing, a diffusive oscillatory instability called overstability may take place. In both of these cases, instability is present because heat diffuses more rapidly than a typical solute, i.e., the ratio τ of the solute diffusivity to thermal diffusivity (or inverse Lewis number) satisfies $\tau <1$. In this paper, we concentrate on the nonlinear properties of the turbulent state arising in the former case and suppose that the destabilizing effect is provided by an unstable salt stratification as frequently occurs in the oceans, focusing on the case $\tau \ll 1$.

Linear stability properties of the salt-finger regime were first determined by Stern [4]. As the magnitude of the unstable mode grows, secondary instabilities are triggered [5] and, finally, lead to a statistically steady state [6,7,8]. However, the physical processes behind the presence of such an equilibrium state are still incompletely understood. Stern [9] has suggested that a collective instability can lead to the disruption of the salt-finger field and lead to saturation and identified a dimensionless quantity, now called the Stern number, that can be used as a criterion for saturation. Others have supposed that saturation will arise once the growth rate of a secondary instability, such as that identified in [5], becomes comparable to that of the salt-finger instability itself [8,10]. More recently Shen [6] has suggested that finger collision plays an essential role because of the generation of large vertical gradients and, hence, strong dissipation. In an interesting paper, Paparella and von Hardenberg [11] proposed a saturation mechanism based on finger clustering that leads to the generation of large-scale structures and, hence, to saturation at the large scale.

Except for weakly nonlinear approaches, e.g., [7,12], most researchers employ numerical studies of the primitive equations consisting of the Navier–Stokes equation coupled to scalar equations for heat and salinity transport. Such studies, particularly in three dimensions, have proven invaluable and identified a number of novel processes, e.g., [13]. To shed light on some of these, we propose here a systematic procedure that leads to a simplified set of equations that is easier to study, both theoretically and numerically. In the following, we refer to these, following [14], as reduced equations. The procedure described below has been used before [15,16,17] and focuses on the double limit of a large density ratio, where the stratification is dominated by thermal contributions, and an asymptotically small diffusivity ratio, where salinity diffuses much more slowly than heat. The latter is motivated by the small value of the diffusivity ratio in the ocean while large density ratios are observed in fingering layers in experiments [18] and in some oceanic measurements [19], and is a limit beyond the current capability of direct numerical simulations. Two reduced models are derived depending on the order of the Schmidt number $\mathrm{Sc}$, which captures the relative importance of momentum and salinity diffusion. When these two effects are comparable, $\mathrm{Sc}=O(1)$, a reduced model that bears similarity to the Rayleigh–Bénard configuration is derived. This model is applicable to astrophysical objects, e.g., [10,17,20]. In the oceanic parameter regime with fast momentum diffusion, $\mathrm{Sc}\gg 1$, we derive a reduced model with a prognostic-diagnostic form, which we refer to as the inertia-free salt convection model. This model is the main topic of this paper.

It turns out that this model already appears in the work of Radko and Stern [16]. However, beyond one three-dimensional simulation in a very small aspect ratio domain, these authors did not investigate the properties of this model. We believe in fact that the model merits far greater attention and that it correctly captures the saturation of the salt-finger instability and the resulting transport properties in a physically-relevant regime. In this paper, we therefore study the properties of this model at some length, albeit in two dimensions. To do so, we employ a doubly-periodic setting as appropriate for oceanic and astrophysical applications with distant boundaries, in contrast to earlier studies of salt-finger convection in vertically-bounded domains, e.g., [21,22,23,24,25]. In this formulation, elevator modes consisting of vertically-invariant upward and downward moving fingers are exact nonlinear solutions, and we believe that analogous modes are crucial even in vertically-bounded domains because of their efficiency in extracting energy from the salinity field. We mention that the two-dimensional model we study is expected to capture the dynamics of the three-dimensional system for $O(1)$ Prandtl numbers and sufficiently large density ratios [26], but not for low Prandtl numbers for which simulations of the two-dimensional and three-dimensional primitive equations lead to qualitatively different predictions [27].

The structure of this paper is as follows. In Section 2, we formulate the salt-finger convection problem, followed in Section 3 by a derivation of two distinct reduced models valid for finite and infinite Schmidt numbers, respectively. The time evolution and the properties of the saturated state of the latter are studied in Section 4 and Section 5 through a combination of numerical simulation and scaling analysis. Our results are summarized in Section 6. Two appendices, one summarizing the dependence of the results on the domain size and the other an analytical approximation to one of the saturation regimes, complete the paper.

2. Formulation

2.1. Setup and Nondimensionalization

We consider a two-dimensional Boussinesq fluid of infinite extent in the horizontal and vertical directions, denoted by x and z, respectively. An initial uniform density stratification is generated by means of a stabilizing background temperature gradient ${\beta }_T>0$ and a destabilizing background salinity gradient ${\beta }_S>0$. We introduce the thermal and saline expansion coefficients ${\alpha }_T>0$ and ${\alpha }_S>0$ so that the density of the fluid $\rho (T,S)$ at temperature T and salinity S is given by $\rho (T,S)={\rho }_0(1-{\alpha }_T(T-T_0)+{\alpha }_S(S-S_0))$, where ${\rho }_0=\rho (T_0,S_0)$ is a reference state. The remaining material properties of the fluid are specified by the thermal diffusivity ${\kappa }_T$, the solute diffusivity ${\kappa }_S$ and the viscosity ν, all of which are taken to be constant. Acceleration due to gravity is denoted by g.

The equations are made dimensionless using the natural salt-finger width d given by:

$$d={\left(\frac{\nu {\kappa }_T}{g{\alpha }_T{\beta }_T}\right)}^{1/4}.$$

(1)

This length defines the characteristic temperature $\mathcal{T}$ and salinity $\mathcal{S}$ scales:

$$\begin{array}{cc}\hfill \mathcal{T}& ={\beta }_{T}d,\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill \mathcal{S}& ={\beta }_{S}d.\hfill \end{array}$$

(2b)

These scales are used to nondimensionalize the temperature and salinity profiles $T_{\mathrm{total}}$ and $S_{\mathrm{total}}$, which are decomposed into the initial background state and a fluctuating component indicated by a tilde:

$$\begin{array}{c}{T}_{\mathrm{total}}=z+\tilde{T},\end{array}$$

(3a)

$$\begin{array}{c}{S}_{\mathrm{total}}=z+\tilde{S}.\end{array}$$

(3b)

The incompressible velocity field is expressed in terms of a stream function ψ,

$$\mathit{u}\equiv u{\mathit{e}}_x+w{\mathit{e}}_z=\left(\begin{array}{c}-{\partial }_z\psi \\ {\partial }_x\psi \end{array}\right),$$

(4)

where $(u,w)$ and ψ are made dimensionless using ${\kappa }_S/d$ and ${\kappa }_S$, respectively. Finally, time is made dimensionless using the salinity diffusion timescale:

$$t_{\mathrm{salt}}=\frac{d^2}{{\kappa }_S}.$$

(5)

2.2. Nondimensional Equations

The governing dimensionless equations for fluctuations around the conduction state $\psi =0$, $T_{\mathrm{total}}=S_{\mathrm{total}}=z$ are:

$$\begin{array}{cc}\hfill \frac{\tau }{Pr}\left[\frac{\partial }{\partial t}{\nabla }^2\psi +\mathcal{J}\left(\psi ,{\nabla }^2\psi \right)\right]& =\frac{1}{\tau }\frac{\partial \tilde{T}}{\partial x}-\frac{1}{\tau {\mathrm{R}}_{\rho }}\frac{\partial \tilde{S}}{\partial x}+{\nabla }^4\psi ,\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}\tilde{T}+\mathcal{J}\left(\psi ,\tilde{T}\right)+\frac{\partial \psi }{\partial x}& ={\tau }^{-1}{\nabla }^2\tilde{T},\hfill \end{array}$$

(6b)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}\tilde{S}+\mathcal{J}\left(\psi ,\tilde{S}\right)+\frac{\partial \psi }{\partial x}& ={\nabla }^2\tilde{S}.\hfill \end{array}$$

(6c)

Three dimensionless ratios enter this set of equations: the Prandtl number Pr, the inverse Lewis number τ and the density ratio ${\mathrm{R}}_{\rho }$ specifying the relative contribution of the temperature and salinity to the background density gradient. A fourth dimensionless ratio, the Schmidt number $\mathrm{Sc}$, can be introduced as an alternative to the diffusivity ratio. These quantities are defined by:

$$Pr=\frac{\nu }{{\kappa }_T},\tau =\frac{{\kappa }_S}{{\kappa }_T},\mathrm{Sc}=\frac{\nu }{{\kappa }_S},{\mathrm{R}}_{\rho }=\frac{{\alpha }_T{\beta }_T}{{\alpha }_S{\beta }_S}.$$

(7)

In this paper, we adopt doubly-periodic boundary conditions for ψ, $\tilde{T}$ and $\tilde{S}$ and consider a statically-stable configuration with ${\mathrm{R}}_{\rho }>1$.

2.3. Connection with an Alternative Description in Terms of Flux Rayleigh Numbers

The above formulation can be connected to the classical description of doubly-diffusive convection in terms of the temperature and salinity Rayleigh numbers, ${\mathrm{Ra}}_T$ and ${\mathrm{Ra}}_S$, respectively. These are defined for a layer of vertical extent H, viz.,

$${\mathrm{Ra}}_T=\frac{g{\alpha }_T{\beta }_T{H}^4}{\nu {\kappa }_T}\mathrm{and}{\mathrm{Ra}}_S=\frac{g{\alpha }_S{\beta }_S{H}^4}{\nu {\kappa }_S},$$

(8)

implying that:

$${\mathrm{R}}_{\rho }\tau =\frac{{\mathrm{Ra}}_T}{{\mathrm{Ra}}_S}.$$

(9)

Large ${\mathrm{Ra}}_T$ flows correspond to slender fingers (at least in the linear regime) as the aspect ratio $d/H$ of the fingers is given by:

$$\frac{d}{H}={\mathrm{Ra}}_T^{-1/4}.$$

(10)

2.4. Linear Stability Analysis

The growth of salt fingers can be understood by analyzing the linear stability of the conduction state with respect to normal mode perturbations of the form $(\psi (x,z,t),\tilde{T}(x,z,t),\tilde{S}(x,z,t))=\Re \left\{({\psi }_e,T_e,S_e){\mathrm{e}}^{\lambda t+\mathrm{i}(kx+mz)}\right\}$. The growth rate λ obeys the dispersion relation:

$${\lambda }^3+{|K|}^2\left(1+\frac{1}{\tau }+\frac{Pr}{\tau }\right){\lambda }^2+\left[{|K|}^4\left(\frac{1}{\tau }+\frac{Pr}{\tau }+\frac{Pr}{{\tau }^2}\right)+\frac{Pr{k}^2}{{\tau }^2{|K|}^2}\left(1-\frac{1}{{\mathrm{R}}_{\rho }}\right)\right]\lambda +\frac{Pr}{{\tau }^2}\left[{|K|}^6+k^2\left(1-\frac{1}{\tau {\mathrm{R}}_{\rho }}\right)\right]=0,$$

(11)

where ${|K|}^2\equiv k^2+m^2$. For a Rayleigh–Taylor stable layer such that ${\mathrm{R}}_{\rho }>1$, the cubic, quadratic and linear coefficients of this dispersion relation are positive. An exponentially-growing instability is present if and only if there exists a wavenumber k such that the constant term is negative. This condition is equivalent to:

$$\frac{1}{\tau {\mathrm{R}}_{\rho }}>1.$$

(12)

The salt-finger instability is therefore present in the regime (e.g., [3]):

$$1<{\mathrm{R}}_{\rho }<\frac{1}{\tau }.$$

(13)

This condition defines the regime of interest for the remainder of this paper.

The strongest linear instability is associated with vertically-invariant modes with $m=0$, called elevator modes in the literature. With periodic boundary conditions in the vertical, these modes are exact solutions of the system (6) and so play an important role in the dynamics. In a vertically-bounded domain, elevator modes are not permitted, but their effect persists, at least when the domain has a large vertical extent, $H\gg d$.

3. Reduced Models with a Small Diffusivity Ratio

3.1. Scaling of Parameters and Variables

In this section, we focus on the small diffusivity ratio limit $\tau \ll 1$ of the system (6), i.e., ${\kappa }_S\ll {\kappa }_T$. Two cases are considered, distinguished by the magnitude of the Schmidt number $\mathrm{Sc}$. For a Schmidt number of order one, where ${\kappa }_S$ ∼ ν, we obtain the small Prandtl number regime $\mathrm{Pr}=O(\tau )$ of astrophysical relevance. This results in a modified Rayleigh–Bénard system with salinity-driven instability and rapidly-diffusing temperature, as described in Section 3.2. The large Schmidt number regime, where ${\kappa }_S\ll \nu$ and $\mathrm{Pr}\gg \tau$, is relevant for oceanic thermohaline flows where $\mathrm{Sc}\approx 700$ and $\tau \approx 0.01$. This results in a model where inertial forces are small and salinity is the only slowly-diffusing quantity:

$$\begin{array}{c}{\kappa }_S\ll {\kappa }_T,\end{array}$$

(14a)

$$\begin{array}{c}{\kappa }_S\ll \nu .\end{array}$$

(14b)

We refer to this system, considered in Section 3.3, as the inertia-free salt convection model.

In all cases, we suppose, moreover, that the density ratio ${\mathrm{R}}_{\rho }$ is large and comparable to ${\tau }^{-1}$. This assumption is supported by both laboratory experiments in the fingering domain [28] and oceanic measurements [19] where ${\mathrm{R}}_{\rho }$ can reach $O({10}^2)$. Furthermore, in the astrophysical context, e.g., for stars along the red giant branch (RGB), ${\mathrm{R}}_{\rho }$ can be as large as $O({10}^6)$, which is of the same order of magnitude as the corresponding inverse Prandtl number [29]. Hence, we define the order one parameter:

$$\mathrm{Ra}=\frac{1}{{\mathrm{R}}_{\rho }\tau }=\frac{{\mathrm{Ra}}_S}{{\mathrm{Ra}}_T}.$$

(15)

This parameter is analogous to a Rayleigh number for our system, as we demonstrate below, and we refer to it as the Rayleigh ratio. The parameter $\mathcal{R}$,

$$\mathcal{R}\equiv \mathrm{Ra}-1=\frac{1}{{\mathrm{R}}_{\rho }\tau }-1,$$

(16)

measures the supercriticality of the system and will therefore prove to be useful as well.

We now discuss the scaling of the fluid variables, starting with the stream function, which remains untouched. The temperature fluctuations $\tilde{T}$ are scaled with τ and become an asymptotically small quantity in the limit $\tau \downarrow 0$. In this limit, the salinity fluctuations $\tilde{S}$ remain finite, but are scaled for convenience by the order one factor ${\mathrm{Ra}}^{-1}$:

$$\begin{array}{c}{T}_{\mathrm{total}}=z+\tau T,\end{array}$$

(17a)

$$\begin{array}{c}{S}_{\mathrm{total}}=z+{\mathrm{Ra}}^{-1}S,\end{array}$$

(17b)

where T and S are now order one quantities. This formulation permits a feedback on the background salinity profile, but excludes any leading order modification of the background temperature profile.

3.2. Modified Rayleigh–Bénard Convection Model: $\mathrm{Sc}=O(1)$

In the limit of small τ and $\mathrm{Sc}=O(1)$, with the scaling $\tilde{T}=\tau T$, we obtain a distinguished regime described by the following set of equations:

$$\begin{array}{cc}\hfill \frac{1}{\mathrm{Sc}}\left[\frac{\partial }{\partial t}{\nabla }^2\psi +\mathcal{J}\left(\psi ,{\nabla }^2\psi \right)\right]& =-\frac{\partial S}{\partial x}+\left({\nabla }^4+{\Delta }^{-1}{\partial }_x^2\right)\psi ,\hfill \end{array}$$

(18a)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}S+\mathcal{J}\left(\psi ,S\right)+\mathrm{Ra}\frac{\partial \psi }{\partial x}& ={\nabla }^{2}S.\hfill \end{array}$$

(18b)

Here, the inverse Laplacian in (18a) comes from solving

$$\frac{\partial \psi }{\partial x}={\nabla }^{2}T$$

(19)

for the temperature. This equation describes the leading order balance in the temperature equation, Equation (6b), in the above-mentioned parameter regime and has previously been applied to both salt-finger convection [17] and temperature-stratified shear flow [30] within the same small Péclet number approximation. Equations (18) describe a modified Rayleigh–Bénard convection (mRBC) system with large-scale vorticity dissipation given by ${\Delta }^{-1}{\partial }_x^2\psi$ in addition to the usual small-scale vorticity dissipation via ${\nabla }^4\psi$. This new source of dissipation is brought about by the stabilizing temperature field and serves to select an intrinsic finite-scale optimal mode, in contrast to the standard RBC system where the scale of the optimal mode is determined by the domain height. The resulting model applies under conditions prevailing in astrophysics, for instance RGB stars, characterized by extremely small parameters $Pr=O(\tau )$ ∼ ${10}^{-6}$. This regime appears to be accessible in two-dimensional direct numerical simulations of the primitive doubly-diffusive equations [29], but not in three-dimensional simulations, which can reach, at present, only $Pr=O(\tau )$ ∼ ${10}^{-2}$ [20].

3.3. Inertia-Free Salt Convection Model: $\mathrm{Sc}\gg 1$

In this regime, characterized by Condition (14b), viscosity is greater than the diffusivity of the salinity, and as a result, inertial modes are suppressed. This regime is relevant to the oceans and yields a further reduction of Equations (18) to a model that is now first order in time:

$$\begin{array}{c}\left({\partial }_x^2+{\nabla }^6\right)\psi ={\partial }_x{\nabla }^{2}S,\end{array}$$

(20a)

$$\begin{array}{c}\frac{\partial }{\partial t}S+\mathcal{J}\left(\psi ,S\right)+\mathrm{Ra}\frac{\partial \psi }{\partial x}={\nabla }^{2}S.\end{array}$$

(20b)

The associated temperature field T is slaved to the stream function ψ as described by (19). In the following, we refer to these equations as the inertia-free salt convection (IFSC) model. This model is specified by a single external parameter $\mathrm{Ra}$, i.e., in this regime, the properties of salt-finger convection do not depend on the density ratio ${\mathrm{R}}_{\rho }$ and the diffusivity ratio τ separately, but on their combination $\mathrm{Ra}=1/({\mathrm{R}}_{\rho }\tau )$. This finding is similar to that of Lignières [30] and Jones [31], who showed that the effective parameter in stably-stratified shear flow is the parameter combination $\mathrm{Pe}\mathrm{Ri}$, where $\mathrm{Pe}\ll 1$ is the Péclet number and $\mathrm{Ri}\gg 1$ is the Richardson number, as well as Knobloch and Spruit [32], who showed that the stability of rotating stars is governed by the parameter combination $\mathrm{Pr}\mathrm{R}$, where $\mathrm{Pr}\ll 1$ is the Prandtl number and $\mathrm{R}\gg 1$ is a suitably-defined buoyancy number.

The IFSC model is the simplest model of salt-finger convection that can be derived and depends on the single intrinsic parameter $\mathrm{Ra}$ only. The coupled form of this model is reminiscent of other geophysical models, such as those describing β-plane vorticity dynamics [33], Bénard convection in the large Prandtl number limit [34] and β-convection [35]. However, the diagnostic equation, Equation (20a), exhibits two distinguishing features: (i) in contrast to the aforementioned models, the salinity and stream function are necessarily out of phase, a crucial fact for linear instability; (ii) the existence of a complex salinity-stream function relation where the two terms ${\partial }_x^2$ and ${\nabla }^6$ dominate on different horizontal scales. Consequently, (20a) implies the potential for dynamics on multiple spatial scales, a fact we believe provides the key to explaining the observed transitions in the saturated state as $\mathrm{Ra}$ increases; see Section 5.1.

In view of its simplicity and its strong relevance to oceanic flows, we focus in this paper on the IFSC model (20) and leave the mRBC model (18) for future study.

4. IFSC Model: From Linear Instability to the Saturated State

In this section, we document a scenario that leads from vanishingly small perturbations of the purely conductive state to a convective saturated state in the IFSC model. During the first stage, a linear instability occurs and preferentially amplifies vertically-invariant structures (elevator modes). These structures are then subject to a secondary instability, resulting in undulations in the vertical direction. Finally, the system saturates and exhibits structures of a much shorter vertical extent than observed during the first two phases.

4.1. Linear Stability

Substituting a normal mode ansatz of the form ${\mathrm{e}}^{\lambda t+\mathrm{i}(kx+mz)}$ into Equations (20) linearized around $S=\psi =0$, we obtain the growth rate λ as a function of $\mathrm{Ra}$:

$$\lambda =\mathrm{Ra}\frac{k^2{|K|}^2}{k^2+{|K|}^6}-{|K|}^2,$$

(21)

where, as before, ${|K|}^2=k^2+m^2$. It follows that the threshold for instability, determined in Equation (12), is preserved in the IFSC model and occurs at:

$${\mathrm{Ra}}_c=1.$$

(22)

Thus, the quantity $\mathcal{R}\equiv \mathrm{Ra}-1$ measures the distance from the threshold, i.e., the supercriticality of the system.

We show the linear growth rate $\lambda (k,m)$ in Figure 1 for (a) $\mathrm{Ra}=1.1$ and (b) $\mathrm{Ra}=10$. Figure 1 highlights the anisotropic dependence of λ on k and m: for fixed m, not too large, there exists an optimal k that maximizes the growth rate, whereas when k is fixed, an optimal m exists for small k, but for large k, the growth rate decreases monotonically with increasing m. For the elevator modes with zero vertical wavenumber, $m=0$, the band of unstable horizontal wavenumbers is well captured by the supercriticality $\mathcal{R}$:

$$0<k^4<\mathcal{R}.$$

(23)

As already mentioned, these modes are exact solutions of the nonlinear system (20) with periodic boundary conditions in the vertical and so play a potentially important role in the nonlinear regime. We present the growth rates of the elevator modes in Figure 2 for several values of $\mathrm{Ra}$.

In the absence of boundaries, the fastest growing mode is an elevator mode, such that:

$$k_{\mathrm{opt}}^4=\frac{1}{2}\left(-2-\mathrm{Ra}+\sqrt{{\mathrm{Ra}}^2+8\mathrm{Ra}}\right)\mathrm{and}{m}_{\mathrm{opt}}=0,$$

(24)

with growth rate given by:

$${\lambda }_{\mathrm{opt}}=\sqrt{\frac{1}{2}\left(-2-\mathrm{Ra}+\sqrt{{\mathrm{Ra}}^2+8\mathrm{Ra}}\right)}{\displaystyle \frac{3\mathrm{Ra}-\sqrt{{\mathrm{Ra}}^2+8\mathrm{Ra}}}{\sqrt{{\mathrm{Ra}}^2+8\mathrm{Ra}}-\mathrm{Ra}}}.$$

(25)

In Figure 3, we re-plot these results in a log-log plot, indicating the asymptotic limits valid for $\mathcal{R}\ll 1$ and $\mathcal{R}\to \infty$:

$$\begin{array}{cc}\hfill \underset{\mathcal{R}\to 0}{lim}{k}_{\mathrm{opt}}& =\frac{1}{3^{1/4}}{\mathcal{R}}^{1/4},\underset{\mathrm{Ra}\to \infty }{lim}{k}_{\mathrm{opt}}=1;\hfill \end{array}$$

(26a)

$$\begin{array}{cc}\hfill \underset{\mathcal{R}\to 0}{lim}{\lambda }_{\mathrm{opt}}& =\frac{2}{3^{3/2}}{\mathcal{R}}^{3/2},\underset{\mathrm{Ra}\to \infty }{lim}{\lambda }_{\mathrm{opt}}=\frac{\mathrm{Ra}}{2}.\hfill \end{array}$$

(26b)

When stress-free, fixed temperature and fixed salinity boundary conditions are applied on horizontal boundaries separated vertically by a gap $2\pi$, the marginal stability curve corresponds to $m_{\mathrm{opt}}=1$ and yields:

$${\mathrm{Ra}}_{\mathrm{mar}}=1+\frac{{\left(k^2+1\right)}^3}{k^2},$$

(27)

a result highly reminiscent of the marginal stability curve of the classic Bénard problem with similar boundary conditions: ${\mathrm{Ra}}_{\mathrm{mar}}^{(\mathrm{RBC})}={\left(k^2+1\right)}^3/k^2$ [36].

We conclude this section with an important energy argument. In a doubly-periodic domain with period $L_x$ in the horizontal direction and $L_z$ in the vertical direction, we define the spatially-averaged salinity potential energy $E_S$ as:

$$E_S=\frac{1}{L_x{L}_z}\int \frac{1}{2}{S}^2(x,z,t)\mathrm{d}x\mathrm{d}z.$$

(28)

For $\mathit{k}\in \left(2\pi n_x/L_x,2\pi n_z/L_z\right)$ with integers $n_x$ and $n_z$, one defines the Fourier series ${\widehat{S}}_{\mathit{k}}$:

$$\begin{array}{c}{\widehat{S}}_{\mathit{k}}(t)=\frac{1}{L_x{L}_z}\int S(x,z,t){\mathrm{e}}^{\mathrm{i}\mathit{k}·\mathit{r}}{\mathrm{d}}^2\mathit{r},\end{array}$$

(29a)

$$\begin{array}{c}S(x,z,t)=\sum _{\mathit{k}}{\widehat{S}}_{\mathit{k}}(t){\mathrm{e}}^{-\mathrm{i}\mathit{k}·\mathit{r}},\end{array}$$

(29b)

where $\mathit{r}=(x,z)$, such that the average energy can be expressed as the sum of contributions in spectral space, as a consequence of Parseval’s theorem:

$$E_S=\frac{1}{2}\sum _{\mathit{k}}{\left|{\widehat{S}}_{\mathit{k}}\right|}^2.$$

(30)

Thus, the energy contained in wavenumber $\mathit{k}$ is ${\widehat{E}}_S(\mathit{k})=\frac{1}{2}{\left|{\widehat{S}}_{\mathit{k}}\right|}^2$.

A spectral energy budget is readily obtained by taking a Fourier transform of Equation (20b), multiplying each Fourier component by the complex conjugate ${\widehat{S}}_{\mathit{k}}^{*}$ and summing over $\mathit{k}$, yielding:

$$\frac{\mathrm{d}}{\mathrm{d}t}{E}_S=\sum _{\mathit{k}}\lambda (\mathit{k}){\left|{\widehat{S}}_{\mathit{k}}\right|}^2,$$

(31)

where λ is the growth rate (21). This energy equation implies that in the saturated state, the energy input from unstable modes compensates energy dissipation by damped modes, with the advection term transferring energy between them.

We define the spatially-averaged vertical salinity flux $F_S$ as:

$$F_S=\frac{1}{L_x{L}_z}\int {\psi }_{x}S\mathrm{d}x\mathrm{d}z.$$

(32)

This flux can be expressed as:

$$F_S=\sum _{\mathit{k}}{\widehat{F}}_S(\mathit{k}),$$

(33)

where:

$${\widehat{F}}_S(\mathit{k})=\frac{\mathrm{i}k}{2}\left({\widehat{\psi }}_{\mathit{k}}{\widehat{S}}_{\mathit{k}}^{*}-{\widehat{\psi }}_{\mathit{k}}^{*}{\widehat{S}}_{\mathit{k}}\right)=-\frac{k^2{\left|K\right|}^2}{k^2+{\left|K\right|}^6}{\left|{\widehat{S}}_{\mathit{k}}\right|}^2,$$

(34)

using the diagnostic relation (20a). Thus, $F_S<0$, as expected.

4.2. Secondary Instability

At some point, the growing elevator modes will trigger secondary instabilities that may be responsible for the break up of the exponentially-growing elevator modes and the saturation of the instability. Figure 4a shows the salinity field obtained from a simulation of the IFSC model (20) in a $2{\ell }_{\mathrm{opt}}\times 5{\ell }_{\mathrm{opt}}$ domain with an initial condition that is a combination of the optimal mode and superposed small amplitude random noise. Here, ${\ell }_{\mathrm{opt}}\equiv 2\pi /k_{\mathrm{opt}}$. The figure shows that, at $t=180$, a secondary instability sets in with a vertical wavenumber comparable to that of the optimal mode. As an indication of what secondary instabilities may be present, we assume a quasi-static approximation where the salt fingers are considered steady. This is a reasonable approximation in the sense that when the secondary instability is important, its growth rate should be larger than that of the salt fingers. We analyze the stability of this state within both the IFSC model (20) and the modified RBC system (18) by means of single vertical mode Floquet theory [5], which considers perturbations of the form:

$$\left(\begin{array}{c}\psi \\ S\end{array}\right)={\mathrm{e}}^{\alpha t+\mathrm{i}k\left(px+qz\right)}\sum _{n=-N}^{n=N}\left(\begin{array}{c}{\psi }_n\\ \mathrm{i}{S}_n\end{array}\right){\mathrm{e}}^{\mathrm{i}nkx}.$$

(35)

Comparisons are made with the full system (6) with $\tau =0.01$, $\mathrm{Ra}=1.1$ and $\mathrm{Sc}=100$, which corresponds to $\mathrm{Pr}=1$ and $R_{\rho }={0.011}^{-1}\approx 90.91$. Given that τ is an asymptotically small parameter, the corresponding mRBC model requires $\mathrm{Ra}=1.1$ and $\mathrm{Sc}=100$, while the corresponding IFSC model only requires $\mathrm{Ra}=1.1$.

We observe in time-stepping simulations of the IFSC model that the secondary instability of the elevator modes is triggered at a salinity amplitude of approximately $S_e\approx 6.5$. This amplitude is selected as the common amplitude for the Floquet theory in all models. For the primitive equations, this implies a renormalization by a factor of ${\mathrm{Ra}}^{-1}$ (see Equation (17b)).

In Figure 4b–d, we show a contour plot of the growth rate α of the secondary instability as a function of the Floquet wavenumbers p and q. The observed match between the three panels demonstrates the validity of the reduced models. Moreover, the prediction of the single vertical mode Floquet theory with $S_{\mathrm{e}}=6.5$ that the largest growth rate occurs for $p=0$ and $q=1.05$ is in excellent agreement with the onset of the observed secondary instability. For comparison, we mention that for $Pr=7$, ${\mathrm{R}}_{\rho }=2$ and $\tau =1/24$, Stern and Simeonov [37] found that the secondary instability sets in with vertical wavenumber $q=0.8$, expressed in units of the optimal wavenumber. These results confirm the validity of the single vertical mode Floquet theory in describing the onset of the secondary instability.

4.3. Towards Saturation

The saturated states of the IFSC model (20) depend on the domain size: when the domain is small, a steady state is reached; in an intermediate-sized domain, chaotic states are observed; for sufficiently large domains, the system reaches statistically steady states involving a broad range of scales. The results reported below are carried out in domains of sufficiently large size that the results are independent of the domain size used. For comparison, we summarize the corresponding small-domain and intermediate-domain results in Appendix A.

In Figure 5, we show the results obtained for $\mathrm{Ra}=1.1$ in a domain of size $32\times 32$ in units of the optimal wavelength ${\ell }_{\mathrm{opt}}=14.86$ with an initial condition in the form of a small amplitude random field. We ran the simulation to $t=4800$, which is large compared to the characteristic time ${\lambda }_{\mathrm{opt}}^{-1}=84.84$ for the growth of the optimal mode, and observed (Figure 5a) that a statistically steady state is finally reached. Three stages toward saturation can be identified: (i) Dominance of salt fingers, $t\lesssim 138$. This is the stage corresponding to peak energy and flux generation. Figure 5b shows the salinity field in the vicinity of this peak $(t\approx 138)$; long finger structures are observed. (ii) Secondary instabilities, $138\lesssim t\lesssim 500$. The salt fingers cannot grow without bound owing to the onset of secondary instability, whose growth rate increases with the amplitude of the growing fingers, as discussed in detail in the previous section. (iii) Saturated state, $t\gtrsim 500$. This regime is characterized by collisions of upward and downward fingers [6]. Figure 5c provides a snapshot of the final statistically steady salinity field. Finger-like structures can still be recognized, but these are of small vertical extent and quite different from those present in Stage (i). It appears that it is a combination of these two effects—secondary instability and salt-finger collision (see Section 5.2)—that is responsible for the observed statistically steady saturated state. In an alternative argument, Radko [38] (see also [8]) suggests that saturation occurs when the growth rates of the linear fingering instability and the secondary instability become comparable. Such a balance implies a balance between linear and nonlinear terms on the same horizontal scale, but as discussed in Section 6, this balance is not always satisfied.

In the next section, we characterize the large-domain saturated state in more detail.

5. Saturated States in the Inertia-Free Salt Convection Model

In large domains (i.e., for large thermal Rayleigh numbers ${\mathrm{Ra}}_T$, cf. Equation (10)), the IFSC model (20) reaches domain size-independent statistically steady states, whose statistical properties depend on the Rayleigh ratio $\mathrm{Ra}$. In this section, we study this $\mathrm{Ra}$ dependence, since this is crucial for a sound parameterization of small-scale salt-fingering convection and its effect on large-scale fields. Global quantities—available potential energy and flux—are studied in Section 5.1, where we identify two distinct regimes according to their dependence on $\mathrm{Ra}$. The corresponding probability density functions are shown in Section 5.2.

Figure 6 shows the spectral density of kinetic energy (first column), salinity potential energy (second column) and salinity flux (third column) in the statistically stationary state as a function of the horizontal wavenumber k and the vertical wavenumber m (top two rows). The third row shows the corresponding integrated spectra along the vertical (plain lines) and horizontal (dashed lines) directions, plotted as a function of the horizontal and vertical wavenumber, respectively. One distinguishing feature of the vertically-integrated spectra is that they peak around the optimal wavenumber $k_{\mathrm{opt}}$. Away from the peak, the vertically-integrated spectra increase weakly for wavenumbers below $k_{\mathrm{opt}}$, but fall off rapidly for wavenumbers larger than this energy injection wavenumber. Similar tendencies are also observed in the horizontally-integrated spectra. In both cases this fall-off is more rapid for smaller $\mathrm{Ra}$ ($\mathrm{Ra}=1.1$) than for larger $\mathrm{Ra}$ ($\mathrm{Ra}=5$). Of particular interest is the fact that energy is transferred to ever larger scales and in fact appears to peak at these scales. We surmise that this is a consequence of the patchiness of the saturated state (Figure 5c), a property that is also observed in three-dimensional simulations of the primitive equations [11]. Figure 7 compares the salinity fields in the stationary state for $\mathrm{Ra}=1.1$ (Figure 7a) and $\mathrm{Ra}=5$ (Figure 7b). The fingers are clearly visible, as are collisions between descending salt fingers and rising fresh-water fingers. Large-scale patchy structures with fingers of finite length are also observed, indicating the presence of multiscale dynamics that is characteristic of Regime II; cf. Section 5.1.2.

5.1. Regimes

We consider domains of size $32{\ell }_{\mathrm{opt}}\times 32{\ell }_{\mathrm{opt}}$ where ${\ell }_{\mathrm{opt}}\equiv 2\pi /k_{\mathrm{opt}}$ is the optimal wavelength and depends on the $\mathrm{Ra}$ value used. In all of the simulations, the initial conditions are taken to be a small amplitude random field, and the model equations are integrated for a sufficiently long time that a statistically stationary state is reached. All averages are calculated after discarding an initial transient.

In Figure 8, we show the dependence of the energy $E_S\equiv (1/2)\int S^2\mathrm{d}x\mathrm{d}z/\int \mathrm{d}x\mathrm{d}z$ and the salinity flux $F_S\equiv -\int {\psi }_{x}S\mathrm{d}x\mathrm{d}z/\int \mathrm{d}x\mathrm{d}z$ on the parameter $\mathrm{Ra}$, ranging from $1.02$–10. In this figure, three intervals of power-law dependence on the supercriticality $\mathcal{R}\equiv \mathrm{Ra}-1$ are seen. The first of these corresponds to Regime I, while the latter two correspond to different sub-limits of Regime II.

The dynamics of the IFSC model (20) are driven by the linear salt-finger instability, which has the distinctive feature that the optimal (horizontal) wavenumber is finite because of the coexistence of large- and small-scale damping stemming from the stabilizing temperature and viscosity, respectively. In the following, we therefore assume that the characteristic horizontal scale of the system is determined by the optimal wavenumber $k_{\mathrm{opt}}$ and confirm this assumption in Figure 9, which shows the correspondence between the peak in the energy spectrum obtained from numerical simulations (Figure 6) and expression (24) for the optimal wavenumber. Here, the constant ratio ${\mathrm{e}}^{0.15}\approx 1.16$ between $k_{\mathrm{finger}}$ and $k_{\mathrm{opt}}$ confirms the validity of our assumption.

In Figure 10, we compare the dominant balances in the prognostic equation with $\mathrm{Ra}=1.02$ and $1.7$. Here, $NL=\mathcal{J}\left(\psi ,S\right)$ and $LI=\mathrm{Ra}{\psi }_x-{\nabla }^{2}S$ denote nonlinear advection and the linear part of the prognostic equation. In the left panel (Regime I), we find a balance between the time derivative, nonlinear advection and linear instability, while the right panel shows that Regime II is characterized by a balance between the time derivative and nonlinear advection only, with the small difference between these two terms balancing the remaining linear term $LI$. In Figure 11, we compare the corresponding dominant balances in the diagnostic relation at the same parameter values, $\mathrm{Ra}=1.02$ and $1.7$. In the left panel (Regime I), we find a dominant balance between ${\partial }_x^2\psi$ and ${\nabla }^2{\partial }_{x}S$, while in the right panel, ${\nabla }^6\psi$ becomes comparable to or larger than ${\partial }_x^2\psi$ and so participates fully in the dominant balance.

The different dominant balances in the prognostic and diagnostic equations characterizing Regimes I and II indicate that these regimes are fundamentally different. Specifically, consideration of the regime implications on the inertial-free diagnostic balance indicates that Regime I is dominated by a balance between the thermals and salinity buoyancy forces on the characteristic spatial scale of the fingers; recall that in Regime I, $k_{\mathrm{opt}}\ll 1$ (Figure 3). In this regime viscous dissipation is only important at scales smaller than the finger scale and takes place primarily at the interface between fingers. Viscous dissipation gains importance in the inertia-free balance in Regime II. This occurs as a direct consequence of: (i) the reduced separation between the fingering scale $2\pi /k_{\mathrm{opt}}$∼1 and the viscous dissipation scale; and (ii) the increased collisions between rising and descending fingers. In the following, we include in the term collision any process that leads to a local inversion of the salinity gradient, such as might occur when a descending finger develops an overdense bulbous head.

Based on the expression (24) for the optimal wavenumber and the dominant balances identified above, we can calculate the parameter dependences in both regimes.

5.1.1. Regime I

Regime I is a small supercriticality regime with $\mathcal{R}=\mathrm{Ra}-1\ll 1$, where the optimal scale is of $O({\mathcal{R}}^{-1/4})$ (cf. Equation (26a)) and, hence, large. Here, by large scale we mean that the dominant balance in the diagnostic equation, Equation (20a), is:

$${\partial }_x^2\psi \sim {\nabla }^2{\partial }_{x}S,$$

(36)

a conclusion confirmed in the left panel of Figure 11. This statement implies an inviscid inertia-free balance between the thermal and salinity buoyancy forces. In this regime, the appropriate balance in the prognostic equation, Equation (20b), is between the nonlinear and linear terms (Figure 10):

$$\mathcal{J}\left(\psi ,S\right)\sim -\mathrm{Ra}{\psi }_x+{\nabla }^{2}S\sim \lambda S,$$

(37)

where the last relation comes from the linear stability problem. Here, λ ∼ ${\mathcal{R}}^{3/2}$ as obtained from Equation (26b) and confirmed in Figure 3a.

We assume power-law dependence on $\mathcal{R}$ of the field magnitudes:

$$\psi \sim {\mathcal{R}}^a\mathrm{and}S\sim {\mathcal{R}}^b.$$

(38)

With the assumption that the order of magnitude of both the horizontal and vertical derivatives is determined by the optimal wavenumber, i.e.,

$${\partial }_x\sim {\partial }_z\sim {\mathcal{R}}^{1/4}$$

(39)

(Equation (26a)) we substitute (38) into (36) and (37) to obtain:

$$\begin{array}{cc}\hfill 1/2+a& =3/4+b,\hfill \end{array}$$

(40a)

$$\begin{array}{cc}\hfill 1/2+a+b& =3/2+b,\hfill \end{array}$$

(40b)

implying that

$$a=1\mathrm{and}b=3/4.$$

(41)

These exponents are used to generate the Regime I black line in both panels of Figure 8. We observe a good match between the theoretical and numerical results.

5.1.2. Regime II

In Regime II, the supercriticality $\mathcal{R}$ becomes larger than that in Regime I, and therefore, in the diagnostic equation, Equation (20a), a full balance obtains (cf. Figure 11):

$$({\partial }_x^2+{\nabla }^6)\psi \sim {\nabla }^2{\partial }_{x}S.$$

(42)

However, in contrast to Regime I, the leading order prognostic balance now involves the time derivative and the nonlinear terms only, leaving the linear term subdominant. In this regime, we suppose that energy is transferred from the small energy injection scale, which peaks at the optimal wavelength of the fingers, to large scales, a scenario inspired by the observed increase in energy with increasing scale (Figure 6) and the multiscale operator present in the diagnostic equation. This balance in energy transfer follows the same phenomenology as the energy cascade in isotropic turbulence, where the dissipative effect is also subdominant in the prognostic momentum equation, but determines the energy flux. However, our system differs from isotropic turbulence by the presence of two distinct dominant scales in the energy transfer that participate in the nonlocal interactions in spectral space.

The large-scale field must satisfy the diagnostic relation with negligible small-scale damping:

$${\partial }_x^2{\psi }^{(l)}\sim {\nabla }^2{\partial }_x{S}^{(l)},$$

(43)

where the superscript $(l)$ indicates the large-scale field. As to the prognostic balance, we assume that the large-scale field reaches a distinguished regime, such that the nonlinear and linear terms balance:

$$\mathcal{J}\left({\psi }^{(l)},S^{(l)}\right)\sim -\mathrm{Ra}{\psi }_x^{(l)}+{\nabla }^2{S}^{(l)}\approx \mathcal{R}{\nabla }^2{S}^{(l)}.$$

(44)

Here, the last relation takes advantage of the large scale of the fields and the definition of the supercriticality $\mathcal{R}$.

As in Regime I, we continue to assume that the scales of the fields are controlled by the linear theory optimal mode. However, since $k_{\mathrm{opt}}$ is no longer small, we need to keep the full expression for the optimal wavenumber in Regime II:

$${\partial }_x\sim {\partial }_z\sim k_{\mathrm{opt}}.$$

(45)

Using (45), we can now calculate the $\mathrm{Ra}$ dependence of the large-scale stream function and salinity from (43) and (44):

$${\psi }^{(l)}\sim \mathcal{R}\mathrm{and}{S}^{(l)}\sim \frac{\mathcal{R}}{k_{\mathrm{opt}}}.$$

(46)

We express the above-mentioned energy transfer picture in terms of a balance between the advection of potential energy ($S^2$) through small- to large-scale interaction and small-scale energy pumping:

$$\mathcal{J}\left({\psi }^{(\alpha )},S^{(\beta )}{S}^{(\gamma )}\right)\sim \mathrm{Ra}{\psi }_{x}S,$$

(47)

where $\mathcal{J}\left({\psi }^{(\alpha )},S^{(\beta )}{S}^{(\gamma )}\right)$ is the larger of $\mathcal{J}\left({\psi }^{(l)},SS\right)$, $\mathcal{J}\left(\psi ,S^{(l)}S\right)$, $\mathcal{J}\left({\psi }^{(l)},S^{(l)}S\right)$ and $\mathcal{J}\left(\psi ,S^{(l)}{S}^{(l)}\right)$, and the quantities without the superscript $(l)$ denote small-scale fields. We need to calculate the four possibilities and check that the computed exponentials ensure the correct small- to large-scale interaction, i.e., that the chosen small- to large-scale advection term used to balance the energy input is the largest of the four. After going through these four cases, we conclude that the right term is $\mathcal{J}\left(\psi ,S^{(l)}{S}^{(l)}\right)$. This choice of interaction term is confirmed in Figure 6: the salinity spectrum peaks at both large and small scales, while the flux spectrum only peaks at the small scale, implying that advection is dominated by the small-scale velocity.

The resulting balance:

$$\mathcal{J}\left(\psi ,S^{(l)}{S}^{(l)}\right)\sim \mathrm{Ra}{\psi }_{x}S$$

(48)

gives us the $\mathrm{Ra}$ dependence of S:

$$S\sim \frac{k_{\mathrm{opt}}}{\mathrm{Ra}}{\left(S^{(l)}\right)}^2\sim \frac{{\mathcal{R}}^2}{k_{\mathrm{opt}}\mathrm{Ra}}.$$

(49)

Substituting (49) into (42) now yields:

$$\psi \sim \frac{k_{\mathrm{opt}}}{1+C{k}_{\mathrm{opt}}^4}S\sim \frac{{\mathcal{R}}^2}{(1+C{k}_{\mathrm{opt}}^4)\mathrm{Ra}},$$

(50)

where the quantity C measures the degree of isotropy. Empirically, we find that $C=2$ is a good match, for reasons that remain to be explored. Expressions (49) and (50) are used in Figure 8 to generate the red curves.

Finally, making use of (49) and (50), we can obtain the scaling in the two subcases of Regime II: when $\mathcal{R}$ is comparable to unity, $k_{\mathrm{opt}}\sim {\mathcal{R}}^{1/4}$ and $\mathrm{Ra}=1+\mathcal{R}\sim {\mathcal{R}}^{1/4}$ (cf. Appendix B), and expressions (49) and (50) can be approximated by a power-law dependence on $\mathcal{R}$:

$$\psi \sim S\sim {\mathcal{R}}^{3/2}.$$

(51)

When $\mathrm{Ra}$ is large, Equations (49) and (50) become:

$$\psi \sim S\sim {\mathcal{R}}^1.$$

(52)

We refer to these two regimes as II${}_1$ and II${}_2$, respectively. Expressions (51) and (52) are used to generate the four black lines in Figure 8 when $\mathcal{R}$ is not small.

We summarize the scaling results for Regimes I and II in

Table 1. Exponents of $\mathcal{R}\equiv \mathrm{Ra}-1$ in the power-law scalings of S, ψ and the salinity flux $F_S$ in Regimes I and Sub-regimes II${}_1$ and II${}_2$ for small and large $\mathcal{R}$, respectively.

Regime	S	ψ	${\mathit{F}}_{\mathit{S}}$
Regime I	$3/4$	1	2
Regime II${}_1$	$3/2$	$3/2$	$13/4$
Regime II${}_2$	1	1	2

.

5.2. Probability Density Functions

To obtain a detailed understanding of the statistics of the saturated fields, we now turn to the study of the probability density functions (pdfs) of the different fields. In Figure 12 and Figure 13, we show the pdfs of the salinity and its derivatives and the pdfs of the stream function and its derivatives in both Regime I ($\mathrm{Ra}=1.02$) and Regime II ($\mathrm{Ra}=1.7$). For $\mathrm{Ra}=1.02$, the pdfs are obtained from data over $t\in$(12,000, 13,000) with a step size $\Delta t=100$, while for $\mathrm{Ra}=1.7$, we show one snapshot at $t=791$. All pdfs presented are normalized by their variance.

We observe that most fields obey Gaussian or quasi-Gaussian statistics. For more detailed assessment, we fit the pdfs to a general stretched-Gaussian form:

$$P(a)=P_{\max }{exp}^{-{\left|a/\beta \right|}^{\alpha }},$$

(53)

where α, β and $P_{\max }$ are constants and $\alpha =2$ corresponds to the Gaussian distribution. The parameter β is adjusted to obtain a probability density function with variance of one, and $P_{\max }$ normalizes the distribution:

$$P_{\max }=\frac{\alpha }{2\beta \Gamma (1/\alpha )}\mathrm{and}{\beta }^2=\frac{\Gamma (1/\alpha )}{\Gamma (3/\alpha )},$$

(54)

where Γ is the gamma function. From the second relation above, we calculate β in

Table 2. Curve-fitting parameters for the probability density functions in Figure 12 and Figure 13. Here the $“+”$ and $“-”$ symbols denote the positive and negative arguments of the pdfs, respectively.

Field	$\mathit{Ra}=1.02$		$\mathit{Ra}=1.7$
	α	β	α	β
S	1.8	1.33	1.8	1.33
$S_x$	2	1.41	2	1.41
$S_z^{+}$	1.5	1.43	0.75	0.5
$S_z^{-}$	1.7	1	1	0.45
ψ	2	1.41	2	1.41
w	2	1.41	2	1.41
u	2	1.41	2	1.41

for the pdfs of S, $S_x$, ψ, w and u. This form has been utilized in the context of Rayleigh–Bénard convection in [39].

In addition to the above general form, there are several other interesting features that can be learned from these figures. In the right panel of Figure 12, the pdfs of $S_z$ are not symmetric with respect to their peaks, indicating obvious skewness. This asymmetry originates from the collisions between positive (rising) and negative (descending) fingers: collisions occur where $S_z$ is positive and tend to increase the salinity gradient. Nontrivial skewness is a common feature of convection where coherent structures, such as fingers or plumes, exist and has been studied in Rayleigh–Bénard convection to identify plume generation [40]. In the same way, enhanced positive tails in the pdfs of $S_z$ for larger $\mathrm{Ra}$ indicate the increasing importance of finger collision.

However, there is no such asymmetry in ψ because ${\psi }_z$ is the horizontal velocity, which preserves on average the reflection symmetry $u\to -u$ of the system. The pdfs of $S_x$, ψ, w and u are all Gaussian (cf. Figure 12 and Figure 13), which is very useful for parameterization. Since Section 5.1 develops a theory for the dependence of the variance of the different fields on the parameter $\mathcal{R}$ (equivalently $\mathrm{Ra}$), the statistics of the quantities with Gaussian distribution are fully understood, and these can be used to parameterize the diffusivity induced by salt-fingering convection as $\mathrm{Ra}$ varies.

6. Discussion and Conclusions

In this paper, we systematically derived two reduced models for salt-fingering convection in the limit of small diffusivity and large density ratios, $\tau \ll 1$ and $R_{\rho }=O({\tau }^{-1})$. In the modified RBC model (18), the stabilizing temperature plays the role of large-scale damping, which when combined with the small-scale viscosity leads to a finite optimal wavenumber corresponding to the fastest linear growth rate. This model may find application in the astrophysical context where $\mathrm{Sc}=O(1)$. In the parameter regime corresponding to the oceans where $\mathrm{Sc}$ is large, we obtained a simpler model, the inertia-free salt convection model (20), which has a prognostic-diagnostic form, and we studied this model in detail. In large domains, corresponding to large thermal and solutal Rayleigh numbers, we checked that this model preserves the linear and secondary instabilities present in the primitive equations and that it captures the three stages—finger dominance, finger disruption and saturation—that have been observed to lead to a statistically steady saturated state. Both secondary instability and finger collision were found to be responsible for the presence of this state, thereby extending the conclusion reached by Shen [6] that finger collisions are important at $O(1)$ density ratios to the limit of large density ratios.

The properties of the saturated state depend on the Rayleigh ratio $\mathrm{Ra}$. We identified two parameter regimes characterized by different balances in both the diagnostic and prognostic equations, corresponding to weak and strong supercriticalities, leading to distinct saturated states. In Regime I, where the width of the salt fingers scales as $k_{\mathrm{opt}}^{-1}$∼${\mathcal{R}}^{-1/4}\ll 1$, the temperature and salinity buoyancy effects dominate the diagnostic equation, and viscous dissipation is negligible. In the prognostic equation, a balance between linear instability and nonlinear advection indicates that the growth rates of primary and secondary instabilities are of the same order (cf. [10,38]). In this regime, we find power-law dependences on the supercriticality of the saturated state. Regime II occurs at larger supercriticality. Even though this indicates stronger driving forces, viscous dissipation enters the diagnostic balance. This counterintuitive fact results from the smaller salt-finger scale at larger supercriticalities. In addition to the small salt-finger scale, a large scale is also present. This scale is fed by an energy transfer from the small finger scale at which energy is injected into the system, and it is this energy transfer that characterizes Regime II. The result relies on the dependence of the energy injection scale on the Rayleigh ratio, leading to two sub-limits with distinct power-law dependence.

The presence of a saturated state even in the presence of exact elevator mode solutions that grow exponentially forever requires that our simulations are initialized with broad-band initial conditions, cf. [41,42]. In this case, the system is unable to ‘find’ the growing elevator mode state, resulting in a saturated state exhibiting a broad range of spatial scales. Figure 5a shows an example of the transient generation of elevator modes followed by their disruption via secondary instability. The resulting statistically steady state is characterized by pdfs that have in general a stretched Gaussian form. In particular, the horizontal and vertical velocities have Gaussian distributions and can therefore be used to calculate the diffusivity of turbulent salt-finger convection, which can be further developed into a parameterization of transport processes. In this connection, we mention that in two-dimensional Rayleigh–Bénard convection, such Gaussian pdfs are universal [43]. Moreover, passive tracers in such a system have identical statistics [43]. In salt-finger convection, it is the salinity field that drives convection, but a tracer like the temperature is not passive since it contributes to buoyancy. As a result, one might expect departures from Gaussian pdfs, and indeed, such departures have been detected [44]. However, in both of our reduced models, the temperature field is slaved to the stream function, and the resulting models resemble the equations governing two-dimensional Rayleigh–Bénard convection. We conjecture that this is the reason why we find Gaussian velocity pdfs in our computations.

Our IFSC model (20) is linked to the weakly nonlinear model derived by Radko [7] (see Equations (8) and (11) of [7]). Even though different parameter regimes are considered in the two models—our reduced model for $\tau \to 0$ captures a fully-nonlinear regime where dynamics are not confined to the onset of instability ($\mathrm{Ra}$ ranges from 1–∞), while Radko’s [7] model relaxes the small τ assumption, but is restricted to dynamics near the onset: the two models match in the overlapping regime $\tau \ll 1$ and $\mathrm{Ra}=1+ϵR$ with $ϵ\ll 1$ and $R=O(1)$. In this regime, the optimal scale from Equation (24) is of $O({ϵ}^{-1/4})$, suggesting the rescaling:

$${\partial }_t\to {ϵ}^{9/4}{\partial }_t,\nabla \to {ϵ}^{1/4}\nabla ,\psi \to ϵ\psi \mathrm{and}S\to {ϵ}^{3/4}S.$$

(55)

With this rescaling, the leading order contribution to (20) reads:

$${\partial }_{t}S+\mathcal{J}\left({\partial }_x^{-1}{\nabla }^{2}S,S\right)+\left(R-{\partial }_x^{-2}{\nabla }^6\right){\nabla }^{2}S=0,$$

(56)

which is identical to the model derived by Radko [7] in the limit of small τ.

In our simulations, we found that the characteristic length of the fingers increases as $\mathrm{Ra}$ increases (cf. Figure 5 and Figure 7), in contrast to a result of Merryfield and Grinder (cf. Figure 8 in [45]). Evidently, salt-finger convection exhibits different limiting parameter regimes: in our model, τ and ${\mathrm{R}}_{\rho }^{-1}$ are constrained to be of the same order and to be small, while in the simulations of Merryfield and Grinder, the value of τ is fixed.

We may define a salinity Nusselt number by:

$${\mathrm{Nu}}_S=-\frac{F_S}{\mathrm{Ra}},$$

(57)

and use the Regime II scalings to obtain ${\mathrm{Nu}}_S\sim \mathrm{Ra}$ in the limit of large $\mathrm{Ra}$. This result is to be compared with the asymptotic result ${\mathrm{Nu}}_S\sim {\mathrm{Ra}}^{1/3}$ obtained by Yang et al. [23] for a vertically-bounded domain. Evidently, the presence of elevator modes permitted in the doubly-periodic domain leads to a very substantial enhancement of the salinity flux. In view of (10), our thermal Rayleigh number ${\mathrm{Ra}}_T$ is of order ${10}^7$, while the salinity Rayleigh number ${\mathrm{Ra}}_S$ is of order ${10}^8$ (see Equation (15)). These values are comparable to those used in [23], although our density ratio is much larger than that in [23].

Salt-finger convection often results in the formation of salinity staircases [11,38,46,47], and these have been extensively studied in oceanographic measurements, numerical simulations and theories. In our IFSC model, we have not observed the formation of staircases. This may be because our model filters out gravity waves, which are believed to be important for staircase formation through collective instability [5,9]. The γ-instability mechanism proposed by Radko [38] provides an alternative explanation, but requires a nonmonotonic dependence of the flux ratio on the mean density gradient. In our reduced model, the background temperature gradient remains at leading order constant. Consequently, within the model, this condition can be interpreted as requiring a nonmonotonic dependence of the salt flux on the mean salinity gradient, a situation that never arises. Thus, neither mechanism for staircase formation is present. Having said this, Brown et al. [10] and Garaud et al. [48] recently discovered that staircases can form in three dimensions even when the flux increases monotonically with the density gradient, although their mechanism also relies on gravity wave generation. In summary, the possibility that the model (in two or three dimensions) captures staircase formation must be kept open, particularly at large $\mathrm{Ra}$, since oceanic staircases are observed in regions with density ratios that are in fact quite small.

In a future publication, we plan to study the large $\mathrm{Ra}$ regime of the model and to discuss the extent to which the conclusions of this paper carry over to three dimensions.