2. Formulation

To explore the circulation patterns in an effectively unbounded doubly stratified fluid, the total temperature and salinity fields $(T_{tot}^\ast ,S_{tot}^\ast )$ are separated into the background state $\,(T_{bg}^\ast ,S_{bg}^\ast ),$ representing uniform vertical gradients, and a departure $({T^\ast },{S^\ast })$ from them:

(2.1)\begin{equation}\left. {\begin{array}{*{20}{@{}c@{}}} {T_{tot}^\ast= T_{bg}^\ast+ {T^\ast } = {A_T}{z^\ast } + {A_{T0}} + {T^\ast },}\\ {S_{tot}^\ast= S_{bg}^\ast+ {S^\ast } = {A_S}{z^\ast } + {A_{S0}} + {S^\ast },} \end{array}} \right\}\end{equation}

where $({A_T},{A_{T0}},{A_S},{A_{S0}})$ are constants and the asterisks denote dimensional field variables. Our focus is on finger-favourable stratification, where $\partial T_{bg}^\ast{/}\partial {z^\ast } = {A_T} > 0$ and $\partial S_{bg}^\ast{/}\partial {z^\ast } = {A_S} > 0$. We ignore planetary rotation, compressibility and the nonlinearity of the equation of state and express the governing Boussinesq equations in terms of perturbations $({T^\ast },{S^\ast })$:

(2.2)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {\dfrac{{\partial {T^\ast }}}{{\partial {t^\ast }}} + {\boldsymbol{v}^\ast }\boldsymbol{\cdot }\boldsymbol{\nabla }{T^\ast } + {A_T}{w^\ast } = {k_T}{\nabla^2}{T^\ast },}\\ {\dfrac{{\partial {S^\ast }}}{{\partial {t^\ast }}} + {\boldsymbol{v}^\ast }\boldsymbol{\cdot }\boldsymbol{\nabla }{S^\ast } + {A_S}{w^\ast } = {k_S}{\nabla^2}{S^\ast },}\\ {\dfrac{{\partial {\boldsymbol{v}^\ast }}}{{\partial {t^\ast }}} + {\boldsymbol{v}^\ast }\boldsymbol{\cdot }\boldsymbol{\nabla }{\boldsymbol{v}^\ast } =- \dfrac{1}{{\rho_0^\ast }}\boldsymbol{\nabla }{p^\ast } + g(\alpha {T^\ast } - \beta {S^\ast })\boldsymbol{k} + \upsilon {\nabla^2}{\boldsymbol{v}^\ast },}\\ {\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}^\ast } = 0,} \end{array}} \right\}\end{equation}

where ${\boldsymbol{v}^\ast } = ({u^\ast },{v^\ast },{w^\ast })$ is the velocity, $\boldsymbol{k}$ is the vertical unit vector, ${p^ \ast }$ is the dynamic pressure, g is gravity, $(\alpha ,\beta )$ are the thermal expansion and haline contraction coefficients and $\rho _0^\ast $ is the reference density.

We adopt the traditional double-diffusive non-dimensionalization (e.g. Radko Reference Radko2013) based on scales of individual salt fingers. Thus, $l = {({k_T}\upsilon /g\alpha {A_T})^{1/4}}$, ${k_T}/l$, ${l^2}/{k_T}$ and $\rho _0^\ast \upsilon {k_T}/{l^2}$ are used as the units of length, velocity, time and pressure, respectively, whereas temperature and salinity are non-dimensionalized as follows:

(2.3a,b)\begin{equation}{T^\ast } = {A_T}l\,T,\quad {S^\ast } = \frac{\alpha }{\beta }{A_T}l\,S.\end{equation}

After non-dimensionalization, the governing equations are reduced to

(2.4)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {\dfrac{{\partial T}}{{\partial t}} + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }T + w = {\nabla^2}T,}\\ {\dfrac{{\partial S}}{{\partial t}} + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }S + \dfrac{w}{{{R_\rho }}} = \tau {\nabla^2}S,}\\ {\dfrac{1}{{Pr}}\left( {\dfrac{\partial }{{\partial t}}\boldsymbol{v} + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla v}} \right) =- \boldsymbol{\nabla }p + (T - S)\boldsymbol{k} + {\nabla^2}\boldsymbol{v},}\\ {\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v} = 0,} \end{array}} \right\}\end{equation}

where ${R_\rho } = \alpha {A_T}/\beta {A_S}$ is the background density ratio, $\tau = {k_S}/{k_T}$ is the diffusivity ratio and $Pr = \upsilon /{k_T}$ is the Prandtl number. Some calculations in this study are performed in two dimensions (x, z), which reduces the Boussinesq system (2.4) to

(2.5)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {\dfrac{{\partial T}}{{\partial t}} + J(\psi ,T) + \dfrac{{\partial \psi }}{{\partial x}} = {\nabla^2}T,}\\ {\dfrac{{\partial S}}{{\partial t}} + J(\psi ,S) + \dfrac{1}{{{R_\rho }}}\dfrac{{\partial \psi }}{{\partial x}} = \tau {\nabla^2}S,}\\ {\dfrac{\partial }{{\partial t}}{\nabla^2}\psi + J(\psi ,{\nabla^2}\psi ) = Pr\left[ {\dfrac{\partial }{{\partial x}}(T - S) + {\nabla^4}\psi } \right]\!,} \end{array}} \right\}\end{equation}

where $\psi $ is the streamfunction associated with the velocity field $(u,w) = ( - \partial \psi /\partial z, \partial \psi /\partial x)$, and $J(a,b) \equiv (\partial a/\partial x)(\partial b/\partial z) - (\partial a/\partial z)(\partial b/\partial x)$ is the Jacobian.

3. Sub-microscale dynamics of double-diffusive convection

To illustrate the key features of fingering convection, we present a typical 2-D DNS. The governing equations (2.5) are integrated in time using the Fourier-based spectral model (e.g. Stern et al. Reference Stern, Radko and Simeonov2001; Radko Reference Radko2019a) with periodic boundary conditions assumed for $(\psi ,T,S)$ in each direction. The periodicity of perturbations reflects the effectively unbounded character of double diffusion in the ocean. Salt fingers are most common at depths of several hundred metres and therefore are unaffected by the surface and sea-floor boundaries. The unbounded model of double diffusion has been validated by oceanic observations. The most recent and comprehensive analysis performed for a wide range of environmental conditions (Brown & Radko Reference Brown and Radko2024) unambiguously confirmed the model's ability to represent field measurements and the lack of systematic biases.

The calculation in figure 1 was performed with the oceanographically relevant (e.g. Radko Reference Radko2013) values of governing parameters $(\textit{Pr} ,\tau ,{R_\rho }) = (7,0.01,2)$ on the computational domain of size ${L_x} \times {L_z} = 50 \times 50\,,$ resolved by $({N_x},{N_z}) = (2048,2048)$ grid points. The experiment was initiated by the random low-amplitude distribution of $(\psi ,T,S)$, introduced to seed primary instabilities. The experiment is extended for $t = 1000$ time units – dimensionally equivalent to 10 days for oceanic scales. The first stage in the development of instability is marked by the spontaneous emergence of vertically oriented salt fingers. The temperature and salinity patterns, shown in figure 1(a,b) at $t = 20$, are relatively regular and vary on comparable spatial scales. However, by $t = 50$, primary salt fingers develop secondary instabilities (e.g. Radko & Smith Reference Radko and Smith2012) and the flow field becomes more turbulent and disorganized (figure 1c,d). The statistical equilibration of salt fingers is completed by $t = 100$ (figure 1e,f). Afterward, the quasi-steady state is maintained indefinitely.

Figure 1. Salt finger DNS performed with $(Pr,\tau ,{R_\rho }) = (7,0.01,2)$. The instantaneous patterns of temperature (left panels) and salinity (right panels) are shown for (a,b) $t = 20$, (c,d) $t = 50$ and (e,f) $t = 100$.

A salient feature of the resulting solutions is the emergence of fine structures in the salinity field that are not reflected in the corresponding temperature and velocity patterns. To better illustrate this phenomenon, we present (figure 2a,b) a magnified view of T and S over a small area of

(3.1)\begin{equation}\varOmega = \{ - 2 < x < 2,\; - 2 < z < 2\} .\end{equation}

The sub-microscale salinity patterns appear in the form of narrow fronts and their width is controlled by the molecular dissipation of salt. The latter assumption is readily confirmed by reproducing the simulation in figure 1 with a much lower diffusivity ratio of $\tau = 5 \times {10^{ - 4}}$ (figure 2c,d). To fully resolve the ‘salinity’ spectrum, we employ a finer mesh with $({N_x},{N_z}) = (4096,4096)$ grid points. The resulting salinity pattern (figure 2d) is characterized by even narrower and more tightly packed fronts than in the baseline experiment (figure 2b). This dynamics conforms to the view expressed by Batchelor (Reference Batchelor1959), who attributed the fine features of weakly diffusive tracers to larger-scale strain.

Figure 2. The magnified view of temperature (a) and salinity (b) in the small area $\varOmega $ for the state in figure 1(e,f). Panels (c) and (d) show the corresponding patterns realized in the analogous experiment performed with $\tau = 5 \times {10^{ - 4}}$ at $t = 100$. Note the emergence of numerous sub-microscale fronts of salinity in (d).

Further insight into the origin of these sub-microscale structures is provided by the salinity variance equation:

(3.2)\begin{equation}\frac{\partial }{{\partial t}}\langle {S^2}\rangle + \frac{2}{{{R_\rho }}}\langle wS\rangle =- \chi ,\end{equation}

where angle brackets denote spatial averages and $\chi \equiv 2\tau \langle {|{\boldsymbol{\nabla }S} |^2}\rangle $. Equation (3.2) was obtained by multiplying the salinity equation in (2.5) by S and averaging the result in x and z. In statistically steady configurations, as considered in our study, the variance production term $\varPi = (2/{R_\rho })\langle wS\rangle$ is balanced by the molecular dissipation (–χ). This so-called production–dissipation balance (Osborn & Cox Reference Osborn and Cox1972) concisely describes the link between the large-scale stirring of the salinity field represented by the Π-term and the molecular dissipation on small scales (χ-term). It is also instructive to examine the production–dissipation balance for very low diffusivity ratios. In the limit $\tau \to 0$, the magnitude of the production term approaches a finite value (e.g. Radko Reference Radko2008), demanding that $|\chi |= O(1)$. This, in turn, requires salinity gradients $|{\boldsymbol{\nabla }S} |\sim \sqrt {\chi /2\tau } = O({\tau ^{ - 0.5}})$ to increase without bound with decreasing $\tau $. The amplification of salinity gradients for small diffusivity ratios is spectacularly manifested by the emergence of numerous sub-microscale fronts, clearly visible in figure 2(d).

The tendency for the formation of sub-microscale salinity patterns can also be quantified (figure 3) by examining the dissipation spectra $D(\kappa )$ that satisfy

(3.3)\begin{equation}\chi = \int {D(\kappa )\,\textrm{d}\kappa } ,\end{equation}

where $\kappa $ is the 2-D wavenumber. The spectra were computed for the baseline experiment performed with $\tau = 0.01$ (figure 3a) and for the low-diffusivity experiment based on $\tau = 5 \times {10^{ - 4}}$ (figure 3b). In both cases, spectra were averaged in time over the quasi-equilibrium interval $100 < t < 1000\,.$ The DNS results in figure 3 are compared with the canonical Batchelor (Reference Batchelor1959) spectrum, which in our non-dimensional units takes the form

(3.4)\begin{equation}{D_B}(\kappa ) = {C_B}\tau \chi \sqrt {\frac{{Pr}}{\varepsilon }} \kappa \,\textrm{exp}( - {C_B}{(\kappa {\eta _B})^2}),\end{equation}

where $\varepsilon $ is the mean non-dimensional dissipation rate of turbulent kinetic energy, ${\eta _B} = {(Pr\,{\tau ^2}{\varepsilon ^{ - 1}})^{1/4}},$ and ${C_B}$ is the Batchelor constant. The latter was evaluated from the best fit of (3.4) to the numerical data, resulting in ${C_B} = 4.01$ and ${C_B} = 3.95$ for the calculations in (a) and (b), respectively. The estimates of ${C_B}$ for temperature dissipation by oceanic turbulence (e.g. Dillon & Caldwell Reference Dillon and Caldwell1980; Oakey Reference Oakey1982) suggest that the Batchelor constant is in the range of 3.4–4.1. The consistency of these values with the salinity-based ${C_B}$ in double-diffusive experiments – which differ in many ways from the more conventional turbulence problems – underscores the remarkable universality of Batchelor's theory of scalar dissipation.

Figure 3. The dissipation spectra for the baseline ($\tau = 0.01$) and low-diffusivity ($\tau = 5 \times {10^{ - 4}}$) simulations are shown by solid curves in (a) and (b), respectively. The dashed curves represent the corresponding Batchelor (Reference Batchelor1959) spectra computed as the best fits to the DNS data.

It should be noted, however, that while the DNS spectra conform to the Batchelor model for relatively low wavenumbers, they exhibit more gentle spectral decay at high ones. Such shallow drop-off is better represented by Kraichnan's (Reference Kraichnan1974) model, which accounts for the fluctuations of the strain rate. Overall, however, the general relevance of Batchelor–Kraichnan ideas to the foregoing simulations is undeniable. Of particular interest to our discussion is the link between scalar dissipation and the strain rate espoused by these theories. To further explore this connection, we examine the statistics of the local dissipation ${\chi _{loc}} \equiv 2\tau {|{\boldsymbol{\nabla }S} |^2}$ as a function of the strain rate (s). The strain rate is defined, in both 2-D and 3-D systems, as

(3.5a,b)\begin{equation}s = \sqrt {2{s_{ij}}{s_{ij}}} ,\quad {s_{ij}} = \frac{1}{2}\left( {\frac{{\partial {v_i}}}{{\partial {x_j}}} + \frac{{\partial {v_j}}}{{\partial {x_i}}}} \right)\!,\end{equation}

where

(3.6a,b)\begin{equation}({v_1},{v_2},{v_3}) = (u,v,w),\quad ({x_1},{x_2},{x_3}) = (x,y,z).\end{equation}

For each value of strain (s), we estimate the mean dissipation (${\chi _s}$) by averaging ${\chi _{loc}}$ over all locations where the strain is in the range $[s - 0.5\varDelta s,\;s + 0.5\varDelta s]$. In the following calculation, we use $\varDelta s = 0.02\,\textrm{max}(s)$ and average the results in time over the quasi-equilibrium interval $100 < t < 1000$.

The ${\chi _s}(s)$ patterns obtained in this manner for the baseline and low-diffusivity experiments are shown in figure 4. Both simulations reveal a rapid monotonic increase in the dissipation rate with increasing strain. To illustrate the physical mechanisms governing the evolution of sub-microscale variability, we present explicit analytical solutions (Appendix A) that underscore the link between finger-scale strain and the emergence of much finer salinity patterns. We believe that this link captures the essence of the dissipation of salinity in double-diffusive systems, and it should be reflected in any physics-based subgrid closure model, such as described below.

Figure 4. The mean dissipation rate realized for a given strain is presented for the baseline experiment performed with $\tau = 0.01$ (red curve) and for the low-diffusivity simulation performed with $\tau = 5 \times {10^{ - 4}}$ (blue curve).

4. The SMF model

The simulations of fingering convection reveal the presence of salinity features that are much smaller than those found in the corresponding temperature and momentum fields. Such disparity of scales is the Achilles heel of double-diffusive DNS, which severely limits the scope of numerical investigations in this area and motivates the present effort to develop a subgrid closure for salinity. The proposed algorithm makes no attempt to parameterize momentum and temperature, which are fully resolved, and is referred to as the SMF model.

Our first attempt represented the subgrid transfer of salinity as a Fickian diffusive process with spatially variable eddy diffusivity K:

(4.1)\begin{equation}\frac{{\partial \bar{S}}}{{\partial t}} + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\bar{S} + \frac{w}{{{R_\rho }}} = \boldsymbol{\nabla }\boldsymbol{\cdot }(K\,\boldsymbol{\nabla }\bar{S}) + \frac{\partial }{{\partial z}}\left( {\frac{K}{{{R_\rho }}}} \right) + \tau {\nabla ^2}\bar{S},\end{equation}

where $\bar{S}$ is the filtered salinity field, and diffusivity K is proportional to the large-scale strain. This subgrid model is analogous to the widely used Smagorinsky (Reference Smagorinsky1963) scheme. The key difference is that the original Smagorinsky closure represents the subgrid transfer of momentum. We, on the other hand, are interested solely in modelling the sub-microscale dynamics of salinity.

Closure (4.1) exhibited a reasonably rapid convergence to the corresponding DNS-based solutions with increasing resolution. However, it was subsequently noticed that its bi-harmonic counterpart

(4.2)\begin{equation}\frac{{\partial \bar{S}}}{{\partial t}} + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\bar{S} + \frac{w}{{{R_\rho }}} =- \boldsymbol{\nabla }\boldsymbol{\cdot }({K_2}\,\boldsymbol{\nabla }{\bar{S}_2}) + \tau {\bar{S}_2},\quad {\bar{S}_2} = {\nabla ^2}\bar{S}\end{equation}

offers a consistently more accurate representation of key mixing characteristics. The coefficient ${K_2}$ in (4.2) is parameterized as

(4.3)\begin{equation}{K_2} = C{\delta ^4}s,\end{equation}

where $\delta = {L_x}/{N_x} = {L_y}/{N_y} = {L_z}/{N_z}$ is the grid spacing, and C is an adjustable parameter. The superior performance of the bi-harmonic Smagorinsky closure is attributed to its selective suppression of sub-microscale components. The bi-harmonic forcing term in (4.2) has a minimal direct impact on finger-scale components and therefore the model is less invasive and, consequently, more precise. The impact of closure (4.2) on sub-microscale salinity patterns is explored in greater detail in Appendix A, where we use theoretical arguments to rationalize its ability to effectively control grid-scale processes.

An interesting question concerns the choice of the coefficient C in (4.3). A straightforward option – and the one utilized in the Smagorinsky (Reference Smagorinsky1963) closure – is a selection of a constant numerical value for C that is applied to all systems considered. More involved are the dynamic models (e.g. Khani & Waite Reference Khani and Waite2015) in which the calibration of subgrid schemes is based on the properties of a particular flow. In the present study, we consider both constant-C and adaptive-C models (§ 5) and find that their accuracy and convergence rates are generally comparable. The selection of the optimal value of C in the adaptive-C closure is based on the requirement to fully resolve the dissipation spectrum of the filtered salinity variance $\bar{D}(\kappa )$. Specifically, we insist that the maximal dissipation ${D_{max }} = \textrm{ma}{\textrm{x}_\kappa }\bar{D}\textrm{(}\kappa \textrm{)}$ significantly exceeds the high-wavenumber dissipation ${D_{high}} = \textrm{ma}{\textrm{x}_{\kappa > 0.5{\kappa _0}}}\bar{D}(\kappa )$, where ${\kappa _0}$ is the largest resolved wavenumber. At the same time, we wish to avoid excessive damping of the salinity field. To achieve a reasonable balance between these competing requirements, we assume the target ratio of $D{}_{max}/{D_{high}} = 30$. In the course of simulations, C is gradually increased (decreased) when the ratio $D{}_{max}/{D_{high}}$ drifts below (above) its target value.

5. Validation: 2-D simulations

All subsequent SMF runs are based on the bi-harmonic Smagorinsky model (4.2) of the salinity field. We start by analysing its performance in two dimensions. An obvious benefit of 2-D simulations is their efficiency, which permits a more systematic exploration of the parameter space. It should also be noted that externally induced large-scale flows, which are ubiquitous in the ocean, favour the formation of salt sheets aligned in the direction of the background shear (Linden Reference Linden1974; Kimura & Smyth Reference Kimura and Smyth2007; Radko et al. Reference Radko, Ball, Colosi and Flanagan2015). In such cases, salt fingers become effectively two-dimensional. This tendency implies that analyses based on 2-D models may be more oceanographically relevant.

To illustrate the efficacy of the proposed subgrid closure, we now reproduce the earlier simulations (§ 3) with the SMF model. Equation (4.2) is used to integrate salinity, the original equations in (2.5) are used for temperature and vorticity, and the previously employed spectral code is altered only by including the subgrid forcing term. The key benefit of the SMF model is its ability to perform integrations on much coarser grids. For instance, figure 5 presents the analogue of the simulation in figure 1 performed with only $({N_x},{N_z}) = (128,128)$ grid points using the adaptive-C model. This mesh can barely represent the temperature and momentum fields, and most of the original salinity dissipation spectrum (figure 3) is unresolved. Nevertheless, the SMF model offers the description of fingering convection (figure 5) that is consistent with the corresponding fully resolved simulation (cf. figure 1). By $t = 20$, vertically elongated patterns form in both temperature and salinity (figure 5a,b). These primary salt fingers gradually amplify and, by $t = 50$ (figure 5c,d), develop secondary instabilities, which lead to the statistical equilibration of the flow field. After $t = 100$ (figure 5e,f), the system maintains a quasi-steady state.

Figure 5. The SMF analogue of the simulation in figure 1. The instantaneous patterns of temperature (left panels) and filtered salinity (right panels) are shown for (a,b) $t = 20$, (c,d) $t = 50$ and (e,f) $t = 100$.

A more quantitative assessment of the SMF model is afforded by the analysis of kinetic energy (E), as well as turbulent fluxes of heat (${F_T}$), salt (${F_S}$) and density (${F_\rho }$):

(4.4a–d)\begin{equation}{F_T} = \langle wT\rangle ,\quad {F_S} = \langle wS\rangle ,\quad {F_\rho } = \langle w(S - T)\rangle ,\quad E = {\textstyle{1 \over 2}}\langle {v_i}{v_i}\rangle .\end{equation}

The corresponding expression of the salt flux in terms of filtered salinity is given by

(4.5)\begin{equation}{F_S} = \langle w\bar{S}\rangle - \left\langle {{K_2}\frac{\partial }{{\partial z}}{\nabla^2}\bar{S}} \right\rangle .\end{equation}

The first term on the right-hand side of (4.5) represents the resolved salt flux, which greatly exceeds the second (parameterized) component in all SMF simulations. The corresponding density flux ${F_\rho } = {F_S} - {F_T}$ is evaluated using (4.5). The prediction of quantities (4.4) is one of the key objectives of the double-diffusive convection theory. Therefore, the utility of any reduced-dynamics double-diffusive model is ultimately determined by its ability to accurately evaluate (4.4).

To that end, in figure 6 we plot the mean values of $({F_T},{F_S},{F_\rho },E)$ for a series of SMF runs performed with $(Pr,\tau ,{R_\rho }) = (7,0.01,2)$ in which we systematically vary the resolution ($\delta $). These simulations were performed with $N \equiv {N_x} = {N_z} = 64,$ 96, 128 192, 256, 384, 512, 768 and 1024. We average the diagnostic quantities (4.4) in time over the quasi-equilibrium interval $100 < t < 1000\,.$ Both adaptive-C and constant-C models are presented, along with the estimates based on the fully resolved DNS in figure 1. The adaptive-C experiments indicate the optimal values of the Smagorinsky coefficient do not exceed C = 0.11 regardless of the resolution. This finding guides the calibration of the fixed-C model, and C = 0.1 is used for the latter. The adaptive-C and fixed-C experiments are mutually consistent, revealing the monotonic convergence of quantities (4.4) to their DNS-based counterparts with decreasing $\delta $.

Figure 6. Convergence analysis of the SMF model. The magnitudes of heat flux (a), salt flux (b), density flux (c) and energy (d) are plotted as functions of the resolution $(\delta )$. The simulations performed with the adaptive-C model are shown in blue and the constant-C model (C = 0.1) are indicated by the red curves. The dashed lines represent fully resolved DNS. The error bars reflect the uncertainties in the estimates of averages associated with the finite length of records for the SMF model. The shaded regions represent the analogous uncertainties for the DNS model.

While the SMF model substantially underestimates the integral quantities (4.4) for $N = 64$ and $N = 96$, its accuracy improves considerably in better-resolved runs. For instance, the differences between the SMF-based and DNS-based estimates of the salinity (temperature) flux reduce to less than 15% (5%) for N ≥ 128 $(\delta \le 0.39)$. It should be noted that all quantities in (4.4) are characterized by substantial temporal variability. The uncertainty of our estimates of the temporal averages that can be attributed to this variability is assessed as follows. We compute a series of averages over smaller intervals $(\varDelta {t_{av}} = 450)$ within the equilibrium range $100 < t < 1000$ and then use the root-mean-square variability of the resulting estimates as a measure of error. The length-of-record uncertainty, which is also indicated in figure 6, turns out to be particularly large for the temperature fluxes. The estimates of ${F_T}$ offered by the SMF model for N ≥ 256 $(\delta \le 0.195)$ are statistically undistinguishable from their DNS-based counterparts.

Figure 7 presents the spectra of kinetic energy $\hat{E}(\kappa )$, where

(4.6)\begin{equation}E = \int {\hat{E}(\kappa )\,\textrm{d}\kappa } .\end{equation}

The energy spectra are shown for the fully resolved simulation (figure 1) and its SMF counterpart performed with N = 128 (figure 5). The two spectra are mutually consistent for relatively low wavenumbers $(\kappa \lesssim 3)$ but the SMF model appears to be more aggressive in suppressing high-wavenumber harmonics. It should be emphasized, however, that both models reveal the lack of any substantial velocity patterns on scales that are much less than the typical size of salt fingers $(\kappa \gtrsim 10)\,.$ These diagnostics support the principal tenet of the SMF model, which assumes that sub-microscale dynamics of velocity and temperature is inconsequential.

Figure 7. The spectrum of kinetic energy for the fully resolved DNS (solid curve) and its SMF counterpart (dashed curve).

Figure 8 presents the convergence analysis for the low-diffusivity case with $\tau = 5 \times {10^{ - 4}}$. Once again, we observe the rapid reduction in the model error with decreasing $\delta $. For $N = 196$, the SMF experiments underestimate the vertical T and S transport by less than 2% and 12%, respectively. The errors of such magnitude can be tolerated in most double-diffusive applications, particularly considering the dramatic gain in efficiency relative to the corresponding fully resolved DNS – recall, for instance, that the experiment in figure 2(c,d) was performed with $N = 4096$. It is also possible that the relatively poor performance of the SMF model for N = 64 and N = 96 is partially caused by the incompletely resolved temperature and momentum fields, rather than by the deficiencies of the subgrid closure itself.

Figure 8. The same as in figure 6 but for the low-diffusivity SMF experiments performed with $\tau = 5 \times {10^{ - 4}}$. All simulations employ the adaptive-C closure.

6. Validation: 3-D simulations

Our next step is the validation of the SMF model in three dimensions. It is well known that cross-scale energy transport is fundamentally different in 2-D and 3-D turbulence (e.g. Kraichnan Reference Kraichnan1967). Therefore, despite its impressive performance in two dimensions, the ability of the SMF model to represent 3-D systems should not be taken for granted.

To establish a baseline for the validation, we perform a set of 3-D DNS. The size of the computational domain used in these simulations is ${L_x} \times {L_y} \times {L_z} = 30 \times 30 \times 30$, and it is resolved by $({N_x},{N_y},{N_z}) = (800,800,800)$ grid points. Our first example (figures 9 and 10) is a DNS performed with $(Pr,\tau ,{R_\rho }) = (7,0.01,2)$. Figure 9 presents a typical salinity field realized in the quasi-equilibrated evolutionary stage. In agreement with earlier studies (e.g. Radko & Smith Reference Radko and Smith2012; Radko et al. 2015), 3-D salt fingers (figure 9) proved to be structurally similar to – but more intense than – their 2-D counterparts (figure 1). In figure 10, we plot horizontal and vertical cross-sections of temperature (a,c) and salinity (b,d) for the state in figure 9 taken across the centre of the computational domain. The (x,z) sections are characterized by the emergence of vertically elongated filaments, whereas the (x,y) patterns are statistically isotropic. Particularly relevant for our discussion is the abundance of sub-microscale features in both horizontal and vertical sections of salinity – features that are not reflected in the corresponding temperature fields.

Figure 9. The 3-D DNS performed with $(Pr,\tau ,{R_\rho }) = (7,0.01,2)$. Presented is the instantaneous pattern of salinity at t = 100.

Figure 10. Horizontal (a,b) and vertical (c,d) sections of temperature (a,c) and salinity (b,d) for the state in figure 9.

To determine whether the DNS results (figures 9 and 10) can be adequately reproduced using the SMF model, we turn to simulations based on (4.2) with the adaptive-C closure. The mesh used for these runs is systematically refined by increasing the number of grid points in each direction: $N \equiv {N_x} = {N_y} = {N_z} = 64$, 96, 128 192, 256 and 384. Other parameters match those of the simulation in figure 9. Figures 11 and 12 present the SMF experiment performed with N = 128. Despite the dramatic difference in resolution, the finger-scale processes in figures 9 and 11 are structurally similar. The fully resolved DNS (figure 9) exhibits a wider range of salinity scales, but the typical magnitudes of salinity perturbations in the two simulations are consistent. The comparison of the T–S cross-sections in figures 10 and 12 carries the same message. The SMF model filters the sub-microscale components of salinity in a way that does not significantly alter the dynamics of larger scales of motion.

Figure 11. The SMF experiment performed with N = 128 and $(Pr,\tau ,{R_\rho }) = (7,0.01,2)$. Presented is an instantaneous pattern of salinity at t = 100. Note the qualitative similarity of the flow field to its fully resolved counterpart in figure 9.

Figure 12. Horizontal (a,b) and vertical (c,d) sections of temperature (a,c) and salinity (b,d) for the state in figure 11.

To be more quantitative in assessing the performance of the SMF model, we now examine its convergence. In figure 13, we plot the quasi-equilibrium values of $({F_T},{F_S},{F_\rho },E)$ as functions of the resolution $\textrm{(}\delta \textrm{)}$ along with the corresponding prediction based on the fully resolved simulation (figures 9 and 10). These diagnostics reveal an impressive skill of the SMF model in predicting key integral characteristics of the flow field. In all cases, the model error rapidly decreases with increasing resolution. Particularly impressive is the model performance for very low N. The errors of the predicted salt and temperature fluxes for the smallest grid with N = 64 $\textrm{(}\delta = 0.469\textrm{)}$ are only 4% and 20%, respectively. It is interesting to note that the estimates of salt fluxes by the SMF model are consistently better in three dimensions (figure 13) than in two dimensions (figure 6). For instance, the accuracy of 4% for F_S in two dimensions is attained only for $\delta \le 0.2\ \textrm{(}N \ge 256\textrm{)}$. The superior performance of the SMF model in three dimensions can be attributed to the downscale (upscale) transfer of energy in 3-D (2-D) turbulence (Kraichnan Reference Kraichnan1967). The inverse 2-D cascade implies that inaccuracies in the representation of sub-microscale salinity patterns could spread upscale and contaminate larger scales of motion. On the other hand, the direct 3-D cascade effectively constrains the model errors to the high-wavenumber part of the spectrum.

Figure 13. Convergence analysis of the 3-D SMF model performed with ${R_\rho } = 2$. The magnitudes of heat flux (a), salt flux (b), density flux (c) and energy (d) are plotted as functions of the resolution $(\delta )$.

Finally, to explore the relevant parameter space, we perform two additional convergence studies for ${R_\rho } = 1.5$ and ${R_\rho } = 3$, while keeping $(Pr,\tau ) = (7,0.01)$. The density ratio of ${R_\rho } = 1.5$ leads to some of the more intense salt fingering in the ocean, such as realized in the Caribbean Sea, whereas ${R_\rho } = 3$ represents its relatively mild form, more common in the Pacific thermocline. The results (figures 14 and 15) are consistent with our earlier findings (cf. figure 13). The key flow properties (4.4) systematically approach their DNS-based counterparts when $\delta $ is decreased. In the most violent regime (figure 14), the SMF model is slightly less accurate than in its less energetic counterparts (figures 13 and 15), particularly in terms of predicting the density flux. However, in all cases, adequate estimates are obtained even when the number of points in each dimension is less by an order of magnitude than demanded by fully resolved DNS. Considering the associated increase in the time step controlled by the CFL (Courant–Friedrichs–Lewy) condition, the computational savings brought by the 3-D SMF model are of the order of 10⁴.

Figure 14. The same as in figure 13 but for the experiments performed with ${R_\rho } = 1.5$.

Figure 15. The same as in figures 13 and 14 but for the experiments performed with ${R_\rho } = 3$.

A pragmatic question can be raised at this point regarding resolution guidelines for the SMF model. Based on all simulations in this study, two-dimensional and three-dimensional, we suggest that the grid spacing of ${\delta _{cr}} = 0.2$ is sufficient to capture the key processes at play. This criterion is equivalent to resolving the nominal Stern (Reference Stern1960) finger scale $l = {({k_T}\upsilon /g\alpha {A_T})^{1/4}}$ by at least five grid points. Admittedly, our ad hoc recommendation is based on the Fourier spectral model and may require modification for other numerical algorithms.