Crucial role of obliquely propagating gravity waves in the quasi-biennial oscillation dynamics

​​​​​​​In climate modelling, the reality of simulated flows in the middle atmosphere is largely affected by the model's representation of gravity wave processes that are unresolved, while these processes are usually simplified to facilitate computations. The simplification commonly applied in existing climate models is to neglect wave propagation in horizontal direction and time. Here we use a model that fully represents the propagation of unresolved waves in all directions, thereby elucidating its dynamical effect upon the most important climate mode in the tropical stratosphere, i.e. the quasi-biennial oscillation. Our simulation shows that the waves in the equatorial stratosphere, which are known to drive this climate mode, can originate far away from the Equator in the troposphere. The waves propagating obliquely toward the Equator are found to play a huge role in the phase progression of the quasi-biennial oscillation as well as in its penetration into the lower stratosphere. Such waves will require further attention, given that current climate models are struggling to simulate the quasi-biennial oscillation down to the lower stratosphere, which may be needed to reproduce its observed impacts on the surface climate.

1 Introduction

Atmosphere models simulate flows on scales bounded by their resolution, while the effects of smaller-scale unresolved processes on the simulated flows are taken into account by additional formulations, so-called parameterizations, based on our knowledge of such processes. In climate modelling, atmospheric gravity waves (GWs), an internal wave mode with horizontal wavelengths of about 1–1000 km, are subject to parameterization. They play a pivotal role in large-scale circulations and their variability in the stratosphere and above (Fritts and Alexander, 2003; Kim et al., 2003). Their most important process in this regard is transporting momentum from the troposphere to upper layers through wave propagation, and therefore GW parameterizations primarily are to represent this process. As a simplification, existing GW parameterizations conventionally consider the wave propagation to be purely vertical and steady in time (e.g. Lindzen, 1981; Warner and McIntyre, 1999; Hines, 1997; Scinocca, 2003), while in the real atmosphere, the propagation is oblique and transient. The effects of this usual simplification on modelled atmospheric circulations and climate variability are however not well known.

The quasi-biennial oscillation (QBO) (Ebdon and Veryard, 1961; Baldwin et al., 2001) is the prominent climate mode of the tropical stratosphere. It is characterized by persistent alternations of the flow direction between easterly and westerly, which are driven by momentum transported primarily by GWs (e.g. Dunkerton, 1997; Kawatani et al., 2010; Ern et al., 2014; Kim and Chun, 2015). This oscillation also propagates downward to the tropopause layer and has a broad impact on atmospheric circulations such as the stratospheric polar vortex (Holton and Tan, 1980), extratropical surface climate (Marshall and Scaife, 2009; Gray et al., 2018), and tropical convection (Gray et al., 2018; Haynes et al., 2021; Yoo and Son, 2016). The atmospheric modelling community has strived to reproduce the QBO in climate simulations and seasonal predictions (e.g. Butchart et al., 2018; Richter et al., 2020; Coy et al., 2022). Currently, many climate models are able to simulate this oscillation with reasonable periods, using GW parameterizations tuned to supply the required momentum forcing. However, the models exhibit a common bias, i.e. a significant underestimation of the QBO easterly magnitude in the lower stratosphere (Bushell et al., 2022; Anstey et al., 2022a). Probably related to this deficiency, climate models could not properly reproduce the aforementioned tropospheric impacts of the QBO (Anstey et al., 2022b; Martin et al., 2023). Moreover, the simulated QBO shows large deviations among models in its spatial structure and future evolution (Richter et al., 2020, 2022). This discrepancy as well as the common bias in the current models may reflect a lack of our knowledge in detailed dynamics of the QBO.

Here we perform a climate simulation of the QBO using a unique GW parameterization, the Multi-Scale Gravity Wave Model (MS-GWaM; see Sect. 2.2), newly developed to represent the three-dimensional and transient wave propagation (referred to as the 3d-TR experiment). The simulation result is compared to a control experiment in which the conventional simplification of GW parameterization (purely vertical and steady propagation) is applied (the 1d-ST experiment). Our results, for the first time, present the role of obliquely propagating GWs in the QBO dynamics that has been veiled by the usual simplification of existing parameterizations. These waves are found to provide momentum forcing required especially for the descent and amplification of the easterly QBO phase in the lower stratosphere, where the aforementioned common bias of climate models exists.

2 Methods

2.1 Experimental design

All the experiments use a common setup, except for the use of simplifications in the GW parameterization. The ICOsahedral Non-hydrostatic model (ICON) (Zängl et al., 2015), the German operational modelling system for numerical weather and climate predictions, is used (version 2.6.4-nwp5). For the study, we replace its original non-orographic GW parameterization (Scinocca, 2003; Orr et al., 2010) with the newly developed three-dimensional transient parameterization MS-GWaM (Sect. 2.2). In addition, a fourth-order vertical damping of divergence is implemented (added to the horizontal damping of divergence established in ICON) instead of using the existing second-order background vertical diffusion. Suppressing the latter is found to be beneficial in simulating the QBO with less artificial vertical damping in the stratosphere.

The experiments are performed with climatological-mean annual-cycle forcing (e.g. ozone and sea-surface temperature) for recent decades, for the purpose of simulating mean characteristics of the QBO over its cycles (rather than capturing its variations among the cycles). Each simulation is for 20 years after about 2 years of a spin-up period. We use a horizontal grid spacing of ∼ 160 km (20 480 horizontal grid cells) with 180 vertical layers up to an altitude of 120 km. A sponge-layer damping is applied from 85 km upward. The vertical grid spacing is constantly 400 m from the mid-troposphere to the mid-stratosphere (36 km) and slowly increases above (reaching ∼ 1.2 km at the sponge-layer bottom).

The experiments of the study differ only in the GW parameterization, one fully representing the three-dimensional, transient wave propagation (3d-TR experiment) and another applying the conventional simplifications that have been used in climate models, i.e. representing only the vertical propagation with the steady-state assumption (1d-ST experiment). Additionally, an experiment with the one-dimensional but transient parameterization (1d-TR experiment), as an intermediate-level simplification, is also performed and briefly explained in Appendix B. The different treatments in the wave propagation modelling in these experiments are described in Sect. 2.2 and 2.3. It should be noted that, for the wave generation in the troposphere, the parameterized wave spectra are virtually the same for all the experiments in the climatological mean (refer to Fig. A1), and therefore any differences in the simulated QBO between the experiments are due entirely to the wave propagation modelling.

2.2 Gravity-wave parameterization: three-dimensional

A GW parameterization that models three-dimensional transient wave dynamics, MS-GWaM, has recently been developed using a Lagrangian ray-tracing approach and implemented in ICON (Bölöni et al., 2021; Kim et al., 2021; Voelker et al., 2023). Its detailed theoretical basis can be found in Achatz (2022) and Achatz et al. (2023). Below we briefly describe its governing equations for modelling the wave propagation.

For GWs at a position x and time t, their frequencies ω and wavenumbers k obey the following dispersion relation:

$$\begin{array}{}\text{(1)}& \mathit{\omega }=\mathit{U}\cdot \mathit{k}+\sqrt{{\displaystyle \frac{N^{\mathrm{2}}|{\mathit{k}}_{\mathrm{h}}{|}^{\mathrm{2}}+f^{\mathrm{2}}(k_z^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})}{|\mathit{k}{|}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}}}}\equiv \mathrm{\Omega }\left(\mathit{k},\mathit{x},t\right),\end{array}$$

with k_h and k_z being respectively the horizontal and vertical components of k, where f is the Coriolis parameter, and all the flow variables, i.e. horizontal wind U, Brunt–Väisälä frequency N, and pseudo-incompressible scale-height parameter Γ⁻¹, are functions of (x,t). Ω is defined as a function that expresses ω in terms of $(\mathit{k},\mathit{x},t)$ given the dispersion relation. The equations for modelling wave propagation consist of the ray equations

$$\begin{array}{}\text{(2)}& \left(\dot{\mathit{x}},\dot{\mathit{k}}\right)=\left({\mathrm{\nabla }}_{\mathit{k}}\mathrm{\Omega },-{\mathrm{\nabla }}_{\mathit{x}}\mathrm{\Omega }\right)\end{array}$$

to predict the position and wavenumber changes following GW rays as well as the equation for wave-action density $\mathcal{N}(\mathit{k},\mathit{x},t)$ in the six-dimensional phase space spanned by x and k:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{D}\mathcal{N}}{\mathrm{D}t}}=\left({\displaystyle \frac{\partial }{\partial t}}+\dot{\mathit{x}}\cdot {\mathrm{\nabla }}_{\mathit{x}}+\dot{\mathit{k}}\cdot {\mathrm{\nabla }}_{\mathit{k}}\right)\mathcal{N}=\mathcal{S}.\end{array}$$

The wave-action density is conserved in that space up to the source or sink 𝒮 arising from wave generation or dissipation.

In the parameterization, the wave-action field is discretized spatially and spectrally into finite volumes in the phase space (so-called ray volumes), and Eqs. (2) and (3) are solved for each ray volume in a Lagrangian manner. From the predicted 𝒩 field, all the fields that are required to calculate the wave effects on the model flow (such as momentum fluxes and forcing presented in Sect. 3.2 and 3.3) can be derived. Details of the discretization and the calculation of wave effects as well as the wave dissipation modelling can be found in Voelker et al. (2023). In the 3d-TR experiment, we use about 40 000 ray volumes per model grid column and time at most for accurate modelling.

The tropical source of waves taken into account by the parameterization is cumulus convection, which is also parameterized, independently, by ICON's subgrid cumulus scheme (Bechtold et al., 2008). The formulation of convectively generated GW spectra and its implementation in our parameterization for the source of 𝒩 generally follow Song and Chun (2005) and Kim et al. (2021) respectively. A difference exists in the present implementation compared to that work, as documented in Appendix A.

2.3 Gravity-wave parameterization: one-dimensional

The one-dimensional transient parameterization (Bölöni et al., 2021; Kim et al., 2021), which neglects the horizontal propagation, uses the same equations and methods as those described in Sect. 2.2, except for applying ${\dot{\mathit{x}}}_{\mathrm{h}}={\dot{\mathit{k}}}_{\mathrm{h}}=\mathrm{0}$ to the equations (where x_h denotes the horizontal position of a wave). We use the same number of ray volumes in the 1d-TR experiment as in the 3d-TR experiment (∼ 40 000 per model-grid column and time at most).

From the one-dimensional equations, the steady-state approximation is further applied in the 1d-ST experiment, neglecting local time derivatives. Denoting the vertical group velocity $c_{gz}=\dot{z}$ (with z being the vertical coordinate) and using a general property of rays in phase space (${\mathrm{\nabla }}_{\mathit{x}}\cdot \dot{\mathit{x}}+{\mathrm{\nabla }}_{\mathit{k}}\cdot \dot{\mathit{k}}=\mathrm{0}$), Eq. (3) reduces to a diagnostic equation

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial }{\partial z}}\left\{c_{gz}\mathcal{N}\right\}=\left\{\mathcal{S}\right\}\end{array}$$

by integration over k_z for a given k_h at z, where {⋅} denotes the integral. The widely used equation form in conventional GW parameterizations, which is also used in our 1d-ST experiment, is obtained accordingly as

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial \left\{{\mathcal{F}}_{\mathrm{p}}\right\}}{\partial z}}={\mathit{S}}_{\mathrm{p}}\end{array}$$

by defining pseudo-momentum 𝓟=k_h𝒩 with its vertical flux 𝓕_p=c_gz𝓟, where S_p={k_h𝒮} is the source or sink of pseudo-momentum. Therefore, the parameterization with the one-dimensional steady-state approximation reduces to modelling wave sources and sinks at every horizontal position and time.

2.4 Reanalysis data

As a reference for the tropical stratospheric wind field, the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Reanalysis (ERA-Interim, Dee et al., 2011b) is utilized. Its good agreement with observations regarding the QBO has been well documented (e.g. SPARC, 2022). The native model-level product of the 6-hourly zonal wind is used, which essentially employs the pressure vertical coordinates in the stratosphere. For comparison to the ICON simulation results, we convert the pressure coordinates to the altitude coordinates by approximation using a scale height of 6.3 km and a reference altitude of 18.6 km at 70 hPa.

3 Results

3.1 Modelled structure of the QBO

The vertical and latitudinal profiles of the QBO winds in the simulations are shown in Figs. 1 and 2 respectively, along with those in ERA-Interim (ERA). In the vertical profiles (Fig. 1), a couple of differences are found between the two experiments: (i) periods of the oscillation are much longer in 1d-ST (36–48 months) than those in 3d-TR (24 months), and (ii) the downward propagation of easterly phases is less pronounced in 1d-ST, exhibiting slower descents and weaker easterly amplitudes between ∼ 27 and 19 km. Westerly phases, on the other hand, show comparable speeds of descent between the experiments until the descents halt, while afterwards they are prolonged at ∼ 21 km in 1d-ST until the easterly phases above penetrate down to this altitude. The contrast in the simulated QBO periods therefore results from the different speeds of easterly-phase progression. Compared to ERA, the periods and peak amplitudes of the QBO are overall well reproduced in 3d-TR, while the easterly jets tend to be a bit weaker at 21–24 km.

Figure 1 Time series of vertical profiles of the tropical stratospheric zonal winds averaged over 5° N–5° S in the two experiments respectively using the three-dimensional transient gravity-wave parameterization (3d-TR) and using the parameterization simplified by the conventional (one-dimensional steady-state) approximation (1d-ST), along with those in ERA-Interim (ERA) for 20 years. The winds have been averaged monthly and zonally. The simulations are designed to represent the climate of recent decades around the year 2000, and accordingly the time series in ERA are plotted for the period centred on the decade of the 2000s (1996–2015).

Figure 2 As in Fig. 1 but for the time series of latitudinal profiles at 24 km altitude.

The latitudinal profiles of the winds exhibit another notable difference between the experiments. As found above, the easterly QBO phases penetrate well down to altitudes below 27 km in 3d-TR. Accordingly, the wind structure with alternating directions around the Equator is reproduced at 24 km, in agreement with that in ERA (Fig. 2). In 1d-ST, in contrast, as the equatorial QBO easterlies are too weak, peak easterlies appear 10–20° off the Equator in the summer hemisphere (i.e. the Southern Hemisphere at the beginning of each year and the Northern Hemisphere 6 months later). Furthermore, their magnitudes are overestimated by ∼ 10 m s⁻¹ compared to those in ERA and 3d-TR at the same locations. The result in Figs. 1 and 2 demonstrates that the simplified representation of GW propagation can lead to different latitudinal and vertical structures of the tropical stratospheric flow in climate simulations.

3.2 Oblique propagation of gravity waves

For an interpretation of the above findings, we first examine GW propagation in 3d-TR. Since the major differences in the QBO characteristics between the two experiments are associated with the easterly-phase descents (Fig. 1), we focus on the easterly momentum carried by GWs, which is responsible for these descents. Figure 3 presents horizontal fields of upward fluxes of easterly momentum ($-\mathit{\rho }\langle w^{\prime }{u}^{\prime }\rangle$ due to waves with ${\widehat{c}}_{\mathit{\lambda }}<\mathrm{0}$, where ρ is the density, the bracketed term is the covariance of the wave perturbations of vertical and zonal winds, and ${\widehat{c}}_{\mathit{\lambda }}$ is the intrinsic zonal phase speed) at altitudes of 14 and 24 km (filled and open contours respectively). Here, only the GWs generated by tropical convection occurring in a 1 h time window on a day are taken into account as an example. These target waves have been decomposed into three groups in Fig. 3 based on their initial horizontal wavelengths (i.e. wavelengths at generation). The flux fields are integrated in time, so that they would be approximately conserved up to wave dissipation between the two altitudes if the wave propagation were purely vertical. Therefore, changes in the horizontal distribution of the fluxes with altitude indicate oblique propagation of the waves, possibly along with the wave dissipation effect. In particular, the waves with horizontal wavelengths larger than 300 km observed over Africa at 14 km altitude are found to propagate southwestward by up to about 15° until they reach 24 km altitude (Fig. 3). In contrast, waves with wavelengths smaller than 300 km travel much less in horizontal directions (≲ 5°), which are mostly westward.

Figure 3 Horizontal fields of time-integrated upward fluxes of easterly momentum (contoured at 0.5 mPa h) due to gravity waves parameterized in the 3d-TR experiment at two altitudes, 14 km (blue, filled) and 24 km (red, open), for comparison. Only the waves that are generated during a certain time window (for 1 h on a day in June) are taken into account to trace the given waves' displacement, and they are decomposed based on the horizontal wavelengths at their generation (λ₀): λ₀< 100 km, 100 km $\le {\mathit{\lambda }}_{\mathrm{0}}<$ 300 km, and 300 km $\le {\mathit{\lambda }}_{\mathrm{0}}<$ 1000 km (from top to bottom). The fluxes are integrated over a period long enough (4 d) to cover the entire wave propagation up to the 24 km altitude.

Such equatorward-slanted propagation over considerable distances as seen for the case in Fig. 3 occurs preferentially at a particular phase of the QBO but persistently in every QBO cycle throughout the 20-year simulation period. Figure 4 shows zonally averaged upward fluxes of easterly momentum due to GWs (shading) along with zonal winds (contours), composited for each QBO phase during the 20 years. As the QBO in 3d-TR has regular 2-year periods (Fig. 1), we define its phases simply by eight consecutive 3-month periods for each QBO cycle, such that the nominal first phase (Phase 1) corresponds to the period with the maximum easterly wind being located at about 31 km (∼ 10 hPa) altitude (refer to the zonal-wind fields in Fig. 4). Phase 1 corresponds to February–April of every other year in 3d-TR.

Figure 4 Composite mean of zonally averaged easterly-momentum fluxes due to gravity waves (shading) along with zonal winds (blue contours with dashed lines for easterly winds and solid lines for zero and westerly winds, at intervals of 5 m s⁻¹) in the 3d-TR experiment, for each QBO phase (Phase 1 to Phase 8). In this experiment, the QBO phases are defined by consecutive 3-month periods, with the first one (Phase 1) being the period when the maximum easterly wind is located at ∼ 31 km altitude (∼ 10 hPa). Phase 1 corresponds to February–April of every other year (note that the periods of the QBO in the 3d-TR experiment are 2 years regularly).

In general, the easterly-momentum fluxes in the tropical upper troposphere (∼ 15 km) are broadly distributed with latitude, and their maxima are often located off the Equator (Fig. 4) following the seasonal dependence of convection. In the stratosphere, the momentum fluxes tend to decrease with altitude in westward sheared layers throughout the QBO cycle, due to wave dissipation. Oblique propagation of waves is manifested during Phases 2–3 by a meridionally slanted structure of the fluxes. In particular, the equatorward propagation can be identified, originating from around 10° N in the upper troposphere (as also observed for the example given in Fig. 3). In the slanted structure, the increase in flux with altitude up to ∼22 km over the Equator is attributed to equatorward wave propagation, while the flux decrease at latitudes higher than ∼ 8° N additionally involves wave dissipation due to critical-level filtering by the easterly flow with magnitudes of ∼ 10 m s⁻¹.

The propagation path of waves is controlled by their ambient wind structure, which the QBO modulates, as well as by their own characteristics (Lighthill, 1978). Our simulation shows that waves carrying easterly momentum (waves with ${\widehat{c}}_{\mathit{\lambda }}<\mathrm{0}$) tend to propagate obliquely toward the Equator when the ambient flow is weakly easterly in the upper troposphere to lower stratosphere, as during Phases 2–3 presented in Fig. 4. This condition is satisfied when the QBO easterly is maximal in the middle stratosphere (Fig. 4). On the other hand, the waves propagate more vertically in westerly ambient flows (e.g. during Phases 8–1) or tend to dissipate in vertically sheared flows when the easterly QBO phase has descended to the lower stratosphere (Phases 5–7). These behaviours are qualitatively consistent with the theory that GWs are modulated to propagate more vertically than horizontally where the ambient flow velocity backs away from the phase velocities of the waves, i.e. where the intrinsic phase velocities are increased (e.g. Lighthill, 1978). It is found from further investigations that GWs travelling long distances toward the Equator in the lower stratosphere generally have horizontal wavelengths larger than about 300 km (not shown).

The westerly-momentum fluxes ($\mathit{\rho }\langle w^{\prime }{u}^{\prime }\rangle$ due to waves with ${\widehat{c}}_{\mathit{\lambda }}>\mathrm{0}$) composited for the eight QBO phases are shown in Fig. 5 for completeness. Here, the fluxes peaking off the Equator do not exhibit the slanted structure. This may be attributed to the easterly ambient flows in the stratosphere off the Equator, which modulate the waves that carry westerly momentum to propagate more vertically.

Figure 5 As in Fig. 4 but for westerly-momentum fluxes.

3.3 Effect of oblique wave propagation on the QBO

The persistent occurrences of the equatorward propagation during Phases 2–3 suggest that it may robustly play a role in the QBO dynamics. Here we investigate the GW forcing of the QBO around Phases 2–3 in 3d-TR and 1d-ST. Phase 2 in 3d-TR corresponds to May–July of every other year where the easterly-maximum wind is located at around 28 km (see Sect. 3.2 and Fig. 4). Consistently, Phase 2 in 1d-ST is defined as the 3-month period with the easterly maximum being located around this altitude in each QBO cycle. For comparison to 3d-TR, only those cycles where Phase 2 corresponds to May–July (three cycles out of five in 1d-ST) are considered in the following composite analysis. Figure 6 shows the easterly-momentum fluxes and zonal-wind forcing due to GWs (shading and green contour respectively) during Phase 2 (centre panels), along with those over the consecutive 3-month periods before and after Phase 2 (left and right panels respectively), in 3d-TR and 1d-ST (upper and lower panels respectively).

Figure 6 As in Fig. 4 but for Phase 2 in the 3d-TR and 1d-ST experiments (panels b and e respectively), with zonal-mean zonal momentum forcing due to gravity waves (green contour, at −0.2 m s⁻¹ d⁻¹) being superimposed. In the 1d-ST experiment, Phase 2 is defined as the 3-month period where the easterly maximum is located at about 28 km in each QBO cycle, but only the cycles in which Phase 2 corresponds to May–July are composited here, so that the features in the same season and phase are compared between the two experiments. The consecutive 3-month periods before and after Phase 2 are shown in panels (a), (d), (c), and (f) respectively in each experiment (which, in the 3d-TR experiment, are Phases 1 and 3, as in Fig. 4). The numbers of the composited QBO cycles are 10 and 3 in the 3d-TR and 1d-ST experiments respectively.

In 1d-ST, by construction, the wave propagation is purely vertical, and the momentum fluxes only decrease with altitude where waves dissipate. The GW forcing of zonal winds typically occurs where the vertical gradient of the flux is large. In 3d-TR, 3 months before Phase 2, the overall distribution of the momentum fluxes and forcing is similar to that in 1d-ST (Fig. 6a and d). However, during Phase 2 when the equatorward wave propagation is manifested, the momentum fluxes at about 24 km altitude exhibit their maximum around the Equator, and they strongly dissipate higher up due to the large shear associated with the equatorial QBO jet (Fig. 6b). This induces substantially large easterly-momentum forcing below the easterly-maximum altitude, thereby leading to the descent of the easterly maximum afterwards (see Fig. 6c). This behaviour is in strong contrast with 1d-ST, where the momentum forcing occurs off the Equator with a weaker magnitude during Phase 2 and therefore the easterly descent is much slower. This result demonstrates that the descent of the easterly QBO phase is largely affected by the wave propagation path, explaining the differences in the speed of the descent and vertical penetration between 3d-TR and 1d-ST shown in Figs. 1 and 2.

GWs can also influence large-scale tropical waves explicitly resolved by the model (e.g. Kim and Achatz, 2021). Therefore, changes in the QBO between 3d-TR and 1d-ST could, in part, be attributed to the changes in the resolved waves due to the different representations of GW propagation. However, we confirmed that the resolved wave forcing (measured by the Eliassen–Palm flux divergence) is an order of magnitude smaller than the zonal-mean GW forcing in the westward sheared layer in 3d-TR (not shown). This implies a relatively minor impact of the resolved waves on the descent of the easterly QBO phase in our model, compared to the GW forcing shown in Fig. 6. On the other hand, in the eastward sheared layer, equatorial Kelvin waves induce westerly-momentum forcing with a magnitude comparable to that by GWs. Interaction of Kelvin waves with obliquely propagating GWs may merit further investigation, although our results do not show a distinct difference in the descent of the westerly phases, on average, between 3d-TR and 1d-ST (Fig. 1).

3.4 Comparison to the one-dimensional parameterization being tuned

Our results show that, via oblique propagation, waves that originate off the Equator provide the equatorial stratospheric flow with momentum which significantly accelerates the QBO. In climate modelling with conventional one-dimensional GW parameterizations, a practical and general approach to accelerate the QBO has been to empirically enhance the magnitude of the momentum flux of waves at their launch locations over the Equator so that the required momentum can be supplemented above in the stratosphere. While a reasonable timescale of the oscillation could be acquired by this approach, the spatial structures of the modelled flows should be examined in comparison to those resulting from the realistic oblique wave propagation.

Following the approach, we repeat the 1d-ST simulation but with GW fluxes increased by 50 % (empirically determined) at launch locations and compare its result (Fig. 7) to 3d-TR. The periods of the QBO in this experiment are modelled to be 2–3 years as intended due to fast descents of easterly phases (Fig. 7). However, the easterly phases descend less in depth, with shorter phase durations than those in 3d-TR, while westerlies are too strong above ∼ 23 km (see Fig. 1). In the summer hemisphere in the lower stratosphere, the excessive easterly bias found in 1d-ST still remains with similar magnitudes (see Fig. 2). The discrepancy of these results from 3d-TR reflects the fact that the oblique equatorward propagation of waves in 3d-TR occurs and accelerates the QBO preferentially during the easterly-descending phase of the oscillation, whereas the simple tuning in the one-dimensional parameterization accelerates or amplifies the entire phases and also over-accelerates flows off the Equator (e.g. in the summer hemisphere) to model reasonable periods of the equatorial oscillation. Given this physical reason, it is convincing that such a discrepancy would remain even if another climate model was used for the current study, although some quantitative details would change. In addition, to mimic the effects of realistic wave propagation somehow using a one-dimensional parameterization, its tuning will need to be designed in a sophisticated way, based on the understanding of actual processes of GWs.

Figure 7 Time series of vertical profiles (a) and 24 km latitudinal profiles (b) of the tropical stratospheric zonal winds in the experiment using 1d-ST parameterization but tuned by raising the launching fluxes of gravity waves by 50 % in order to obtain realistic periods of the oscillation (2–3 years).

4 Discussion

Although not presented in this study, we note that ICON with its original GW parameterization (Scinocca, 2003; Orr et al., 2010), which uses a spatiotemporally uniform, prescribed wave spectrum and a different wave-dissipation scheme from that used here, simulated the QBO with generally weak amplitudes when the same experimental setup as described in Sect. 2.1 was used.

The obliquely propagating waves that significantly affect the QBO in 3d-TR have horizontal wavelengths of 300–1000 km with variable vertical wavelengths down to ∼ 1 km. Waves on these scales are subject to parameterization, as they are not fully resolved by current climate models due to the limitation in horizontal and vertical resolutions as well as the difficulty in properly generating the wave source (multi-scale convection, such as mesoscale convective systems). In our simulation, the waves on those scales account for only about 10 % of the parameterized GW spectrum in the tropics (see Fig. A1 for the spectrum). Given their large effects on the QBO (Figs. 1 and 2) even with the relatively small contribution to the spectrum, quantitative observational investigations of them will be required to better understand and model the QBO. It may still be improbable to explicitly capture three-dimensional GW propagation using current measurement techniques. Nonetheless, a recent observational campaign (Haase et al., 2018) produced statistics showing that a substantial portion of tropical GWs detected in the lowermost stratosphere (∼ 20 km) had their sources at far horizontal distances (∼ 10 °) in the troposphere (Corcos et al., 2021), which supports our simulation result of oblique propagation.

It is especially in the easterly-descending QBO phase in the lower stratosphere where the effect of obliquely propagating waves is large in our simulation (Fig. 6), but this effect could be even larger depending on the quantitative details of the waves. The oblique wave propagation process is therefore a strong hint of the aforementioned common model bias of the lower-stratospheric QBO easterlies, which needs to be corrected to reproduce the observed downward impact of the QBO on the surface climate (Anstey et al., 2022b). Finally, it should be highlighted that the QBO projection on a changing climate, which was not robustly simulated among models and/or GW parameterizations (Richter et al., 2022; Schirber et al., 2015), may be more reliable using a three-dimensional GW parameterization because the wave propagation features vary depending on flow structures under the changing climate.

Appendix A: Parameterized wave spectrum at generation

As documented in Sect. 2.2, the formulation and implementation of convectively generated GW spectra follow Song and Chun (2005) and Kim et al. (2021). In the present implementation, however, a notable difference exists from that work. While for the horizontal and temporal scales of convective latent heating (δ_h, δ_t), which are preset parameters used in the source formulation, a single scale set has been taken there (5 km and 20 min for the horizontal and temporal scales respectively), here a distinctly larger-scale set (100 km, 12 h) is used in addition in order to take the multi-scale nature of tropical convection into account. The latter scale is chosen as a representative scale of the convective heating distribution in mesoscale convective systems that are unresolved by climate models (e.g. Tao and Moncrieff, 2009), and it is found to be important to generate waves that have wavelengths larger than ∼ 300 km in our simulations (refer to Trinh et al., 2016). The calculated spectrum at wave generation using those two scale sets, averaged over the tropics for the whole simulation period, is presented in Fig. A1.

Figure A1 Horizontal-wavelength spectrum of vertical fluxes of absolute horizontal momentum due to gravity waves launched at 15° N–15° S over the 20-year simulation period in the 3d-TR (black solid line), 1d-TR (pink dotted line), and 1d-ST experiments (sky-blue solid line).

Appendix B: 1d-TR experiment

While the conventional simplification applied in 1d-ST consists of the two approximations (one-dimensional and steady-state propagation), the impact of oblique wave propagation examined in Sect. 3 should also be confirmed by exclusively applying the one-dimensional simplification but with transient GW parameterization. An additional experiment performed with this simplification (1d-TR, Fig. B1) shows qualitatively similar results to 1d-ST, also exhibiting too long periods of the oscillation (36–48 months) with slow downward penetration of the easterly phase, together with the excessive easterly bias at the 10–20° latitudes of the summer hemisphere.

Figure B1 Time series of vertical profiles (a) and 24 km latitudinal profiles (b) of the tropical stratospheric zonal winds in the experiment using the transient gravity-wave parameterization simplified by the one-dimensional approximation (1d-TR).