2.1. The unbounded gradient model and the interface model

Depending on the different initial and boundary conditions in simulations of the salt-fingering system, the evolution of the fingering field can vary considerably. In fact, three different representations (the unbounded gradient model, the interface model and the bounded gradient model) have been employed by different researchers to discuss different aspects of salt-fingering instabilities and these have been summarized in Radko (Reference Radko2013). The unbounded gradient model and interface model, which we will employ herein to calibrate and test the turbulent flux laws, will be briefly reviewed in this subsection in order to fix ideas.

The nonlinear hydrodynamic field equations for the heat-salt doubly diffusive system are

(2.1a)\begin{gather} \frac{\textrm{D} u}{\textrm{D}t} = -\frac{1}{\rho_0}\frac{\partial p'}{\partial x} + \nu\nabla^2u, \end{gather}

(2.1b)\begin{gather}\frac{\textrm{D} v}{\textrm{D}t} = -\frac{1}{\rho_0}\frac{\partial p'}{\partial y} + \nu\nabla^2v, \end{gather}

(2.1c)\begin{gather}\frac{\textrm{D} w}{\textrm{D}t}= -\frac{1}{\rho_0}\frac{\partial p'}{\partial z} - g\frac{\rho-\rho_0}{\rho_0}+\nu\nabla^2w, \end{gather}

(2.1d)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

(2.1e)\begin{gather}\frac{\textrm{D} \varTheta}{\textrm{D}t} = \kappa_\theta\nabla^2\varTheta, \end{gather}

(2.1f)\begin{gather}\frac{\textrm{D} S}{\textrm{D}t} = \kappa_s\nabla^2S, \end{gather}

(2.1g)\begin{gather}\rho = \rho_0(1 -\alpha (\varTheta-\varTheta_0) + \beta (S-S_0)). \end{gather}

In these field equations the velocity vector field is $\boldsymbol {u}(x,y,z,t)=\{u,v,w\}$, the potential temperature field is $\varTheta (x,y,z,t)$ and the salinity field is $S(x,y,z,t)$. In the linear equation of state (2.1g), the thermal expansion rate $\alpha$ and salinity contraction coefficient $\beta$ are both assumed to be constant. In this system, $p'(x,y,z,t)$ is the deviation of the pressure from a background state of hydrostatic balance, and $\nu$,$\kappa _\theta$ and $\kappa _s$ are kinematic viscosity and molecular diffusivity for heat and for salt, respectively, in the advection–diffusion equations ((2.1e) and (2.1f)).

Equation (2.1) are the governing equations for both the unbounded gradient model and the interface model. The major distinction between these models concerns the different initial and boundary conditions to be applied. The unbounded gradient model is depicting salt-fingering fields in a homogeneous environment with constant background temperature gradient $\varTheta _{z0}$ and salinity gradient $S_{z0}$, and so the initial conditions are

(2.2a)\begin{gather} \varTheta(x,y,z)= \varTheta_{z0}z, \end{gather}

(2.2b)\begin{gather}S(x,y,z)= S_{z0}z, \end{gather}

in which $\varTheta_{z0} = {\rm \Delta} \varTheta /L$ and $S_{z0} = {\rm \Delta} S / L$, as illustrated in figure 1(a). The boundary conditions for velocity as well as for the perturbation fields $\varTheta '$ and $S'$ (defined as deviations from the horizontally uniform-gradient profiles) are periodic in all three coordinate directions. Specifically, the vertical boundary conditions for $\varTheta$ and S are

(2.3a)\begin{gather} \varTheta(x,y,z=L,t)= \varTheta(x,y,z=0,t)+{\rm \Delta} \varTheta, \end{gather}

(2.3b)\begin{gather}S(x,y,z=L,t)= S(x,y,z=0,t)+{\rm \Delta} S. \end{gather}

By applying these vertical boundary conditions, we may treat the system as effectively consisting of an infinitely large domain comprised of an infinite number of such sub-systems with the average temperature (salinity) gradient fixed to the values of ${\rm \Delta} \varTheta /L$ (${\rm \Delta} S/L)$. Since the mean scalar gradient is everywhere the same, the unbounded gradient model is a suitable basis for analysis of equilibrium properties of salt-fingering fields. However, the special boundary condition (2.3) is somewhat artificial since it prevents the average background gradient from being influenced by the fluxes. This effect, on the other hand, can be clearly captured in the interface model that we proceed to discuss in what follows.

Figure 1. Illustration of the settings for the unbounded gradient model (a) and the interface model (b). In both parts of the figure, we use a solid line to represent the initial profile for temperature and the dashed line to represent the initial profile for salinity. Here ${\rm \Delta} \varTheta$, ${\rm \Delta} S$ and $L$ are the characteristic scales for potential temperature and salinity and length for each model that are used to non-dimensionalize each system.

Although there are a variety of different possible definitions for the interface model, for present purposes, we choose to assume hyperbolic tangent forms for the initial conditions for both temperature and salinity fields as in figure 1(b), namely:

(2.4a)\begin{gather} \varTheta(x,y,z)= {\rm \Delta} \varTheta \tanh\left(\frac{z}{L}\right), \end{gather}

(2.4b)\begin{gather}S(x,y,z)= {\rm \Delta} S \tanh\left(\frac{z}{L}\right). \end{gather}

In the formulation of this model we also assume periodic boundary conditions in the horizontal coordinates for both of the scalar fields ($\varTheta$, $S$) and velocity vector field $\boldsymbol {u}$. In the vertical coordinate direction we will assume free-slip, no-flux boundary conditions for both top and bottom boundaries of the domain.

In figure 1 we have also shown the length scale $L$, temperature scale ${\rm \Delta} \varTheta$ and salinity scale ${\rm \Delta} S$ that we will use to non-dimensionalize the field equations describing this system: in the unbounded gradient model, ${\rm \Delta} \varTheta$ and ${\rm \Delta} S$ are the corresponding scalar variations across the domain and $L$ is the vertical height of the domain; whereas in the interface model, ${\rm \Delta} \varTheta$ and ${\rm \Delta} S$ are half the variations of the temperature and salinity fields across the interface and $L$ is the measure of interface depth so that ${\rm \Delta} \varTheta /L ({\rm \Delta} S/L)$ represent the temperature/salinity gradients at the midpoint of the interface.

The original Boussinesq equations (2.1) with the incompressibility condition (2.1d) can then be non-dimensionalized assuming length scale $L$, time scale $L^2/\kappa _\theta$ and scalar scales ${\rm \Delta} \varTheta$, ${\rm \Delta} S$ to obtain

(2.5a)\begin{gather} \frac{\textrm{D} u_\ast}{\textrm{D} t_\ast} = -\frac{\partial p_\ast}{\partial x_\ast} + Pr\nabla^2_\ast u_\ast, \end{gather}

(2.5b)\begin{gather}\frac{\textrm{D} v_\ast}{\textrm{D} t_\ast} = -\frac{\partial p_\ast}{\partial y_\ast} + Pr\nabla^2_\ast v_\ast, \end{gather}

(2.5c)\begin{gather}\frac{\textrm{D} w_\ast}{\textrm{D} t_\ast}= -\frac{\partial p_\ast}{\partial z_\ast} +Ra Pr \varTheta_\ast -Ra Pr S_\ast/R_\rho+ Pr\nabla_\ast^2w_\ast , \end{gather}

(2.5d)\begin{gather}\boldsymbol{\nabla_\ast}\boldsymbol{\cdot} \boldsymbol{u_\ast} = 0, \end{gather}

(2.5e)\begin{gather}\frac{\textrm{D} \varTheta_\ast}{\textrm{D} t_\ast} = \nabla^2_\ast \varTheta_\ast, \end{gather}

(2.5f)\begin{gather}\frac{\textrm{D} S_\ast}{\textrm{D} t_\ast} = \tau \nabla^2_\ast S_\ast, \end{gather}

in which we have employed asterisks to denote non-dimensional variables. The non-dimensional transformation for each variable as well as the definition of non-dimensional parameters that appear in this non-dimensional system are listed in table 1. The Prandtl number and the diffusivity ratio are fixed to typical oceanographic values ($Pr=7$ and $\tau =0.01$) in this paper. The thermal Rayleigh number is essentially controlling the size of the system relative to the scale of the finger width, so that the volume-averaged physical quantities of the system are insensitive to variation of $Ra$ as long as $Ra$ is high enough (no less than order $10^7$ to guarantee that the system contains sufficient fingers so that horizontal averages are meaningful). The dominant non-dimensional number of the system is then the bulk density ratio $R_\rho$, which measures the degree of compensation between temperature and salinity gradients in terms of their effects on density stratification. The system will be susceptible to salt-fingering instability when $R_\rho$ satisfies $1<R_\rho <1/\tau$ (for both models), as revealed by linear stability analysis (e.g. Walin Reference Walin1964; Baines & Gill Reference Baines and Gill1969). The values of $Ra$ and $R_\rho$ employed in our numerical simulations for both the unbounded gradient model and the interface model will be discussed in detail in § 3 of this paper.

Table 1. Non-dimensional transformation and definition of non-dimensional numbers.

2.2. Energy budget equations and irreversible fluxes

In a salt-fingering-favourable system the unstably stratified salinity field releases its associated potential energy to create the kinetic energy required for mixing. This leads to an ‘up-gradient’ density flux that might be seen as somewhat counterintuitive. In this subsection we will discuss the energy budgets in the salt-fingering system from a somewhat novel perspective that separates the potential energy into a linear combination of those associated with temperature and salinity separately (this can also be applied to the diffusive convection counterpart to the salt-fingering system that we will discuss elsewhere). Application of this methodology will enable us to discuss the ‘true’ irreversible mixing in each of these two fields.

We begin by defining the average potential energy $PE$ of the system in a certain volume $V$ as

(2.6)\begin{equation} PE=\langle \rho g z\rangle=\frac{1}{V}\int \rho g z \, \textrm{d} V, \end{equation}

where $\langle \cdot \rangle$ represents volume averages over the system of interest.

Due to our assumption of a linear equation of state (2.1g), we may further decompose the potential energy $PE$ into a temperature component $PE_\varTheta$ and a salinity component $PE_S$ as

(2.7a)\begin{align} PE &= \langle \rho g z \rangle, \end{align}

(2.7b)\begin{align}&= \langle \rho_0 g z \rangle-\langle \rho_0\alpha (\varTheta-\varTheta_0) g z \rangle+\langle \rho_0\beta(S-S_0) g z \rangle, \end{align}

(2.7c)\begin{align}&= \langle (1+\alpha\varTheta_0-\beta S_0) \rho_0 g z \rangle-\rho_0\alpha g \langle \varTheta z \rangle+ \rho_0 \beta g \langle S z \rangle, \end{align}

(2.7d)\begin{align}&\equiv PE_0+PE_\varTheta+PE_S, \end{align}

where $PE_\varTheta$ and $PE_S$ denote the potential energy contributed by the temperature and salinity stratification of the fluid, separately. Equation (2.7) shows that the total potential energy can be viewed as a simple linear combination of $PE_\varTheta$ and $PE_S$, a result that depends only upon the assumption of a linear equation of state. Here $PE_0$ is just a constant that may be ignored as it can play no role in the understanding of mixing processes.

The evolution of $PE$, $PE_\varTheta$, as well as $PE_S$ in a closed system subject to the Boussinesq approximation satisfy the simple evolution equations

(2.8a)\begin{gather} \frac{\textrm{d} PE}{\textrm{d} t}= \mathcal{D}_{p}+\mathcal{H}, \end{gather}

(2.8b)\begin{gather}\frac{\textrm{d} PE_\varTheta}{\textrm{d} t}= \mathcal{D}_{p\varTheta}+\mathcal{H}_\varTheta, \end{gather}

(2.8c)\begin{gather}\frac{\textrm{d} PE_S}{\textrm{d} t}= \mathcal{D}_{pS}+\mathcal{H}_S, \end{gather}

in which

(2.9a)\begin{gather} \mathcal{H}_\varTheta= -\rho_0 g \alpha \langle \varTheta'(x,y,z,t) w'(x,y,z,t) \rangle, \end{gather}

(2.9b)\begin{gather}\mathcal{H}_S= \rho_0 g \beta \langle S' (x,y,z,t)w'(x,y,z,t) \rangle, \end{gather}

(2.9c)\begin{gather}\mathcal{H}= \mathcal{H}_\varTheta+\mathcal{H}_S \end{gather}

(2.9d)\begin{gather}\qquad\qquad\qquad\qquad = g\langle \rho'(x,y,z,t)w'(x,y,z,t) \rangle, \end{gather}

(2.9e)\begin{gather}\mathcal{D}_{p\varTheta}= \rho_0g\alpha \kappa_\varTheta \left\langle \frac{\partial \bar{\varTheta}(z)}{\partial z} \right\rangle, \end{gather}

(2.9f)\begin{gather}\mathcal{D}_{pS}= -\rho_0g\beta\kappa_S \left\langle \frac{\partial \bar{S}(z)}{\partial z} \right\rangle, \end{gather}

(2.9g)\begin{gather}\mathcal{D}_p= \mathcal{D}_{p\varTheta}+\mathcal{D}_{pS}. \end{gather}

In the above equations the overbar $\bar {f}(z)$ represents the horizontal average of field $f$ and $f'$ represents the deviation from it, i.e. $f'(x,y,z,t)=f(x,y,z,t)-\bar {f}(z)$. Here $\mathcal {H}$ represents the total buoyancy flux which is comprised of contributions from heat $\mathcal {H}_\varTheta$ and salt $\mathcal {H}_S$; $\mathcal {D}_p$ represents the irreversible energy transfer associated with molecular diffusivity of both diffusing species ($\mathcal {D}_{p\varTheta }$ and $\mathcal {D}_{pS}$) in the motionless background that cannot be influenced by the turbulent motions.

In a system that supports salt-fingering instability, there is warm salty water lying above cold fresh water, leading upon the development of instability involving both heat flux $F_\varTheta \equiv \langle w'\varTheta ' \rangle$ and salt flux $F_S\equiv \langle w'S'\rangle$ directed downwards, and, thus, $\mathcal {H}_\varTheta >0$,$\mathcal {H}_S<0$ from ((2.9b), (2.9c)). Meanwhile, positive vertical gradients for both $\varTheta$ and $S$ will lead to $\mathcal {D}_{p\varTheta }>0$ and $\mathcal {D}_{pS}<0$. This suggests a continuously increasing temperature potential energy $\textrm {d} PE_\varTheta / \textrm {d} t>0$ as well as a continuously decreasing salinity potential energy $\textrm {d}PE_S/\textrm {d}t<0$ (from (2.8b) and (2.8c)), which implies that the salinity field is releasing potential energy whereas the temperature field is gaining potential energy in the salt-fingering field. The ratio of these buoyancy fluxes ($\gamma =|\mathcal {H}_\varTheta |/|\mathcal {H}_S|=\alpha \langle w'\varTheta ' \rangle /\beta \langle w'S' \rangle$) is always smaller than 1 due to more effective salt exchange than heat exchange in the fingering system (Radko Reference Radko2013). This fact ensures that $\mathcal {H}=\mathcal {H}_\varTheta +\mathcal {H}_S<0$ and, thus, the system is continuously releasing potential energy. These fluxes, following the initial onset of instability, eventually impact the average kinetic energy of the system which is defined as $KE\equiv \langle 1/2 \rho _0 \boldsymbol {u}^2 \rangle$, whose evolution in a closed system can be shown to be

(2.10)\begin{equation} \frac{\textrm{d} KE}{\textrm{d} t}=-\mathcal{H}-\epsilon, \end{equation}

in which the viscous dissipation $\epsilon$ is $2\rho _o \nu \langle e_{ij}e_{ij} \rangle$. These two terms, which are of opposite sign, balance each other in an equilibrium (statistical equilibrium) state of the system.

The above definitions for heat and salt derived buoyancy flux ($\mathcal {H}_\varTheta$, $\mathcal {H}_S$) and flux ratio $\gamma$ are the quantities that have been employed in the discussion of double diffusion processes (see, e.g. Shen Reference Shen1995; Traxler et al. Reference Traxler, Stellmach, Garaud, Radko and Brummell2011; Radko & Smith Reference Radko and Smith2012). However, these definitions clearly fail to distinguish the differences between reversible stirring and irreversible mixing, both of which will be contributing to the transfer of energy between $PE$ and $KE$. As in the shear-driven stratified turbulence case (see Winters et al. Reference Winters, Lombard, Riley and D'Asaro1995; Caulfield & Peltier Reference Caulfield and Peltier2000; Peltier & Caulfield Reference Peltier and Caulfield2003), stirring refers to the temporary exchange between $PE$ and $KE$ caused by large-scale fluid motion which is reversible in principle and cannot contribute to irreversible mixing which occurs at the molecular scale. The critical step in the recognition of whether irreversible mixing occurs is to monitor whether the background potential energy is evolving as the turbulence continues. As illustrated in figure 2 of Caulfield & Peltier (Reference Caulfield and Peltier2000) and figure 2 of Peltier & Caulfield (Reference Peltier and Caulfield2003), the background potential energy (also known as the minimum potential energy) is the potential energy that would be associated with the adiabatic redistribution of fluid parcels in such a way that the density $\rho (x,y,z,t)$ is monotonically decreasing with height and, therefore, statically stable. The definition of the background potential energy guarantees that it will not be influenced by the macroscopic fluid motions (the stirring processes), rather, the sorted fluid elements will only lead to an increase of background potential energy in time when molecular diffusion is irreversibly changing the densities of neighbouring fluid parcels.

In our analysis of the salt-fingering system we will employ similar ideas to distinguish between irreversible mixing and reversible stirring. However, as density is influenced by both temperature and salinity in salt-fingering convection, we will need to sort both scalar fields separately rather than sorting the density field. Otherwise, the temperature and salinity information would be lost in the sorting process and it would then be impossible to determine distinct heat and salinity fluxes. How this sorting procedure should be performed will depend strongly upon the particular flavour of doubly diffusive turbulence being analysed: in the salt-fingering problem, since the salinity field is unstably stratified which supplies energy to both the potential energy associated with the temperature field and the kinetic energy reservoir, we need to adiabatically rearrange the salinity profile into a state such that the salinity is monotonically increasing with height (corresponding to a top heavy structure). In this way, irreversible mixing in the salinity field will lead to a monotonic decrease of the background salinity related potential energy, which is consistent with the basic energy flows in the salt-fingering field. For the temperature field, this must also be adiabatically rearranged to a state in which temperature is monotonically increasing with height. (The sorting directions will be just opposite for both the temperature and the salinity field in a diffusive convection system, which we will address elsewhere.) Since temperature is negatively correlated with density, we essentially sort the temperature to a bottom heavy structure just as in the shear-driven turbulence case. Irreversible mixing of the temperature field will lead to an increase in the corresponding background temperature potential energy, which is then also consistent with the direction of the energy flow in the salt-fingering problem.

Explicitly, we firstly rearrange all the fluid elements adiabatically to a configuration in which warmer elements are always on top of the colder elements. In this way the temperature of the sorted configuration will only depend on the vertical coordinate, which can be denoted as $\varTheta _\star (z,t)$. Then we do the sorting of the salinity field to ensure that the saltier elements are always above the fresher elements, which will give us a sorted salinity configuration $S_\star (z,t)$.

The background potential energy for temperature $BPE_\varTheta$ and the background potential energy for salinity $BPE_S$ can then be defined based on these two profiles as

(2.11a)\begin{gather} BPE_\varTheta = -\rho_0\alpha g \langle \varTheta_\star(z) z\rangle, \end{gather}

(2.11b)\begin{gather}BPE_S = \rho_0\beta g \langle S_\star(z) z\rangle, \end{gather}

(2.11c)\begin{gather}BPE = BPE_\varTheta+BPE_S. \end{gather}

In all the possible rearrangements of the fluid elements, $\varTheta _\star (z)$ is the configuration which has the minimum $PE_\varTheta$ of all the configurations (thus, $BPE_\varTheta \leq PE_\varTheta$), and $S_\star (z)$ is the configuration which has the maximum $PE_S$ of all the configurations (thus, $BPE_S \geq PE_S$). As the temperature diffusion term $\kappa _\theta \nabla ^2 \varTheta$ in (2.1e) mixes the warmer water columns with the colder water columns, the total mass of the sorted profile $\varTheta _\star (z)$ stays the same (in a closed system) but the centre of gravity is lifted, leading to the increase of $BPE_\varTheta$. Meanwhile, macroscopic fluid motion (stirring) only changes the relative positions of fluid parcels, which will not influence the configurations after sorting and, thus, will not lead to changes of $BPE_\varTheta$. Thus, the rate of change of $BPE_\varTheta$ will reflect the strength of irreversible mixing in the temperature field. The same analysis applies for the salinity field.

The above analysis can be justified simply by investigation of the evolution equations of $BPE_\varTheta$ and $BPE_S$. Winters et al. (Reference Winters, Lombard, Riley and D'Asaro1995) derived the evolution equation for the background potential energy in the shear-driven case, which can be directly applied to the evolution for $BPE_\varTheta$ and $BPE_S$ in the double diffusion case of interest to us here since both temperature and salinity are governed by the same advection–diffusion equations as that which governs density in the shear-driven stratified turbulence case discussed in Winters et al. (Reference Winters, Lombard, Riley and D'Asaro1995). In a closed system, these evolution equations can be shown to be

(2.12a)\begin{gather} \frac{\textrm{d}}{\textrm{d} t}BPE_\varTheta= \mathcal{M}_\varTheta+D_{p\varTheta}, \end{gather}

(2.12b)\begin{gather}\frac{\textrm{d}}{\textrm{d} t}BPE_S= \mathcal{M}_S+D_{pS}, \end{gather}

(2.12c)\begin{gather}\mathcal{M}_\varTheta+\mathcal{D}_{p\varTheta}= \kappa_\theta \rho_0 \alpha g\left\langle \frac{\textrm{d} z_{\theta\star}}{\textrm{d} \varTheta}|\boldsymbol{\nabla} \varTheta|^2 \right\rangle, \end{gather}

(2.12d)\begin{gather}\mathcal{M}_S+\mathcal{D}_{pS}= -\kappa_s \rho_0 \beta g \left\langle \frac{\textrm{d} z_{s\star}}{\textrm{d} S}|\boldsymbol{\nabla} S|^2 \right\rangle. \end{gather}

Since both temperature and salinity profiles are sorted to increase upwards, it is clear that $\mathcal {M}_\varTheta +\mathcal {D}_{p\varTheta }>0$ and $\mathcal {M}_S+\mathcal {D}_{pS}<0$ from ((2.12c) and (2.12d)). Thus, we have confirmed our analysis above that $BPE_\varTheta$ will always be a monotonically increasing function with time while $BPE_S$ will be a monotonically decreasing function with time. The distinctions between $\mathcal {M}_\varTheta$($\mathcal {M}_S$) and $D_{p\varTheta }$ ($D_{pS}$) were introduced in Caulfield & Peltier (Reference Caulfield and Peltier2000) and Peltier & Caulfield (Reference Peltier and Caulfield2003), in which $D_{p\varTheta }$ ($D_{pS}$) characterizes energy transfer in the motionless background, and $\mathcal {M}_\varTheta$($\mathcal {M}_S$) represents the part of the irreversible fluxes induced by motion which produces a much larger iso-scalar surface area for mixing to occur.

Distinct available potential energies may then be defined as the difference between the potential energies and background potential energies, namely:

(2.13a)\begin{equation} APE_\varTheta = PE_\varTheta-BPE_\varTheta, \end{equation}

(2.13b)\begin{equation}APE_S = PE_S-BPE_S, \end{equation}

(2.13c)\begin{equation}\frac{\textrm{d}}{\textrm{d} t}APE_\varTheta = \mathcal{H}_\varTheta-\mathcal{M}_\varTheta, \end{equation}

(2.13d)\begin{equation}\frac{\textrm{d}}{\textrm{d} t}APE_S = \mathcal{H}_S-\mathcal{M}_S. \end{equation}

Analysing the energy flows in the available potential energy reservoirs will lead to a qualitatively useful depiction of energy flows between the different energy reservoirs that are active in salt-fingering turbulence, as illustrated in figure 2. Concerning the temperature field, energy is being transferred from the $KE$ reservoir into the $PE_\varTheta$ reservoir at a rate represented by $\mathcal {H}_\varTheta$. A portion of the energy flux is effected through irreversible processes (mixing) $\mathcal {M}_\varTheta$ and this is eventually transferred into the background potential energy $BPE_\varTheta$, while the remaining portion of this energy flux $\mathcal {H}_\varTheta -\mathcal {M}_\varTheta$ enters into the $APE_\varTheta$ reservoir. The energies stored in $APE_\varTheta$ is ‘available’, meaning that it may be released back to the $KE$ reservoir through a negative flux $\mathcal {H}_\varTheta $. On the other hand, for the salinity field, energy is transferred from the $PE_S$ reservoir into the $KE$ reservoir. The irreversible flux $\mathcal {M}_S$ transports energy from background potential energy $BPE_S$ to the $APE_S$ irreversibly, which leads to the monotonic decrease of $BPE_S$. This loss of the energy in $APE_S$ to KE, on the other hand, is reversible and can potentially be returned back to $APE_S$ through a positive $\mathcal {H}_S$.

Figure 2. Exchange of energy fluxes between different energy reservoirs, graphics for (2.10), (2.12) and (2.13).

2.3. Evaluating irreversible salt-fingering fluxes in two different models

We will first discuss the application of the above formalism to the interface model. As the interface model is strongly inhomogeneous in the vertical direction, it is expected to be useful to investigate the variation of the irreversible fluxes of heat and salt in the vertical direction across the finite depth of the interface. The depth-dependent irreversible heat and salt fluxes are defined as

(2.14a)\begin{gather} F^{irr}_\varTheta(z_{\theta \star}) = -\kappa_\theta\frac{\textrm{d} z_{\theta\star}}{\textrm{d} \varTheta}\overline{|\boldsymbol{\nabla} \varTheta'(x,y,z)|^2}_{z_{\theta \star}}, \end{gather}

(2.14b)\begin{gather}F^{irr}_S(z_{s \star}) = -\kappa_s\frac{\textrm{d} z_{s\star}}{\textrm{d} S}\overline{|\boldsymbol{\nabla} S'(x,y,z)|^2}_{z_{s \star}}, \end{gather}

in which the irreversible fluxes have depth-dependence on the sorted coordinate $z_{\theta \star }$ ($z_{s \star }$) and the representation $|\, f(x,y,z) |_{z_{\theta \star } (z_{s \star })}$ refers to the average of the field $f(x,y,z)$ with the same $z_{\theta \star } (z_{s \star })$. Such depth-dependent definitions are consistent with the bulk-averaged definitions (2.12) (as discussed in Salehipour & Peltier (Reference Salehipour and Peltier2015)). Although the above definitions of irreversible heat flux (salt flux) are closely related to the previously mentioned heat (salt)-induced buoyancy flux through $\mathcal {M}_\varTheta (z_{\theta \star })=-g\alpha F^{irr}_\varTheta (z_{\theta \star })$ ($\mathcal {M}_S(z_{s \star })=g\beta F^{irr}_S(z_{s \star })$), we prefer to use these more straightforward quantities of $F^{irr}_\varTheta (z_{\theta \star })$ and $F^{irr}_S(z_{\theta \star })$ to be compared with $F_\varTheta (z)$ and $F_s(z)$ in what follows.

The calculation of irreversible fluxes in the unbounded gradient model, on the other hand, is not nearly so straightforward. As we have discussed in § 2.1, the unbounded gradient model essentially describes an infinitely large, homogeneous system that is defined by applying the periodic boundary conditions in all three coordinate directions. It is for this reason that we have to consider fluid elements that are not ‘directly’ contained in the simulation domain when the time-dependent adiabatic rearrangements of fluid parcels are constructed. Consider the typical iso-temperature contour (non-dimensionalized) in the sketch shown in figure 3(a), the periodic vertical boundary condition applied allows for the occurrence of the intersections of the iso-temperature contour with the top and bottom boundaries. If we sort the fluid elements only in the simulation domain, we would obtain the result shown in figure 3(b), i.e. the sorted profiles will have a larger vertical gradient near the vertical boundaries, which is pathological since the artifacts near the vertical boundaries have been artificially introduced. Thus, we sort the fluid elements in the new domain shown in figure 3(c) which contains essentially the same fluid elements as the original domain and retains the same boundary conditions (repeatable vertically) as before. Sorting in this newly defined domain precisely avoids the inhomogeneous-gradient problem illustrated in figure 3(b) and is applied to correctly sort the fluid parcels involved in the salt-fingering turbulent flows to be discussed in what follows.

Figure 3. Illustration of the incorrect method of sorting fluid elements (a,b) and the correct method of fluid-element sorting (c,d) to calculate irreversible fluxes in the unbounded gradient model. Panels (a) and (c) illustrate a typical iso-temperature contour (or surface in three dimensions) in the unbounded gradient model, (b) shows the sorted temperature when the sort region is depicted by the thick line in (a), while (d) shows the sorted temperature when the sort region is depicted by the thick line in (c).

2.4. An alternative definition of background potential energies in double-diffusive systems

In this subsection we will discuss an alternative possibility for extending the idea of background potential energy from a single-component system to a doubly diffusive system. This approach has been discussed in detail in the recent work of Middleton & Taylor (Reference Middleton and Taylor2020) and we will present and discuss their main results using our notation in this subsection. These results will then be compared with our results presented in the previous subsection in (2.12)–(2.14), on the basis of which we will argue that the definitions of background potential energies proposed in the current paper are distinctly preferable for application to the study of irreversible mixing in doubly diffusive systems.

As we have discussed above, the background potential energy in a single-component system is associated with the stable rearrangement of the fluid parcels based on their densities. Middleton & Taylor (Reference Middleton and Taylor2020) directly extends this definition of background potential energy to the double-diffusive system, as

(2.15)\begin{equation} BPE^{MT}= g \langle \rho z_{\rho \star} \rangle, \end{equation}

where $z_{\rho \star }(x,y,z,t)$ represents the vertical position in the sorted profile of the fluid parcel at $(x,y,z,t)$ and $\rho (x,y,z,t)$ is determined by the temperature and salinity through the linear equation of state (2.1g). In a closed system the evolution of $BPE^{MT}$ can be derived following the same strategy as in Winters et al. (Reference Winters, Lombard, Riley and D'Asaro1995) to be

(2.16)\begin{align} \frac{\textrm{d}}{\textrm{d} t}BPE^{MT} &=-\kappa_\theta g \left\langle \frac{\textrm{d} z_{\rho \star}}{\textrm{d} \rho} |\boldsymbol{\nabla} \rho_\theta |^2 \right\rangle- \kappa_s g \left\langle \frac{\textrm{d} z_{\rho \star}}{\textrm{d} \rho} |\boldsymbol{\nabla} \rho_s |^2 \right\rangle \nonumber\\ &\quad -(\kappa_\theta+\kappa_s)g \left\langle \frac{\textrm{d} z_{\rho \star}}{\textrm{d} \rho} (\boldsymbol{\nabla} \rho_s \boldsymbol{\cdot} \boldsymbol{\nabla} \rho_\theta) \right\rangle, \end{align}

where $\rho _\theta \equiv -\rho _0 \alpha (\varTheta -\varTheta _0)$ and $\rho _s \equiv \rho _0 \beta (S-S_0)$ are the components of density that are determined by the temperature field and salinity field separately. Equation (2.16) is essentially equivalent to equation (2.19) of Middleton & Taylor (Reference Middleton and Taylor2020) in a closed system. In order to make a clearer comparison, we rewrite the evolution equation for the background potential energies (2.12c) using $\rho _\theta$ and $\rho _S$ as

(2.17)\begin{align} \frac{\textrm{d}}{\textrm{d} t}BPE &= \frac{\textrm{d}}{\textrm{d} t}BPE_\varTheta+\frac{\textrm{d}}{\textrm{d} t}BPE_S\nonumber\\ & =-\kappa_\theta g \left\langle \frac{\textrm{d} z_{\theta \star}}{\textrm{d} \rho_\theta} |\boldsymbol{\nabla} \rho_\theta |^2 \right\rangle- \kappa_s g \left\langle \frac{\textrm{d} z_{s \star}}{\textrm{d} \rho_s} |\boldsymbol{\nabla} \rho_s |^2 \right\rangle. \end{align}

The first two terms on the right-hand side of (2.16) have very similar forms to those describing our irreversible buoyancy fluxes defined in (2.17). In fact, there are only two differences between these two forms: firstly, the derivatives of the sorted vertical coordinate are computed with respect to different scalars in (2.16) and (2.17); secondly, the two terms in (2.16) are both positive, while in (2.17) the temperature term is positive whereas the salinity term is negative (in a salt-fingering-favourable configuration). Both of these distinctions are suggesting that the terms in (2.16) cannot be regarded as representations of irreversible mixing processes. Firstly, the appearance of total density $\rho$ in the first two terms of (2.16) renders them not independent: different salinity configurations will change the value of the temperature term and vice versa. However, the salinity dissipation process is clearly independent of the heat dissipation process, making the correlations in (2.16) difficult to interpret. Secondly, salinity mixing in a salt-fingering-favourable environment mixes the upper layer salty water with lower layer fresh water, which would lead to a decrease rather than increase in potential energy. This further fact makes it difficult to interpret the second term in (2.16) as representing salinity mixing.

The third term in (2.16) is negative in a salt-fingering-favourable environment, since thermal and saline density components compensate each other to make $\boldsymbol {\nabla }\rho _\theta \boldsymbol {\cdot } \boldsymbol {\nabla }\rho _s <0$. The magnitude of this negative term can be larger than the first two positive terms in double-diffusive fields (as discussed in Middleton & Taylor (Reference Middleton and Taylor2020)), allowing $BPE^{MT}$ to be released even in a closed system. This property contradicts the original physical meaning of the background potential energy that was defined in a single-component system, which can only increase in a closed system due to irreversible mixing processes. Our methodology described in §§ 2.2 and 2.3, on the other hand, perfectly retains the original physical meaning of background potential energies and provides intuitively reasonable forms for the irreversible fluxes in double-diffusive fields. It therefore seems clear that the formalism we have developed for the separation of irreversible and reversible fluxes in doubly diffusive convection is to be preferred. The price to be paid for its use is that at each timestep in the evolution of the DNS, two full sorts of the individual fluid parcels must be performed. Since the algorithm we employ to perform the required sorts is highly parallelizable (as will be discussed in the next section), however, this fact is not an impediment to the use of this more accurate methodology.

4.1. Direct numerical simulation results for the unbounded gradient model

4.1.1. Stages of time-dependent evolution

The evolution of salt fingers in the unbounded gradient model has been thoroughly studied and is well understood (e.g. Shen Reference Shen1995; Traxler et al. Reference Traxler, Stellmach, Garaud, Radko and Brummell2011). There are generally three stages in this evolutionary process, which can be separated on the basis of the evolution of volume averaged traditionally defined fluxes $F_\varTheta \equiv \langle w'\varTheta ' \rangle$, $F_S\equiv \langle w'S'\rangle$, an example of which is shown in the dashed curves in figure 4(a,b). In the first stage, the fastest growing mode of linear salt-fingering instability experiences exponential growth, which leads to the exponential growth of both vertical heat and salt fluxes as shown in figure 4(a,b). In the second stage, these very large vertical fluxes are eventually suppressed by the development of secondary zig-zag instabilities (as originally discussed by Holyer (Reference Holyer1984)). This zig-zag instability introduces vertical shear of horizontal velocity which acts so as to introduce local tilt along the axis of the original long and slender salt fingers (the first slightly tilted fingers can be seen in figure 4d) and acts so as to break the individual fingers into patches (figure 4e). Following this zig-zag instability induced breakup of the initial fingers, the system then evolves into the third stage of activity in which the turbulent fields of both temperature and salinity have been statistically homogenized as depicted in figure 4(f).

Figure 4. Evolution of non-dimensional heat diffusivities (a) and salt diffusivities (b) as well as flux ratio (c) in a DNS with $R_\rho =2$, both traditional fluxes (or flux ratio) and irreversible fluxes (or flux ratio) are shown in dashed and solid lines, separately. Salinity fields are shown in (d)–(f) at non-dimensional times 0.00507, 0.00512, 0.00766 separately to illustrate the twisting of salt fingers under the action of the secondary zig-zag instabilities.

4.1.2. Results for the irreversible fluxes

As discussed in § 2.2, only the irreversible component of the vertical fluxes may contribute to true mixing in the salt-fingering system and, therefore, to the eddy diffusivity which would be employed in a large-scale model that was unable to resolve the turbulence directly. These fluxes can be calculated in our DNS in the unbounded gradient model following the procedures discussed in § 2.3. The typical evolution of both the irreversible heat flux $F^{irr}_\varTheta$ and the irreversible salt flux $F^{irr}_S$ are shown in the solid lines in figures 4(a) and 4(b) separately. Here $F^{irr}_\varTheta$ follows almost the same trajectory as that for its traditionally defined counterpart $F_\varTheta$: $F^{irr}_\varTheta$ experiences exponential growth in the initial phase of salt finger development, but decreases in strength once the zig-zag instability begins to twist the fingers and eventually equilibrates to very nearly the same value as the equilibrium $F_\varTheta$. However, the evolution of $F^{irr}_S$ differs distinctly from that of $F_S$. In the initial phase of finger growth, $F^{irr}_S$ grows at nearly the same growth rate as $F_S$ but at any particular time remains smaller by more than an order of magnitude. This initial stage is followed by an even faster rate of growth that begins when $F_S$ begins to decline. Here $F^{irr}_S$ ascends to its peak value and then evolves towards the statistical equilibrium strength characteristic of stage 3 in which it has nearly the same value as $F_S$.

The above evolutionary trajectory of the fluxes can be understood by recognizing the physical meaning of irreversible fluxes – they represent the flux contributed by diffusion acting on sufficiently small scales so as to modify the temperature/salinity of the turbulently displaced fluid parcels and their neighbours. The traditionally defined fluxes $F_S$ and $F_\varTheta$, however, only capture the extent to which fluid displaced from the high temperature/salinity environment is inserted into the low temperature/salinity environment. The difference between the strengths of these differently defined fluxes are closely related to rates of molecular diffusion. In the limiting case in which it may be assumed that the tracer diffuses infinitely quickly (horizontally), any displacement of a fluid parcel from a region in which the background tracer field has a low value into a region in which the background tracer field has a higher value will inevitably lead to irreversible mixing at the small scales with the same strength. In this case there will therefore be no distinctions between different definitions of the fluxes. Relaxing this constraint of infinite diffusivity slightly will lead to the case of the temperature field illustrated above: although the absolute value of molecular temperature diffusivity is small (around $1.7\times 10^{-7} \ \textrm {m}^2\ \textrm {s}^{-1}$), it is not a small value in the micro-scale system we are studying considering that we have non-dimensionalized the system in such a way that $\kappa _{\theta *}=1$ (shown in table 1). In this circumstance, any displaced fluid parcel that has a different temperature from that of its surroundings will take some time to mix temperatures locally. It is for this reason that a small phase delay can be observed in figure 4(a): $F^{irr}_\varTheta$ follows almost the same trend as $F_\varTheta$, it simply evolves somewhat more slowly. In other words, the displaced fluid parcels are eventually mixed in temperature with their surroundings quite efficiently.

The molecular diffusivity for salinity, on the other hand, is much smaller than temperature diffusivity ($\kappa _s=0.01\kappa _\theta$ under oceanographic conditions and for the current simulations which are tuned to that environment), which leads to a more complicated picture than that for temperature discussed above. Quantitative analysis can be performed by considering the previously discussed different stages of salt finger development (stage 1) illustrated in figure 4(d), the fluid parcels can travel a very great distance vertically through the surroundings with much higher/lower salinity while leaving its original salinity essentially unchanged due to the extremely small value of the salinity diffusivity. In this process, only vanishingly small irreversible mixing for the salinity field $F^{irr}_S$ occurs even though $F_S$ is large (both $F^{irr}_S$ and $F_S$ are exponentially growing). In stage 2, however, the zig-zag instability begins to twist the fingering structures and to advect fluid parcels with extreme salinities into closer contact with their neighbours. This forces these fluid parcels to exchange salinity with neighbouring parcels (as can be seen in figure 4e), therefore, the more fingering structures are twisted and destroyed, the stronger the irreversible mixing will be. This explains why $F^{irr}_S$ dramatically increases right at the time when the zig-zag instability begins to suppress fingering growth as shown in figure 4(b). This is one of the most fundamental conclusions of this paper: while the traditionally defined salt flux provides a useful diagnostic of the salt-fingering generation mechanism, the irreversible salt flux is capturing the extent to which fingers are being destroyed. The final statistical equilibrium stage of salt-fingering turbulence can best be understood as a state in which there exists a balance between finger generation and disruption (the nature of this balance is captured in the theoretical model of Radko & Smith (Reference Radko and Smith2012)) which suggests that the equilibrium state is obtained when the growth rate of first-order fingering instability is balanced with the growth rate of the zig-zag instability). As the growth rate of fingers being generated at a given time will determine the number of fingers that are available for disruption at a later time, there is an expected phase delay between $F^{irr}_S$ and $F_S$ demonstrated in the high-resolution inset of figure 4(b).

After comparing the evolution of irreversible heat and salt fluxes in the unbounded gradient model above, it is also informative to discuss their ratios: the ratio of irreversible fluxes is defined as

(4.1)\begin{equation} \gamma^{irr}\equiv\frac{|\mathcal{M_\varTheta}|}{|\mathcal{M_S}|}=\frac{|\alpha F^{irr}_\varTheta|}{|\beta F^{irr}_S|}. \end{equation}

The evolution of $\gamma ^{irr}$ is compared with the traditional definition $\gamma$ in figure 4(c). The traditional $\gamma$ in the salt finger growth stage (around 0.584) is a slight overestimation of the value of $\gamma$ in the equilibrium stage (around 0.546), a trend that has been discussed in detail in earlier work of Traxler et al. (Reference Traxler, Stellmach, Garaud, Radko and Brummell2011). However, the ratio of irreversible fluxes $\gamma ^{irr}$ is extremely high in the early stage of evolution of the turbulent flow simply due to the much higher temperature diffusivity than salinity diffusivity. Following the transition stage 2, $\gamma ^{irr}$ equilibrates at a level that is nearly the same at equilibrium as that for $\gamma$. However, the variations about the equilibrium are significantly stronger than those characteristic of the traditional ratio $\gamma$, which suggests a larger error bar is necessary for characterization of the numerically determined irreversible flux ratio. This effect is not due to larger variations in $F^{irr}_\varTheta$ ($F^{irr}_S$) than $F_\varTheta$ ($F_S$) (since they are characterized by fluctuations of similar magnitude as shown in figure 4a,b), but is simply due to the phase delays mentioned above: while irreversible temperature flux remains in phase with the finger generation process at a given time, the irreversible salinity flux is more nearly in phase with the finger disruption mechanism at a given time. Thus, the value of the irreversible flux ratio at a given time actually involves the state of the system that existed at a significantly earlier time, and this inevitably leads to stronger variations with time of $\gamma ^{irr}$ about its equilibrium value.

4.1.3. Equilibrium fluxes and a revised scheme for the parameterization of salt-fingering turbulence in global ocean models

As discussed above, the salt-fingering turbulence in the unbounded gradient model eventually evolves into an equilibrium state which is statistically stable. This makes it possible for us to obtain useful estimates of the diapycnal diffusivities based on the measured fluxes. Specifically, the diapycnal diffusivities can be calculated by substituting the equilibrium fluxes and the constant background gradients into the definitions, for both traditionally defined fluxes and the irreversible fluxes, as

(4.2a)\begin{gather} K_\varTheta = -\frac{F_{\varTheta}}{\varTheta_{z0}}, \end{gather}

(4.2b)\begin{gather}K_S = -\frac{F_S}{S_{z0}}, \end{gather}

(4.2c)\begin{gather}K^{irr}_\varTheta = -\frac{F^{irr}_\varTheta}{\varTheta_{z0}}, \end{gather}

(4.2d)\begin{gather}K^{irr}_S = -\frac{F^{irr}_S}{S_{z0}}. \end{gather}

The equilibrium irreversible diapycnal diffusivity (averaged over a finite-time period in the equilibrium stage) is strongly dependent upon $R_\rho$, which is the only non-dimensional parameter that varies significantly in oceanographic environments. In table 3 we show explicitly the equilibrium irreversible diapycnal diffusivities for heat and salt ($K^{irr}_\varTheta$, $K^{irr}_S$) as well as the irreversible flux ratio $\gamma ^{irr}$ for each of 10 different values of $R_\rho$. These irreversible diapycnal diffusivities are consistent with the traditionally defined diapycnal diffusivities obtained in previous three-dimensional simulations with similar resolutions (Traxler et al. Reference Traxler, Stellmach, Garaud, Radko and Brummell2011; Radko & Smith Reference Radko and Smith2012). The functional dependence of $K^{irr}_\varTheta$ and $\gamma ^{irr}$ on $R_\rho$ can be well fit using the specific forms:

(4.3a)\begin{gather} K^{irr}_\varTheta = \frac{a}{(R_\rho-1)^bR_\rho^c}, \quad (a=78.09, b=0.52, c=0.87), \end{gather}

(4.3b)\begin{gather}\gamma^{irr} = aR_\rho^3+bR_\rho^2+cR_\rho+d,\quad (a=-0.00068, b=0.0163, c=-0.125,d =0.77), \end{gather}

(4.3c)\begin{gather}K^{irr}_S = \frac{K^{irr}_\varTheta}{\gamma^{irr}} R_\rho. \end{gather}

Table 3. Equilibrium $K^{irr}_\varTheta$, $K^{irr}_S$, as well as $\gamma ^{irr}$ averaged over the statistical equilibrium state of simulation numbers 1–10.

Both the data and our fitting results are displayed in figure 5. The expression for $K^{irr}_\varTheta$ we employ is based upon the assumption that $K^{irr}_\varTheta$ goes to infinity as $R_\rho$ approaches 1 and $K^{irr}_\varTheta$ goes to 0 as $R_\rho$ goes to infinity (when salt fingers are extremely weak). We use a third-order polynomial form for $\gamma ^{irr}$ in order to capture the transition from a decreasing $\gamma ^{irr}$ to increasing $\gamma ^{irr}$, as shown in figure 5(c). It is worth noting that these expressions differ from the expressions suggested by Radko & Smith (Reference Radko and Smith2012), in which simulations have been performed over a much more restricted range of $R_\rho$ ($1<R_\rho <3$) to which the similarity theories described in Radko (Reference Radko2008) could be applied.

Figure 5. Averaged $K^{irr}_\varTheta$ (a), $K^{irr}_S$ (b) and $\gamma ^{irr}$ (c) over the statistical equilibrium state as a function of $R_\rho$ from unbounded model simulations (simulation numbers 1–10). Numerical fitting of (4.3) is shown in each figure to be compared with the data.

4.2. Direct numerical simulation results for the interface model

4.2.1. Time evolution of simulations for the interface model

As in the case of the unbounded gradient model, we will begin our discussion of the interface model by first providing an overall description of the evolution of the system.

Figures 6(a) and 6(b) display the evolution of temperature induced buoyancy flux $\mathcal {H}_\varTheta$ and minus salinity induced buoyancy flux $-\mathcal {H}_S$ at each depth on a logarithmic scale from our simulation with initial density ratio $R_{\rho 0}=2$. The fluxes firstly develop around the midpoint of the interface where the background gradients are steepest and then extend to the physical boundaries of the computational domain, a process that is accompanied by the growth and expansion of the fingering-affected region as visualized in figure 6(c,d). The fingers appear to be more abundant and narrower than those depicted in figure 4 due to the much larger domain in which we have elected to simulate salt-fingering turbulence in the interface model.

Figure 6. (a,b) Contour plot of $\mathcal {H}_\varTheta$ and $\mathcal {H}_S$ in the simulation case 11. (c,d) Snapshots of salinity fields taken at two time slices of the simulation, note that the axis scales in (c) and (d) have been adjusted for better visualization purposes.

It is clear that at the time the simulation is arrested the fluxes are still extending vertically in space and have begun to be confined by the upper and lower boundaries of the domain. Continuing the evolution of the system beyond this time would inevitably introduce effects due to boundary reflections which would render the results unphysical. Although the system as a whole is not in statistical equilibrium at this time, a quasi-equilibrium state has already been established across the region of strong background gradients that define the interface (say the region between $-2<z_\ast <2$). In this region, salt-fingering turbulence (illustrated in figure 6b) will be recognized as having the same characteristics as the equilibrium fingering field in the unbounded gradient model (see figure 4f).

4.2.2. Reversible and irreversible buoyancy fluxes in the interface model

Since the interface model is describing a closed system with well-defined potential energies, we are in a position to make detailed tests of the accuracy of satisfaction of the energy budget equations discussed in § 2.2 which we proceed to perform in this section.

We first illustrate the evolution of temperature/salinity potential energies (solid lines) and temperature/salinity background potential energies (dashed lines) in figure 7(a,b). In accord with the discussion in § 2.2, $PE_\varTheta$ is always larger than $BPE_\varTheta$ (the minimum potential energy state for temperature), the difference between which corresponds to the available potential energy $APE_\varTheta$ which is able to be released back to the kinetic energy reservoir. In figure 7(b) the salinity potential energies are decreasing, with $BPE_S$ (the maximum potential energy for salinity) being always larger than $PE_S$. The differences between $PE_S$ and $BPE_S$ correspond to the quantity of potential energy that may be ‘returned’ from the kinetic energy reservoir. The evolution of $PE_S (PE_\varTheta )$ and $BPE_S(BPE_\varTheta )$, as discussed in § 2.2, should follow (2.8) and (2.12) if the system is a closed system. This is tested in figure 7(c,d), where the time derivative of potential energies are plotted as solid lines to compare with their predicted contributions $\mathcal {H}+D_p$ and $\mathcal {M}+D_p$ calculated from their definitions in (2.9) and (2.12). Solid and dashed lines agree well which each other in figure 7(c,d), indicating that the above equations have been obeyed quite well in our numerical simulations. Such agreement also reveals that the current resolution implemented in our DNS is high enough to precisely capture the critical influence of small-scale diffusion.

Figure 7. Time evolution of temperature potential energy (a) and salinity potential energy (b) of simulation case 11. Test of the potential energy budget equations for $PE_\varTheta$, $BPE_\varTheta$(c), as well as $PE_S$ and $BPE_S$ (d).

Although it may seem that the irreversible fluxes are always weaker than the traditionally defined fluxes within the simulation time (shown in figure 7c,d), it is not the case if we compare the irreversible fluxes and traditional fluxes within a finite region near the centreline shown in figure 8 (we choose to focus upon the range $-2<z_\ast <2$ in this case). The evolution of fluxes in this region strongly resembles the evolution process in the unbounded gradient model: there are also evident phase delays in the irreversible fluxes (more evident phase delays in the salinity fluxes than temperature) in the early evolution of the flow, after which both definitions tend to agree with each other once the system enters a statistically steady state. This demonstrates that the system has already reached a quasi-equilibrium state across the width of the transition layer ($-2<z_\ast <2$) (such agreement only occurs in the statistical equilibrium state as shown in our analysis of results from the unbounded gradient model). The irreversible fluxes are generally smaller than the traditional fluxes while salt fingers continue to grow in the outer regions of our domain, averaging over the entire domain will still lead to smaller irreversible fluxes as seen in figure 7(c,d).

Figure 8. Comparison of the evolution of $F_\varTheta$ and $F^{irr}_\varTheta$ (a) as well as $F_S$ and $F^{irr}_S$ (b) averaged over the centreline region $-2<z_\ast <2$.

The above discussions rely upon an important assumption that lies at the core of the present paper, namely that the evolution of salt-fingering turbulence in a local region in the interface model can be understood on the basis of our simulations of salt-fingering turbulence in the unbounded gradient model. Different regions (with different vertical positions) may have different growth rates or different fingering widths, but, we hypothesize, they always follow the same path of turbulence evolution that we have shown to obtain in the unbounded gradient model at the same density ratio. If this hypothesis is correct it both suggests that the local irreversible fluxes will be smaller than traditional fluxes in the finger growth stage (with a phase delay), and that the local equilibrium stage can be characterized by coincidence between the reversible and irreversible fluxes. This hypothesis will be directly tested and established as correct in the next subsection where we will apply the parameterization of diapycnal diffusivity obtained through analysis of data from the unbounded gradient model to the interface model.

4.2.3. Effective diapycnal diffusivity for the interface model

The manner in which we calculate this important quantity for the interface model is necessarily distinctly different than the simple method employed to compute diapycnal diffusivity for the unbounded gradient model. First, the interface model is inhomogeneous vertically and there is always stronger flux near the centre of the interface compared with boundary regions, as clearly shown in figure 6(a). Secondly, no equilibrium state exists in this run-down experiment, and, thus, both the fluxes and the related diapycnal diffusivity are essentially time-dependent. Previous researchers (e.g. Kimura et al. Reference Kimura, Smyth and Kunze2011) computed averages of the fluxes over the interface region as well as a time average over periods during which the buoyancy Reynolds number did not vary significantly to obtain a value for a flux which depends primarily on the non-dimensional parameters of the initial profiles. They further employed the functional dependence so determined to obtain a representation of diapycnal diffusivity dependent upon both $R_\rho$ and $Ri$ in circumstances in which the system is also influenced by shear. In the current work, however, we argue that averaging over a finite region near the centre of the interface will inevitably introduce large errors in the inference of diapycnal diffusivity by implying that there are large variances in both the local density ratios and local fluxes in the averaging region. Explicitly, the density ratio and diapycnal diffusivity as a function of depth defined by the horizontal average of the temperature and salinity fields are as follows:

(4.4a)\begin{gather} R_\rho(z) = \frac{\alpha\bar{\varTheta}_z}{\beta \bar{S}_z}, \end{gather}

(4.4b)\begin{gather}K_\varTheta(z) = -\frac{F_\varTheta(z)}{\bar{\varTheta}_z}, \end{gather}

(4.4c)\begin{gather}K_S(z) = -\frac{F_S(z)}{\bar{S}_z}. \end{gather}

It will be observed that the above formulae are based on the traditionally defined diapycnal diffusivity. As discussed in § 4.2.2, the irreversible fluxes equilibrate to the same value as the traditionally defined fluxes in the statistical equilibrium state. Since our focus will be on the equilibrium properties in this subsection, such distinctions will be irrelevant in the following discussion. In figure 9(a) we show the vertical variations of $R_\rho (z)$ at the centreline region for our simulation number 11. Although $R_\rho (z)$ is initialized to be $R_{\rho 0}=2$ for all depths, the system evolves into a state in which $R_\rho$ is depth dependent. The difference between local $R_\rho (z)$ and the initially specified global value is especially significant near the centre of the interface model, with $R_\rho (z_\ast =0)$ being much larger than 2 and $R_\rho (z_\ast =-2)$, $R_\rho (z_\ast =2)$ much smaller than 2. This pattern does not rely on the time chosen as demonstrated in figure 9(b) in which we have plotted the evolution of $R_\rho (z)$ at each of these three positions as a function of time. As the system approaches the statistical equilibrium state appropriate to a specific depth, the $R_\rho (z)$ also attains a quasi-equilibrium value. In figure 9(b) the equilibrium value of $R_\rho$ is shown to be larger than $R_{\rho 0}$ near the centre of the interface model but smaller than $R_{\rho 0}$ at $z_\ast =2,-2$ whose heights are displaced from the centreline. This pattern arises as the first-order response to the spatial distribution of fluxes shown in figure 6. As strong heat and salt fluxes in the central region of the interface support a strong flux of both heat and salt downwards across the centreline, both temperature gradient $\varTheta _z$ and salinity gradient $S_z$ decrease in the vicinity of $z_\ast =0$. Since this effect is stronger for salinity than for temperature (salt fluxes are stronger in a salt-fingering process), $R_\rho =\alpha \bar {\varTheta }_z / \beta {\bar{S}_z}$ at the centreline will increase because of the more prominent decrease in the denominator than the numerator. For the $z_\ast =2$ region (the same argument applies for $z_\ast =-2$), however, these regions will experience an increase in both $S_z$ and $\varTheta _z$, as both salinity and temperature below $z_\ast =2$ are delivered downwards across the centreline, while salinity and temperature above $z_\ast =2$ remain relatively unchanged (fluxes are much smaller away from the centreline). This effect is also stronger in the salinity field than in the temperature field, leading to a decrease in $R_\rho$ at $z_\ast =2$ (the denominator is increasing more strongly). The same mechanism leads to a density inversion zone (with $R_\rho <1$) in regions further removed from the centre of the transition region (not shown). Such regions have been previously reported and discussed in the two-dimensional numerical simulations of Shen (Reference Shen1989), Shen & Veronis (Reference Shen and Veronis1997) and Singh & Srinivasan (Reference Singh and Srinivasan2014), but further discussion of these regions is beyond the scope of the present paper.

Figure 9. (a) Local $R_\rho$ near the centreline region at $t_\ast =0.098$ for simulation number 11. (b) The time-dependence of local $R_\rho$ at $z_\ast =2, z_\ast =0$ and $z_\ast =-2$.

These large variations with $R_\rho$ as a function of depth indicate that, instead of averaging over a particular depth range (as in Kimura et al. (Reference Kimura, Smyth and Kunze2011)) in which $R_\rho$ varies so intensely (the diapycnal diffusivities and flux ratio evaluated using this method are discussed in appendix B), it would appear to be a better strategy to treat each depth separately, especially insofar as diapycnal diffusivity is concerned. These diapycnal diffusivities do vary significantly as a function of depth, as shown by the solid lines in figure 10. Specifically, we plot in figure 10(a–c) the evolution of $K_\theta$,$K_s$ (both non-dimensionalized by dividing by the molecular diffusivity $\kappa _\theta$) and the flux ratio $\gamma$ at several different depths. It is apparent that the diffusivities at different local depths tend to fluctuate about different values. It is then a well-motivated assumption to make that the diapycnal diffusivities at depth $z$ are controlled by the local density ratio $R_\rho (z)$. To test this hypothesis, we apply the parametrization schemes developed on the basis of results from the unbounded gradient model (4.3) to the $R_\rho (t)$ data at different depths and also plot these data on figure 10(a–c) separately as the dashed curves. These dashed curves should be understood as representing the equilibrium diapycnal diffusivities of the turbulence at a particular value of $R_\rho$ at the corresponding time. It will be observed that there is a generally good match between the dashed curve and the solid curve. Specifically, the dashed curve gives a more accurate prediction at $z_\ast =0$, where $R_\rho$ is in a relatively stationary environment which better resembles the settings of the unbounded gradient model. At $z_\ast =2$ and $z_\ast =-2$, both the solid and dashed curves have large variations due to strong local flux fluctuations, however, the two curves are varying around similar equilibrium levels, which reveals their connections. These results supports the validity of our hypothesis that we can accurately treat the local (depth dependent) environment of the interface model as a particular realization of the unbounded gradient model from which local fluxes for the interface model can be determined from the local value of $R_\rho$.

Figure 10. Time evolution of diapycnal diffusivities for temperature $K_\theta$ (a) and salinity $K_S$ (b), as well as flux ratio $\gamma$(c) at $z_\ast =2, z_\ast =0$ and $z_\ast =-2$. In each figure, the solid line represents data from the interface model and the dashed line represents the prediction of $K_\theta$, $K_s$, as well as $\gamma$ calculated by applying our new parametrization method to the time evolution of local $R_\rho$ shown in figure 9(b).