2. Numerical methods

In this paper, we examine SRI amplitude modulations linked to spiral pattern changes in a Taylor–Couette system subjected to heating from above and cooling from below, which establishes a stable density stratification along the axial ( $z$ ) direction (i.e. with less dense fluid on top and more dense fluid as we move towards the bottom of the cylindrical cavity along its axial (z)-axis).

We will use the DNS solver developed by Abide et al. (Reference Abide, Viazzo, Raspo and Randriamampianina2018), considering the same physical model, numerical methods and experimental parameters as we used previously (Seelig, Harlander & Gellert Reference Seelig, Harlander and Gellert2018; Meletti et al. Reference Meletti, Abide, Viazzo, Krebs and Harlander2020, Reference Meletti, Abide, Viazzo and Harlander2023), so that the results can be consistent with previous findings. The code employs a fourth-order accurate spatial discretisation and high-performance computing (HPC) (Abide et al. Reference Abide, Binous and Zeghmati2017, Reference Abide, Viazzo, Raspo and Randriamampianina2018).

The physical model adopts a Taylor–Couette flow configuration filled with an incompressible fluid that is warmed on top and cooled at the bottom, producing a stable density stratification in the axial direction. The physical model is schematically presented in figure 1(a).

This is the same configuration used to study the SRI experimentally by Meletti et al. (Reference Meletti, Abide, Viazzo, Krebs and Harlander2020).

The fluid properties considered are those of the M5 silicon oil used in the experiments presented by Seelig et al. (Reference Seelig, Harlander and Gellert2018) and Meletti et al. (Reference Meletti, Abide, Viazzo, Krebs and Harlander2020), with kinematic viscosity $\nu = 5\times 10^{-6}\,{\mathrm m}^2\,{\mathrm s}^{-1}$ ; specific weight $\rho = 923\,{\mathrm{kg}}\,{\mathrm m}^{-3}$ ; coefficient of thermal expansion $\alpha =1.08\times 10^{-3}\,{\mathrm K}^{-1}$ ; thermal conductivity $k=0.133\, {\mathrm W}\,({\mathrm{K m})^{-1}}$ ; specific heat $c_p = 1630 \,{\mathrm J}\,({\mathrm{kg\,K}})^{-1}$ ; and Prandtl number $Pr = 57$ .

The vertical temperature gradient imposed is of $\Delta T/\Delta z = {5.71}\,{\textrm {K m}^{-1}}$ and ensures the stable axial density stratification. The geometrical parameters of the TC cylindrical cavity (also presented in the bottom part of table 1) are: inner cylinder radius $r_{\textit{in}}=75\,{\mathrm{mm}}$ , outer cylinder radius $r_{\textit{out}}=145\,{\mathrm{mm}}$ , gap size $d=(r_{\textit{out}}-r_{\textit{in}})=70\,{\mathrm{mm}}$ ; cylinders height $H=700\,{\mathrm{mm}}$ ; aspect ratio $\varGamma =10$ ; radii ratio $\eta \approx 0.52$ .

Table 1.

Flow parameters and non-dimensional numbers used in this study.

Notes: $r_{\textit{in}}=75\ \mathrm{mm}$ , $r_{\textit{out}}=145\ \mathrm{mm}$ , $d=r_{\textit{out}}-r_{\textit{in}}=70\,\mathrm{mm}$ , $\eta =r_{\textit{in}}/r_{\textit{out}}\approx 0.52$ , $H=700\ \mathrm{mm}$ , $\varGamma = H/d=10$ , $\alpha =1.08\times 10^{-3}\ \mathrm{K^{-1}}$ , $\partial _z T = 5.71\ \mathrm{K\,m^{-1}}$ , $\rho =923\ \mathrm{kg\,m^{-3}}$ , $k=0.133\ \mathrm{W\,m^{-1}K^{-1}}$ , $c_p=1630\ \mathrm{J\,kg^{-1}K^{-1}}$ .

The numerical code applied here solves the Navier–Stokes aligns using the Boussinesq approximation to incorporate buoyancy forces that, in a rotating frame with angular velocity $\boldsymbol{\varOmega }$ , and in non-dimensional form, read:

(2.1) \begin{align} \frac {\partial \tilde {\boldsymbol{u}}}{\partial \tilde {t} } + Ro ( \tilde {\boldsymbol{u}}\boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} ) \tilde {\boldsymbol{u}} + 2\boldsymbol{k} \times \tilde {\boldsymbol{u}} = -\tilde {\boldsymbol{\nabla }} \tilde {p} + \frac {1}{Fr^2Ro}\tilde {T}\boldsymbol{k} + E \tilde {\boldsymbol{\nabla} }^2 \tilde {\boldsymbol{u}}, \end{align}

where $\boldsymbol{k} \equiv \boldsymbol{\varOmega }/|\boldsymbol{\varOmega }|$ , $\tilde {\boldsymbol{\nabla }}$ is the non-dimensional gradient operator. $Ro$ is the Rossby number $Ro \equiv U/\varOmega L$ , $E=\nu /\varOmega L^2$ is the Ekman number and $ \textit{Fr} \equiv \varOmega /N$ is the Froude number, with $N$ being the buoyancy frequency, here defined as $N=\sqrt {\alpha \boldsymbol{ g} \partial T/\partial z}$ , where $\boldsymbol{ g}$ is the gravity acceleration. In our results, we consider the Reynolds number based on the inner cylinder rotation, as $ \textit{Re} = \varOmega _{\textit{in}} r_{\textit{in}} d/\nu$ . These flow parameters are presented in table 1.

For a direct comparison of experimental results (Meletti et al. Reference Meletti, Abide, Viazzo, Krebs and Harlander2020), the Navier–Stokes aligns under the Boussinesq approximation in our code are solved in the dimensional form and in a non-rotating (lab) frame, reading:

(2.2) \begin{align} \left \{ \begin{array}{ll} \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{u}} = 0 & \text{ in } D, \\[12pt]\displaystyle \partial _t \boldsymbol{u} + \frac {1}{2}\left [(\boldsymbol{ u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{ u} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{ u}\boldsymbol{ u}) \right ] = - \boldsymbol{\nabla }p + \nu \Delta \boldsymbol{ u} +\boldsymbol{F} & \text{ in } D,\\[12pt]\displaystyle \partial _t T + \frac {1}{2}\left [(\boldsymbol{ u} \boldsymbol{\cdot }\boldsymbol{\nabla }) T+ \boldsymbol{\nabla }(\boldsymbol{ u} T) \right ] = \kappa \boldsymbol{\nabla} ^2 T & \text{ in } D.\\ \end{array} \right . \end{align}

Here, $D$ denotes the computational domain, $\kappa$ is the fluid’s thermal conductivity, $p$ is the pressure, $T$ represents the temperature field, and $\boldsymbol{ u} = (u_r, u_\phi , u_z)$ denotes the velocity vector field in the radial, azimuthal and axial directions, respectively. The buoyancy force, denoted by $\boldsymbol{F}$ , results from density variations and is defined as

(2.3) \begin{align} \boldsymbol{F} = \alpha \boldsymbol{ g} \frac {\partial T}{\partial z}\Delta z . \end{align}

All analysis will consider intermediate Reynolds numbers of $ \textit{Re} = 400$ , with rotation ratio between inner and outer cylinders $\mu = \varOmega _{\textit{out}}/\varOmega _{\textit{in}} = 0.35$ . This value of $\mu$ is considered so that we impose a velocity profile that is almost Keplerian (slightly slower), to focus on accretion disk applications (see Visser & Dullemond Reference Visser and Dullemond2010; Lyra & Umurhan Reference Lyra and Umurhan2019).

In our simulations, velocity boundary conditions are applied at the cylinder walls, with the top and bottom boundaries rotating at the same angular speed as the outer wall. To impose the temperature gradient, the temperature at the top and bottom lids are prescribed as adiabatic Dirichlet boundary conditions. Assuming negligible heat loss through the lateral walls compared with the thermal forcing at the lids, adiabatic boundary conditions are applied laterally, at the inner and outer cylinder’s walls. Simulations begin with small, randomly distributed white noise perturbations to trigger instability development. Time discretisation is conducted using the methods proposed by Hugues & Randriamampianina (Reference Hugues and Randriamampianina1998). Spatial discretisation is enhanced using spectral Fourier methods in the azimuthal direction and a fourth-order compact finite difference scheme in the radial and axial directions, as described by Abide & Viazzo (Reference Abide and Viazzo2005). High-performance computing techniques, as detailed by Abide, Binous & Zeghmati (Reference Abide, Binous and Zeghmati2017), facilitate significant reductions in simulation time through parallelisation. Following Lopez & Marques (Reference Lopez and Marques2020) and Lopez & Marques (Reference Lopez and Marques2022), centrifugal buoyancy effects did not present a strong influence in our simulations once the instability was already established and they were therefore omitted from the numerical model. The simulations employ a grid resolution of $32 \times 64 \times 200$ in the $\phi \times r \times z$ directions.

3. Axial modes activation and spiral components separation

In this section, we will present the interaction of the SRI with the base flow and how this is related to the spiral pattern changes. We will further discuss how this is connected to the activation of different axial wavenumbers. Moreover, we discuss how the oscillating behaviour of the spiral axial propagation can be separated into two components, each of them related to weak nonlinear features.

The initial detailed comparison of experimental and numerical data by Meletti et al. (Reference Meletti, Abide, Viazzo, Krebs and Harlander2020) identified low-frequency amplitude modulations, as those shown here in figure 2. Similar patterns were later obtained numerically by Lopez & Marques (Reference Lopez and Marques2022) for a short cylinder geometry. These amplitude modulations occur across all velocity components ( $u_\phi$ , $u_r$ and $u_z$ ) and the temperature $T$ , affecting thus the temporal variations of the circulation and the stratification.

Figure 3.

Comparison of the time average axial velocity profiles with time averages taken during an upward travelling spiral period, during a downward travelling spiral period and during the transition from an upward to a downward spiral period. The black dashed line shows the time average over four full periods of amplitude modulations, revealing upward, downward and standing spiral patterns. The results are from numerical simulation performed with $ \textit{Re}=400$ , $\mu =0.35$ and $\Delta T/\Delta z = {5.7\,}{\textrm {K m}^{-1}}$ at a fixed radial position $r = r_{\textit{in}}+d/3$ and $\phi =0$ . The inset in each image shows values of $f-\overline {f}^t$ , i.e. the mean values during the upward, downward and transition spiral travelling regime minus the average during the full simulation period.

Figure 3 shows the time-averaged velocity profiles in the axial direction at a fixed radial position $r = r_{\textit{in}} + d/2$ . The inset in each figure shows $u-\overline {u}^t$ , i.e. the velocity minus the total mean value (which is represented by the black dashed line). The time-averaged azimuthal velocity ( $\overline {u_\phi }$ ), displayed in figure 3(a), shows that during the upward or downwards travelling regime (when the velocity amplitude of is larger), the azimuthal flow is slightly accelerated with respect to the mean values or compared with the $u_\phi$ during the spiral transition between these two regimes. The radial velocity profiles ( $\overline {u_r}$ ) presented in figure 3(b), instead, do not change significantly during the periods when the spirals travel either upwards or downwards, showing no great differences from the overall time average. However, figure 3(c) illustrates how the time-averaged axial velocity ( $\overline {u_z}$ ) is affected by the instability. The mean axial flow becomes positive during upward spiral propagation and negative during downward spiral propagation. During the transition between the two propagation patterns, the average axial velocity approaches zero. These variations in the time-averaged axial velocity are important because they show that the spiral propagation affects the mean flow. This connection highlights the interaction between the base flow and the emerging instabilities. Figure 3(d) shows the temperature profiles (and their fluctuations). In these cases, we see some difference in temperature fluctuations, of approximately 0.5 °C, especially near the cavity lids, which are considered to be small. Note that near the top and bottom boundaries in figure 3, there is a strong radial inward flow caused by Ekman circulation, with a strong upward flow near the bottom wall and a downward flow near the top wall (shown in the radial velocity presented in figure 3 c). These Ekman effects are qualitatively similar to those reported by Lopez & Marques (Reference Lopez and Marques2020). Despite these Ekman effects, simulations with periodic boundary conditions at the top and bottom lids (not presented here) also exhibited amplitude modulations associated with pattern changes (see Meletti et al. Reference Meletti, Abide, Viazzo and Harlander2023), showing that the presence of lids and Ekman effects is not the reason for the modulations and pattern transitions to occur.

To better understand the slow axial mean flow oscillations observed in our SRI data, we plot the time series of the eddy momentum flux divergence and the axial mean flow in figure 4. More precisely, at a fixed azimuthal position $\phi =0$ and at $r=r_{\textit{in}}+d/2$ , we take the axial mean of the momentum flux $u'_r u'_z$ (primes denote the eddy components $u'=u-\overline {u}^t$ , i.e. the deviations from the mean) at each time step and compute the derivative of this mean flux with respect to $r$ . In wave–mean flow interaction theory, this eddy flux divergence is the forcing term for the axial mean flow. Clearly, there is a correlation between the axial mean flow and the flux divergence, which suggests that the eddies in the form of the spiral waves force the mean flow oscillations, as will be discussed in more detail in § 5. Note that the time-series were filtered with a 15 min moving average to highlight the low-frequency variations of the mean flow. This time window was chosen arbitrarily. Other periods were also tried for the windowing and the qualitative features observed do not change, considering our purposes here, i.e, using it as a low-pass filter to highlight the low-frequency features.

Figure 4.

Time series of the wave momentum flux divergence $\partial _r (\langle u_r^\prime u_z^\prime \rangle _z)$ and the axial mean flow $\langle u_z\rangle _z$ . The values are averaged in the axial direction and taken at a fixed radial position $r=r_{\textit{in}} + d/2$ and $\phi =0$ . For wave–mean flow interactions, the momentum flux divergence is the forcing of the mean axial flow. The amplitude of both time series was normalised by the maximum amplitude in the established low-frequency region (for $t\gt {100}\,{\min }$ ).

Figure 5.

Two-dimensional power spectra of $u_\phi (r)$ , computed from the space–time diagrams shown in figure 2(b–d) when the spiral is propagating (a) upwards; (b) during transition from an upward to a downward propagation regime and (c) downwards. The $x$ -axis is the frequency in $\mathrm{Hz}$ , while the $y$ -axis represents axial wavenumber $k$ (axial modes). The spectra amplitudes are normalised by the maximum amplitude value of the upward (and downward) propagating spirals $P_{0,max}$ . The spectra are obtained in a frame of reference fixed in the laboratory. During the transition, the maximum amplitude of the spectra was half of the maximum amplitude found while the spiral was travelling upward or downward ( $P/P_{0,max} = 0.5$ ).

To investigate if there are changes in the SRI frequency and in the axial wavenumbers during these different spiral propagating regimes, we computed the 2-D-FFT obtained from $u_\phi$ space–time diagrams in the axial direction (presented in figure 2). Figure 5 shows the results obtained at moments when the spiral is travelling upwards (figure 5 a, on the left-hand side), downwards (figure 5 c, on the right-hand side), and during the transition from the upward to the downward propagation (figure 5 b, middle image). The diagrams in figure 5 present the frequencies $f$ on the $x$ -axis and the axial wavenumber $k$ in the $y$ -axis. We note that the peaks have the same SRI frequencies of $f=0.032$ Hz independently of the spiral direction of propagation, but they change their wavenumbers from positive to negative values, with different modes being activated and suppressed. During the upward travelling spiral, the stronger wavenumber activated in figure 5(a) is $k_{up}=4$ , with maximum amplitude $P_{0,max}$ , and the downward mode, which has wavenumber $k_{down}=-4$ , is weakly activated. The wavenumber activated with $P_{0,max}$ in figure 5(c), related to the downward propagating spiral, is $k_{down}=-4$ , while $k_{up}=4$ has smaller amplitude. In figure 5(b), both modes $k_{up}=4$ and $k_{down}=-4$ are activated, each with approximately half of the maximum amplitude of the spiral travelling upward or downward ( $P_{0,max}/2$ ).

To better understand the spiral oscillatory behaviour and its associated mode activation during the different spiral patterns observed, we separated the up- and downward propagating spirals by applying the RT to the full signal, capturing phases during dominating up, downward and mixed spiral propagation. The RT is a Fourier-like technique to select wave components with different directions of propagation. The RT is particularly suited for finding individual waves that compose noisy or irregular fields (Almar et al. Reference Almar, Michallet, Cienfuegos, Bonneton, Tissier and Ruessink2014). These techniques are interesting for evaluating the results directly using the data obtained, without knowing the wave’s dispersion relation and without the necessity of an analytical model of the SRI. A brief description of the RTs can be found in Appendix A.

Figure 6.

Separation of upward and downward axial travelling components in $u_\phi$ space–time diagram using the RT. $u_\phi$ has been taken at $(r,\phi )=(r_{\textit{in}}+d/2,0)$ . Panels (a), (c), (e) show time intervals when the spiral is travelling upwards. Panels (b), (d), (f) show time intervals when the spiral is travelling downwards. Panels (a), (b) show the full space–time diagram minus the mean flow (computed using the full time signal). Panels (c–f) show the separated upward and downward travelling components obtained using the RT (applied to the full signal, but it filters out the mean flow when it is applied). Simulation performed with $ \textit{Re}=400$ , $\mu = 0.35$ , $\Delta T/\Delta z \approx {5.71}{\textrm {K m}^{-1}}$ and $H={700}\,\textrm {mm}$ .

Figure 6 then shows the separation of the upward and downward components of $u_\phi$ on a space–time diagram, while the spiral is travelling upward (left-hand side) and downward (right-hand side) using the RT. The results are from the simulation with $ \textit{Re}=400$ , $\mu =0.35$ and $\Delta T/\Delta z \approx {5.71}\,{\textrm {K m}^{-1}}$ presented in figure 2. After the two wave fields have been separated, we again computed the 2-D-FFT spectra from the corresponding space–time diagrams (not shown here). We observed that the 2-D-FFT of the downward spiral component is exactly the same as the bottom wavenumber in the full 2-D-FFT spectrum presented in figure 5, but the positive frequency is no longer observed (as expected). The upward wavenumber in figure 5 was also captured in the spectrum of the upward component of the upward-travelling spiral (figure 6 e), but the downward travelling spiral component is removed. Thus, on a simplified model, the wavenumbers we observe in figures 5 can be associated with a superposition of one upward and one downward spiral of axial wavenumber $k=4$ and $k=-4$ propagating in time with the SRI frequency (see Meletti et al. Reference Meletti, Abide, Viazzo, Krebs and Harlander2020 for more details on the SRI frequency measurements and values). In other words, while the spiral is propagating downwards, we see in figure 6 that the spiral travelling upwards is suppressed, reaching smaller amplitudes and a more vertical inclination (figure 6 f). The same occurs with the downward component when the spiral is propagating upwards (figure 6 c). This approach fully confirms the results shown in figure 5, namely that each separated mode is indeed associated with the spiral components travelling upwards and downwards without any changes in the frequency, but with changes in their wavenumbers. This implies that a linear superposition of the two spirals should explain some part of the amplitude modulation. As an extension of the analysis, we also examined the RT on a more complex case involving a taller cavity with $H = {2800}\,\textrm {mm}$ (four times taller than the configuration considered here). In this set-up, the flow exhibits more complicated spiral patterns, but all the conclusions we can obtain from those results are consistent with those presented here, i.e. we observe similar regions being activated/deactivated depending on whether the spiral is moving upwards or downwards in a given time. These additional results are presented in Appendix B.

Figure 7(a) shows the time series envelope of separated upward and downward spiral components. Note that this is not the full SRI velocity, which presents faster oscillations (see Meletti et al. Reference Meletti, Abide, Viazzo, Krebs and Harlander2020), but just the envelope capturing the velocity amplitude modulations. The upward travelling time series was arbitrarily dislocated in the vertical axis for better visualisation (originally, both time series were on top of each other). Note that, also here, when the amplitude of the downward spiral is enhanced, the upward spiral amplitude becomes smaller (and vice versa), showing a phase shift of the up-downward components beating. When the power spectrum of the $u_\phi$ amplitude envelope is compared with the separated upward and downward components obtained by using the RT in figure 7(b), the same (low) frequencies are observed.

Figure 7.

(a) $u_\phi ^\prime$ upward and downward spiral components amplitude modulations taken at $(r,\phi ,z)=(r_{\textit{in}}+d/2,0,H/2)$ . The upward travelling time series was arbitrarily dislocated in the vertical axis for better visualisation (originally, both time series were on top of each other). (b) Power spectra obtained from the amplitude envelopes of the $u_r$ time series, and of their separated upward and downward components, obtained using the RT to separate the signals. The time series is obtained from numerical simulations with $ \textit{Re}=400$ , $\mu =0.35$ and $\Delta T/\Delta z \approx {5.71}\,{\textrm {K m}^{-1}}$ .

Note that we can observe harmonics in the envelope spectra presented in figure 7(b), and also on each of the separated up and downward spiral components, which suggests a weak nonlinear interaction and not simply a linear interaction of two waves travelling with different frequencies. Note also that the amplitudes of each separated spiral component in the power spectra differ from those obtained from the full velocity signal. A similar behaviour is obtained from the $u_\phi$ , $u_r$ and $u_z$ time series, i.e. with peaks corresponding to the same frequencies in $\textrm {Hz}$ , but with different amplitudes $P$ and $P_{max}$ . However, most importantly, we would like to highlight that, when the Radon filtering is applied, we observe a clear phase shift in the modulations of each separated component when we look at the time series in figure 7(a). This phase shift indicates that the harmonics observed in the spectra are not simply modulated in isolation; rather, they suggest a nonlinear interaction, causing the modulations. The fact that these modulations cannot be fully explained by linear interactions alone implies that a more complex dynamic is at play, where the interaction between the spirals and the mean flow results in the alternating strengthening and weakening of the SRI spiral components. This interaction, in turn, leads to the observed alternation between upward and downward propagating modes, each with varying amplitudes. In the following sections, we will investigate how these two individually modulated components could lead to the pattern formations we observed.

4. Toy model: wave-like spiral propagation

In this section, we introduce a toy model consisting of two waves travelling in opposite axial directions to demonstrate how the linear superposition of the SRI spirals can lead to pattern changes. This approach is inspired by the separation of the spirals using the RT technique. We observed in the previous sections that each spiral (travelling upward and downward) is individually modulated and phase-shifted relative to one another. The reasons behind the modulation of each spiral will be discussed later in this paper. For now, we assume that the spirals behave as two individual plane waves, described by the following aligns:

(4.1) \begin{align} \begin{aligned} & \text{wave}_1 = A_1 \cos \left ( \left ( m_1x + l_1y + k_1 z \right ) - \omega t \right ), & \\ & \text{wave}_2 = A_2 \cos \left ( \left ( m_2x + l_2y + k_2 z \right ) - \omega t \right ), & \\ & u_{toy} = \text{wave}_1 + \text{wave}_2 , & \\ \end{aligned} \end{align}

with $0 \leqslant x,y,z \leqslant 2 \pi$ , and amplitudes $A_1$ and $A_2$ . The values of $0\lt x,y,z\lt 2\pi$ result from normalising the cavity height (e.g. ${0}\,\textrm {mm} \lt z \lt {700}\,\textrm {mm}$ ). The amplitude of each plane wave is modulated, and they must be out of phase to achieve the inclined spirals with constructive and destructive interference while they propagate. In the toy model, sinusoidal amplitude modulations $A_1$ and $A_2$ are considered out of phase with an angle $\theta$ , written as

(4.2) \begin{align} \begin{aligned} & A_1 = A \sin(\omega _A t) , & \\ & A_2 = A \sin(\omega _A t+\theta ) , & \\ \end{aligned} \end{align}

where $A$ is a given real value, and $\omega _A \lt \lt \omega$ is the amplitude modulation of each wave, here considered to be the same for $A_1$ and $A_2$ .

Figure 8.

$u_{toy}$ space–time diagram of the toy model composed of two plane waves with sinusoidal amplitude modulations with $\omega _A=7\times 10^{-4}$ , out of phase by an angle $\theta =\pi /3$ , and travelling in opposite axial directions with wavenumbers of $\text{wave}_1$ and $\text{wave}_2$ respectively $(m_1,l_1,k_1)= (1,1,4)$ and $(m_2,l_2,k_2)=(1,1,-4)$ . The frequency $\omega = 0.03$ and the maximum amplitude of each wave is $A=10$ . Data taken at $(x,y)=(0,\pi )$ .

Figure 9.

Snapshots with different spiral patterns in the $\phi$ – $z$ cross-section comparing $u_\phi ^\prime = u_\phi -\overline {u_\phi }$ obtained from(a–c) numerical simulations fixed at a radial position $r\approx r_{\textit{in}}+d/3$ and(d,e,f) the toy model. Panels (a), (d) show moments when the spirals are travelling downwards; panels (b), (e) show the transition; and panels (c), (f) show spirals travelling upwards. The simulations were performed with $ \textit{Re}=400$ , $\mu =0.35$ and $\Delta T/\Delta z \approx {5.71}\,{\textrm {K m}^{-1}}$ . The toy model consists of two plane waves at $y=\pi$ with frequency $\omega =0.001$ , and wavenumbers $(m_1,l_1,k_1)= (1,1,4)$ and $(m_2,l_2,k_2)=(1,1,-4)$ , wave amplitude $A=3$ mm s⁻¹, and modulation frequency $\omega _A=0.01$ with $\theta =\pi /2$ phase difference.

Figure 8 shows the space–time diagram obtained with this toy model, with wavenumbers in the azimuthal, radial and axial directions ( $m=1,l=1,k=4$ ) and ( $m=1,l=1, k=-4$ ), similar to those previously observed in our simulations (with $ \textit{Re}=400$ , $\mu =0.35$ , $\Delta T/\Delta z \approx {5.71}\,{\textrm {K m}^{-1}}$ and $H={700}\,\textrm {mm}$ ). It is possible to see that the linear superposition of both upward and downward waves, each individually modulated and travelling out of phase, could lead to the final spiral pattern transitions observed in the previous sections. The amplitude modulations of the waves presented in figure 8 are out of phase with an angle $\theta =\pi /3$ , but other different phase shifts (and different wavenumbers) produce similar pattern changes.

The toy model and the numerical simulation also show good qualitative agreement when we compare snapshots, i.e. looking at the space structures at given times while the spirals are travelling upward, downward and during the transition. This comparison can be seen in figure 9, where we can see that the results are similar, except for the fact that the spirals in the simulations are confined in a slightly smaller region, due to Ekman effects. The good qualitative agreement of this simplified toy model with the numerical simulations, as well as the possibility of reproducing the spiral pattern changes previously investigated using similar wavenumbers and frequencies obtained from numerical simulations and experimental measurements, shows that this linear superposition of the spirals can drive the spiral pattern changes. In this case, each up and downward component should be interacting with the mean flow that, at times, provides more energy to the upward travelling spiral, and at other times, provides more energy to the downward component. The reason why each individual spiral is modulated will be interpreted as a QBO-like mechanism in the following section, inspired by the fact that the spiral propagation direction affects the mean flow structure (see figure 3). Therefore, each spiral component must be interacting with the mean flow to achieve its individual (phase-shifted) modulation, so that the linear superposition of these two modulated components will lead to the reversal of the spiral direction as presented here in this toy model.

5. Axial mean flow interpretation considering inertial wave interactions

Based on previous observations of mean flow variations associated with the oscillatory spiral behaviour and given that the toy model introduced in the previous section successfully reproduces these pattern changes, we will now interpret the low-frequency amplitude modulation of the SRI as a QBO-like phenomenon. This interpretation is suggested by figures 3(c) and 4, as we will discuss in the following.

We would like to emphasise that employing a QBO-like model to interpret our SRI data is a simplified approach intended to understand the mechanisms responsible for the amplitude modulations we observed. This analogy is motivated by certain dynamical similarities between the SRI and the QBO.

5.2. Derivation of the QBO-like aligns applied to the SRI reversals

To consider these modulations as a QBO-like phenomenon occurring in stratified rotating shear flow, our approach follows the ideas presented by Holton & Lindzen (Reference Holton and Lindzen1972), Plumb (Reference Plumb1977), Seelig & Harlander (Reference Seelig and Harlander2015) and Renaud & Venaille (Reference Renaud and Venaille2020), but focuses on inertial waves due to the weak stratification in the SRI simulation. Instead of considering the classical non-rotating but stratified QBO scenario, where the base flow in the azimuthal direction produces a QBO in the azimuthal–radial plane, we consider the base flow in the axial direction generating the QBO in the axial–radial plane. In other words, we will show here how the wave interactions in the axial–radial plane can drive an oscillating axial mean flow. Since, for the flow discussed here, the ratio between the Coriolis parameter $f=2 \varOmega _{\textit{in}}$ and the buoyancy frequency $N$ is larger than 1 ( $f/N=3$ ), we focus here on inertial waves and not on internal gravity waves as in the model by Holton & Lindzen (Reference Holton and Lindzen1972) and Plumb (Reference Plumb1977). However, we will see that the aligns are analogous to aligns for the atmospheric QBO in the equatorial stratosphere, driven by internal gravity waves. The analysis is local Cartesian, as presented in figure 10; that is, we neglect curvature. For further simplification, the diffusion terms in the Plumb (Reference Plumb1977) model are replaced by viscous drag terms with a drag coefficient $\varsigma$ . Note that we neglect the azimuthal $x$ -dependency in the inertial wave-based QBO model by assuming that this dependency is weak. Only then can we derive a model that is isomorphic to the classical 2-D QBO model. This is in analogy to the internal gravity wave model of Lindzen & Holton (Reference Lindzen and Holton1968), Holton & Lindzen (Reference Holton and Lindzen1972) and Plumb (Reference Plumb1977), for which the meridional dependency had been skipped. Then, the momentum aligns in a frame co-rotating with $\varOmega _{\textit{in}}$ read

(5.1) \begin{align} \partial _t u + v\partial _y u + w\partial _z u = -2\varOmega _{\textit{in}} w -\varsigma u, \\[-28pt] \nonumber \end{align}

(5.2) \begin{align} \partial _tv+v\partial _yv+w\partial _zv = -\frac {1}{\rho _0} \partial _y p -\varsigma v, \\[-28pt] \nonumber \end{align}

(5.3) \begin{align} \partial _tw+v\partial _yw+w\partial _zw = -\frac {1}{\rho _0}\partial _z p + 2\varOmega _{\textit{in}} u - \varsigma w, \\[0pt] \nonumber \end{align}

where $(y,z)$ are the axial and radial directions, respectively, and $(u,v,w)$ are the azimuthal, axial and radial velocity components (see figure 10). Here, $p$ is the generalised pressure that includes centrifugal effects. Note that, in contrast to the strictly 2-D internal wave QBO model, for inertial waves, the flow is 3-D even though all dependent variables depend only on $t$ , $y$ and $z$ .

Figure 10.

Local Cartesian coordinate system used for (5.1)–(5.3). Waves propagating in the $z$ – $y$ -plane are sketched. Note that in such a ‘planetary model’, the rotation vector and the gravity vector are perpendicular to each other. In contrast, these vectors would be parallel and in the opposite direction in the SRI-cylinder. In our model, we neglected the explicit gravity term for simplicity.

Defining the streamfunction $\psi =\psi (y,z,t)$ and a zonal momentum $\tilde u$ as

(5.4) \begin{eqnarray} v=-\partial _z \psi , \quad w=\partial _y\psi , \quad \tilde {u}=2\varOmega _{\textit{in}} u, \end{eqnarray}

this system can be reduced to

(5.5) \begin{eqnarray} \partial _t\boldsymbol{\nabla} ^2\psi +J(\psi ,\boldsymbol{\nabla} ^2\psi ) & = & \partial _y\tilde {u}-\varsigma \boldsymbol{\nabla} ^2\psi , \\[-28pt] \nonumber\end{eqnarray}

(5.6) \begin{eqnarray} \partial _t\tilde {u}+J(\psi ,\tilde {u}) & = & -4\varOmega ^2_{\textit{in}}\partial _y\psi -\varsigma \tilde {u} . \\[0pt] \nonumber \end{eqnarray}

By introducing the streamfunction in the azimuthal–radial plane,

(5.7) \begin{eqnarray} \psi = \psi (x,z,t), \quad u = -\partial _z\psi , \quad w = \partial _x \psi , \end{eqnarray}

and with a gravity vector pointing in the negative $z$ -direction, Holton & Lindzen (Reference Holton and Lindzen1972) and Plumb (Reference Plumb1977) derived the following align for the case of internal gravity waves:

(5.8) \begin{eqnarray} \partial _t\boldsymbol{\nabla} ^2\psi +J(\psi ,\boldsymbol{\nabla} ^2\psi ) &=& \partial _xb+\nu \boldsymbol{\nabla} ^4\psi, \\[-28pt] \nonumber \end{eqnarray}

(5.9) \begin{eqnarray} \partial _{t}b+J(\psi ,b) & = & -N^{2}\partial _{x}\psi -\mu b, \\[0pt] \nonumber\end{eqnarray}

where $b=-g \rho '/\rho _0$ is buoyancy, $\rho '$ and $\rho _0$ is the density perturbation and the mean density, $\nu$ is the kinematic viscosity, and $N^2=-(g/\rho _0) \text{d} \bar \rho /\text{d}z$ is the square of the buoyancy frequency, where $\bar \rho$ is the background density with a linear dependency on $z$ . Here, there is no rotation, and in contrast to the previous case, the waves propagate in the azimuthal–radial plane (the $x$ – $z$ -plane). Except the linear friction term in (5.5), the (5.8), (5.9) and (5.5), (5.6) are mathematically isomorphic. Friction terms are also used, e.g. in simple models for vorticity Ekman pumping (Pedlosky Reference Pedlosky1987).

To derive a ‘QBO model’ from (5.5) and (5.6), we can proceed in the same way as for the gravity waves. First, a mean state $\bar \psi$ is defined, the streamfunction is written as $\psi (y,z,t) = \bar \psi (z)+\epsilon \psi '(y,z,t)$ and the aligns are linearised about the mean streamfunction $\bar \psi$ , where $\psi '(y,z,t)$ is the fluctuation around the mean value and $\epsilon$ is a small perturbation parameter, $\epsilon \ll 1$ . The perturbation is expanded as $\psi '(y,z,t)=\phi (z) \exp (\mathrm{i}l(y-ct))$ , where $l$ is the axial wave number and $c$ the axial phase speed. We obtain a linear align for $\phi (z)$ , and its solution gives the perturbation velocities $v'$ and $w'$ . Then, assuming periodic boundary conditions in the axial $y$ -direction, (5.2) is averaged over this direction to obtain

(5.10) \begin{align} \partial _{t} \bar v +\varsigma \bar v=-\epsilon ^2 \partial _{z} (\overline {v'w'}), \end{align}

where $\bar v=-\partial _z \bar \psi$ is the slowly varying and $z$ -dependent mean flow in the axial $y$ -direction. Note that the forcing term on the right-hand side of (5.10) is analogous to the eddy forcing term we computed for the SRI simulations and plotted in figure 4. Note further that even though the SRI-cylinder is bounded by a top and bottom lid, for deriving (5.10), we considered periodicity in the axial direction. This simplification is motivated by findings of Meletti et al. (Reference Meletti, Abide, Viazzo and Harlander2023) showing that the slow axial flow modulations also occur under axially periodic boundary conditions. To draw an analogy with Plumb (Reference Plumb1977), we will now focus here on the QBO-like phenomenon induced by a single boundary (inner cylinder wall).

5.3. QBO-like oscillations induced by the inner boundary

Plumb (Reference Plumb1977) found for the atmospheric QBO case,

(5.11) \begin{align} \partial _{t} \bar u -\nu \partial _{zz} \bar u = -\epsilon ^2 \partial _{z} (\overline {u'w'}), \end{align}

i.e. a slowly varying mean flow in the $x$ -direction driven by waves in the $x$ – $z$ -plane. Again, we see a strong analogy between (5.10) and (5.11). Holton & Lindzen (Reference Holton and Lindzen1972) and Plumb (Reference Plumb1977) were able to ‘parametrise’ the wave momentum flux as

(5.12) \begin{align} (\overline {u'w'}) \sim \exp \left ( -s \int _0^z \frac {N \mu }{k (\bar u-c)^2}\,\text{d}z' \right ), \end{align}

where $s$ is the sign of the $z$ -component of the group velocity, $k$ is the wavenumber in the $x$ -direction and $c=\omega /k$ is the phase velocity, where $\omega$ is the wave frequency. In analogy, we can write for the inertial wave case,

(5.13) \begin{align} (\overline {v'w'}) \sim \exp \left ( -s \int _0^z \frac {2 \varOmega \zeta }{l (\bar v-c)^2}\,\text{d}z' \right ). \end{align}

Considering two waves propagating towards the positive and negative $y$ -direction, as proposed in § 4, the mean flow align can be written in the universal non-dimensional form,

(5.14) \begin{align} \partial _{t} \hat {\bar v} + Re^{-1} \hat {\bar v}=-\frac {\partial }{\partial \hat z} \left ( \exp \left ( - \int _0^{\hat z} \frac {1}{\tilde {l} (\hat {\bar v}-1)^2}\, \text{d}\hat z' \right ) - \exp \left ( - \int _0^{\hat z} \frac {1}{\tilde {l} (\hat {\bar v}+1)^2}\, \text{d} \hat z' \right )\right ), \end{align}

where $ \textit{Re}^{-1}=\varsigma /\varOmega _Q$ , comparing the damping frequency with the QBO-like frequency $\varOmega _Q$ . Note that for the internal gravity wave forcing, Renaud & Venaille (Reference Renaud and Venaille2020) scaled the Reynolds number as $ \textit{Re}^{-1}=\omega ^2 h^2/(2 \nu \gamma )$ , where $\omega$ and $h$ is the wave’s frequency and amplitude,and $\nu$ and $\gamma$ the kinematic viscosity and the buoyancy damping constant, respectively. In contrast to the classical QBO case, as previously mentioned, the waves in the SRI case are not forced by topography but by a flow instability. Thus, the wave amplitude depends on flow parameters like the damping constant and the angular velocity of the inner cylinder, and is also controlled by a nonlinear saturation of the growing waves, making it more difficult to estimate and include its value in the model. To bring (5.10) in the non-dimensional form (5.14), we then used the scaling $\hat t=\varOmega _Q t, \hat {\bar v}=\bar v/c, \hat z=z \varOmega _Q/c, \hat w'=\epsilon w'/c, \hat v'=\epsilon v'/c$ , and for the axial wavenumber $\hat l=l c/(2 \varsigma )$ and $\tilde {l}=\hat l \varOmega _q/\varOmega$ . Here, $c$ is the typical wave phase speed in the axial $y$ -direction. Note that these scalings lead to a formal similarity between (5.14) and (37) of Renaud & Venaille (Reference Renaud and Venaille2020). However, we highlight that, in the SRI context, the Reynolds number $ \textit{Re}=\varOmega _Q/\varsigma$ cannot be directly treated as an externally controlled parameter spanning from 0 to infinity, as in the Holton–Lindzen–Plumb framework of Renaud & Venaille (Reference Renaud and Venaille2020). In SRI experiments, $ \textit{Re}$ is a non-trivial function that depends on other control parameters, such as the angular frequency of the modulations $\varOmega _Q$ and friction $\varsigma$ . As a result, these parameters also influence the wave amplitude in a non-trivial way. In the regimes explored here, increasing $ \textit{Re}$ (interpreted as the forcing parameter in Plumb’s model) does not lead to a quasi-periodic route to chaos, as observed by Renaud & Venaille (Reference Renaud and Venaille2020), but may instead promote flow stabilisation, as described by Rüdiger et al. (Reference Rüdiger, Seelig, Schultz, Gellert, Egbers and Harlander2017) and Meletti et al. (Reference Meletti, Abide, Viazzo, Krebs and Harlander2020).

To give an impression of the range of the non-dimensional parameters, we can use dimensional values from the numerical simulations. The mean velocity is, depending on the $z$ -position, $\bar v \lessapprox 2 \times 10^{-3}$ m s $^{-1}$ , $(\omega , \varOmega , \varOmega _Q) \approx (2 \times 10^{-1}, 4 \times 10^{-1}, 2 \times 10^{-3})$ rad s $^{-1}$ , and the wavelength is $\lambda \approx 0.14$ m implying $l \approx 45$ rad m $^{-1}$ , which gives an axial wave phase speed of $c \approx 4.4 \times 10^{-3}$ ms $^{-1}$ . The solution of the non-dimensional version of (5.14) is controlled by $ \textit{Re}^{-1}$ and $\tilde {l}$ , where the former determines the flow regime (stable or oscillatory) and the amplitude, and the latter controls the period of the oscillations. We used a damping rate close to the angular frequency of the slow oscillations $\varsigma =2 \varOmega _Q$ to be in the oscillatory regime. With the values above, $\tilde {l}$ is of order 1.

Note that using (5.14) but with only a single-wave forcing on the right-hand side, the mean flow would converge to a steady state with a narrow boundary layer at the inner cylinder wall and a constant flow above this layer. Adding a second wave, the mean flow becomes unstable and starts to oscillate. More details can be found from Holton & Lindzen (Reference Holton and Lindzen1972), Plumb (Reference Plumb1977) and Renaud & Venaille (Reference Renaud and Venaille2020). Note that the mechanism described here has been reproduced in a laboratory experiment with internal gravity waves by Plumb & McEwan (Reference Plumb and McEwan1978). For a recent review, see also Semin & Pétrelis (Reference Semin and Pétrelis2024).

Figure 11.

(a) Radial distribution of the axial velocity averaged in the axial direction and in time $( \overline {\langle u_z \rangle _z}^t )$ at $(r,\phi )=(r_{\textit{in}}+d/3,0)$ considering three different time intervals: when the spiral is travelling upwards, downwards and during the transition. (b) Axial and time mean of the eddy momentum flux $( \overline {\langle u_r^\prime u_z^\prime \rangle _z}^t )$ .

Assuming wave forcing at $z=0$ , we solved (5.14) numerically using the fourth-order Runge–Kutta scheme for the time derivative and the trapezoidal rule for the integrals. To connect (5.14) more tightly to the simulation, we used the simulated axial mean flow profile $\overline {\left \langle u_z \right \rangle _z}^t$ as an initial condition. The profile considered is the one shown in figure 11(a) as the blue solid line. We further show the momentum flux taken from the SRI simulation in figure 11(b). The slopes of the curves shown give the forcing of the low-frequency oscillations. We see that at the inner part of the gap (left-hand side in figures 11 a and 11 b), the sign of the forcing, with positive/negative slope, differs from that at the outer part, with negative/positive slopes. For (5.14) of the QBO-like model, the sign differs between a region $0 \lessapprox z \lessapprox 0.1\boldsymbol{\cdot }\max (z)$ and $z \gtrapprox 0.1\boldsymbol{\cdot }\max (z)$ . This leads to a downward propagation of the QBO-like signal and usually different flow directions close to the bottom and further up (see, e.g. figure 3a of Renaud & Venaille (Reference Renaud and Venaille2020)).

The result of the integration is shown in figure 12(a), where we plotted the axial velocity ${\overline v}(t,z)$ (corresponding to $=\left \langle u_z \right \rangle _z(t,r)$ in the DNS) for $ \textit{Re}^{-1}=2$ and $\tilde {l}=1$ . It can be seen that over the entire period considered, the axial flow direction oscillates. The largest amplitudes can be found moving towards $z=0$ in the coordinate system presented in figure 10, which corresponds to the area near $r=r_{\textit{in}}$ in the DNS results. The amplitude weakens away from the inner boundary. Comparing the mean flow and the wave forcing of the QBO model (not shown) with the corresponding time-series of the SRI simulations shown in figure 4, we find qualitative agreements concerning the shape and the phase lag of the curves. This suggests that the mechanism we describe here using the simple QBO model is consistent with the data shown in figure 4 where we found a correlation between $\partial _{r} \langle u'_z u'_r \rangle _z$ and $\langle u_z \rangle _z$ , corresponding to $\partial _{z} (\overline {v'w'})$ and $\bar v$ of the QBO model. We hence conclude that the wave–mean flow interaction of the QBO model is also a relevant process for the oscillating axial mean flow in the SRI simulation.

Figure 12.

Comparison of space–time diagrams using a non-dimensional time $^*$ . The time is normalised for better comparison, such that the modulations of the QBO model and the DNS are the same. (a) QBO model; (b) SRI axial velocity $u_z$ modulations at mid-height position $z=H/2$ ; (c) QBO-like model considering two boundaries by mirroring the results of panel (a) at the top and phase shifting them ${53}{^\circ }$ . By considering two boundaries, the results of the QBO-like model align more closely with the low-pass-filtered SRI data shown in panel (b). (a) Space–time diagram of ${\overline v}(t,z)$ of the QBO-like model in the radial direction, (b) SRI space–time diagram of $u_z(t,r)$ at $\phi =0,z=H/2$ using a low pass filter to highlight the low frequencies, (c) Space–time diagram when superposing two QBO solutions at the opposing boundaries and applying a phase difference between the solutions.

5.4. Comparison of the SRI data with the QBO-like model

A limitation of the model sketched is the neglect of the radially bounded nature of our SRI system, where waves near the inner and outer boundaries strongly influence the flow. This effect is clear when we compare the $t$ – $r$ -Hovmöller diagrams of the QBO-like model in figure 12(a) with the low-frequency modulations of $u_z$ shown for the SRI simulation in figure 12(b). The modulation amplitude is large at both radial boundaries since both are source regions for waves. A first step to adapt (5.14) to the SRI simulation is to solve it for a wave source at $z=0$ and at $z=1$ , and then superpose both solutions. For the SRI, the waves at the inner and outer boundary are phase shifted, as can be seen, e.g. in figure 6 of Yavneh et al. (Reference Yavneh, McWilliams and Molemaker2001) or figure 2 of Molemaker et al. (Reference Molemaker, McWilliams and Yavneh2001). In fact, for an instability due to wave resonance, a phase shift between the resonating waves is rather typical. Hence, we applied a ${53}{^\circ }$ phase shift to the waves at $z=1$ . Note that we applied the phase shift heuristically, lacking a theory for the consequences of this shift to the low-frequency oscillations. It is, however, empirically motivated by SRI data, where there is a phase shift between the wave patterns near the inner and outer cylinders. Thus, in this framework, we do not solve the aligns at the outer boundary; instead, we simply mirror the inner-cylinder solution to the outer wall and heuristically apply a phase shift, which is then added linearly to the QBO-like solution (presented in figure 12 a) to obtain a qualitative representation of the outer cylinder’s contribution. The result of superposing the forcing at both boundaries and applying a phase shift is shown in figure 12(c), where the space–time structure of the modified QBO model becomes more consistent with the bounded SRI results. Despite these limitations, the general features of figures 12(b) and 12(c) are not very different: both systems display QBO-like oscillations, and in both cases, boundary influences extend only until a finite distance from the boundaries.

Finally, in figure 13, we show a comparison of a time series of a mirrored QBO-like result similar to figure 12(c) with the corresponding SRI time series. Here, the comparison is done in dimensional units using $ \textit{Re}^{-1}=2$ and $\tilde {l}=1.35$ . The time series is taken both at mid-gap ( $r=r_{\textit{in}}+d/2$ ) and, for the DNS, at mid-height position ( $z=H/2$ ). In dimensional units, the agreement between the time series is very good, although the frequency of the QBO-like model is somewhat slightly higher.

In summary, we note that the real SRI configuration, as captured by DNS, is naturally more complex than the one modelled with the reduced QBO-like model. A quantitative analysis of the SRI wave–mean flow interactions would require an appropriately formulated numerical model, starting from (2.2). Therefore, we do not attempt to justify the realism of the inertial wave–mean flow model in modelling the SRI spiral pattern changes. Rather, we use this simplified model to show how it produces similar results that support an interpretation of the spiral propagation changes in the axial direction with QBO-like features, illustrating how the low-frequency oscillations in the SRI can be interpreted as resulting from mean-flow/instability interactions.

Figure 13.

Comparison of dimensional time series between the mirrored QBO-like model (5.14) and the SRI axial velocity averaged in the azimuthal direction $\langle u_z\rangle _\phi$ . The time series have been taken at mid-gap ( $r=r_{\textit{in}}+d/2$ ) position and the SRI time series also at mid-height position ( $z=H/2$ ). For the simulation, we used $ \textit{Re}^{-1}=2$ and $\tilde {l}=1.35$ . The scaling is given in the text.