2 Formulation

The temperature and salinity fields $(T_{tot}^{\ast },S_{tot}^{\ast })$ are separated into linear vertical background profiles $(T_{bg}^{\ast },S_{bg}^{\ast })$ and a departure $(T^{\ast },S^{\ast })$ from them:

where $(A_{T},A_{T0},A_{S},A_{S0})$ are constants. The asterisks hereafter denote dimensional quantities and the subscripts ‘ $tot$ ’ represent the total field variables. Our focus is on the diffusive case, in which both temperature and salinity decrease upward: $\unicode[STIX]{x2202}T_{bg}^{\ast }/\unicode[STIX]{x2202}z^{\ast }=A_{T}<0$ and $\unicode[STIX]{x2202}S_{bg}^{\ast }/\unicode[STIX]{x2202}z^{\ast }=A_{S}<0$ . In the present version of the theory, we ignore planetary rotation, compressibility and the nonlinearity of the equation of state, which reduces the governing Boussinesq equations of motion to

where $\boldsymbol{v}_{tot}^{\ast }=(u_{tot}^{\ast },v_{tot}^{\ast },w_{tot}^{\ast })$ is the velocity, $p_{tot}^{\ast }$ is the dynamic pressure, $g$ is gravity and $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ are the thermal expansion and haline contraction coefficients.

Consider the basic state (denoted by overbars) representing unidirectional shear flow with a sinusoidal velocity profile:

It is assumed that the shear (2.3) is maintained against viscous dissipation by the externally imposed pressure gradient force

For convenience, the governing equations are non-dimensionalized using $H$ as the unit of length, $U_{0}$ as the unit of velocity and $H/U_{0}$ as the unit of time. The non-dimensionalization is implemented as follows:

which transforms the governing equations to

where

is the background density ratio,

is the Richardson number, $N_{b}$ is the buoyancy frequency, $Pe=(U_{0}H)/k_{T}$ is the Peclet number, $\unicode[STIX]{x1D70F}=k_{S}/k_{T}$ is the diffusivity ratio and $Pr=\unicode[STIX]{x1D708}/k_{T}$ is the Prandtl number. It should be noted that the Richardson number is based on the maximal shear $(2\unicode[STIX]{x03C0}U_{0}/H)$ . While this convention leads to a rather cumbersome coefficient of the buoyancy term in the vertical momentum equation, it represents the most physical definition of the Richardson number. The Peclet number can also be expressed in terms of $Ri$ and $N_{b}$ as follows:

Representative values of the Richardson number and the buoyancy frequency in the diffusive regions of the World Ocean are $Ri\sim 1-10$ , $N_{b}\sim 10^{-3}s^{-1}$ . The relevant vertical scale of the basic shear that could trigger thermohaline–shear instabilities is highly uncertain. If we assume that $H\sim 1m$ , then (2.9) implies $Pe\sim 10^{3}$ , which provides very rough guidance for the exploration of the parameter space. The basic state (2.3), (2.4) reduces in non-dimensional units to

The non-dimensional perturbation velocity and pressure are defined accordingly as

3 Direct numerical simulations

Preliminary insights into thermohaline–shear destabilization were first generated from numerical solutions. The governing equations (2.6) were integrated in time using the dealiased pseudospectral method (e.g. Stern, Radko & Simeonov Reference Stern, Radko and Simeonov2001; Stellmach et al. Reference Stellmach, Traxler, Garaud, Brummell and Radko2011) with periodic boundary conditions applied to the perturbation fields in each spatial direction. Both two- and three-dimensional DNS were performed. Double-diffusive DNS are computationally expensive, particularly for thermohaline layering experiments. In this case, it becomes necessary to resolve a wide range of dynamically significant scales – from the anticipated step height, which greatly exceeds the scale of heat dissipation, to the salinity dissipation scale, which is smaller than the heat scale by a factor of $\sqrt{\unicode[STIX]{x1D70F}}$ . In this study we consider oceanographically relevant parameters $(Pr,\unicode[STIX]{x1D70F})=(10,0.01)$ , and such a small diffusivity ratio restricts the integration periods accessible by three-dimensional DNS. Two-dimensional DNS are certainly less demanding computationally and can be used to simulate relatively weak instabilities evolving on longer time scales. Furthermore, as will be shown here, they adequately capture key features of the nonlinear evolution of thermohaline–shear instabilities – staircase formation and subsequent layer-merging events.

Figure 1 presents a typical 3D simulation. The chosen parameters $(Ri,R_{\unicode[STIX]{x1D70C}},Pe)=\left(0.5,2,10^{4}\right)$ define a system that is individually stable with respect to dynamic and to diffusive instabilities. This calculation was initiated by the basic state (2.10) perturbed by small-amplitude random initial $(T,S)$ distributions. The computational domain $(L_{x},L_{y},L_{z})=(2,1,1)$ , which corresponds to a dimensional size of $6~\text{m}\times 3~\text{m}\times 3~\text{m}$ , was resolved by a numerical mesh with $768\times 384\times 384$ grid points. The evolution of this system is illustrated by a sequence of temperature perturbation snapshots (figure 1 a–f). It should be noted that the perturbation initially (figure 1 a) attains its largest values in the low-shear part of the flow $(z\approx 3/4)$ – one of the many features distinguishing this phenomenon from Kelvin–Helmholtz instability. This pattern gradually amplifies in time, leading to a substantial nonlinear modification of the vertical stratification and shear. As a result, the Richardson number is reduced below the critical value of $Ri_{cr}=1/4$ in selected locations and becomes susceptible to Kelvin–Helmholtz instabilities. These secondary instabilities, clearly visible in figure 1(b), tend to be localized to high-shear $(z\approx 0,1/2,1)$ regions. The ensuing mixing produces a series of well-defined horizontally uniform layers (figure 1 c). While the layers formed initially are relatively thin, they merge sequentially, until there is only one well-defined layer left by $t=20\,902$ (figure 1 f).

Figure 1. Three-dimensional DNS of layering in the dynamically and diffusively stable regime $(R_{\unicode[STIX]{x1D70C}},Ri,Pe,\unicode[STIX]{x1D70F})=(2,0.5,10^{4},0.01)$ . The instantaneous temperature perturbations are shown in (a–f) at $t=4012,5007,6010,7054,15\,014$ and 20 902 respectively.

The layering dynamics is revealed more clearly by the plots (figure 2) of the total horizontally averaged density profiles, which include the contribution from linear background stratification:

where $\langle \cdots \rangle$ denotes averaging in $(x,y)$ . These profiles demonstrate the formation of relatively homogeneous layers separated by high-gradient interfaces – structures that are suggestive of the observed thermohaline staircases. The transition from uniformly stratified to layered systems is slow, with a growth rate of $\unicode[STIX]{x1D706}\sim 1.5\times 10^{-3}$ , which corresponds to a dimensional time scale of several days.

Figure 2. The horizontally averaged density profiles for the simulation in figure 1. One can note the formation of relatively small homogeneous layers separated by high-gradient interfaces (c) and their subsequent merger (d–f).

During the early stages of the experiment (figure 1 a,b), the perturbation is relatively uniform in the $y$ -direction, which suggests that the most rapidly growing instability mode is quasi-two-dimensional and therefore the linear phase can be adequately represented by the 2D model. This proposition was readily confirmed by reproducing the foregoing simulation in two dimensions $(x,z)$ and observing similar evolutionary patterns. For instance, figure 3 presents an analogous two-dimensional experiment for $Ri=1$ (all other parameters were kept the same). The increase in the Richardson number results in a much lower growth rate $(\unicode[STIX]{x1D706}\sim 2\times 10^{-4})$ , which places this case beyond the current capabilities of 3D DNS. Nevertheless, this system evolves in a very similar manner – the compact perturbation localized in the low-shear zone gradually amplifies (figure 3 a–c) and eventually transforms the uniform stratification into a well-defined staircase. As previously (cf. figure 2), the layers that are formed first are not steady but undergo a series of merging events, until there is only one step left within the limits of the computational domain. The tendency for layers to merge is a generic property of thermohaline staircases, which is realized in both diffusive and fingering configurations, in two and in three dimensions (e.g. Radko et al. Reference Radko, Bulters, Flanagan and Campin2014a ,Reference Radko, Flanagan, Stellmach and Timmermans b ).

Figure 3. Two-dimensional DNS of layering for $(R_{\unicode[STIX]{x1D70C}},Ri,Pe,\unicode[STIX]{x1D70F})=(2,1,10^{4},0.01)$ . The instantaneous temperature perturbations are shown in the left panels and the horizontally averaged density profiles are on the right for $t=30\,575,49\,020,52\,118$ and 55 823 in (a–d) respectively.

The relatively slow evolution of the perturbations in figures 1–3 is consistent with our interpretation that these instabilities are being driven by the molecular diffusion of density components. Since the foregoing simulations are only weakly dissipative $(Pe=10^{4})$ , it is not surprising that the corresponding non-dimensional growth rates are also very low $(\unicode[STIX]{x1D706}\sim 10^{-3}{-}10^{-4})$ . The reduction in the Peclet number, however, dramatically increases both the growth rates and the parameter range of thermohaline–shear instability. Figure 4, for instance, presents the simulation performed for $(Pe,Ri)=(100,10)$ . The perturbation in this case amplifies much faster $(\unicode[STIX]{x1D706}\sim 0.02)$ , and the Richardson number is a factor of 40 greater than the critical value of $Ri_{cr}=1/4$ suggested by classical theory for non-dissipative shear flows (Richardson Reference Richardson1920; Howard Reference Howard1961; Miles Reference Miles1961).

Figure 4. Two-dimensional DNS in the low- $Pe$ regime: $(R_{\unicode[STIX]{x1D70C}},Ri,Pe,\unicode[STIX]{x1D70F})=(2,10,10^{2},0.01)$ . The instantaneous temperature perturbations are shown in the left panels and the horizontally averaged density profiles are on the right for $t=108,311,500$ and 600 in (a–d) respectively.

Another notable difference between the low- $Pe$ (figure 4) and high- $Pe$ (figure 3) simulations is the much larger ratio of the vertical/horizontal wavenumbers $(m/k)$ in the unstable modes. This ratio was $m/k\sim 2$ in high- $Pe$ experiments, but the most unstable solutions were found for $m/k\sim 50$ in low- $Pe$ cases, a finding that is supported and rationalized by the following (§ 4) linear stability analysis. However, there are a number of similarities in the evolution of the low- $Pe$ system in figure 4 and high- $Pe$ systems (e.g. figure 3) as well. The transition starts with the exponential growth of unstable perturbations (figure 4 a,b). When instabilities amplify to a level sufficient to substantially modify the vertical density and velocity gradients, they produce isolated regions that are susceptible to Kelvin–Helmholtz instability (figure 4 c). These secondary instabilities manifest themselves in the form of relatively narrow eddies, with $O(1)$ lateral extent, which are clearly visible in the late stages of the experiment (figure 4 c,d). The ensuing intense mixing generates patches of relatively uniform density separated by interfaces (figure 4 d). Since the system in figure 4 is much more dissipative than our previous examples, the staircase is not as well defined as in figures 1–3. The interfaces are diffused considerably by molecular dissipation and the mixed layers are not as homogeneous as in high- $Pe$ simulations. Nevertheless, we expect that in nature, the initiation of layering by means of thermohaline–shear instability represents only the first evolutionary stage. Layer-merging events (e.g. Radko et al. Reference Radko, Bulters, Flanagan and Campin2014a ,Reference Radko, Flanagan, Stellmach and Timmermans b ) can systematically increase the step heights in the staircase, concurrently increasing the effective Peclet number, and eventually produce a series of thick well-defined layers.

Figure 5. Two-dimensional DNS in the low- $Pe$ regime: $(R_{\unicode[STIX]{x1D70C}},Ri,Pe,\unicode[STIX]{x1D70F})=(2,10,10^{2},0.01)$ . The vertical extent of the computational domain is increased relative to that in figure 4 to $L_{z}=64$ . The instantaneous temperature perturbations are shown for $t=1009,2005,3500$ and 4040 in (a–d) respectively.

Figure 6. (a) Space–time diagram of the mean density gradient for the experiment in figure 5. (b) Space–time diagram for the experiment in which $Ri=100$ ; other parameters are the same as in (a).

Numerical support for the layer-merging hypothesis is provided by the simulation in figure 5, which is identical in all respects to the calculation in figure 4, except that the vertical size of the computational domain is now $L_{z}=64$ . As expected, the dominant vertical scale of instability modes that emerge first matches the wavelength of the basic shear (figure 5 a). These modes produce a system of horizontal layers. In time, however, layers merge sequentially and their number reduces from $N=64$ at $t=1009$ (figure 5 a) to $N=3$ at $t=4040$ (figure 5 d). The merging pattern is revealed very clearly in the space–time diagram of the horizontally averaged density gradient $-\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{mean}/\unicode[STIX]{x2202}z$ (figure 6 a). Radko (Reference Radko2007) suggested that layered systems can exhibit two merging patterns: $B$ and $H$ modes. $B$ -mergers represent the strengthening of stronger interfaces at the expense of weaker interfaces, which gradually erode. $H$ -mergers occur when interfaces drift vertically and coalesce. Figure 6(a) indicates that both $H$ - and $B$ -mergers are realized, although $H$ -mergers are more common. This feature is somewhat unusual since non-sheared double-diffusive systems are generally more susceptible to $B$ -mergers (e.g. Radko et al. Reference Radko, Bulters, Flanagan and Campin2014a ,Reference Radko, Flanagan, Stellmach and Timmermans b ). Our finding that the average layer thickness greatly increases in time suggests that the effective step height is not directly controlled by the vertical shear scale $(H)$ . This is an encouraging result. It opens prospects for developing theoretical models of the equilibrium-layer thickness without the explicit knowledge of the relevant wavelength of the basic shear, a quantity that is poorly constrained by observations or models. It should be noted that the non-dimensional layer thickness of $h\sim 20$ realized at the final stages of this experiment (figure 5 d) can be converted to dimensional values using (2.9), resulting in

which is broadly consistent with step sizes of diffusive staircases observed in the Arctic and Southern Oceans (e.g. Kelley et al. Reference Kelley, Fernando, Gargett, Tanny and Ozsoy2003).

Figure 6(b) presents an analogous experiment in which the Richardson number is increased to $Ri=100$ . The outcome is qualitatively similar. The spontaneous transition from uniform to stepped stratification starts with the formation of relatively thin layers, the scales of which are set by the periodicity of the background shear. These layers then undergo a series of mergers, which systematically increases the average layer thickness in the staircase. The time scale of the transition is, however, significantly longer for $Ri=100$ than for $Ri=10$ . The initial layering does not occur until $t=6000$ – which corresponds to a dimensional time scale of years – whereas in figure 6(a) $(Ri=10)$ layers are already visible at $t=500$ .

4 Linear stability analysis

The preliminary experiments (§ 3) revealed that dynamically and double-diffusively stable flows in the ocean can be unstable with respect to mixed thermohaline–shear instability, and its nonlinear evolution can ultimately produce layered stratification. Therefore, this section offers a more systematic description of thermohaline–shear instability by means of a Floquet-based linear stability analysis.

To examine stability properties of doubly stratified shear flows, the governing equations are first linearized with respect to the basic state (2.10), resulting in

where the velocity and pressure perturbations are defined in (2.11). The solution is sought in the following form, suggested by the Floquet theory:

where $\unicode[STIX]{x1D706}$ is the growth rate, $k$ and $l$ are the horizontal wavenumbers, and $m_{f}$ is the Floquet coefficient, which controls the fundamental wavelength in $z$ . When (4.2) is substituted into the linear system (4.1) and the individual Fourier components are collected, the stability problem reduces to the matrix form

where

and $\unicode[STIX]{x1D63C}$ is the $(12N+2)\times (12N+2)$ matrix, whose elements are functions of $k,l,m_{f},Pe,R_{\unicode[STIX]{x1D70C}},Ri,N$ .

For each set of governing parameters, the eigenvalue with the maximum real part determines the growth rate of the most rapidly amplifying mode. The spatial pattern of the perturbation is then reconstructed from the corresponding eigenvector $\unicode[STIX]{x1D743}$ :

where $\tilde{T}$ represents the temperature distribution in the fastest growing mode. The patterns of $S$ , $\boldsymbol{v}$ and $p$ in the normal modes are reconstructed in a similar manner. Extensive experimentation with the model indicated that the most rapid growth occurs when the Floquet factor is zero. This means that the fastest growing mode of instability has the same periodicity as that of the shear, and $m_{f}=0$ will be used in all subsequent calculations. Likewise, testing the performance of the model for different resolutions revealed that $N=100$ is sufficient for an accurate representation of the thermohaline–shear instabilities in the considered parameter range, and this value will be employed hereafter.

Figure 7 presents some typical experiments in which $Ri$ and $Pe$ are varied while other parameters are held constant: $(R_{\unicode[STIX]{x1D70C}},Pr,\unicode[STIX]{x1D70F})=(2,10,10^{-2})$ . The growth rate $\text{Re}(\unicode[STIX]{x1D706})$ is plotted as a function of $(k,l)$ for $(Ri,Pe)=(1,10^{4}),(1,10^{2}),(10,10^{2})$ in (a–c) respectively. In all cases, the system is stable with respect to dynamic and diffusive instabilities individually. Yet, it is apparently destabilized by the interaction between weak shear and two-component stratification (figure 7 a–c). The maximal growth rate rapidly increases with decreasing $Pe$ . For instance, the growth rate in the high- $Pe$ regime (figure 7 a) is less by two orders of magnitude than in the low- $Pe$ case (figure 7 b) for the same $Ri$ . The decrease in $Pe$ also shifts the (non-dimensional) wavelength of the fastest growing mode to much larger values. For the same $Pe$ , a reduction in shear (larger $Ri$ ) results in considerably lower growth rates (cf. figure 7 b,c). It is important to emphasize that in all cases considered, the largest growth rates were always attained for $l=0$ , which means that the wave fronts of the most rapidly amplifying perturbations are aligned perpendicular to the direction of the basic current. All of these features are consistent with the observed evolution of thermohaline–shear instabilities in DNS (§ 3).

Figure 7. The growth rates $\text{Re}(\unicode[STIX]{x1D706})$ predicted by the Floquet-based stability analysis are plotted as a function of the horizontal wavenumbers $(k,l)$ . Calculations are presented for $(Ri,Pe)=(1,10^{4}),(1,10^{2}),(10,10^{2})$ in (a–c) respectively. Only positive growth rates are shown.

In figure 8, we plot the representative spatial pattern of the most rapidly amplifying mode for $(Ri,Pe)=(1,10^{4})$ , which was reconstructed using (4.5). The amplitude of this mode was normalized by the maximal value of the temperature perturbation, and the resulting solution is denoted by $(\hat{T},{\hat{S}},\hat{u} ,{\hat{w}})$ . As observed in high- $Pe$ DNS (figures 1, 3), the unstable perturbations in the temperature and salinity fields tend to be localized to the low-shear part of the flow. Because of the symmetry of our system, there are two modes corresponding to the same (largest) growth rate: one is centred at $z=1/4$ and another, the mirror image of the first, at $z=3/4$ . In figure 9, the same calculation is performed for a low- $Pe$ case: $(Ri,Pe)=(10,10^{2})$ . In this regime, the structure of the fastest growing normal mode is substantially different. The temperature perturbation is relatively uniform at all $z$ -levels (figure 9 a). Salinity, on the other hand, tends to be localized in the high-shear zones $(z\approx 0,1/2,1)$ . In both low- $Pe$ and high- $Pe$ regimes (figures 8 and 9), the amplitude of the horizontal (vertical) velocity is largest in the high-shear (low-shear) regions.

Figure 8. The most rapidly amplifying normal mode identified by a Floquet-based stability analysis for $(R_{\unicode[STIX]{x1D70C}},Ri,Pe,\unicode[STIX]{x1D70F})=(2,1,10^{4},0.01)$ . The perturbation is normalized by the maximal temperature value. The patterns of $T$ , $S$ , $u$ and $w$ are shown in (a–d) respectively.

Figure 9. The same as figure 8 but for the low- $Pe$ case $(R_{\unicode[STIX]{x1D70C}},Ri,Pe,\unicode[STIX]{x1D70F})=(2,10,10^{2},0.01)$ .

Figure 10 presents the plot of maximal growth rates, in logarithmic coordinates, as a function of $(Pe,Ri)$ for various density ratios $(R_{\unicode[STIX]{x1D70C}})$ . A reduction in $R_{\unicode[STIX]{x1D70C}}$ has a destabilizing effect on the system. The parameter space occupied by unstable modes tends to be very wide at low density ratios ( $R_{\unicode[STIX]{x1D70C}}=1.5$ and $R_{\unicode[STIX]{x1D70C}}=2$ in figure 10 a,b) and contracts with the increasing $R_{\unicode[STIX]{x1D70C}}$ ( $R_{\unicode[STIX]{x1D70C}}=10$ and $R_{\unicode[STIX]{x1D70C}}=50$ in figure 10 c,d). Likewise, the thermohaline–shear instability is sensitive to the diffusivity ratio (figure 11). For $\unicode[STIX]{x1D70F}=0.005$ (figure 11 a) and $\unicode[STIX]{x1D70F}=0.1$ (figure 11 b), unstable modes occupy most of the parameter space explored in figure 11. Increasing $\unicode[STIX]{x1D70F}$ tends to stabilize the system. For $\unicode[STIX]{x1D70F}=1$ , the dynamics of our system simplifies to that of an effectively one-component fluid, and the maximal unstable Richardson numbers reduce to the rather unremarkable interval of $Ri<1$ . While being of marginal relevance for our primary objective (i.e. explaining the origin of diffusive staircases), it is interesting to examine the case in figure 11(d). This configuration $(\unicode[STIX]{x1D70F}=1.5)$ represents a fluid stratified in the fingering sense (the slower diffusing component is unstably stratified) but linearly stable with respect to double diffusion $(R_{\unicode[STIX]{x1D70C}}>\unicode[STIX]{x1D70F})$ . Curiously, in this regime we find no evidence of new types of instability driven by a combination of shear and different molecular diffusivities of density components. Thus, thermohaline–shear destabilization appears to be an exclusively diffusive phenomenon.

Figure 10. The decimal logarithm of the maximal growth rate $\log _{10}(\text{Re}(\unicode[STIX]{x1D706}))$ is plotted as a function of $(Ri,Pe)$ for $R_{\unicode[STIX]{x1D70C}}=1.5,2,10$ and 50 in (a–d) respectively for $\unicode[STIX]{x1D70F}=0.01$ . Only the growth rates exceeding $10^{-5}$ are shown. One can note the reduction in the intensity of thermohaline–shear instability with increasing density ratio.

Figure 11. The decimal logarithm of the maximal growth rate $\log _{10}(\text{Re}(\unicode[STIX]{x1D706}))$ is plotted as a function of $(Ri,Pe)$ for $\unicode[STIX]{x1D70F}=0.005,0.1,1$ and 1.5 in (a–d) respectively for $R_{\unicode[STIX]{x1D70C}}=2$ . Only the growth rates exceeding $10^{-5}$ are shown. One can note the reduction in the intensity of thermohaline–shear instability with increasing diffusivity ratio.

5 Physical interpretation and energetics

In this section, an attempt is made to generate some insight into the dynamics of thermohaline–shear instability. One of the key questions in this regard is the ultimate source of energy required to amplify the unstable perturbations. The answer to this question is by no means obvious. The source of energy of Kelvin–Helmholtz instability is the kinetic energy of the basic shear, whereas double diffusion is driven by the release of potential energy stored in the unstably stratified density component. To identify the energy source for thermohaline–shear instabilities, the (non-dimensional) perturbation energy equation is derived as follows. First, we subtract the temperature and salinity equations in the linearized system (4.1), multiply the result by $(S-T)$ , and integrate the result over one fundamental wavelength in $x$ , $y$ and $z$ to obtain the density variance equation:

Similarly, the kinetic energy equation is obtained by multiplying the $x$ , $y$ and $z$ momentum equations by $u$ , $v$ and $w$ respectively, adding them, and integrating the result:

Next, the density variance equation (5.1) is multiplied by $(4\unicode[STIX]{x03C0}^{2}Ri)/(R_{\unicode[STIX]{x1D70C}}-1)^{2}$ and added to the kinetic energy equation (5.2), which results in the total perturbation energy equation:

where

The $P_{dd}$ term in (5.3) represents the production of perturbation energy by diffusion of heat and salt, $P_{sh}$ is interpreted as the production of energy by shear, and $P_{visc}$ is the viscous dissipation.

Finally, the individual terms in (5.3) were evaluated for the fastest growing normal modes identified using the linear analysis in § 4. Of particular interest is the relative contribution of double-diffusive, shear-driven and viscous energy production/dissipation mechanisms. The answer turned out to be regime-dependent. Figure 12 shows the variation in the energy balance for a series of Floquet-based calculations in which $Pe$ was systematically increased, whereas all other governing parameters were kept constant: $(R_{\unicode[STIX]{x1D70C}},Ri,\unicode[STIX]{x1D70F})=(2,1,0.01)$ . The relative contribution of various processes is quantified by plotting

as a function of $Pe$ . These diagnostics highlight the stark differences between the low- $Pe$ and high- $Pe$ regimes, which were evident in the DNS (§ 3) and in the linear stability analyses (§ 4). For high Peclet numbers $(Pe>500)$ , the dominant balance is between the energy production by basic shear and its viscous dissipation. Nevertheless, double diffusion still accounts for a surprisingly large fraction of energy production $(P_{dd\,rel}=0.55{-}0.65)$ . The reason for such a high double-diffusive contribution is that much of the perturbation kinetic energy produced by shear is immediately dissipated by viscosity.

Figure 12. The energy production/dissipation balances for various values of $Pe$ . One can note the existence of two dynamically dissimilar regimes: (i) the high- $Pe$ regime $(Pe>500)$ in which the dominant balance is between the production of the perturbation energy by shear and its viscous dissipation and (ii) the low- $Pe$ regime $(Pe<500)$ in which the energy released by the unstable basic temperature gradient is transferred directly into the perturbation field, resulting in its amplification.

The dynamics of low- $Pe$ systems is dissimilar. The perturbation energy is produced almost entirely by double diffusion $(P_{dd\,rel}=1.05{-}1.1)$ . The shear-driven contribution is small $(P_{sh\,rel}=0.01{-}0.05)$ and so is the viscous dissipation $(-P_{visc\,rel}=0.05{-}0.15)$ . This result, however, should be interpreted with caution. Even though shear-induced effects are not strongly reflected in the energetics, it would be erroneous to conclude that they are largely irrelevant for the dynamics of thermohaline systems at low $Pe$ . On the contrary, the release of potential energy by double diffusion is only possible due to the catalytic role of shear. In the absence of shear, the smooth thermohaline stratification in the oceanographically relevant parameter regime would be linearly stable.

The schematic diagram in figure 13 attempts to rationalize the destabilizing tendency of shear. Before delving into the dynamics of the interaction between shear and double diffusion, let us recall the conventional interpretation of oscillatory diffusive instability. A particle initially located at the equilibrium level $(A)$ is slightly displaced (figure 13) upward into a colder and fresher region $(A_{1})$ . In the absence of diffusion and friction, the particle would oscillate at the buoyancy frequency $(N_{b})$ . However, since heat diffuses faster than salt, the particle quickly adjusts its temperature to its new environment, but largely retains its salinity. Thus, the particle is now denser than the surrounding fluid and the buoyancy force is downward. Therefore, as the particle moves downward, it gains extra kinetic energy due to the buoyancy forcing associated with the diffusion of heat. This energy gain on the way down is slightly higher than the corresponding loss of kinetic energy on the way up (the particle is colder on the downward part of the cycle) and therefore there is a small net gain of energy in each cycle. This residual gain can, in principle, drive amplifying (overstable) oscillations. However, the oscillatory nature of diffusive instability makes it highly vulnerable to friction, and therefore its parameter range (1.1) is extremely narrow.

Figure 13. Schematic diagram illustrating the mechanism of thermohaline–shear instability. While diffusive instability is inherently oscillatory, and therefore highly susceptible to viscous damping, the presence of vertical shear allows for direct instability modes.

Let us now consider the effects of vertical shear, as indicated in figure 13. In this configuration, possibilities exist not only for oscillatory but also for direct modes of instability. After the particle is displaced upward (from point $A$ to $A_{1}$ ) and its buoyancy is reduced, it can be advected by the shear into the downward-moving region $(B_{1})$ . As a result, the downward motion there is further reinforced by buoyancy forcing. Likewise, as the particle is displaced below the equilibrium level $(B_{2})$ and it becomes lighter than the surrounding fluid, it can be advected by shear into the upward-moving region $(A_{2})$ . Thus, while individual Lagrangian particles in this configuration still inevitably oscillate, there is now a possibility for the direct amplification of the velocity field at any given point in space. This dynamics results in the direct – and therefore much more efficient – transfer of energy from thermal stratification into amplifying perturbations. In the direct mode, the buoyancy forcing acts continuously to amplify the perturbation, rather than alternating between reinforcing and opposing action. The absence of fast oscillations in the Eulerian velocity field implies that the proposed mechanism is also less susceptible to viscous damping. Therefore, the parameter range for thermohaline–shear instability is much wider than for diffusive convection (figures 10 and 11), and the instability conditions are readily met in the ocean.