The sub-microscale dynamics of double-diffusive convection

Abstract

This study investigates the dynamics of fingering convection on scales much smaller than the typical size of individual salt fingers. On such scales, salinity patterns exhibit the spontaneous emergence of sharp fronts induced by finger-scale strain. In contrast, velocity and temperature fields are largely devoid of sub-microscale variability, which is attributed to the rapid molecular dissipation of heat and momentum. The presence of fine salinity structures fundamentally limits the efficiency of direct numerical simulations (DNS) of double-diffusive processes. In the oceanographic context, the computational cost of resolving sub-microscale salinity features exceeds that of temperature-only DNS by up to four orders of magnitude, severely restricting the types of double-diffusive systems that can be studied numerically. To address this complication, we introduce the sub-microscale filtering (SMF) algorithm, which resolves temperature and velocity while parameterizing the sub-microscale dynamics of salinity. The proposed closure draws inspiration from the Smagorinsky scheme, which represents unresolved processes by the downgradient strain-dependent momentum flux. The SMF model is successfully validated through fully resolved simulations.

1. Introduction

The term double-diffusive convection represents a wide range of processes that are driven by unequal diffusivities of properties affecting fluid density. Following the seminal discovery of primary double-diffusive instabilities (Stern 1960), investigations of multi-component flows evolved into a distinct and vibrant branch of fluid mechanics. Double-diffusive phenomena have been studied most extensively in the oceanographic context (Schmitt 1994; Kelley et al. 2003; Radko 2013), where the density of seawater is controlled by its temperature and salinity. Since heat diffuses approximately a hundred times faster than salt, double-diffusive processes in the ocean tend to be intense and widespread (You 2002). In several critical regions, they dominate the vertical mixing of water masses (Schmitt et al. 2005; Guthrie, Fer, & Morison 2015), motivating in-depth studies of double-diffusive dynamics and transport. The fingering mode of double-diffusive convection is particularly common in the main thermocline, a stratified region extending from the base of the surface mixed layer to depths of approximately one kilometre. In the Atlantic Ocean, at least 90% of the main thermocline is fingering favourable (Schmitt 1994). In addition to oceanic systems, double-diffusive convection is also realized – and plays a significant role – in numerous astrophysical, geophysical and engineering applications (e.g. Turner 1985; Radko 2013; Garaud 2018).

In the spirit of the times, much recent progress in understanding double-diffusive convection has been brought by numerical modelling (e.g. Smyth & Kimura 2011; Hieronymus & Carpenter 2016; Yang, Verzicco & Lohse 2016, 2020; Ouillon et al. 2020; Ma & Peltier 2021; Brown & Radko 2022). However, the principal limitation of such efforts lies in the prohibitive computational cost of modelling many realistic configurations. Double-diffusive processes are initiated by molecular diffusion and therefore primary instabilities in seawater operate on spatial scales of centimetres, placing them in the category of oceanic microstructure. While interesting in their own right, an even broader geophysical significance of primary double-diffusive instabilities lies in their ability to induce much larger flows. This dynamics is exemplified by salt fingers, the vertically elongated microscale filaments that form in the ocean regions where temperature and salinity increase upward. Salt fingers frequently induce secondary circulation patterns, which include lateral intrusions (Stern 1967; Ruddick & Kerr 2003; Ruddick & Richards 2003), collective instability waves (Stern 1969; Stern, Radko & Simeonov 2001) and thermohaline staircases (Stern & Turner 1969; Stellmach et al. 2011). The spatial extents of these secondary double-diffusive patterns exceed the scale of individual fingers by at least 2–3 orders of magnitude, which presents a fundamental obstacle for numerical models. Despite continuous advancements in high-performance computing, microstructure-resolving simulations of double-diffusive processes are still extremely demanding. Even with a major investment of computational resources, users are often forced to make painful concessions in terms of the realism of model configurations and chosen parameters. Several models attempted to parameterize, rather than resolve, double-diffusive microstructure (e.g. Stern & Simeonov 2002; Simeonov & Stern 2004; Radko & Sisti 2017). The latter approach, however, suffers from uncertainties in the formulation of such parameterizations and is expected to provide only qualitative descriptions of the processes at play.

Yet another modelling challenge – and the one that the present study strives to address – stems from the difference in the molecular diffusivities of heat and salt, ${k_T}$ and ${k_S}$, respectively. The dissipation scale of salinity is less than the dissipation scale of temperature by a factor of $\sqrt {{k_T}k_S^{ - 1}} \sim 10$ (e.g. Radko 2008), which conforms to the Batchelor (1959) scaling of scalar dissipation. Thus, the requirement to resolve sub-microscale salinity features increases the computational cost by a factor of 10³ in two-dimensional (2-D) simulations and 10⁴ in three dimensions relative to the investment in equivalent temperature-only direct numerical simulations (DNS). Unfortunately, this complication makes several fundamental double-diffusive problems numerically inaccessible even in two dimensions. Note that the difference in diffusivities captures the very essence of double diffusion and therefore it must be reflected in any meaningful simulations. Many earlier numerical studies reduced the associated computational burden by considering the salt/heat diffusivity ratios $\tau \equiv {k_S}k_T^{ - 1}$ that are much larger than the oceanic value of $\tau \sim 0.01$ (e.g. Stern et al. 2001; Kimura, Smyth, & Kunze 2011; Brown & Radko 2022). While acceptable for preliminary explorations, concerns regarding the ability of such models to accurately represent oceanic processes call for a more cardinal solution to the problem. Another approach involves using a much finer mesh for salinity than for temperature and velocity (e.g. Kimura et al 2011; Smyth & Kimura 2011). However, since the fine mesh must still resolve the salinity dissipation scale, this technique reduces the computational cost by less than a factor of five. Thus, while the advent of double-grid DNS represents a welcome addition to the modelling toolbox, many double-diffusive problems remain numerically inaccessible.

This study addresses the twin challenges of explaining the sub-microscale dynamics of double-diffusive convection and offering a physics-based representation of key processes that can be implemented in coarse-resolution numerical models. The formation of the sub-microscale salinity patterns is attributed to the strain by the finger-scale velocity field, which stretches the salinity filaments and systematically reduces their effective width. The resulting direct cascade of salinity variance is arrested only at scales that are small enough to be effectively damped by the molecular dissipation of salt. This chain of events is consistent with the classical description of the small-scale behaviour of weakly diffusive scalars in the presence of turbulence by Batchelor (1959). To capture this dynamics, we develop a strain-dependent subgrid salinity model. This closure is inspired by the pioneering large eddy simulation (LES) approach of Smagorinsky (1963), where an analogous scheme was implemented to parameterize the subgrid momentum transfer. In contrast, the present algorithm targets salinity and makes no attempt to parameterize momentum and temperature, which are fully resolved. This sub-microscale filtering (SMF) model is used to simulate fingering convection in the oceanographically relevant regime and is validated by the corresponding salinity-resolving experiments.

The material is organized as follows. The governing equations are described in § 2. Section 3 presents the physical arguments and supporting simulations that establish the link between the sub-microscale salinity dissipation and the finger-scale velocity strain. Section 4 introduces the strain-based subgrid model of salinity. We validate the SMF model using fully resolved 2-D and 3-D simulations in §§ 5 and 6, respectively. The results are summarized, and conclusions are drawn, in § 7.

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 2013) 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. 2001; Radko 2019a) 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 2024) 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 2013) 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 2012) 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 (1959), 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 1972) 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 2008), 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 (1959) 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 1980; Oakey 1982) 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 (1959) 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 (1974) 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 (1963) 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 (1963) 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 2015) 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 1974; Kimura & Smyth 2007; Radko et al. 2015). 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 1967). 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 2012; 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 1967). 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 (1960) 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.

7. Discussion

This study explores the sub-microscale dynamics of fingering convection. We are interested in flow features with scales that are much less than the typical size of individual salt fingers. On such scales, the slower diffusing property (salinity in the oceanographic context) develops narrow fronts that play a prominent role in the dissipation of the salinity variance. The emergence of these salinity structures stands in stark contrast with the lack of sub-microscale patterns of velocity and the faster diffusing density component (temperature in the ocean), which is attributed to the rapid molecular dissipation of heat and momentum. The very dissimilar evolution of salinity is not surprising, given that molecular viscosity and temperature diffusivity exceed the diffusivity of salt by two and three orders of magnitude, respectively. It is shown that the formation of sub-microscale salinity patterns is controlled by the finger-scale strain. Both shear and normal components of strain tend to systematically stretch and tighten salinity fronts until they are narrow enough to be affected by weak molecular dissipation of salt.

The spontaneous development of sub-microscale salinity patterns presents a major obstacle to numerical modelling of double-diffusive phenomena. The failure to represent these structures has catastrophic consequences for the stability and accuracy of models, rendering them largely unusable. At the same time, the requirement to fully resolve all relevant scales makes double-diffusive DNS prohibitively expensive. The computational cost can exceed that of the equivalent temperature-only simulations by three or four orders of magnitude, placing several double-diffusive problems beyond the reach of even the most advanced supercomputers.

Motivated by the principal limitation of the traditional DNS as a modelling tool for many double-diffusive problems, we develop a model of two-component flows that fully resolves velocity and temperature but parameterizes the subgrid transfer of salinity. Based on the analysis of sub-microscale dynamics (§ 3), we consider the strain-dependent subgrid closure for salinity. The proposed algorithm is reminiscent of the treatment of the momentum transfer in the pioneering model of Smagorinsky (1963). The convergence analysis indicates the proposed model (SMF) performs well for fingering convection in the oceanographically relevant regime. It produces representative flow patterns and accurate estimates of mixing even when the number of points in each dimension is less by an order of magnitude than required by DNS. It is interesting to note that the attempts (not shown) to implement explicit spectral cutoff filters for salinity have led to highly erratic model behaviour and were eventually abandoned. Such a contrast in terms of reliability and fidelity with the presented model is attributed to its robust, physics-based foundation. Of particular significance is the link between the salinity dissipation and the finger-scale strain, exploited by the proposed algorithm. The connection between strain and dissipation is the crux of the classical theories of scalar transfer in turbulent systems (Batchelor 1959; Kraichnan 1974) and is shown here to control the dynamics of double-diffusive convection as well.

While the developments reported here appear to be highly promising, they should be viewed just as a first step in the search for the optimal SMF scheme. The first LES turbulence models inspired a multitude of enhancements in several branches of computational fluid dynamics – the efforts reviewed recently by Moser, Haering & Yalla (2021). By the same token, we hope that the presented algorithm will place the double-diffusive community on the path toward continuously evolving SMF modelling. The essential physics of sub-microscale processes and traditional turbulence problems are fundamentally different in many ways. As a result, the design of SMF models must be guided by principles that are not reflected in conventional LES algorithms. Nevertheless, such efforts will eventually make it possible to address grand-challenge problems in our field, such as the realistic modelling of thermohaline staircases and double-diffusive interleaving. The natural next step should involve the application of SMF models to diffusive convection problems, where cold and fresh water masses overly the relatively warm and salty ones.

Another major item on our wish list is the expansion of this analysis beyond the heat–salt systems. In this regard, we note that some double-diffusive applications are even more challenging computationally than the oceanic case (Schmitt 1983; Radko 2013). For instance, the diffusivity of silicate melt components can be extremely low ($\tau \sim {10^{ - 3}} \unicode{x2013} {10^{ - 8}}$) and so is the compositional diffusivity in stars ($\tau \sim {10^{ - 6}} \unicode{x2013} {10^{ - 8}}$). As a result, the disparity of dissipation scales of density components is dramatic, and the fully resolved DNS are not only computationally expensive but currently impossible. For such systems, SMF modelling may be the only resort in the foreseeable future.