A reduced model for salt-finger convection in the small diffusivity ratio limit

A simple model of nonlinear salt-finger convection in two dimensions is derived and studied. The model is valid in the limit of small solute to heat diffusivity ratio and 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\'enard 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 the energy input through linear instability and the 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 solute and heat Rayleigh numbers. The spectra and probability density functions are also computed.


I. 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 : for example, heat and salt contribute to the stratification in the oceanic setting, heat and chemical composition in the 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 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 τ < 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 τ 1.
Linear stability properties of the salt-finger regime were first determined by Stern 3 . As the magnitude of the unstable mode grows secondary instabilities are triggered 4 , and finally lead to a statistically steady state 5,6 . However, the physical processes behind the presence of such an equilibrium state are still incompletely understood. Stern 7 has suggested that a a) Electronic address: j.h.xie@berkeley.edu b) Electronic address: benjamin.miquel@colorado.edu c) Electronic address: julien@colorado.edu d) Electronic address: knobloch@physics.berkeley.edu 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 4 , becomes comparable to that of the salt-finger instability itself 8,9 . More recently Shen 5 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 10 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, eg. 6,11 , 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 proved invaluable and identified a number of novel processes, eg. 12 . To shed light on some of these we propose here a systematic procedure that leads to a simplified set of equations that are easier to study, both theoretically and numerically. In the following we refer to these, following 13 , as reduced equations. The procedure described below has been used before 14,15,16 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 17 and in some oceanic measurements 18 , and is a limit beyond the current capability of direct numerical simulations. Two reduced models are derived depending on the order of Schmidt number Sc, which captures the relative strength of momentum and salinity diffusion. When these two effects are comparable, Sc = O(1), a reduced model that bears similarity to the Rayleigh-Bnard configuration is derived. This model is applicable to astrophysical objects, eg. 9,16,19 . In the oceanic parameter regime with fast momentum diffusion, Sc 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 15 . 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, eg. 20,21,22,23,24 . In this formulation elevator modes consisting of vertically invariant spatial dynamics 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.
The structure of this paper is as follows. In §II we formulate the salt-finger convection problem followed in §III 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 §IV and §V through a combination of numerical simulation and scaling analysis. Our results are summarized in §VI. 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.

A. 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 β T > 0 and a destabilizing background salinity gradient β S > 0. We introduce the thermal and saline expansion coefficients α T > 0 and α S > 0 so that the density of the fluid ρ(T, S) at temperature T and salinity S is given by The remaining material properties of the fluid are specified by the thermal diffusivity κ T , the solute diffusivity κ S and the viscosity ν, all of which are taken to be constant. The acceleration due to gravity is denoted by g.
The equations are made dimensionless using the natural salt finger width d given by This length defines the characteristic temperature T and salinity S scales: These scales are used to nondimensionalize the temperature and salinity profiles T total and S total which are decomposed into the initial background state and a fluctuating component indicated by a tilde: The incompressible velocity field is expressed in terms of a streamfunction ψ, where (u, w) and ψ are made dimensionless using κ S /d and κ S , respectively. Finally, time is made dimensionless using the salinity diffusion timescale:

B. Nondimensional equations
The governing dimensionless equations for fluctuations around the conduction state ψ = Three dimensionless ratios enter this set of equations: the Prandtl number Pr, the inverse Lewis number τ and the density ratio R ρ specifying of the relative contribution of the temperature and salinity to the background density gradient. A fourth dimensionless ratio, the Schmidt number Sc, can be introduced as an alternative to the diffusivity ratio. These quantities are defined by In this paper, we adopt doubly periodic boundary conditions for ψ,T andS, and consider a statically stable configuration with R ρ > 1.

numbers
The above formulation can be connected to the classical description of doubly diffusive convection in terms of the temperature and salinity Rayleigh numbers, Ra T and Ra S , respectively. These are defined for a layer of vertical extent H, viz., implying that Large 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

D. 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 ψ(x, z, t),T (x, z, t),S(x, z, t) = (ψ e , T e , S e )e λt+i(kx+mz) . The growth rate λ obeys the dispersion relation where |K| 2 ≡ k 2 + m 2 . For a Rayleigh-Taylor stable layer such that R ρ > 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 The salt-finger instability is therefore present in the regime 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 d.

A. Scalings of parameters and variables
In this section we focus on the small diffusivity ratio limit τ 1 of the system (6), i.e.
We refer to this system, considered in §III.C, as the inertia-free salt convection model.
In all cases we suppose, moreover, that the density ratio R ρ is large and comparable to τ −1 . This assumption is supported by both laboratory experiments in the fingering domain 25 and oceanic measurements 18 where R ρ can reach O(10 2 ). Also, in the astrophysics context, e.g. for stars along the red giant branch (RGB), R ρ can be as large as O(10 6 ), which is of the same order as the corresponding Prandtl number 26 . Hence, we define the order one parameter: This parameter is analogous to a Rayleigh number for our system, as we demonstrate below.
The parameter R, 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 streamfunction which remains untouched. The temperature fluctuationsT are scaled with τ and become an asymptotically small quantity in the limit τ ↓ 0. In this limit the salinity fluctuationsS remain finite but are scaled for convenience by the order one factor Ra: 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.

B. Modified Rayleigh-Bénard convection model: Sc = O(1)
In the limit of small τ and Sc = O(1), with the scalingT = τ T , we obtain a distinguished regime described by the following set of equations: 1 Sc Once ψ is determined the scaled temperature T is obtained from Equations (18)  In this regime, characterized by condition (14b), viscosity is larger 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 Eqs. (III.B) to a model which is now first order in time: In the following we refer to these equations as the inertia-free salt convection (IFSC) model.
As before, the temperature T is determined from equation (19). 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.

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.

A. Linear stability
Substituting a normal mode ansatz of the form e λt+i(kx+mz) into Eqs. (20) linearized around S = ψ = 0 we obtain the growth rate λ as a function of Ra: where, as before, |K| 2 = k 2 + m 2 . It follows that the threshold for instability, determined in Eq. (12), is preserved in the IFSC model and occurs at Thus the quantity R ≡ Ra−1 measures the distance from threshold, i.e., the supercriticality of the system.
We show the linear growth rate λ(k, m) in  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 R: As already mentioned, these modes are exact solutions of the nonlinear equations (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 Fig. 2 for several values of Ra.
In the absence of boundaries the fastest growing mode is an elevator mode such that  with growth rate given by In Fig. 3 we replot these results in a log-log plot, indicating the asymptotic limits valid for R 1 and R → ∞: When stress-free, fixed temperature and fixed salinity boundaries conditions are applied on horizontal boundaries separated vertically by a gap 2π, the marginal stability curve corresponds to m opt = 1 and yields a result highly reminiscent of the marginal stability curve of the classic Bénard problem with similar boundary conditions: 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: For k ∈ (2πn x /L x , 2πn z /L z ) with integers n x and n z , one defines the Fourier seriesŜ k : where 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: 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 This flux can be expressed as whereF using the diagnostic relation (11). Thus F S < 0, as expected.

B. 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 4(a) shows the salinity field obtained from a simulation of the IFSC model (20) in a 2 opt × 5 opt domain with an initial condition that is a combination of the optimal mode and superposed small amplitude random noise. Here opt ≡ 2π/k 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 4 , which considers perturbations of the form Comparisons are made with the full system (6) (6), (c) the modified RBC equations (18) and (d) the IFSC model (20). Wavenumbers are normalized by the optimal wavenumber k opt given by (24).
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 Ra (see Eq. (17)).
In Fig. 4 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 §VI this balance is not always satisfied.
In the next section we characterize the large-domain saturated state in more detail.

MODEL
In large domains (i.e., for large thermal Rayleigh numbers Ra T , cf. Eq. (10)), the IFSC model (20) reaches domain size-independent statistically steady states, whose statistical properties depend on the Rayleigh ratio Ra. In this section we study this 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 §V.A, where we identify two distinct regimes according to their dependence on Ra. The corresponding probability density functions are shown in §V.B. 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 vertical and horizontal wavenumber, respectively. One distinguishing feature of the vertically integrated spectra is that they peak around the optimal wavenumber k opt . Away from the peak, the vertically integrated spectra increase weakly for wavenumbers below k opt , but fall off rapidly for wavenumbers larger than this energy injection wavenumber. Similar tendencies are also observed in the horizontally integrated spectra. At smaller Ra (Ra = 1.1) this fall-off is more rapid with increasing horizontal wavenumber than with increasing vertical wavenumber, a situation that reverses for larger Ra (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 (Fig. 5(c)), a property that is also observed in three-dimensional simulations of the primitive equations 10 . Figure   7 compares the salinity fields in the stationary state for Ra = 1.1 (Fig. 7(a)) and Ra = 5 ( Fig. 7(b)). 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. §V.A.2.

A. Regimes
We consider domains of size 32 opt ×32 opt where opt ≡ 2π/k opt is the optimal wavelength and depends on the Ra value used. In all 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 Fig. 8 we show the dependence of the energy E S ≡ (1/2) S 2 dxdz/ dxdz and the salinity flux F S ≡ − ψ x Sdxdz/ dxdz on the parameter Ra, ranging from 1.02 to 10. In this figure, three intervals of power-law dependence on the supercriticality R ≡ Ra − 1 are seen. The first of these corresponds to Regime I, while the latter two correspond to different sublimits 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 opt , and confirm this assumption in Fig. 9, which shows the correspondence between the peak in the energy spectrum obtained from numerical simulations (Fig. 6) and expression (24) for the optimal wavenumber. Here the constant ratio e 0.15 ≈ 1.16 between k finger and k opt confirms the validity of our assumption.
In Fig. 10  I) we find a dominant balance between ∂ 2 x ψ and ∇ 2 ∂ x S, while in the right panel ∇ 6 ψ becomes comparable to or larger than ∂ 2 x ψ 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 opt 1 (Fig. 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π/k opt ∼ 1 and the viscous dissipation scale, and (ii) the increased collisions between rising and descending fingers.
Based on the expression (24) for the optimal wavenumber and the dominant balances identified above, we can calculate the parameter dependences in both regimes.

Regime I
Regime I is a small supercriticality regime with R = Ra − 1 1, where the optimal scale is of O(R −1/4 ) (cf. Eq. (24)) and hence large. Here, by large scale we mean that the dominant balance in the diagnostic equation (20a) is a conclusion confirmed in the left panel of Fig. 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 (20b) is between the nonlinear and linear terms ( Fig. 10): where the last relation comes from the linear stability problem. Here λ ∼ R 3/2 as obtained from Eq. (25) and confirmed in Fig. 3.
We assume power-law dependence on R of the field magnitudes: With the assumption that the order of magnitude of both the horizontal and vertical derivatives is determined by the optimal wavenumber, i.e., (see Eq. (26a)) we substitute (38) into (36) and (37) to obtain implying that a = 1 and b = 3/4.
These exponents are used to generate the Regime I black line in both panels of Fig. 8. We observe a good match between the theoretical and numerical results.

Regime II
In Regime II, the supercriticality R becomes larger than that in Regime I and therefore in the diagnostic equation (20a) a full balance obtains (cf. Fig. 10): 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 (Fig. 6)  The large-scale field must satisfy the diagnostic relation with negligible small-scale damping: 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 term balance: Here the last relation takes advantage of the large scale of the fields and the definition of the supercriticality 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 opt is no longer small we need to keep the full expression for the optimal wavenumber in Regime II: Using (45) we can now calculate the Ra-dependence of the large-scale streamfunction and salinity from (43) and (44): 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: where J ψ (α) , S (β) S (γ) is the larger of J ψ (l) , SS , J ψ, S (l) S , J ψ (l) , S (l) S and J ψ, S (l) S (l) , 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 J ψ, S (l) S (l) . This choice of interaction term is confirmed in Fig. 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 gives us the Ra dependence of S: Substituting (49) into (42) now yields 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 Fig. 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 R is comparable to unity, k opt ∼ R 1/4 and Ra = 1+R ∼ R 1/4 (cf. Appendix B), and expressions (49) and (50) can be approximated by a power-law dependence on R: When Ra is large, Eqs. (49) and (50) become 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 Fig. 8 when R is not small.
We summarize the scaling results for regimes I and II in Table I.

B. 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 Fig. 12 and   13 we show the pdfs of the salinity and its derivatives and the pdfs of streamfunction and its derivatives of both Regime I (Ra = 1.02) and Regime II (Ra = 1.7). For Ra = 1.02, the pdfs are obtained from data over t ∈ (12000, 13000) with step size ∆t = 100, while for 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: where α, β and P max are constants and α = 2 corresponds to the Gaussian distribution. The parameter β is adjusted to obtain a probability density function with variance 1, and P max normalizes the distribution: where Γ is the gamma function. From the second relation above we calculate β in Table II for the pdfs of S, S x , ψ, w and u. This form has been utilized in the context of Rayleigh-Bénard convection in 33 .
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 Fig. 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 34 . However, there is no such asymmetry in ψ because ψ z is the horizontal velocity, which preserves on average the reflection symmetry u → −u of the system. The pdfs of S x , ψ, w and u are all Gaussian (c.f. Fig. 12 and 13), which is very useful for parameterization. Since §V.A develops a theory for the dependence of the variance of the different fields on the parameter R (equivalently Ra), the statistics of the quantities with Gaussian distribution are fully understood, and these can be used parameterize the diffusivity induced by salt-fingering convection as Ra varies.  Table II.  Table II.

VI. DISCUSSION
In this paper, we systematically derived two reduced models for salt-fingering convection in the limit of small diffusivity and large density ratios, τ 1 and R ρ = O(τ −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 Sc = O(1). In the parameter regime corresponding to the oceans where 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, 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. 35,36 . 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 5(a) 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 37 . Moreover, passive tracers in such a system have identical statistics 37 . 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 38 . However, in both our reduced models the temperature field is slaved to the streamfunction 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 6 (see Eqs. (8) and (11) With this rescaling the leading order contribution to (20) reads which is identical to the model derived by Radko 6 in the limit of small τ .
In our simulations we found that the characteristic length of the fingers increases as Ra increases (cf. Figs. 5 and 7), in contrast to a result of Merryfield and Grinder (cf. Fig.   8 in 39 ). Evidently salt-finger convection exhibits different limiting parameter regimes: in our model τ and R −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 and use the Regime II scalings to obtain Nu S ∼Ra in the limit of large Ra. This result is to be compared with the asymptotic result Nu S ∼Ra 1/3 obtained by Yang et al. 22 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 Ra T is of order 10 7 while the salinity Rayleigh number Ra S is of order 10 8 (see (15)). These values are comparable to those used in 22 although our density ratio is much larger than that in 22 .
Salt-finger convection often results in the formation of salinity staircases 10,32,40,41 , 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 4,7 . The γ-instability mechanism proposed by Radko 32 provides an alternative explanation but requires a nonmonotonic dependence of the salt flux on the mean salinity gradient. Since in our system this dependence is always monotonic neither mechanism for staircase formation is present. However, Brown et al. 9 recently discovered that staircases can form even when the flux increases monotonically with the mean salinity gradient raising the possibility that the Rayleigh ratio used in our simulations is insufficiently large to generate this interesting state.
In a future publication we plan to discuss the extent to which the conclusions of this paper carry over to three dimensions.
Here the maximum and minimum are taken over the interval Ra − 1 ∈ [e −2 , 1] corresponding to Regime II 1 (cf. Fig. 8). 1/20, which are also negligible. Therefore we take the optimal choice of the exponent α = 1/4 and use it to generate the black lines for Regime II 1 in both panels of Fig. 8. We need to mention that the power-law expressions for Regime II 1 are only numerical approximations which are not asymptotic.