3.1. Asymmetric density anomalies across a diffusive interface due to nonlinearities in the equation of state

In this subsection, we use a one-dimensional model of molecular diffusion to show how the nonlinear nature of the equation of state affects the density anomalies that form within the layers either side of a diffusive interface. Consider a tracer in one dimension, $C(z, t)$ , with molecular diffusivity $\kappa _{C}$ , evolving according to the diffusion equation

(3.1) \begin{equation} \frac {\partial C}{\partial t} = \kappa _{C}\frac {\partial ^2 C}{\partial z^2}. \end{equation}

From an initial two-layer state with a step change occurring at $z = a$ , concentration $C^{*}$ in the deeper layer and a $\Delta C$ difference across the step change, we have the initial conditions

(3.2) \begin{equation} C(z, 0) = \begin{cases} C^{*} + \Delta C &\text{for } \ z \gt a, \\ C^{*}& \text{for } \ z \lt a. \end{cases} \end{equation}

An analytic solution to the initial value problem posed by (3.1) and (3.2) can be found using an error function:

(3.3) \begin{equation} C(z, t) = C^{*} + \frac {\Delta C}{2}\left (1 + \mathrm{erf}\left (\frac {z - a}{\sqrt {4\kappa _{C}t}}\right )\right )\!. \end{equation}

Using (3.3) we can model the evolution of two tracers with different diffusivities from an initial step change to investigate the impact of nonlinearity on diffusive convection. The two tracers of interest in the ocean are salinity, with molecular diffusivity $\kappa _{S}$ , and temperature, with molecular diffusivity $\kappa _{\varTheta }$ . For salinity and temperature, we set the initial state as in experiments I and IV, but use the nominative values $\kappa _{S} = {1\times {10}^{-9}}\ \textrm{m}^{2}\ \textrm{s}^{-1}$ and $\kappa _{\varTheta } = {1\times {10}^{-7}}\ \textrm{m}^{2}\ \textrm{s}^{-1}$ in the ocean (de Lavergne et al. Reference de Lavergne, Groeskamp, Zika and Johnson2022). We use dimensional quantities here, and when reporting our DNS results below, so as to more readily make comparisons between the linear and nonlinear equation of state, and to observations in § 4. Figure 1(a) shows the initial conditions (3.2) as solid lines, and the profiles at time $t^{*} = {83}\,\textrm{min}$ (dashed lines) computed using (3.3). The exact choice of $t^{*}$ here is not important provided that the turbulent boundary layers have formed either side of the established diffusive interface.

Figure 1.

(a) Initial salinity and temperature (solid blue and red, respectively), and salinity and temperature at $t = t^{*}$ (dashed blue and red, respectively) computed using (3.3), with $S^{*} = 34.7$ , $\Delta S = {-0.12}\ {\textrm {g kg}^{-1}}$ , $\kappa _{S} = {1\times {10}^{-9}}\ {\textrm{m}^2\ \textrm {s}^{-1}}$ and $\varTheta ^{*} = {0.5}\,^\circ {\textrm {C}}$ , $\Delta \varTheta = {-2}\,^\circ {\textrm {C}}$ , $\kappa _{\varTheta } = {1\times {10}^{-7}}\ {\textrm{m}^2\ \textrm {s}^{-1}}$ . The initial density ratio is $R_{\rho } = 1.05$ . (b) Initial density anomaly (solid grey) and density anomalies at $t = t^{*}$ , calculated from the evolved salinity and temperature profiles in (a), using $\rho _{\textit{linear}}$ (dashed blue) and $\rho _{\textit{nonlinear}}$ (dashed orange). In both cases, the density is potential density referenced to 0 dbar, and is presented as an anomaly from the initial mean density.

In figure 1(b), we compare the density evolution of the salinity and temperature profiles in figure 1(a) between linear and nonlinear equations of state. The linear equation of state used is (2.6), and the nonlinear equation of state $\rho _{\textit{nonlinear}} (S,\varTheta ,p )$ is TEOS-10 (IOC et al. 2010) provided by the Gibbs Seawater Oceanographic Toolbox (McDougall & Barker Reference McDougall and Barker2011). To compare the density calculated using the linear and nonlinear equations of state, we plot the potential density (referenced to 0 dbar) anomaly from the initial mean potential density, $\sigma _{0}' = \sigma _{0} - \langle \sigma _{0}^{\textit{initial}} \rangle _{z}$ , in figure 1(b).

The main feature of figure 1(b) is that with both linear and nonlinear equations of state, density anomalies form, and grow, either side of the diffusive interface due to diffusive convection, which can occur as a result of the two orders of magnitude difference in molecular diffusivities. The regions where the density anomalies form are known as unstable diffusive boundary layers, and are what cause turbulent mixing within the layers, while the transport across the diffusive interface is set by the molecular diffusivity of salinity, the slower diffusing tracer (Carpenter et al. Reference Carpenter, Sommer and Wüest2012). The effect of nonlinearity in figure 1(b) is to cause an asymmetry between the density anomalies within each layer. Note that the density difference between the maximum and minimum density formed by double diffusion is the same with both the linear and nonlinear equations of state. The effect of cabbeling is to shift these extrema, leading to the midpoint density not being at exactly the same height as the midpoint salinity and temperature, as is the case when the equation of state is linear.

The reason for the difference in density anomalies between the linear and nonlinear equations of state in figure 1(b) is due to the difference in the geometry of isopycnals in salinity–temperature space. In figure 2, the evolved salinity and temperature profiles (dashed lines in figure 1 a) are mapped into salinity–temperature space. In figures 2(a,b), isopycnals are plotted at the initial shallow layer, minimum, initial deeper layer, and maximum density anomalies. In figure 2(a), $\rho _{\textit{linear}}$ from (2.6) is used to calculate density. The linear isopycnals in figure 2(a) create an equal gap between the initial density of the shallow layer and the minimum density, and the initial density of the deeper layer and the maximum density. This results in equal density anomalies either side of the diffusive interface when the equation of state is linear. In figure 2(b), the curvature of the isopycnals in salinity–temperature space creates a larger gap between the density of the initial deep water and the density maximum than the gap between the density of the initial shallow water and the density minimum. This creates the unequal density anomalies within each layer when the equation of state is nonlinear. This effect was noted by McDougall (Reference McDougall1981b ), and again by Özgökmen & Esenkov (Reference Zika, Sloyan and McDougall1998) and Castro-Folker & Stastna (Reference Castro-Folker and Stastna2023), who suggested that the larger gravitational instability that forms in the deeper layer drives more turbulence in the deeper layer.

Figure 2.

The salinity and temperature profiles from figure 1(a) in salinity–temperature (S–T) space, with the density contours calculated using (a) $\rho _{\textit{linear}}$ and (b) $\rho _{\textit{nonlinear}}$ . The colour bar is categorical, with only the minimum, shallower, deeper and maximum density anomalies shown.

3.2. Diagnosing the asymmetric density anomalies

The results of the one-dimensional model, (3.1) and (3.2) from § 3.1, show that with the nonlinear equation of state, there is an increase in the density anomaly within the deeper layer, and a decrease in the density anomaly within the shallower layer, when compared to a linear equation of state where the density anomalies are equal either side of a diffusive interface (see figures 1 and 2). To diagnose the asymmetry caused by the effects of the nonlinear equation of state on molecular diffusion, we define an asymmetry parameter

(3.4) \begin{equation} \delta _{\textit{mol}} = \frac {|\rho _{\textit{min}} - \rho _{\textit{shallow}}|}{|\rho _{\textit{max}} - \rho _{\textit{deep}}|}, \end{equation}

where $\rho _{\textit{shallow}}$ and $\rho _{\textit{deep}}$ are the densities computed in the shallow and deep layers, respectively, away from the diffusive interface, while $\rho _{\textit{min}}$ and $\rho _{\textit{max}}$ are the density extrema formed by diffusive convection. Finding $\rho _{\textit{shallow}}$ and $\rho _{\textit{deep}}$ can be done as either an average density over the well-mixed part of the layers or, as we do with the one-dimensional model and the initial conditions in our experiments, as the density at the shallowest and deepest waters in the domain.

The asymmetry parameter $\delta _{\textit{mol}}$ is a measure of how the nonlinear equation of state impacts instabilities that arise because of the differing molecular diffusivities of salinity and temperature. The faster diffusing temperature causes the thickness of the temperature interface to grow faster than the salinity interface, i.e. $h_{\varTheta } / h_{S} \gt 1$ , where $h$ is thickness (Carpenter et al. Reference Carpenter, Sommer and Wüest2012), which, irrespective of the equation of state, forms instabilities within the layers either side of the diffusive core. The asymmetry parameter quantifies if the differing thermal expansion and haline contraction coefficients within the layers impact the instabilities that form within each layer when the equation of state is nonlinear. Therefore, the asymmetry parameter does not determine if instabilities occur due to differing diffusivities; rather, it takes into account how these instabilities change when the equation of state is nonlinear.

Two non-dimensional parameters that will influence the asymmetry parameter in (3.4) are the diffusivity ratio $\tau$ , which sets the steepness of the change between $\rho _{\textit{min}}$ and $\rho _{\textit{max}}$ in $S{-}\varTheta$ space, and the density ratio $R_{\rho }$ . A density ratio range $1 \lt R_{\rho } \lt 10$ , especially the lower end of this range, is thought to allow the spontaneous formation of thermohaline staircases in the diffusive convection regime (Kelley et al. Reference Kelley, Fernando, Gargett, Tanny and Özsoy2003), so we restrict ourselves to this range. Oscillatory instabilities can develop when the density ratio is linearly unstable, i.e. in the range $1 \lt R_{\rho } \lt 1.14$ , possibly leading to subcritical diffusive convection (Veronis Reference Urakawa and Hasumi1965). However, the oscillatory modes can be damped by viscous drag, meaning that the formation of diffusive interfaces is the more common outcome in the ocean (Radko Reference Radko2013). Recall that we use (2.7), from McDougall (Reference McDougall1981b ), to compute $R_{\rho }$ , which takes into account the different values of the thermal expansion and haline contraction coefficient within each layer either side of an interface when the equation of state is nonlinear. Nominative values for the molecular diffusivities of salinity and temperature are $\kappa _{S} \approx {1.4\times {10}^{-9}}\,{\,\textrm{m}^{2}\,\textrm{s}^{-1}}$ and $\kappa _{\varTheta } \approx {1.4\times {10}^{-7}}\,{\,\textrm{m}^{2}\,\textrm{s}^{-1}}$ (de Lavergne et al. Reference de Lavergne, Groeskamp, Zika and Johnson2022), leading to $\tau \approx 0.01$ in the ocean. We consider the effect that the diffusivity ratio has on the asymmetry parameter as modelling studies, this one included, employ $\tau \gt 0.01$ to decrease the Schmidt number, which reduces the computational cost of DNS experiments while still capturing double diffusive effects.

Figure 3.

Variation of the asymmetry parameter $\delta _{\textit{mol}}$ in (3.4), at fixed $R_{\rho }$ values with varying $\tau$ . Here, $\delta _{\textit{mol}}$ is calculated from a density profile computed from salinity and temperature profiles evolved using the one-dimensional model (3.3).

Figure 3 shows how the asymmetry parameter $\delta _{\textit{mol}}$ , calculated from salinity and temperature profiles evolved using our one-dimensional model, varies with the diffusivity ratio $\tau$ at three fixed $R_{\rho }$ values. For the linear equation of state, $\delta _{\textit{mol}} \approx 1$ at all times, which is expected because the density anomalies that form within the layers separated by a diffusive interface are equal. For the nonlinear equation of state, $R_{\rho }$ appears to have little influence on the magnitude of the asymmetry. As $\tau$ approaches 1, figure 3 shows a decrease in $\delta _{\textit{mol}}$ . This occurs because in the limit of isotropic turbulent mixing, $\tau \approx 1$ , the cabbeling process, triggered by diapycnal mixing, only forms a density anomaly between the mixed water and the deeper water, i.e. a density anomaly only forms on the deeper side of the interface (Bisits et al. Reference Bisits2025). This results in the numerator of (3.4) being zero when $\tau \approx 1$ . As $\tau$ increases from 0.01 to 1, the density anomalies in both the shallow and deep layers decrease; however, from § 3.1, the density anomaly in the shallow layer is smaller than the density anomaly in the deeper layer, resulting in the decrease of the asymmetry parameter. According to our one-dimensional model, taking $\tau = 0.1$ instead of $\tau = 0.01$ for the DNS experiments in this study leads to an approximately 15 % relative decrease in $\delta _{\textit{mol}}$ .

Figure 4.

(a,d) The density ratio, where values outside the range $1 \lt R_{\rho } \lt 10$ have been removed. (b,e) The asymmetry parameter $\delta _{\textit{mol}}$ calculated from a density profile generated from salinity and temperature profiles evolved using the one-dimensional model (3.3). (c, f) The cabbeling parameter (2.8); when $\delta _{\textit{turb}} \gt 1$ , cabbeling can produce denser water upon mixing. For all calculations, we set $\tau = 0.1$ to match what we use in our simulations. The initial conditions for the experiments listed in table 2 are shown, as well as observed values of $\Delta \theta$ and $\Delta S$ across diffusive interfaces from the Arctic Ocean (Timmermans et al. Reference Timmermans, Garrett and Carmack2008), and the maximum observed $\Delta \varTheta$ across diffusive interfaces in the Baltic Sea (Umlauf et al. Reference Turner2018), Lake Banyoles (Sánchez & Roget Reference Sommer, Carpenter, Schmid, Lueck, Schurter and Wüest2007) and the Weddell Sea (Robertson et al. Reference Radko1995).

We now consider how changes in the temperature and salinity differences across an interface, $\Delta \varTheta$ and $\Delta S$ , respectively, impact the asymmetry parameter $\delta _{\textit{mol}}$ , and try to determine the relative importance of the nonlinear equation of state, using the cabbeling parameter (2.8), at a diffusive interface. In figures 4(a,d), we show the density ratio, in figures 4(b,e) the asymmetry parameter, and in figures 4(c, f) the cabbeling parameter. With a linear equation of state, we see no asymmetry between the density anomalies within each layer (see figure 4(b), where $\delta _{\textit{mol}}=1$ everywhere), and unsurprisingly, no instability to cabbeling ( $\delta _{\textit{turb}} = 1$ indicates no cabbeling instability; see figure 4 c). For the nonlinear equation of state, the main takeaway from figure 4(e) is that $\delta _{\textit{mol}}$ is largely determined by changes to the temperature difference. This is consistent with cabbeling having an influence at diffusive interfaces, as cabbeling is caused by changes to the thermal expansion coefficient with temperature. We also note that the density asymmetry is still present even at larger $R_{\rho }$ values and when the two layers are stable to cabbeling (i.e. $\delta _{\textit{turb}} = 1$ ; see figure 4 f). This suggests that when turbulence is low enough to allow the onset of diffusive convection, the primary effect of the nonlinear equation of state is to create asymmetry between the density anomalies, and that the temperature difference across the interface is the primary driver of the asymmetry, rather than the importance of the nonlinear equation of state as determined by $\delta _{\textit{turb}}$ .

3.3. The effect of asymmetric density anomalies in DNS experiments

3.3.1. Initial evolution of DNS experiments

Figure 5 shows the initial evolution of a vertical slice of the salinity and temperature as an anomaly from the initial conditions for experiment I (figures 5 a–d), linear equation of state, and experiment IV (figures 5 e–h), nonlinear equation of state. In figures 5(a,e), we can see at $t = {3}\,\textrm{min}$ , i.e. the snapshot before turbulence is triggered, the faster diffusing temperature has created a temperature interface that is thicker than the salinity interface. The density anomalies formed by the differing diffusivities then generate turbulence within the layers either side of the diffusive interface with both a linear and nonlinear equation of state (see figures 5(b–d) and 5(f–h), respectively). This shows that the onset of convection when double diffusion is the primary mixing process is largely unaffected by whether the equation of state is linear or nonlinear.

Figure 5.

Evolution of a vertical slice of the salinity anomaly $S'$ on the $x{-}z$ face, the temperature anomaly $\varTheta '$ on the $y{-}z$ face, and the vertically averaged vertical velocity $\langle w \rangle _{z}$ on the $x{-}y$ face at $z = 0$ . The $x{-}z$ vertical slice of the salinity anomaly is taken at $y = {-2.5}\ \textrm{cm}$ , and the $y{-}z$ vertical slice of the temperature anomaly is taken at $x = {-2.5}\ \textrm{cm}$ . The labels in (a) and (e) indicate which field is plotted on which face. (a–d) Experiment I, which uses a linear equation of state. (e–h) Corresponding experiment IV, which uses the nonlinear equation of state.

With the linear equation of state, the tracer anomalies appear to have spread the same vertical distance from the interface in both the shallow and deep layers at $t = {6}\,\textrm{min}$ (figure 5 b). The equal vertical spreading of salinity and temperature continues, and at $t = {24}\,\textrm{min}$ (figure 5 d), the magnitude of the salinity and temperature anomalies either side of the interface look to be equal. Qualitatively, the other experiments in table 2 with a linear equation of state, experiments II and III, are similar, though the turbulence is sustained for less time due to the smaller instability, from the smaller $\Delta \varTheta$ , being mixed away more quickly.

At $t = {6}\,\textrm{min}$ with a nonlinear equation of state (figure 5 f), the temperature anomaly extends further into the deeper layer than the shallower layer, and there is a larger salinity plume in the deeper layer versus the shallower layer. By time $t = {24}\,\textrm{min}$ with a nonlinear equation of state (figure 5 h), there is an accumulation of salinity and temperature anomalies in the shallower layer just above the interface, with no such accumulation occurring below the interface in the deeper layer. This suggests that the tracer anomalies are being mixed away more efficiently in the deeper layer than in the shallower layer. Experiments V and VI are qualitatively similar to IV, though the accumulation of salinity and temperature above the interface is less prominent because of the smaller initial temperature differences.

Figures 6(a–c) compare the horizontally averaged profiles of the density anomaly prior to the onset of convection for experiments with the same initial conditions but a different equation of state. In figures 6(a) and 6(b), we see that the experiments that use the nonlinear equation of state, IV and V, have larger density anomalies in the deeper layer when compared to the equivalent experiments with a linear equation of state, I and II, respectively. The difference between the density anomalies in experiments III and VI (figure 6 c) is much less pronounced, agreeing with our theoretical one-dimensional model in § 3.2 that the influence of nonlinearities in the equation of state on the diffusive convection instability start to become significant at larger temperature differences.

In table 3, we compare the value of the asymmetry parameter in the DNS experiments computed at $t={3}\,\textrm{min}$ , which is prior to the onset of turbulent mixing within each layer, to that predicted by our one-dimensional molecular diffusion model. For experiments I–VI, table 3 shows that there is good agreement between the asymmetry parameter in the DNS experiments and our one-dimensional model – see Appendix D for further details, and how the profiles evolve in salinity–temperature space. This demonstrates that the asymmetry parameter is able to predict how the nonlinear nature of the equation of state can affect the onset of diffusive convection. Further, we suggest that it is this asymmetry in the density anomalies that drives the differing behaviour between the initial evolution of experiments I and IV described above and seen in figure 5. We now investigate how the experiments evolve in time, focusing on comparisons between those that begin from the same initial conditions but use a different equation of state, and whether the value of the asymmetry parameter provides a relative measure of how nonlinearities in the equation of state affect diffusive convection.

Table 3.

Comparison between the asymmetry parameter calculated using the one-dimensional molecular diffusion model, (3.1) and (3.2), from the initial conditions of each experiment, and the value computed from the DNS experiments at time $t={3}\,\textrm{min}$ (prior to the onset of turbulence).

Table 3. Long description

A table comparing asymmetry parameters from different models and experiments. The table has six columns labeled I, II, III, IV, V, and VI, and two rows labeled One-dimensional model and DNS experiments. The values in the table are as follows: Row 1: One-dimensional model, I, 1.0; II, 1.0; III, 1.0; IV, 0.591; V, 0.801; VI, 0.902. Row 2: DNS experiments, I, 0.995; II, 1.011; III, 0.979; IV, 0.632; V, 0.792; VI, 0.917.

3.3.2. Interface evolution available potential energy

Section 3.3.1 shows how the nonlinear equation of state forms an asymmetry between the density anomalies prior to the onset of turbulence within each layer. Now we consider the effect that the asymmetry has on the turbulence within each layer, the volume integrated $\textit{APE}$ , and the diffusive interface as the simulations run down. Figures 6(d–i), show Hovmoller plots of the density anomaly $\sigma _{0}'$ for all experiments in table 2, with the colour bar limits taken from the experiments with a linear equation of state. Note that the density anomaly from the mean in experiment III is shifted by ${1\times {10}^{-3}}\ {\textrm {kg m}^{-3}}$ , so the initial state is more closely aligned to experiment VI. This constant shift is only added to account for the thermal expansion and haline contraction coefficients used in the linear equation of state not being calculated near the midpoint values of salinity and temperature between the layers of experiment III (which results in the initial density anomaly not being centred about zero).

Figure 6.

(a–c) Horizontally averaged profiles of the potential density anomaly from the initial mean density at $t = {3}\,\textrm{min}$ . The dashed lines are from the experiments that include the nonlinear equation of state, and the solid lines are from the corresponding experiments with a linear equation of state. (d–i) Hovmoller plots of the horizontally averaged density as an anomaly from the mean initial density for all experiments: (d–f) experiments I–III, using a linear equation of state; (g–i) experiments IV–VI, using the nonlinear equation of state. Note that the limits on the colour bar for each row have been set to the extrema of the density anomaly from the experiments with a linear equation of state (d–f) to capture the asymmetry in the density anomalies that form with the nonlinear equation of state.

Figure 7.

(a) Time series of the height of the salinity interface, calculated using (3.5) for all experiments listed in table 2. (b) Time series of the density ratio calculated using (2.7) from the mean salinity and temperature over $-0.2\ \text{m} $ to ${-0.1}\ \textrm{m} $ in the shallower layer, and $-0.4\ \text{m} $ to ${-0.3}\ \textrm{m} $ in the deeper layer. (c) Plot of $\textit{APE}$ relative to initial $\textit{APE}$ , i.e. $\mathcal{APE}$ from (3.6), for all experiments. (d) Time derivative of $\mathcal{APE}$ for experiments I and IV. Experiments with a linear equation of state are solid lines, and experiments with a nonlinear equation of state are dashed lines.

A noticeable effect when comparing figures 6(d) and 6(g), experiments I and IV, is an upward migration of the diffusive interface when the equation of state is nonlinear. To track the interface height in our experiments, we reshape and sort the three-dimensional salinity field into a monotonic background profile $\hat {S}(z^{*}, t)$ , where $z^{*} \in [0, L_{z}]$ and the salinity is sorted from saltiest to freshest, equally partitioned with the total number of elements within the three-dimensional domain (Winters et al. Reference Wang, Chong, Stevens, Verzicco and Lohse1995; Carpenter et al. Reference Carpenter, Sommer and Wüest2012). Inverting $\hat {S}(z^{*}, t)$ and finding the height of the midpoint salinity, $\hat {S}_{\textit{mid}}$ , gives us a method to calculate the salinity interface height:

(3.5) \begin{equation} z^{*} \big(\hat {S}_{\textit{mid}}, t\big) = \int _{\hat {S} \lt \hat {S}_{\textit{mid}}}\mathrm{d}z^{*}. \end{equation}

Note an increase in $z^{*}(\hat {S}_{\textit{mid}}, t)$ indicates upward migration of the salinity interface.

Figure 7(a) shows (3.5) plotted against time. We can see that with the nonlinear equation of state, experiments IV–VI, there is an upward migration of the diffusive interface before it plateaus at latter stages. In contrast, with a linear equation of state, the turbulence causes small fluctuations to the interface height at the start of the simulations, but once subsided, the interface appears to be approximately at its initial height $z^{*} = {0.25}\ \textrm{m}$ . This interface migration was observed by McDougall (Reference McDougall1981b ) in laboratory experiments, and attributed to more active turbulence in the deeper layer. McDougall (Reference McDougall1981b ) suggested that the interface migration occurred only when $R_{\rho } \lt 1.6$ . In figure 7(a), we can see that the interface is still migrating in experiment IV, though at a slower rate, when $R_{\rho } \gt 1.6$ – see figure 7(b) for the time series of the density ratio that is calculated using (2.7) from the mean salinity and temperature over $-0.2\ \text{m} $ to ${-0.1}\ \textrm{m} $ in the shallower layer, and $-0.4\ \text{m} $ to ${-0.3}\ \textrm{m} $ in the deeper layer. Further, figures 4(d,e), our one-dimensional model shows that the asymmetry parameter can still be significant when $R_{\rho } \gt 1.6$ provided that there is a relatively large temperature difference between the two layers. These results suggest that interface migration can occur at higher values of $R_{\rho }$ than suggested by McDougall (Reference McDougall1981b ) provided that there is a temperature difference sufficient to create asymmetric density anomalies either side of the interface.

To drive the turbulence, the kinetic energy reservoir needs a source of $\textit{APE}$ . In experiments I–VI, the only source for the kinetic energy is $\textit{APE}$ . Middleton & Taylor (Reference Middleton and Taylor2020) have shown that double diffusion creates a pathway from the $\textit{BPE}$ reservoir to the $\textit{APE}$ reservoir, which is not the case when the diffusivities are equal. The addition of the nonlinearity of the equation of state, which can create $\textit{APE}$ due to cabbeling’s densification upon mixing, also provides a possible source of $\textit{APE}$ (Bisits et al. Reference Bisits2025). We also note that there is very little difference between the density ratio throughout the simulations that have the same initial conditions but use a different equation of state – see figure 7(b). This suggests that the overall stability of the layers is not likely a contributor to any differences in $\textit{APE}$ between experiments with a linear equation of state and their counterparts with a nonlinear equation of state. Figure 7(c) shows the volume integrated $\textit{APE}$ referenced to the initial $\textit{APE}$ ,

(3.6) \begin{equation} \mathcal{APE} = \frac {\textit{AP}E - \textit{APE}_{t=0}}{\textit{APE}_{t=0}}, \end{equation}

for all experiments, and figure 7(d) shows the time derivative of $\mathcal{APE}$ for experiments I and IV. In all experiments, there is an initial generation of $\textit{APE}$ that is then gradually eroded, with intermittent production, for the remainder of the experiment.

The experiments with a nonlinear equation of state do produce marginally more $\textit{APE}$ than the corresponding experiments with a linear equation of state. The time derivative of $\textit{APE}$ in experiments I and IV shows similar behaviour, with the rate of $\textit{APE}$ production at times $t \approx {15}{-}{25}\,\textrm{min}$ and $t \approx {35}{-}{45}\,\textrm{min}$ slightly higher in experiment IV, but again it is not a substantial difference between the linear and nonlinear equations of state. For cabbeling to produce $\textit{APE}$ , the value of the cabbeling parameter must fall in the narrow red range in figure 4(e), which corresponds to density ratios that are quite close to 1. In the experiments with larger initial Rayleigh numbers, such as I and IV, the density ratio increases quickly, moving out of the range where cabbeling is able to generate $\textit{APE}$ upon mixing. Comparing the evolution of $\textit{APE}$ for experiments III and VI, we do see that there is noticeably more $\textit{APE}$ produced when the equation of state is nonlinear. An explanation for this is that the slower increase in the density ratio allows a longer window in the unstable to cabbeling regime where cabbeling can produce $\textit{APE}$ upon mixing. These results indicate that diffusive convection is the primary mechanism for generating the $\textit{APE}$ required to drive turbulence in our experiments, with nonlinearities contributing little to $\textit{APE}$ generation except in situations where the system can remain unstable to cabbeling for longer.

We suggest that rather than generating significantly more $\textit{APE}$ due to cabbeling, the effect of the nonlinear equation of state is to redistribute the $\textit{APE}$ so that there is more $\textit{APE}$ in the deeper layer versus the shallower layer because of the asymmetric density anomalies. This is consistent with what we see when comparing figures 6(d) and 6(g), and to lesser extent when comparing figures 6(e) and 6(h), and figures 6(f) and 6(i), where larger density anomalies reach deeper into the lower layer when the equation of state is nonlinear, indicating that there is more turbulence in the deeper layer, which agrees with McDougall (Reference McDougall1981b ) and Castro-Folker & Stastna (Reference Castro-Folker and Stastna2023). Our results show that smaller values of the initial asymmetry parameter, i.e. increased asymmetry in the density anomalies between the shallow and deep layers, provide a good indicator of when effects of the nonlinear nature of the equation of state are likely to influence diffusive interfaces. Further, it is only when the initial value of the asymmetry parameter is less than 1 that we see interface migration and larger density anomalies in the deeper layer than the shallow layer, indicating that the cause of these effects is the asymmetric density anomalies that form due to the nonlinear equation of state.

3.4. Sustaining diffusive interfaces

A characteristic of thermohaline staircases, and the diffusive interfaces that make up these staircases, is their spatial and temporal coherence (Radko Reference Radko2013). Our simulations do not capture this because the no-flux vertical boundary conditions mean that the system runs down. Ideally, we would like to explore how nonlinearity impacts the spontaneous formation and maintenance of thermohaline staircases by evolving salinity and temperature anomalies about increasing linear vertical gradients in a triply periodic domain (see e.g. Flanagan et al. (Reference Flanagan, Radko, Shaw and Stanton2014) and Ma & Peltier (Reference Ma and Peltier2022), who do this with a linear equation of state). However, with the nonlinear equation of state, setting background linear gradients in salinity and temperature does not result in a linear background density gradient. Further, there is a discontinuity in the density gradient between the top and bottom of the domain. In our attempts, this discontinuity led to an accumulation of denser water at the bottom of the domain and lighter water at the top of the domain, which prevented any staircase formation. Running the same set-up with a linear equation of state did lead to spontaneous staircase formation as in Flanagan et al. (Reference Flanagan, Radko, Shaw and Stanton2014) and Ma & Peltier (Reference Ma and Peltier2022), indicating that further work needs to be done to extend the unbounded gradient model (Radko Reference Radko2013) when the equation of state is nonlinear. This led us to assess the impact of nonlinearity on a sustained interface using a different approach.

From our ‘run-down’ simulations in table 2, we calculate a time series of the molecular fluxes of salinity and temperature across the interface from the change in salinity and temperature content below the midpoint salinity and temperature iso-surfaces as in Carpenter et al. (Reference Carpenter, Sommer and Wüest2012). Then, with the goal of reaching a quasi-steady state for $\Delta S$ , $\Delta \varTheta$ and $R_{\rho }$ , we set constant salinity and temperature flux boundary conditions at the top and the bottom of the domain for a chosen $R_{\rho }$ . This means that we add any heat/salt lost from the deeper layer to the shallower layer into the bottom of the domain using constant heat/salt flux. Similarly for the shallower layer, the flux at the top of the domain compensates for the heat/salt gained from transport across the diffusive interface.

We carried out the procedure described in the preceding paragraph for experiments II and V to produce two more simulations, $\mathrm{{II}}'$ and $\mathrm{{V}}'$ , that are re-runs of II and V with constant flux boundary conditions for salinity and temperature added to the top and bottom of the domain. Note that in experiments $\mathrm{{II}}'$ and $\mathrm{{V}}'$ , the vertical boundary conditions for momentum remain insulating, and horizontally the boundary conditions for momentum, salinity and temperature remain periodic. We aimed to maintain a constant value $R_{\rho } \approx 1.3$ . Calculating the salinity and temperature fluxes across the diffusive interface for experiment II when $R_{\rho } \approx 1.3$ gave us $F_{S}^{\textit{{II}}} \approx {4.77\times {10}^{-8}}\ {\textrm {g kg}^{-1}\ \textrm{m s}^{-1}}$ and $F_{\varTheta }^{\textit{{II}}} \approx {4.5\times {10}^{-6}}\,{^\circ \textrm {C m s}^{-1}}$ . These values are used as constant flux boundary conditions for salinity and temperature at the top and bottom of the domain in $\mathrm{{II}}'$ . From experiment V, we found $F_{S}^{\textit{{V}}} \approx {1.272\times {10}^{-7}}\ {\textrm {g kg}^{-1}\ \textrm{m s}^{-1}}$ and $F_{\varTheta }^{\textit{{V}}} \approx {6.4\times {10}^{-6}}\,{^\circ \textrm {C m s}^{-1}}$ when $R_{\rho } \approx 1.3$ . To achieve the quasi-steady state with the nonlinear equation of state, we also needed to reduce $F_{S}^{\textit{{V}}}$ by half at the top of the domain because of the smaller density anomaly driving less turbulence in the shallower layer. Figure 8 shows the salinity and temperature differences between the shallow and deep layers, and $R_{\rho }$ for $\mathrm{{II}}'$ and $\mathrm{{V}}'$ (figures 8(a,b), respectively). In both experiments, the salinity and temperature differences and $R_{\rho }$ achieve a quasi-steady state after $t \approx {250}\,\textrm{min}$ , with $R_{\rho } \approx 1.27$ and $R_{\rho } \approx 1.33$ in $\mathrm{{II}}'$ and $\mathrm{{V}}'$ , respectively. For experiment $\mathrm{{V}}'$ , we also find $\delta _{\textit{turb}} \approx 1.3$ when the system is in a quasi-steady state showing that effects of the nonlinear equation of state will be important. This quasi-steady state allows us to observe how nonlinearities impact the diffusive interface for a longer period of time.

Figure 8.

Runs of experiments $\mathrm{{II}}'$ and $\mathrm{{V}}'$ (solid lines). The temperature and salinity in the shallower layer are the averages over $z \in [{-0.2}, {-0.1}]\ \textrm{m}$ , and in the deeper layer over $z \in [{-0.4}, {-0.3}]\ \textrm{m}$ ; then $\Delta \varTheta$ and $\Delta S$ are the differences between these two averages, normalised by their initial values. For comparison, the salinity interface height from the ‘run-down’ experiments II and V is also shown in (c) (dashed lines).

In figure 8(c), we see that the upward interface migration continues in experiment $\mathrm{{V}}'$ , whereas in experiment $\mathrm{{II}}'$ , the interface remains approximately at its starting height. From the slope of the linear fit to the salinity interface height for $t = {250}{-}{420}\,\textrm{min}$ , when we have a quasi-steady state, we estimate the entrainment velocity $U_{e}$ , which is a measure of the rate of mechanical entrainment of fluid into the deeper layer from the shallower layer used in McDougall (Reference McDougall1981b ). The product of the entrainment velocity and the temperature difference across an interface, $U_{e}\,\Delta \varTheta$ , is then the entrainment flux of temperature across the interface. A non-dimensional entrainment parameter $E$ can then be defined as the ratio of the entrainment flux of temperature and the temperature flux across the interface:

(3.7) \begin{equation} E = \frac {U_{e}\,\Delta \varTheta }{F_{\varTheta }^{\textit{interface}}}. \end{equation}

See Appendix C for how (3.7) relates to the entrainment parameter defined in McDougall (Reference McDougall1981b ). When $E = 0$ , which is only possible for $U_{e} = 0$ when there is an interface separated by $\Delta \varTheta$ , there is no net entrainment of fluid from the shallower layer into the deeper layer, resulting in symmetric temperature flux across the interface, and no upward migration of the interface. Computing (3.7) for experiment $\mathrm{{V}}'$ using $\Delta \varTheta = 0.7$ , the average value during $t = {250}{-}{420}\,\textrm{min}$ , and $F_{\varTheta } = F_{\varTheta }^{\textit{{V}}}$ , we find $E \approx 0.03$ . This means that in experiment $\mathrm{{V}}'$ , the entrainment flux of temperature is 3 % of the temperature flux across the diffusive interface in the quasi-steady period $t = {250}{-}{420}\,\textrm{min}$ . For $\mathrm{{II}}'$ , we find $E \approx {7.8\times {10}^{-4}}$ during $t = {250}{-}{420}\,\textrm{min}$ , showing that the entrainment flux is negligible when compared to the temperature flux. The non-zero value for $E$ in experiment $\mathrm{{II}}'$ can be explained by the interface height experiencing perturbations around its initial height, as we see in figure 8(c), due to the turbulence that double diffusion drives. Experiments II $'$ and V $'$ demonstrate that the effects of the nonlinear equation of state seen in experiments I–VI, upward interface migration and more turbulent mixing in the deeper layer driven by asymmetric density anomalies, can persist when a sufficiently large temperature difference is maintained.

Our results show that the density anomalies that form from the differing molecular diffusivities of salinity and temperature are influenced by the nonlinear equation of state. The presence of differing density anomalies can be detected using the asymmetry parameter $\delta _{\textit{mol}}$ . The asymmetry parameter can be estimated using our one-dimensional model provided that the differences in salinity and temperature across an interface, and the salinity and temperature in the deeper layer, are known. As the asymmetry parameter decreases from 1, where the density anomalies either side of the diffusive interface are equal, to zero, the asymmetry in density anomalies increases. This provides a relative measure of how significant the effect of nonlinearities in the equation of state are on diffusive convection. The temperature difference, rather than the density ratio as previously thought, appears to be the primary mechanism for creating the asymmetry, suggesting that cabbeling is playing a role at diffusive interfaces. Figure 4(e) shows that for $ {-0.25} \lt \Delta \varTheta \lt {0}\,^\circ {\textrm {C}}$ , $\delta _{\textit{mol}} \gt 0.95$ , suggesting that the asymmetry in density anomalies will be subtle until the temperature difference across a diffusive interface is larger than $|\Delta \varTheta | \gt {0.25}\,^\circ {\textrm {C}}$ .

The most noticeable effect of the nonlinear equation of state in our DNS experiments is an upward migration of the diffusive interface. The most significant interface migration in our run-down experiments occurs in experiment IV, which has the smallest value of $\delta _{\textit{mol}}$ . The extent of the interface migration reduces as $\delta _{\textit{mol}}$ increases in experiments V and VI, agreeing with our one-dimensional theoretical model that the magnitude of $\delta _{\textit{mol}}$ provides a relative measure of how nonlinearities in the equation of state influence diffusive convection. As experiment V $'$ shows, the interface migration can persist, which may impact the stability of thermohaline staircases. When the equation of state is nonlinear, we propose that the interface migration is caused by the asymmetry in density anomalies shifting the distribution of $\textit{APE}$ to drive more turbulence in the deeper layer. That is, whilst there is not notably more total $\textit{APE}$ produced when the equation of state is nonlinear, the deeper layer has more $\textit{APE}$ , leading to more turbulence in the deeper layer than the shallower layer. The implications of our results on diffusive interfaces in the ocean are discussed in § 4.