Effects of phase boundary and shear on diffusive instability

Abstract

We study the diffusive instability subject to a solid–liquid interface and a Kolmogorov flow using modal, non-modal analyses and energy analysis. It is found that the phase boundary has different effects on the threshold of diffusive convection for weak and strong salinity stratification. In the context of oceanography where the salinity Rayleigh number R_S is very high, the ice–water interface has negligible influence on the onset of diffusive convection. In the presence of shear, the diffusive convection for R_S < 10⁶ tends to be inhibited and with the increase of shear intensity, the oscillatory, steady convective and Kelvin–Helmholtz instabilities will be successively dominant. For R_S > 10⁶, the shear has a destabilizing effect on the diffusive convection, due to the generation of thermohaline-sheared instability found by Radko (J. Fluid Mech., vol. 805, 2016, pp. 147–170). Non-modal analysis indicates that for realistic parameters of high-latitude oceans, with the transition of the ultimate energy source of instability from the density gradient to background current, the thermohaline-shear instability is expected to transition from oscillatory to steady instability. The initial transient amplification due to double diffusion has also been studied, which is due to the generation of overstable instability at the initial phase. For R_S < 10⁶, the optimal initial condition to achieve the maximum transient growth favours longitudinal rolls. For thermohaline-shear instability, however, it favours transverse rolls and specifically, in oscillatory thermohaline-shear instability, the transient amplification can be enhanced by the shear by one order of magnitude, thus having important influence on the stability of the system.

1. Introduction

Double-diffusive convection represents a mixing process driven by two fluid components that diffuse at different rates. It is of great significance in oceanography, geophysics, astrophysics and engineering (e.g. Turner 1985; Schmitt 1994; Radko 2013; Garaud 2018). Regarding the double-diffusive convection in oceanography, the temperature field constitutes the faster diffuser while the salinity field the slower diffuser (Schmitt 1994). Away from the freezing point, the density of water is a decreasing function of temperature but an increasing function of salinity. In low-latitude regions, warm and salty waters are generally located above cold and fresh waters (Radko 2013), which makes the (sub)tropical oceans susceptible to the well-known fingering convection, i.e. salt fingers. In contrast, seawater in polar regions is stratified such that both salinity and temperature decrease upward, thus prone to diffusive convection, viz. oscillatory double-diffusive convection (ODDC). The onset conditions for instability in double-diffusive convection have been well understood in classical cases (e.g. Veronis 1965; Baines & Gill 1969; Balmforth et al. 2006). Generally, the linear instability of double-diffusive convection can take place in the form of either steady or oscillatory convection, depending on the detailed sets of parameters (Veronis 1965; Balmforth et al. 2006).

Over the past decades, sea ice in the Arctic Ocean has been observed to decrease at least by 40% in volume (Rothrock, Yu & Maykut 1999; Wadhams & Davis 2000; Deser & Teng 2008). Solar radiation and basal melting are the two main driving factors of the dramatic sea ice loss (e.g. Perovich et al. 2014). Regarding the basal melting, the thermohaline stratification lying between the cold, fresh meltwater and the warm, saline Atlantic intruding water was quite stable in the past, so that the warm Atlantic intruding water, which carries enough heat to melt all the ice in the Arctic, cannot effectively transfer its heat to the lower surface of the ice (Turner 2010). However, recent measurement of the Eurasian Basin has indicated that the Atlantic warm tongue is shoaling and the stability of the stratification layer is weakened, which allows more heat to penetrate through the halocline and enter in the lower surface of sea ice, thereby leading to the so-called ‘Atlantification’ of the Arctic Ocean (Polyakov et al. 2017). In fact, basal melting plays an even more important role in the melting of the Antarctica ice shelves, which have lost more than 50% of their mass by the basal double-diffusive process (Rignot et al. 2013; Begeman et al. 2018). As mentioned before, high-latitude oceans are generally stratified in the diffusive regime. For instance, a borehole oceanographic observation on the Ross Ice Shelf indicates that near the grounding zones, the value of the Turner angle, which is a common measure of the stratification pattern, is mainly within (−90°, −45°), thereby susceptible to diffusive convection (Begeman et al. 2018). So far, however, instability properties of the subglacial double-diffusive convection in various ocean conditions, e.g. buoyancy by double diffusion, ice melting and background current, have not been well understood. Specifically, while previous studies have demonstrated the non-trivial effect of a solid–liquid phase boundary on the stability of Rayleigh–Bénard convection (Davis, Müller & Dietsche 1984; Kim, Lee & Choi 2008; Vasil & Proctor 2011; Toppaladoddi & Wettlaufer 2019), this effect has not been examined in a subglacial diffusive convection system. In addition, despite of the systematic examination of the combined effects of the shear and double diffusion (i.e. the thermohaline-shear instability) from the perspective of asymptotic stability analysis and direct numerical simulation (DNS) (Radko 2016; Brown & Radko 2019; Radko 2019), the transient behaviour and the nature (i.e. oscillatory, steady and dynamic) of this newly found instability are still left to be resolved. These two points are reviewed individually below.

Regarding the effects of melting/freezing on the buoyancy-driven instability, let us first have a look at the results of investigations on single-component cases (that is, with the heat transport only). Davis et al. (1984) have studied the bifurcation and pattern formation of single-component buoyancy-driven convection underneath a steady solid–liquid interface. They considered a steady base state by assuming that the heat entering and leaving the interface balances each other. It is found that at the onset of flow motion, the solid–liquid interface behaves like a planar rigid solid of imperfect thermal conductivity. With the increase of solid depth, the conductivity effectively decreases so that the critical Rayleigh number Ra_c and critical wavenumber k_c both decrease. For sufficiently large solid thickness, (Ra_c, k_c) asymptote to (1493, 2.815). Recently, this work has been extended by Toppaladoddi & Wettlaufer (2019) to study the effect of low-Péclet-number Couette flow. It is found that the imposed shear tends to stabilize the system by inhibiting the vertical heat flux. Vasil & Proctor (2011) studied the instability of Rayleigh–Bénard convection subject to a different solid–liquid interface. They considered a Stefan system with an isothermal solid, in which the liquid layer will slowly deepen as time progresses. At the large-Stefan-number limit in which the melting is very slow, the phase interface acts as a fixed-flux boundary, which results in the convection threshold decreasing from (1707.763, 3.116) (Chandrasekhar 1961) to (1295.78, 2.552). The decrease of convection threshold due to the phase boundary has also been observed by Kim et al. (2008) in a fast-melting system.

In two-component cases, few attempts have been made to study the double-diffusive instability in the presence of an ice–water interface. In places such as Greenland (e.g. McDougall 1983; McPhee 2008; Chandler et al. 2013) and the Arctic (e.g. Polashenski, Perovich & Courville 2012; Polashenski et al. 2017), strong solar radiation in summer is able to cause the seasonal melting/freezing at the surface of ice shelves. Martin & Kauffman (1974) studied the dynamics of under-ice melt ponds resulting from the aforementioned process. They focused on the ice growth (the so-called false bottom) in the system, which was found to go through three different stages: nucleation, lateral growth and upward migration due to the temperature and salinity gradient. The formation of a false bottom is closely correlated with the double diffusion. Recently, Hirata, Goyeau & Gobin (2012) investigated the onset of double-diffusive convection in the under-ice melt ponds. They considered the water density inversion and showed that the depth of the ice matrix plays a key role in the stability of the under-ice melt layer. Particularly, increasing the depth of the ice layer will decrease the convection threshold, which is consistent with the single-component cases (Davis et al. 1984; Toppaladoddi & Wettlaufer 2019). However, the depth of the melt layer has a marginal influence, due to the compensation between stabilizing salinity gradient and destabilizing temperature gradient which simultaneously change when increasing the liquid depth.

Regarding the effect of background current on double-diffusive instability, much attention has been paid to salt fingers in the early stage. Linden (1974) studied the effect of a steady Couette flow on fingering convection between two planes. Results showed that the shear affects the fingering instability in a way similar to thermally stratified shear flow. For instance, longitudinal rolls are the preferred mode and unaffected by the shear. In addition, the instabilities of both transverse and oblique rolls tend to be inhibited by the shear since the mean flow can absorb energy from the spanwise perturbations. Recently, Konopliv, Lesshafft & Meiburg (2018) performed a transient growth analysis of a similar problem. Their results showed that even though the Couette flow is vertically unbounded, it can also damp the fingering instability. Most importantly, they detected the transient growth of fingering convection which exists regardless of whether there is a shear or not. This transient growth is attributed to a mechanism relying on the density stratification, different from the Orr mechanism or lift-up mechanism which relies on shear (Jerome, Chomaz & Huerre 2012). It should be noted that Linden (1974) and Konopliv et al. (2018) were not concerned with the shear-driven instability in their studies. In fact, the basic shear itself may yield dynamic instability such as Kelvin–Helmholtz instability, which can interact with double-diffusive stratification and thus lead to the formation of thermohaline staircases (Symth & Kimura 2007). Radko & Stern (2011) studied the fingering instability in the presence of a shear with sinusoidal velocity profile (namely the Kolmogorov flow). By assuming triply periodic boundary conditions, they identified two kinds of instability that dominate the system. For small values of Richardson number R_i, the classical Kelvin–Helmholtz instability is the fastest growing mode; for large values of R_i, however, the Kelvin–Helmholtz instability will be superseded by the thermohaline-shear instability, a newly found oscillatory instability that bears strong resemblance to the collective instability observed in thermally stratified Kolmogorov flow (Balmforth & Young 2002).

Recently, to understand the thermohaline layering of diffusively stable oceans in high-latitude regions, Radko (2016) investigated the diffusive instability in the presence of the Kolmogorov flow. Both results of DNS and linear stability analysis indicated that the interaction between the shear and density stratification can surprisingly destabilize the system that is individually stable to Kelvin–Helmholtz and diffusive instabilities. It is clear that as with Kelvin–Helmholtz instability, the fastest growing mode in the thermohaline-shear instability is the transverse rolls, which is in contrast to the fingering cases (Linden 1974; Symth & Kimura 2007; Radko & Stern 2011; Konopliv et al. 2018). Radko (2016) performed an energy analysis to suggest that there exist two regimes of the energy balance of thermohaline-shear instability. Different from the low-Pe regime (Pe is the Péclet number) in which the energy production by shear is negligible compared with that by temperature gradient, the dominant energy balance in the high-Pe regime is between the energy produced by shear and dissipated by viscosity. It is noted that in this seminal work, there are three questions remaining to be resolved, which partially motivates us to further study the sheared diffusive convection by analysing the energy budgets with the help of non-modal analysis.

First, in the high-Pe regime, the author attributed the main energy source of amplifying the disturbance to the unstable temperature gradient, which is definitely true in the low-Pe regime only. This might be because in the energy analysis, the author fixed all controlling parameters except Pe, which leads to the overestimation of energy contribution of heat diffusion (especially in the high-Pe regime). Note that in the thermohaline-shear instability, the critical temperature gradient is varied when varying Pe (specifically, we will show that in fact the critical Rayleigh number is a rapidly decreasing function of Pe).

Second, Radko (2016) rationalized the catalytic role of shear in the energy transfer from the density gradient to the perturbed hydrodynamic field by hypothesizing that the existence of a shear increases the residual energy gain of a particle in an oscillation cycle. However, this hypothesis requires that the thermohaline-shear instability is within the oscillatory rather than steady instability regime, which will be addressed in this work.

Finally, the analysis of Radko (2016) focused only on the asymptotic growth of the perturbation, while as mentioned before, it has been confirmed by Konopliv et al. (2018) that the initial transient growth can occur in double-diffusive convection even in the absence of a shear. Considering that the short-term disturbance behaviour can sometimes overdraw the long-term one determined by the eigenvalue analysis (Schmid & Henningson 2001), it is necessary to inspect the initial transient growth of thermohaline-shear instability.

Nevertheless, as a plausible mechanism of the generation of staircases in the Arctic, thermohaline-shear instability has been discussed by Garaud (2017) in wider geophysical prospects and also studied numerically with a linear velocity profile (Yang et al. 2022). In addition, to consider the effects of internal gravity waves, this instability has been investigated in the case of shear with time-dependent magnitude (Brown & Radko 2019; Radko 2019). Specifically, a temporal record of the quadratic perturbation norm (not the non-modal analysis) by Radko (2019) indicated that the overall exponential growth is modulated by small-scale oscillations. Since the small-scale oscillations acquire the angular frequency twice of the time-dependent mean shear, it is believed that this oscillatory behaviour is inherited from the base flow so that it cannot help to identify whether the thermohaline-shear instability is oscillatory or steady.

As shown above, few works have studied the effects of an ice–water interface on the linear instability of diffusive convection. However, previous studies of effects of shear on diffusive instability focused mainly on the thermohaline-shear instability, so a more comprehensive investigation in this regard is necessary. As for the thermohaline-shear instability, there are still several essential problems left to be resolved. This leads to the aims of the present study: (i) investigate the diffusive instability in the presence of an ice–water interface; (ii) study the transition from oscillatory to steady, and to dynamic instability; and (iii) examine the transient growth of diffusive convection in the absence/presence of shear. To this end, in addition to the modal stability analysis, non-modal stability analysis and energy analysis are also adopted to improve our understanding of the problems. This paper is organized as follows. In § 2, we formulate the problem. In § 3, we derive the linearized equations and provide a systematic code validation. The results are present in § 4 and conclusions are drawn in § 5.

2. Problem formulation

This work extends on that of Toppaladoddi & Wettlaufer (2019) to account for salinity stratification. As sketched in figure 1, we consider an idealized system which consists of an infinite horizontal liquid layer of thickness h ₀ sandwiched between a lower fixed plane and an upper solid layer. The upper solid layer corresponds to the ice shelf of thickness d ₀. A Cartesian coordinate system (x, y, z) is centred on the bottom of the channel with the y-axis pointing vertically upwards. The base of the water is held at the temperature T_b higher than that at the upper surface of ice shelf T_u. The melting/freezing temperature of ice, i.e. T_i, is such that T_u < T_i < T_b, which allows the phase interface to move upwards or downwards, depending on the heat fluxes towards and away from the interface. When the heat flux across the liquid layer is larger than those across the ice layer, the system undergoes melting, which produces fresh water near the ice–ocean interface so that the ambient salinity is stratified as with the temperature in the liquid. The configuration considered here, that both the temperature and salinity are increasing functions of the water depth, is based on the observations of Begeman et al. (2018). The salinities of lower and upper surfaces of the liquid layer are denoted by S_b and S_i. Following Hewitt (2020), we assume the ice void of salinity for convenience. In addition, as with Radko (2016), we impose a Kolmogorov flow in the liquid layer to model the effect of ambient current, with U ₀ being the maximum of the velocity profile.

Figure 1. Schematic description of the problem.

2.1 Governing equations

In the liquid, the governing equations for sheared double-diffusive process are (e.g. Radko 2013)

(2.1a)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U} = 0,\end{gather}

(2.1b)\begin{gather}\frac{{\partial \boldsymbol{U}}}{{\partial t}} + \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla U} =- \frac{{\boldsymbol{\nabla }P}}{{{\rho _0}}} + \upsilon {\nabla ^2}\boldsymbol{U} + g\frac{{\rho - {\rho _0}}}{{{\rho _0}}}\boldsymbol{y},\end{gather}

(2.1c)\begin{gather}\frac{{\partial T}}{{\partial t}} + \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla }T = {k_T}{\nabla ^2}T,\end{gather}

(2.1d)\begin{gather}\frac{{\partial S}}{{\partial t}} + \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla }S = {k_S}{\nabla ^2}S,\end{gather}

where U = (U, V, W) is the velocity field, P the dynamic pressure, g the acceleration due to gravity, ρ (ρ ₀) the (reference) density, υ the kinematic viscosity and y the unit vector in the y direction. Here, (T, S) represent the temperature and salinity, with their molecular diffusivities denoted by (k_T, k_S). We assume the linear equation of state for simplicity (Radko 2013):

(2.2)\begin{equation}\rho = {\rho _0}[1 + {\alpha _S}(S - {S_0}) - {\alpha _T}(T - {T_0})],\end{equation}

where (α_T, α_S) are the constant expansion/contraction coefficients, and (T ₀, S ₀) are the reference temperature and salinity.

The heat transport in the solid is governed by thermal diffusion, viz.

(2.3)\begin{equation}\frac{{\partial {T_s}}}{{\partial t}} = {k_T}{\nabla ^2}{T_s},\end{equation}

where T_s is the temperature in the solid. Note that the thermal diffusivities in the solid and the liquid are identical for simplicity.

The ice–water interface is at the temperature of melting/freezing, and its evolution is governed by the Stefan condition, i.e.

(2.4a,b)\begin{equation}T = {T_s} = {T_i},\quad {\rho _0}L\frac{{\partial h}}{{\partial t}} = \boldsymbol{n}\boldsymbol{\cdot }\kappa (\nabla {T_s} - \nabla T),\end{equation}

where L is the latent heat of fusion, h = h(x, z) the height of the interface, n the unit vector normal to the interface and κ the thermal conductivity.

In addition to the Stefan condition, the remaining boundary conditions are

(2.5a)\begin{gather}\boldsymbol{U} = {\textbf{0}},\quad T = {T_b},\quad S = {S_b}\quad \textrm{at}\;y = 0,\end{gather}

(2.5b)\begin{gather}\boldsymbol{U} = {\textbf{0}},\quad S = {S_i}\quad \textrm{for}\;y = {h_0},\end{gather}

(2.5c)\begin{gather}{T_s} = {T_u}\quad {\rm at}\;y = {h_0} + {d_0}.\end{gather}

To place the above formulations in dimensionless form, we introduce the following scales:

(2.6)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {x,y,z \sim {h_0},\quad u,v,w \sim {U_0},\quad t \sim \dfrac{{{h_0}}}{{{U_0}}},}\\ {T - {T_0} \sim \varDelta T = {T_b} - {T_i},\quad S - {S_0} \sim \varDelta S = {S_b} - {S_i},\quad P \sim \rho U_0^2.} \end{array}} \right\}\end{equation}

Here we have T ₀ = T_i. These scalings give rise to the following dimensionless numbers, viz.

(2.7)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {Ra = \dfrac{{g{\alpha_T}\varDelta Th_0^3}}{{\upsilon {k_T}}},\quad {R_S} = \dfrac{{g{\alpha_S}\varDelta Sh_0^3}}{{\upsilon {k_S}}},\quad Re = \dfrac{{{U_0}{h_0}}}{\upsilon },}\\ {Pr = \dfrac{\upsilon }{{{\alpha_T}}},\quad \beta = \dfrac{{{k_S}}}{{{k_T}}},\quad St = \dfrac{L}{{{c_p}\varDelta T}},\quad \varLambda = \dfrac{{{T_i} - {T_u}}}{{\varDelta T}}.} \end{array}} \right\}\end{equation}

Here, Ra and R_S are the thermal and the salinity Rayleigh numbers which measure the buoyancy due to temperature and salinity gradients, respectively. The Reynolds number Re represents the strength of the shear flow. The Prandtl number Pr compares kinematic viscosity to thermal diffusivity. Additionally, β is the diffusivity ratio of salt and heat. Since we are concerned here with the diffusive convection of high-latitude oceans, we will consider Pr and β to be fixed at Pr = 10 and β = 0.01 throughout all calculations, unless otherwise specified. The Stefan number St denotes the ratio between the latent heat and specific heat (c_p), and $\varLambda$ is the ratio of the temperature differences across the liquid and the ice layers. Maintaining the pre-scaled notation, the governing equations become

(2.8a)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U} = 0,\end{gather}

(2.8b)\begin{gather}\frac{{\partial \boldsymbol{U}}}{{\partial t}} + \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla U} =- \boldsymbol{\nabla }P + \frac{1}{{Re}}{\nabla ^2}\boldsymbol{U} + \frac{{Ra}}{{Pr \cdot R{e^2}}}\boldsymbol{\cdot }T\boldsymbol{\cdot }\boldsymbol{y} - \frac{{{R_S} \cdot \beta }}{{Pr \cdot R{e^2}}}\boldsymbol{\cdot }S\boldsymbol{\cdot }\boldsymbol{y},\end{gather}

(2.8c)\begin{gather}\frac{{\partial T}}{{\partial t}} + \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla }T = \frac{1}{{Pr \cdot Re}}{\nabla ^2}T,\end{gather}

(2.8d)\begin{gather}\frac{{\partial S}}{{\partial t}} + \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla }S = \frac{\beta }{{Pr \cdot Re}}{\nabla ^2}S,\end{gather}

(2.8e)\begin{gather}\frac{{\partial {T_s}}}{{\partial t}} = \frac{1}{{Pr \cdot Re}}{\nabla ^2}{T_s},\end{gather}

with accompanying boundary conditions

(2.9a)\begin{gather}\boldsymbol{U} = {\textbf{0}},\quad S = 1,\quad T = 1\quad \textrm{at}\;y = 0,\end{gather}

(2.9b)\begin{gather}\boldsymbol{U} = {\textbf{0}},\quad S = \textrm{0},\quad T = {T_s} = 0,\quad St\frac{{\partial h}}{{\partial t}} = \boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{\nabla }{T_s} - \boldsymbol{\nabla }T),\quad \textrm{for}\;y = 1,\end{gather}

(2.9c)\begin{gather}{T_s} =- \varLambda \quad \textrm{at}\;y = 1 + {d_0}.\end{gather}

3. Linear stability equations

In this section, we derive the stability equations for infinitesimal disturbances in the present system. We first focus on the dynamics in the liquid. Based on the linear stability theory (Schmid & Henningson 2001), all variables are separated into the basic state and a departure from it, i.e.

(3.1)\begin{equation}\boldsymbol{U} = \bar{\boldsymbol{U}} + \boldsymbol{u},\quad P = \bar{P} + p,\quad T = \bar{T} + \theta ,\quad {T_s} = \overline {{T_s}} + {\theta _s},\quad S = \bar{S} + s,\end{equation}

where u = (u, v, w), p, θ, θ_s and s represent the velocity vector, pressure, temperature and salinity in the perturbed state, while $\bar{\boldsymbol{U}}$, $\bar{P}$, $\bar{T}$, $\overline {{T_s}}$ and $\bar{S}$ are the corresponding base-state variables. Inserting the decompositions into (2.8a–d), subtracting the equation for the basic state and omitting the nonlinear terms, we obtain

(3.2a)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0,\end{gather}

(3.2b)\begin{gather}\frac{{\partial \boldsymbol{u}}}{{\partial t}} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\bar{\boldsymbol{U}} + (\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} =- \boldsymbol{\nabla }p + \frac{1}{{Re}}{\nabla ^2}\boldsymbol{u} + \frac{{Ra}}{{R{e^2} \cdot Pr}}\boldsymbol{\cdot }\theta \boldsymbol{\cdot }\boldsymbol{y} - \frac{{{R_S} \cdot \beta }}{{R{e^2} \cdot Pr}}\boldsymbol{\cdot }s\boldsymbol{\cdot }\boldsymbol{y},\end{gather}

(3.2c)\begin{gather}\frac{{\partial \theta }}{{\partial t}} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\bar{T} + (\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{\nabla })\theta = \frac{1}{{Re \cdot Pr}}{\nabla ^2}\theta ,\end{gather}

(3.2d)\begin{gather}\frac{{\partial s}}{{\partial t}} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\bar{S} + (\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{\nabla })s = \frac{\beta }{{Re \cdot Pr}}{\nabla ^2}s.\end{gather}

The boundary conditions for the velocity and salinity fluctuations read (the Stefan condition will be discussed later)

(3.3)\begin{equation}\boldsymbol{u} = 0\quad \textrm{and}\quad s = 0\quad \textrm{for}\;y = 0\;\textrm{and}\;y = 1.\end{equation}

3.1 Base-state solutions

In this study, we consider a steady and horizontally homogeneous base state in which the ice–water interface is unmovable unless the vertical flow motion sets in. This steady state is reasonable due to the fact that, given that the solid is conducting, any melting process will eventually reach an equilibrium state, as long as the heat fluxes entering and leaving the interface balance (e.g. Purseed, Favier & Duchemin 2020). Therefore, in the base state, we have (Toppaladoddi & Wettlaufer 2019)

(3.4)\begin{equation}\frac{{\textrm{d}\bar{T}}}{{\textrm{d}y}} = \frac{{\textrm{d}\overline {{T_s}} }}{{\textrm{d}y}}\quad \textrm{at}\;y = 1.\end{equation}

This leads to

(3.5)\begin{equation}{d_0} = \varLambda ,\end{equation}

which indicates that the thicknesses of the liquid and solid layers in the equilibrium state are determined by the values of T_b, T_i and T_u (Davis et al. 1984).

The base-state temperature profile is given by (Davis et al. 1984; Toppaladoddi & Wettlaufer 2019)

(3.6a,b)\begin{equation}\bar{T}(y) = 1 - y,\quad 0 \le y \le 1\quad \textrm{and}\quad \overline {{T_s}} (y) = 1 - y,\quad 1 \le y \le 1 + {d_0}.\end{equation}

Regarding the base state of salinity, we have

(3.6c)\begin{equation}\bar{S}(y) = 1 - y,\quad 0 \le y \le 1.\end{equation}

For the base flow, we are mainly interested in the hydrostatic base state and Kolmogorov flow (Radko 2016), which are respectively given by

(3.6e)\begin{gather}\bar{\boldsymbol{U}}(y) = 0\quad \textrm{for}\;\textrm{hydrostatic}\;\textrm{state},\end{gather}

(3.6f)\begin{gather}\bar{U}(y) =- \textrm{sin}(2{\rm \pi} y),\quad \bar{V} = \bar{W} = 0\quad \textrm{for}\;\textrm{Kolmogorov}\;\textrm{flow}\;\textrm{with}\;0 \le y \le 1.\end{gather}

3.2 Matrix representation

After the given base state is substituted in the linearized equations and the pressure term is eliminated from the obtained equations, the linearized system can be rewritten as

(3.7a)\begin{gather}\begin{array}{ccccc} \dfrac{{\partial {\nabla ^2}v}}{{\partial t}} & = \left( { - \bar{U}\dfrac{\partial }{{\partial x}}{\nabla^2} + \dfrac{{{\textrm{d}^2}\bar{U}}}{{\textrm{d}{y^2}}}\dfrac{\partial }{{\partial x}} + \dfrac{1}{{Re}}{\nabla^4}} \right)v + \dfrac{{Ra}}{{R{e^2} \cdot Pr}}\left( {{\nabla^2} - \dfrac{{{\partial^2}}}{{\partial {y^2}}}} \right)\theta \\ & \quad - \dfrac{{{R_S} \cdot \beta }}{{R{e^2} \cdot Pr}}\left( {{\nabla^2} - \dfrac{{{\partial^2}}}{{\partial {y^2}}}} \right)s, \end{array}\end{gather}

(3.7b)\begin{gather}\frac{{\partial \eta }}{{\partial t}} =- \frac{{\textrm{d}\bar{U}}}{{\textrm{d}y}}\frac{\partial }{{\partial z}}v + \left( { - \bar{U}\frac{\partial }{{\partial x}} + \frac{1}{{Re}}{\nabla^2}} \right)\eta ,\end{gather}

(3.7c)\begin{gather}\frac{{\partial \theta }}{{\partial t}} =- \frac{{\textrm{d}\bar{T}}}{{\textrm{d}y}}v + \left( { - \bar{U}\frac{\partial }{{\partial x}} + \frac{1}{{Re \cdot Pr}}{\nabla^2}} \right)\theta ,\end{gather}

(3.7d)\begin{gather}\frac{{\partial s}}{{\partial t}} =- \frac{{\textrm{d}\bar{S}}}{{\textrm{d}y}}v + \left( { - \bar{U}\frac{\partial }{{\partial x}} + \frac{\beta }{{Re \cdot Pr}}{\nabla^2}} \right)s,\end{gather}

where η is the vorticity in the y-direction, defined by

(3.8)\begin{equation}\eta = \frac{{\partial u}}{{\partial z}} - \frac{{\partial w}}{{\partial x}}.\end{equation}

For clarity, the linearized system (3.7) can be recast in matrix notation, viz.

(3.9)\begin{equation}\boldsymbol{\mathsf{A}}\frac{{\partial \boldsymbol{\zeta }}}{{\partial t}} = \boldsymbol{\mathsf{B}}\zeta \quad \Rightarrow \quad \frac{{\partial \boldsymbol{\zeta }}}{{\partial t}} = \boldsymbol{\mathsf{L}}\zeta ,\end{equation}

where ζ = (v, η, θ, s)^T and $\boldsymbol{\mathsf{L}}$ = $\boldsymbol{\mathsf{A}}$⁻¹$\boldsymbol{\mathsf{B}}$ is the linearized Navier–Stokes operator.

We seek solutions of (3.9) in the following wave-like form:

(3.10)\begin{equation}f(x,y,z,t) = {f^\ast }(y)\ \textrm{exp}(\textrm{i}kx + \textrm{i}mz - \textrm{i}\omega t),\end{equation}

where f corresponds to (v, η, θ, s), f* = (v*, η*, θ*, s*) is the shape function, k and m the streamwise and spanwise wavenumbers, and ω the circular frequency of the perturbation. Note that ω is complex valued, with the real part ω_r being the phase speed while the imaginary part ω_i the growth rate. Substituting (3.10) into (3.9), we arrive at an eigenvalue problem,

(3.11)\begin{equation}- \textrm{i}\omega {\boldsymbol{\zeta }^\ast } = \boldsymbol{\mathsf{L}}{\boldsymbol{\zeta }^\ast },\end{equation}

where ζ* is the eigenvector of the eigenvalue −iω.

The boundary conditions for the amplitudes of the perturbed velocity, vorticity and salinity are

(3.12a)\begin{equation}{v^\ast } = \frac{{\partial {v^\ast }}}{{\partial y}} = {\eta ^\ast } = {s^\ast } = 0\quad \textrm{for}\;y = 0\;\textrm{and}\;y = 1.\end{equation}

The boundary condition for the perturbed temperature at the interface is derived from the Stefan condition (2.9b), which has been discussed by Davis et al. (1984). It is found that the interface acts as an imperfect thermal conductor such that

(3.12b)\begin{equation}\textrm{D}{\theta ^\ast } + {\theta ^\ast }\gamma \ \textrm{coth}(\gamma \varLambda ) = 0\quad \textrm{for}\;y = 1,\end{equation}

where γ = (k ²+m ²)^1/2 is the total horizontal wavenumber. This condition is consistent with that of Nield (1968) who studied the buoyancy-driven instability underlying a solid layer of finite conductivity and finite thickness. Therefore, we in fact consider the presence of an ice–water interface which is steady before the onset of convective rolls, instead of a melting or freezing boundary. However, while the modified Stefan condition (3.12b) is independent of St, we still assume that St is large enough following Toppaladoddi & Wettlaufer (2019) and the possibility of large enough St has been discussed by Vasil & Proctor (2011). The condition for θ* at the lower boundary is

(3.12c)\begin{equation}{\theta ^\ast } = 0\quad \textrm{for}\;y = 0.\end{equation}

Regarding the heat diffusion in the solid, linearizing (2.8e) as with the fluid system leads to

(3.13)\begin{equation}({\textrm{D}^2} - {\gamma ^2})\theta _s^\ast= 0,\end{equation}

where θ_s* is the amplitude of the perturbed temperature in the solid. The corresponding boundary conditions are given by

(3.14a,b)\begin{equation}\theta _s^\ast= \theta \quad \textrm{at}\;y = 1\quad \textrm{and}\quad \theta _s^\ast= 0\quad \textrm{at}\;y = 1 + \varLambda .\end{equation}

Davis et al. (1984) have given the analytical solution to (3.13), i.e.

(3.15)\begin{equation}\theta _s^\ast= \frac{{\textrm{sinh}[\gamma (1 + \varLambda - y)]}}{{\textrm{sinh}(\gamma \varLambda )}} \cdot {\theta ^\ast }(y = 1)\quad \textrm{for}\;1 \le y \le 1 + \varLambda ,\end{equation}

which is easily solved once we determine the perturbed liquid temperature at the interface.

3.3 Energy norm

The modal stability analysis introduced in § 3.2 predicts the conditions for the exponentially growing instability of a linear system, but fails to examine its transient behaviour in a finite time interval. The transient behaviour is because the non-orthogonality of the eigenvectors in the presence of a basic shear might cause a non-trivial transient energy growth before the eventually exponential growth or decay of the disturbance (Schmid & Henningson 2001). To give insight into the transient behaviour of the present system, we resort to the non-modal stability analysis (Schmid 2007). To this end, we define the dimensionless perturbation energy as follows:

(3.16)\begin{equation}\int_\varOmega {E\,\textrm{d}V} = \int_\varOmega {({e_k} + {e_\theta } + {e_S})\,\textrm{d}V} = \int_\varOmega {[{c_1}({u^2} + {v^2} + {w^2}) + {c_2}{\theta ^2} + {c_3}{s^2}]\,\textrm{d}V,}\end{equation}

where E represents the total perturbation energy, e_k the perturbed kinetic energy, and e_θ and e_S the potential energies of the disturbance with respective to temperature and salinity fields, respectively (Konopliv et al. 2018). Since the system is horizontally periodic, the control volume $\varOmega$ is defined as 0 ≤ y ≤ 1 so that dV = dy. The definitions of perturbation energy components on the right-hand side of (3.16) follow Konopliv et al. (2018), with the accompanying coefficients given by c ₁ = 1, c ₂ = Ra/(Re ²Pr) and c ₃ = R_Sβ/(Re ²Pr).

Since we describe the system in terms of v−η−θ−s, (3.16) is rewritten as

(3.17a)\begin{equation}\int_\varOmega {E\,\textrm{d}V} = \tfrac{1}{4}\int {{\boldsymbol{\zeta }^{{\ast ^{\dagger} }}}\,\textrm{diag}({\boldsymbol{\mathsf{M}}_1},{\boldsymbol{\mathsf{M}}_2},{\boldsymbol{\mathsf{M}}_3},{\boldsymbol{\mathsf{M}}_4}){\boldsymbol{\zeta }^\mathrm{\ast }}\,\textrm{d}y} = \tfrac{1}{4}\int {{\boldsymbol{\zeta }^{{\ast ^{\dagger} }}}\boldsymbol{\mathsf{M}}{\boldsymbol{\zeta }^\mathrm{\ast }}\,\textrm{d}y} ,\end{equation}

with

(3.17b–e)\begin{equation}{\boldsymbol{\mathsf{M}}_1} = \boldsymbol{\mathsf{I}} + \frac{1}{{{\gamma ^2}}}\boldsymbol{\mathsf{D}}_1^{\dagger} {\boldsymbol{\mathsf{D}}_1},\quad {\boldsymbol{\mathsf{M}}_2} = \frac{1}{{{\gamma ^2}}}\boldsymbol{\mathsf{I}},\quad {\boldsymbol{\mathsf{M}}_3} = \frac{{Ra}}{{R{e^2}Pr}}\boldsymbol{\mathsf{I}},\quad {\boldsymbol{\mathsf{M}}_4} = \frac{{{R_S}\beta }}{{R{e^2}Pr}}\boldsymbol{\mathsf{I}},\end{equation}

where the superscript ^† denotes the complex conjugate and $\boldsymbol{\mathsf{D}}$₁ the first-order derivative with respective to y. Since the weight matrix is positive definite, we apply a Cholesky decomposition to it, i.e. $M=\boldsymbol{\mathsf{F}}$^†$\boldsymbol{\mathsf{F}}$ with ξ* = $\boldsymbol{\mathsf{F}}$ζ*, which leads to (Zhang et al. 2015)

(3.18)\begin{equation}\int_\varOmega {E\,\textrm{d}V} = \tfrac{1}{2}\int {{\boldsymbol{\zeta }^{{\ast ^{\dagger} }}}\boldsymbol{\mathsf{M}}{\boldsymbol{\zeta }^\mathrm{\ast }}\,\textrm{d}y} = \tfrac{1}{2}\int {{{\boldsymbol{\zeta} }^{{\ast ^{\dagger} }}}{\boldsymbol{\mathsf{F}}^{\dagger} }\boldsymbol{\mathsf{F}}{\boldsymbol{\zeta }^\mathrm{\ast }}\,\textrm{d}y} = \tfrac{1}{2}\int {{\boldsymbol{\xi }^{{\ast ^{\dagger} }}}{\boldsymbol{\xi }^\mathrm{\ast }}\,\textrm{d}y} ={\parallel} {\boldsymbol{\xi }^\mathrm{\ast }}_2\parallel ,\end{equation}

where ||ξ*|| is the L ₂-norm of ξ*. The transient growth is defined as the maximum possible growth for all initial conditions $\boldsymbol{\xi }_0^\ast $, i.e.

(3.19)\begin{equation}G(t) = \mathop {\max }\limits_{\boldsymbol{\xi }_0^\ast } \frac{{||{\boldsymbol{\xi }^\ast }(t)|{|_2}}}{{||{\boldsymbol{\xi }^\ast }(0)|{|_2}}}.\end{equation}

3.4 Code validation

In this study, we solve the numerical problem (3.11) using the spectral method based on collocation points as the root of Chebyshev polynomials (Weideman & Reddy 2000). This method has also been used in our previous work (Zhang et al. 2015).

First, to test the performance of our stability code for various grid resolutions, we consider a case in the diffusive instability regime. Without loss of generality, the combined effects of phase boundary and low-Re Kolmogorov flow are considered. As shown in table 1, N = 250 is sufficient to converge the most unstable mode at Ra = 1.15 × 1707.763, R_S = 10 × 1707.763, Re = 0.1, k = 3.2, m = 0 and $\varLambda$ = 1.

Table 1. Resolution check for diffusive convection in the presence of phase boundary and Kolmogorov flow at Ra = 1.15 × 1707.763, R_S = 10 × 1707.763, Re = 0.1, k = 3.2, m = 0, $\varLambda$ = 1. In the table, ω_r and ω_i are the real and imaginary components of the frequency of the most unstable mode.

Second, we verify our model by examining the convection threshold of Rayleigh–Bénard convection in the presence of a solid–liquid interface. In this case, the basic shear and salinity stratification are neglected. As shown in figure 2, the critical Rayleigh number Ra_c and the critical wavenumber k_c are the decreasing functions of the thickness of the solid layer $\varLambda$. Our calculations are in good agreement with those of Davis et al. (1984) and Toppaladoddi & Wettlaufer (2019).

Figure 2. Linear stability of Rayleigh–Bénard convection subject to phase boundary. (a) Critical Rayleigh number Ra_c as a function of $\varLambda$ and (b) critical wavenumber k_c as a function of $\varLambda$.

Finally, we also reproduce the marginal stability curves of double-diffusive convection for validation in the absence of the phase boundary and shear. For clarity, here we assume that $R_T^\ast $ (respectively $R_S^\ast $) represents the controlling parameter of the component that is unstably (respectively stably) stratified. As was done by Balmforth et al. (2006), we summarize the linear stability of double-diffusive convection from Veronis (1965) as follows:

(3.20a)\begin{gather}\textrm{steady:}\quad R_T^\ast= R_S^\ast+ R{a_{c0}}\quad \textrm{if}\;{\sigma ^2}\frac{{{\beta ^\ast } + Pr}}{{1 + Pr}} < {\beta ^\ast } < 1\;\textrm{or}\;{\beta ^\ast } > 1,\end{gather}

(3.20b)\begin{gather}\textrm{oscillatory:}\quad R_T^\ast= ({\beta ^\ast } + Pr\textrm{)}\left( {\frac{{{\beta^\ast }R_S^\ast }}{{1 + Pr}} + \frac{{R{a_{c0}}(1 + {\beta^\ast })}}{{Pr}}} \right)\notag\\ \textrm{if}\;{\sigma ^2}\frac{{{\beta ^\ast } + Pr}}{{1 + Pr}} > {\beta ^\ast }\;\textrm{and}\;{\beta ^\ast } < 1.\end{gather}

Here, ${\sigma ^2} = R_S^\ast{/}R_T^\ast $ and β* is the diffusivity ratio between the stably stratified component and unstably stratified component. Specifically, β* > 1 corresponds to the fingering convection in which $R_T^\ast $ represents R_S, whereas β* < 1 corresponds to the diffusive convection in which $R_T^\ast $ represents Ra. It is clear that fingering instability always arises as a steady type of motion whereas in diffusive cases, both steady convective instability and oscillatory instability may appear, with the transition occurring at σ ² = β*(1 + Pr)/(β* + Pr), as shown in (3.20). In addition, assuming that the boundary is perfectly conducting, one has Ra_{c 0} = 1707.763 for the no-slip condition while Ra_{c 0} = 657.511 for the stress-free boundary (Chandrasekhar 1961). To highlight the effect of the boundary condition, we consider the no-slip condition in this study. As shown in figure 3, our results are consistent with Veronis (1965).

Figure 3. Linear stability boundaries of (a) fingering convection and (b) oscillatory convection with no-slip and perfectly conducting boundaries. Note that R_{Sc 0} = Ra_{c 0} = 1707.763.

4. Results of stability analysis

4.1 The classical diffusive convection

While the linear instability of classical diffusive convection (i.e. that in the absence of the phase boundary and shear) has been well understood, it is necessary to study this problem from the perspective of transient energy evolution before proceeding to the investigation on the effects of phase boundary and shear. As mentioned before, diffusive convection can bifurcate into steady convection for weak salinity stratification or into overstable motions for intense salinity stratification. These two typical cases are illustrated in figures 4 and 5, in which we show the eigenspectra in the modal analysis, the transient energy growths in the non-modal analysis and the perturbed salinity fields. Note that when we consider the hydrostatic base state, we set m = 0. In the steady instability regime, all the instability modes are located at the core mode which corresponds to a zero phase speed. The transient energy evolves exponentially with time after an initial rapid growth. The reasons for the initial rapid growth has been examined by Konopliv et al. (2018) in the fingering regime, which is due to the interaction between double diffusion and velocity fluctuation field, rather than due to the Orr mechanism or lift-up mechanism. The situation in the oscillatory instability regime is dissimilar. As shown in figure 5(a), there exist two most unstable modes lying on the either side of the core mode. They acquire the identical growth rate but the opposite phase speed. This is because the overstable motions, a standing oscillation in the form of internal gravity waves, are characterized by the periodic transitions between the clockwise and counterclockwise circulation patterns (Veronis 1965; Radko 2013). One of the most unstable modes shown in the eigenspectrum results in the clockwise rolls while its counterpart results in the counterclockwise rolls. The periodic transition of overstable motions is evidenced by the oscillation of the transient energy evolution shown in figure 5(b). In this regard, Schmid & Langre (2003) have also observed oscillating transient evolution in the analysis of a generic case of coupled-mode flutter, and Predtechensky et al. (1994) have experimentally observed the periodic behaviour of double-diffusive convection in a Hele-Shaw cell. It is expected that the perturbed salinity field shown in figure 5(c) is distorted unidirectionally at a given instant, which is dissimilar to that in the steady instability regime shown in figure 4(c). Other notable feature distinguishing overstable motions from steady convection are that the initial transient amplification seems to disappear here, which we will discuss in detail in § 4.2.2.

Figure 4. Diffusive convection in steady instability regime: (a) eigenspectrum; (b) transient energy growth; (c) perturbed salinity field. The parameters are Ra = 1.005Rc, R_S = 0.005Rc, k = γ_c, m = 0, Re = 0 and $\varLambda$ = 0. Hereafter, Rc = 1707.763, γ_c = 3.116 and λ = 2${\rm \pi}$/k. In panel (b), two additional cases with Ra = 1.0048Rc and 1.0052Rc are shown for comparison.

Figure 5. As with figure 4 but in the oscillatory instability regime. The parameters are Ra = 1.0209Rc, R_S = Rc, k = γ_c, m = 0.

To give an insight into the intrinsic transition of overstable motions between the clockwise and counterclockwise patterns, we perform a transient energy analysis in a cycle, as shown in figure 6, and plot a sequence of accompanying perturbed salinity profiles, as shown in figure 7. The transient energy analysis is explained in Appendix A. As shown in figure 6, there are four energy components, which are VT (VS) for the energy exchange between the hydrodynamic velocity field and the temperature (salinity) field, VD for the viscous dissipation and Prod for the energy transfer between the perturbed hydrodynamic field and the base flow.

Figure 6. Transient energy growth and transient energy budget versus time. The parameters are as with figure 5. VT (VS) is the energy transfer between the velocity fluctuation field and temperature (salinity) field, Prod is the energy production by the basic shear and VD is the viscous dissipation.

Figure 7. Time evolution of the perturbed salinity fields in a cycle. The corresponding times stamps of the plots are specified in figure 6.

It is noted that in contrast to VS, VT is always positive, indicating that the temperature stratification always releases its potential energy into the hydrodynamic field, regardless of whether the fluid parcel is moving upwards or downwards. This is because the high diffusivity of heat allows the fluid parcel to quickly adjust its temperature to the environment. As with the transient energy G, the transient energy budgets also exhibit periodic oscillation, except that Prod is always zero as we consider a hydrostatic base state here. The dominant energy balance is between the heat diffusion and viscous dissipation. However, it is the energy exchange between the hydrodynamic field and salinity stratification (rather than temperature stratification) that acts in phase with G. This is because the overstable motions benefit from the low diffusivity of salt (Veronis 1965). A cycle of overstable oscillation clearly contains two oscillations of the transient energy G; one corresponds to the clockwise cell while the other to the counterclockwise cell. As shown in figure 6, the first oscillation of G (let us assume it is clockwise) starts with the increase of VS which is initially negative. When VS becomes positive, namely the salinity stratification starts to transform its potential energy into kinetic energy, the transition from clockwise to counterclockwise cells occurs, which is visualized in figure 7(c–g). After this transition, the newly generated counterclockwise cell increases in intensity until the peak (labelled as ‘h’ in figure 6) is reached. As the energy production by salt diffusion is reduced and finally the salinity field begins to absorb energy from the hydrodynamic field, the intensity of the counterclockwise rolls is weakened, thereby leading to another transition from counterclockwise to clockwise cell. The second transition is visualized in figure 7(k–o).

4.2 The effects of phase boundary

Regarding the effects of an ice–water interface on diffusive convection, we first have a look at the stability boundary shown in figure 8. Davis et al. (1984) pointed out that the effect of a steady solid–liquid interface on the buoyancy-driven instability in its underlying liquid layer is determined by the temperature difference between the liquid and solid layers, i.e. $\varLambda$, which is equal to the solid thickness in the equilibrium state. Since the solid–liquid interface becomes a worse thermal conductor as $\varLambda$ increases, the critical Rayleigh number Ra_c and the critical wavenumber γ_c which mark the onset of flow motion are decreasing functions of $\varLambda$ (Davis et al. 1984).

Figure 8. Effects of phase boundary on the convection threshold of diffusive convection: (a–c) low-R_S cases; (d) and (e) high-R_S cases. In panel (b), VT (VS) is the energy transfers between the velocity fluctuation field and temperature (salinity) field, Prod is the energy production by the basic shear and VD is the viscous dissipation. In panel (e), ΔRa_c denotes the difference between Ra_c for conducting boundary and phase boundary. Here, Re = 0 and m = 0.

In double-diffusive cases, the effect of a phase boundary on the convection threshold is more complex, which depends not only on $\varLambda$, but also on R_S. When R_S is relatively low, as shown in figure 8(a and b), (Rac, γ_c) also decreases with larger values of $\varLambda$, which is qualitatively consistent with the results of Hirata et al. (2012) who studied the onset of double-diffusive convection in an under-ice melt pond. The decrease of Ra_c can be explained by a modal energy analysis shown in figure 8(c). The derivation of perturbation energy equation is given in Appendix A. The notation Prod, VD, VT and VS represent the same energy components as the transient energy analysis. We highlight the energy difference between $\varLambda$ > 0 and $\varLambda$ = 0, and to isolate the effects of $\varLambda$, all the parameters except $\varLambda$ are kept constant. It is readily seen that as $\varLambda$ increases, energy production by heat diffusion (VT) is increased. The increasing potential energy released from the unstable temperature stratification is more than that collected and stored by the stable salinity stratification. A part of this difference is dissipated by viscosity, while the remaining is transferred into the perturbation kinetic energy, driving the amplification of overstable oscillations.

In high-R_S cases, however, we find that the presence of an ice–water interface tends to stabilize the system, despite that this effect is too small to make a real difference since Ra_c is very high in this situation. Therefore, it is deduced that in oceanic systems where R_S is relatively high, the existence of an ice–water interface has negligible effects on the linear stability of the subglacial diffusive convection. In other contexts, where R_S is relatively low (for example, the solidification of binary alloy), the existence of a phase boundary can change the heat flux at the interface so that the convection threshold will be decreased in a certain degree.

It is interesting to study the transition between steady and oscillatory instability regimes. In figure 9, we plot the magnified view of figure 8(a) in the low-R_S regime. According to (3.20) (Veronis 1965), the transition from steady convective instability to oscillatory instability occurs at R_S/Ra((β + Pr)/(1 + Pr)) = β, which is simplified to be Ra = 91R_S for β = 0.01 and Pr = 10. In the case of a perfectly conducting boundary, our calculation follows this prediction well, as shown by figure 9(a). However, in the presence of a phase boundary, the transition is modified and the parameter range in which the conductive state is linearly unstable to steady convection is expanded, as shown in figure 9(b).

Figure 9. Magnified view of figure 8: (a) case with a conducting boundary; (b) case with a phase boundary.

The effects of phase boundary on transient energy growth in the oscillatory instability regime are illustrated in figure 10. It is seen that the ice–water interface does not change the oscillatory nature of overstable motions. However, an originally decaying oscillatory mode with a perfectly conducting boundary will become linearly unstable when considering an ice–water interface, which is consistent with the result of the modal analysis in low-R_S cases where the phase boundary can decrease the threshold of diffusive convection. Furthermore, it is found that with the increase of ice thickness, the frequency of overstable oscillation will be decreased, as shown in figure 10(c).

Figure 10. (a,b) Effects of phase boundary on transient energy growth; (c) frequency of overstable oscillation f ₁ as a function of $\varLambda$. Here, k = γ_c and m = 0.

4.3 The effects of shear

Now we turn to examine the effects of Kolmogorov flow on the diffusive instability. To isolate the effects of the shear, we set $\varLambda$ = 0. Before investigating the thermohaline-shear instability identified by Radko (2016) from the perspective of instability nature (oscillatory or steady) and initial transient growth, we first consider cases of R_S < 10⁶ for which three different instabilities are expected to occur for low (~0.1), moderate (~1) and high Re (~10).

4.3.1 Three instability regimes for R_S < 10⁶

When the shear intensity is weak, as shown in figure 11(a), the location of the most unstable modes in the eigenspectrum is analogous to the unsheared overstable motions and, compared with the remaining modes, they travel with the largest absolute phase speed. Moreover, as shown in the inserted perturbed salinity fields in figure 11(a), these two instability modes are triggered at y = 2/5 and 3/5, which is different from 1/4 and 3/4 predicted by Radko (2016). This might be because Radko (2016) considered the instability of transverse rolls with the vertical direction being periodic, while we consider here the instability of oblique rolls with the vertical boundary being no-slip. Nevertheless, it is believed that these two modes pertain to the oscillatory instability instead of the steady convective or Kelvin–Helmholtz instability, which is evidenced by the oscillatory behaviour of transient energy growth G and the accompanying energy budget E shown in figure 11(b and c). However, as shown in figure 11(b), the frequency of overstable oscillation is decreased for larger values of Re, which implies that the shear is prone to inhibiting the overstable oscillation here.

Figure 11. Oscillatory diffusive instability in the presence of Kolmogorov flow: (a) eigenspectra for a growing and a decaying case; (b) transient energy growths versus time; (c) transient energy budget versus time for Re = 0.04. The remaining parameters are Ra = 1.027Rc, R_S = Rc, $\varLambda$ = 0, k = γ_c × (1/3)^0.5, m = γ_c × (2/3)^0.5.

It should be clarified that the oscillatory behaviour identified here is dissimilar to that of Radko (2019). While Radko (2019) also found that the exponential growth of the quadratic perturbation norm is modulated on small time scales, this modulation in fact has an angular frequency twice that of the background shear which has a time-dependent amplitude. Therefore, we argue that the oscillation of Radko (2019) should be attributed to the time-dependent nature of the basic shear, whereas the oscillation shown in our study is inherited from the oscillatory nature of the overstable motions.

The diffusive instability in the presence of Kolmogorov flow of moderate intensity is displayed in figure 12. Unlike the low-Re cases, as shown in figure 12(a), there exists only one most unstable mode, located on the core-mode branch which has zero phase speed. This suggests that the instability is triggered at y = 1/2 where the basic velocity is zero but the shear rate is the largest. Since the Kelvin–Helmholtz instability by the Kolmogorov flow is negligible for Re = 1, the instability shown here is still buoyancy driven. Specifically, here the diffusive instability is steady rather than oscillatory, which is evidenced by the smooth growth of G shown in figure 12(b). The transient energy budget shown in figure 12(c) suggests that in the presence of moderate-Re Kolmogorov flow, the energy productions by the basic shear and salt diffusion are negligible in comparison to the energy balance between the heat diffusion and viscous dissipation. In this situation, the diffusive convection behaves like a single-component flow (e.g. Rayleigh–Bénard convection), despite that the circulation patterns in the present case are affected by the shear. The perturbed salinity field shown in figure 13 indicates that the system develops initially as the multi-cells pattern which will not persist for long. As time progresses, the system evolves to a stable circulation pattern, consisting of large-size cells localized at y = 1/2, and the perturbed energy starts to grows/decays exponentially with time.

Figure 12. As with figure 11 but in the steady diffusive instability regime. The two cases in panel (b) correspond to those in panel (a) with Re = 1. The transient energy budget E_k shown in panel (c) corresponds to the case of Ra = 1.2Rc in panel (a). Other parameters are the same as in figure 11.

Figure 13. Time evolution of salinity field for the case of Ra = 1.2Rc in figure 12.

Once the intensity of the basic shear is sufficiently strong, the fluid system might be susceptible to the dynamic instability induced by the basic shear itself, which refers to Kelvin–Helmholtz instability for Kolmogorov flow. Figure 14 illustrates this situation. Looking at figure 14(a), one finds that, while the original Kelvin–Helmholtz instability modes remain and are still the most unstable, the presence of double diffusion generates a large number of new modes. These new modes lead to the initial fast growth of G, as shown in figure 14(b). This transient amplification is due to the overstable motions generated at the initial instant, which is evidenced by the oscillation of transient energy budget E_k shown in figure 14(c). From figure 14(c), it is seen that the energy production by the basic shear is always positive, which is consistent with the idea that the high-Re Kolmogorov flow is a destabilizing forcing in the present system. In the long term, this forcing is dominant and the kinetic energy transferred from the basic shear into the perturbed hydrodynamic field is mainly compensated by the viscous dissipation. In the short term, however, the potential energy released/stored by the double-diffusive stratification is considerable and has the propensity to govern the initial dynamics. This is visualized in figure 15 where we display the time evolution of the perturbed salinity field. At t = 10, a pair of convective rolls arise in the upper region of the domain while the middle region is nearly uniform, implying that the oscillatory instability is triggered faster than the Kelvin–Helmholtz instability. As time progresses, this pair of rolls become distorted and elongated while simultaneously, near the inflection point of basic velocity profile (i.e. y = 1/2), the Kelvin–Helmholtz vortices set in. Subsequently, the overstable motions decrease in amplitude and gradually disappear, whereas the emerging Kelvin–Helmholtz streaks are intensified and eventually become dominant. Therefore, in the asymptotic stage, one observes in figure 14(b) a smooth growth of G similar to the case without double diffusion and in figure 15(d), the typical Kelvin–Helmholtz pattern.

Figure 14. Kelvin–Helmholtz instability subject to double diffusion: (a) eigenspectra in the presence/absence of double diffusion; (b) accompanying transient energy growths versus time; (c) transient energy budget for the case with double diffusion. The parameters are Re = 50.6, k = 3.262, m = 0. For the case with double diffusion, we have Ra = 1984, R_S = 100Rc.

Figure 15. Time evolution of salinity field for the case with double diffusion in figure 14.

Preceding discussions characterize the properties of three instabilities that can arise in the present system, i.e. the oscillatory, steady and Kelvin–Helmholtz instabilities. To determine the boundaries distinguishing them, we present the stability maps for various shear intensity in figure 16. Regarding the transition from oscillatory to steady convective instability, it is found that there exist two paths of transition in the present system. (i) For the low-Re case, as shown in figure 16(a), the system first bifurcates into overstable motions, which, as with in the nonsheared case (Balmforth et al. 2006), will be superseded by the steady convection once Ra exceeds a critical value. (ii) Likewise, increasing the shear intensity, one will find a critical Re above which the steady convection is the only linear instability. This critical Re is located at the so-called codimension-two point (Predtechensky et al. 1994), which connects the two paths of the transition from oscillatory to steady instability. The codimension-two points of the neutral curves for various R_S constitute a curve which grows with increasing R_S; see figure 16(b). One notable difference between the two paths is that the first path is valid for rolls with arbitrary orientation, whereas the second path is only valid for rolls with a transverse component since the shear has no effect on the longitudinal rolls; see the vertical curve in figure 16(b). In other words, in the presence of shear, longitudinal rolls are always oscillatory as long as (Ra, R_S) are lying in the oscillatory instability regime determined by (3.20b), whereas the transverse and oblique rolls might be either oscillatory or steady, depending on both the buoyancy and the shear intensity.

Figure 16. (a) Boundaries of oscillatory and steady instabilities of transverse rolls for R_S = 100Rc; (b) critical Re at codimension-two point for various R_S, with TRs being the transverse rolls (m = 0), ORs the oblique rolls (we consider k = m) and LRs the longitudinal rolls (k = 0); (c) transition from diffusive instability to Kelvin–Helmholtz (KH) instability for R_S = 100Rc. The vertical dashed line marks the transition of transverse rolls from the diffusive instability to the mixed (diffusive and Kelvin–Helmholtz) instability region.

The transition from buoyancy-driven instability to dynamic instability is illustrated in figure 16(c). Since the modal stability of longitudinal rolls are unaffected by the shear, longitudinal rolls are always unstable to the diffusive instability. For perturbations with a transverse component, the transition from diffusive instability to dynamic instability is observed. As shown in figure 16(c), there exist four distinct regions in the Ra–Re plane which correspond to the linearly stable region, diffusive instability region, Kelvin–Helmholtz instability region, and the mixed diffusive and Kelvin–Helmholtz instability region. In the diffusive instability region, buoyancy is the only destabilizing factor and likewise, high-Re Kolmogorov flow is the only destabilizing factor in the Kelvin–Helmholtz instability region. In the mixed region, both buoyancy and shear will contribute to the destabilization of the system. Since the Kelvin–Helmholtz instability is two-dimensional with m = 0, the transition to Kelvin–Helmholtz instability occurs at a higher Re for oblique rolls in comparison to the transverse rolls.

4.3.2 Transient growth for R_S < 10⁶

We are now able to explain a contradiction arising in § 4.3.1 which is concerned with the initial transient growth of G. As mentioned before, Konopliv et al. (2018) evaluated the initial transient growth of the fingering convection and attributed it to a mechanism relying on the existence of density gradients. While we also detect this rapid growth in steady instability regime in the absence of shear and in the Kelvin–Helmholtz instability regime, it appears to be eliminated in low- and moderate-Re cases shown in figures 11 and 12.

Let us recall the transient growth mechanism of Konopliv et al. (2018), i.e. ‘The fast initial growth occurs when θ* and s* initially have opposite signs, so they act in phase on the density field… As a result, v* grows strongly, which in turn will cause θ* and s* to move in the same direction until they eventually are of the same sign and sufficiently larger for the transient growth to decay.’ This description is highly consistent with our observation in figure 6 where we show the transient energy budget of overstable oscillation, despite that in Konopliv et al. (2018), no stable oscillation in the long term but only an initial fast growth is observed. This indicates that the so-called transient growth mechanism of double-diffusive convection by Konopliv et al. (2018) is in fact the mechanism responsible for the generation of the overstable motions. In other words, in a double-diffusive system, the optimal initial condition to achieve the maximum transient growth must be susceptible to the overstable motions, despite that unlike the diffusive configuration, this periodic flow might not persist too long in fingering convection. Therefore, it is expected that in the presence of shear, the optimal initial condition must be the longitudinal rolls since the shear tends to inhibit the oscillatory instability with a transverse component, as discussed in last section. This is evidenced by figure 17(a) where we show the contours of G as in the k-m plane in moderate-Re case. Now it is easy to deduce that the initial transient growth in figures 11 and 12 being insignificant is because we consider oblique rolls there. Furthermore, a conclusion can be drawn that regardless of the existence of a shear and of the instability type (oscillatory, steady or dynamic), the initial transient growth is exactly the first growth of G in a cycle of overstable oscillation, i.e. the growing stage in figure 6(a–h). In the oscillatory instability regime, particularly, this transient growth as well as the subsequent decay will repeat themselves over and over again, thereby leading to the so-called overstable motions. The effects of (Ra, R_S, Re) on the transient growth is shown in figure 17(b). It is seen that the larger values of Ra or smaller values of R_S result in a larger G_m (G_m = max(G) in a decaying or neutral case) since they correspond to a stronger buoyancy force. However, a more intense shear intensity also yields a larger G_m, which might be due to the lift-up mechanism, despite that this mechanism cannot induce the initial transient growth in pure Kolmogorov flow.

Figure 17. (a) Dependence of G_m on the horizontal wavenumber (k, m) with G_m = max(G), Ra = 2385, R_S = 40Rc and Re = 1; (b) G_m as a function of m; (c and d) G_m/max(G_m) as a function of m where G_m = max(G) for a decaying case. In panel (d), Ra = 1.7Rc and R_S = 80Rc.

It is noted that in many other sheared complex fluid systems, the longitudinal rolls also achieve the maximum transient growth, e.g. Jerome et al. (2012) for Rayleigh–Bénard–Poiseuille/Couette flow and Zhang et al. (2015) for electro-convection with Poiseuille flow. However, transient growth in these systems is subsequently different than our present system since they are due to the Orr mechanism or lift-up mechanism and relies on the existence of high-Re shear. In modal analysis, the longitudinal rolls are unaffected by the background shear, which is valid for the present system as with the aforementioned sheared convective flows. However, the non-modal results indicate that the buoyancy or Coulomb force has non-negligible effects on the initial transient growth by the Orr mechanism or lift-up mechanism. Specifically, in the lift-up mechanism, the optimal wavenumber to achieve the maximum transient growth is independent of shear intensity, but is influenced by the strength of the force that drives the convective motions (Jerome et al. 2012; Zhang et al. 2015). In the transient growth mechanism for a double-diffusive system (let us call it an overstable mechanism), more complex behaviours are expected. Regarding the effects of buoyancy on the optimal wavenumber, as shown in figure 17(c), it is seen that increasing Ra or equivalently decreasing R_S results in a larger optimal wavenumber. As for the effect of shear, it is found that the low- (not shown here) and moderate-Re Kolmogorov flow has negligible effects on the optimal initial condition whereas in the case of the high-Re Kolmogorov flow, a larger optimal wavenumber is required to achieve to the maximum transient growth; see figure 17(d).

4.3.3 Thermohaline-shear instability

In this section, we consider a more realistic situation and investigate the effects of shear on diffusive instability in the context of oceanography where Ra and R_S are much higher than those considered in §§ 4.3.1 and 4.3.2. The typical parameters of diffusive convection in high-latitude oceans are 1.5 < R_ρ < 10 and 1 < Ri < 10, Pe ~ 10³, where R_ρ is the background density ratio, Ri the Richardson number and Pe the Péclet number (Radko 2016). After a simple transformation from the term of R_ρ–Ri–Pe formulation to the term of Ra–R_S–Re formulation, we deduce that 10⁵ < Ra < 10⁶, 1.5 × 10⁷ < R_S < 10⁹, Re ~ 100. Here we considered R_S ~ 10⁸, while Re and Ra are independent and dependent variables, respectively. Our aim is to identify whether the thermohaline-shear instability found by Radko (2016) is oscillatory or steady and to discuss whether the transition between different instability regimes is realistically expected in the ocean.

In figure 18, we show the stability diagram for the sheared diffusive convection as a function of Re for R_S = 10⁸, τ = 0.01 and Pr = 10. In contrast to figure 16, here the critical Ra is a decreasing function of Re. This implies that at each point of the stability boundary within 0 < Re < 140.4, the system is individually stable with respect to diffusive instability and to Kelvin–Helmholtz instability, but unstable to the superposed instability. This superposed instability is precisely the thermohaline-shear instability identified by Radko (2016), which favours the transverse disturbances. The non-modal analysis suggests that the thermohaline-shear instability is oscillatory for Re lower than 81, above which the instability is steady. The transition from oscillatory instability to steady instability occurs at a codimension-two point, which has been specified in figure 18(a). After collecting all the codimension-two points in a wide range of R_S, we obtain the critical Re and critical Ra marking the transition from oscillatory to steady instability, as shown in figure 18(b). It is seen that the critical Re and critical Ra are the increasing functions of R_S. Moreover, further increasing Re, the critical Ra will become negative after Re exceeds 140.4, marking the transition to the Kelvin–Helmholtz instability regime.

Figure 18. Transition from oscillatory to steady instability in high-R_S case: (a) critical Ra as a function of Re for R_S = 10⁸; (b) critical Re and critical Ra for codimension-two point as a function of R_S; (c) asymptotic energy budget along the stability boundary in panel (a).

Here, we perform an asymptotic energy analysis for thermohaline-shear instability. Different from Radko (2016), in which only Pe is varied and the remaining parameters are all fixed, we proceed along the linear instability boundary, along which both Ra and Re are varied. It is seen from figure 18(c) that there exist two different regimes in which the energy budgets are entirely different, and the transition between these two regimes coincides with the codimension-two point that distinguishes the oscillatory and steady instabilities. In the oscillatory instability regime, the ultimate energy resource to amplify disturbances is the heat diffusion. The energy production by the shear is negligible or even negative, but the presence of the background shear promotes the release of energy from the temperature gradient into the perturbed hydrodynamic field. This is consistent with the hypothesis of Radko (2016) that rationalizes the catalytic effects of shear on the transfer of potential energy into kinetic energy. In the steady instability regime, however, the basic shear directly produces energy to the perturbed hydrodynamic field as with in many dynamic instabilities such as Kelvin–Helmholtz instability or Tollmien–Schlichting instability. In this sense, different from its oscillatory counterpart which involves complex interaction between shear and double diffusion (i.e. the catalytic effect), the steady thermohaline-shear instability resembles more a simple superposition of instabilities due to buoyancy and shear. Clearly, the condition for the steady thermohaline-shear instability to occur is more restrictive than the oscillatory one, since it requires a higher intensity and special velocity profile (containing inflection points) of the background shear. Therefore, we argue that the generation of thermohaline staircases in the Arctic Ocean is more likely due to the mechanism inducing the oscillatory instead of steady thermohaline-shear instability. This argument is consistent with the result of Radko (2019) that the thermohaline-shear instability can be triggered by the interaction between double diffusion and an unbounded, time-dependent Couette flow.

Since the codimension-two point occurring in figure 18 is within the parameter ranges of the real oceanic systems and the sheared transverse rolls in this situation are the least stable, the transition from oscillatory to steady instability discussed here is realistically possible. However, the transition into Kelvin–Helmholtz instability regime is realistically impossible, because it requires a sufficient shear intensity which is not satisfied in oceanic conditions.

It is necessary to investigate when the transient growth by double diffusion becomes significant in destabilizing the present system. To this end, we examine the effects of R_S and Re on G_m, considering transverse rolls. Different from the cases in figure 17, it is obvious that the optimal initial condition for thermohaline-shear instability favours the transverse rolls. As shown in figure 19(a), with the increase of R_S, G_m first undergoes nonmonotonic variation and eventually approaches a fixed value that is larger than the low-R_S case but is still of order 10. This is due to the fact that the oscillatory convection is the low-amplitude motion, independent of the strength of the stabilizing force which refers to R_S here. When considering the presence of a base flow, it is found that G_m is increased by two orders of amplitude and a sudden jump of G_m at the codimension-two point is observed. Specifically, in the oscillatory instability regime where G_m is an increasing function of Re, the transient energy returns to the low-amplitude oscillation after the initial fast growth, see the case of Re = 40 shown in figure 19(b). In the steady instability regime where G_m decreases with larger values of Re, however, the transient energy of a marginally stable case in the long term maintains a level identical to G_m. Therefore, it is deduced that the initial transient growth has important influence on the oscillatory thermohaline-shear instability. The reasons are twofold: (i) from the perspective of modal stability, transverse rolls are the least stable in thermohaline-shear instability and (ii) from the perspective of non-modal stability, the presence of shear makes the amplitude of initial transient growth an order higher of the oscillatory convection. In contrast, the transient growth by double-diffusion has negligible effect on the steady instability which, in the marginally stable case, is of the amplitude identical to the transient growth. In addition to the energy resource of instabilities, this is another principle difference between the steady thermohaline-shear instability and the oscillatory one.

Figure 19. Effects of R_S and Re on the transient growth by double-diffusion. In panel (a), the results for Re = 0 are obtained by setting Re = 1 but considering a hydrostatic base flow. All the cases in panels (a) and (b) are marginally stable, with the corresponding neutral stability curves shown in figures 8(d) and 18(a).

5. Conclusion

In this paper, we have studied the effects of a solid–liquid interface and a Kolmogorov flow on the diffusive instability, using modal and non-modal stability analyses as well as energy analysis. In the absence of the phase boundary and shear, we showed that, with the overstable motions transitioning from the clockwise to counterclockwise circulation patterns, the transient energy undergoes periodic oscillation, driven by the salinity field which alternates between releasing and storing potential energy. For steady convective instability, however, a smooth evolution of transient energy is observed.

Regarding the effects of the phase boundary, the change of critical Rayleigh number Ra_c depends on the thickness of the solid phase and the salinity Rayleigh number R_S. Specifically, for relatively small R_S (R_S < 10⁶), Ra_c is decreased with larger values of solid thickness, due to the increasing energy transfer from the temperature gradient to the perturbed hydrodynamic field. In the context of ice melting in polar regions where R_S is relatively large, the presence of an ice–water interface tends to stabilize the diffusive instability. However, in the absence of basic shear, since the buoyancy force required to drive the convective flow for large R_S is very strong, this stabilizing effect is in fact trivial. However, while the existence of a phase boundary does not change the oscillatory nature of overstable motions, it makes the underlying fluid more susceptible to steady convection instead of the oscillatory convection.

In the presence of shear, three instabilities could be dominant for relatively small R_S at different values of Reynolds number Re, which are the oscillatory instability for Re ~ 0.1, steady convective instability for Re ~ 1 and Kelvin–Helmholtz instability for Re > 50. The properties of these instabilities have been studied from the perspective of transient energy evolution. Specifically, in the Kelvin–Helmholtz instability regime, it is found that the dynamic instability develops more slowly than the oscillatory instability. Therefore, in the early phase, the overstable motions will arise in the low-shear region, and yields an initial rapid growth and the short-term periodic oscillation of the transient energy. As time progresses, the buoyancy-driven patterns will be elongated and gradually disappear, while at the same time, Kelvin–Helmholtz streaks will be generated in the inflection point of the basic velocity profile and eventually the dynamic instability becomes dominant.

For relatively large R_S, we examined the thermohaline-shear instability identified by Radko (2016). By inspecting the temporal evolution of disturbances, we found that the thermohaline-shear instability could transition from the oscillatory instability regime to the steady one, and this transition is realistically possible in the context of oceanography. The energy analysis along the neutral stability curve of thermohaline-shear instability suggested that in the oscillatory instability regime, the heat diffusion and viscous dissipation are dominant in energy balance, with energy production by shear being negligible. This result confirmed that the hypothesis proposed by Radko (2016) to explain the catalytic role of shear is reasonable when the thermohaline-shear instability is oscillatory. In the steady instability regime, however, the mean shear makes direct and considerable contribution to the energy to intensify instabilities. Therefore, the steady thermohaline-shear instability is more like a direct superposition of instabilities by shear and buoyancy. We have also investigated the mechanism of the initial transient growth driven by double-diffusion (Konopliv et al. 2018). Results imply that the optimal initial condition to achieve the maximum transient growth must be susceptible to the overstable motions. For relatively small R_S, since the shear inhibits all disturbances with a transverse component, the optimal initial condition favours the longitudinal rolls. For relatively large R_S when the thermohaline-shear instability is induced, the optimal initial condition favours the longitudinal rolls, and the initial transient growth, driven by double diffusion and enhanced by the basic shear, might have significant impact on oscillatory thermohaline-shear instability.

In the future, it is necessary to numerically/experimentally investigate whether a thermohaline-shear system that is stable to the oscillatory instability predicted by modal analysis could be destabilized by the initial transient growth. We intend to investigate the subcritical bifurcation of double-diffusive convection subject to shear by weakly nonlinear instability properties and the transition from convective instability to absolute instability of the present system.