# Quasi-biennial oscillation: laboratory experiments

Benoît Semin, François Pétrelis

► To cite this version:

Benoît Semin, François Pétrelis. Quasi-biennial oscillation: laboratory experiments. Comptes Rendus. Physique, 2024, 25 (S3), pp.1-25. 10.5802/crphys.195. hal-04767464

HAL Id: hal-04767464

https://hal.science/hal-04767464v1

Submitted on 5 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
Physique

Benoît Semin and François Pétrelis

**Quasi-biennial oscillation: laboratory experiments**

Published online: 30 September 2024

**Part of Special Issue:** Geophysical and astrophysical fluid dynamics in the laboratory

**Guest editors:** Stephan Fauve (Laboratoire de Physique de l'ENS, CNRS, PSL Research University, Sorbonne Université, Université Paris Cité, Paris, France) and Michael Le Bars (CNRS, Aix Marseille Univ, Centrale Marseille, IRPHE, Marseille, France)

https://doi.org/10.5802/crphys.195

<figure><img src="image_2.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Review article / *Article de synthèse*

Geophysical and astrophysical fluid dynamics in the laboratory /
*Dynamique des fluides géophysiques et astrophysiques au*
*laboratoire*

Quasi-biennial oscillation: laboratory
experiments

*Oscillation quasi-biennale : expériences de laboratoire*

Benoît Semin*<sup>*,a</sup>* and François Pétrelis<sup>b</sup>

<sup>a</sup> PMMH, CNRS, ESPCI Paris, PSL University, Sorbonne Université, Université Paris-Cité, F-75005, Paris, France

<sup>b</sup> LPENS, CNRS, ENS Paris, PSL University, Sorbonne Université, Université Paris-Cité, F-75005, Paris, France

*E-mails*: benoit.semin@espci.fr (B. Semin), petrelis@lps.ens.fr (F. Petrelis)

**Abstract.** The quasi-biennial oscillation (QBO) is an oscillation of the wind in the equatorial stratosphere. This wind is a mean flow induced by atmospheric waves, including internal gravity waves, which explain that the period (28 month) is not linked to any astrophysical forcing. This oscillation has only been reproduced in 3 laboratory experiments, which share a similar geometry. We present the details of our experimental set-up, and we explain which improvements allowed us to obtain quantitative measurements during long times. We show experimentally the feedback of the mean flow on the waves, which is one of the key ingredient of the oscillation. The details of the analytical resolution of the 1D model of Plumb and McEwan are given. We compare experimental results to analytical and numerical results, and found a qualitative agreement. The period decreases when the forcing increases, and the amplitude of the mean flow is not monotonic with respect to height and displays two local maxima as a function of height close to the threshold. The bifurcation is always a Hopf one, but can be subcritical or supercritical depending on the dominant dissipation mechanism of the mean flow which can be tuned experimentally by changing the Brunt–Väisälä frequency. We argue that an investigation of the bifurcation in general circulation models (GCM) is of interest to better understand the evolution of the QBO due to climate change.

**Résumé.** L'oscillation quasi-biennale (QBO) est une oscillation du vent dans la stratosphère équatoriale. Ce vent est un écoulement moyen engendré par des ondes atmosphériques, notamment des ondes internes de gravité, ce qui explique que la période (28 mois) n'est liée à aucun forçage astrophysique. Cette oscillation n'a été reproduite qu'à 3 reprises par des expériences analogues de laboratoire, et ces 3 expériences ont été réalisées dans la même géométrie. Nous présentons en détail notre dispositif expérimental, et en particulier les améliorations qui nous ont permis d'obtenir des mesures quantitatives sur de longues durées. Nous montrons expérimentalement la rétroaction du l'écoulement moyen sur les ondes, qui est l'un des ingrédients clés de la QBO. Nous donnons les détails de la résolution analytique du modèle 1D de Plumb et McEwan. Nous comparons les résultats expérimentaux à ceux obtenus analytiquement et numériquement, et nous

<footnote>*Corresponding author</footnote>
montrons qu'ils sont en accord qualitatif. La période de l'écoulement moyen diminue lorsque le forçage aug-
mente, et l'amplitude d'écoulement moyen en fonction de la hauteur n'est pas monotone et présente au voi-
sinage du seuil deux maxima locaux. La bifurcation est toujours une bifurcation de Hopf, mais elle peut être
sous-critique ou supercritique en fonction du mécanisme dominant de dissipation de l'écoulement moyen,
qui peut être changé expérimentalement en modifiant la fréquence de Brunt-Väisälä. Nous concluons en ex-
pliquant qu'une étude de la bifurcation dans les modèles de circulation générale (GCM) serait intéressante
pour mieux comprendre l'évolution de la QBO liée au changement climatique.

**Keywords.** Fluid mechanics, internal gravity waves, wave-mean flow interaction, instability, non-linear physics.

**Mots-clés.** Mécanique des fluides, ondes internes de gravité, interaction ondes-écoulement moyen, instabilité, physique non-linéaire.

Manuscript received 27 March 2024, revised 15 July 2024, accepted 5 August 2024.

# 1. Introduction

The quasi-biennial oscillation (QBO) is the periodic reversal of the wind in the equatorial strato-
sphere [1] (the stratosphere is the second layer of the atmosphere, between 16 km to 50 km at
the equator). The period of the oscillation is 28 months on average, and is not directly linked
to any astrophysical forcing. This wind is known to be generated by a broad spectrum of at-
mospheric waves, including small scale internal gravity waves and planetary waves like Kelvin
waves, Rossby-gravity waves and inertia-gravity waves [1]. The interaction between the wind and
the waves leads to the periodic reversal of the wind.

A key feature of the QBO is its downward propagation without loss of amplitude [1]. The
position where the mean flow changes sign drifts towards the location from where the waves are
emitted, i.e. the troposphere (the troposphere is the bottom layer of the atmosphere, between
0 km to 16 km at the equator).

The QBO has been mainly studied through observations on Earth, but similar wind oscillations
exist on other planets like Saturn [2] or Jupiter [3, 4]. Internal waves can also be sustained in stably
stratified layers in stars [5], and models of the Sun tachocline also observe periodic reversals of a
mean flow similarly to the QBO mechanism [6]. Overall this mechanism is robust and generic.

The QBO used to be one of the most predictable fluctuations of the large-scale circulation of
the Earth's atmosphere [7]. However, two disruptions of the QBO have been observed recently
in 2015/16 and 2019/20, after more than 50 years without observed disruption [8]. These
disruptions are thought to be due to strong extra-equatorial waves forcing.

A key interest of the QBO is the correlation between phenomena outside the tropical strato-
sphere, the so-called teleconnections [7, 9]. The QBO influences the convection in the tropo-
sphere close to the equator [10] and modulates the mixing and the transport of chemical con-
stituents in the atmosphere, including ozone among others [11]. Even if this wind oscillation is
located close to the equator, it has some influence in the extratropics. This wind oscillation is
correlated with hurricanes in North America [1] and affects winter conditions in Europe [12]. It
also influences the Madden-Julian oscillation [13, 14], an eastward propagation of large regions
of both enhanced and suppressed tropical rainfall, observed over the Indian and Pacific oceans.

The QBO has been simulated using general circulation models (GCM) for more than 20 years
now [15]. It used to be difficult for the models used for climate change forecast to produce a QBO.
For instance, most of the models included in CMIP5 (Coupled Model Intercomparison Project
Phase 5) do not reproduce the QBO [16]. Now, many climate models simulate QBO-like oscilla-
tions, but often with an amplitude weaker than the observed one in the lower stratosphere [7].
Moreover, these models rely on the parametrization of small-scale tropical waves which cannot
be resolved [7, 17], and the uncertainty on the response of the QBO to climate change is at least partly due to uncertainties on internal gravity waves [18].

Simple 1D models have been developed when the basics mechanisms of the QBO have been understood [19, 20]. These models have inspired the analog experiments of the QBO (see [21] and discussion below), and reciprocally the experiments have shown the relevance of the physical ingredients used in these simple models. We will discuss this aspect later in the present article. These simple models are still useful to extract some basics mechanisms, as for example in [22] where the successive bifurcations of the Plumb and McEwan model are investigated and where the existence of a second bifurcation is shown, leading to quasi-periodic trajectories. In [23] the effect of the boundary condition is investigated. The excitation with a spectrum of waves has been studied in [24], where the authors show that the spreading of the wave energy across a wide frequency range leads to more regular oscillations.

Couston et al. 2018 [25] have shown that an oscillating flow can emerge from internal gravity waves generated by turbulent convective motion, using direct numerical simulations (DNS), i.e. without parametrization. A 2D model coupling the troposphere and the stratosphere and taking into account the phase change of water has been used to investigate the effect of the QBO on precipitations [26]. A comparison between 2D and 3D direct numerical simulations reproducing the experiments has been performed for different forcing parameters [27].

Three analog experiments of the QBO have been performed: the first has been performed by Plumb and McEwan in 1978 [21], the second by Otobe et al. in 1998 [28] and the third by Semin et al. in 2018 [29]. These 3 experiments share a similar geometry: a linearly stratified aqueous solution is placed between two concentric vertical cylinders, of size about 0.5 m. The forcing is achieved using 16 membranes placed at the bottom [21] or top [28, 29] of the fluid. Each membrane oscillates in opposite phase with its two neighbors, so that the forcing is a pattern of azimuthal standing waves, which is equivalent to two counter-propagating waves. The choice of the geometry, size and number of membranes has thus been done by Plumb and McEwan in 1978 [21].

In the first experiment by Plumb and McEwan [21], the flow was measured using tracers. Since this study dates back before particle image velocimetry (PIV) became widely available, the particles were tracked manually from "Cine film". The average velocity on a cycle was measured. Tracers generated by electrolysis on a vertical wire were also used. This first experiment has validated the wave/mean flow theory, the instability of mean flow, the reversal of the mean flow and the origin of the downward propagation. This was an essential validation of a mechanism which is far from being obvious. A stability diagram is displayed in this article, which shows a very qualitative agreement with the theory. Only measurements over one period are presented.

The second experiment by Otobe et al. [28] has confirmed the results of Plumb and McEwan [21]. The mean flow is measured with particles, and the waves with a Moiré technique. The distortion of the waves by the mean flow is presented. The difficulty to obtain reproducible results is emphasized: the authors state for example "The period of the oscillation is 1 to 2 hours, but there is no apparent correlation with the experimental parameters".

The third experiment has been done by Semin et al. [29]. We will present this experiment in more detail in the present article. This is a quantitative experiment, where both the waves and the mean flow are measured using PIV. Continuous slow injection of new (and heavier) fluid at the bottom of the tank and slow withdrawal of fluid at the top prevented the weakening of the stratification at the top where the membranes move, and thus allows very long measurements. We have determined the nature of the bifurcation. The bifurcation is always a Hopf bifurcation, but can be supercritical or subcritical depending on the main dissipation of the mean flow.

An experiment close to the analog QBO experiments has been performed by Léard et al. [30] (see also [31]). The authors take advantage of the non-monotonic change of the density of water
around 4 °C to excite directly the waves in a stratified layer by the convection in the lower layer. A mean flow appears in the upper stratified layer, but with a dynamics different from the QBO, explained by the authors by the difference in Prandtl numbers between water and air. Another experiment has been reported by Léard et al. 2021 [32]. The waves in a stratified fluid are forced by 12 jets in a upper unstratified layer. A mean flow of very low intensity is obtained in this set-up where the forcing is turbulent and statistically symmetric, contrary to the previously reported experiments of mean flow generation by internal waves [33, 34]. However, no reversal has been obtained, probably due at least partly to the mixing at the interface between the two layers. These experiments justify a posteriori the use of membranes to force the fluid in the analog QBO experiments.

In this article, we present in Section 2 the details of our experimental set-up, including many tricks not discussed in our previous articles in [29, 35]. The details of the analytical resolution of the Plumb and McEwan model is given in Section 3 and in Appendix B. Numerical simulations are also discussed. The experimental results are presented in Section 4.

## 2. Experimental set-up

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic view of the experimental set-up.</figcaption></figure>

A schematic view of the experimental set-up is shown in Figure 1 and a picture in Figure 2 left. It has been already described in [35]. It is made of two transparent concentric vertical cylinders. The outer diameter of the inner cylinder is 365 mm, the inner diameter of the outer cylinder is 600 mm, and the height is $H = 410$ mm.

The gap $h$ between these two cylinders is filled with a linearly density-stratified solution of NaCl or $MgCl_2$ in water, obtained using the double-bucket method [28]. We use NaCl for water softener (Laugil) and magnesium chloride ice melt (France Déneigement). For NaCl, we use 62 L of tap water in bucket 1, and 48.5 L of tap water and 15 kg of NaCl in bucket 2. For $MgCl_2$, we use 55.9 L of tap water in bucket 1, and 19.5 L of tap water and 36.4 kg of $MgCl_2$ in bucket 2 ($MgCl_2$ is already partially hydrated when we buy it). Moreover 40 mL of chlorine bleach is added in each bucket to avoid algea proliferation. During the injection, bucket 1 is constantly mixed using a VWR VOS overhead stirrer, as required by the double-bucket method. The density profile is obtained from conductivity, measured using a conductivity probe (Mettler-Toledo SevenGo Pro) positioned at different heights using a linear manual translation stage. We have used a pycnometer to measure the density of solution that we used to calibrate the conductivity.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left: picture of the experimental set-up. The green arrows indicate the motion of the membranes, and the red and blue arrow the direction of the mean flow. Right: picture of a membrane.</figcaption></figure>

Note that the conductivity of MgCl₂ is not monotonic with respect to density: a given value of conductivity may correspond to two values of density. We use the stably stratified property of the fluid to select the right value. The stratification is quantified by the Brunt–Väisälä frequency, $N = \sqrt{-(g/\rho_0)d\rho_0/dz}$, where g is the acceleration of gravity, $\rho_0$ is the background density at the top and z is the vertical axis. The typical value of N is 1.5–2.2 rad·s⁻¹, the highest values are obtained using MgCl₂.

The fluid motion is forced using 16 silicone membranes (Profiles Market, thickness 0.3 mm). The membranes are in contact with the top of the fluid and are fixed using nylon and stainless steel screws on a laser-cut PMMA annulus (see Figure 2 right). We used Boehm hollow punches to make the fixation holes in the membranes. However, it happened several times that a membrane broke due to tearing, which then induced mixing at the top and unreliable experimental results, which were then discarded. The silicone membranes do not degrade in salty water. A gap of a few millimeters is located between the annulus and the cylinders, so that air can be removed from the membrane by pushing gently with a plastic stick.

Each membrane can be individually moved up and down by a linear stepper motor (Nanotec L2818-OV15601 with SMCI36 controller), which is computer controlled using Labview. The membranes are fixed to the motor stick using two rubber plates (typical dimensions 45 mm, in black in Figure 2 right) which are sandwiched between two transparent and rigid PMMA plates. The motion of each motor is a sine function, and two neighboring membranes are driven by motors in opposition of phase so that the forcing is a standing wave whose azimuthal wavelength $\lambda_x \approx 200$ mm is twice the curvilinear distance between 2 motors. We note M the amplitude of the motion of the membranes and we report here on measurements performed with a forcing period $T_f = 15$ s. The length of the stick of the motor number m is:

$$h(t) = M \sin(\omega t + m\pi) + h_{0,m} \quad (1)$$

where $\omega = 2\pi/T_f$ and $h_{0,m}$ is the initial value of the stick length.
The fluid velocity is measured using particle image velocimetry (PIV). Since the velocities are slow, we use a low-cost laser sheet (Apinex, 5 mW) and an industrial camera (Allied Vision Stingray) at 10 or 12 frames per second. The particles are imaged in the whole height on a width of about 30 mm. Movies of 900 images are acquired every 100 s or 1500 images every 250 s. We made the seeding particles using wax and alumina powder, to obtain particles with a range of average densities similar to the range of densities in the fluid. The protocol to make the beads is detailed in Appendix A. The beads are mixed with hand washing liquid and water (50%) before being gently poured at the top of the fluid, to avoid air bubble accumulation at their surface. We use OpenPIV [36] on Matlab to determine the velocity field from two successive images. To take into account the small lens effect due to the cylinder, a different calibration length is used on the vertical and horizontal axis. The relaxation time of the particle is $t_b = \rho_b(2r_b)^2/(18\eta)$, where $\rho_b$ is the density of the particle (slightly above that of water), $r_b \sim 0.5$ mm is the radius of the particle and $\eta$ the dynamic viscosity of water. An estimate of this time is thus $t_b \sim 0.06$ s $\ll T_f$: despite their large size, the particles are expected to follow the flow. The typical interrogation windows has a width of 64 pixels and a height of 32 pixels (corresponding to 11 mm), and the overlap between the interrogation windows is 50%. In each interrogation window, the velocity as a function of time is fitted using a sine function to obtain the amplitude and the phase of the wave (oscillation at the forcing angular frequency $\omega$) and the mean flow. We have checked that the mean flow is azimuthal and only depends on the height $z$, except very close to the boundaries.

One of the main improvements between our set-up and the previous ones is the continuous injection of fluid at the bottom and removal at the top using a peristaltic pump (Ismatec Reglo ICC 4 channels). Even if the membrane motion at the top tends to slowly mix the fluid, this ensures that the stratification remains linear (see Figure 3). No measurement of mean flow lasting more than 2 periods had been reported before our experiments. In the first experiments that we performed, without pump, we also observed only 1 or 2 oscillations of the mean flow, and the measured stratification at the top was very weak at the end of the experiment. Since the stratification affects the emission and propagation of waves in the system, the weakening of the stratification is consistent with the impossibility to generate a mean flow any more. To avoid such a weakening of the stratification due to mixing, salty fluid is continuously injected at the bottom (around 9 mL·min⁻¹) and removed at the top of the experiment in 3 points: between the annulus of membranes and the cylinder (around 21 mL·min⁻¹) and on top of two membranes (around 14 mL·min⁻¹). This configuration ensures that no liquid accumulates at the top. The typical vertical velocity associated with this flow is 1 µm·s⁻¹, which is significantly smaller than the other velocities involved in the experiment. This modification allows us to perform quantitative measurements in a nearly steady state (variation of N less than 0.1 rad·s⁻¹ during a run of 90 000 s, see Figure 3). However, since the saturated concentration of the salt is finite, this system imposes a limit on the duration of the experiment. When the saturated fluid is injected, the experiment has to be stopped.

# 3. Theoretical approach

## 3.1. *Description of the model*

The aim of this section is to present the model of [21, 35], which is a minimal model of the QBO. In the next section, we present several analytical results obtained on this model. Some of them were mentioned in the supplementary material of [29].

In this model, a mean flow is generated by two counter-propagating internal gravity waves. These two waves propagate through a linearly stratified fluid, characterized by a Brunt-Väisälä frequency, $N = \sqrt{-(g/\rho_0)\frac{d\rho_0}{dz}}$, where g is the acceleration of gravity, $\rho_0$ is the background
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Stratification before (red dots) and after (blue diamonds) a 90000 s long experiment. The orange dashed line is a linear fit of the red dots, and the cyan dotted line is a linear fit of the blue diamonds. Stratification using NaCl, $T = 15$ s, $M = 13.5$ mm. The Brunt-Väisälä frequency computed from the linear fits is $N = 1.48 \text{ rad} \cdot \text{s}^{-1}$ before the experiment, and $N = 1.40 \text{ rad} \cdot \text{s}^{-1}$ after.</figcaption></figure>

density and $z$ is the upwards vertical axis. The angular frequency of the waves is $\omega$, one of the wave is propagating in the azimuthal direction $+x$, and the other in the direction $-x$; the horizontal wavevector is written $k_x$. These waves verify the linear dispersion relation for internal gravity waves:

$$ \omega^2 = N^2(\sin(\theta))^2, \qquad (2) $$

where $\theta$ is the angle between the wavevector $k$ and the vertical. Due to the structure of the internal gravity waves, the horizontal wave component $u'$ and the vertical component $v'$ verify:

$$ v' = u' \tan(\theta), \qquad (3) $$

The equation of the model are obtained using the Navier-Stokes equation, the Boussinesq approximation and the Wentzel-Kramers-Brillouin approximation. As in the original work of Plumb and McEwan [21], we consider upward propagating waves in the calculation. The following dimensionless equation is obtained for the mean flow $u$:

$$ \frac{\partial u}{\partial t} = -\frac{\partial D}{\partial z} + \Lambda_1 \frac{\partial^2 u}{\partial z^2} - \Lambda_2 u, \qquad (4) $$

where:

$$ D = F \left( c_1 \exp \left[ - \int_0^z \left( \frac{\alpha_1}{(u-c_1)^2} + \frac{\alpha_2}{(u-c_1)^4} \right) dz' \right] + c_2 \exp \left[ - \int_0^z \left( \frac{\alpha_1}{(u-c_2)^2} + \frac{\alpha_2}{(u-c_2)^4} \right) dz' \right] \right), \quad (5) $$

traces back to the forcing of the mean flow by the waves. $v$ is the kinematic viscosity, $c_1 = 1$ and $c_2 = -1$, $\Lambda_1, \Lambda_2, \alpha_1, \alpha_2$ are positive, $\alpha_1 + \alpha_2 = 1$ and $F$ is the control parameter. An experiment is characterised entirely by 3 dimensionless parameters, for instance $\Lambda_1, \Lambda_2$, and $\alpha_1$. However, we keep $F$ as a free parameter here because we will expand $F$ close to the threshold to apply the perturbation method. The boundary conditions are:

$$ u(0) = 0, \qquad (6) $$

$$ u(+\infty) = 0, \qquad (7) $$

and corresponds to a no slip boundary condition at $z = 0$ like in [21] and a vanishing velocity at infinity. We chose the latter condition to simplify the analysis and because the mean flow is actually localised close to the forcing boundary.
To make the link between the dimensionless variables and the physical variables, we note in the following set of equation and only in this set the dimensional variables with a ~:

$$
\begin{align*}
\tilde{F} &= \langle \tilde{u}' \tilde{v}' \rangle, \\
\tilde{c} &= \frac{\tilde{\omega}}{\tilde{k}_x}, \\
\tilde{d} &= \left( \frac{\tilde{N}\tilde{\gamma}}{\tilde{k}_x \tilde{c}^2} + \frac{\tilde{N}^3 \tilde{v}}{\tilde{k}_x \tilde{c}^4} \right)^{-1}, \\
\tilde{z} &= \frac{\tilde{z}}{\tilde{d}}, \\
u &= \frac{\tilde{u}}{\tilde{c}}, \\
\alpha_2 &= \frac{\tilde{N}^3 \tilde{v} \tilde{d}}{\tilde{k}_x \tilde{c}^4}, \\
\Lambda_1 &= \frac{\tilde{v} \tilde{c}}{\tilde{F}(0) \tilde{d}}, \\
\Lambda_2 &= \frac{\tilde{\gamma} \tilde{c} \tilde{d}}{\tilde{F}(0)},
\end{align*}
$$

where $\tilde{F}$ the momentum flux, $\tilde{c}$ the horizontal phase velocity, $\tilde{d}$ the vertical dissipation length, $\tilde{v}$ is the kinematic viscosity, $\tilde{\gamma}$ the dissipation due to the walls, $\tilde{F}(0)$ the momentum flux at $\tilde{z} = 0$.

From the dimensionless solution of equation (4), it is possible to obtain the dimensional variables: the dimensionless lengths should be multiplied by $\tilde{d}$ to obtain the dimensional lengths, the dimensionless times should be multiplied by $\tilde{c}\tilde{d}/\tilde{F}(0)$ and the dimensionless mean flow by $\tilde{c}$.

## 3.2. Linear and weakly non-linear solutions of the equations

For the dimensionless system of equations (4) and (5), with the boundary conditions (6) and (7), we have obtained an explicit solution of the linear problem from which an equation for the threshold is obtained. We have also derived the amplitude equation in the vicinity of the threshold. Details of the calculation are given in Appendix B and we sum up the results here.

We introduce the stream function $\psi$

$$ \psi = \int_{0}^{z} u \, dz \qquad (8) $$

so that

$$ u = \frac{\partial \psi}{\partial z} \qquad (9) $$

We are considering eigenmodes such that:

$$ \frac{\partial \psi}{\partial t} = \mu \psi \qquad (10) $$

with $\mu \in \mathbb{C}$.

We define the parameters

$$ a = \frac{F \times 2(2\alpha_1 + 4\alpha_2)}{\Lambda_1}, \qquad (11) $$

$$ b = \frac{\mu + \Lambda_2}{\Lambda_1}, \qquad (12) $$

that appear in the expression of $\psi$

$$ \psi(z) = \psi_0(z) (1 + K f(z)). \qquad (13) $$
$\psi_0$, K and f are given by

$$ \psi_0 = J(2\sqrt{b}, 2\sqrt{a}\exp(-z/2)), \quad (14) $$

$$ K = \frac{-1}{f(0)} = \frac{-1}{\int_{z_0}^{0} \frac{\int_{+\infty}^{z} \psi_0(z') dz'}{(\psi_0(z))^2} dz}. \quad (15) $$

$$ f(z) = \int_{z_0}^{z} \frac{dz' \int_{+\infty}^{z'} \psi_0(z'') dz''}{(\psi_0(z'))^2}, \quad (16) $$

where J is the Bessel function of the first kind, of order $2\sqrt{b}$, and of argument $2\sqrt{a}\exp(-z/2)$. $z_0$ is an arbitrary positive constant; the final solution does not depend on the value chosen for $z_0$.

Using the boundary conditions and properties of the Bessel functions, we obtain an integral equation that relates the parameters of the problem and the growth rate

$$ \int_0^1 \left[ f(2\sqrt{b}-1, 2\sqrt{a}v) + f(2\sqrt{b}+1, 2\sqrt{a}v) \right] dv = 0. \quad (17) $$

We have solved equation (17) numerically using a Matlab minimization routine and Bessel functions of complex order and argument [37]. We fix the value of $\alpha_1$ (and thus $\alpha_2 = 1 - \alpha_1$). We set $F = 1$ which amounts to changing $\Lambda_i$ into $\tilde{\Lambda}_i = \Lambda_i/F$ and dropping the $\tilde{\cdot}$ to ease notation. We then compute $a$ from equation (11). Using equation (17), we obtain the corresponding value for $b$ at the threshold. Since $\mu$ is imaginary at the threshold and $\Lambda_2$ is real, we obtain from equation (12) the value of $\Lambda_2 = \Lambda_1 \Re(b)$ and $\mu = \Lambda_1 \Im(b)$. The eigenmodes are obtained using equation (13) with equations (16), (15) and (14), and their time evolution are given by multiplying the eigenmodes by $\exp(\mu t)$.

We have also derived the amplitude equation for the unstable mode in the vicinity of its onset of instability. We note

$$ F = F_0 + \varepsilon F_1 \quad (18) $$

where $F_0$ is the value at the instability threshold ($\Lambda_1$ and $\Lambda_2$ being fixed) and $\varepsilon$ a small parameter (dimensionless distance to the threshold).

We introduce a slow timescale $T = \varepsilon t$ and write $u$ as:

$$ u = \sqrt{\varepsilon} (u_0 + \varepsilon u_1 + \varepsilon^2 u_2 + \dots) \quad (19) $$

At order $\varepsilon^{1/2}$, we recover the linear problem so that $u_0$ is the product of the eigenmode of the linear problem $u$ with a slowly varying amplitude:

$$ u_0(t, T, z) = \frac{1}{2} (A(T)u(t, z) + A^*(T)u^*(t, z)), \quad (20) $$

where $u$ is given by equation (9) and (13).

At order $\varepsilon^{3/2}$, using the solvability condition, we obtain the amplitude equation for $A$

$$ \frac{\partial A}{\partial T} = \alpha A F_1 + \beta |A|^2 A, \quad (21) $$

where explicit expressions for $\alpha$ and $\beta$ are given in the annex.

This is the amplitude equation of a pitchfork Hopf bifurcation. Indeed, this is constrained by the oscillatory nature of the unstable mode that breaks invariance under shift in the origin of time.

The coefficients can be calculated numerically using the expressions given in Appendix B. In particular, we have calculated the following expression using a Matlab routine

$$ S = \frac{\Re(\beta)}{\Re(\alpha)}, \quad (22) $$

where $\Re$ denotes the real part. The nature of the bifurcation is determined by the sign of $S$: if $S > 0$, the bifurcation is subcritical, and it is supercritical if $S < 0$.
The value of S is determined by a and b through equation (14) and the equations mentioned after. Since $\mu \in i\mathbb{R}$ at the threshold, the value of b at the threshold only depends on the mean flow dissipation through $\Lambda_1$ and $\Lambda_2$ (see equation (12)). The wave dissipation, i.e. the values of $\alpha_1$ and $\alpha_2$, only changes the value at the threshold of F. The nature of the bifurcation thus depends on the dissipation of the mean flow, whereas the nature of the dissipation of the waves only changes a prefactor which has always the same sign.

### 3.3. Numerical simulations

We have solved numerically the system of equations (4) and (5) using a centered second-order finite difference method with grid size $\delta z = 10^{-2}$ (or $\delta z = 4 \times 10^{-3}$ for the stability diagram). The total length is 4 in dimensionless unit, i.e. $4d$ in physical units: this is sufficient to mimic an infinite domain, since the waves are almost entirely damped for $z = 4$. Since the optimization of the numerical simulations is not the core of the present article, we have used an explicit Euler method in time with $\delta t = 10^{-4}$ (or $\delta t = 10^{-5}$ for the stability diagram), which is simpler than the schemes used by [38] and [25]. The boundary conditions are no-slip at both the bottom and upper boundary, according to equations (6) and (7). This numerical scheme is implemented in Matlab. The equations and the boundary conditions are not exactly the same as the ones in the works of [38] and [25], who only consider bulk dissipation of the mean flow an a stress-free boundary condition at the upper boundary.

The numerical system can be solved using explicit Euler method, is numerically rather stable, and rapid to solve using a laptop (a few minutes with a 2 GHz processor laptop, using a Matlab script). We point out that this system can be used as a non-trivial example of numerical resolution at the undergraduate level.

## 4. Experimental results

### 4.1. Interaction between wave and mean flow

The mechanism of the QBO is based on the interaction between the waves and the mean flow. We investigate this interaction in this section.

The characteristic of the waves and the mean flow 2 min after the beginning of the forcing at $M = 14.5$ mm is shown in Figures 4 (a,b,c). The waves are established, but the mean flow is still zero. This is in agreement with the study of a single wave in the same set-up [35], and confirms that the characteristic time-scale of the mean flow is much larger than the one of the waves.

The initial phase of the wave is shown in Figures 4 (b). This phase corresponds to the expected stationary wave in the x direction. The measured vertical wavelength $\lambda_z = 46$ mm is in reasonable agreement with the expected value $\lambda_z = \lambda_x \tan(\theta) \approx 40$ mm. The local amplitude of the wave is given by Figure 4 (a). The amplitude of the wave is not constant in the x direction, as expected for the sum of two counter-propagative waves. The wave amplitude does not vanish, which is linked to the choice of the position of the field of view, and to the fact that it is shorter than half the wavelength in the x direction $\lambda_x/2 = 100$ mm. By averaging the amplitude in the x direction, we obtain Figure 5. The amplitude of the wave decreases with height. This is due to the dissipation of the wave, mainly due to the dissipation by viscosity in the bulk: $\alpha_2 = 0.965$ where $\alpha_2$ is the dimensionless quantification of the bulk dissipation whereas $\alpha_1 = 0.035$ is the dissipation due to the wall. We can fit the curve with the function:

$$u = \mu_0 \exp(-z/(2d_e)) \quad (23)$$

where $d_e$ is the experimental decay length and $\mu_0$ is related to the amplitude of the wave at the top through $u_0 = \mu_0/q$, where $q = 1.6$ in our experiment takes into account that the measured
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Measurements of the wave and mean flow for an experiment with $M = 14.5$ mm with $\text{MgCl}_2$ ($N = 2.16 \text{ rad} \cdot \text{s}^{-1}$). The first row of the figure corresponds to the beginning of the forcing, when the mean flow is still small (sub-figures a,b,c). The second row corresponds to a time at which the mean flow is established. We chose a time at which the mean flow is negative everywhere for the simplicity of the interpretation (sub-figures d,e,f). The first column corresponds to the root-mean square value of the horizontal wave velocity $\langle\sqrt{\langle u' \rangle^2}\rangle$ (sub-figures a,d). The second column corresponds to the wave phase, obtained from the horizontal wave velocity $u'$ (sub-figures b,e). The third column corresponds to the mean flow (sub-figures c,f).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Red dots: wave root-mean-square value averaged in the *x* direction, from Figure 4 (a), i.e. *M* = 14.5 mm with MgCl₂, no mean flow. Dashed line: exponential fit, with *μ*₀ = 5.45 mm·s⁻¹, *d*₀ = 81 mm.</figcaption></figure>

amplitude is linked to the sum of two counter-propagating waves, measured in a finite field of view. The fit is good, indicating that the model is satisfactory. The value of *d*₀ from the fit is *d*₀ = 81 mm, is in rather good agreement with the theoretical value *d* = 95 mm, obtained using *γ* = 10⁻³ s (see [35]), *ν* = 10⁻⁶ m²·s⁻¹ and the measured value of *N*. The small discrepancy is probably due to the difference in viscosity, the solution of MgCl₂ being more viscous than water or NaCl solution. We did not perform a quantitative analysis because the concentration and thus the viscosity is not homogeneous in the solution.

We now consider the same forcing parameters but after 75 min in Figures 4 (d,e,f). The mean flow is well established at that time. We have chosen a snapshot at which the mean flow has the same sign everywhere so that the qualitative explanations are easier. The phase displayed in Figure 4 (d) is different from the one when there is no mean flow (Figure 4 (b): there is a kink in the phase at the top. At the bottom, the phase is tilted with an angle 10°, where the theoretical angle for a single wave is $\theta = \arcsin(\omega/N) = 11°$. This figure shows that there is a feedback of the mean flow on the wave, and that one of the wave is dominant at the bottom, the other being of much smaller amplitude. The wave at the bottom has a positive horizontal wave velocity (see for instance [39] for the basic geometry of internal gravity waves). This feedback is consistent with the denominator in equation (5), and consistent with the qualitative explanation of the QBO [1] for which the wave which has a horizontal phase velocity of the same sign as the mean flow accelerates the mean flow more efficiently than the wave which has a horizontal phase velocity of the opposite sign, and is thus more attenuated. The modification of the waves by the mean flow has been shown in [28] using the Moiré technique and numerically in [27], but our method is more direct because we extract the wave phase. The feedback of the mean flow on the wave is also seen by comparing the amplitude of the wave without and with a mean flow in Figures 4 (a) and (d). In the latter case, the amplitude is lower everywhere but very close to the membrane. This is consistent with an increased attenuation of one of the waves due to the mean flow. The feedback of the mean flow on the waves for a single wave has been already shown experimentally in the geometry of a dipole and a wave beam in [40] and in a similar geometry in our previous work [35]. It has also been shown with two waves in [28], but with a more qualitative measurement of the waves.

## 4.2. Characteristics of the mean flow

An example of a space-time diagram of the mean flow, known as Hovmöller diagram in the atmospheric sciences, is shown in Figure 6. In this example, the mean flow is zero during a very long time and then suddenly appears without any apparent cause. This is consistent with
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Example of mean flow as a function of time for a forcing amplitude *M* = 13.5 mm, a period *T*<sub>f</sub> = 15 s and *N* = 1.5 rad·s<sup>-1</sup>, stratification using NaCl.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Mean flow as a function of time at *z* = 320 mm for the same experiment as the one shown in Figure 6: *M* = 13.5 mm, *T*<sub>f</sub> = 15 s, *N* = 1.5 rad·s<sup>-1</sup>, NaCl. Blue points: experimental data, blue dashed line: curve linking the experimental points, black line: sine fit.</figcaption></figure>

a subcritical transition to the QBO for this value of *N*, as discussed in [29] and in Section 4.3.
More than 9 periods of the mean flow are observed in this example, which is much higher than
in any previous experiments [21, 28]. As mentioned in Section 2, this has been achieved by the
continuous injection of fluid at a very small velocity using an peristaltic pump.

The value of the mean flow at a given height for the same experiment is shown in Figure 7.
This figure shows that the mean flow is actually quite regular, and can be fitted at a given height
by a sine function during more than 9 periods. This sinusoidal behaviour of the mean flow is
consistent with the weakly nonlinear methods used in the theoretical approach of Section 3.

Since the mean flow is sinusoidal in time, it is possible to plot the relative phase of the mean
flow as a function of height, as done in Figure 8 left. The phase difference between the top and
the bottom is close to $\pi$ in the field of view. This is consistent with the usual view of the QBO and
the theoretical calculations where the mean flow at the top and at the bottom are in opposition
of phase.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left: relative phase of the mean flow, $M = 12$ mm, $\text{MgCl}_2$. The horizontal dashed line corresponds to $\pi$. Right: amplitude of the mean flow as a function of height $z$. Blue dots: experiment, $M = 12$ mm (threshold at $M = 11.3$ mm), $\text{MgCl}_2$; magenta line and diamond symbols: numerical resolution of the 1D model, $M = 10.0$ mm (threshold at $M = 9.1$ mm), orange line: analytical linear profile, height coefficient chosen to match the data.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Left: period of the mean flow as a function of the forcing amplitude $M$, for $\text{MgCl}_2$. Blue points: experiment, the forcing is changed every 8000 s to reach a quasi-steady state, and the period is measured in the last 4000 s. See [29] for the space-time diagram. Magenta line and diamond symbols: numerical simulation. Right: amplitude as a function of $M$. Same color code. Height $z = 320$ mm for the experimental data, i.e. about $d_e$ from the forcing, and distance $d$ for the simulations. The blue dashed line and the violet dashed dotted line are of the form $\propto \sqrt{M-M_t}$, where $M_t$ is the threshold, which is different for the experiments and the simulations.</figcaption></figure>

The amplitude of the mean flow as a function of height is displayed in Figure 8 right. Experimental results, numerical and analytical results are shown in this figure. The chosen parameters are slightly different, due to the difference in the thresholds between the experiment and the numerical simulation. The value of $d$ for scaling the simulation is taken from Figure 5, and the value of $F$ is extrapolated from the value of $\mu_0$ obtained from this figure, and by using the fact that $F$ is quadratic in $M$. Since the wave velocity varies by about 20% across the width [35], and that the previous hypotheses may not be perfectly valid, we cannot expect a perfect quantitative
agreement. However, Figure 8 right shows that the values of the velocities and the shape of the mean flow profile are in qualitative agreement. In particular, the variation of the mean flow with height is not monotonic, leading to a local minimum in the amplitude. Experimentally, this effect is robust and has been observed irrespective of the stratification. Numerically, it is observed provided that the distance to the threshold is not too large. This shape is also observed in the linear analytical shape.

The period as a function of the forcing *M* is shown in Figure 9 left. The set of experiments used includes the one used in Figure 8, and the numerical simulations are performed using the same parameters. The period is finite at the threshold and decreases when *M* increases. The experiments and the simulations are in qualitative agreement. The period of the mean flow, about *T*<sub>QBO</sub> = 2500 s is much larger than the wave period *T*<sub>f</sub> = 15 s as expected. Note that obtaining a reasonable value of the mean flow period is not easy experimentally. It is achieved successfully here thanks to the continuous injection of new fluid, as shown by the failure of [28] to obtain meaningful values for this period.

The mean flow amplitude as a function of the forcing *M* is shown in Figure 9 right. The value of the thresholds between experiment and theory are different but close. In both cases, the transition is continuous and the first point can be fitted by a form $\propto \sqrt{M-M_t}$. This behavior and the finite value of the period at the threshold shows that the bifurcation is a supercritical Hopf bifurcation in that experiment. The nature of the bifurcation will be investigated in more details in the next section.

4.3. *Stability diagram*

The stability diagram is shown in Figure 10. The bifurcation is subcritical at small values of Λ₁, and supercritical for large values of Λ₁. The tricritical point is at about Λ₁ = 0.12 and Λ₂ = 1.88. The values of the thresholds are in agreement to within 3% between the numerical and analytical solutions. The nature of the bifurcation is roughly in agreement. The behaviour of the bifurcation for the two different salts depends on the ratio Λ₂/Λ₁ = γd²/v. This ratio depends on N, and in our experiments of the salt because we use the more soluble and more expensive MgCl₂ only to achieve high values of N. The two lines displayed in Figure 10 correspond to two experiments, one with NaCl and a low value of N and one with MgCl₂ and a larger value of N. During an experiment with MgCl₂ where the amplitude of the forcing is changed, the ratio Λ₂/Λ₁ changed slightly due to the change in the stratification (see Section 2) and to the fact that the kinematic viscosity of a MgCl₂ solution depends on its density. The ratio γ/v is independent of the viscosity, but d decreases when the viscosity is increased. More precisely, it is possible to evaluate the value of Λᵢ for the point where there is no mean flow, because the wave amplitude can be fitted by an exponential (see an example in Figure 5). This is the case for the largest forcing amplitude (plain diamonds), for which the mean flow is initially zero. It is also possible for the points at which there is no more mean flow (empty diamonds). In between, it is possible to interpolate the values of Λᵢ for the other values of the forcing for which there is a mean flow, assuming that F(0) is quadratic in the forcing amplitude and that the viscosity increases linearly with time (the latter hypothesis is reasonable because the change is small). The dashed lines in Figure 10 correspond to a guide to the eye of where the experimental points should be without viscosity effects.

For MgCl₂, we have a rather good agreement between the experimental threshold and the theoretical one, without any adjustable parameter. This is not the case with NaCl. The reason is not clear, but we can note that the decay length *d* for NaCl is different from the calculated one.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Stability diagram, from analytical solutions ($\alpha_2 = 1$), simulations, with experimental values. The solution is stable (i.e. no mean flow) for large values of $\Lambda_1$ and $\Lambda_2$, and unstable (oscillating mean flow) for small values of $\Lambda_1$ and $\Lambda_2$. Red line and red circles: analytical values for the threshold, supercritical bifurcation; black circle: analytical value of the tricritical point; blue line and blue circles: analytical values for the linear threshold, subcritical bifurcation. Red squares: threshold from numerical simulations, case where the bifurcation is supercritical; blue squares: threshold from numerical simulations, case where the bifurcation is subcritical. Dashed orange line: line corresponding approximately to $MgCl_2$ during an experiment where the wave flux is changed; dashed violet line: same for NaCl. Plain orange diamond: experimental values where a mean flow is observed for $MgCl_2$; empty orange diamonds: same but without mean flow; orange asterisk: experimental values for which a mean flow is observed, with value of forcing interpolated. Violet diamonds: same as orange diamonds, but for NaCl.</figcaption></figure>

**5. Conclusion**

The quasi-biennial oscillation (QBO) is a the periodic reversal of the wind in the equatorial stratosphere, due to the forcing by atmospheric waves. The period, of about 28 months, is not linked to any external astrophysical forcing, which is a remarkable feature [41].

Three analog experiments have been performed in 1978 by Plumb and McEwan [21], in 1998 by Otobe et al. [28] and in 2018 by Semin et al. [29]. The geometry is similar in all three experiments: a linearly stratified aqueous solution between two cylinders is forced using 16 membranes. Simple 1D model [20] and a careful consideration of the dimensionless numbers of the problem have lead to the choice of this particular set-up.

The first experiment [21] has provided a proof of the concept of the wave induced mean flow oscillation. The last experiment [29] has shown that the bifurcation between a state without mean flow and a state with a mean flow depends experimentally on the value of the Brunt–Väisälä frequency. The model shows that the nature of the bifurcation depends on the main dissipation mechanism of the mean flow, i.e. in the bulk or in the boundary layer close to the wall of the cylinders. This main dissipation mechanism changes in our experiments when the Brunt–Väisälä frequency changes.

In this article, we have presented the details of the experimental set-up discussed in [29, 35], including many experimental tricks.

We also give the details of the analytical resolution of the 1D model of Plumb and McEwan [21]. This model describe the generation of a mean flow by two counter-propagating internal gravity waves. We have presented an analytical solution of this model up to the third order, in order
to determine the threshold, the eigenvectors and the nature of the bifurcation. Two dissipation mechanisms are taken into account for the waves and the mean flow: one linked to the effect of viscosity in the bulk, described by the dimensionless coefficient $\Lambda_1$ for the mean flow and $\alpha_2$ for the waves, and one for the dissipation at the lateral walls, described by the dimensionless coefficient $\Lambda_2$ for the mean flow and $\alpha_1$ for the waves. We have solved analytically this model, and we have obtained the value of the threshold as a function of the dimensionless parameters. The nature of the dissipation of the waves changes only quantitatively the threshold. In contrast, the nature of the dissipation of the mean flow plays a major role in the bifurcation: the bifurcation is supercritical for large values of $\Lambda_1$ and subcritical for small values of $\Lambda_1$.

It is remarkable that the model of Plumb and McEwan can be solved analytically. However, as shown in this article, the resolution is long. On the other hand, the numerical resolution of this model is simple and rapid: we think that it would be a good example for advanced undergraduate students.

We have measured simultaneously the waves and the mean flow using PIV. We have shown an example of feedback of the mean flow on the waves. The experimental results are consistent with the mechanism of the QBO: the attenuation of a wave (and its momentum transfer to the mean flow) depends on the relative sign of the mean flow and the horizontal phase velocity of the wave. In particular, at the special time where the mean flow has the same sign everywhere, one of the waves is filtered out and only the other wave is measured at the bottom.

We have also shown an example of a mean flow as a function of time for a fixed forcing. The example we have chosen correspond to a subcritical bifurcation: there is no mean flow during a very long time, and then it appears due to an unknown perturbation. The space-time diagram also shows that the oscillation of the mean flow is very stable, and last more that 9 periods. This is significantly more than the duration of the oscillation of the mean flow reported in previous works [21, 28]. This is due to the use of a pump to continuously inject fluid at the bottom and remove fluid at the top, which keep the stratification linear even close to the membranes and the injection of waves efficient.

The experimental mean flow at a given height can be well fitted by a sinusoidal fit. This is consistent with the weakly nonlinear approach used in the analytical model and is used to obtain the relative phase and the amplitude of the mean flow as a function of height $z$. The amplitude profile of the mean flow is qualitatively similar for the experimental, analytical and numerical results. In both cases, the amplitude is non-monotonous. We recover the main characteristics of the mean flow profile: a phase difference of $\pi$ between the top and the bottom, consistent with the usual qualitative description of the QBO [1]. The variation of the amplitude of the mean flow with height is not monotonic, displaying a local minimum in theory, simulations and experiments.

We have shown a set of values of the period of the mean flow as a function of the forcing. The experiments and the simulations are in good qualitative agreement and show that this period decreases when the amplitude of the forcing increases.

We also discuss in detail one of the main result of [29]: the bifurcation is always a Hopf one: when a mean flow is generated, it is always oscillating at a finite frequency, but the transition can be subcritical (sudden, possible hysteresis) or supercritical (continuous). The experimental results for the bifurcation are in agreement with the analytical results, and we have been able to find two different configurations for which the bifurcation is subcritical (small value of $N$, NaCl) or supercritical (large value of $N$, MgCl$_2$).

One of the perspective for the analog experiments is to study the effect of the ratio of the dissipation length to the height of water, and to force the membrane with two frequencies. The goal is to reproduce experimentally the numerical results of [42] about the robustness of the periodic behavior of the QBO.
Another perspective for the analog experiments is to investigate a less regular forcing. It is surely feasible to use a similar set-up with membrane, and apply an irregular forcing, or at least use only a fraction of the motors. Obtaining a QBO without membrane is very difficult and has not been achieved yet. No oscillation has been observed in numerical simulations of the experiments for random forcing in the case of perfectly smooth membrane [27]. Even with a very careful optimisation of the forcing and of the stratification it is not clear that such an experiment will work.

The geophysics perspective is the use of the tools of nonlinear physics, in particular bifurcation theory, for the QBO in GCM (general circulation model). These tools are becoming (again) very popular in atmospheric science to study the effects of climate change, and the bifurcations are one of the types of the so-called "tipping points" [43-45]. The QBO is expected to decrease in amplitude due to climate change [7]. The experiments have shown that the bifurcation can be subcritical or supercritical, depending mainly on the effective dissipation form of the mean flow. A first question in GCM will be to identify the dissipation of the main flow, which is neither molecular viscosity nor friction on the wall like in the experiment, but probably involves turbulence and maybe the coupling of the velocity field with other fields. An even more relevant question that can be addressed using GCM is the nature of the bifurcation: will the QBO disappear, and if yes continuously or suddenly? A first approach would be to plot the amplitude and period as a function of the CO₂ concentration or the wave forcing and search for the variations of the mean flow close to the onset of instability.

## Appendix A. Fabrication of the seeding beads

To make the seeding beads, we use the following ingredients:

* *Paraffin wax:* We used paraffin wax with a solidification point at 52–54 °C, supplier VWR, reference 26157.291

* Aluminum oxide (alumina) 3 µm powder. Supplier Lam Plan reference A003 P105K

* *Hand washing liquid:* We use Green Lime Hand Washing Liquid 500 ml from the brand "Arbre Vert"

It is very probable not necessary to use exactly the same ingredients to obtain beads, we just mention the ones that we use. We are not even sure that the hand washing liquid is necessary.

The protocol that we use is given below.

We mix 31 g of paraffin, 6 g of alumina and one droplet of washing-up liquid in a borosilicate glass beaker with a teflon stirring bar. Due to the size of alumina particles, we use a FFP2 mask for safety.

The beaker should be placed in a boiling water bath ("bain-marie") and stirred. To do this, we put a crystallizing dish with water on a magnetic hotplate stirrer. The beaker with wax was placed inside the water, close enough to the surface of the hotplate so that the bar can stir the mixture.

When the wax is liquid, we poured the mixture in a plastic weighing boat to cool it down. It is possible to use a spoon to take the alumina which remains at the bottom of the beaker.

Then we let the alumina particles sediment. We started to stir again when the wax is very viscous but not yet solid.

When the wax is cold but still malleable by hand, we made a sphere with the mixture. We waited then one night so that the mixture is hard.

Then we used a chopper to make the particles (Moulinette Mini Chopper from the brand Moulinex).

It is important to stop after a few seconds (5-20) to avoid any melting of the wax.

The last step is to sieve the particles. The amount of beads around 1 mm is about 5–10 g.
## Appendix B. Linear and weakly nonlinear solution

Equations (4) and (5) are investigated analytically close to the threshold in the following.

We will determine the linear threshold, the eigenmodes and where the bifurcation is super-critical or sub-critical in the parameter space.

This will be achieved using a weakly non-linear expansion up to the third order in the amplitude of the mean flow $u$.

At this order, the operator $-\partial D/\partial z$ writes (see equation (5)):

$$ -\frac{\partial D}{\partial z} = F(Lu + \Phi(u)) \qquad (24) $$

where $L$ is the linear operator:

$$ Lu = -\frac{\partial}{\partial z} (\exp(-z) \ell(u)) = \frac{\partial}{\partial z} \left( 2(2\alpha_1 + 4\alpha_2) \exp(-z) \int_0^z u \, dz \right) \qquad (25) $$

$\ell(u)$ is defined by

$$ \ell(u) = -2(2\alpha_1 + 4\alpha_2) \int_0^z u \, dz' \qquad (26) $$

$\Phi$ is the non-linear operator.

$$ \Phi(u) = -\frac{\partial}{\partial z} (\exp(-z) \xi(u)) \qquad (27) $$

and $\xi(u)$ is given by

$$ \xi(u) = -2(4\alpha_1 + 20\alpha_2) \int_0^z u^3 \, dz' + 2(2\alpha_1 + 4\alpha_2)(3\alpha_1 + 10\alpha_2) \left( \int_0^z u \, dz' \right) \left( \int_0^z u^2 \, dz' \right) - \frac{1}{3}(2\alpha_1 + 4\alpha_2)^3 \left( \int_0^z u \, dz' \right)^3 \qquad (28) $$

### B.1. Linear order

At linear order, the model (4) reduces to:

$$ \frac{\partial u}{\partial t} = F \frac{\partial}{\partial z} \left( 2(2\alpha_1 + 4\alpha_2) \exp(-z) \int_0^z u \, dz \right) + \Lambda_1 \frac{\partial^2 u}{\partial z^2} - \Lambda_2 u \qquad (29) $$

We introduce the variable $\psi$ like in [20]:

$$ \psi = \int_0^z u \, dz \qquad (30) $$

Then (29) becomes:

$$ \frac{\partial^2 \psi}{\partial t \partial z} = F \frac{\partial}{\partial z} \left( 2(2\alpha_1 + 4\alpha_2) \exp(-z) \psi \right) + \Lambda_1 \frac{\partial^3 \psi}{\partial z^3} - \Lambda_2 \frac{\partial \psi}{\partial z} \qquad (31) $$

Using the linearity of the equations, we use complex variables, so that the solution is the real part of the corresponding complex variable. To simplify the notations, we use the same variable for the complex and real variables. We are considering eigenmodes of (31) so that:

$$ \frac{\partial \psi}{\partial t} = \mu \psi \qquad (32) $$

$\mu \in \mathbb{C}$. At the threshold, the real part of the growth rate vanishes and thus $\mu \in i\mathbb{R}$.

We then integrate (31):

$$ \mu \psi = F \times 2(2\alpha_1 + 4\alpha_2) \exp(-z) \psi + \Lambda_1 \frac{\partial^2 \psi}{\partial z^2} - \Lambda_2 \psi + K_1 \qquad (33) $$

where $K_1$ is an unknown complex number. This equation can be rewritten:

$$ \frac{\partial^2 \psi}{\partial z^2} + (a \exp(-z) - b) \psi = K, \qquad (34) $$
with:

$$a = \frac{F \times 2(2\alpha_1 + 4\alpha_2)}{\Lambda_1}, \qquad (35)$$

$$b = \frac{\mu + \Lambda_2}{\Lambda_1}, \qquad (36)$$

and K is unknown.

We first solve the homogeneous equation:

$$\frac{\partial^2 \psi}{\partial z^2} + (a \exp(-z) - b) \psi = 0. \qquad (37)$$

The solutions are:

$$\psi = J\left(\pm 2\sqrt{b}, 2\sqrt{a}\exp(-z/2)\right), \qquad (38)$$

where J is the Bessel function of the first kind, of order $\pm 2\sqrt{b}$ and of argument $2\sqrt{a}\exp(-z/2)$.

The function $J(-2\sqrt{b}, 2\sqrt{a}\exp(-z/2))$ does not have a finite limit when $z \to +\infty$, and we disregard this solution due to the boundary condition (see equation (7)).

The solution of (37) is thus:

$$\psi_0 = J\left(2\sqrt{b}, 2\sqrt{a}\exp(-z/2)\right), \qquad (39)$$

up to a multiplicative constant, which is irrelevant for the linear problem.

We then solve (34), and we search a particular solution using the method of the variation of parameter:

$$\psi(z) = f_1(z)\psi_0(z). \qquad (40)$$

Using equation (34) and that $\psi_0$ is solution of the homogeneous equation, we obtain:

$$2 \frac{df_1}{dz} \frac{d\psi_0}{dz} + f_1 \frac{d^2\psi_0}{dz^2} = K. \qquad (41)$$

We multiply the previous equation by $\psi_0$:

$$\frac{d}{dz} \left( \psi_0^2 \frac{df_1}{dz} \right) = K\psi_0, \qquad (42)$$

which gives:

$$\psi_0^2 \frac{df_1}{dz} = \int_{+\infty}^{z} K \psi_0(z') dz', \qquad (43)$$

assuming that $f_1$, $\frac{df_1}{dz}$ vanish when $z \to +\infty$ and using that $\psi_0$ vanishes when $z \to +\infty$.

We write $f_1 = Kf$, so that we have:

$$f(z) = \int_{z_0}^{z} \frac{dz' \int_{+\infty}^{z'} \psi_0(z'') dz''}{(\psi_0(z'))^2}, \qquad (44)$$

where $z_0$ is an arbitrary positive constant; the final solution does not depend on the value chosen for $z_0$.

So the solution of (34) is thus:

$$\psi(z) = \psi_0(z)(1 + Kf(z)), \qquad (45)$$

where $\psi_0$ is given by (39) and $f$ by (44).

We can also write the solution for $u$:

$$u = \frac{\partial}{\partial z} \left( \psi_0(z) (1 + K f(z)) \right), \qquad (46)$$
Using (7), (30), we have the following boundary conditions:

$$ \psi(0) = 0, \tag{47} $$

$$ \frac{d\psi}{dz}(0) = 0, \tag{48} $$

$$ \frac{d\psi}{dz}(+\infty) = 0, \tag{49} $$

Using (47) in (45) gives:

$$ 0 = \psi_0(0)(1 + Kf(0)). \tag{50} $$

So that:

$$ K = \frac{-1}{f(0)} = \frac{-1}{\int_{z_0}^{0} \frac{\int_{+\infty}^{z} \psi_0(z') dz'}{(\psi_0(z))^2} dz}. \tag{51} $$

Using (48) in (45) gives:

$$ 0 = \psi_0'(0)(1 + Kf(0)) + \psi_0(0)(Kf'(0)) = \psi_0(0)(Kf'(0)). \tag{52} $$

We thus have $f'(0) = 0$, i.e.:

$$ \int_0^{+\infty} \psi_0(z) \, dz = 0 \tag{53} $$

i.e.:

$$ \int_0^{+\infty} J(2\sqrt{b}, 2\sqrt{a}\exp(-z/2)) \, dz = 0 \tag{54} $$

This equation allows us to compute the threshold of the instability.

Using the new variable

$$ v = \exp(-z/2), \tag{55} $$

equation (54) gives:

$$ \int_0^1 \frac{J(2\sqrt{b}, 2\sqrt{a}v)}{v} \, dv = 0. \tag{56} $$

Using the following property of Bessel functions:

$$ \frac{2\alpha}{y} J(\alpha, y) = J(\alpha - 1, y) + J(\alpha + 1, y), \tag{57} $$

we finally obtain equation (17).

## B.2. Non-linear order

Using equation (24), equation (4) is considered up to the third order in the amplitude of $u$:

$$ \frac{\partial u}{\partial t} = F(Lu + \Phi(u)) + \Lambda_1 \frac{\partial^2 u}{\partial z^2} - \Lambda_2 u \tag{58} $$

More precisely, we note

$$ F = F_0 + \epsilon F_1 \tag{59} $$

where $F_0$ is the value at the instability threshold ($\Lambda_1$ and $\Lambda_2$ being fixed) and $\epsilon$ a small parameter (dimensionless distance to the threshold).

We introduce a slow timescale $T = \epsilon t$ which correspond to the slow variation of the amplitude of the oscillation. Using the multiple scale expansion, $\frac{\partial u}{\partial t}$ should be replaced by:

$$ \frac{\partial u}{\partial t} + \frac{\partial u}{\partial T} \frac{\partial T}{\partial t} = \frac{\partial u}{\partial t} + \epsilon \frac{\partial u}{\partial T} \tag{60} $$

We choose the expansion for $u$:

$$ u = \sqrt{\epsilon} (u_0 + \epsilon u_1 + \epsilon^2 u_2 + \dots) \tag{61} $$
At order $\epsilon^{1/2}$, equation (58) gives:

$$ \frac{\partial u_0}{\partial t} = F_0 L u_0 + \Lambda_1 \frac{\partial^2 u_0}{\partial z^2} - \Lambda_2 u_0 \quad (62) $$

This is the linear order (see equation (29)), which has already been solved. The solution is proportional to

$$ u = \frac{\partial \psi}{\partial z} \quad (63) $$

where $\psi$ is given by equations (45), (51), (44), and (39).

We write $u_0$ as a product of the eigenmode of the linear problem $u$ with a slowly varying amplitude:

$$ u_0(t, T, z) = \frac{1}{2} (A(T)u(t, z) + A^*(T)u^*(t, z)), \quad (64) $$

where $u$ is given by equation (63).

At order $\epsilon^{3/2}$, equation (58) gives:

$$ \frac{\partial u_1}{\partial t} + \frac{\partial u_0}{\partial T} = F_0 L u_1 + F_1 L u_0 + \Phi(u_0) + \Lambda_1 \frac{\partial^2 u_1}{\partial z^2} - \Lambda_2 u_1 \quad (65) $$

which can be rewritten:

$$ \underbrace{\frac{\partial u_1}{\partial t} - F_0 L u_1 - \Lambda_1 \frac{\partial^2 u_1}{\partial z^2} + \Lambda_2 u_1}_{\mathcal{L} u_1} = \underbrace{- \frac{\partial u_0}{\partial T} + F_1 L u_0 + \Phi(u_0)}_{\mathcal{H} u_0} \quad (66) $$

The next step is to calculate the adjoint operator of $\mathcal{L}$. We introduce the (hermitian) scalar product:

$$ \langle a | b \rangle = \lim_{T \to +\infty} \frac{2}{T} \int_{-T}^{T} \int_{0}^{+\infty} a^* b \, dt \, dz' \quad (67) $$

where $a^*$ is the complex conjugate of $a$.

The fields verify the following properties: they do not diverge when $t \to \pm\infty$ and they vanish when $x=0$ or $+\infty$.

After integration by parts we obtain for the adjoint operator:

$$ \mathcal{L}^\dagger c = + \frac{\partial c}{\partial t} + F_0 \times 2(2\alpha_1 + 4\alpha_2) \int_0^z \exp(-z') \frac{\partial c}{\partial z} dz' + \Lambda_1 \frac{\partial^2 c}{\partial z^2} - \Lambda_2 c \quad (68) $$

We now search a vector $q$ in the kernel of $\mathcal{L}^\dagger$:

$$ \mathcal{L}^\dagger q = 0 \quad (69) $$

We thus have:

$$ \frac{\partial}{\partial z} (\mathcal{L}^\dagger q) = 0 \quad (70) $$

We define $\varphi$ by:

$$ \varphi = \frac{\partial q}{\partial z} \quad (71) $$

We search modes proportional to $\exp(\kappa t)$ so that

$$ \frac{\partial \varphi}{\partial t} = \kappa \varphi \quad (72) $$

We thus obtain:

$$ \frac{\partial^2 \varphi}{\partial z^2} + (a \exp(-z) - B) \varphi = 0 \quad (73) $$

with:

$$ a = \frac{F_0 \times 2(2\alpha_1 + 4\alpha_2)}{\Lambda_1} \quad (74) $$

$$ B = \frac{-\kappa + \Lambda_2}{\Lambda_1} \quad (74) $$
This equation is nearly identical to the one solved for the direct problem (equation (34)). We thus
have:

$$
\varphi = J(2\sqrt{B}, 2\sqrt{a}\exp(-z/2)) \tag{75}
$$

so that we have finally obtained an expression of an element of the kernel of the adjoint of the
linear operator:

$$
q = \int_{0}^{z} J(2\sqrt{B}, 2\sqrt{a}\exp(-z'/2)) dz' \quad (76)
$$

From equation (66), the solvability condition then writes:

$$
\langle q | \mathcal{L} u_1 \rangle = \left\langle q \left| -\frac{\partial u_0}{\partial T} + F_1 L u_0 + \Phi(u_0) \right\rangle \right| = 0 \quad (77)
$$

Using the expression of $u_0$, see eq. (64), we obtain the amplitude equation:

$$
\frac{\partial A}{\partial T} = A F_1 \frac{\langle q | Lu \rangle}{\langle q | u \rangle} + |A|^2 A \frac{\langle q | \Phi(\frac{1}{2}(u(t,z) + u^*(t,z))) \rangle}{\langle q | u \rangle}, \quad (78)
$$

from which we obtain the expression of the coefficients $\alpha$ and $\beta$ as defined in Eq. (21):

$$
\alpha = \frac{\langle q | Lu \rangle}{\langle q | u \rangle}, \tag{79}
$$

$$
\beta = \frac{\langle q | \Phi(\frac{1}{2}(u(t,z) + u^*(t,z))) \rangle}{\langle q | u \rangle}. \quad (80)
$$

The different coefficients are given by

$$
\begin{align}
\langle q | u \rangle &= \int_0^{+\infty} \left( \int_0^z J(2\sqrt{B}, 2\sqrt{a}\exp(-z'/2)) dz' \right)^* \left( \frac{1}{2} \frac{\partial \psi}{\partial z} \right) dz, \tag{81} \\
&= -\frac{1}{2} \int_0^{+\infty} J(2\sqrt{B}, 2\sqrt{a}\exp(-z'/2))^* \psi dz. \nonumber
\end{align}
$$

where $\psi$ is given by eq (45), and

$$
\begin{align}
\frac{1}{2(2\alpha_1 + 4\alpha_2)} \langle q | Lu \rangle &= \int_0^{+\infty} \left( \int_0^z J(2\sqrt{B}, 2\sqrt{a}\exp(-z'/2)) dz' \right)^* \left( \frac{1}{2} \frac{\partial}{\partial z} (\exp(-z)\psi(z)) \right) dz \nonumber \\
&= -\frac{1}{2} \int_0^{+\infty} J(2\sqrt{B}, 2\sqrt{a}\exp(-z'/2))^* \exp(-z) \psi dz \tag{82}
\end{align}
$$

and for the non-linear term by

$$
\left\langle q \left| \Phi \left( \frac{1}{2} (u(t,z) + u^*(t,z)) \right) \right\rangle \right\rangle = \left\langle q \left| \frac{\partial}{\partial t} \left( \exp(-z) \xi(u) \right) \right\rangle \right\rangle = Q_1 + Q_2 + Q_3, \quad (83)
$$

with:

$$
Q_1 = 2(4\alpha_1 + 20\alpha_2) \int_0^{+\infty} dz J(2\sqrt{B}, 2\sqrt{a}\exp(-z/2))^* \exp(-z) \int_0^z \frac{3}{8} |u(z')|^2 u(z') dz',
$$

$$
Q_2 = -2(2\alpha_1 + 4\alpha_2)(3\alpha_1 + 10\alpha_2) \int_0^{+\infty} dz J(2\sqrt{B}, 2\sqrt{a}\exp(-z/2))^* \exp(-z) \\
\frac{1}{8} \left( \left( \int_0^z u(z')^* dz' \right) \left( \int_0^z u(z')^2 dz' \right) + 2 \left( \int_0^z u(z') dz' \right) \left( \int_0^z |u(z')|^2 dz' \right) \right),
$$

$$
Q_3 = \frac{1}{3} (2\alpha_1 + 4\alpha_2)^3 \int_0^{+\infty} dz J(2\sqrt{B}, 2\sqrt{a}\exp(-z/2))^* \exp(-z) \frac{3}{8} |\psi(z)|^2 \psi(z).
$$
Acknowledgements

We gratefully acknowledge J. da Silva Quintas, C. Goncalves, E. Nicolau, C. Herrmann, and L. Bonnet for their technical support. We thank our collaborators S. Fauve, N. Garroum and G. Facchini for discussions and with whom [29] or [35] were written.

Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

References

[1] M. P. Baldwin, L. J. Gray, T. J. Dunkerton, et al., "The quasi-biennial oscillation", *Rev. Geophys.* **39** (2001), pp. 179–229.

[2] T. Fouchet, S. Guerlet, D. F. Strobel, A. A. Simon-Miller, B. Bézard and F. M. Flasar, "An equatorial oscillation in Saturn's middle atmosphere", *Nature* **453** (2008), pp. 200–202.

[3] X. Li and P. L. Read, "A mechanistic model of the quasi-quadrennial oscillation in Jupiter's stratosphere", *Planet. Space Sci.* **48** (2000), no. 7-8, pp. 637–669.

[4] A. P. Showman, X. Tan and X. Zhang, "Atmospheric Circulation of Brown Dwarfs and Jupiter- and Saturn-like Planets: Zonal Jets, Long-term Variability, and QBO-type Oscillations", *Astrophys. J.* **883** (2019), no. 1, article no. 4.

[5] E.-J. Kim and K. B. MacGregor, "Gravity Wave-driven Flows in the Solar Tachocline", *Astrophys. J.* **556** (2001), no. 2, pp. L117–L120.

[6] C. Pinçon, K. Belkacem and M.-J. Goupil, "Implications of the generation of internal gravity waves by penetrative convection for the internal rotation evolution of low-mass stars", *EPJ Web Conf.* **160** (2017), article no. 02002.

[7] J. A. Anstey, S. M. Osprey, J. Alexander, et al., "Impacts, processes and projections of the quasi-biennial oscillation", *Nat. Rev. Earth Environ.* **3** (2022), no. 9, pp. 588–603.

[8] Y. Wang, J. Rao, Y. Lu, Z. Ju, J. Yang and J. Luo, "A revisit and comparison of the quasi-biennial oscillation (QBO) disruption events in 2015/16 and 2019/20", *Atmos. Res.* **294** (2023), article no. 106970.

[9] N. Butchart, "The stratosphere: a review of the dynamics and variability", *Weather Clim. Dynam.* **3** (2022), no. 4, pp. 1237–1272.

[10] C. C. Collimore, D. W. Martin, M. H. Hitchman, A. Huesmann and D. E. Waliser, "On The Relationship between the QBO and Tropical Deep Convection", *J. Climate* **16** (2003), no. 15, pp. 2552–2568.

[11] K. P. Bowman, "Global Patterns of the Quasi-biennial Oscillation in Total Ozone", *J. Atmos. Sci.* **46** (1989), no. 21, pp. 3328–3343.

[12] A. G. Marshall and A. A. Scaife, "Impact of the QBO on surface winter climate", *J. Geophys. Res.* **114** (2009), article no. D18110.

[13] Z. Martin, S.-W. Son, A. Butler, H. Hendon, H. Kim, A. Sobel, S. Yoden and C. Zhang, "The influence of the quasi-biennial oscillation on the Madden-Julian oscillation", *Nat. Rev. Earth Environ.* **2** (2021), no. 7, pp. 477–489.

[14] L. L. Hood and C. A. Hoopes, "Arctic Sea Ice Loss, Long-Term Trends in Extratropical Wave Forcing, and the Observed Strengthening of the QBO-MJO Connection", *J. Geophys. Res. Atmos.* **128** (2023), no. 24, article no. e2023JD039501.

[15] M. Takahashi, "Simulation of the Quasi-Biennial Oscillation in a general circulation model", *Geophys. Res. Lett.* **26** (1999), no. 9, pp. 1307–1310.

[16] N. Butchart, J. A. Anstey, K. Hamilton, et al., "Overview of experiment design and comparison of models participating in phase 1 of the SPARC Quasi-Biennial Oscillation initiative (QBOi)", *Geosci. Model Dev.* **11** (2018), no. 3, pp. 1009–1032.

[17] F. Lott and L. Guez, "A stochastic parameterization of the gravity waves due to convection and its impact on the equatorial stratosphere", *J. Geophys. Res. Atmos.* **118** (2013), pp. 8897–8909.

[18] U. Achatz, M. J. Alexander, E. Becker, et al., "Atmospheric Gravity Waves: Processes and Parameterization", *J. Atmos. Sci.* **81** (2024), no. 2, pp. 237–262.

[19] R. S. Lindzen and J. R. Holton, "A theory of the Quasi-biennial Oscillation", *J. Atmos. Sci.* **25** (1968), pp. 1095–1107.

[20] R. A. Plumb, "The interaction of two internal waves with the mean flow: implications for the theory of the quasi-biennial oscillations", *J. Atmos. Sci.* **34** (1977), pp. 1847–1858.

[21] R. A. Plumb and A. D. McEwan, "The instability of a Forced Standing Wave in a Viscous Stratified Fluid : A Laboratory Analogue of the Quasi-Biennial Oscillation", *J. Atmos. Sci.* **35** (1978), pp. 1827–1839.
[22] A. Renaud, L.-P. Nadeau and A. Venaille, "Periodicity Disruption of a Model Quasi-biennial Oscillation of Equatorial Winds", Phys. Rev. Lett. 122 (2019), no. 21, article no. 214504.

[23] A. Renaud and A. Venaille, "On the Holton-Lindzen-Plumb model for mean flow reversals in stratified fluids", Q. J. R. Meteorol. Soc. 146 (2020), no. 732, pp. 2981-2997.

[24] P. Léard, D. Lecoanet and M. Le Bars, "Multimodal Excitation to Model the Quasi-biennial Oscillation", Phys. Rev. Lett. 125 (2020), no. 23, article no. 234501.

[25] L.-A. Couston, D. Lecoanet, B. Favier and M. Le Bars, "Order Out of Chaos: Slowly Reversing Mean Flows Emerge from Turbulently Generated Internal Waves", Phys. Rev. Lett. 120 (2018), no. 24, article no. 244505.

[26] H.-H. Bui, E. Nishimoto and S. Yoden, "Downward Influence of QBO-Like Oscillation on Moist Convection in a Two-Dimensional Minimal Model Framework", J. Atmos. Sci. 74 (2017), no. 11, pp. 3635-3655.

[27] N. P. Wedi and P. K. Smolarkiewicz, "Direct numerical simulations of the Plumb-McEwan laboratory analog of the QBO", J. Atmos. Sci. 63 (2006), pp. 3226-3252.

[28] N. Otobe, S. Sakai, S. Yoden and M. Shiotani, "Visualization and WKB analysis of the internal gravity wave in the QBO experiment", Nagare multimedia 17 (1998). preprint, https://www2.nagare.or.jp/mm/98/otobe/index.htm.

[29] B. Semin, N. Garroum, F. Pétrélis and S. Fauve, "Nonlinear saturation of the large scale flow in a laboratory model of the quasibiennial oscillation", Phys. Rev. Lett. 121 (2018), no. 13, article no. 134502.

[30] P. Léard, B. Favier, P. Le Gal and M. Le Bars, "Coupled convection and internal gravity waves excited in water around its density maximum at 4C", Phys. Rev. Fluids 5 (2020), article no. 024801.

[31] M. Le Bars, L.-A. Couston, B. Favier, P. Léard, D. Lecoanet and P. Le Gal, "Fluid dynamics of a mixed convective/stably stratified system—A review of some recent works", Comptes Rendus. Physique 21 (2020), no. 2, pp. 151-164.

[32] P. Léard, L. Margaillan, T. Raymond, M. Rouby and M. Le Bars, "Large-scale flow driven by turbulently generated internal gravity waves", Phys. Rev. Fluids 6 (2021), no. 8, article no. 084801.

[33] B. King, H. P. Zhang and H. L. Swinney, "Tidal flow over three-dimensional topography in a stratified fluid", Phys. Fluids 21 (2009), article no. 116601.

[34] G. Bordes, A. Venaille, S. Joubaud, P. Odier and T. Dauxois, "Experimental observation of a strong mean flow induced by internal gravity waves", Phys. Fluids 24 (2012), article no. 086602.

[35] B. Semin, G. Facchini, F. Pétrélis and S. Fauve, "Generation of a mean flow by an internal wave", Phys. Fluids 28 (2016), article no. 096601.

[36] H. Ben-Gida, R. Gurka and A. Liberzon, "OpenPIV-MATLAB – An open-source software for particle image velocimetry; test case: Birds' aerodynamics", SoftwareX 12 (2020), article no. 100585.

[37] H. Cai, Bessel Function of Complex Order and Argument. MATLAB Central File Exchange, accessed: 2024-02-26, https://math.stackexchange.com/questions/585240/bessel-function-with-complex-argument.

[38] S. Yoden and J. Holton, "A new look at equatorial quasi-biennial oscillation Models", J. Atmos. Sci. 45 (1988), pp. 2703-2717.

[39] T. Dauxois, S. Joubaud, P. Odier and A. Venaille, "Instabilities of Internal Gravity Wave Beams", Ann. Rev. Fluid Mech. 50 (2018), no. 1, pp. 131-156.

[40] R. Godoy-Diana, J.-M. Chomaz and C. Donnadieu, "Internal gravity waves in a dipolar wind: a wave-vortex interaction experiment in a stratified fluid", J. Fluid Mech. 548 (2006), pp. 281-308.

[41] K. Hamilton, "Dynamics of the Tropical Middle Atmosphere: A Tutorial Review", Atmosphere-Ocean 36 (1998), no. 4, pp. 319-354.

[42] X. Chartrand, L.-P. Nadeau and A. Venaille, "Recovering Quasi-Biennial Oscillations from Chaos", J. Atmos. Sci. 81 (2024), pp. 1213-1224.

[43] S. Bathiany, M. Scheffer, E. H. Van Nes, M. S. Williamson and T. M. Lenton, "Abrupt climate change in an oscillating world", Sci. Rep. 8 (2018), no. 1, article no. 5040.

[44] N. Wunderling, A. von der Heydt, Y. Aksenov, et al., "Climate tipping point interactions and cascades: A review", Earth Syst. Dyn. 15 (2024), pp. 41-74.

[45] Intergovernmental Panel on Climate Change (IPCC), Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, 2023.