3.1 Change in potential energy density ${\rm\Delta}\mathscr{E}_{p}$

In an inviscid flow with positive temperature gradients in both the upward vertical and horizontal directions, a particle moving from a colder higher region to a hotter lower region (while also traversing a horizontal distance ${\rm\Delta}L$ ) experiences a change in potential energy density which under the Boussinesq approximation (Lappa Reference Lappa2012) is

where ${\it\xi}$ is the angle from the horizontal of the path taken by the particle, and ${\it\delta}=\tan ^{-1}(-\partial _{r}{\it\theta}/\partial _{z}{\it\theta})$ is the inclination of the lines of constant temperature to the horizontal. When $0<{\it\xi}<{\it\delta}$ , ${\rm\Delta}\mathscr{E}_{p}<0$ and potential energy is released, which may amplify a perturbation by converting potential energy to kinetic energy, leading to baroclinic instability. Differentiating with respect to ${\it\xi}$ demonstrates that the maximum energy transfer occurs when ${\it\xi}={\it\delta}/2$ . This optimal orientation relates the horizontal displacement to a fixed particle path length through ${\rm\Delta}L=L_{p}\cos ({\it\delta}/2)$ , and hence equation (3.8) can be expressed for change in potential energy normalised by path length squared as

where ${\rm\Delta}\mathscr{E}_{p}/{L_{p}}^{2}$ is normalised by ${\it\rho}_{0}{\it\nu}{\it\kappa}/R^{4}$ . Equation (3.9) will be plotted on axisymmetric base flow fields to identify those regions where the necessary condition for baroclinic instability, ${\rm\Delta}\mathscr{E}_{p}<0$ , is satisfied. Interpretation of modes will be aided by considering the alignment of disturbance mode structures with regions of negative ${\rm\Delta}\mathscr{E}_{p}$ , while instability structures confined to regions of ${\rm\Delta}\mathscr{E}_{p}>0$ would suggest a non-baroclinic origin. As a necessary but not sufficient criterion for baroclinic instability, (3.9) may yield ${\rm\Delta}\mathscr{E}_{p}<0$ in flows where baroclinic instability is not observed, and its inviscid and local nature may lead to disparity against the global stability solutions for the viscous flows obtained within this study. A powerful alternative for assessing the capacity for a flow to be baroclinically active is the available potential energy (Winters & Barkan Reference Winters and Barkan2013), which is described in the subsection to follow.

3.2 Available potential energy $E_{a}$ and density $\mathscr{E}_{a}$

The development of baroclinic instability is contingent on potential energy in the fluid being available for conversion to kinetic energy. The elaboration to follow draws from Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995), Winters & Barkan (Reference Winters and Barkan2013) and references therein, and in this subsection quantities are expressed in their dimensional form. To begin, at any instant in time a body of fluid contains a total potential energy

where $\mathscr{V}$ is the fluid volume. If we imagine permitting the parcels of fluid to adiabatically (i.e. without heat or molecular mass transfer Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995) settle to their equilibrium heights, then a minimum potential energy state is achieved, known as the background potential energy $E_{b}\equiv g\int {\it\rho}z_{\ast }\,\text{d}\mathscr{V}$ . Here $z_{\ast }(\boldsymbol{x},t)$ is the equilibrium height of fluid parcel with density ${\it\rho}$ , and the one-dimensional function ${\it\rho}(z_{\ast })$ is independent of the spatial distribution of fluid parcels.

It follows then that the available potential energy is the difference between potential and background potential energies, i.e.

Winters & Barkan (Reference Winters and Barkan2013) made the further contribution of defining the available potential energy density $\mathscr{E}_{a}(\boldsymbol{x},t)$ , a positive-definite field that integrates to the total available potential energy, i.e.

The available potential energy density is zero where fluid parcels are at their equilibrium height, and larger values indicate greater displacement from equilibrium. It is formally defined as

For convenience (3.13)–(3.14) are recast in terms of temperature consistent with (3.1)–(3.3), yielding

When presented in this study, $\mathscr{E}_{a}$ is given in dimensionless form, scaled by ${\it\rho}_{0}{\it\nu}{\it\kappa}/gR^{2}$ .

Winters & Young (Reference Winters and Young2009), Winters & Barkan (Reference Winters and Barkan2013) describe horizontal convection as being driven by the conversion of available potential energy to kinetic energy, emphasising the importance of the $\mathscr{E}_{a}$ distribution and buoyancy fluxes in these flows. Buoyancy relates to temperature through a linear equation of state approximation, $b={\it\alpha}g{\it\theta}$ .

3.3 Linear stability analysis

The potential for non-axisymmetric three-dimensional instability developing and significantly altering its characteristics on an axisymmetric base flow motivates an application of linear stability analysis. In this analysis, first an axisymmetric solution to equations (3.1)–(3.3) is obtained by fixing $\partial /\partial {\it\phi}=0$ . Then the linearised governing equations are obtained by decomposing the velocity, pressure and temperature into the sum of an axisymmetric field ( $\bar{\boldsymbol{u}}$ , $\bar{p}$ , $\bar{{\it\theta}}$ ) and a small non-axisymmetric perturbation ( $\boldsymbol{u}^{\prime }$ , $p^{\prime }$ , ${\it\theta}^{\prime }$ ). The perturbation field is constructed as a single complex Fourier mode of an azimuthal expansion of the flow field and the wavenumber of the perturbation is a parameter in the stability analysis. Substituting these into (3.1)–(3.3) and retaining only terms of order of the perturbation field yields the linearised Navier–Stokesand thermodynamic equations

where the linear advection terms are $\boldsymbol{L}_{\boldsymbol{u}}^{\prime }=-[(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}^{\prime }+(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\bar{\boldsymbol{u}}]$ and $L_{{\it\theta}}^{\prime }=-[(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}){\it\theta}^{\prime }+(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\bar{{\it\theta}}]$ .

An operator $\mathscr{A}(T)$ is defined as being equivalent to time integration of a perturbation field by equations (3.16–3.18) over some time interval $T$ . Linear stability of a base flow of a given $(\mathit{Ra},Q,\mathit{Pr})$ and azimuthal wavenumber ${\it\beta}$ is dictated by the leading eigenmode of $\mathscr{A}(T)$ , which is obtained using an implicitly restarted Arnoldi method in conjunction with time integration of the base flow and one or more perturbation fields (Cogan, Ryan & Sheard Reference Cogan, Ryan and Sheard2011; Vo et al. Reference Vo, Montabone and Sheard2015). The leading complex eigenvalue $|{\it\mu}|$ is related to the exponential growth rate ${\it\sigma}$ of the corresponding eigenmode (corresponding to the most asymptotically unstable perturbation field) via ${\it\sigma}=\log |{\it\mu}|/T$ . In the case of steady-state base flows (all cases considered in this study reach a steady-state axisymmetric equilibrium), $T$ may be chosen arbitrarily. If $|{\it\mu}|>1$ for any ${\it\beta}$ , then ${\it\sigma}>0$ and the flow is unstable, whereas if $|{\it\mu}|<1$ for all ${\it\beta}$ , then the axisymmetric flow is asymptotically stable to all infinitesimal disturbances. The azimuthal wavenumber relates to the azimuthal wavelength of the instability through ${\it\beta}=2{\rm\pi}/{\it\lambda}$ .

3.4 Azimuthally averaged perturbation kinetic energy

To complement the linear stability analysis, consideration is given to energy conversion from the base flow to growing disturbances via analysis of the volume integrated energetics of the predicted eigenmodes (Gill Reference Gill1982; Ménesguen, McWilliams & Molemaker Reference Ménesguen, McWilliams and Molemaker2012). The fluctuation kinetic energy per unit mass is defined as $k^{\,\prime }=({u_{z}^{\prime }}^{2}+{u_{r}^{\prime }}^{2}+{u_{{\it\phi}}^{\prime }}^{2})/2$ , and an equation describing its evolution is obtained by taking the dot product of $\boldsymbol{u}^{\prime }$ with (3.17), given that $\partial _{t}k^{\,\prime }=(1/2)\partial _{t}(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime })=\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\partial _{t}\boldsymbol{u}^{\prime }$ . The azimuthal average of the equation is then taken (denoted by overbars), resulting in the azimuthally averaged perturbation kinetic energy equation

where

and $\unicode[STIX]{x1D634}_{ij}^{\,\prime }$ is the instantaneous rate-of-strain tensor, which in cylindrical ( $z$ , $r$ , ${\it\phi}$ ) coordinates is written

For brevity we have adopted tensor notation to represent the double-dot-product of the rate-of-strain tensor (i.e. $\unicode[STIX]{x1D634}_{ij}^{\,\prime }\unicode[STIX]{x1D634}_{ij}^{\,\prime }$ represents the sum of the square of each element in $\unicode[STIX]{x1D634}_{ij}^{\,\prime }$ ).

In (3.19), the left-hand side is the rate of change of azimuthally averaged perturbation kinetic energy due to time dependence in the perturbation field, and the first group of terms (in square brackets) on the right-hand side describes the rate of change of azimuthally averaged perturbation kinetic energy due to its advection on the base flow in the $r$ – ${\it\phi}$ plane. Omitted from the right-hand side are the divergence terms, which describe transport of azimuthally averaged kinetic energy on the $r$ – ${\it\phi}$ plane due to pressure fluctuations and viscous stresses, but which identically contribute nothing to the overall change in the integral of $\overline{k^{\,\prime }}$ over the flow domain. The next group of terms on the right-hand side describe the rate of $\overline{k^{\,\prime }}$ production from shear in the base flow. Following this is a term describing the rate of $\overline{k^{\,\prime }}$ dissipation due to viscous stresses, and finally a perturbation buoyancy flux term. In general, this last term would be written $-(4\mathit{Ra}^{1/5}/\mathit{Pr}\,Q^{2})[{\hat{g}}_{r}\overline{u_{r}^{\prime }{\it\theta}^{\prime }}+{\hat{g}}_{{\it\phi}}\overline{u_{{\it\phi}}^{\prime }{\it\theta}^{\prime }}+{\hat{g}}_{z}\overline{u_{z}^{\prime }{\it\theta}^{\prime }}]$ , with ${\hat{g}}_{z}$ , ${\hat{g}}_{r}$ and ${\hat{g}}_{{\it\phi}}$ being components of the unit vector in the direction of gravity, but in (3.19) only the axial component is included as gravity acts solely in the negative axial (downward) direction in the present study.

The perturbation buoyancy flux describes the reversible exchange of energy between potential energy in the axisymmetric base flow and kinetic energy in the perturbation field. This may be demonstrated by considering the evolution of potential energy. Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) combine their buoyancy transport equation and the incompressibility constraint to obtain an evolution equation for total potential energy; here instead the local potential energy evolution is considered. Defining a potential energy density $\mathscr{E}_{p}={\it\rho}z$ such that (3.10) becomes $E_{p}\equiv g\int \mathscr{E}_{p}\,\text{d}\mathscr{V}$ , which has axisymmetric mean and fluctuating perturbation components $\bar{\mathscr{E}}_{p}$ and $\mathscr{E}_{p}^{\prime }$ , respectively, appendix A details the derivation of the evolution equation for azimuthally averaged potential energy density in the mean flow. In dimensionless form this equation is written

where a Froude number $\mathit{Fr}=R{\it\Omega}^{2}/g$ has been introduced. The first two terms on the right-hand side of (3.22) describe advective transport of potential energy density by the mean flow and perturbation, respectively. The third term describes vertical advective buoyancy fluxes in the mean flow and perturbation field and the final term describes the effect of thermal dissipation in the mean flow. When comparing the prefactors to the dimensionless azimuthally averaged vertical advective buoyancy flux terms in the perturbation kinetic energy equation (3.19) and (3.22), the Froude number captures the ratio of the quantities used to scale kinetic and potential energy in the present study. Notice that in (3.22) the perturbation buoyancy flux term is negative, whereas the corresponding term is positive in (3.19), reflecting the reversible exchange role of these terms. Crucially, vertical advective buoyancy flux facilitates a conversion of potential energy in the axisymmetric mean flow to kinetic energy in a non-axisymmetric perturbation.

The power of equation (3.19) is revealed when we consider its application to the predicted linear instability modes of the axisymmetric base flows. These modes necessarily exhibit exponential growth (at a rate ${\it\sigma}$ ), and it is straightforward to show then that

where total perturbation kinetic energy $E_{k}=\int _{{\it\Omega}}\overline{k^{\,\prime }}\,\text{d}\mathscr{V}$ . Hence when normalised by $2E_{k}$ , (3.19) equates to the exponential growth rate of eigenmodes arising from the linear stability analysis, and therefore, the integrals of each term normalised by $2E_{k}$ provide insight into their individual contributions to the overall growth or decay of each instability mode. Furthermore, the spatial variation of each term on the $z$ – $r$ plane will provide insight into the local features of the underlying axisymmetric velocity and temperature fields that most actively contribute to the instability growth, facilitating a quantitative basis for the classification of instability modes reported herein.

3.5 Spatial and temporal discretisation

The governing equations (3.1)–(3.3) are solved in cylindrical coordinates using a nodal spectral element method in space, and a third-order scheme based on backwards differentiation is employed for time integration (Karniadakis, Israeli & Orszag Reference Karniadakis, Israeli and Orszag1991). The same scheme is also employed for solution of the linearised equations (3.16)–(3.18). Lagrangian tensor-product polynomial shape functions are imposed upon each macro-element. The polynomial order is varied to control spatial resolution and interpolated at the Gauss–Lobatto–Legendre quadrature points. The cylindrical formulation of the solver employed here has been validated in previous studies (Sheard & Ryan Reference Sheard and Ryan2007; Sheard Reference Sheard2009) and was used recently to study axisymmetric rotating radial horizontal convection by Hussam et al. (Reference Hussam, Tsai and Sheard2014).

The computational domain in the $z$ – $r$ plane is discretised into quadrilateral elements. A rectangular mesh comprising $1320$ elements was constructed to discretise the meridional semi-plane. To resolve the flow accurately, the grid size was much smaller in the vicinity of the side and bottom walls, particularly the heated boundary, with coarser grid spacing in the interior. A grid resolution study was undertaken to determine a suitably accurate element polynomial degree. The test was performed at the upper end of the parameter range of this study, $\mathit{Ra}=10^{9}$ and $Q=60$ as a stringent test of the mesh resolution. Three parameters are monitored for convergence: domain integral of temperature ( ${\it\theta}_{avg}$ ) and square of velocity magnitude ( $L_{2}$ ), and the growth rate of the leading eigenmode obtained from the linear stability analysis at the dominant azimuthal wavenumber ( ${\it\sigma}$ ). As shown in table 1, ${\it\theta}_{avg}$ and $L_{2}$ are highly resolved: ${\it\theta}_{avg}$ errors between polynomial degree $N_{p}$ and $N_{p}+1$ decrease from $3.0\times 10^{-2}\,\%$ down to $4.4\times 10^{-3}\,\%$ from $N_{p}=4$ to $6$ , while $L_{2}$ errors decrease from $4.8\times 10^{-4}\,\%$ down to $1.6\times 10^{-4}\,\%$ . However, the higher resolution required to resolve the instability eigenmodes is reflected in the ${\it\sigma}$ case, with errors decreasing from $2.8\,\%$ at $N_{p}=4$ down to $3.7\times 10^{-2}\,\%$ . Hence $N_{p}=7$ is used hereafter to constrain errors to less than $O(0.1\,\%)$ .

For time integration of equations (3.2)–(3.3), the advection/convection terms are concurrently solved explicitly, followed by a projection of the velocity field onto a divergence-free space, and finally implicit solves for velocity components and temperature. This procedure extends the backwards differentiation algorithm of Karniadakis et al. (Reference Karniadakis, Israeli and Orszag1991) to the coupled Boussinesq equations. The temperature transport formulation of the present code has been validated in studies on buoyancy-driven flows (Sheard & King Reference Sheard and King2011; Hussam et al. Reference Hussam, Tsai and Sheard2014).

Table 1.

Convergence of the integral temperature magnitude throughout the domain ${\it\theta}_{avg}$ , integral velocity magnitude $L_{2}$ and maximum growth rate ${\it\sigma}$ for different polynomial degree $N_{p}$ at $\mathit{Ra}=10^{9}$ and $Q=60$ .

5.1 Eigenmode energetics along instability mode branches

Along each of the mode branches revealed in figures 4 and 6, the energetics of the azimuthally averaged eigenvector fields were analysed by evaluating each term of equation (3.19). In every case considered in this work, the same four terms consistently had the largest magnitudes across the $12$ calculated terms: the next largest term was at most $5\,\%$ of the dominant term in all cases, with a mean of $2.1\,\%$ . As an example, table 2 displays the computed integral summations of each term for the leading instability mode at $Q=5$ and $k=12$ , with the dominant terms highlighted.

Table 2.

Volume integrated values of terms in (3.19) for the leading instability eigenmode with ${\it\beta}=12$ for a flow with $Q=5$ and $\mathit{Ra}=10^{9}$ . The integrated term values are scaled by $R^{5}{\it\Omega}^{3}$ . For each term, the two right-most columns respectively include the percentage contributions of each term to the sum and range ( $\max {-}\min$ ) of the terms contributing to $\partial _{t}\overline{k^{\,\prime }}$ . Highlighted in bold are the four terms of largest absolute value contribution.

For brevity, nomenclature is assigned to each of the four dominant terms, as per table 3. The first two terms describe the production of perturbation kinetic energy in the $r$ - ${\it\theta}$ plane due respectively to horizontal shear in, and azimuthal rotation of, the base flow, and are denoted by $P_{1}$ and $P_{2}$ . The third term, denoted by $D$ , describes the viscous dissipation of perturbation kinetic energy and the fourth term, denoted by $B$ , describes the transfer from perturbation potential energy to kinetic energy. Note that due to the absence of the integration operators $\langle \cdots \,\rangle$ from the terms in table 3 compared to the corresponding terms in table 2, $P_{1}$ , $P_{2}$ , $D$ and $B$ each represent a spatially varying field on the $z$ – $r$ plane.

Figure 7.

(a–d) Plots of the respective percent-of-range contributions of $\langle P_{1}\rangle$ , $\langle P_{2}\rangle$ , $\langle D\rangle$ and $\langle B\rangle$ against $Q$ , for each of the mode branches identified in figure 4.

Table 3.

Nomenclature assigned to the four dominant terms from (3.19) and table 2 in terms of their absolute contributions to the rate of change of azimuthally averaged perturbation kinetic energy.

Percent-of-range contributions of each of the four dominant terms were evaluated along the loci of maximum growth rate for the mode branches from figures 4 and 6. These are plotted in figure 7. Several different instability behaviours are suggested by the markedly different proportional contributions of the four dominant terms in these plots:

1. (1) At both $Q\gtrsim 30$ (M5) and $Q\approx 18$ (M4), the two production terms $\langle P_{1}\rangle$ and $\langle P_{2}\rangle$ are weakly negative and positive, respectively (magnitudes being less than approximately $6\,\%$ ), and sum to approximately zero. The dissipation term $\langle D\rangle$ is strongly negative ( ${\approx}-47\,\%$ ), while the buoyancy term $\langle B\rangle$ is slightly stronger ( ${\approx}52\,\%$ ) and positive.

2. (2) The M1 and M3 modes at smaller $Q$ values exhibit similar proportions to M4 and M5 for the $\langle D\rangle$ and $\langle B\rangle$ terms, but are notably different in the production terms. In addition to being stronger (M1 and M2 having production term magnitudes up to $12\,\%$ and $18\,\%$ , respectively), the terms are opposite-signed, with $\langle P_{1}\rangle$ predominantly positive and $\langle P_{2}\rangle$ negative.

3. (3) The M2 mode is different again: production terms $\langle P_{1}\rangle$ and $\langle P_{2}\rangle$ are strong ( ${\approx}50\,\%$ ), with $\langle P_{1}\rangle$ positive and $\langle P_{2}\rangle$ negative. This mode exhibits the weakest proportions of dissipation (between $-7\,\%$ and $-14\,\%$ ) and buoyancy (between approximately $12\,\%$ and $25\,\%$ ) amongst all mode branches.

4. (4) The M1 mode exhibits the most significant change in the relative contributions of the terms over intermediate $Q$ values $5\lesssim Q\lesssim 15$ . At low $Q$ , the mode exhibits modest production term contributions, a strong negative dissipation contribution, and a stronger positive buoyancy contribution, while at higher $Q$ , the distribution is similar to the M2 mode.

In order to confirm if the categorisation of instability modes based on their energetics is valid, attention is turned to the eigenvector fields in the subsections to follow.

5.2 Eigenmode structure and dominant energetics fields

The first mode branch to be considered is the M1 branch. In contrast to the other mode branches described herein, the energetics analysis in figure 7 shows that the M1 branch experiences a significant change in the proportional contributions of individual terms to the evolution of azimuthally averaged perturbation kinetic energy as a function of $Q$ : when $Q\lesssim 5$ the M1 branch has energy evolution contributions consistent with the M3 branch; the shift in the distribution of contributions occurs over $5\lesssim Q\lesssim 15$ ; and for $Q\gtrsim 15$ the M1 branch has energy evolution contributions resembling the M2 branch. Figure 8 illustrates the eigenmode structure and the integral contributions of all $12$ terms from table 3 for the M1 branch at values of $Q$ capturing each of the aforementioned behaviours.

Figure 8.

(a–i) Structure of the M1 eigenmode; contour plots of perturbation temperature, azimuthally averaged perturbation kinetic energy evolution viscous dissipation ( $D$ ) and available potential energy conversion fields ( $B$ ) at $Q$ and ${\it\beta}$ as shown. Arbitrary contour levels equispaced about zero are plotted to elucidate the perturbation field structure. For each $Q$ – ${\it\beta}$ combination, the same contour levels are used for $D$ and $B$ fields. (j) Plots the integral contributions of each of the terms contributing to the azimuthally averaged perturbation kinetic energy equation, for each of the $Q$ – ${\it\beta}$ combinations in (a–i). White, light grey and dark grey bars represent (a–c)  $Q=2,{\it\beta}=5^{\circ }$ , (d–f)  $Q=10,{\it\beta}=15^{\circ }$ , and (g–i)  $Q=20,{\it\beta}=20^{\circ }$ , respectively. From left to right, the plot follows the order of terms shown in table 2, while the four dominant terms, $\langle P1\rangle$ , $\langle P2\rangle$ , $\langle D\rangle$ and $\langle B\rangle$ are labelled. The first two terms arise from advection of $\overline{k^{\,\prime }}$ on the mean flow, the next $8$ terms arise from production of $\overline{k^{\,\prime }}$ due to velocity shear, as labelled. The contributions to $\partial \overline{k^{\,\prime }}/\partial t$ are scaled by ${\it\nu}^{3}/R^{4}$ .

The first case, $Q=2$ and ${\it\beta}=5$ , demonstrates that the eigenmode is concentrated adjacent to the hot outer part of the forcing boundary and up the outer side wall (occupying approximately the outer radial $25\,\%$ of the base and lower $40\,\%$ of the side wall). This is reflected by the ${\it\theta}^{\prime }$ , $D$ and $B$ fields shown in figure 8(a–i). The dissipation field $D$ is strongest at the boundaries in the vicinity of the junction between base and side wall. This behaviour is consistent with dissipation due to viscous shear in perturbation velocity boundary layers in this vicinity. In terms of the $B$ field, regions of high production extend from the hot corner. At $Q=2$ , the mode is only very weakly unstable (ref. the small positive peak growth rate in figure 5). Correspondingly, in figure 8(j) the integrals of the two dominant terms, $D$ and $B$ , are seen to be approximately equal and opposite.

In the region of shifting mode behaviour ( $Q=10$ and ${\it\beta}=15$ , figure 8 d–f), the eigenmode is found to become narrower and extend further in the vertical direction. In notable contrast to the $Q=2$ case, here the perturbation field structure is displaced from the side wall, the dissipation field $D$ is again strongest at the bottom boundary at the same location as the temperature field features and the $B$ field is dominated by a region of strong energy production extending vertically from the base through much of the interior. At $Q=20$ , the eigenmode structure is positioned at approximately $80\,\%$ of the tank radius from the axis. The $B$ field exhibits a similar dominant vertical structure to the $Q=10$ case, but the structure is broader.

In figure 8(j) the production term contributions $\langle P1\rangle$ and $\langle P2\rangle$ are consistently positive and negative, respectively, and approximately cancel. While not shown in this paper, the authors consistently found the fields for the $P1$ and $P2$ production terms to be similarly structured, but opposite signed. Hence the $P1$ and $P2$ fields typically sum to approximately zero both locally and globally. With increasing $Q$ , the M1 branch exhibits a strong increase in the magnitude of $\langle P1\rangle$ and $\langle P2\rangle$ . At $Q=2$ , the $\langle D\rangle$ and $\langle B\rangle$ contributions have a significantly greater magnitude than $\langle P1\rangle$ and $\langle P2\rangle$ , but by $Q=20$ this is reversed.

The M2 branch is represented at $Q=10$ with dominant wavenumber ${\it\beta}=5$ in figure 9. In figure 9(a–f), fields ${\it\theta}^{\prime }$ , $D$ and $B$ are plotted, along with the axial, radial and azimuthal components of vorticity in the eigenvector field ( ${\it\omega}_{z}^{\prime }$ , ${\it\omega}_{r}^{\prime }$ and ${\it\omega}_{{\it\phi}}^{\prime }$ , respectively). Strong similarities are observed between the M2 branch eigenmode visualised here and the M1 branch at $Q=20$ in figure 8 (moreover, the energetics data in figure 7 indicates that the behaviour at $Q=20$ will be maintained over $15\lesssim Q\lesssim 30$ ). However, here the M2 eigenmode structures are approximately three times broader than the $Q=20$ M1 eigenmode, likely reflecting the lower dominant wavenumbers of the M2 branch (here ${\it\beta}=5$ ) in comparison to the M2 branch beyond $Q\approx 15$ ( ${\it\beta}\gtrsim 15$ ).

Figure 9.

Representative structure of the M2 eigenmode at $Q=10$ and ${\it\beta}=5$ . (a–f) Plots fields of ${\it\theta}^{\prime }$ , $D$ , $B$ and perturbation vorticity components ${\it\omega}_{z}^{\prime }$ , ${\it\omega}_{r}^{\prime }$ and ${\it\omega}_{{\it\phi}}^{\prime }$ . For ${\it\theta}^{\prime }$ , $D$ and $B$ , contour levels and shading are as per figure 8. The same set of arbitrary contour levels are used for each vorticity component. (g) Plots contributions to $\partial _{t}\overline{k^{\,\prime }}$ as per figure 8.

The perturbation kinetic energy contributions plotted in figure 9(b) demonstrate that the production terms sum to a weakly negative contribution. The gain from $\langle B\rangle$ is more than twice the deficit due to viscous dissipation through $\langle D\rangle$ , reflecting the importance of $\langle B\rangle$ to the instability growth. The spatial distribution of $B$ reveals that this conversion of perturbation available potential energy to kinetic energy is strongest in a broad region extending almost the full depth of the enclosure at approximately $50\,\%$ of the radial distance from the axis. The vorticity fields reveal alternating-sign axial vorticity structures extending from the top to the base of the enclosure. The radial and azimuthal vorticity fields are both strongest in a thin layer adjacent to the bottom wall, which the $D$ field in figure 9(a–f) demonstrates is responsible for significant viscous dissipation, but because it is localised to a slender region, its integral $\langle D\rangle$ is modest, and does not inhibit instability growth. Revisiting the axisymmetric base flow for this case with $Q=10$ depicted in figure 2(g–i), the region occupied by the eigenmode structures has inclined isotherms with thermal gradient directed diagonally upward and radially outward. Under rotation, this flow pattern is susceptible to baroclinic instability, and the change in potential energy field ${\rm\Delta}\mathscr{E}_{p}$ in figure 2(g–i) is indeed found to be strongly negative in the region occupied by the eigenmode structures, supporting the assertion that this instability mode has a baroclinic origin. Furthermore, the available potential energy density field $\mathscr{E}_{a}$ for the axisymmetric base flow demonstrates that the flow has accumulated a significant amount of potential energy that cannot be released into the axisymmetric flow, but is available for release to the non-axisymmetric (three-dimensional) modes. Barkan et al. (Reference Barkan, Winters and Smith2013) demonstrated using numerical simulation that two-dimensional rotating horizontal convection in a rectangular enclosure accumulated significantly more available potential energy than its three-dimensional counterpart. The present results demonstrate a global linear instability mode with baroclinic features that may facilitate a conversion of accumulated axisymmetric available potential energy to perturbation kinetic energy in a non-axisymmetric flow state.

Figure 10.

Representative structure of the M5 eigenmode at $Q=50$ and ${\it\beta}=25$ . Contour levels and shading are as per figure 9.

A representative eigenmode from the M5 branch is displayed in figure 10 for $Q=50$ and ${\it\beta}=25$ . Figure 10(g) demonstrates the unusual feature amongst these mode branches that $\langle P1\rangle <0$ and $\langle P2\rangle >0$ which is contrary to all but the fleeting M4 branch. The eigenmode structures are concentrated in a narrow region adjacent to the side wall. The radial vorticity component is small compared to the axial and azimuthal components, both of which are prominent.

Figure 11.

Ratio of the sum of production terms to the perturbation vertical advective buoyancy flux for each of the mode branches identified in figure 4. Symbols are as per figure 7.

In addition to thermal instability mechanisms, consideration was also given to the possibility that hydrodynamic shear-layer instabilities may play a role in this system. For instance, the development of swirling flow about a vertical axis of rotation may invite shear-layer instabilities on horizontal planes. The Rayleigh–Kuo criterion (Rayleigh Reference Rayleigh1879; Kuo Reference Kuo1949) is a necessary but not sufficient condition for instability that requires the horizontal gradient of absolute vorticity to change sign somewhere within the domain. The base flows obtained in § 4 were tested at various depths across the investigated range of $Q$ against the Rayleigh–Kuo criterion. All $Q\geqslant 4$ satisfied this criterion at radial positions $0.85\lesssim r<1$ . This does appear to correspond to the location of eigenmode structures in this case. However, it would be expected that an instability mechanism of this kind would present in the azimuthally averaged perturbation kinetic energy evolution analysis through a significant net contribution of the $P1$ and $P2$ terms, as these describe the production of perturbation kinetic energy due to perturbation velocity on the horizontal plane and the radial variation in azimuthal velocity (ref. table 3). No such net contribution from these production terms is found for this mode nor any of the mode branches elucidated in this study. The dominance of thermal instability mechanisms in this system is discussed in detail in the section to follow.

5.3 The dominant mechanism of perturbation kinetic energy growth

It is apparent from the plots of the contributions of terms of $\partial \overline{k^{\,\prime }}/\partial t$ that the two terms largest in magnitude, $P1$ and $P2$ , approximately offset each other, and that the other production terms are significantly smaller in magnitude. Hydrodynamic instabilities, such as shear-layer instability, would present through these production terms. On the other hand, it was shown in § 3.4 and (3.22) that buoyancy-driven instabilities, such as baroclinic or Rayleigh–Bénard type instabilities, will present through vertical advective buoyancy flux (term $B$ ). Insight into the dominant mechanism driving perturbation kinetic energy growth (and hence growth of the underlying instability modes) is gleaned from consideration of the relative contributions of the sum total of the production terms to the buoyancy flux term. Figure 11 plots this ratio against $Q$ for each mode branch. The largest value of this ratio is found to be approximately $5.5\,\%$ for the M3 mode. As the buoyancy flux contribution was always positive, any negative values of the ratio $\sum \langle P\rangle /\langle B\rangle$ correspond to negative net contributions from the production terms. Modes M1 and M2 both exhibit this behaviour. The M4 and M5 branches are approximately zero, corresponding to a near-zero contribution to instability growth by the production terms in these cases. Given that the production terms contribute at best approximately $5\,\%$ of the contribution of vertical advective perturbation buoyancy flux, it is apparent that the instability modes identified in this study consistently have a thermal origin, being driven by the conversion of available potential energy in the axisymmetric mean flow to kinetic energy in the perturbation field.

Recently, Tsai et al. (Reference Tsai, Hussam, Fouras and Sheard2016) confirmed via a local one-dimensional linear stability analysis that the instability leading to unsteady flow in planar horizontal convection had a thermal origin – specifically of a Rayleigh–Bénard kind at the hotter end of the heated bottom boundary. Their study also demonstrated the near insensitivity of the instability to velocity shear in the boundary layer adjacent to the heated boundary. This is consistent with the observations reported herein of negligible net contribution to perturbation kinetic energy growth from the production terms, relative to the contribution from perturbation buoyancy flux, across all instability branches.

The role played by the axisymmetric mean thermal field towards the exchange of vertical advective perturbation buoyancy flux is now considered. Appendix B describes the derivation of an evolution equation for azimuthally averaged buoyancy flux,

Briefly, this equation is obtained by premultiplying the perturbation thermodynamic equation (3.18) by the vertical perturbation velocity $u_{z}^{\prime }$ , rearranging and azimuthally averaging. The first two terms on the right-hand side (5.1a ) describe vertical buoyancy flux production by buoyancy flux through vertical and radial gradients of axial velocity in the base flow. The contributions arising from the perturbation pressure gradient (5.1b ) and both thermal and viscous dissipation (5.1c ) follow. The transport of the buoyancy flux field on the $z$ – $r$ plane is described by (5.1d ), and the production of buoyancy flux through interaction of the perturbation velocity field and thermal gradients in the axisymmetric mean flow is described by (5.1e ).

Volume integration of the two thermal production terms (5.1e ) for each of the instability modes identified in this study will reveal the nature of the instability mechanism – whether either or both of the radial (horizontal) and vertical thermal gradients actively drive the exchange of perturbation buoyancy flux feeding the instability growth – and the corresponding distributions of these production terms on the $z$ – $r$ plane will elucidate the regions of the flow producing this activity.

It is pertinent to discuss these two terms in the context of thermal instability: the change of potential energy criterion outlined in § 3.1 demonstrated that under an upward and radially outward thermal gradient, movement of a fluid element from a higher cooler position to a lower hotter position could incite baroclinic instability via exchange of potential energy in the background flow to kinetic energy in a disturbance. With respect to the second term of (5.1e ), downward (negative $u_{z}^{\prime }$ ) and radially outward (positive $u_{r}^{\prime }$ ) motion from a cooler to a hotter region is facilitated by a positive $\partial _{r}\bar{{\it\theta}}$ . Gathering these returns the form of the second term, $-u_{z}^{\prime }u_{r}^{\prime }\partial _{r}\bar{{\it\theta}}$ , and hence eigenmodes exhibiting regions with a positive value for this term may indicate a baroclinic instability mechanism. Consider the first term of (5.1e ) in the context of Rayleigh–Bénard instability: upward or downward vertical motion of a fluid element (i.e. positive ${u_{z}^{\prime }}^{2}$ ) through an adverse vertical thermal gradient (negative $\partial _{z}\bar{{\it\theta}}$  ) causes a reduction in potential energy, facilitating instability growth. Hence positive values of the first term, $-{u_{z}^{\prime }}^{2}\partial _{z}\bar{{\it\theta}}$ may be representative of a Rayleigh–Bénard instability mechanism.

Figure 12.

Contour plots showing the production of azimuthally averaged vertical advective buoyancy flux by spatial gradients in the axisymmetric thermal field, for select cases from figures 8 to 10. Solid and dashed contour lines separately depict $-\overline{{u_{z}^{\prime }}^{2}}\,\partial _{z}\bar{{\it\theta}}$ and $-\overline{u_{z}^{\prime }u_{r}^{\prime }}\,\partial _{r}\bar{{\it\theta}}$ , respectively, and their sum is shown by the flooded contours. Only positive values are plotted, with darker shading representing larger values. In each frame, contour levels are equispaced between zero and an arbitrarily maximum magnitude. (a) M1: $Q=2$ , ${\it\beta}=5^{\circ }$ , (b) M2: $Q=10$ , ${\it\beta}=5^{\circ }$ , (c) M1: $Q=20$ , ${\it\beta}=20^{\circ }$ , (d) M5: $Q=50$ , ${\it\beta}=25^{\circ }$ .

Figure 13.

Growth rate as a function of azimuthal wavenumber ${\it\beta}/Q$ for $10^{6}\leqslant \mathit{Ra}\leqslant 10^{9}$ for $Q=10$ . The highlighted symbols show the interpolated maximum growth rate for each of the M1, M2 and M3 mode branches.

Figure 12 plots these thermal mean flow buoyancy flux production terms for four cases from figures 8 to 10. These show the various dominant modes with increasing $Q$ . From figure 12(a), it is apparent that the M1 mode is initially an instability apparently produced by an adverse thermal gradient adjacent to the bottom wall at the outer hot end, and adjacent to the enclosure side wall. The local Rayleigh number across the thermal boundary layer at the radial location where the left-most local maximum of buoyancy flux production occurs was calculated to be $\mathit{Ra}_{{\it\theta}}=738$ . This is somewhat lower than the accepted critical value for Rayleigh–Bénard convection ( $\mathit{Ra}_{crit.}=1707.76$ ; Reid & Harris Reference Reid and Harris1958), which is consistent with the stable axisymmetric flows found at this Rayleigh number. Hence figure 12(a) may depict production of perturbation buoyancy flux associated with a thermal instability of the side-wall jet.

As per earlier observations, in figure 12(b) the M2 branch is found to be similar to the high- $Q$ M1 branch, having characteristics consistent with baroclinic instability.

Finally, the M5 branch at $Q=50$ with ${\it\beta}=25$ shown in figure 12(d) has a thin vertical region adjacent to the lower half of the side wall exhibiting strong perturbation buoyancy flux production by the locally adverse thermal gradient. In this region the distribution of these production terms is similar to those of the M1 branch at $Q=2$ with ${\it\beta}=5$ in figure 12(a), though the absence of an additional region of positive values of the vertical production term above the bottom boundary likely explains the less evident similarities between these cases in the eigenmodes (cf. figures 8(a–c) and 10).

The next section addresses the Rayleigh number dependence of the flow stability.

5.4 Rayleigh number dependence

The Rayleigh number dependence of the maximum growth rate and the associated disturbances is now considered. Figure 13 shows the predicted growth rate against ${\it\beta}/Q$ for $Q=10$ . The case of $Q=10$ is considered as it was particularly active to instability, exhibiting each of the M1, M2 and M3 branches at $\mathit{Ra}=10^{9}$ , permitting the Rayleigh number dependence of several of the mode branches to be examined. Decreasing Rayleigh number results in a decrease in the number of instability modes from three to two modes, with the suppression of the thermal M3 mode. This is not surprising given that thermal instability is expected to be strongly dependent on Rayleigh number, so reducing Rayleigh number suppresses the M3 mode. It is found that decreasing Rayleigh number from $\mathit{Ra}=10^{9}$ to $10^{7}$ leads to a significant reduction in the maximum growth rate of the dominant instability, the M2 mode, and a switch of the dominant instability mode to the M5 mode. In fact, the maximum growth rate decreases by $46\,\%$ and $90\,\%$ from $\mathit{Ra}=10^{9}$ down to $\mathit{Ra}=10^{8}$ and $10^{7}$ , respectively. However, the peak wavenumber where the maximum growth rate occurs is approximately independent of Rayleigh number. For the M2 mode, instability emerges with peak wavenumber at ${\it\beta}/Q=0.55$ regardless of Rayleigh number, while it emerges with peak wavenumber ranging from $1.3$ to $1.5$ for the M3-mode instability. The suppression of the M2 mode identified earlier as having features consistent with baroclinic instability with decreasing Rayleigh number illuminates the role of the base flow in determining the sensitivity of the flow to baroclinic instability. It was earlier shown that for $\mathit{Ra}=10^{9}$ , baroclinic instability is supported through much of regime II, but is suppressed in regime I at $Q=60$ . Taking the $Q_{\text{I-II}}$ criterion from Hussam et al. (Reference Hussam, Tsai and Sheard2014), it is found that from $\mathit{Ra}=10^{8}$ to $10^{7}$ to $10^{6}$ , the threshold rotation parameter decreases from $Q_{\text{I-II}}=29$ to $16$ to $8.6$ . Hence at $Q=10$ the Rayleigh numbers producing baroclinic instability correspond to those in which $Q_{\text{I-II}}>10$ (i.e. lying within regime II). $\mathit{Ra}=10^{6}$ is stable at $Q=10$ , and lies within the very strongly rotating regime I.

This is a significant outcome from this study: following the definitions of Hignett et al. (Reference Hignett, Ibbetson and Killworth1981) baroclinic instability is supported in rotating horizontal convection under strong rotation, but not very strong rotation. It is the strong rotation regime that is most relevant to Earth’s oceans.