3.1. Flow visualisations

A first impression of the 3-D spatial organisation of the wake vortices is provided by the isosurfaces of the vertical component of vorticity ($\omega _z=(\boldsymbol {\nabla } \times \boldsymbol {u})_z$) in figure 2(a–c) for cases BuInf, Bu25 and Bu1, respectively. At moderately strong stratification ($Fr\le O(0.1)$), vertical overturning motions are suppressed to approximately one order of magnitude smaller than horizontal flows, with the latter well represented by $\omega _z$.

Figure 2. Visualisation of the isosurfaces of $| \omega _z | D/U_{\infty }=1.5, 1.0, 3.0$ for BuInf, Bu25 and Bu1 (top to bottom), respectively. Red and blue indicate positive (anticyclonic since $f_{c}$ is negative) and negative (cyclonic) vorticity, respectively. Note the vertical axis is normalised with the height $h$ of the hill, which is approximately $0.3$ times the base diameter $D$. The three snapshots (a–c) are taken at time instants $tU_{\infty }/D=252.0, 221.6, 152.0$, respectively.

In the near-wake region ($x/D<3$), the $\omega _z$ isosurface is space filling and multiscale, indicating greater turbulence intensity than further downstream. But after $x/D=3$, the vortices quickly organise into spatially compact coherent structures that persist to the end of the computational domain with little change. Rotation clearly influences the shape and size of the structures. For case BuInf in figure 2(a), the structures are slanted ‘tongues’ with horizontal dimensions greater than their height. As the strength of rotation is increased, the horizontal size shrinks and the height increases as in figure 2(b) for Bu25. With further increase of rotation, the vortex structures become aligned with the vertical axis, appearing as tall ‘columns’ that extend from the flat bottom to the hill height, in figure 2(c).

The wake vortices are composed of cyclones (rotation is in the same direction as that of the frame) and anticyclones. It is found in figure 2(b,c) that, with the presence of rotation, the cyclones (negative vortices in blue) are different from anticyclones (positive vortices in red) with regards to size, shape and vorticity distribution, as will be elaborated in later sections. In figure 2(c), cyclones are taller and thinner than anticyclones and, as will be shown, have stronger interior vorticity.

Figure 3(a–c) shows the instantaneous $\omega _z$ at several heights for all three cases. In all three cases, a pattern of Q2-D Kármán vortices is distinct in all planes at $z/h=0.12, 0.25, 0.75$. The flow is akin to a vortex wake rather than its unstratified counterpart of a 3-D turbulent wake (Garcia-Villalba et al. Reference Garcia-Villalba, Li, Rodi and Leschziner2009).

Figure 3. (a–c) Instantaneous and (d) time-averaged $\omega _z$ at horizontal planes at $z/h=0.75, 0.50, 0.12$ (from top to bottom in each panel). (a,d) Cases BuInf, (b) Bu25, (c) Bu1. In (a–c), the snapshots are taken at non-dimensional time instants $tU_{\infty }D= 77.9, 123.9,63.3$ counted after the fully developed state is reached, respectively.

The mean vertical vorticity ($\langle \omega _z \rangle$) of case BuInf is shown in figure 3(d). For the lower two planes, the mean looks very similar to those obtained in low Reynolds number 2-D cylinder wakes, such as in Barkley (Reference Barkley2006) and Mittal (Reference Mittal2008), both at $Re_D=100$. Such similarity supports the Q2-D feature of the hill wake at $Fr=0.15$. At the same time, there is a notable difference: the flow in each horizontal place that cuts the hill does not represent an independent 2-D flow around a cylinder with the local hill diameter. Among different planes at different heights in figure 3(d), the length of the attached shear layer (with dark colours), is approximately a constant, regardless of the variation of the local hill diameter. This is consistent with the fact that the VS frequency is a global constant which will be discussed in detail in the next section, as the shear layer length is correlated to the shedding frequency (Williamson & Brown Reference Williamson and Brown1998).

In the SPOD analysis of the next section, the vertical vorticity on different horizontal planes ($\omega _z(x,y,t;z)$) is chosen as the quantity of interest since the large-scale wake structures involve predominantly Q2-D vortical motion (figures 2–3), the shedding and evolution processes of which are well represented by $\omega _z$. Moreover, the focus of this work is on the influence of increasing rotation strength on the horizontal motions, which tends to constrain the flow to be around the vertical axis at the large scales and significantly enhances $\omega _z$, as will be shown later. In the rotating cases Bu25 and Bu1, the stability of the anticyclones will also be studied (§ 5) with $\omega _z$ being one of the most important metrics, hence we apply $\omega _z$ rather than other vortex identification criteria for overall consistency.

As seen in figures 2–3, $\omega _z$ structures exhibit dissimilar vertical organisation at different levels of rotation, although Kármán vortices are present in each horizontal plane. Owing to stratification, a vertical length scale of $U_{\infty }/N = Frh = 0.15h$ is introduced to the flow, and whether vortex structures remain coherent over vertical distances larger than $U_{\infty }/N$ regardless of rotation (suggested affirmatively by figure 2) needs quantitative investigation. To that end, we perform SPOD on vertically offset horizontal planes at $z/h=0.12, 0.25,0.50,0.75$ ($N_x\times N_y\times N_z=1538\times 1280\times 1$) and the vertical centre plane at $y=0$ ($N_x\times N_y \times N_z=1538\times 1\times 320$) to allow the choice of different dominant (VS) frequencies at different heights by the flow and avoid imposing a priori a global frequency.

3.2. Spectral proper orthogonal decomposition and its numerical implementation

Proper orthogonal decomposition (POD) is a matrix-factorisation-based modal decomposition of complex systems introduced into fluid mechanics by Bakewell & Lumley (Reference Bakewell and Lumley1967) and Lumley (Reference Lumley1967, Reference Lumley1970).

Consider a statistically stationary square-integrable multi-variable signal $\boldsymbol {q}(\boldsymbol {x},t) \in \mathcal {L}_{{{\boldsymbol{\mathsf{W}}}}}^2(\boldsymbol {\varOmega })$ with zero mean. Here, $\mathcal {L}^2_{{{\boldsymbol{\mathsf{W}}}}}$ is the Hilbert space equipped with a weighted inner product

(3.1)\begin{equation} (\boldsymbol{q}_1, \boldsymbol{q}_2)_{{{\boldsymbol{\mathsf{W}}}}} = \int_{\boldsymbol{\varOmega}} \boldsymbol{q}_2^{H}(\boldsymbol{x}) {{\boldsymbol{\mathsf{W}}}}(\boldsymbol{x}) \boldsymbol{q}_1(\boldsymbol{x})\,{\rm d}\kern0.07em\boldsymbol{x} ,\end{equation}

on a bounded domain $\boldsymbol {\varOmega }$, and $({\cdot })^{H}$ denotes Hermitian transpose. The weight matrix ${{\boldsymbol{\mathsf{W}}}}$ is Hermitian positive–definite and the weighted $2$-norm is defined as $\|\boldsymbol {q}_1\|_{{{\boldsymbol{\mathsf{W}}}}}=(\boldsymbol {q}_1,\boldsymbol {q}_1)_{{{\boldsymbol{\mathsf{W}}}}}^{1/2}$. The symbol $\langle {\cdot } \rangle _{E}$ denotes the ensemble average over all realisations, and it is equivalent to the time average under ergodicity. The goal of POD is to find an empirical function $\boldsymbol {\psi }(\boldsymbol {x})$ that solves the optimisation problem

(3.2)\begin{equation} \boldsymbol{\psi}(\boldsymbol{x}) = \arg \max_{\|\boldsymbol{\psi}\|_{{{\boldsymbol{\mathsf{W}}}}}=1} \langle (\boldsymbol{q},\boldsymbol{\psi})_{{{\boldsymbol{\mathsf{W}}}}} \rangle_{E}, \end{equation}

which defines $\boldsymbol {\psi }(\boldsymbol {x})$ as the function on which the projection of $\boldsymbol {q}(\boldsymbol {x},t)$ is maximised in the sense of the least squares. Since $\mathcal {L}^2_{{{\boldsymbol{\mathsf{W}}}}}$ is an infinite-dimensional space, a practical way to obtain the empirical function $\boldsymbol {\psi }(\boldsymbol {x})$ is to approximate it within a finite-dimensional space spanned by $\{\boldsymbol {\psi }^{(i)}\}_{i=1}^{M}$, where $\boldsymbol {\psi }^{(i)}$ is the $i$th spatial mode and $M$ is the order of truncation. It was shown in Holmes et al. (Reference Holmes, Lumley, Berkooz and Rowley2012) that the optimisation problem (3.2) can be converted to a Fredholm eigenvalue problem as

(3.3)\begin{equation} \mathcal{R} \boldsymbol{\psi}^{(i)}(\boldsymbol{x}) = \int_{\boldsymbol{\varOmega}} {{\boldsymbol{\mathsf{R}}}}(\boldsymbol{x},\boldsymbol{x}'){{\boldsymbol{\mathsf{W}}}}(\boldsymbol{x}') \boldsymbol{\psi}^{(i)}(\boldsymbol{x}')\,{\rm d}\kern0.07em\boldsymbol{x}' = \lambda^{(i)} \boldsymbol{\psi}^{(i)}(\boldsymbol{x}), \end{equation}

where $\mathcal {R}$ is a linear operator and ${{\boldsymbol{\mathsf{R}}}}(\boldsymbol {x},\boldsymbol {x}')= \langle \boldsymbol {q}(\boldsymbol {x}) \boldsymbol {q}^{H}(\boldsymbol {x}') \rangle _{E}$ is the two-point correlation tensor. Since ${{\boldsymbol{\mathsf{R}}}}$ is Hermitian positive–definite, its eigenvalues $\lambda ^{(i)}$ are real positive that each represents a fraction of the fluctuation energy, and the eigenvectors $\{\boldsymbol {\psi }^{(i)}\}_{i=1}^{M}$ form an orthogonal basis under the inner product (3.1).

In the framework of SPOD, the eigenvalue problem (3.3) is cast as

(3.4)\begin{equation} \mathcal{R} \boldsymbol{\psi}^{(i)}(\boldsymbol{x},t) = \int_{\boldsymbol{\varOmega}} \int_{-\infty}^{\infty} {{\boldsymbol{\mathsf{R}}}}(\boldsymbol{x},\boldsymbol{x}',t,t') {{\boldsymbol{\mathsf{W}}}}(\boldsymbol{x}') \boldsymbol{\psi}^{(i)}(\boldsymbol{x}',t')\,{\rm d}\kern0.07em\boldsymbol{x}'\,{\rm d} t' = \lambda^{(i)} \boldsymbol{\psi}^{(i)}(\boldsymbol{x},t), \end{equation}

and ${{\boldsymbol{\mathsf{R}}}}(\boldsymbol {x},\boldsymbol {x}',t,t') = \langle \boldsymbol {q}(\boldsymbol {x},t) \boldsymbol {q}^{H}(\boldsymbol {x}',t') \rangle _{E}$ is the two-point, two-time correlation tensor. With time homogeneity, it reduces to ${{\boldsymbol{\mathsf{R}}}}(\boldsymbol {x},\boldsymbol {x}',\tau )$ as a function of $\tau =t-t'$, and is the Fourier transform pair of the spectral density tensor ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x},\boldsymbol {x}',f) = \langle \hat {\boldsymbol {q}}(\boldsymbol {x},f) \hat {\boldsymbol {q}}^{H}(\boldsymbol {x}',f) \rangle _{E}$:

(3.5)\begin{equation} {{\boldsymbol{\mathsf{R}}}}(\boldsymbol{x},\boldsymbol{x}',\tau) = \int_{-\infty}^{\infty} {{\boldsymbol{\mathsf{S}}}}(\boldsymbol{x},\boldsymbol{x}',f) \exp({{\rm i} 2 {\rm \pi}f \tau})\, {\rm d} f. \end{equation}

Hence, $\boldsymbol {\phi }^{(i)}(\boldsymbol {x},f) = \boldsymbol {\psi }^{(i)} (\boldsymbol {x},t) \exp ({-{\rm i} 2 {\rm \pi}f \tau })$ will be the corresponding eigenmodes of the following eigenvalue problem:

(3.6)\begin{equation} \mathcal{S} \boldsymbol{\phi}^{(i)}(\boldsymbol{x},f) = \int_{\boldsymbol{\varOmega}} {{\boldsymbol{\mathsf{S}}}}(\boldsymbol{x},\boldsymbol{x}',f) {{\boldsymbol{\mathsf{W}}}}(\boldsymbol{x}') \boldsymbol{\phi}^{(i)}(\boldsymbol{x}',f) \, {\rm d}\kern0.07em\boldsymbol{x}' = \lambda^{(i)}(f) \boldsymbol{\phi}^{(i)}(\boldsymbol{x},f), \end{equation}

which will be solved separately for each frequency. Here, $\hat {\boldsymbol {q}}(\boldsymbol {x},f)$ denotes the Fourier mode of ${\boldsymbol {q}}(\boldsymbol {x},t)$ at frequency $f$, and can be represented by the eigenfunctions $\boldsymbol {\phi }^{(i)}(\boldsymbol {x},f)$ as

(3.7)\begin{equation} \hat{\boldsymbol{q}}(\boldsymbol{x},f) = \sum_{i=1}^{\infty} \sqrt{\lambda^{(i)}(f)} \boldsymbol{\phi}^{(i)}(\boldsymbol{x},f). \end{equation}

In the case of SPOD, the weighted inner product in (3.1)–(3.2) will be a space–time integral

(3.8)\begin{equation} (\boldsymbol{q}_1, \boldsymbol{q}_2)_{{{\boldsymbol{\mathsf{W}}}}} = \int_{\boldsymbol{\varOmega}} \int_{-\infty}^{\infty} \boldsymbol{q}_2^{H}(\boldsymbol{x},t) {{\boldsymbol{\mathsf{W}}}}(\boldsymbol{x}) \boldsymbol{q}_1(\boldsymbol{x},t) \, {\rm d}\kern0.07em\boldsymbol{x}\, {\rm d} t. \end{equation}

We note that, at the same frequency $f$, different eigenvectors $\boldsymbol {\phi }^{(i)}(\boldsymbol {x},f), \boldsymbol {\phi }^{(j)}(\boldsymbol {x},f)$ are orthogonal under the spatial inner product (3.1) due to the symmetric positive definiteness of ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x},\boldsymbol {x}',f)$. But eigenvectors $\boldsymbol {\phi }^{(i)}(\boldsymbol {x},f_1), \boldsymbol {\phi }^{(i)}(\boldsymbol {x},f_2)$ at the same rank ($i$) associated with different frequencies are not necessarily orthogonal under the space-only inner product.

Our numerical implementation is similar to those in Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018) and Schmidt & Colonius (Reference Schmidt and Colonius2020). Data are sampled into blocks of sequenced snapshots (shown below is the $l$th block)

(3.9)\begin{equation} {{\boldsymbol{\mathsf{Q}}}}^{(l)}= [{{\boldsymbol{\mathsf{q}}}}_1^{(l)}, {{\boldsymbol{\mathsf{q}}}}_2^{(l)},\ldots, {{\boldsymbol{\mathsf{q}}}}_{N_{FFT}}^{(l)}] \in \mathbb{R}^{N_d \times N_{FFT}}, \end{equation}

where each column ${{\boldsymbol{\mathsf{q}}}}_i^{(l)}$ is one snapshot. The total number of snapshots in one block (ensemble) is $N_{FFT}$, where fast Fourier transform (FFT) will be performed over each block. The degree of freedom of one snapshot is $N_d = N_x \times N_y \times N_z \times N_{var}$, where $N_{var}$ is the dimension of the vector $\boldsymbol {q}(\boldsymbol {x},t)$. In this study, we apply SPOD on the vertical component of vorticity ($\omega _z$, hence $N_{var}=1$) in horizontal and vertical two-dimensional cross-sections of the flow (with either $N_z=1$ or $N_y=1$, respectively).

A discrete Fourier transform (DFT) is performed on each block ${{\boldsymbol{\mathsf{Q}}}}^{(l)}$ to yield

(3.10)\begin{equation} \hat{{{\boldsymbol{\mathsf{Q}}}}}^{(l)} = [\hat{{{\boldsymbol{\mathsf{q}}}}}_1^{(l)}, \hat{{{\boldsymbol{\mathsf{q}}}}}_2^{(l)},\ldots, \hat{{{\boldsymbol{\mathsf{q}}}}}_{N_{FFT}}^{(l)}] \in \mathbb{C}^{N_d \times N_{FFT}}. \end{equation}

Then, the Fourier modes are sorted according to frequency (labelled as the $k$th discrete frequency) to form

(3.11)\begin{equation} \hat{{{\boldsymbol{\mathsf{Q}}}}}_{k} = [\hat{{{\boldsymbol{\mathsf{q}}}}}_{k}^{(1)}, \hat{{{\boldsymbol{\mathsf{q}}}}}_{k}^{(2)},\ldots, \hat{{{\boldsymbol{\mathsf{q}}}}}_{k}^{(N_{blk})}] \in \mathbb{C}^{N_d \times N_{blk}}. \end{equation}

The sampled spectral density at the $k$th frequency is then ${{\boldsymbol{\mathsf{S}}}}_k = \hat {{{\boldsymbol{\mathsf{Q}}}}}_{k} \hat {{{\boldsymbol{\mathsf{Q}}}}}_{k}^{H} /(N_{blk}-1)$ and the discrete form of the eigenvalue problem (3.6) is

(3.12)\begin{equation} {{\boldsymbol{\mathsf{S}}}}_k {{\boldsymbol{\mathsf{W}}}} {\boldsymbol{\varPhi}}_k = {\boldsymbol{\varPhi}}_k {\boldsymbol{\varLambda}}_k, \end{equation}

with the weight matrix ${{\boldsymbol{\mathsf{W}}}}$ containing the weights of numerical quadrature at each grid point. In practice, (3.12) is typically solved with the method of snapshots of Sirovich (Reference Sirovich1987) by replacing $\boldsymbol { \varPhi }_k = \hat {{{\boldsymbol{\mathsf{Q}}}}}_{k} \boldsymbol {\varPsi }_k$ such that (3.12) becomes an equivalent eigenvalue problem

(3.13)\begin{equation} \frac{1}{N_{blk}-1} \hat{{{\boldsymbol{\mathsf{Q}}}}}_{k}^{H} {{\boldsymbol{\mathsf{W}}}} \hat{{{\boldsymbol{\mathsf{Q}}}}}_{k} {\boldsymbol \varPsi}_k = {\boldsymbol \varPsi}_k \boldsymbol{\varLambda}_k, \end{equation}

that has a much smaller dimension when $N_{blk} \ll N_d$ is true. Hence, the eigenmodes of ${{\boldsymbol{\mathsf{S}}}}_k$ are recovered as $\tilde {\boldsymbol { \varPhi }}_k = \hat {{{\boldsymbol{\mathsf{Q}}}}}_{k} \boldsymbol {\varPsi }_k \boldsymbol {\varLambda }_k^{-1/2}$ such that the eigenvalue decomposition is

(3.14)\begin{equation} {{\boldsymbol{\mathsf{S}}}}_k =\tilde{\boldsymbol{ \varPhi}}_k \boldsymbol{\varLambda}_k \tilde{\boldsymbol{ \varPhi}}_k^{H} = \sum_{i=1}^{N_{blk}} \lambda_{k}^{(i)} \tilde{\boldsymbol{\phi}}_k^{(i)} (\tilde{\boldsymbol{\phi}}_k^{(i)})^{H}. \end{equation}

The physical meaning of the spatial modes $\tilde {\boldsymbol { \varPhi }}_k(\boldsymbol {x})$ can be interpreted as either the eigenvector of the spectral density tensor ${{\boldsymbol{\mathsf{S}}}}_k$ or the left singular vector of the Fourier mode $\hat {{{\boldsymbol{\mathsf{q}}}}}_{k}$, at a discrete frequency $f_k$.

The SPOD method takes advantage of extracting spatial modes that evolve at a single frequency from a time-resolved database, as in table 1. It was applied to analyse stratified wakes by Nidhan et al. (Reference Nidhan, Chongsiripinyo, Schmidt and Sarkar2020) and Nidhan, Schmidt & Sarkar (Reference Nidhan, Schmidt and Sarkar2022), who showed that it can successfully extract the large-scale VS motions and the associated characteristic frequency. In oceanography applications, Zeiden et al. (Reference Zeiden, MacKinnon, Alford, Rudnick, Voet and Wijesekera2021) applied a similar approach called empirical orthogonal functions (EOFs) therein to the flow past an island and successfully separated the vortical modes and the tidal modes. But in their case, the EOFs are fit to three mooring points instead of the bulk of the flow, hence they are different from ours where the eigenmodes will be emphasised as global modes.

When converting the time series into Fourier modes in (3.10), Welch's method (Welch Reference Welch1967) is used to reduce the variance of the spectrum, with $N_{FFT} = 512$ snapshots in each block, and a Hamming window on each block to enforce periodicity. An overlap ratio of $50\,\%$ ($N_{ovlp}=256$) between two sequential blocks is chosen to offset the effect of low weights near the edges of the window. We end up with 13 blocks and the ensemble average will be taken over all blocks to obtain SPOD eigenspectra. The convergence of the method is checked via reducing the frequency resolution to $N_{FFT} = 256$, or reducing the total number of snapshots from $N_t=4000$ to 3000, and 2000. A high confidence level is obtained for the largest six eigenvalues as well as the sum of all eigenvalues, at each frequency. This follows the fact that for general eigenvalue-revealing algorithms, large eigenvalues converge faster. And it is noted that, in the present wakes, converging high-rank SPOD eigenvalues with much smaller magnitudes is still challenging even with $N_t=4000$ snapshots. In this paper, no particular analysis will be conducted for higher than the sixth eigenvalue at any frequency.

In this database, a constant maximal Courant–Friedrichs–Lewy (CFL) number is kept during the simulation, which results in an uneven (but almost uniform) time spacing of snapshots. To obtain uniformly spaced data for DFT, a piecewise cubic Hermite interpolation (PCHIP) is performed in time.

3.3. The SPOD eigenspectra and vortex shedding frequencies

Figure 4(a–c) shows the global vertical enstrophy spectra $S_{\omega _z \omega _z}(f)$ at different 2-D planes for BuInf, Bu25 and Bu1, respectively. The spectral density at a discrete frequency $f_k$ is computed by summing over all SPOD eigenvalues at this frequency

(3.15)\begin{equation} S_{\omega_z \omega_z}(f_k) = \sum_i^{N_{blk}} \lambda^{(i)}(f_k) = {\rm tr}({{\boldsymbol{\mathsf{S}}}}_k {{\boldsymbol{\mathsf{W}}}}) = \int_{\boldsymbol{\varOmega}} {{\boldsymbol{\mathsf{S}}}}(\boldsymbol{x},\boldsymbol{x},f_k)\,{\rm d}\kern0.07em\boldsymbol{x}, \end{equation}

and is the spectral density of the area-integrated squared $\omega _z$. It is independent of SPOD.

Figure 4. (a–c) Global power spectra $S_{\omega _z \omega _z}$ of $\omega_z$ as a function of Strouhal number ($St=f\kern 0.06em D/U_{\infty }$). Spectra are shown at four horizontal planes $z/h=0.12,0.25,0.50,0.75$, and the central vertical plane at $y=0$. The VS Strouhal number is marked as $St_{VS}$ as well as the VS harmonics. The values take $St_{VS}=0.264, 0.249$ and $0.266$ in (a–c), for cases BuInf, Bu25 and Bu1, respectively. (d) Decay of the largest eigenvalue at the VS frequency $f_{VS}$ and the harmonics (indexed by $n=f/f_{VS}$) for the horizontal plane $z/h=0.25$ (results are similar for other locations). Colour codes same as in (a–c).

For all three cases and all planes examined, the spectra display strong harmonic spikes, with the strongest peak defined as the VS frequency ($\,f_{VS}$) in each plane. The VS frequency is independent of the vertical location of the four horizontal planes, and is also the same in the central vertical plane, even though performing SPOD on separate planes allows the freedom of selecting different frequencies. This indicates that, for each case, $f_{VS}$ is a global constant (for the heights examined), and the VS modes are 3-D global modes. The VS Strouhal number is $St_{VS} = f_{VS} D/U_{\infty } = 0.264,0.249,0.266$ for BuInf, Bu25 and Bu1, respectively. It is noted that, since the characteristic frequency of the global mode ($\,f_{VS}$) does not depend on the height or the local hill diameter, normalising it with a different length scale than the hill base diameter $D$ would just make $St_{VS}$ different by a scalar multiple. However, as will be discussed later, the numerical values of $St_{VS}$ using $D$ correspond well to that of VS from a circular cylinder. Also, since $D$ determines the size of the largest scales in the flow, it is a natural choice for the normalising scale. In the eigenspectra, the successive peaks are harmonics of the VS frequency at $2St_{VS}, 3St_{VS}$ and so on.

Perfect et al. (Reference Perfect, Kumar and Riley2018) found that whether the VS frequency is a global constant or is controlled by the local hill diameter depends on a non-dimensional parameter: the Burger number ($Bu$). The Burger number characterises the relative importance of two counteracting mechanisms for vertical coupling: rotation and stratification. When $Bu$ is small (rotation is strong), the bulk of the flow adjusts to be around the axis of rotation quickly, and geostrophic balance is established, where the vertical variation is minimised. As $Bu$ is increased, the vertical intercommunication is progressively weakened by stratification. Perfect et al. (Reference Perfect, Kumar and Riley2018) use the diameter at half-height to be the characteristic horizontal scale so that their values of Rossby number ($Ro^*$) and Burger number ($Bu^*$) are related to our values by $Ro^* = 2 Ro$ and $Bu^* = 4 Bu$. They performed a number of simulations and found that ${Bu}^{*}_{cri}=5.5$, equivalently ${Bu}_{cri} = 1.4$, is a regime-separation criterion below which the rotation is strong enough to couple different layers to form vertically aligned vortices. Also, when ${Bu^*} > 12$, equivalently ${Bu} > 3$, they found stratification to be more prominent so that vortices are shed at different layers relatively independently. Turning to the present results, in each of the cases at $Fr=0.15$ whose $Ro$ span unity to infinity, the VS frequency is independent of height, indicating a potential disagreement between Perfect et al. (Reference Perfect, Kumar and Riley2018) and our results. On the other hand, for Froude numbers similar to $Fr=0.15$ in Perfect et al. (Reference Perfect, Kumar and Riley2018), almost all the cases (see figure 4, therein) were labelled as ‘vertically coupled shedding’ and had strong rotation with $Ro$ between 0.025 and 0.25. Only their weakest rotation case with $Ro = 0.25$ was labelled as ‘vertically decoupled shedding’. It is also worth noting that, in the ROMS simulations, vertical motions and pressure correlations are quite approximate (especially in the near wake where VS is accompanied by small-scale turbulence) since the momentum equation in that direction is reduced to a hydrostatic balance, and the pressure field might play an important role in coupling VS.

Nevertheless, the fact that the modes extracted by SPOD are global modes, and they evolve at the same frequency, implies that the large-scale vortices at $Fr = 0.15$ are horizontally and vertically coherent, as opposed to the asymptotic limit of vertically uncoupled layered dynamics in strongly stratified flows. Our findings indicate that the stratification of $Fr=0.15$ is not strong enough to vertically decouple the vortex dynamics in the wake of the conical seamount, regardless of the presence of rotation, and inclusion of $Fr$ is necessary in addition to $Bu$. The pure Burger number criterion has the limitation that, for instance, with no rotation and small stratification, $Bu$ will be far larger than $Bu_{cri}$ but the VS frequency can still be a global constant. The $Fr$ that marks the transition from vertically coupled to decoupled VS in non-rotating hill wakes is subject to future research.

Boyer et al. (Reference Boyer, Davies, Holland, Biolley and Honji1987) studied experimentally the wake behind a conical obstacle in linearly stratified rotating flows. Their parameters spanned $0.08< Fr<0.28$, $0.06< Ro<0.4$ and for three Reynolds numbers $Re_D=380,760,1140$. Based on the measurement on single horizontal cross-sections, they found the VS Strouhal number to be only a weak function of both $Re$ and $Ro$. The robustness of the VS frequency to rotation strength is also observed numerically in this work, in a similar $Fr$–$Ro$ regime but at a turbulent Reynolds number $Re_D=10\,000$. Even though their Strouhal numbers are measured as the vortex advection velocity divided by the mean separation of two same sign vortices, and were in the range $0.20< St<0.35$, our VS Strouhal numbers still present a good quantitative agreement with theirs. Moreover, we interpret the VS frequency revealed by SPOD as the characteristic frequency of the most energetic global mode, which also agrees with single-point frequency measurements in the intermediate wake ($x/D>3$, not shown).

The values of $St_{VS}$ for all three cases are close to $St=0.2665$, which is the $Re \rightarrow \infty$ asymptote of the $St$–$Re$ relationship in low Reynolds number 2-D cylinder wakes proposed by Williamson & Brown (Reference Williamson and Brown1998). Note that their relation $St= 0.2665\unicode{x2013}1.018\sqrt {Re}$ is given for $50< Re_D<180$, which is before the transition to 3-D wakes. This transition happens at around $Re=188.5$ according to a global Floquet instability of the periodic wake (Barkley & Henderson Reference Barkley and Henderson1996). As a result, the $St$–$Re$ relationship experiences a discontinuity as a sudden jump of $St$ during this transition (Fey, König & Eckelmann Reference Fey, König and Eckelmann1998; Williamson & Brown Reference Williamson and Brown1998), and the asymptote of $St=0.2665$ is not reached in 3-D cylinder wakes. Fey et al. (Reference Fey, König and Eckelmann1998) showed that the maximum VS Strouhal number in a cylinder wake of about $St=0.21$ is reached right before the onset of the Kelvin–Helmholtz instability in the shear layer. We interpret the observed values of VS frequencies in our hill wakes as a saturation of the Q2-D VS, which will not be observed in 3-D cylinder wakes at this Reynolds number and above. This interpretation is consistent with the finding in Boyer et al. (Reference Boyer, Davies, Holland, Biolley and Honji1987) that the robustness of $f_{VS}$ to rotation rate is not affected by the Reynolds number, in their low-to-moderate Reynolds number experiments.

The enstrophy distribution among different eigenvalues is an important measure of the complexity of the system. Figure 5 shows the sum of all eigenvalues, and also the first to the sixth eigenvalues (from dark to light) according to their absolute value. For frequencies at or close to the VS frequency and its harmonics, the first eigenvalue accounts for most of the enstrophy, and is an order of magnitude larger than the second eigenvalue. However, for larger frequencies ($St>2$), the dominance of the leading eigenvalues is lost, as the degree of freedom for small-scale motions is increased. The significance of low rankness in the spectra is therefore twofold. There are strong harmonic spikes in the SPOD eigenspectra in figure 4(a–c), implying that a great portion of enstrophy is contained in the large-scale VS motions. For the VS shedding frequency and its harmonics, the first two eigenvalues contribute the most of the enstrophy.

Figure 5. The SPOD eigenspectra: (a) BuInf, (b) Bu25 and (c) Bu1. Horizontal plane at $z/h=0.25$ is shown (results are similar for other locations). From top to bottom are the summation of all eigenvalues (darkest, spectral density as in figure 4a–c), and the first to the sixth eigenvalues (from dark to light: $i=1,2,3,4,5,6$), as a function of Strouhal number ($St=f\kern 0.06em D/U_{\infty }$). The difference between the first and the second eigenvalues is filled with colour.

Moreover, the enstrophy distribution among the harmonics is shown in figure 4(d). The decay of the first eigenvalue $\lambda ^{(1)}$ as a function of the integer harmonic index $n=f/f_{VS}$ is plotted. Similar to figure 5, one horizontal plane $z/h=0.25$ is chosen, but the results are qualitatively similar for the other three selected horizontal planes. For all three Burger numbers, $\lambda ^{(1)}$ decays approximately as a power progression as $n^{-3}$, which is slower than the geometric decay of distinct POD eigenvalues in 2-D cylinder wakes (Noack et al. Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003). Nevertheless, higher harmonics contain a progressively smaller amount of enstrophy. Additionally, by comparing the absolute value of $\lambda ^{(1)}$ in figure 4(d) and $\sum _i \lambda ^{(i)}$ in figure 5, it is generally true that vertical enstrophy is increased as the rotation rate increases. This is an effect of system rotation that promotes vortical motion around the rotation axis as will be discussed in more detail in § 5.

In all, the large scales of the vortical motion can be well represented by a few characteristic frequencies and the associated leading modes. Such low-rank behaviour suggests the possibility of reduced-order modelling of stratified hill wakes even at large Reynolds numbers.

3.4. The SPOD eigenfunctions and large-scale global modes

Apart from the temporal characteristics uncovered by SPOD eigenspectra, the SPOD eigenmodes represent energetic flow structures that are coherent in space and time and, therefore, dynamically important. The dominant modes in each case are found to be the VS modes, which are three-dimensionally coherent and characterised by the VS frequency ($\,f_{VS}$) and its superhamonics.

The real parts of the leading SPOD modes (corresponding to the largest eigenvalue at each frequency) are plotted in figures 6 and 7, for two frequencies $f_{VS}$ and $2f_{VS}$ and on various 2-D planes. The eigenfunctions are free up to a scalar multiple and are normalised by the largest modulus of the spatial mode in the same plane. The interpretation of the magnitude of the eigenmodes (or the darkness of the colours) as the intensity of the structures is only meaningful when the comparison is made within the same plane.

Figure 6. Leading SPOD eigenmodes corresponding to $f_{VS}$ (a,c) and $2f_{VS}$ (b,d), for cases BuInf (a,b) and Bu25 (c,d). The plotted quantity is the real part of SPOD eigenfunction on each plane, normalised by the maximum value in the plane. In each figure, the top plane is the vertical plane, $y=0$, and the lower planes are, from the second row to bottom, horizontal planes at $z/h=0.75, 0.50, 0.12$.

Figure 7. For case Bu1. Caption same as in figure 6.

For case BuInf (no rotation), the VS mode in figure 6(a) is reminiscent of the marginally stable global modes found in the linear stability analysis of low Reynolds number cylinder wakes (Barkley & Henderson Reference Barkley and Henderson1996; Mittal Reference Mittal2008). However, the phases of the global mode are different in different planes, with the phases of higher planes leading. This is clearer in the vertical plane ($\kern 1.5pt y=0$), which shows the tilting and elongation of the structures as an effect of stratification. Note that the vertical coordinate is normalised by the hill height $h$ which is approximately 0.3 times the base diameter; in un-normalised coordinates, the angle of the slanted ‘tongues’ is very shallow. The tilting angle from the horizontal of the structures in the top panel of figure 6(a) ranges from approximately $8^{\circ }$ close to the bottom wall to $2^{\circ }$ near the top, with an average of approximately $4^{\circ }$, which is almost constant during the evolution. Although the VS structures are not vertical, they still evolve cohesively at $f_{VS}$ and experience little change during the advection. In the later § 4, the lack of change of the tilting angle will be shown to be a result of a roughly constant advection velocity of the vortices. A similar slanted mode as in figure 6(a) was found in a stratified sphere wake by Chongsiripinyo, Pal & Sarkar (Reference Chongsiripinyo, Pal and Sarkar2017) (referred as ‘surfboards’ therein), but whether it is a common feature of stratified bluff-body wakes is subject to future research.

Moreover, the standing-wave-like spatial modes of $f_{VS}$ in figure 6(a) are symmetric about $y=0$, representing the advection of the perturbation vorticity $\omega '_z=\omega _z - \langle \omega _z \rangle$. The projection of $\omega '_z$ on the $f_{VS}$ mode has a period of $T_{VS}$ and the streamwise wavelength of the mode is interpreted as the average spacing of the Kármán vortices. The appearance of the symmetric $f_{VS}$ mode on the antisymmetric mean flow (see figure 3d), in analogy to low Reynolds number 2-D cylinder wakes (Kumar & Mittal Reference Kumar and Mittal2006), can be viewed as an accompaniment to the symmetry-breaking bifurcation from a steady wake to periodic shedding.

The modes of $2 f_{VS}$ in figure 6(b) are antisymmetric, with a wavelength approximately half of that in the $f_{VS}$ mode. A result of this reflection antisymmetry is that the eigenspectrum at the vertical plane (black dashed line in figure 4(a) does not exhibit a peak at $2f_{VS}$, unlike the spectra at horizontal planes. The antisymmetry implies a zero magnitude of the $2f_{VS}$ mode at the centreline ($\kern 1.5pt y=0)$, hence the top row in figure 4(b) does not imply any dynamical importance.

The imaginary parts of the eigenmodes are not shown, which are different from the real parts by a streamwise phase shift of ${\rm \pi} /2$, and are therefore representing the same VS dynamics. The phase shift between two standing-wave-like eigenmodes with the same streamwise wavelength (such as the real and imaginary parts of the $f_{VS}$ mode here) is essential for them to accommodate one travelling-wave-like structure (such as the advecting VS motion) as noted by Rempfer & Fasel (Reference Rempfer and Fasel1994).

For cases Bu25 and Bu1, the reflectional symmetry is broken by the Coriolis force and, unlike BuInf, the peak at $2f_{VS}$ is observed in the eigenspectra at the vertical plane (black dashed lines). The rightward ($-y$) shift of the vortex wake is consistent with the direction of the Ekman veering near the bottom by the unbalanced mean pressure gradient. This asymmetry can also be found in the SPOD modes in figures 6(c,d) and 7(a,b), where the cyclones (on the left side of the bulk of the wake, looking from above) are preferred and amplified, as compared with the anticyclones. However, in the vertical mode (first row in figure 6c,d), slanted structures are still statistically significant, even though figure 2(b) shows that, for case Bu25, some anticyclones have already become columnar. This suggests that the stratification in Bu25 is still dominating, and rotation has not been able to modify the spatial organisation of the structures. However, the presence of rotation significantly alters the structure inside the anticyclones, which will be analysed in § 5.

For case Bu1 (strong rotation), the global modes in figure 7(a,b) show excellent vertical alignment, agreeing with the shed ‘columns’ in figure 2(c). In the vertical plane (see the first row in figure 7a), both cyclonic and anticyclonic structures extend slightly higher than the obstacle peak, whereas the structures in cases Bu25 and BuInf are limited to below the obstacle peak, signifying that the obstacle's range of influence is vertically increased by increasing rotation. We note that those ‘columns’ are different from the typical Taylor column in rotating flow over obstacles. A Taylor column has an infinite height on the top of a finite object that generates it, whereas the columnar global mode in case Bu1 has a finite height. In figure 7(a,b) it can be seen that the parts of the global mode near the hill ($0< x/D<1$) are smaller and slanted, and they eventually become organised into tall ‘columns’ during their advection downstream as rotation is experienced, instead of being columnar at the generation.

Finally, the relation between the global VS modes found in the present stratified wakes and the global modes in 2-D Kármán wakes is worth further exploration. The VS modes observed in each 2-D horizontal plane of the non-rotating case (BuInf, see figures 6–7) are visually similar to those observed in 2-D Kármán wakes, and the associated VS Strouhal number ($St_{VS}$) is close to the $Re\rightarrow \infty$ limit of the $St$–$Re$ relation of 2-D cylinder wakes. It is worth pointing out that strong stratification is crucial for the similarity between the present high Reynolds number wakes and 2-D Kármán wakes. Alternating shedding vortices are not seen at a Reynolds number similar to $Re_D=10\,000$ in unstratified cylinder wakes (Fey et al. Reference Fey, König and Eckelmann1998) and stratified sphere wakes at $Fr>O(0.5)$ (Pal et al. Reference Pal, Sarkar, Posa and Balaras2016), where the attached shear layer will have lost its stability.

On the other hand, the slanted, forward-tilting vertical structures seen in the $y=0$ plane of BuInf case are suggestive of vertical coherence of horizontal vortical motion, which is absent for 2-D flows. As a part of the global mode, the VS frequency and spatial mode in each horizontal plane depends rather on the hill base diameter than on the local hill diameter at the same height – another difference from typical 2-D Kármán wakes.

The change of the titling angle of the global mode as $Ro$ decreases until the vortex structures are fully upright at $Ro=0.15$ adds an additional aspect – the influence of rotation. Even though the vertical alignment of the structures is varied, the dominant modes are still the VS modes with the shedding frequencies being similar in each case, indicating that the VS motions at present $Fr=O(0.1)$ are relatively robust to rotation. Slanted vortical structures as in BuInf were also observed in stratified non-rotating wakes past a sphere by Pal et al. (Reference Pal, Sarkar, Posa and Balaras2016) and Chongsiripinyo et al. (Reference Chongsiripinyo, Pal and Sarkar2017). Whether these structures are common in strongly stratified wakes past a body of revolution with vertically varying diameter and how these structures are modified by further increase in stratification deserves future study.

5.1. Probability distribution of vorticity

The p.d.f.s of $|\omega _z|$ for positive vorticity (solid line) and negative vorticity (dashed line) are shown separately and on the planes $z/h=0.12$ (figure 11a) and $0.50$ (figure 11b). Each p.d.f. is normalised such that the area under each line is equal. In the p.d.f., $|\omega _z|$ is normalised with convective units ($D$ and $U_{\infty }$) for consistency between rotating and non-rotating cases (where $f_{c}$ is zero in BuInf), but the Coriolis frequency ($|\,f_{c}|$) will also be used for normalisation in cases Bu25 and Bu1. For case BuInf, symmetry is achieved for vorticity of all magnitudes, as expected. Comparing cases Bu25 and Bu1, there is an increasing asymmetry between the p.d.f. of positive and negative vorticities, as the rotation strength is increased.

Figure 11. The p.d.f. of $|\omega _z|$ at (a) $z/h=0.12$, and (b) $z/h=0.50$. Solid lines indicate positive vorticity (AVs), and dashed lines indicate negative vorticity (CVs). Vorticity is presented in convective units and the dotted vertical reference lines represent the magnitudes of the non-dimensionalised Coriolis frequency $|\,f_{c}| D/U_{\infty }$, which equals $1/Ro$ and takes the value of $4/3$ for case Bu25 and $20/3$ for case Bu1, respectively.

Consider the Bu1 case (blue lines). In plane $z/h=0.12$ shown in figure 11(a), there is a local peak for cyclonic vorticity (negative $\omega _z$) and one for anticyclonic vorticity (positive $\omega _z$), both being slightly larger than the system rotation rate and close to $1.1|\,f_{c}|$. In plane $z/h=0.50$ shown in figure 11(b), there is a local peak only for anticyclonic (positive) vorticity at $0.7 |\,f_{c}|$.

For case Bu25, the p.d.f. of cyclonic vorticity at $z/h=0.12$ shown in figure 11(a) also has a peak above $|\,f_{c}|$ (with peak relative vorticity $3.1 |\,f_{c}|$, which is significantly greater that in case Bu1). On the other hand, anticyclonic vorticity does not show any observable peak near $|\,f_{c}|$.

The local peaks in p.d.f.s are interpreted as the values that are more commonly found in the flow compared with their neighbours, instead of the most intense ones. In the next section, we will show the existence and persistence of intense anticyclones ($\omega _z > |\,f_{c}|$) in both Bu1 and Bu25 and discuss their stability.

5.2. Vorticity conditioned to individually tracked vortices

The vorticity distribution inside wake vortices is essential to the understanding of their kinematics, the idealisation and modelling of such vortical wakes and the role of vortex stability. The previous section on the vorticity p.d.f., while giving an overall statistical view of cyclonic/anticyclonic vorticity, does not reveal the properties of individual coherent wake vortices – the focus of this section. The identification of vortex centres in § 4 is leveraged to quantify the vorticity conditioned to the identified vortex centres and thus reveal the flow inside and around the vortices. As elaborated below, rotation significantly impacts the downstream evolution of the vortex-conditioned distribution of $\omega _z$ and, notably, the difference between AV and CV is larger for $Bu = 25$ than for $Bu =1$.

Figure 12 shows profiles of the average of $\omega _z(\tilde {x})$, conditioned to instantaneous individual vortex centres. Here, $(\tilde {x},\tilde {y})=(x-x_c,y-y_c)$ is a new set of horizontal coordinates fixed to individual vortex centres. Since wakes are spatially developing, the results are presented at various values of $x_c$ to diagnose the downstream evolution of vortex-conditioned properties. Vortices with centres less than $2D$ apart are assumed to possess similar properties and are grouped for a regional average. For example, the group located at $x_c/D=4$ represents vortex centres that fall in the section of $3< x/D<5$ and so forth. Each vortex centre in the same group is shifted to $\tilde {x} = \tilde {y}=0$ before the statistics are gathered. For each group, more than 2000 vortices are available for the ensemble average.

Figure 12. Conditionally averaged vorticity distribution around vortex centres, $\omega _zD/U_{\infty }$ for positive vortices (AVs, solid) and $-\omega _zD/U_{\infty }$ for negative vortices (CVs, dashed). (a,b) BuInf, (c,d) Bu25 and (e, f) Bu1. Horizontal planes shown are at $z/h=0.12$ (a,c,e) in the left column and at $z/h=0.50$ (b,d, f) in the right column. The horizontal dotted lines in (c–f) indicate the case-dependent absolute magnitude of the Coriolis frequency ($|\,f_{c}|$). In panels (c–f), the strongest CVs have peak magnitudes of $\omega _z/|\,f_{c}|=5.1, 3.0, 2.2, 2.2$, and the strongest AVs have peaks $\omega _z/|\,f_{c}|=-2.4, -1.2, -1.3, -0.8$, respectively.

For case BuInf shown in figure 12(a,b), vorticity profiles are Gaussian-like except close to the edges. The peak magnitude decays and the width grows as a result of diffusion. Otherwise, no significant change is present when $x$ increases. In figure 12(a) at $z/h=0.12$, the magnitude of $\omega _z$ near $x/D=\pm 1$ is close to zero – the asymptotic far-field condition for isolated vortices. However, in figure 12(b) at $z/h=0.50$, the vorticity can change sign when $\tilde {x}/D$ varies between $0$ and $-1$, becoming substantially negative for the group $x_c/D=4$ but less so further downstream. A similar sign change is not observed when $\tilde {x}$ varies between $0$ and $1$. As seen in the middle and bottom rows of figure 3(a), wake vortices of both signs at $z/h=0.50$ are more spatially diffuse than at $z/h=0.12$ and are almost side by side in the region around $x/D=4$, making possible the sign change of vorticity.

For case Bu25 shown in figure 12(c,d), CVs and AVs are substantially different. The CVs are stronger and more compact, while the AVs are weaker and wider. The latter is likely because of the cyclo-geostrophic instability associated with the AVs upstream before the conditional statistics are gathered. In figure 12(c), the vorticity distribution for $x_c/D=4$ displays short-wavelength wiggles and represents active instability of the AVs as can also be seen in the bottom row of figure 3(b). We emphasise that the short-wavelength wiggles are in a profile that is obtained by averaging over an ensemble of approximately 2000 members and are, thus, statistically significant. This is later confirmed in § 5.3, where the generalised Rayleigh discriminant on the left side of the aforementioned AVs is shown to lie beyond the stability limit.

Figure 12(d) shows differentiated behaviour between CVs and AVs when the vortex edge is approached. The conditionally averaged vorticity in AVs (positive $\omega _z$) tends to change sign toward the right end, while that around CVs (negative $\omega _z$) changes sign toward the left end. This indicates a preferred configuration of vortex dipoles with a positive vortex on the left and a negative vortex on the right, as can also be seen in the middle row of figure 3(b). At the same time, the dipoles are quite asymmetric; the positive $\omega _z$ AV is weaker and more spatially diffuse than its partner, suggesting that the weaker AVs are more susceptible to the induced motion of the CVs. It is in case Bu25 that the asymmetry between CVs and AVs is most pronounced. As shown in figure 12(d), the ratio of peak magnitudes between cyclones and anticyclones is approximately four on average, and increases as downstream distance is increased. In combination with the fact that the absolute magnitude of the anticyclones is weak (order $D/U_{\infty }$), they are more easily influenced by adjacent cyclones during mutual interaction.

Strong rotation favours the formation of coherent vortices. For case Bu1 shown in figure 12(e, f), both CVs and AVs are strongest in units of $U_{\infty }/D$ in all three cases. In the plane at $z/h=0.12$, both CVs and AVs have magnitude greater than $|\,f_{c}|$ while only CVs exceed $|\,f_{c}|$ in the plane at $z/h=0.50$. Moreover, both CVs and AVs undergo little change in terms of vorticity magnitude and distribution during their advection – a key difference from the other two Burger numbers. This is in agreement with the mean flow characteristics to be discussed in § 6 that the streamwise change in the momentum wake of Bu1 is slower than in the other two cases.

In terms of the strength of the AVs, the vorticity magnitude in the vortex core exceeds $|\,f_{c}|$ at all streamwise locations in the plane at $z/h=0.12$, in both cases Bu25 and Bu1. Thus, the absolute vorticity stability criterion is not conclusive to the behaviour of AVs since, contrary to that stability criterion, AVs with $|\omega _z| > f$ are found to advect in a stable manner.

Table 3 summarises the properties of CVs and AVs, in terms of peak intensity and vortex size. The peak intensity $\omega _{z,p}$ is the maximum of $\omega _z$ of each curve in figure 12, and the vortex size $l_{0.05}$ is the horizontal distance within which the magnitude of $\omega _z$ is greater than $0.05 \omega _{z,p}$. It can be seen that, at the same spatial location, vortex intensity is generally higher in Bu1 in convective units ($U_{\infty }/D$). Moreover, the sizes of both CVs and AVs are consistently smaller in case Bu1. In combination with greater peak vorticity, velocity gradients are much larger in vortices in case Bu1 than in Bu25. In terms of vorticity in units of $f_{c}$, Bu25 stands out with the largest magnitude of $\omega _z/|\,f_{c}|$.

Table 3. Summary of the properties of CVs and AVs, at two horizontal planes $z/h=0.12, 0.50$ and two streamwise locations of vortex centres $x_c/D=4, 10$. The peak vertical vorticity $\omega _{z,p}$ is given in convective ($U_{\infty }/D$) as well as rotation ($\,f_{c}$) units. Vortex sizes are characterised by $l_{0.05}$, which is the diameter at which the intensity of $\omega _z$ decays to $0.05\omega _{z,p}$.

5.3. Stability of anticyclones

Vortex stability is an important question since it is related to the generation of turbulence and small-scale motions. Here, we assess the ability of various criteria in the literature to identify the stability characteristics of the advecting wake vortices of the present simulations. It will be shown that more recent criteria that account for stratification and viscous dissipation constitute a significant improvement over earlier attempts.

The study of the centrifugal (inertial) instability of swirling flows can be traced back to Rayleigh (Reference Rayleigh1917), who showed that the flow could become unstable if the squared angular momentum decreases radially. The Rayleigh criterion is a necessary and sufficient condition for the instability of inviscid columnar vortices subject to 3-D axisymmetric perturbations (Drazin Reference Drazin2002). In a rotating frame, the Rayleigh criterion for inertial instability is equivalent to the existence of a region with a negative product of Coriolis frequency and absolute vorticity (Holton Reference Holton1972)

(5.1)\begin{equation} {f_{c}}({\omega_z (r) + f_{c} }) < 0, \end{equation}

or, in terms of the local Rossby number (note $f_{c}$ can take either sign),

(5.2)\begin{equation} {Ro_{\rm L} (r)} = \frac{\omega_z (r) }{f_{c}} <{-}1 . \end{equation}

The criterion for inertial instability, (5.1), is widely used and implies that anticyclones with vorticity magnitude exceeding $|\,f_c|$ are unlikely. However, it assumes a sheared parallel flow as the base state (Holton Reference Holton1972). Kloosterziel & Van Heijst (Reference Kloosterziel and Van Heijst1991) found in their experiment that CVs could be unstable too, contrary to (5.2). With the inclusion of the centrifugal term, a generalised Rayleigh discriminant (Kloosterziel & Van Heijst Reference Kloosterziel and Van Heijst1991; Mutabazi, Normand & Wesfreid Reference Mutabazi, Normand and Wesfreid1992) for centrifugal instability in rotating vortical flows was established as

(5.3)\begin{equation} \chi(r) = [\omega_z(r)+f_{c}]\left[2\frac{u_{\theta}(r)}{r}+f_{c}\right] < 0, \end{equation}

where $\omega _z+f_{c}$ is interpreted as the absolute vorticity and $u_{\theta } + f_{c} r/2$ as the absolute velocity. It implies instability when the absolute velocity and absolute vorticity are of opposite signs. Equation (5.3) assumes a specific form (axisymmetric) of perturbations. Nevertheless, (5.3) is good enough in many circumstances since axisymmetric perturbations are, in general, more unstable than non-axisymmetric ones (Billant & Gallaire Reference Billant and Gallaire2005).

While the former two criteria (5.1) and (5.3) concern the stability to 2-D perturbations, the inclusion of effects in the third dimension, such as stable stratification, finite dissipation and the vertical variation of the vortex base flow, is necessary for actual geophysical flows and will complicate the determination of stability. With stratification that suppresses small wavenumbers and finite vertical dissipation that suppresses large wavenumbers, the range of unstable vertical wavenumbers is reduced and (5.3) overestimates the unstable region (Lazar, Stegner & Heifetz Reference Lazar, Stegner and Heifetz2013b). To reduce the predicted unstable region, Lazar et al. (Reference Lazar, Stegner and Heifetz2013b) proposed a new criterion for the marginally stable Burger number

(5.4)\begin{equation} \sqrt{Bu_{v}} = \left(\frac{3}{8|a_0|}\right)^{3/2} \frac{1}{\sqrt{Ek}} \frac{(|2Ro_{v}+1|)^{7/4}}{|Ro_{v}|}, \end{equation}

where

(5.5a–c)\begin{equation} Ro_{v}= \frac{V_{max}}{f_{c} r_{max}};\quad Bu_{v} = \left(\frac{Nh}{f_{c} r_{max}}\right)^2;\quad Ek = \frac{\nu}{{|\,f_{c}|} h^2} ,\end{equation}

are the vortex Rossby number, the vortex Burger number and the vertical Ekman number. Here, $V_{max}$ and $r_{max}$ refer to the peak magnitude of the azimuthal velocity and its location, respectively. Positive and negative $Ro_{v}$ represent cyclones and anticyclones, respectively. The constant $a_0=-2.338$ is the first zero of the Airy function. In the parameter space of $Bu_{v}$–$Ro_{v}$, (5.4) at each constant $Ek$ corresponds to a stability curve that separates regimes that are stable and unstable to axisymmetric perturbations.

More recently, Yim, Stegner & Billant (Reference Yim, Stegner and Billant2019) pointed out that the most unstable azimuthal wavenumber of the centrifugal mode is not necessarily $m=0$ (axisymmetric), but depends on $Bu_{v}$ and will have an impact on the determination of stability. Accordingly, they suggested the use of the stability curve given as

(5.6)\begin{equation} \sqrt{Bu_{v}} = \frac{0.23}{\sqrt{Ek}} \frac{(Ro_{v}+0.3)^2}{\sqrt{|Ro_{v}|}}, \end{equation}

which considered the dependence of the most unstable azimuthal wavenumbers on ${Bu_{v}}$. It was also found that the stability of AVs is insensitive to the vertical variation of the initially axisymmetric vortex base flow. In this section, we will compare the absolute vorticity criterion (5.1), the generalised Rayleigh discriminant (5.3) as well as the new criteria (5.4) and (5.6), and assess their ability to predict the stability of the wake vortices in the simulations.

As was shown in figure 12(c,e), anticyclones at $z/h=0.12$ are observed to possess a large region with $|\omega _z|>|\,f_{c}|$, although they advect as stable vortices for a significant range of downstream evolution distance. Hence, the absolute vorticity condition (5.1) is not sufficient for instability.

Figure 13 shows the conditionally averaged generalised Rayleigh discriminant $\chi (\tilde {x})$ in (5.3) as a function of streamwise distance from the vortex centre. It can be seen that the unstable region is reduced significantly compared with (5.1). In case Bu25 and at $z/h=0.12$ (figure 13a), at roughly $x/D=4$ (the darkest green line), there is a small region, located near the left edges of the AVs, that is found to be unstable, but otherwise, all regions have at least marginally stable $\chi$. Furthermore, the vortices tend to evolve to a more stable state. For Bu1, the AVs in the $z/h=0.50$ plane (figure 13b) have stable discriminant $\chi$, while the peripheries of the AVs in the plane $z/h=0.12$ have marginally unstable $\chi$ which does not experience noticeable change during the advection. The statistics of $\chi$ of the AVs in the plane at $z/h=0.25$ are similar to those at $z/h=0.12$ and are not shown. Consistent with the instantaneous snapshots in figure 3 and the statistics in figure 12, AV instability is only observed for Bu25 at $z/h=0.12$, and is captured properly by the generalised Rayleigh discriminant (5.3). On the other hand, in case Bu1 at $z/h=0.12$ (figure 13c), where the discriminant is marginally unstable at the periphery of the vortices, no actual change to the vorticity profile of the AVs is observed (see figure 12e). Hence, a sufficient condition for stability requires other considerations, e.g. stratification and dissipation (Lazar et al. Reference Lazar, Stegner and Heifetz2013b; Yim et al. Reference Yim, Stegner and Billant2019).

Figure 13. Conditionally averaged generalised Rayleigh discriminant $\chi /f_{c}^2$ for AVs in cases Bu25 (a,b) and Bu1 (c,d). Horizontal planes shown are at $z/h=0.12$ (a,c), and at $z/h=0.50$ (b,d). Streamwise locations, $x/D = 4, 6, 8$ and 10 are shown in dark to light colours. The dashed line in each figure marks the stability criterion of $\chi = 0$.

Prior to applying (5.4) and (5.6), the vortex sizes and shapes in terms of $V_{max}$ and $r_{max}$ are required. The radial direction is substituted by the streamwise ($x$) direction and the azimuthal velocity component by the spanwise velocity ($v$). The peak velocity $V_{max}$ is defined as the maximum azimuthal (herein transverse) velocity and the corresponding peak location $r_{max}$ is interpreted as vortex radius. It is noted that for both criteria (5.4) and (5.6) applied here, 2-D vortex profiles subject to 3-D perturbations are assumed, and the vertical variability of the vortex structure is not considered.

Figure 14(a) shows $V_{max}$ and $r_{max}$ for Bu 25 and Bu1 at various streamwise locations. It can be seen that, for both Bu1 and Bu25, the vortex radii agree reasonably well with the radial location of the least stable $\chi$ (figure 13), consistent with theoretical analysis and experimental observation that the edges of vortices are the least stable regions (Kloosterziel & Van Heijst Reference Kloosterziel and Van Heijst1991; Carnevale et al. Reference Carnevale, Briscolin, Kloosterziel and Vallis1997; Lazar et al. Reference Lazar, Stegner and Heifetz2013b; Yim, Billant & Ménesguen Reference Yim, Billant and Ménesguen2016). Moreover, AVs in Bu25 have a much larger radii as well as variability during the evolution, compared with Bu1. The AVs in Bu1, which have greater $V_{max}$ (over twice stronger than Bu25) and smaller vortex radii, have the largest average vorticity.

Figure 14. Vortex properties: (a) conditionally averaged maximum azimuthal velocity $V_{max}$ and corresponding location $r_{max}$ and (b) square root of the vortex Burger number $\sqrt {Bu_{v}}$ and the vortex Rossby number $Ro_{v}$. Bu25 is shown in green and Bu1 in blue. Colours from dark to light indicate planes at $z/h=0.12, 0.25, 0.50$. Each circle marks one of the locations from $x/D=4,6,8,10$. The circles are connected by lines following the order of $x$-locations and the last location ($x/D=10$) is marked with the largest circle. In (b), dashed lines denote the stability curves (5.4) by Lazar et al. (Reference Lazar, Stegner and Heifetz2013b) and dotted lines (5.6) by Yim et al. (Reference Yim, Stegner and Billant2019) for cases Bu1 (blue) and Bu25 (green). The Ekman numbers are $Ek=|Ro|/Re_D (D/h)^2 =$ $8.33 \times 10^{-4}$ and $1.67 \times 10^{-4}$ for cases Bu25 and Bu1, respectively. The left side of the neutral stability curves is stable.

Figure 14(b) shows the evolution of AVs, on average, in the $Bu_{v}$–$Ro_{v}$ parameter space. The AVs in both Bu25 and Bu1 are characterised by vortex Rossby numbers of $O(0.5)$. The stability curves (5.4) and (5.6) are also plotted for cases Bu1 and Bu25. The left side of a stability curve is the stable region, and vice versa. It can be seen that all AVs in both cases fall on the stable side, and they all tend to evolve to a more stable state (lower $|Ro_{v}|$ and further away from the stability limit). The stability results are in agreement with our observation that there is no apparent sign of instability of AVs except for at $z/h=0.12$ in Bu25, where AVs are still not independently distinguishable from the turbulent near wake. The more conservative determination of cyclo-geostrophic instability in the present wakes utilising (5.4) and (5.6) as compared with (5.3) also confirms the point of view in Lazar et al. (Reference Lazar, Stegner and Heifetz2013b) and Yim et al. (Reference Yim, Stegner and Billant2019) that, in real geophysical environments, stratification and vertical dissipation will further shrink the range of unstable vertical wavenumbers from the low- and high-wavenumber end, respectively, and lead to a greater range of stability.