2.1. A parametrization of turbulent diapycnal diffusivities

Calibration of an accurate parametrization for mixing in a stratified turbulent flow has always been recognized as a fundamental problem in theoretical fluid dynamics (see the recent reviews of Peltier & Caulfield (Reference Peltier and Caulfield2003), Gregg et al. (Reference Gregg, D'Asaro, Riley and Kunze2018) and Caulfield (Reference Caulfield2021)). Owing to the accumulating data from numerical simulations, laboratory experiments and field measurements, new parametrizations underpinned by improved physical understanding have emerged in recent decades (see e.g. Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Bouffard & Boegman Reference Bouffard and Boegman2013; Salehipour et al. Reference Salehipour, Peltier, Whalen and MacKinnon2016). These revised parametrizations are possible replacements for that of Osborn (Reference Osborn1980), which was based upon an analysis of the turbulent kinetic energy equation and which suggested that the mixing efficiency should be fixed to a value near  0.2.

Among these different schemes for the parametrization of turbulent diffusivity, that of Bouffard & Boegman (Reference Bouffard and Boegman2013) will be of particular interest in the present context. It is apparently unique in that it stresses the importance of the Prandtl and Schmidt numbers as well as a measure of turbulence intensity by the buoyancy Reynolds number. It contains no explicit dependence upon shear as in that of Salehipour et al. (Reference Salehipour, Peltier, Whalen and MacKinnon2016).

Specifically, the Bouffard & Boegman (Reference Bouffard and Boegman2013) parametrization for diapycnal diffusivity is based upon the following functional fits to the datasets employed to constrain it. For the salinity and temperature contributions to the diapycnal diffusivities $K_\rho$ as a function of $Re_b$ and $Pr$, distinctly different fits are provided. The data on which these functional fits were based include those from both laboratory experiments (e.g. Jackson & Rehmann Reference Jackson and Rehmann2003; Rehmann & Koseff Reference Rehmann and Koseff2004) and numerical simulations (e.g. Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Smyth, Nash & Moum Reference Smyth, Nash and Moum2005). Their parametrization is as follows:

(2.1a)$$\begin{gather} K^{BB}_\rho(Re_b, Pr)=\kappa, \quad \text{for} \ Re_b<10^{{2}/{3}}Pr^{-{1}/{2}}, \end{gather}$$

(2.1b)$$\begin{gather}K^{BB}_\rho(Re_b, Pr)=\frac{0.1}{Pr^{{1}/{4}}}\nu Re_b^{{3}/{2}}, \quad \text{for} \ 10^{{2}/{3}}Pr^{-{1}/{2}}< Re_b<(3 \ln \sqrt {Pr})^2, \end{gather}$$

(2.1c)$$\begin{gather}K^{BB}_\rho(Re_b, Pr)= 0.2 \nu Re_b, \quad \text{for} \ (3 \ln \sqrt {Pr})^2< Re_b<100, \end{gather}$$

(2.1d)$$\begin{gather}K^{BB}_\rho(Re_b, Pr)= 2 \nu Re_b^{{1}/{2}}, \quad \text{for} \ Re_b>100. \end{gather}$$

The regions of different dependence upon $Re_b$ from (2.1a) to (2.1d) are described as the ‘molecular-diffusion-driven regime’, the ‘buoyancy-controlled regime’, ‘the transitional regime’ and, at highest buoyancy Reynolds number, the ‘energetic regime’. In figure 1 we have plotted the dependence of $K_\rho$ as a function of $Re_b$ for the Prandtl number for temperature, $Pr=7$, and the Prandtl number (or Schmidt number) for salinity, $Pr=700$. It is clear on the basis of figure 1 that the difference in $Pr$ is associated with significantly different dependences of diapycnal diffusivities as a function of $Re_b$ when the buoyancy Reynolds number is moderate ($Re_b<100$). In this range, the low diffusivity renders the salinity much more difficult to mix than temperature. There exists a slight discontinuity between the buoyancy-controlled regime and the transitional regime in (2.1) at $Pr=7$ but it does not have any significant effect in our following analysis. It is also important to note that the transitional regime as described by the Osborn (Reference Osborn1980) relationship, assumed as basis for the parametrization, does not exist in the parametrization for $Pr=700$.

Figure 1. Prediction of $K_s(Re_b)$ and $K_\theta (Re_b)$ based the parametrization of Bouffard & Boegman (Reference Bouffard and Boegman2013).

2.2. A linear stability criterion for the staircase onset

Given the parametrization of diapycnal diffusivity for both heat and salt, we are in a position to analyse their implications for the stability of a water column that is stratified in both temperature and salinity. We will perform this assessment in the context of the similar mean-field model as that recently employed to investigate layer formation in the salt-fingering regime (Ma & Peltier Reference Ma and Peltier2021a, hereafter MP). The governing mean-field equations for the large-scale temperature field $\varTheta (z)$ and salinity field $S(z)$ are as follows:

(2.2a)$$\begin{gather} \frac{\partial \varTheta }{\partial t}={-}\frac{\partial}{\partial z} F_\varTheta= \frac{\partial}{\partial z} \left( K_\varTheta(Re_b) \frac{\partial \varTheta}{\partial z}\right), \end{gather}$$

(2.2b)$$\begin{gather}\frac{\partial S }{\partial t}={-}\frac{\partial}{\partial z} F_S= \frac{\partial}{\partial z} \left(K_S(Re_b) \frac{\partial S}{\partial z}\right). \end{gather}$$

Equations (2.2) are coupled because both diffusivities are dependent on $Re_b=\varepsilon /(\nu N^2)$. This control parameter of the model will be perturbed by instability growth through its dependence on $N^2$, which will be modified by the variations in $\varTheta$ and $S$. In our model, $\varTheta$ and $S$ are defined in density units so that the equation of state can be written as $\rho =\rho _0+S-\varTheta$.

In the above equations, we have assumed that $K_\varTheta$ and $K_S$ are independent of the gradient Richardson number $Ri_g$. The reasonableness of this assumption is based upon the fact that the Richardson number is usually large in the Arctic and therefore its effect on the turbulent level is relatively minor compared with the low- to mid-latitude ocean. We will furthermore assume $K_\varTheta$ and $K_S$ to be independent of the density ratio $R_\rho$. In work that will be published elsewhere, we will demonstrate this independence explicitly using direct numerical simulations of Kelvin–Helmholtz billows in the diffusive-convection environment.

It is also important to understand that a multicomponent fluid may still behave, insofar as the characteristics of its turbulence are concerned, as if it were a single-component fluid. The condition under which this is possible is that the turbulent diffusivities of the two components are the same. In this case the advection–diffusion equations for the two species combine into a single such equation for density. As shown by Smyth et al. (Reference Smyth, Nash and Moum2005) and Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) as well as our figure 1, this occurs if and only if the intensity of the turbulence is sufficiently high as measured by the buoyancy Reynolds number to be higher than 100. It is therefore clear on this basis alone that a theory based upon a layer formation criterion for a single-component fluid should not apply to the formation of Arctic staircases since the experimental measurements in field experiments (e.g. Chanona, Waterman & Gratton Reference Chanona, Waterman and Gratton2018; Dosser et al. Reference Dosser, Chanona, Waterman, Shibley and Timmermans2021) have established that the buoyancy Reynolds numbers in the Arctic Ocean are usually smaller than 100, with the exception of some hot spots near the continental shelf.

In order to investigate the conditions that might support the existence of a layering mode of instability, we assume furthermore that the background stratification is characterized by uniform gradients and the one-dimensional fields can be decomposed into a combination of background fields $\bar {\varTheta } =-\varTheta _{z0}z$ and $\bar {S}= -S_{z0}z$ and weak perturbations $\varTheta '$ and $S'$ as

(2.3a)$$\begin{gather} \varTheta(z)=\bar{\varTheta}(z)+\varTheta'(z), \end{gather}$$

(2.3b)$$\begin{gather}S(z)=\bar{S}(z)+S'(z). \end{gather}$$

Assuming an initial state characterized by a constant viscous dissipation $\varepsilon _0$, the buoyancy Reynolds number $Re_b$ will be perturbed from its background value $\overline {Re}_b$ by the amount $Re_b'$, which is given by

(2.4)\begin{equation} Re_b'=\frac{\partial Re_b}{\partial \rho_z}\frac{\partial \rho'}{\partial z}=\frac{\rho_0}{\nu g} \frac{\varepsilon_0}{ \left(\dfrac{\partial \bar{\rho}}{\partial z}\right)^2} \frac{\partial \rho' }{\partial z} ={-} \overline{Re}_b \frac{\dfrac{\partial S' }{\partial z} -\dfrac{\partial \varTheta' }{\partial z}}{\dfrac{\partial\bar{\rho}} {\partial z}}, \end{equation}

which follows by employing a Taylor expansion to determine the perturbation to the squared buoyancy frequency that appears in the definition of the buoyancy Reynolds number. We then expand $K_\varTheta$ and $K_S$ in (2.2) in $Re_b$ keeping only the first-order terms to obtain

(2.5a)$$\begin{gather} \frac{\partial \varTheta'}{\partial t}=\left.-\frac{\partial K_\theta} {\partial Re_b} \right|_{\overline{Re}_b} \overline{Re}_b \frac{1}{R_\rho-1} \left(\frac{\partial^2 S'}{\partial z^2} - \frac{\partial^2 \varTheta'}{\partial z^2}\right)+K_\theta |_{\overline{Re}_b} \frac{\partial^2 \varTheta'}{\partial z^2}, \end{gather}$$

(2.5b)$$\begin{gather}\frac{\partial S'}{\partial t}=\left.-\frac{\partial K_s}{\partial Re_b} \right|_{\overline{Re}_b} \overline{Re}_b \frac{R_\rho}{R_\rho-1} \left(\frac{\partial^2 S'}{\partial z^2} - \frac{\partial^2 \varTheta'}{\partial z^2}\right) +K_s |_{\overline{Re}_b} \frac{\partial^2 S'}{\partial z^2}. \end{gather}$$

In order to analyse the stability of this linear system, we expand the perturbation terms in normal modes as $(\varTheta ', S')=(\hat {\varTheta },\hat {S})\exp (\lambda t)\exp (\textrm {i} k z)$. Upon substitution into (2.5), the system is reduced to an eigenvalue problem for a $2 \times 2$ matrix whose eigenvalues $\lambda$ are roots of the quadratic equation

(2.6)\begin{align} &\lambda^2+k^2\left(K_\theta+K_s+\frac{\partial K_s}{\partial Re_b} Re_b \frac{R_\rho}{R_\rho-1}- \frac{\partial K_\theta}{\partial Re_b} Re_b \frac{1}{R_\rho-1}\right)\lambda \nonumber\\ &\quad+\,k^4\left(K_\theta K_s+\frac{\partial K_\theta}{\partial Re_b} K_s Re_b \frac{1}{R_\rho-1} - \frac{\partial K_s}{\partial Re_b} K_\theta Re_b \frac{R_\rho}{R_\rho -1}\right) =0. \end{align}

A sufficient condition for an eigenvalue $\lambda$ with positive real part is that the third term in (2.6) be negative definite, namely,

(2.7)\begin{equation} K_\theta K_s+\frac{\partial K_\theta}{\partial Re_b} K_s Re_b \frac{1}{R_\rho-1} - \frac{\partial K_s}{\partial Re_b} K_\theta Re_b \frac{R_\rho}{R_\rho -1}<0. \end{equation}

By further writing the diapycnal diffusivity in the form of flux coefficient as in $K_\theta =\nu \varGamma _\theta Re_b$ and $K_s=\nu \varGamma _s Re_b$, the above condition for instability is further simplified to

(2.8)\begin{equation} \frac{\partial \varGamma_s}{\partial Re_b}\frac{1}{\varGamma_s} >\frac{\partial \varGamma_\theta}{\partial Re_b} \frac{1}{\varGamma_\theta} \frac{1}{R_\rho}. \end{equation}

By assuming that $K_S$ and $K_\theta$ have different power-law dependences on $Re_b$ as $K_S\sim Re_b^{\beta _s}$ and $K_\varTheta \sim Re_b^{\beta _\theta }$, the criterion of (2.8) will become

(2.9)\begin{equation} \beta_s-1>\frac{\beta_\theta-1}{R_\rho}. \end{equation}

By comparing (2.9) with the power-law dependences of Bouffard & Boegman (Reference Bouffard and Boegman2013) shown in (2.1), we will immediately arrive at the conclusion that the system will be unstable as long as the value of $Re_b$ remains in the buoyancy-controlled regime of salinity, namely,

(2.10)\begin{equation} 0.18< Re_b<97. \end{equation}

Although the instability criterion does not depend explicitly upon $R_\rho$, the non-dimensional growth rate (evaluated by numerically solving the quadratic equation for the linear instability) is dependent on both $Re_b$ and $R_\rho$, as shown in figure 2. The growth rate is mainly controlled by $Re_b$ and reaches its maximum value at $Re_b$ sufficiently close to but smaller than the upper stability boundary at $Re_b=97$.

Figure 2. Filled contour plot of growth rate (real part of the non-dimensional eigenvalue $\lambda /(k^2 \nu )$) as a function of $Re_b$ and $R_\rho$ evaluated by solving the quadratic equation (2.6) with the implementation of Bouffard & Boegman 's (Reference Bouffard and Boegman2013) parametrization in (2.1).

The physical interpretation of our linear stability analysis can be understood by taking the limit $R_\rho \to \infty$ in the above formulae. In this limit, the contribution from the temperature field vanishes and the salinity is equivalent to density, reducing (2.5) to the form

(2.11)\begin{align} \frac{\partial \rho'}{\partial t}&=\left.\left(\left.-\frac{\partial K_\rho}{\partial Re_b} \right|_{\overline{Re}_b} \overline{Re}_b +K_\rho\right)\right|_{\overline{Re}_b} \frac{\partial^2 \rho'}{\partial z^2} \nonumber\\ &\equiv K_\rho^{eff}\frac{\partial^2 \rho'}{\partial z^2}. \end{align}

In the above limiting form, our instability criterion is equivalent to the condition that the diffusion coefficient is $K_\rho ^{eff}<0$, initially proposed in Phillips (Reference Phillips1972) and extended in the recent work of TZ. The condition that the diffusion coefficient in this diffusion equation is negative is just, as expected, the result for a fluid in which density is determined by only a single diffusing species. When the possible influence of shear on the condition for the formation of layers is eliminated from equation (2.14) of TZ, the same result is obtained. However, TZ's criterion is not applicable to the formation of Arctic staircases, as $R_\rho$ is everywhere finite, and the typical value is 2–7 (Timmermans et al. Reference Timmermans, Toole, Krishfield and Winsor2008). Therefore, our more general analysis based on a two-component system described in (2.5) has to be applied in the Arctic Ocean.

Our criterion (2.10) basically states two important things concerning the staircase formation process in the polar ocean. First, whether or not a staircase forms does not depend on the density ratio $R_\rho$, as it does in the salt-fingering staircase formation process. This is consistent with the analysis of a large dataset documenting staircase measurements in the Canada Basin of the Arctic Ocean in Shibley et al. (Reference Shibley, Timmermans, Carpenter and Toole2017), where the presence of staircases is observed for a wide range of values of $R_\rho$ measured to range from 2 to 9 (demonstrated in figures 5d and 7a of Shibley et al. Reference Shibley, Timmermans, Carpenter and Toole2017). This is distinctly different from the salt-fingering staircases, where no staircases have been found in places where the salt-fingering density ratio is larger than 2 (Radko Reference Radko2013). Second, the criterion (2.10) implies that the staircases only form in the range of intermediate turbulence strength as determined by the buoyancy Reynolds number. This is also consistent with the field measurements, from which accumulating evidence has established that strong turbulence may prevent staircases from forming at all (Shaw & Stanton Reference Shaw and Stanton2014; Guthrie, Fer & Morison Reference Guthrie, Fer and Morison2017; Shibley & Timmermans Reference Shibley and Timmermans2019). We will return to a detailed discussion of these measurements in § 4.

4.1. Spatial variability of staircases

The most abundant oceanographic measurements in the Arctic are provided using high-resolution Ice-Tethered Profilers (ITPs; Krishfield et al. Reference Krishfield, Toole, Proshutinsky and Timmermans2008; Toole et al. Reference Toole, Krishfield, Timmermans and Proshutinsky2011) that have been employed to sample over the entire central Eurasian and Canada Basins. The recent work of Shibley et al. (Reference Shibley, Timmermans, Carpenter and Toole2017) provided detailed analysis of the ITP measurements over the years 2004–2013 to document the spatial variability of the presence of staircases in the entire Arctic region. As shown in their figure 7(a) most of the staircases are observed in the Canada Basin (staircases are present in 80 % of profiles) and they are generally not observed in the Eurasian Basin (staircases are present in 24 % of profiles) (Shibley et al. Reference Shibley, Timmermans, Carpenter and Toole2017).

We turn next to a discussion of the observations of mixing as reported in the most recent work of Dosser et al. (Reference Dosser, Chanona, Waterman, Shibley and Timmermans2021), which uses the fine-scale parametrization of Polzin, Toole & Schmitt (Reference Polzin, Toole and Schmitt1995) to estimate the dissipation rate by analysing the ITP data that sampled the central Arctic Ocean over the years 2002–2019. They compiled and evaluated the distribution of buoyancy Reynolds number $Re_b$ for the Canadian Basin in summer, which is shown in their figure 4(a). Although the turbulence in the Arctic Ocean is an inherently patchy and intermittent process, as shown in the log-normal distribution of $\varepsilon$ in their figure 7(a), the sampled $Re_b$ lies almost entirely within the range of our criterion $0.18< Re_b<97$. Thus the layering mode is expected to keep growing in time although $Re_b$ may continue to vary, until the staircase forms in the Canada Basin. For the Eurasian Basin, however, Dosser et al. (Reference Dosser, Chanona, Waterman, Shibley and Timmermans2021) have shown that, while the geometric mean values of $\varepsilon$ are almost identical to those in the Canada Basin and Eurasian Basin, monthly averaged $N^2$ is approximately three times smaller in the latter (see figure S2 of the supporting information of Dosser et al. (Reference Dosser, Chanona, Waterman, Shibley and Timmermans2021)), resulting in a buoyancy Reynolds number whose mean value is three times larger than the values in the Canada Basin. This means that a larger portion of the tail of the log-normal distribution of $Re_b$ will be high enough to be stable to our criterion, in which case the staircase formation process would be inoperative with $Re_b>97$. This qualitatively explains why most of the staircases are found in the Canada Basin.

4.2. Temporal variability of staircases

The best example of the temporal variability of diffusive-convection staircases is provided by the recent discussion of Guthrie et al. (Reference Guthrie, Fer and Morison2017). These authors found that the Amundsen Basin (of the Eurasian Basin) of the Arctic Ocean was characterized by well-defined staircase structures in 2013 but they were largely absent in the following year of 2014 (also absent in the previous measurements in 2007 and 2008). Through their analysis, they found that the years in which staircases are absence are characterized by consistently higher mixing activity compared with the year 2013 (approximately an order of magnitude higher), although the density ratios $R_\rho$ are almost the same for these years.

This is an entirely expected result based upon the theory advocated herein. The turbulence level in the Eurasian Basin is too high to allow the staircase-forming instability to develop. Our mechanism requires a sufficiently quiescent level of turbulence for the staircase-forming high-Prandtl-number mechanism to function.