2. Problem formulation

We are interested in the dynamics of salt fingers in the simultaneous limits of fast thermal diffusion and weak or moderate instability. We consider fluctuations atop linear background profiles of salinity and potential temperature in the vertical with constant slopes $\beta _S$ and $\beta _T$ , respectively. We assume the flows are slow enough and the layer height small enough to permit the use of the Boussinesq approximation. In this limit, the standard control parameters include the Prandtl number $\textit{Pr} \equiv \nu /\kappa _T$ , the inverse Lewis number $\tau \equiv \kappa _S/\kappa _T$ , and the density ratio $R_\rho \equiv \alpha _T \beta _T / (\alpha _S \beta _S)$ with $\alpha _T\gt 0$ , $\alpha _S\gt 0$ being the respective coefficients of expansion. We consider periodic boundary conditions in all directions, in which case our system is linearly unstable to the salt-finger instability for $1 \lt R_\rho \lt \tau ^{-1}$ (Baines & Gill Reference Baines and Gill1969), with $R_\rho = \tau ^{-1}$ corresponding to marginal diffusive stability and $R_\rho \lt 1$ to an unstably stratified background and hence dynamical instability. In the regime of interest here, it is helpful to introduce the following control parameters:

(2.1) \begin{align} \mathcal{R} \equiv \frac {\alpha _S \beta _S \kappa _T}{\alpha _T \beta _T \kappa _S} = \frac {1}{R_\rho \tau }, \quad \varepsilon \equiv \mathcal{R} - 1 \quad \text{and} \quad Sc \equiv \frac {\nu }{\kappa _S} = \frac {\textit{Pr}}{\tau }, \end{align}

where $\mathcal{R}$ is the Rayleigh ratio (with marginal stability now at $\mathcal{R} = 1$ ), $\varepsilon$ is a measure of supercriticality and $Sc$ is the Schmidt number. In all results presented, we fix $\textit{Pr} = 5$ , $\tau = 0.05$ and thus $Sc = 100$ .

We follow § 3.1 of Xie et al. (Reference Xie, Miquel, Julien and Knobloch2017) in our choice of non-dimensionalisation, taking the characteristic finger width $d = [\kappa_T \nu / (g \alpha_T \beta_T)]^{1/4}$ (with gravitational acceleration $g$ ) as the length scale and the salinity diffusion time $d^2/\kappa _S$ as the time scale. As our temperature scale, we take the background temperature difference across a height $d$ rescaled by $\tau$ , yielding $\tau \beta _T d$ . Similarly, our unit for salinity fluctuations is the background salinity difference across $d$ rescaled by $\mathcal{R}^{-1}$ , yielding $\mathcal{R}^{-1} \beta _S d$ . The governing equations, with hats over dependent and independent variables (and derivatives) to indicate this choice of non-dimensionalisation, are

(2.2) \begin{align} \frac {1}{Sc} \left ( \frac {\partial }{\partial \hat t} + \hat {\boldsymbol{u}} \boldsymbol{\cdot} \hat{\boldsymbol\nabla }\right ) \hat {\boldsymbol{u}} = - \hat {\boldsymbol\nabla }\hat p + \big(\hat T - \hat S \big) \boldsymbol{e}_z + \hat {\nabla} ^2 \hat {\boldsymbol{u}}, \\[-29pt] \nonumber \end{align}

(2.3) \begin{align} \hat {\boldsymbol \nabla } \boldsymbol{\cdot} \hat {\boldsymbol{u}} = 0, \\[-29pt] \nonumber \end{align}

(2.4) \begin{align} \tau \left ( \frac {\partial }{\partial \hat t} + \hat {\boldsymbol{u}} \boldsymbol{\cdot }\hat {\boldsymbol\nabla }\right ) \hat T + \hat w = \hat {\nabla} ^2 \hat T \\[2pt] \nonumber \end{align}

and

(2.5) \begin{align} \left ( \frac {\partial }{\partial \hat t} + \hat {\boldsymbol{u}} \boldsymbol{\cdot }\hat {\boldsymbol \nabla}\right ) \hat S + \mathcal{R}\, \hat w = \hat {\nabla} ^2 \hat S, \end{align}

with unit vector $\boldsymbol{e}_z$ in the $z$ direction, velocity $\hat {\boldsymbol{u}} = (\hat {u}, \hat {v}, \hat w)$ and temperature, salinity and pressure fluctuations $\hat T$ , $\hat S$ and $\hat p$ . Note that $\hat T$ and $\hat S$ represent fluctuations atop the imposed background profiles of potential temperature and salinity, which vary linearly with depth. These equations are identical to those that have been used to study salt fingers in triply periodic boxes elsewhere (see e.g. Radko Reference Radko2013; Garaud Reference Garaud2018), albeit with a different non-dimensionalisation, and thus with different control parameters appearing in different places.

In what follows, we present DNS of the system for different values of $\mathcal{R}$ . These simulations were performed using the pseudospectral-tau method implemented in Dedalus v2 (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). We use periodic boundary conditions with a horizontal domain size of $(L_x, L_y) = (4 \times 2\unicode{x03C0} /\hat {k}_{\textit{opt}}, 2 \times 2\unicode{x03C0} /\hat {k}_{\textit{opt}})$ , where $\hat {k}_{\textit{opt}}(\textit{Pr}, \tau , R_\rho )$ is the horizontal wavenumber at which the linear instability has the largest growth rate for a given set of parameters (for a subset of our parameters, we have checked that the dynamics and time-averaged fluxes change negligibly upon increasing $L_x$ and $L_y$ ). For the parameters of interest, the salt fingers become very extended in $z$ (cf. Fraser et al. Reference Fraser, Reifenstein and Garaud2024), and thus our domain height must be very large to avoid artificial domain-size effects (see e.g. Appendix A.3 of Traxler et al. Reference Traxler, Stellmach, Garaud, Radko and Brummell2011). For most of the cases reported here, we find $L_z = 64 \times 2\unicode{x03C0} /\hat {k}_{\textit{opt}}$ to be sufficient (where we confirm this by comparing against results with $L_z = 128 \times 2\unicode{x03C0} /\hat {k}_{\textit{opt}}$ ), although we find that shorter domains suffice for $\varepsilon \gtrsim 1$ , while taller domains are necessary for $\varepsilon \lesssim 1/8$ . We dealias using the standard 3/2-rule and use a numerical resolution of eight Fourier modes per $2 \unicode{x03C0} /\hat {k}_{\textit{opt}}$ in each direction, although twice this resolution becomes necessary (as verified by convergence checks) for our simulations with the largest values of $\varepsilon$ . Note that, per our dealiasing procedure, nonlinearities are evaluated on a grid with $3/2$ times this resolution. We time step using a semi-implicit, second-order Adams–Bashforth/backwards difference scheme (Dedalus’ ‘SBDF2’ option, (2.8) of Wang & Ruuth Reference Wang and Ruuth2008), with nonlinearities treated explicitly and all other terms treated implicitly, and an advective Courant–Friedrichs–Lewy safety factor of $0.3$ (sometimes $0.15$ for our highest $\varepsilon$ cases). We initialise simulations with small-amplitude noise in $\hat {S}$ , and all other fields set to zero.

3. Trends across $\boldsymbol{\varepsilon}$

Figure 1. Flow velocity snapshots at $y=0$ in the saturated state from simulations of (2.2)–(2.5) with varying supercriticality: $\varepsilon = 1/80$ (a–c), $\varepsilon = 1/10$ (d–f) and $\varepsilon = 1$ (g–i), with time traces of the corresponding salinity flux $|\hat {F}_S|$ (blue solid lines) shown in (j), (k) and (l), respectively, alongside fluxes from different reduced models (orange dashed and red dash-dotted lines, see § 4). All cases exhibit a multiscale and anisotropic flow where fingers with large vertical extent and vertical velocity (compared with horizontal width and velocity) coexist with small-scale isotropic disturbances. Magenta curves (e,f,h,i) show the time-average (over the saturated state) of the horizontal, helical mean flow $\hat{\boldsymbol{U}}_\perp = (\hat{U}(z), \hat{V}(z), 0)$ that becomes a significant feature for $\varepsilon \gtrsim 0.1$ .

Figure 1 shows velocity snapshots from simulations with $\varepsilon = 1/80$ (panels a–c), $1/10$ (panels d–f) and $1$ (panels g–i), i.e. $R_\rho \simeq 19.75$ , $18.18$ and $10$ , illustrating general trends in this regime. In each case, we see highly anisotropic and multiscale dynamics, with vertically elongated, large-amplitude structures (the characteristic salt fingers) in $w$ , and smaller-amplitude, isotropic eddies seen in each velocity field. The separation between the long vertical and the short isotropic scales shrinks as $\varepsilon$ increases.

At very small $\varepsilon$ (see panels a–c), the fingers become vertically invariant ‘elevator modes’ disturbed by isotropic ripples. For moderate supercriticality, $\varepsilon \sim 0.1$ , the fingers are no longer vertically invariant but still very anisotropic, with much larger vertical scale and velocity than in the horizontal. We remark that at larger $\varepsilon$ , self-connecting structures only persist for insufficiently tall domains and lead to bursty and domain height-dependent dynamics. At very small $\varepsilon$ , self-connecting structures persist in even the tallest domains we can reasonably achieve numerically, but they do not drive bursty dynamics or domain height-dependent dynamics.

Figure 2. Relative helicity (see text) of the mean flow (blue) and of the fluctuations about the mean flow (orange; multiplied by $10^4$ to ease comparison) for two values of $\varepsilon$ . At small $\varepsilon$ , the flow is almost maximally helical, and in both cases the fluctuations are almost non-helical, with $H_{\textit{rel}}[\hat {\boldsymbol{u}}'] \sim 10^{-5}$ .

In this regime, a horizontal mean flow $\hat{\boldsymbol{U}}_\perp = (\hat{U}(z), \hat{V}(z), 0)$ develops spontaneously (cf. Liu et al. Reference Liu, Sharma, Julien and Knobloch2024), as shown by the magenta curves in panels (e,f). The two components of $\hat {\boldsymbol{U}}_\perp$ are $\unicode{x03C0} /2$ out of phase in $z$ – i.e. one component passes through $0$ as the other reaches an extremum – resulting in a mean flow with non-zero helicity $H[\hat {\boldsymbol{U}}_\perp ]$ , where

(3.1) \begin{align} H[\hat {\boldsymbol{u}}] = \int \hat {\boldsymbol{u}} \boldsymbol{\cdot } \big(\hat {\boldsymbol \nabla }\times \hat {\boldsymbol{u}}\big) \textrm {d}V \end{align}

defines the helicity of a given flow field $\hat {\boldsymbol{u}}$ . As the mean flow has no vertical component, its helicity arises from the horizontal components of $\hat {\boldsymbol{U}}_\perp$ and $\hat {\boldsymbol \nabla }\times \hat {\boldsymbol{U}}_\perp$ , and it advects the fingers into a corkscrew-like shape (the graphical abstract accompanying this article shows a volume rendering of the vertical velocity for the same $\mathcal{R} = 2$ simulation as shown in figure 1(g–i) and provides a visualisation of this flow). In fact, this mean flow is nearly a maximally helical (Beltrami) flow in that its relative helicity $H_{\textit{rel}}[\hat {\boldsymbol{U}}_\perp ] \equiv H[\hat {\boldsymbol{U}}_\perp ]/[(\int |\hat {\boldsymbol{U}}_\perp |^2\, \textrm {d}V) (\int |\hat {\boldsymbol \nabla }\times \hat {\boldsymbol{U}}_\perp |^2 \textrm {d}V)]^{1/2}$ , calculated using the time- and horizontally averaged flow, is over 0.99 – alternatively, $H_{\textit{rel}}$ of the instantaneous mean flow saturates at roughly $0.8$ , see figure 2. While $H_{\textit{rel}} \gt 0$ for this particular simulation, other random initial conditions lead to $H_{\textit{rel}} \lt 0$ (not shown) with no clear statistical preference between the two signs based on our limited sample. In stark contrast to the strongly helical mean flow, the relative helicity of the fluctuations about the mean ( $\hat {\boldsymbol{u}}' \equiv \hat {\boldsymbol{u}} - \hat {\boldsymbol{U}}_\perp$ ) is roughly $10^{-5}$ . This leads to the remarkable observation that the system spontaneously forms a symmetry-breaking, maximally helical flow from non-helical fluctuations. In fact the generation of such mean flows by salt fingers was already seen in vertically bounded domains by Yang et al. (Reference Yang, Verzicco and Lohse2016) (see their figure 16), albeit without comment (see also the discussion of fingers ‘twisting about one another’ in Schmitt & Lambert Reference Schmitt and Lambert1979). Similar spontaneous mirror symmetry breaking and Beltrami flow formation have also been observed in other systems, including active matter (Słomka & Dunkel Reference Słomka and Dunkel2017; Romeo et al. Reference Romeo, Słomka, Dunkel and Burns2024) and magnetohydrodynamic (Agoua et al. Reference Agoua, Favier, Delache, Briard and Bos2021) turbulence, but without such strikingly non-helical fluctuations.

For yet larger $\varepsilon$ (see figure 1 (g–i)), both the mean flow and the fingers become more vigorous, and the anisotropy of the fingers is less extreme, permitting shorter vertical domains. In this regime, the mean flow is still very helical at each $z$ , but tends to have a shorter vertical scale and its helicity may change sign with $z$ . The helicity is also less stationary in this regime, and for some initial conditions it is observed to change sign with time, similar to the finite lifetimes of different flow patterns observed in convection (see e.g. Ahlers, Cannell & Steinberg Reference Ahlers, Cannell and Steinberg1985; Winchester, Dallas & Howell Reference Winchester, Dallas and Howell2021; Wang, Goluskin & Lohse Reference Wang, Goluskin and Lohse2023).

Both the multiscale aspect of this system and the trends in scale separation with $\varepsilon$ are readily seen in figure 3, which shows spectra of the kinetic energy as a function of the vertical wavenumber $\hat {k}_z$ at $\hat {k}_\perp = \hat {k}_{\textit{opt}}$ , the fastest-growing wavenumber of the linear instability. Two distinguishing features are seen most clearly in the horizontal kinetic energy, which has a local maximum (or otherwise a clear change in the spectrum, in the case of large $\varepsilon$ ) at isotropic scales where $\hat {k}_z\sim \hat {k}_{\textit{opt}}$ , indicated by the black vertical lines, and a local maximum at smaller $\hat {k}_z$ indicated by the red vertical lines. Note that the gap in $\hat {k}_z$ between these two local maxima shrinks as $\varepsilon$ increases, consistent with the observed decrease in scale separation with increasing $\varepsilon$ (figure 1).

Figure 3. Horizontal (blue) and vertical (orange) kinetic energy spectra (time-averaged over the statistically stationary state) versus $\hat {k}_z$ at $\hat {k}_y = 0$ and $\hat {k}_x = \hat {k}_{\textit{opt}}$ . Black lines show $\hat {k}_z = \hat {k}_{\textit{opt}}$ to highlight the small-scale isotropic flow component while the red vertical lines correspond to the secondary peak in the horizontal spectrum to highlight the anisotropic, small $\hat {k}_z$ flow component. The ratio between these two wavenumbers provides one measure of anisotropy shown in figure 4.

4. Asymptotic models

The IFSC model of Xie et al. (Reference Xie, Miquel, Julien and Knobloch2017), in three dimensions (3-D), is appropriate when $\tau \ll 1$ and $Sc \gg 1$ and is described by the equations

(4.1) \begin{align} \hat w = \hat {\nabla} ^2 \hat T, \\[-29pt] \nonumber \end{align}

(4.2) \begin{align} 0 = - \hat {\boldsymbol \nabla }\hat p + (\hat T - \hat S) \boldsymbol{e}_z + \hat {\nabla} ^2 \hat {\boldsymbol{u}}, \\[2pt] \nonumber \end{align}

with (2.3) and (2.5) left unchanged. While it appears that the model should be valid for all order-one $\varepsilon$ , we show in figure 1 that this is not the case: the model does produce dynamics and fluxes consistent with the full system at sufficiently small supercriticality, roughly $\varepsilon \lesssim 1/80$ , but, for $\varepsilon \approx 1/20$ or larger, it produces dynamics that differ qualitatively from solutions of the full equations – while the full system exhibits helical mean flows that disrupt and twist the fingers on large scales, the IFSC model removes the Reynolds stress term from the horizontal mean of (2.2) and thus lacks these flows. Without mean flows to disrupt the long fingers, such fingers drive temporally bursty dynamics that dramatically raise the fluxes, as shown by the IFSC curve in figure 1(k).

Figure 4. Scalings with respect to $\varepsilon$ of several quantities (indicated in the caption for each panel) for the full system, (2.2)–(2.5) (blue dots), the IFSC model, with (2.4) and (2.2) replaced by (4.1) and (4.2) (green diamonds), and the MIFSC model, where (2.2) is replaced instead by (4.14)–(4.15) (orange crosses). Black dashed lines show scalings predicted by the multiscale asymptotic analysis described in the text. The green dashed lines and the two measures of anisotropy are described in the text.

However, the IFSC model can be modified to capture mean flow generation. Our simulations suggest that anisotropy is an essential aspect of a reduced description of salt fingers valid in the regime depicted in figure 1, possibly with a rescaling of the various fields to retain the Reynolds stress at leading order. In order to capture both the elongated fingers shown in figure 1 and the small-scale isotropic fluctuations therein, we employ a multiscale asymptotic analysis inspired by a related approach to turbulence in stably stratified fluids employed by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024) and Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024).

To begin, we note that the growth rate and optimal wavenumber of the linear instability scale as $\varepsilon ^{3/2}$ and $\varepsilon ^{1/4}$ , respectively, for sufficiently small $\varepsilon$ – for the $\textit{Pr}$ , $\tau$ considered here, this scaling is achieved for $\varepsilon \lesssim 1$ (see e.g. figure 3 of Xie et al. Reference Xie, Miquel, Julien and Knobloch2017). It is thus convenient to rescale the independent and dependent variables as follows (cf. Radko Reference Radko2010):

(4.3) \begin{align} \hat {\boldsymbol{x}} = \varepsilon ^{-1/4} \boldsymbol{x}, \quad \hat t = \varepsilon ^{-3/2} t, \quad \hat {\boldsymbol{u}} = \varepsilon ^{5/4} \boldsymbol{u}, \quad \hat p = \varepsilon ^{3/2} p, \quad (\hat T, \hat S) = \varepsilon ^{3/4} (T,S). \end{align}

Equations (2.2)–(2.5) then become

(4.4) \begin{align} \varepsilon \frac {1}{Sc} \left ( \frac {\partial }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{u} = - \boldsymbol{\nabla }p + \varepsilon ^{-1} (T - S) \boldsymbol{e}_z + {\nabla} ^2 \boldsymbol{u}, \\[-29pt] \nonumber \end{align}

(4.5) \begin{align} \varepsilon \tau \left ( \frac {\partial }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) T + w = {\nabla} ^2 T \\[2pt] \nonumber \end{align}

and

(4.6) \begin{align} \varepsilon \left ( \frac {\partial }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) S + \mathcal{R}\, w = {\nabla} ^2 S, \end{align}

with (2.3) left unchanged. When $\varepsilon =O(1)$ but $\tau \ll 1$ and $Sc\gg 1$ , the inertial terms in (4.4)–(4.5) drop out and the resulting equations correspond to the 3-D IFSC model with (4.6) providing the sole prognostic equation. On the other hand, when $\varepsilon \ll 1$ , we may expand all fields as asymptotic series in $\varepsilon$ as $q = \sum _n q_n \varepsilon ^n$ . Inspecting the $z$ component of (4.4) in the limit $\varepsilon \to 0$ shows that $T_0 = S_0$ , i.e. the dynamics are neutrally buoyant in this limit. In the following it is helpful to introduce the buoyancy $b \equiv \varepsilon ^{-1}(T-S)$ and subtract (4.5) from (4.6), yielding

(4.7) \begin{align} \left ( \frac {\partial }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) (S - \tau T) + w = -{\nabla} ^2 b. \end{align}

We now perform a multiscale asymptotic expansion employing the following four steps:

1. (i) Introduce fast and slow scales in $z$ and $t$ and allow each field $q$ to depend on them both, i.e. $q = q(\boldsymbol{x}_\perp , z_s, z_f, t_s, t_f)$ , with $\partial _z \mapsto \alpha ^{-1} \partial _{z_s} + \partial _{z_f}$ and $\partial _t \mapsto \alpha ^{-1} \partial _{t_s} + \partial _{t_f}$ ; here the fast vertical scale $z_f$ is isotropic, and the anisotropy of the slow vertical scale $z_s$ is controlled by $\alpha$ , which we take to be $\alpha = Sc/\varepsilon \gg 1$ (which we note is qualitatively consistent with the trends seen in figures 1 and 3).

2. (ii) Decompose each field $q$ into its fast-averaged and fluctuating components: $q = \left \langle q \right \rangle _f + \tilde {q}$ where $\left \langle \boldsymbol{\cdot }\right \rangle _f$ denotes an average over $z_f$ and $t_f$ , so that $\left \langle q \right \rangle _f$ depends only on $\boldsymbol{x}_\perp$ , $z_s$ and $t_s$ , while $\tilde {q}$ depends on all coordinates. The only exception to this $q = \left \langle q \right \rangle _f + \tilde {q}$ decomposition is for the pressure and horizontal flow fluctuations (defined as $\boldsymbol{u}'_\perp \equiv \boldsymbol{u}_\perp - \boldsymbol{U}_\perp$ ), where we instead take $p = \alpha ^{-1} \left \langle p \right \rangle _f + \tilde {p}$ and $\boldsymbol{u}'_\perp = \alpha ^{-1} \langle \boldsymbol{u}'_\perp \rangle _f + \tilde {\boldsymbol{u}}'_\perp$ . (Note that this is without loss of generality: we are simply taking different non-dimensionalisations for $\left \langle p \right \rangle _f$ and $\tilde {p}$ and likewise for $\langle \boldsymbol{u}'_\perp \rangle _f$ and $\tilde {\boldsymbol{u}}'_\perp$ .)

3. (iii) Take the fast average of each of the governing equations to obtain an equation for $\partial _{t_s} \left \langle q \right \rangle _f$ , and subtract the result from the original equation to obtain an equation for $\partial _t \tilde {q}$ .

4. (iv) Finally, retain only leading-order terms by removing terms that are small in the limits $\varepsilon \to 0$ , $\tau \to 0$ , $Sc \to \infty$ and $\alpha \to \infty$ (with $\alpha \varepsilon /Sc = 1$ ).

An explicit illustration of these steps is provided in, for example, Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) and Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024).

The resulting equations are

(4.8) \begin{align} \left \langle w \right \rangle _f = {\nabla} _\perp ^2 \left \langle T \right \rangle _f \quad \text{and} \quad \tilde {w} = {\nabla} ^2_f \tilde {T}, \\[-12pt] \nonumber \end{align}

(4.9) \begin{align} \left \langle \boldsymbol{U}_\perp \right \rangle _f \boldsymbol{\cdot }\boldsymbol{\nabla} _\perp \left \langle S \right \rangle _f + \big \langle \tilde {\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla} _f \tilde {S} \big \rangle _f + \left \langle w \right \rangle _f = - {\nabla} _\perp ^2 \left \langle b \right \rangle _f, \\[-29pt] \nonumber \end{align}

(4.10) \begin{align} \left [ \frac {\partial }{\partial t_f} + \left \langle \boldsymbol{U}_\perp \right \rangle _f \boldsymbol{\cdot }\boldsymbol{\nabla} _\perp + \left \langle w \right \rangle _f \frac {\partial }{\partial z_f} \right ] \tilde {S} + \tilde {\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla} _f \tilde {S} & - \big \langle \tilde {\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla} _f \tilde {S} \big \rangle _f \\ \nonumber & + \tilde {\boldsymbol{u}}_\perp \boldsymbol{\cdot }\boldsymbol{\nabla} _\perp \left \langle S \right \rangle _f + \tilde {w} = - {\nabla} _f^2 \tilde {b}, \\[-29pt] \nonumber \end{align}

(4.11) \begin{align} 0 = - \boldsymbol{\nabla} _\perp \left \langle p \right \rangle _f + \left \langle b \right \rangle _f \boldsymbol{e}_z + {\nabla} _\perp ^2 \langle \boldsymbol{u}' \rangle _f \quad \text{and} \quad 0 = - \boldsymbol{\nabla} _f \tilde {p} + \tilde {b} \boldsymbol{e}_z + {\nabla} _f^2 \tilde {\boldsymbol{u}}', \\[2pt] \nonumber \end{align}

(4.12) \begin{align} \frac {\partial }{\partial t_s} \left \langle \boldsymbol{U}_\perp \right \rangle _f + \frac {\partial }{\partial z_s}\langle \tilde {\boldsymbol{u}}_\perp \tilde {w}\rangle _{f,\perp } = \frac {\partial ^2}{\partial z_s^2} \left \langle \boldsymbol{U}_\perp \right \rangle _f, \\[-29pt] \nonumber \end{align}

(4.13) \begin{align} \boldsymbol{\nabla} _s \boldsymbol{\cdot } \langle \boldsymbol{u}' \rangle _f = 0 \quad \text{and} \quad \boldsymbol{\nabla} _f \boldsymbol{\cdot }\tilde {\boldsymbol{u}}' = 0, \\[6pt] \nonumber \end{align}

where we have introduced $\boldsymbol{\nabla} _s \equiv (\partial _x, \partial _y, \partial _{z_s})$ and $\boldsymbol{\nabla} _f \equiv (\partial _x, \partial _y, \partial _{z_f})$ , with $\langle \boldsymbol{\cdot }\rangle _{f, \perp }$ denoting an average over $\boldsymbol{x}_\perp$ in addition to $z_f$ and $t_f$ . The structure of this system of equations is broadly similar to the IFSC model: the temperature equations reduce to a diagnostic balance between diffusion and advection of the background, the salinity equations retain nonlinearity on fast scales (and so do not yield a quasilinear structure as in the work of Chini et al. Reference Chini, Michel, Julien, Rocha and Caulfield2022), and the momentum equations for fluctuations about the mean flow involve a dominant balance between pressure, buoyancy and viscosity. However, in contrast to the IFSC model, (4.12) retains the Reynolds stress term absent from (4.2). Thus, modifying the IFSC model to capture the leading-order dynamics of the full equations in this limit merely requires retaining the Reynolds stress in the $\boldsymbol{k}_\perp = 0$ component of the momentum equation:

(4.14) \begin{align} \frac {1}{Sc} \left (\frac {\partial }{\partial \hat t} \hat {\boldsymbol{U}}_\perp + \frac {\partial }{\partial \hat z} \langle \hat {\boldsymbol{u}}'_\perp \hat {w}'\rangle _\perp \right ) = \frac {\partial ^2}{\partial \hat z^2} \hat {\boldsymbol{U}}_\perp , \\[-12pt] \nonumber \end{align}

(4.15) \begin{align} 0 = - \hat {\boldsymbol \nabla }\hat p' + \big(\hat T' - \hat S' \big) \boldsymbol{e}_z + \hat {\nabla} ^2 \hat {\boldsymbol{u}}', \\[9pt] \nonumber \end{align}

with temperature given by (4.1), and (2.3) and (2.5) retained in full. Figure 1 shows that this MIFSC model captures the same dynamics as the full equations even for $\varepsilon = O(1)$ .

The rescaling applied to arrive at the multiscale asymptotic system offers a natural suggestion for the scaling of the various fields in this limit. For (2.2)–(2.5), this analysis predicts that both $\hat{S}$ and $\hat{T}$ should scale as $\varepsilon ^{3/4}$ , $\hat {\boldsymbol{U}}_\perp$ and $\hat{w}$ as $\varepsilon ^{5/4}$ , and the scaling of $\hat {\boldsymbol{u}}'_\perp$ should differ between fast and slow scales in $\hat{z}$ , with $ \langle \boldsymbol{u}'_\perp \rangle _f \sim \varepsilon ^{9/4}$ and $\tilde {\boldsymbol{u}}'_\perp \sim \varepsilon ^{5/4}$ . The predicted stronger $\varepsilon$ dependence of $ \langle \boldsymbol{u}'_\perp \rangle _f$ is consistent with figure 3, which shows that the horizontal kinetic energy peaks at large $\hat {k}_z$ for small $\varepsilon$ but at small $\hat {k}_z$ for large $\varepsilon$ .

We present a more quantitative comparison between these predictions and DNS in figure 4 by calculating the root-mean-square of $\hat {\boldsymbol{u}}'_\perp$ , $\hat {\boldsymbol{U}}_\perp$ , $\hat w$ , $\hat S$ and $\hat T$ (denoted by subscript ‘rms’ in figure 4; here, $\hat{\boldsymbol{u}}'_{\perp,rms} \equiv \sqrt{\hat{u}'^{2}_{\textit{rms}}+ \hat{v}'^{2}_{\textit{rms}}} $ etc.), and the volume-averaged fluxes $\hat F_S = \langle \hat w \hat S \rangle$ and $\hat F_T = \langle \hat w \hat T \rangle$ , with the blue dots corresponding to DNS of the full equations, green diamonds to the IFSC model, and orange crosses to the MIFSC model. Anisotropy is quantified by two means: using the wavenumber ratio corresponding to the two spectral peaks in figure 3 (crosses; the lowest $\varepsilon$ values are suppressed because the low- $\hat {k}_z$ peak is difficult to identify and is likely constrained by domain size), and the ratio $\hat w_{\textit{rms}}/\hat {\boldsymbol{u}}'_{\perp , \textit{rms}}$ (dots). In each panel, black dashed lines correspond to the predicted scalings (for $\hat {\boldsymbol{u}}'_{\perp , \textit{rms}}$ , the black line shows the predicted $\tilde {\boldsymbol{u}}'_\perp \sim \varepsilon ^{5/4}$ scaling while the green line shows the $ \langle \boldsymbol{u}'_\perp \rangle _f \sim \varepsilon ^{9/4}$ scaling). The scalings of most of these quantities are consistent with predictions as $\varepsilon \to 0$ , especially $\hat w_{\textit{rms}}$ , $\hat S_{\textit{rms}}$ , $\hat T_{\textit{rms}}$ , $\hat F_S$ and $\hat F_T$ , with the anisotropy measurement inconclusive. In contrast, the rms of the mean flow, $\hat {\boldsymbol{U}}_{\perp , \textit{rms}}$ , follows the predicted $\varepsilon ^{5/4}$ scaling at larger $\varepsilon$ only; at small $\varepsilon$ , its scaling is perhaps more consistent with $\varepsilon ^{5/2}$ (green line), although it is difficult to measure $\hat {\boldsymbol{U}}_{\perp , \textit{rms}}$ accurately at the smallest $\varepsilon$ values as the mean flow is then very weak and evolves very slowly.