1. Introduction
The equatorial stratosphere of Earth is home to zonal winds that change direction roughly every 14 months, a phenomenon known as the quasi-biennial oscillation (QBO). The QBO is fuelled by planetary waves such as Kelvin and Rossby waves, as well as gravity waves that are emitted at the bottom of the stratosphere with a broad range of frequencies and wavelengths, and with a high degree of temporal and spatial intermittency (Baldwin et al. Reference Baldwin2001). Similar oscillating winds are observed on other planets, including Jupiter and Saturn (Dowling Reference Dowling2008; Read Reference Read2018), and may also exist within stars (McIntyre Reference McIntyre1994; Rogers, MacGregor & Glatzmaier Reference Rogers, MacGregor and Glatzmaier2008; Showman, Tan & Zhang Reference Showman, Tan and Zhang2019). The QBO serves as a reference for slow, large-scale mean flows that emerge in stratified and/or rotating fluids under fast, small-scale wave forcing (Vallis Reference Vallis2017). The QBO reversals were long thought robust and regular, compared with the chaotic nature of atmospheric waves, until they showed abrupt phase changes in 2015/2016 (Osprey et al. Reference Osprey, Butchart, Knight, Scaife, Hamilton, Anstey, Schenzinger and Zhang2016; Newman et al. Reference Newman, Coy, Pawson and Lait2016) and 2019/2020 (Anstey et al. Reference Anstey, Banyard, Butchart, Coy, Newman, Osprey and Wright2021), highlighting the need to quantify the variability of QBO-type flows. Today, two mutually non-exclusive sources of variability are proposed: on the one hand, the variability may result from teleconnections with external climatic variability modes (Osprey et al. Reference Osprey, Butchart, Knight, Scaife, Hamilton, Anstey, Schenzinger and Zhang2016; Newman et al. Reference Newman, Coy, Pawson and Lait2016; Fletcher et al. Reference Fletcher, Guerlet, Orton, Cosentino, Fouchet, Irwin, Li, Flasar, Gorius and Morales-Juberías2017); and on the other hand, it may be a form of intrinsic variability, resulting from the stochasticity of the waves generated by the turbulent convection in the troposphere (Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2018; Showman et al. Reference Showman, Tan and Zhang2019).
The first studies that identified waves as the driving force of the QBO were the seminal works of Holton & Lindzen (Lindzen & Holton Reference Lindzen and Holton1968; Holton & Lindzen Reference Holton and Lindzen1972), which described the mechanism of momentum transfer from planetary and gravity waves to horizontal winds. Shortly after, Plumb (Reference Plumb1977) encapsulated the minimal ingredients required for the emergence of a reversing mean flow in an idealised model, subsequently named the Holton–Lindzen–Plumb (HLP) model. Several laboratory experiments have observed wave-driven oscillating mean flows, providing evidence in support of the HLP model (Plumb et al. Reference Plumb and McEwan1978; Semin et al. Reference Semin, Garroum, Pétrélis and Fauve2018). Recently, Renaud, Nadeau & Venaille (Reference Renaud, Nadeau and Venaille2019) and Renaud & Venaille (Reference Renaud and Venaille2020) have shown that oscillating mean flows reproduced by the HLP model under monochromatic wave forcing contain a level of intrinsic variability that increases significantly as the amplitude of the wave forcing increases. In their simulations, the mean flow transitions from a periodic reversal to a bi-stable state and eventually to chaos at higher forcing levels. Subsequently, Léard et al. (Reference Léard, Lecoanet and Le Bars2020) investigated predictions with the HLP model forced by a broadband wave spectrum and showed how the spectral bandwidth influences the transitions towards bi-stability and chaos.
Idealised simulations of a model troposphere–stratosphere configuration recently showed that the chaotic nature of a broadband gravity wave spectrum forcing a QBO-type flow could significantly impact the stability of mean-flow reversals (Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2018). Reproducing a convincingly chaotic atmospheric wave field has been the subject of several stochastic wave parameterisations for global circulation models (Piani, Norton & Stainforth Reference Piani, Norton and Stainforth2004; Eckermann Reference Eckermann2011; Lott, Guez & Maury Reference Lott, Guez and Maury2012; Lott & Guez Reference Lott and Guez2013), although they did not quantify the impact on reversal variability. Very recently, a stochastic forcing was similarly implemented in the context of the HLP model by Ewetola & Esler (Reference Ewetola and Esler2024). In their work, the authors implemented a stochastic wave forcing in the form of a fast monochromatic wave modulated by a slow, random envelope evolving according to a parameterised class of stochastic differential equations (SDE). Importantly, Ewetola & Esler (Reference Ewetola and Esler2024) found that the details of the SDE do not matter; the impact of wave stochasticity on mean-flow statistics depends on a single parameter, which is proportional to the envelope time scale
$\tau$
, as long as
$\tau$
is short compared with the QBO time scale. The study focused primarily on average properties of the QBO, however, leaving unaddressed the question of the variability of the reversals.
In this paper, we perform wave-resolved simulations of the two-dimensional Boussinesq equations, which constitute the starting point for the HLP model, and investigate mean-flow statistics and reversal variability when the wave forcing is monochromatic, polychromatic or stochastic, using a simple Gaussian random field as wave generator. In § 2 we present the model equations and the forcing schemes used in our simulations. In § 3 we present our results and highlight, under Gaussian random forcing, that the amplitude and correlation time of the waves both increase the relative variability of mean-flow reversals. Finally in § 4 we summarise our work, discuss our model limitations and provide guidelines for future work.
2. Methods
2.1. Model equations and simulation details
We investigate the dynamics of a two-dimensional stably stratified fluid with internal gravity waves emitted from the bottom boundary. To this end, we solve the Navier–Stokes equations under the Boussinesq approximation in dimensionless form. Considering a Cartesian coordinate system
$(x,z)$
with
$x$
the horizontal axis and
$z$
the upward vertical axis, the evolution equations for the dimensionless velocity vector
$\boldsymbol{u}$
, buoyancy
$b$
and gravity potential
$\phi$
can be written as
\begin{equation} \begin{cases} \partial _t{\boldsymbol{u}} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{u}} = -\boldsymbol{\nabla }\phi + b\hat {z} + \frac 1{R_w}{\nabla} ^2{\boldsymbol{u}},\\ \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{u}} = 0,\\ \partial _tb + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }) b + \mathcal{N}^2w = \frac 1{R_w \textit{Pr}}{\nabla} ^2b - \varGamma b, \end{cases} \end{equation}
where
$R_w$
is a wave-based Reynolds number,
$\mathcal{N}$
is the dimensionless Brunt–Väisälä frequency,
$\textit{Pr}$
is the Prandtl number and
$\varGamma$
is the dimensionless radiative damping rate. The characteristic time and length scales used to derive (2.1) are the inverse forcing wave frequency
$1/\omega _0$
and the inverse forcing wavenumber
$1/k_0$
, respectively. We use scales that are typical of the waves and not of the mean flow, as wave scales will remain constant in our simulations while mean-flow scales will vary with our parametrisation of the forcing. Thus, the four dimensionless parameters in (2.1) are expressed as
\begin{equation} R_w=\frac {\omega _0}{\nu {k_0}^2},\quad \mathcal{N}=\frac {N}{\omega _0},\quad \textit{Pr}=\frac \nu \kappa ,\quad \varGamma =\frac \gamma {\omega _0}, \end{equation}
where
$\nu$
is the kinematic viscosity,
$N$
is the Brunt–Väisälä frequency,
$\kappa$
is the thermal diffusivity and
$\gamma$
is the radiative damping rate. The four parameters of (2.2) are constant in all our simulations, with prescribed values
$R_w=765$
,
$\mathcal{N}=3.33$
,
$\textit{Pr}=1$
and
$\varGamma ={6.9\times 10^{-2}}$
. These values are typical of the stratosphere forced by short internal gravity waves (Vallis Reference Vallis2017), except for the wave-based Reynolds number
$R_w$
, which has typically values closer to
$50\ 000$
for such waves (Renaud et al. Reference Renaud, Nadeau and Venaille2019). Lower values of
$R_w$
are necessary in our simulations to keep the numerical cost reasonable. By keeping these parameters constant across all simulations, we focus our analysis on the sensitivity of the mean-flow dynamics to the parameters of the bottom boundary condition, which control the temporal content of wave emission, as described in § 2.2.
The domain size in the horizontal direction is
$L=4\pi$
, such that we force two horizontal wavelengths at the bottom boundary. For the domain height we take
$H=20$
. Waves are damped as they propagate upward because of radiative cooling, viscosity and thermal diffusivity. In all simulations, the total damping is dominated by radiative damping, as the radiative damping length
$\varLambda _{\textit{rad}}=\sqrt {\mathcal{N}^2-1}/(\varGamma \mathcal{N}^2)=4.15$
is much smaller than both the viscous and molecular damping lengths
$\varLambda _{{v\textit{isc}}}=R_w/\mathcal{N}^3=20.6$
and
$\varLambda _{\textit{mol}}= \textit{Pr}\varLambda _{{v\textit{isc}}}$
, respectively. As a result, the domain height is roughly 5 times the damping length
$\varLambda \approx \varLambda _{\textit{rad}}$
. We choose periodic boundary conditions in the
$x$
direction. In the
$z$
direction, we use a no-slip, zero-buoyancy condition on top, and we set the pressure gauge by requiring that the horizontal mean of
$\phi$
is zero on the top boundary. We also impose a sponge layer on both velocity and buoyancy on the last tenth of the domain height to ensure that no waves are reflected from the top boundary downward. On the bottom boundary, we impose
$(x,t)$
-dependent wave shapes for the velocity and buoyancy, as described in § 2.2. We solve equations numerically with the pseudospectral open-source code Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). We expand variables in the
$x$
and
$z$
directions using 64 Fourier modes and 256 Chebyshev modes, respectively. The numerical solution is advanced in time using a second-order multistage explicit–implicit Runge–Kutta scheme with the Courant–Friedrichs–Lewy parameter set to 0.5.
In order to speed up the emergence of the mean flow and reduce computing times, we seed all simulations with a vertically sheared jet at the initial time. The initial jet amplitude is equal to one tenth of the amplitude of the expected mean flow, which should be close to the phase velocity
$c=1$
of the forced wave. The initial transient growth then typically lasts between
$5000$
and
$10\ 000$
time units, and we run all simulations until roughly
$t=160\ 000$
. A snapshot of the velocity field from one simulation is shown in figure 1.
Figure 1. Snapshot showing the instantaneous horizontal velocity field
$u(x,z)$
at
$t=30\ 000$
(see colour bar to the right) and the vertical profile of the horizontally averaged velocity
$\overline {u}(z)$
(black line).
2.2. Wave emission from the bottom boundary
All simulations are forced by imposing time-varying horizontal profiles for the velocities
$u$
and
$w$
and the buoyancy
$b$
on the bottom boundary. The forcing profiles are separated in space and time and written in dimensionless form as
\begin{equation} \begin{cases} u(x,z=0,t) = \sqrt {\mathcal{F}}\cos (x)u_0(t),\\ w(x,z=0,t) = \sqrt {\mathcal{F}}\sin (x)w_0(t),\\ b(x,z=0,t) = \sqrt {\mathcal{F}}\sin (x)b_0(t), \end{cases} \end{equation}
where
$\mathcal{F}$
is the total directional horizontal wave momentum flux, i.e. the absolute momentum flux sent both to the left and to the right. The mathematical details of the forcing signals
$u_0$
,
$w_0$
and
$b_0$
are given in Appendix A. From the value
$\mathcal{F}$
for the input momentum flux we build a forcing Reynolds number defined as
\begin{equation} \textit{Re} = \frac {\mathcal{F}\varLambda }{c\nu }, \end{equation}
where
$\varLambda =({\varLambda _{\textit{rad}}}^{-1}+{\varLambda _{{v\textit{isc}}}}^{-1}+{\varLambda _{\textit{mol}}}^{-1})^{-1}$
is the total damping length and
$c=\omega _0/k_0$
is the forcing wave phase velocity. The forcing Reynolds number
$\textit{Re}$
is equal to the effective Reynolds number in Renaud et al. (Reference Renaud, Nadeau and Venaille2019) and is varied between simulations to control the forcing intensity. The precise choice of the cosine/sine horizontal profiles in (2.3) derives from the polarisation relations for linear gravity waves and imposes some constraints on the relations between the Fourier coefficients of
$u_0$
,
$w_0$
and
$b_0$
. By enforcing these relations we make sure to force linear wave modes, with two horizontal wavelengths in the domain width.
We test the influence of the wave emission process on the reversing mean-flow dynamics by considering three different types of forcing wave spectra. We first consider a monochromatic forcing, as it has been used in numerous studies since the derivation of the HLP model (Holton & Lindzen Reference Holton and Lindzen1972; Plumb Reference Plumb1977; Renaud & Venaille Reference Renaud and Venaille2020). Our second forcing scheme builds upon a multimodal, continuous wave spectrum with a finite frequency bandwidth, similar to the polychromatic scheme investigated in Léard et al. (Reference Léard, Lecoanet and Le Bars2020). Finally, our third forcing scheme is stochastic in time and is built using random fields techniques that were recently developed in turbulence modelling (Pereira, Garban & Chevillard Reference Pereira, Garban and Chevillard2016; Chevillard et al. Reference Chevillard, Garban, Rhodes and Vargas2019; Reneuve & Chevillard Reference Reneuve and Chevillard2020).
The main characteristics of the three forcing schemes are illustrated in figure 2. The monochromatic forcing is characterised by a pure sine of constant amplitude in the time domain (top row), which corresponds to a single frequency
$\omega =1$
in Fourier space (bottom row). We recall the ratio of
$\omega _0$
to
$N$
is the same in all simulations with
$\mathcal{N}= {N}/{\omega _0}=3.33$
. The polychromatic forcing is designed to generate a fast-oscillating carrier wave modulated by a slower, periodic envelope of period
$\tau$
. The modulated wave signal, which repeats itself every
$\tau$
time units, is defined in frequency domain as a narrow-band spectrum of bandwidth
$\Delta \omega \sim 1/\tau$
. We set the carrier wave frequency of the polychromatic signal equal to the monochromatic forcing wave frequency and use the modulation period
$\tau$
as free parameter. The stochastic forcing is designed to generate a fast-oscillating carrier wave, with the same carrier wave properties as for the other two forcing schemes. The slowly varying envelope is correlated over the same characteristic time scale
$\tau$
as for the polychromatic forcing. However, the modulation is now random. This choice allows the stochastic forcing signal to be very similar to the polychromatic (but deterministic) forcing, but with the slow envelope slightly evolving each time instead of perfectly repeating itself. In the Fourier domain, this corresponds to a random and noisy spectrum for a single realisation of the signal. By performing ensemble averages on different signal realisations, or by computing multiple periodograms of a single realisation, the random forcing spectrum converges on average towards the same spectrum as the polychromatic forcing. In practice, a single realisation of such a stochastic signal can be computed in Fourier space by multiplying the desired average spectrum with the Fourier coefficients of a Gaussian white noise. The resulting signal in real space is then very similar to the long-time solutions of the SDE recently proposed by Ewetola & Esler (Reference Ewetola and Esler2024).
Figure 2. Time and frequency domain comparison of the (a), (d) monochromatic, (b), (e) polychromatic and (c), ( f) stochastic forcing signals. Arrows in (b), (c) highlight the modulation period
$\tau$
(respectively correlation time) for the polychromatic (respectively stochastic) forcing. Solid black lines in (df) indicate the harmonic content of the momentum flux spectrum
$\hat {\mathcal{F}}= (1/2)\hat {u_0}\hat {w_0}^*$
. Grey lines in ( f) showcase a single realisation of the random forcing spectrum. The dashed red lines indicate the Gaussian spectrum (respectively average spectrum) profile for the polychromatic (respectively stochastic) forcing, for a typical value of
$\tau$
corresponding to
$\Delta \omega =0.1$
(narrower and broader spectra with
$\Delta \omega ={0.01}{-}{0.3}$
were also used, as discussed in the text). The vertical dashed blue lines represent the Brunt–Väisälä frequency.
Each simulation is first characterised by the choice of forcing type, and then by the values of two control parameters: the forcing Reynolds number
$\textit{Re}$
and the correlation time
$\tau$
(monochromatic simulations can be seen as the singular case of polychromatic forcing with
$\tau =+\infty$
). These two specific parameters have been proven to be key factors controlling the phenomenology of the horizontal wind in HLP simulations with monochromatic, polychromatic and stochastic forcings (Renaud et al. Reference Renaud, Nadeau and Venaille2019; Léard et al. Reference Léard, Lecoanet and Le Bars2020; Ewetola & Esler Reference Ewetola and Esler2024). Our study aims to unravel their impact on the mean-flow dynamics in wave-resolving simulations, thus relaxing the constraints of the HLP model.
3. Results
3.1. Comparative phenomenology of monochromatic, polychromatic and stochastic forcings
We first discuss the results of side-to-side simulations that we ran using the three forcing types described in § 2.2. The momentum flux is fixed independently of the forcing type by using a constant forcing Reynolds number
$\textit{Re}=34$
for all runs and the modulation and correlation times for the polychromatic and stochastic forcings are
$\tau =63$
. We first focus on the phenomenology of the mean flows emerging in each case. To this end, we show Hovmöller diagrams of the mean flow for three different runs in figure 3 (a–c). Both monochromatic and polychromatic forcings result in mean flows reversing periodically. The descending shear front is very sharp and reversal periods are visually indistinguishable from each other. The reversing mean flow obtained with a stochastic forcing appears more erratic, even though the momentum flux is the same in all cases. The descending front is more noisy and all mean-flow phases are distinct.
Figure 3. Hovmöller diagrams and probability density functions of the mean-flow reversal times for simulations with (a), (d) monochromatic, (b), (e) polychromatic and (c), ( f) stochastic forcings. The PDFs are computed on a single run for the mono- and polychromatic forcings, including 25 mean-flow phases. The PDF for the stochastic forcing is computed on four different runs adding up to 240 reversals.
We quantify the variability of mean-flow reversals by measuring the probability density functions (PDF) of phase periods
$\mathcal{T}$
. To compute
$\mathcal{T}$
, we first find the altitude
$z_{\textit{max}}$
where the horizontally averaged flow
$\overline {u}(z,t)$
has the greatest temporal variance. Then, we interpolate the times of rising zero crossings of the temporal signal
$\overline {u}(z_{\textit{max}},t)$
. The successive differences of zero-crossing times constitute the ensemble of phase periods
$\mathcal{T}$
.
The PDFs for
$\mathcal{T}$
are shown in figure 3 (d–f). In agreement with our earlier qualitative observation, the PDFs are identical and show a single peak at
$\mathcal{T}=2400$
for the monochromatic and polychromatic forcings. The relative standard deviation is less than 0.001, revealing a strictly periodic mean flow. This result is consistent with earlier studies that showed periodic reversals at moderate Reynolds numbers
$\textit{Re}$
, i.e. below the threshold for the transition to frequency doubling, bistability and chaos, under both monochromatic (Plumb Reference Plumb1977; Renaud & Venaille Reference Renaud and Venaille2020) and polychromatic forcing (Léard et al. Reference Léard, Lecoanet and Le Bars2020). The PDF obtained with stochastic forcing has the same average value
$\mathcal{T}=2400$
as with the other two forcing schemes. However, it is not as peaked as the other two cases and instead shows a broad distribution. Phase periods range from 1800 to 3800, such that the relative standard deviation is 0.11 and hence orders of magnitude higher than with the monochromatic and polychromatic forcings.
We reproduce this comparative experiment with different parameter values, varying the forcing intensity through
$\textit{Re}$
and
$\alpha$
, and the modulation period
$\tau$
(respectively correlation time) of the polychromatic (respectively stochastic) forcing. In all iterations of this experiment, performed in the moderate Reynolds range (i.e. below the threshold for secondary instabilities of the mean flow), the monochromatic and polychromatic forcings always produce perfectly periodic reversing mean flows. The reversals periodicity is thus very robust for those simulations, and the resulting Hovmöller diagrams are all very similar with a few exceptions at extreme parameter values, shown in figure 4, which we will discuss before moving to the case of stochastic forcing in § 3.2.
Figure 4. Peculiar Hovmöller diagrams obtained at extreme parameter values for (a) monochromatic forcing with
$\textit{Re}=68$
, (b) polychromatic forcing with
$\textit{Re}=17$
and
$\tau =630$
and (c) polychromatic forcing with
$\textit{Re}=17$
and
$\tau =21$
.
For monochromatic forcing at the highest Reynolds number value
$\textit{Re}=68$
, we observe the growth of a subharmonic modulation of the mean-flow oscillation, which is in an approximate
$1:3$
ratio considering the theoretical
$\textit{Re}^{-1}$
scaling (Renaud et al. Reference Renaud, Nadeau and Venaille2019) for the main reversal period. This observation is compatible with the frequency-locking regime, which is achieved when
$\textit{Re}$
lies between the critical values for the quasi-periodic and chaotic regimes in Renaud et al. (Reference Renaud, Nadeau and Venaille2019). For polychromatic forcing at the longest correlation time
$\tau =630$
, the wave modulation becomes slow enough to be visible on the Hovmöller diagram in the form of temporal beating at low altitudes. However, the presence of this beating does not seem to affect the periodicity of the reversals, and supplementary simulations (not shown) at slightly different values of
$\tau$
suggest that there is no particular frequency locking between the wave modulation and the mean-flow reversals, which is otherwise suggested by numerical simulations of the HLP model with additional modulation by an annual cycle (Rajendran et al. Reference Rajendran, Moroz, Read and Osprey2016).
Finally, under polychromatic forcing with the shortest correlation time
$\tau =21$
, simulations at the lowest Reynolds number value
$\textit{Re}=17$
do not trigger any oscillating wind, apart from the slow decay of the initial seed, which is invisible in figure 4(c) given the colour scale. This result suggests that, for
$\tau =21$
, the value of
$\textit{Re}=17$
becomes lower than the critical value for the emergence of the mean flow. We interpret this result in light of recent work by Chartrand, Nadeau & Venaille (Reference Chartrand, Nadeau and Venaille2024), which showed that some of the wave momentum flux is unavailable to the mean flow when the waves have multiple frequencies. Specifically, Chartrand et al. (Reference Chartrand, Nadeau and Venaille2024) showed that waves with frequencies higher than a dominant low-frequency wave do not feed the mean flow. In our case, smaller values of
$\tau$
correspond to a broader forcing band in frequency space, meaning more momentum flux is distributed at high frequencies above the centre of the wave packet. Since those frequencies do not contribute to the mean flow, the effective share of the momentum flux that contributes to the mean flow is lower, suggesting a lower effective value for the Reynolds number.
3.2. Mean-flow reversal sensitivity to stochastic parameters
We now investigate the impact of the forcing parameters on the variability of the mean-flow phases when the forcing is stochastic. We vary the input momentum flux through
$\textit{Re}$
and the forcing correlation time scale through
$\tau$
. We ran a total of 16 simulations, considering four different Reynolds numbers, i.e.
$\textit{Re}=17,34,51,68$
, and five different correlation time scales
$\tau =21,63,310,450,630$
. The range of Reynolds numbers was chosen empirically in order to keep the monochromatic simulations in the periodic reversal regime: lower values of
$\textit{Re}$
do not trigger any mean flow, and higher values trigger secondary instabilities, leading to more complex mean-flow structures, as described in Renaud et al. (Reference Renaud, Nadeau and Venaille2019). We base the phase statistics for each point of the parameter space on an aggregate of four distinct simulations, which we run in parallel in order to increase statistical convergence while reducing simulation time. Each of these four runs is subject to a different realisation of the stochastic forcing and is run until
$t \approx 160\ 000$
, thus containing 25–90 mean-flow reversals. In total, our statistics are based on 100–350 mean-flow phases for each point of the parameter space.
The results of this parametric study are shown in figure 5 (we show a subset of 9 representative cases for clarity). The nine interior frames show Hovmöller diagrams of the mean flow over (relatively) short time periods, with
$\textit{Re}$
increasing from top to bottom and
$\tau$
increasing from left to right, while the six exterior frames compare the PDFs for
$\mathcal{T}$
either row-wise or column-wise. For conciseness, we will describe in detail figures 5 (iii) and 5(vi) only, as they most effectively illustrate the impact of the correlation time scale and Reynolds number, respectively, on mean-flow statistics.
Figure 5. Hovmöller diagrams and time-to-reversal PDFs of the mean flows obtained under stochastic forcing with
$\textit{Re}=17,34,68$
(first, second and third row from top) and
$\tau =21,63,630$
first, second and third column from left). The PDFs are shown on separate axes on the fourth row and column: dotted, dashed and solid lines correspond to the first, second and third row (column) diagrams, respectively. No Hovmöller diagrams are shown in (a) and (d) as no mean flow appears in those simulations.
We first see from figure 5 (iii) that an increase in
$\textit{Re}$
leads to a decrease in the average phase duration
$\langle \mathcal{T}\rangle$
. Quantitatively, doubling
$\textit{Re}$
leads to a halving of
$\langle \mathcal{T}\rangle$
, which is a scaling well captured by the expression for the momentum deposit time
$\mathcal{T}_*=c\varLambda /\mathcal{F}\propto 1/\textit{Re}$
(Vallis Reference Vallis2017; Renaud et al. Reference Renaud, Nadeau and Venaille2019) that predicts the reversal time based on dimensional analysis. A second observation is that higher Reynolds numbers are associated with narrower PDFs for
$\mathcal{T}$
. Thus, increasing
$\textit{Re}$
reduces both the mean reversal time and range of phase periods. The relative spread of the PDFs is visible in figure 5 (vi), which shows the impact of the forcing correlation time
$\tau$
for fixed
$\textit{Re}=68$
. All PDFs have approximately the same average value, showing that an increase in
$\tau$
increases the spread of the PDFs both in absolute time units and relative to the mean. The broadening of the PDFs with the correlation time can be well understood by the statistical average argument from Ewetola & Esler (Reference Ewetola and Esler2024): when the forcing correlation time
$\tau$
is very short compared with the average phase duration
$\langle \mathcal{T}\rangle$
, the mean flow receives momentum from many uncorrelated wave packets over
$\langle \mathcal{T}\rangle$
, thus perceiving an ‘averaged’ version of the forcing. On the other hand, when
$\tau$
becomes comparable to
$\langle \mathcal{T}\rangle$
the number of wave packets seen by the mean flow over
$\langle \mathcal{T}\rangle$
decreases significantly, such that it perceives a version of the forcing that is not converged statistically, and exhibits greater variability.
Figures 5(a) and 5(d) are empty as there is no mean flow generated in these simulations. Because those runs correspond to low
$\tau$
, low
$\textit{Re}$
values, we interpret this absence of mean flow by the same spectral spreading argument as for the polychromatic case, as explained in § 3.1.
In order to further investigate the effects of a non-statistically converged wave momentum flux when using stochastic forcing, we show in figure 6 the relative standard deviation of the phase periods
$\sigma /\langle \mathcal{T}\rangle$
as a function of
$\tau /\langle \mathcal{T}\rangle$
, i.e. the ratio of the forcing correlation time to mean reversal time, for each parameter value. We see that the relative standard deviation increases mostly monotonically with
$\tau /\langle \mathcal{T}\rangle$
, suggesting that this ratio is a key parameter controlling the relative variability of reversals in our stochastic simulations. We still observe some degree of spreading around the main trend, suggesting that the rescaling proposed in figure 6 may not fully capture some effects of stochasticity in our simulations. The departure from the main trend is mainly visible for simulations with the highest levels of relative variability (i.e.
$\sigma /\langle \mathcal{T}\rangle \gt 20\%$
, as denoted by the grey area in figure 6), suggesting that this spread could arise both from physical effects and the slow statistical convergence of the standard deviation for our simulation sample. We note that the rescaling significantly reduces the vertical spread. For example, the purple cross and red triangle have a large spread in both
$\tau$
values (
$\tau =450$
and 310, respectively) and
$\sigma$
values (
$\sigma =620$
and
$370$
, respectively, i.e. 50 % difference relative to the mean) but are brought close together at
$\tau /\langle \mathcal{T}\rangle \approx 0.2$
with a reduced vertical spread in
$\sigma /\langle \mathcal{T}\rangle$
(13 % difference relative to mean) after rescaling.
Figure 6. Relative standard deviation
$\sigma /\langle \mathcal{T}\rangle$
of the mean-flow reversal period as a function of the ratio
$\tau /\langle \mathcal{T}\rangle$
, which is the forcing correlation time divided by the mean reversal time. Each symbol represents one point in the
$\tau -\textit{Re}$
parameter space. Symbol colours and styles highlight
$\tau$
and
$\textit{Re}$
values, respectively. The greyed area correspond to values of the relative standard deviation higher than
$20\,\%$
.
3.3. Mean-flow properties’ sensitivity to forcing type
We finally investigate the influence of the forcing type and parameters on the average properties of the reversing mean flow. We focus on the average reversal time
$\langle \mathcal{T}\rangle$
, the peak velocity
$u_{\textit{max}}$
of the mean flow and the altitude value
$z_{\textit{max}}$
at which peak velocity is attained. We report the results in figure 7 for all stochastic simulations, nine comparable polychromatic simulations (i.e. same points in
$\tau -\textit{Re}$
space), and three comparable monochromatic simulations, i.e. with
$\textit{Re}=17,34,68$
. We find very similar results for the average reversal time across all simulations. The Reynolds number is the key parameter controlling
$\langle \mathcal{T}\rangle$
, and we observe qualitative consistency with the HLP model scaling
$\langle \mathcal{T}\rangle \sim 1/\textit{Re}$
. For all forcing types, we find that increasing
$\textit{Re}$
increases the peak velocity
$u_{\textit{max}}$
and decreases
$z_{\textit{max}}$
, in agreement with Renaud & Venaille (Reference Renaud and Venaille2020). However, we observe a significant departure from the asymptotic scalings for single-wave streaming at
$t=+\infty$
(Renaud & Venaille Reference Renaud and Venaille2020) shown by the blacked dotted lines: the mean flow in our simulations features a jet that is systematically weaker and higher than predicted. This discrepancy is not unexpected as the mean flow is always oscillating in our simulations, whereas it is a non-oscillating, steady state in the asymptotic analysis. In fact, we expect that the mean flow in the oscillating state does not have time to reach the peak velocity value and low altitude of the steady single-wave streaming solution before it reverses. Recent deterministic two-dimensional simulations resolving waves (Renaud et al. Reference Renaud, Nadeau and Venaille2019) similarly showed mean flows weaker than in the HLP model, from which the asymptotic scalings are derived, for the same
$\textit{Re}$
values.
Figure 7. Scalings of (a) the average mean-flow reversal period
$\langle \mathcal{T}\rangle$
, (b) peak velocity
$u_{\textit{max}}$
and (c) peak velocity altitude
$z_{\textit{max}}$
with
$\textit{Re}$
. Symbol colours and styles encode the
$\tau$
values and forcing type, respectively. The dashed line in (a) represents a fit of the type
$\langle \mathcal{T}\rangle \propto 1/\textit{Re}$
. Dotted lines in (b) and (c) represent the asymptotic solutions for single-wave streaming at
$t=+\infty$
(Renaud et al. Reference Renaud, Nadeau and Venaille2019).
We also see in figure 7(b,c) that simulations under polychromatic and stochastic forcing generally feature weaker and higher altitude winds than when using monochromatic forcing. We attribute this effect to the same spectral spreading argument from Chartrand et al. (Reference Chartrand, Nadeau and Venaille2024) as already discussed in § 3.1. Indeed, both polychromatic and stochastic simulations have a proportion of the wave momentum flux above the central forcing frequency, resulting in this share not contributing to the mean-flow momentum. This effect is particularly striking for the runs at the lowest value of
$\tau$
(widest temporal spectra), where the momentum flux frequency spreading is the largest.
4. Conclusions
We have shown that internal waves forced stochastically at the bottom of a two-dimensional stably stratified fluid can generate mean flows with a broad distribution of times to reversal, at Reynolds numbers lower than the threshold for secondary instabilities of wave-driven flows. In comparison, both monochromatic and polychromatic forcings, which are deterministic, produce perfectly periodic mean flows for the same Reynolds number with no significant variability.
We have found that the forcing correlation time, relative to the mean reversal period, which is primarily controlled by the Reynolds number, strongly influences the variability of the mean flow. Our results suggest that longer correlation times are associated with a higher relative variability of the reversal period, whereas forcings with a short correlation time produce more regular reversals. We attribute this result to the degree of statistical convergence of the forcing during a single reversal time. Indeed, short correlation times allow for a cumulative contribution of momentum fluxes that is well converged at the scale of the mean flow, leading to more regular reversals. On the contrary, longer correlation times generate fewer emitted wave packets over a characteristic QBO period, resulting in higher variability of the cumulative momentum flux from one reversal to another.
Our study suggests that the forcing type has limited influence on the average properties of the mean flow. In fact, the forcing Reynolds number is the primary parameter determining the average properties of the mean flow. We have found that simulations under monochromatic forcing produce slightly stronger and lower altitudes jets. Polychromatic and stochastic simulations feature a mean flow that is driven by an effective Reynolds number that is lower than the control Reynolds number, because the wave momentum flux distributed at frequencies higher than the dominant, primary forcing wave frequency does not feed the mean flow.
Although our study of a stochastic forcing produced highly variable mean-flow reversals at the most extreme parameters values, we have not observed any reversal anomalies such as the QBO anomaly of 2016. We suggest that this is due to the statistical content of our stochastic forcing, as we have limited ourselves to Gaussian statistics. We expect that intermittent stochastic wave forcing, in the sense of more subtle higher-order statistics and their dependence on spatio-temporal scales, is a promising avenue for the reproduction of reversal anomalies. Thus, future studies should focus on designing stochastic forcing schemes with parameterisable higher-order statistics, in the same spirit as the non-Gaussian features that are modelled in the synthetic turbulence community (Pereira et al. Reference Pereira, Garban and Chevillard2016; Chevillard et al. Reference Chevillard, Garban, Rhodes and Vargas2019; Reneuve & Chevillard Reference Reneuve and Chevillard2020).
Acknowledgements
We thank C. Herbert, A. Venaille and L.-P. Nadeau for fruitful discussions. We gratefully acknowledge support from the CBPsmn (PSMN, Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. The platform operates the SIDUS solution developed by Emmanuel Quemener (Quemener & Corvellec Reference Quemener and Corvellec2013). We thank three anonymous referees for critically reading our manuscript and suggesting substantial improvements.
Funding
This work has received funding from the European Union’s Horizon research and innovation programme under grant agreement IceAblation – 101117317.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Details of the three forcing schemes
A.1. Monochromatic forcing
In the case of monochromatic forcing, the temporal signals in dimensionless form scaled by
$\omega _0$
and
$k_0$
are simply
\begin{equation} \begin{cases} u_0(t) = \sqrt {2(\mathcal{N}^2 - 1)}2\cos (t), \\[7pt] w_0(t) = \sqrt {\frac 2{\mathcal{N}^2 - 1}}2\sin (t), \\[7pt] b_0(t) = \mathcal{N}\sqrt {\frac {2\mathcal{N}^2}{\mathcal{N}^2 - 1}}2\cos (t), \end{cases} \end{equation}
where the pre-factors and cos/sin profiles are chosen to ensure that
$u$
,
$w$
,
$b$
are solutions to the linearised Boussinesq equations on the bottom boundary, as detailed in the Supplemental Material of Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2018). In particular, (2.3) and (A1) ensure that the waves have a positive vertical group velocity and energy propagates upward. Note that the waves represented by (2.3) are standing waves. They can be rewritten as the sum of counter-propagating waves along the
$x$
axis as
\begin{equation} \begin{cases} u(x,z=0,t) = \sqrt {2\mathcal{F}(\mathcal{N}^2 - 1)}\left [\cos (t - x) + \cos (t + x)\right ]\!, \\[7pt] w(x,z=0,t) = \sqrt {\frac {2\mathcal{F}}{\mathcal{N}^2 - 1}}\left [\cos (t - x) - \cos (t + x)\right ]\!, \\[7pt] b(x,z=0,t) = \mathcal{N}\sqrt {\frac {2\mathcal{F}\mathcal{N}^2}{\mathcal{N}^2 - 1}}\left [-\sin (t - x) + \sin (t + x)\right ]\!. \end{cases} \end{equation}
The directional wave momentum flux, say in the positive
$x$
direction, can then be readily computed as
\begin{equation} \mathcal{F}^+ = \langle {u^+w^+}\rangle _{z=0} = 2\mathcal{F}\langle \cos ^2(t - x)\rangle = \mathcal{F}, \end{equation}
with
$u^+$
(respectively
$w^+$
) the horizontal (respectively vertical) velocity of the right-propagating wave, i.e. which is the first term on the right-hand side of the top (respectively middle) line of (A2). Note that the absolute wave momentum fluxes to the right and left are equal such that we may write
$|\mathcal{F}^+|=|\mathcal{F}^-|=\mathcal{F}$
.
A.2. Polychromatic forcing
The polychromatic forcing is simply constructed as the sum of multiple phase-shifted monochromatic waves, i.e. such that
\begin{equation} \begin{cases} u_0(t) = \frac {1}{\mathcal{F}}\sum _{p=0}^{n-1}\sqrt {2\mathcal{F}_p\frac {\mathcal{N}^2 - {\varpi _p}^2}{{\varpi _p}^2}}2\cos ({\varpi _p}t + \varphi _p), \\[7pt] w_0(t) = \frac {1}{\mathcal{F}}\sum _{p=0}^{n-1}\sqrt {2\mathcal{F}_p\frac {{\varpi _p}^2}{\mathcal{N}^2 - {\varpi _p}^2}}2\sin ({\varpi _p}t + \varphi _p), \\[7pt] b_0(t) = \frac {\mathcal{N}}{\mathcal{F}}\sum _{p=0}^{n-1}\sqrt {2\mathcal{F}_p\frac {\mathcal{N}^2}{\mathcal{N}^2 - {\varpi _p}^2}}2\cos ({\varpi _p}t + \varphi _p), \\ \end{cases} \end{equation}
where
$\mathcal{F}_p$
is the directional momentum flux spectrum at frequency
$\varpi _p$
, and the phases
$\varphi _p$
are uniformly distributed between 0 and
$2\pi$
. The sums are performed over reduced frequencies
$\varpi _p = 2\pi p/\tau$
where
$\tau$
is the forcing correlation time, such that the forcing scheme defined by (A4) is
$\tau $
-periodic. The number of forced modes
$n$
is such that
$\varpi _{n-1}\leqslant \mathcal{N}$
and
$\varpi _{n}\gt \mathcal{N}$
.
The wave momentum flux spectrum follows a Gaussian profile centred on
$\varpi =1$
, with cutoffs at both very low frequency and the Brunt–Väisälä frequency, such that
\begin{equation} \mathcal{F}_p = \mathcal{A}\tanh ^4\left (\frac {\varpi _p}{\varpi _m}\right )\exp\! \left [\frac {-(\varpi _p - 1)^2}{2\Delta \varpi ^2}\right ]\exp\! \left (\frac {-\mathcal{N}^2}{\mathcal{N}^2-{\varpi _p}^2}\right )\!, \end{equation}
where
$\varpi _m=0.1$
is the low cutoff frequency,
$\Delta \varpi =2\pi /\tau$
is the spectrum bandwidth and the coefficient
$\mathcal{A}$
is adjusted so that
$\sum _p\mathcal{F}_p=\mathcal{F}$
. It should be noted that the actual choice for
$\mathcal{F}_p$
is mostly independent of the
$\tau$
-periodicity constraint of the forcing time-series. Indeed, it is sufficient for this to build the forcing spectrum as a set of harmonically related components each separated in frequency by
$2\pi /\tau$
, but it is not needed for the gravest harmonics to contain a significant amount of energy, as shown in figure 2. The total directional wave momentum flux is then computed by summing over all frequencies of e.g. the right-propagating wave function, i.e.
\begin{equation} \mathcal{F}^+ = \sum _{p=0}^{n-1}\langle {u_p^+w_p^+}\rangle _{z=0} = \sum _{p=0}^{n-1}2\mathcal{F}_p\langle \cos ^2(\varpi _pt - x)\rangle = \mathcal{F}. \end{equation}
A.3. Stochastic forcing
For the stochastic forcing scheme we model
$u_0$
,
$w_0$
,
$b_0$
using one-dimensional random fields, similar to the ones used in random turbulence modelling (Pereira et al. Reference Pereira, Garban and Chevillard2016; Chevillard et al. Reference Chevillard, Garban, Rhodes and Vargas2019; Reneuve & Chevillard Reference Reneuve and Chevillard2020). The three random fields are the result of linear convolutions involving a Gaussian white noise
$W$
, i.e. such that
\begin{equation} \begin{cases} u_0(t) = \int _{s\in \mathbb{R}}K_u(t-s)W({\mathrm{d}}{s}), \\[5pt] w_0(t) = \int _{s\in \mathbb{R}}K_w(t-s)W({\mathrm{d}}{s}), \\[5pt] b_0(t) = \int _{s\in \mathbb{R}}K_b(t-s)W({\mathrm{d}}{s}), \\ \end{cases} \end{equation}
with
$K_u$
,
$K_w$
and
$K_b$
the convolution kernels, which we can define using inverse Fourier transforms as
\begin{equation} \begin{cases} K_u(t) = \int _{\varpi \in \mathbb{R}}\sqrt {2\hat {\mathcal{F}}(\varpi )\frac {\mathcal{N}^2 - {\varpi }^2}{{\varpi }^2}}e^{i\varpi t}\frac {{\mathrm{d}}\varpi }{2\pi }, \\[5pt] K_w(t) = \int _{\varpi \in \mathbb{R}}\sqrt {2\hat {\mathcal{F}}(\varpi )\frac {{\varpi }^2}{\mathcal{N}^2 - {\varpi }^2}}e^{i\varpi t}\frac {{\mathrm{d}}\varpi }{2\pi }, \\[5pt] K_b(t) = \int _{\varpi \in \mathbb{R}}\mathcal{N}\sqrt {2\hat {\mathcal{F}}(\varpi )\frac {\mathcal{N}^2}{\mathcal{N}^2 - {\varpi }^2}}e^{i\varpi t}\frac {{\mathrm{d}}\varpi }{2\pi }. \\ \end{cases} \end{equation}
The directional wave momentum flux spectrum
$\hat {\mathcal{F}}(\varpi )$
can then be computed using a continuous version of (A5).
In practice, a numerical approximation of (A7) considers a discrete convolution between the sampled kernels and a discrete-time Gaussian white noise, which can be computed efficiently as a product in Fourier space. Using this numerical approximation, the forcing signals must be sampled from a time window larger than the total simulation time to avoid signal periodicity. Here, we typically use a forcing window of 200 000 time units, i.e. roughly twice as long as our simulation time.
In the case of stochastic forcing, both the forcing signals and hence directional wave momentum flux are random processes. However, the statistics of
$u_0$
and
$w_0$
are stationary in time, such that the expected value of
$\mathcal{F}^+$
can be computed, i.e. as
\begin{equation} \mathbb{E}[\mathcal{F}^+] = \frac 12\mathbb{E}[u_0w_0] = \int _{s\in \mathbb{R}}K_u(s)K_w(s){\mathrm{d}}{s} = \int _{\varpi \in \mathbb{R}}\hat {\mathcal{F}}(\varpi )\frac {{\mathrm{d}}\varpi }{2\pi } = \mathcal{F}. \end{equation}
Open access
