Propagation of gravity waves and its effects on pseudomomentum flux in a sudden stratospheric warming event

Effects of realistic propagation of gravity waves (GWs) on distribution of GW pseudomomentum fluxes are explored using a global ray-tracing model for the 2009 sudden stratospheric warming (SSW) event. Four-dimensional (4D; x–z and t) and two-dimensional (2D; z and t) results are compared for various parameterized pseudomomentum fluxes. In ray-tracing equations, refraction due to horizontal wind shear and curvature effects are found important and comparable to one another in magnitude. In the 4D, westward pseudomomentum fluxes are enhanced in the upper troposphere and northern stratosphere due to refraction and curvature effects around fluctuating jet flows. In the northern polar upper mesosphere and lower thermosphere, eastward pseudomomentum fluxes are increased in the 4D. GWs are found to propagate more to the upper atmosphere in the 4D, since horizontal propagation and change in wave numbers due to refraction and curvature effects can make it more possible that GWs elude critical level filtering and saturation in the lower atmosphere. GW focusing effects occur around jet cores, and ray-tube effects appear where the polar stratospheric jets vary substantially in space and time. Enhancement of the structure of zonal wave number 2 in pseudomomentum fluxes in the middle stratosphere begins from the early stage of the SSW evolution. An increase in pseudomomentum fluxes in the upper atmosphere is present even after the onset in the 4D. Significantly enhanced pseudomomentum fluxes, when the polar vortex is disturbed, are related to GWs with small intrinsic group velocity (wave capture), and they would change nonlocally nearby large-scale vortex structures without substantially changing local mean flows.

1 Introduction

Atmospheric gravity waves (GWs) play an important role in the momentum and energy budgets of global circulations in the middle and upper atmosphere. GW pseudomomentum fluxes can induce large-scale momentum forcing, which can substantially change ambient winds, either when transience related to unsteady propagation or dissipation due to breaking or damping occurs (e.g., Fritts and Alexander, 2003; Bühler, 2014).

GWs may also affect global thermal structure through adiabatic vertical motions and heat deposition. GW momentum forcing induces the meridional and vertical mass circulations that contribute to the temperature structure related to the Brewer–Dobson circulation in the stratosphere (e.g., Rosenlof and Holton, 1993; Chun et al., 2011) and can reverse the pole-to-pole radiatively driven latitudinal temperature gradient in the upper mesosphere (e.g., Kim et al., 2003; Smith et al., 2012). Irreversible heat is released when GW momentum forcing is induced (e.g., Becker and Schmitz, 2002; Medvedev and Klaassen, 2003), and it contributes to GW-induced heat deposition.

In general, excitation of GWs is unsteady, and GWs propagate at finite group velocities in the form of localized packets or wave trains. Hence, studies on propagation of GW packets in slowly varying large-scale flows have been carried out using ray-tracing modeling based on the spatial ray theory (e.g., Dunkerton, 1984; Marks and Eckermann, 1995). Hasha et al. (2008) extended the ray theory to spherical geometry. Ribstein et al. (2015) presented more complete formulations in which the magnitude of the three-dimensional (3D) wave number vector is invariant with respect to the Earth's curvature under the deep atmosphere approximation.

However, GW parameterizations (GWPs) for global climate and numerical prediction models have dealt with propagation of GWs under simplifying assumptions that steady GWs propagate instantaneously only in the vertical direction from tropospheric sources to the model top. To consider horizontal and time propagation of GWs, Song and Chun (2008) developed a ray-based GWP for convective GWs for use in the Whole Atmosphere Community Climate Model (WACCM). Senf and Achatz (2011) discussed the validity of the simplifying assumptions of conventional GWPs by computing GW-induced forcing in the spatiotemporally varying large-scale flow associated with the thermal tides in the mesosphere and lower thermosphere (MLT) based on a ray-based method. Kalisch et al. (2014) showed the importance of momentum forcing due to poleward-propagating GWs using a ray-tracing model, and they discussed the implementation of the effects of poleward propagation in global models. Amemiya and Sato (2016) presented a quasi-columnar way of implementing a ray-based GWP in global models, ignoring time propagation of GWs. Yet, it is not clear how ray-based GWPs can be formulated in a way that is consistent with theories on interactions between GW packets and slowly varying mean flows. Moreover, implementation of ray GWPs in models is not straightforward since it requires overcoming the limitations of conventional modeling frameworks where all parameterizations are columnar, and subgrid-scale processes across time steps (e.g., time-propagating GWs) are ignored.

There have been studies to understand the effects of horizontal and transient propagation of GWs on interactions between GWs and slowly varying mean flows. Bühler and McIntyre (2003) presented a theory on wave–mean interaction associated with horizontal refraction of GWs. Bühler and McIntyre (2005) demonstrated a new type of interaction (wave capture) between GWs and horizontally varying vortices using conservation law for the sum of GW pseudomomentum and impulse for GW packets in slowly varying mean flows. Eckermann et al. (2015) showed that horizontal spreading of GWs can be as important as refraction of vertical wave numbers in the wave–mean interaction for orographic GWs. Dunkerton (1981) demonstrated that transient waves with finite vertical group velocities can induce spontaneous mean flow responses such as descent of mean shear layers. Fritts and Dunkerton (1984) and Fritts et al. (2015) explored the roles of self-acceleration of GW phase speeds in wave-induced instabilities and momentum deposition. Muraschko et al. (2015) presented a method based on the phase space Wentzel–Kramers–Brillouin (WKB) theory to accurately compute time–height evolution of wave activity in the time-dependent background flow. Kruse and Smith (2018) confirmed that nondissipative wave–mean interaction due to transience of orographic GWs is nonnegligible over dissipative interactions at initial stages of the generation of orographic GWs. However, despite these various individual efforts, further research is needed to properly consider the effects of horizontal and time propagation of GWs on the interactions between GWs and slowly varying mean flows in representing subgrid-scale GW processes in global models (Plougonven et al., 2020).

Planetary-scale flows in the middle atmosphere, through which GWs propagate, can exhibit substantial spatial inhomogeneity and transience. Substantially disturbed large-scale flows are often found during sudden stratospheric warming (SSW) events in association with large planetary wave (PW) activity (e.g., Albers and Birner, 2014; Song and Chun, 2016), and they may result in substantial changes in horizontal wave numbers and frequencies of propagating GWs. This change in GW spectral properties results in spatiotemporal variations of GW pseudomomentum fluxes.

There have been various modeling studies on the roles of GWs in SSWs. Limpasuvan et al. (2012) demonstrated, using WACCM, that GW momentum forcing is involved in both SSW initiation and recovery from SSWs. Albers and Birner (2014) discussed the roles of GWs in PW resonance before the onset of the 2009 SSW. In recovery phases of SSWs, modeling studies (e.g., Chandran et al., 2013; Limpasuvan et al., 2016) have also reported that combined effects of PWs and GWs are important in the generation and evolution of elevated stratopauses (ESs). However, given that GW refraction and transient propagation cannot be considered in these models with conventional columnar GWPs, there may be limitations in the model-based assessment that is of relative importance between PWs and GWs in evolutions of SSWs and ESs.

Satellite observations have presented evidence of substantial variations of GW activity around SSW onset dates. GW activity is often found to be enhanced in the upper stratosphere before SSW onsets and in high-altitude regions where ESs form in the recovery phase of SSWs (e.g., Yamashita et al., 2013; Thurairajah et al., 2014). These variations of GW activity are also supported by GW-resolving model results for the 2009 SSW (e.g., Yamashita et al., 2010; Limpasuvan et al., 2011). De Wit et al. (2014) demonstrated substantial change in GW momentum fluxes and forcing in the upper mesosphere around the onset of the 2013 SSW using meteor radar observations over Trondheim, Norway. They showed that the magnitude and evolution of estimated GW momentum forcing are comparable to results from WACCM. However, it is unclear how the two estimates of momentum forcing can be similar, even though the modeled upper-mesospheric winds look quite different from the radar observations. This inconsistency may possibly be attributed to the long-distance horizontal propagation of GWs between the lower atmosphere and the upper mesosphere (e.g., Sato et al., 2009; de Wit et al., 2014; Thurairajah et al., 2017).

The present study explores effects of the 4D (x–z, t) propagation of GWs on distributions of pseudomomentum fluxes, a central quantity in GW–mean-flow interaction, for the 2009 SSW. A ray-tracing model for inertia–gravity waves (IGWs) on a sphere, whose prototype was used by Song et al. (2017), is employed to compute trajectories and pseudomomentum fluxes of GWs for specified (time-varying) large-scale flows. Diagnosis of mean flow responses to change in GW pseudomomentum fluxes is not attempted in this study, since slowly varying mean flows are not only modified by GW pseudomomentum but also by the second-order mean pressure fields that can induce mean motions in regions far from localized GW packets (Bühler, 2014). For statistical robustness, ensemble simulations are carried out similar to experiments for stochastic GWPs (e.g., Dunkerton, 1982; Eckermann, 2011). In each simulation, properties of a monochromatic GW packet at a horizontal grid point are randomly drawn from populations for properties of orographic or nonorographic GWs used in GWPs.

The paper is organized as follows. Section 2 presents formulations of the ray-tracing model. Section 3 describes specification of large-scale flow from the ground to the lower thermosphere. Ensembles for parameterized orographic and nonorographic GWs are presented in Sect. 4. In Sect. 5, ray-tracing simulation results for the 2009 SSW are demonstrated by comparing the 4D (x–z, t) and 2D (z, t) results. A summary and discussion is given in the last section.

2 Ray-tracing model for IGWs on sphere

2.1 Kinematic wave theory

A wave packet is defined by a group of phase surfaces over the distance of the order of a dominant wavelength. A ray is a curve for which tangents coincide with a sequence of wave propagation directions (Landau and Lifshitz, 1975).

Kinematic wave theory (Hayes, 1970) relates the ground-based (observed) frequency ω and the 3D wave number k to a variable ψ(r,t), called the phase, as follows:

$$\begin{array}{}\text{(1)}& \mathit{\omega }(\mathit{r},t)=-\partial \mathit{\psi }(\mathit{r},t)/\partial t,\end{array}$$

and

$$\begin{array}{}\text{(2)}& \mathit{k}(\mathit{r},t)={\mathbf{\nabla }}_r\mathit{\psi }(\mathit{r},t),\end{array}$$

where ${\mathrm{\nabla }}_r={\mathit{e}}_{\mathit{\lambda }}/$ $\left(\mathit{r}\cos \mathit{\varphi }\right)\partial /\partial \mathit{\lambda }+$ $({\mathit{e}}_{\mathit{\varphi }}/\mathit{r})\partial /\partial \mathit{\varphi }+$ ${\mathit{e}}_{\mathrm{r}}\partial /\partial \mathit{r}$ in the spherical coordinate system; e_λ, e_ϕ, and e_r are orthogonal unit vectors in the eastward, northward, and radial directions, respectively; λ, ϕ, and r are the longitude, latitude, and radial distance, respectively; r is a position vector; t is time; k can be written as $k{\mathit{e}}_{\mathit{\lambda }}+l{\mathit{e}}_{\mathit{\varphi }}+m{\mathit{e}}_{\mathrm{r}}$; and k, l, and m are zonal, meridional, and vertical wave number components, respectively.

At each (r and t), ω is related to k through a dispersion function (Ω; Bretherton and Garrett, 1968) given by the following:

$$\begin{array}{}\text{(3)}& \mathit{\omega }=\mathrm{\Omega }(k,{\mathrm{\Lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathrm{\Lambda }}_N),\end{array}$$

where Λ_n (n=1, …, N) denotes the properties of the wave propagation medium that vary slowly with respect to phase ψ(r and t).

2.2 Ray-tracing equations

Time evolutions of position, wave number, and observed frequency of a wave packet are described as follows:

$$\begin{array}{}\text{(4)}& \left(\mathrm{d}\mathit{r}/\mathrm{d}t,\mathrm{d}\mathit{k}/\mathrm{d}t\right)=\left({\mathrm{\nabla }}_k\mathrm{\Omega },-{\mathrm{\nabla }}_r\mathrm{\Omega }\right),\end{array}$$

and

$$\begin{array}{}\text{(5)}& \mathrm{d}\mathit{\omega }/\mathrm{d}t=\left(\partial \mathrm{\Omega }/\partial {\mathrm{\Lambda }}_n\right)\partial {\mathrm{\Lambda }}_n/\partial t.\end{array}$$

Here, d∕dt is the time rate of change following the group velocity (c_g) of a wave packet; ∇_k and ∇_r are the partial derivatives with respect to wave numbers and spatial coordinates, respectively; ${\mathrm{\nabla }}_k\mathrm{\Omega }=\partial \mathrm{\Omega }/\partial \mathit{k}=\partial \mathrm{\Omega }/\partial k_i{\mathit{e}}_i=c_{\mathrm{g}i}{\mathit{e}}_i=c_{\mathrm{g}}$, where i (=1, 2, or 3) is the summation index that denotes the zonal, meridional, or radial component in order ${\mathrm{\nabla }}_r\mathrm{\Omega }=\left(\partial \mathrm{\Omega }/\partial {\mathrm{\Lambda }}_n\right){\mathrm{\nabla }}_r{\mathrm{\Lambda }}_n$.

Equation (4) is isomorphic to the Hamilton equation for a physical system characterized by a Hamiltonian denoted by Ω. In deriving the k equation (Eq. 4) from the local time derivative of Eq. (2), a term c_gi(∇_re_i)⋅k appears in which e_i originates from $k_i={\mathit{e}}_i\cdot \mathit{k}={\mathit{e}}_i\cdot (k{\mathit{e}}_{\mathit{\lambda }}+l{\mathit{e}}_{\mathit{\varphi }}+m{\mathit{e}}_{\mathrm{r}})$, but it must become 0 for the form of ray-tracing equations to be independent of choice of coordinate systems (Ribstein et al., 2015). This constraint gives the k equation shown in Eq. (4).

In the computation of Eqs. (4) and (5), component forms are used (see Appendix A1), and the shallow atmosphere approximation ($r=a+z\approx a$, where a is the mean radius of the Earth, and z (≪a) is the height) is applied (Phillips, 1966; Senf and Achatz, 2011). Under this approximation, ∇_r ≈ ∇ = ${\mathit{e}}_{\mathit{\lambda }}/\left(a\cos \mathit{\varphi }\right)\partial /\partial \mathit{\lambda }+({\mathit{e}}_{\mathit{\varphi }}/a)\partial /\partial \mathit{\varphi }+{\mathit{e}}_z\partial /\partial z$, and the magnitude of horizontal wave number $\left|{\mathit{k}}_{\mathrm{h}}\right|$ [$=(k^{\mathrm{2}}+l^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$] is invariant with respect to the Earth's curvature.

2.3 Dispersion relation

The dispersion function Ω is required to compute the ray-tracing equations and is given by the wave dispersion relation.

In the model, the anelastic dispersion relation for IGWs (Marks and Eckermann, 1995) is employed as follows:

$$\begin{array}{}\text{(6)}& {\widehat{\mathit{\omega }}}^{\mathrm{2}}={\left(\mathit{\omega }-k\cdot \mathit{U}\right)}^{\mathrm{2}}={\displaystyle \frac{N^{\mathrm{2}}\left(k^{\mathrm{2}}+l^{\mathrm{2}}\right)+f^{\mathrm{2}}\left(m^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}\right)}{k^{\mathrm{2}}+l^{\mathrm{2}}+m^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}}},\end{array}$$

where $\widehat{\mathit{\omega }}$ (> 0) is the intrinsic frequency; U is the wind vector given by (U, V, and 0); U and V are the zonal and meridional wind components, respectively; N is the static stability; f is the Coriolis parameter; $\mathit{\alpha }=\mathrm{1}/\left(\mathrm{2}H\right)$; and H is the large-scale density ($\stackrel{\mathrm{‾}}{\mathit{\rho }}$) scale height given by $[-(\mathrm{1}/\stackrel{\mathrm{‾}}{\mathit{\rho }})\partial \stackrel{\mathrm{‾}}{\mathit{\rho }}/\partial z{]}^{-\mathrm{1}}$.

The large-scale flow variables (U, V, N², and α²) and f² correspond to Λ_n in Eq. (3), and they are assumed to vary slowly in space and time with respect to GW phases.

2.4 Amplitude equation

Time evolution of wave amplitude in slowly varying large-scale flows is described by the wave action conservation law (Bretherton and Garrett, 1968) as follows:

$$\begin{array}{}\text{(7)}& \partial A/\partial t+\mathrm{\nabla }\cdot \left(c_{\mathrm{g}}A\right)=-A/{\mathit{\tau }}_{\mathrm{dis}},\end{array}$$

where A is the GW action density [$=E/\widehat{\mathit{\omega }}$ (>0)]; E is the phase-averaged GW energy per unit volume; τ_dis (>0) is the wave dissipation timescale (see Appendix A2 for details).

For computation of wave amplitude along a ray, the conservation law (Eq. 7) is changed into an equation for vertical action flux F_A (=c_gzA) after multiplying Eq. (7) by c_gz as follows:

$$\begin{array}{}\text{(8)}& \mathrm{d}{F}_A/\mathrm{d}t-F_A/{\mathit{\tau }}_{\mathrm{def}}=-F_A/{\mathit{\tau }}_{\mathrm{dis}}.\end{array}$$

Here, τ_def is the wave packet deformation timescale (Marks and Eckermann, 1995) at which $\left|F_A\right|$ can increase (decrease) for τ_def>0 (τ_def<0) and is given by the following:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathit{\tau }}_{\mathrm{def}}& =c_{\mathrm{gz}}/\left[\partial c_{\mathrm{gz}}/\partial t+\left(c_{\mathrm{g}\mathit{\lambda }}\partial c_{\mathrm{gz}}/\partial \mathit{\lambda }-c_{\mathrm{gz}}\partial c_{\mathrm{g}\mathit{\lambda }}/\partial \mathit{\lambda }\right)/h_{\mathit{\lambda }}\right.\\ & \left.+\left\{c_{\mathrm{g}\mathit{\varphi }}\cos \mathit{\varphi }\partial c_{\mathrm{gz}}/\partial \mathit{\varphi }-c_{\mathrm{gz}}\partial \left(c_{\mathrm{g}\mathit{\varphi }}\cos \mathit{\varphi }\right)/\partial \mathit{\varphi }\right\}/h_{\mathit{\lambda }}\right],\end{array}\end{array}$$

where c_gλ, c_gϕ, and c_gz are the zonal, meridional, and vertical components of group velocity, respectively; h_λ=acos ϕ.

In the wave–mean interaction, the vertical flux of IGW horizontal pseudomomentum (${\mathit{F}}_{\mathrm{p}}=c_{\mathrm{gz}}{\mathit{p}}_{\mathrm{h}}=c_{\mathrm{gz}}{\mathit{k}}_{\mathrm{h}}A$, wherep_h is the pseudomomentum), rather than the action flux, is a central quantity (Bühler, 2014). Time evolution of F_p along a ray can be obtained by combining results of k_h in Eq. (4) and F_A in Eq. (8). In general, the magnitude and direction of F_p are changed by refraction due to a horizontally varying medium, but $\left|{\mathit{F}}_{\mathrm{p}}\right|$ does not vary owing to curvature terms as $\left|{\mathit{k}}_{\mathrm{h}}\right|$ is invariant with respect to the curvature (Sect. 2.2).

The action conservation equation (Eq. 7) is related to conservation of the angular pseudomomentum of IGWs. Combining the component form equation for k (Eq. A13) multiplied by acos ϕ and the nondissipative form of Eq. (7) gives the angular pseudomomentum conservation law as follows:

$$\begin{array}{}\text{(10)}& \mathrm{d}\mathcal{P}/\mathrm{d}t+\mathcal{P}\mathrm{\nabla }\cdot c_{\mathrm{g}}=\partial \mathcal{P}/\partial t+\mathrm{\nabla }\cdot \left(c_{\mathrm{g}}\mathcal{P}\right)=\mathrm{0},\end{array}$$

where 𝒫 (=kAacos ϕ) is the angular pseudomomentum.

2.5 Dissipation mechanisms

For dissipation of GWs, two separate processes are employed, namely nonlinear wave saturation and molecular diffusion.

Nonlinear saturation is computed by forcing $\left|{\mathit{F}}_{\mathrm{p}}\right|$ not to exceed values for saturated GWs in GW-induced turbulence fields (Lindzen, 1981). The saturation flux (F_p,sat) formulated under the midfrequency approximation (Senf and Achatz, 2011, hereafter SA11) is employed and can be written as follows:

$$\begin{array}{}\text{(11)}& F_{\mathrm{p},\mathrm{sat}}={\mathrm{Fr}}_{\mathrm{c}}^{\mathrm{2}}(\stackrel{\mathrm{‾}}{\mathit{\rho }}/\mathrm{2})\left(\right|k_{\mathrm{h}}|/N)|U_{\mathrm{h}}-c_{\mathrm{p}}{|}^{\mathrm{3}},\end{array}$$

where Fr_c is the critical Froude number (McFarlane, 1987, M87 hereafter) and is set equal to $\mathrm{1}/\sqrt{\mathrm{2}}$ (Hasha et al., 2008), U_h is the horizontal wind parallel to k_h ($U_{\mathrm{h}}=\mathit{U}\cdot {\mathit{k}}_{\mathrm{h}}/\left|{\mathit{k}}_{\mathrm{h}}\right|$), and c_p is the ground-based phase speed ($c_{\mathrm{p}}=\mathit{\omega }/\left|{\mathit{k}}_{\mathrm{h}}\right|)$.

Molecular diffusion is important above the upper mesosphere. Dividing the total GW energy equation by $\widehat{\mathit{\omega }}$ gives the term τ_dis due to viscous damping (Eq. A24). In the model, kinematic viscosity is set equal to thermal diffusivity (i.e., Pr=1, where Pr is the viscous Prandtl number), and thus a complete form of τ_dis Eq. (A26) is simplified as follows:

$$\begin{array}{}\text{(12)}& {\mathit{\tau }}_{\mathrm{dis}}=\mathrm{1}/\left[\mathrm{2}\mathit{\nu }\left(k^{\mathrm{2}}+l^{\mathrm{2}}+m^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}\right)\right],\end{array}$$

where ν is the kinematic molecular viscosity.

Kinematic viscosity (ν) is defined as $\mathit{\mu }/\stackrel{\mathrm{‾}}{\mathit{\rho }}$, and viscosity μ is determined as $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{5}}$ kg m⁻¹ s⁻¹ considering reported values of ν. Vadas and Fritts (2005) employed ν=6.5 m² s⁻¹ at z=90 km, and Pitteway and Hines (1963) suggested ν=4 m² s⁻¹ at the same height. These two values are roughly consistent with the abovementioned value of μ for possible range of $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ at z=90 km.

2.6 Numerical implementation

Time integrations of the ray-tracing equations (Eqs. 4–5 or Eqs. A10–A16) are carried out using the Livermore solver for ordinary differential equation (ODE) with automatic method switching (LSODA) based on the stiffness of an ODE system (Petzold, 1983; Hindmarsh, 1983).

Solutions (λ, ϕ, z, k, l, m, and ω) of the ray-tracing equations proceed in time over a period of multiples of a time step (δt). The time step (δt) is determined as 900 s through some tests (see Figs. S1 and S2 in the Supplement). Large-scale flow variables (U, V, N², and α²) are given at an interval of Δt, and Δt should be a multiple of δt for proper time marching of the LSODA. In this study, Δt=1 h is used.

The solver requires interpolation of Λ_n at space and time locations of a GW packet. For spatial interpolation, a local C¹-continuous tricubic method (Lekien and Marsden, 2005) is employed. The C¹ continuity allows for accurate computation of the equation for k that involves the first-order spatial derivatives of Λ_n. For temporal interpolation, simple linear interpolation is used since Λ_n is assumed to vary linearly during the time interval (Δt) of large-scale variables.

Action flux equation (Eq. 8) is actually an equation for $\left|F_A\right|$. Note that τ_def and τ_dis required for computation of Eq. (8) do not change the sign of F_A.

The dissipation timescale (τ_dis) can be computed along individual rays using ray solutions and large-scale variables at the positions of rays (see Eq. A26). Meanwhile, computation of the deformation timescale (τ_def) requires ray-tube information related to spatiotemporal variations of c_gλ, c_gϕ, and c_gz in the neighborhood of a ray (see Eq. 9). In the model, τ_def is estimated through a gridding method as follows: c_gλ, c_gϕ, and c_gz of rays are recorded and accumulated at vertices of grid cells (ΔλΔϕΔz) that contain ray paths during δt. Then, τ_def is computed at each grid point using a finite difference form of Eq. (9) and gridded group velocity components averaged over overlapped rays. This gridding method is crude compared to the 2D (z–t) phase–space theory (Muraschko et al., 2015), but it is used to estimate, even roughly, ray-tube effects in the 4D (r-t) space.

After τ_def and τ_dis are obtained, Eq. (8) is computed using the Euler method for the time step δt. In case that vertical propagation direction is reversed after δt, the sign of F_A is changed considering the sign of m because $\mathrm{sgn}\left(F_A\right)=\mathrm{sgn}\left(c_{\mathrm{gz}}\right)=-\mathrm{sgn}\left(m\right)$. Further details of numerical implementation are described in Appendix A3.

3 Large-scale atmospheric flow

Time-varying large-scale flows during the 2009 SSW are specified by combining 6 h reanalysis data sets and empirical model results. The reanalysis data are linearly interpolated at a 1 h interval (Δt=1 h). The empirical model results are obtained at the hourly interval using daily 10.7 cm solar flux (F10.7) and 3 h geomagnetic activity (Ap) indices. Hourly whole atmospheric flows for ray modeling are obtained by fitting the third-order B spline curves in the vertical to the time-interpolated hourly reanalysis data and empirical model results in four overlapping layers. Details of the B spline fit can be found in Song et al. (2018).

The four vertical layers are (i) p=10³–1 hPa (z≈0–48 km), (ii) 400–0.1 hPa (7–65 km), (iii) 1–$\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ hPa (48–94 km), and (iv) 5×10⁻³–10⁻⁸ hPa (84–331 km), respectively. Data used in the four layers are as follows (in order of layer altitudes): (i) European Centre for Medium-Range Weather Forecast (ECMWF) Interim (ERA-Interim; Dee et al., 2011) reanalysis, (ii) Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2; Gelaro et al., 2017) reanalysis, (iii) the Advanced-Level Physics High-Altitude (ALPHA) prototype of the Navy Operational Global Atmospheric Prediction System (NOGAPS, collectively NOGAPS–ALPHA; Eckermann et al., 2009) data, and (iv) empirical model results for temperature (NRLMSISE-00; Picone et al., 2002) and geomagnetically quiet time horizontal winds (HWM14; Drob et al., 2015) and disturbed horizontal winds (DWM07; Emmert et al., 2008).

Figure 1 shows latitude–height cross sections of ground-to-space (G2S) zonal wind and temperature at 60^∘ W at 00:00 UTC on 23 January 2009, 1 d before the central date (24 January) of the 2009 SSW. The G2S data demonstrate that vertically smooth whole atmospheric wind and temperature profiles can be constructed by fitting the B spline curves in the vertical direction to the reanalysis data and empirical model results. In the Northern Hemisphere (NH), the polar night jet is already reversed above the lower stratosphere north of 60^∘ N, and the weakened eastward jet is tilted from the midlatitudes towards the Equator in the lower mesosphere. In association with the jet structure, substantial warming is found in the northern (winter) polar stratosphere, and temperature maximum (≈280 K) is as high as that in the summer polar stratopause. In the Southern Hemisphere (SH), typical summertime wind and temperature structure is found, namely easterly flow in the middle atmosphere below the upper mesosphere, warm temperature near the polar stratopause region, and coldest temperature and wind reversals in the polar upper mesosphere.

Figure 1 Latitude–height cross sections of (a) zonal wind and (b) temperature in the ERA-Interim, MERRA-2, NOGAPS–ALPHA, empirical models, and G2S data at 60^∘ W at 00:00 UTC on 23 January 2009. For zonal wind, shading and contour intervals are 2 and 20 m s⁻¹, respectively. Contours for westward winds are plotted in broken lines. For temperature, shading and contour intervals are 5 and 20 K, respectively.

For ray simulations, G2S data at hourly intervals are spatially smoothed using the vertical 1–2–1 smoother and horizontal moving averaging on the spherical surface. Horizontal averaging is done using variables within the area of a spherical cap centered at every latitude–longitude grid point. A spherical cap is defined by an angle between two lines from the sphere center to the center of the spherical cap's surface and to the cap's boundary. The angle is set equal to about 2.7^∘, and for this angle, the area of the spherical cap's surface is equivalent to the area of a circle with a radius of 300 km [π(300 km)²] on a flat surface. Smoothed G2S data are regridded at $\mathrm{2.5}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ horizontal resolution for use in ray simulations.

4 GW ensembles

In this study, orographic and nonorographic GWs are considered separately. Properties of orographic GWs are given by the M87 scheme. Nonorographic GWs are specified based on SA11 and Warner and McIntyre (1996, hereafter WM96).

4.1 Orographic GWs

Orographic GWs (OGWs) are assumed to be excited when low-level horizontal winds are strong, and vertical parcel displacements, due to subgrid topography, are large (M87).

The vertical displacement (h_m) is given by 2σ_h, where σ_h is the standard deviation of the subgrid-scale topography. Reynolds stress due to stationary OGWs (ω=0 and c_p=0) is given by (k_ho∕2) $\min (h_m^{\mathrm{2}},U_{\mathrm{s}}^{\mathrm{2}}/N_{\mathrm{s}}^{\mathrm{2}})$ ρ_sN_sU_s, where k_ho is the horizontal wave number [$=\mathrm{2}\mathit{\pi }/\left(\mathrm{100}\mathrm{km}\right)$], and ρ_s, N_s, and U_s are air density, stability, and horizontal wind magnitude, respectively, averaged within the source layer between the ground and the h_m. Grid- and subgrid-scale topography are obtained by averaging and regridding the National Center for Atmospheric Research (NCAR) Community Earth System Model (CESM) auxiliary data with horizontal resolutions close to $\mathrm{2.5}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$.

Directions of horizontal wave number vectors are set opposite to horizontal wind vectors averaged within source layer. OGWs are launched at the top of the source layers.

4.2 Nonorographic GWs

For nonorographic GWs (NOGWs), three 14-discrete-wave schemes, as in SA11, are considered. One is a modified version of SA11, and the other two are derived from the empirical spectra of WM96.

The empirical GW energy spectrum ($\widehat{E}$) in WM96 is given by a separable function of m, $\widehat{\mathit{\omega }}$, and the azimuth angle (φ; i.e., $\widehat{E}=E_{\mathrm{0}}A\left(m\right)B\left(\widehat{\mathit{\omega }}\right)\mathrm{\Phi }\left(\mathit{\phi }\right)$). The pseudomomentum flux spectrum $\stackrel{\mathrm{‾}}{\mathit{\rho }}{c}_{\mathrm{gz}}\widehat{E}(m,\widehat{\mathit{\omega }},\mathit{\phi })\mathit{k}/\widehat{\mathit{\omega }}$ can be written as a function of k, ω, and φ by multiplying the spectrum by the Jacobian factor ($J=m/\left|\mathit{k}\right|$). This spectrum is discretized to obtain 14-wave schemes through numerical integrations around appropriately chosen k and ω for 14 sets of φ and c_p, as in SA11, in a quiescent atmosphere with a specified stability (0.01 rad s⁻¹). Two 14-wave WM96 schemes are obtained by using two different values (1 and 5∕3) of p in the spectrum $B\left(\widehat{\mathit{\omega }}\right)$ given by $B_{\mathrm{0}}\left(p\right){\widehat{\mathit{\omega }}}^{-p}$. NOGWs are launched at every horizontal grid point at z=6.8 km near 400 hPa.

Figure 2 illustrates angular histograms of spectral properties and Reynolds stress in the three 14-wave NOGW schemes. In these schemes, horizontal propagation directions (φs) and ground-based phase speeds (c_ps) are given for each of 14 GWs, and they are identical to those in SA11. The horizontal wavelength (λ_h) in SA11 ranges from 385 to 596 km. In WM96 with $p=\mathrm{5}/\mathrm{3}$ (WM96a), the range of λ_h is broader (309–782 km) compared to SA11, and in WM96, with p=1 (WM96b), the range is much broader (128–942 km).

Each GW has an identical amount of Reynolds stress in the three schemes. For this, the stresses for GWs with c_p>20 m s⁻¹ are reduced in SA11, and integration intervals of k and ω for the spectra in WM96 are appropriately adjusted. As a result, NOGWs with c_ps smaller (larger) than 20 m s⁻¹ have Reynolds stress to the order of 10⁻³ (10⁻⁵) N m⁻². Reynolds stresses exhibit slightly more westward flux, but they are almost isotropic. The total magnitudes of Reynolds stresses used in this study are $\mathrm{3.6}\times {\mathrm{10}}^{-\mathrm{3}}$ N m⁻² in the eastward, northward, and southward directions and $\mathrm{4.1}\times {\mathrm{10}}^{-\mathrm{3}}$ N m⁻² in the westward direction. These magnitudes are comparable to the total momentum flux of $\mathrm{3.75}\times {\mathrm{10}}^{-\mathrm{3}}$ N m⁻² in each cardinal direction at 450 hPa that is employed in the ECMWF model (Orr et al., 2010). Details of GW properties shown in Fig. 2 are found in Table S1 in the Supplement.

Figure 2 Angular histograms of (a) phase speeds, (b) Reynolds stresses, and (c)–(e) horizonal wavelengths of nonorographic GWs (SA11, WM96a, and WM96b) as a function of propagation directions (φ) at an interval of 45^∘. For wave IDs of 1–8 (9–14), c_p=6.8, 6.8, 10.2, 6.8, 6.8, 6.8, 10.2, and 6.8 m s⁻¹ (32.8, 20.4, 20.4, 32.8, 20.4, and 20.4 m s⁻¹) counterclockwise from due east (φ=0^∘).

4.3 Generation of GW ensembles

In ray simulations, a single GW packet, stochastically chosen from GW source ensembles, is launched at a horizontal grid point where GWs are supposed to be generated. This approach is computationally efficient, allowing for statistical significance tests for differences between ray simulations.

The OGW scheme (M87) launches a single OGW at a horizontal grid point, but NOGW schemes usually specify multiple GWs. Hence, GW ensembles are separately generated for OGWs and NOGWs. For OGWs, the vertical displacement h_m is assumed to be given by s_fσ_h, where the scale factor s_f has a uniform probability distribution between 1 and 3 around its default value 2. For NOGWs, a single GW is selected with uniform chance from 14 discrete waves. GW ensembles are precomputed using a random-number generator for reproducibility.

Figure 3 demonstrates horizontal distributions of zonal OGW and NOGW pseudomomentum fluxes (F_p) at individual launch levels for a particular ensemble member and the ensemble averages of the zonal pseudomomentum fluxes at 00:00 UTC on 20 January 2009, 4 d before the 2009 SSW onset. The stochastic parameters for OGWs and NOGWs are the displacement scale factors (s_f) and wave IDs (1–14), respectively. For each of GW source schemes, 20 ensemble members are generated.

For OGWs, s_f varies randomly between 1 and 3 at grid points where σ_h is nonzero. In the midlatitude tropospheric eastward jet regions, westward OGW F_p is large and may reach about −1 N m⁻² over major mountainous regions, namely the Alps, the Tibetan Plateau, the Rocky Mountains, and the Andes. Large eastward OGW F_p is found in the higher-latitude regions such as Greenland and Antarctica. Zonal F_p in each OGW ensemble member has locally substantial deviations from the ensemble mean (Fig. 3b) in the major mountain areas. Maximum value of the standard deviation from the ensemble mean is about 0.7 N m⁻² and is found in the Tibetan areas.

Figure 3 Longitude–latitude distributions of zonal pseudomomentum fluxes (F_p) for OGW and NOGW_SA11 at 00:00 UTC on 20 January 2009 on the $\mathrm{2.5}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ horizontal grid as follows: (a) zonal OGW F_p above source layers for the first OGW ensemble member, (b) ensemble-averaged zonal OGW F_p, (c) zonal NOGW F_p at z=6.8 km for the first NOGW ensemble member, and (d) ensemble-averaged zonal NOGW F_p at z=6.8 km. OGW F_p is multiplied by an efficiency factor (0.125) as described in Richter et al. (2010).

For NOGWs, wave IDs of 1–14 are randomly spread around the globe. The magnitude of ensemble-averaged zonal NOGW F_p is O(10⁻³ N m⁻²). The horizontal distribution of zonal F_p in each NOGW ensemble member (Fig. 3c) looks noisy, but characteristic structure is more clearly revealed in its ensemble average (Fig. 3d). The sign of zonal NOGW F_p is generally opposite to that of the tropospheric zonal flows. The organized spatial structure of the zonal NOGW F_p only emerges in the ensemble because any single ensemble member is completely stochastic. Meanwhile, the OGW F_p is well organized in a single ensemble member because they are constrained by the OGW source parameterization that depends explicitly on orography.

5 Results

GW ray simulations are carried out for the time period of 25 d from 00:00 UTC on 8 January 2009 to 00:00 UTC on 2 February 2009 for the 20 OGW and NOGW ensemble members.

For each GW ensemble member, two kinds of simulations are carried out, namely four-dimensional (4D; x–z and t) and two-dimensional (2D; z and t) experiments. In the 4D experiments, GW rays propagate horizontally and vertically in spatially varying background media, but in the 2D, rays are allowed to propagate only in the vertical direction. In both the 4D and 2D, GW rays propagate through time-varying flows, and therefore modulations of the observed frequencies of GWs occur. In the 4D cases for M87 and SA11, additional simulations where τ_def=0 in the amplitude equation are carried out to see ray-tube effects. In all simulations, 3 h averaged gridded outputs are generated. GW rays are launched every 3 h, and 3 d old rays are eliminated. These launch intervals and ray lifetimes are chosen considering computational time, the timescale of the large-scale flow, and elapsed time for rays launched in the troposphere to reach the upper mesosphere and lower thermosphere (UMLT; see Figs. S1–S3).

5.1 Zonally averaged GW properties

Figure 4 shows latitude–height distributions of zonal mean zonal wind and ensemble averages of zonal mean zonal F_p of OGWs and three NOGWs in the 4D and 2D experiments at 00:00 UTC on 20 January 2009. The NH polar vortex is not reversed but is much weakened. Hatched areas indicate regions where differences between the 4D and 2D are not statistically significant. The signs of the zonal mean zonal pseudomomentum fluxes below z=80 km are similar overall in the 4D and 2D and seem to be related to the structure of the zonal wind, but they become different from each other above z=80 km in the NH polar regions.

Statistically significant differences between the 4D and 2D are found in five regions. (i) In the latitude–height region of 40–70^∘ N and 40–50 km, westward F_p in the 4D is enhanced about 10 times when compared with the 2D results for OGW, NOGW_SA11, and NOGW_WM96a, and they are about 28 times increased when compared with the 2D results for NOGW_WM96b. (ii) In the NH polar UMLT (70–80^∘ N and 85–100 km), the magnitude of eastward F_p in the 4D is 1.5–5 times larger than that of westward F_p in the 2D for both OGWs and NOGWs. This result implies that the 4D results may better explain the mesospheric cooling around the SSW central dates, compared to the 2D, given that the cooling can be induced by eastward GW momentum forcing in the UMLT. (iii) In the upper troposphere and lower stratosphere (UTLS) above the tropospheric jets, westward F_p of NOGWs is 1.6–2.6 times larger in the 4D. (iv) In the SH mesosphere, eastward F_p of NOGWs in the midlatitudes is reduced by more than half in the 4D. (v) Eastward F_p of OGWs at z=30–80 km around 70^∘ S is enhanced 5–6 times in the 4D. As a result, the magnitude of zonal F_p in the middle atmosphere and its structure in the UMLT can be substantially changed in the 4D where horizontal propagation and refraction are allowed, even though F_p is given identically at source levels in the 4D and 2D experiments.

Additional 4D OGW and NOGW_SA11 experiments (Fig. S4) where no ray-tube effects are considered (τ_def=0) give similar results to those with τ_def≠ 0 shown in Fig. 4, except for the NH upper stratosphere and lower mesosphere (USLM). Statistically significant differences between nonzero and zero τ_def are found in relatively narrow regions around the NH jet core (75^∘ N and 40 km at 00:00 UTC on 20 January 2009; see Fig. S4). It is interesting that ray-tube effects become important in regions near the jet where the large-scale winds vary rapidly in space (see Sect. 5.2). In these regions, the magnitude of westward F_p when τ_def≠0 is reduced by less than half compared to when τ_def=0. These differences are localized compared to the differences between the 4D and 2D experiments, which may indicate that ray-tube effects are relatively limited. However, it should be noted that ray simulations in this study may underestimate ray-tube effects, given that horizontal spread of GW fields that emanate from local sources cannot be considered in the current ray simulations where a single GW ray is launched at a grid point.

Figure 4 Latitude–height cross sections of (a) zonal mean zonal wind (ZU) and (b–i) zonal mean zonal pseudomomentum fluxes (F_p) for OGW and three NOGW schemes in the (top row) 4D and (bottom row) 2D experiments at 00:00 UTC on 20 January 2009. OGW F_p is multiplied by the efficiency factor (0.125). Contour interval of zonal mean zonal wind is 10 m s⁻¹, and negative values are plotted in broken lines. Hatched areas on the pseudomomentum fluxes indicate regions where the paired and two-tailed t test for 20 ensemble members of the 4D and 2D experiments give p values larger than 0.05 (i.e., no statistical significance at the level of 0.05). Here, the p value means there is a probability that mean values in the 4D and 2D experiments would be similar to each other. The hatched areas are identical in the pair of 4D and 2D experiments for a particular GW scheme. A non-parameteric test, such as Wilcoxon's signed rank test, where no probabilistic distribution is assumed, also gives almost similar results to the t test. For these statistical tests, algorithms presented in Boslaugh (2013) are employed.

Figure 5 shows latitude–height distributions of zonal mean zonal wind and ensemble means of zonally averaged zonal wave numbers (k) for OGWs and NOGWs in the 4D and 2D experiments at 00:00 UTC on 20 January 2009. It is clear that the sign of zonal mean k is overall the same as that of the zonal F_p shown in Fig. 4, except for some regions in the NH USLM. This similarity in the signs of k and zonal F_p implies that the zonal mean zonal F_p (Fig. 4) is mostly due to upward-propagating GWs, since the sign of zonal F_p (=kF_A) is determined by the sign of k alone in case that c_gz> 0 and thus F_A> 0. In the NH USLM, however, GWs in the 4D are found to propagate downward in some areas, and positive k at z=30–50 km is related to the GWs with negative c_gz and westward F_p (see Sect. 5.2 for details).

Figure 5 Latitude–height cross sections of (a) zonal mean zonal wind (ZU) and (b–i) zonal mean zonal wave numbers (k) for OGW and three NOGW schemes in the 4D and 2D experiments at 00:00 UTC on 20 January 2009. Contour interval of zonal mean zonal wind is 10 m s⁻¹ and negative values are plotted in broken lines. Hatched areas on the zonal wave numbers indicate regions where the paired and two-tailed t test for the 4D and 2D experiments gives p values larger than 0.05 (i.e., no statistical significance at the level of 0.05).

Statistically significant differences of zonal mean k between the 4D and 2D are also found in similar regions to those of zonal mean zonal F_p shown in Fig. 4. In the NH polar UMLT, the sign of k is reversed between the 4D and 2D, and the magnitude of positive k in the 4D is 1.2–6.3 times larger than that of negative k in the 2D for both OGWs and NOGWs. In the UTLS, negative k of NOGWs is 1.4–2.4 times larger in the 4D. In the SH mesosphere, positive k of NOGWs in the midlatitude regions is reduced by about half in the 4D, and positive k of OGWs around 70^∘ S is enhanced about 1–3 times in the 4D. These changes in k in the 4D with respect to the 2D are roughly similar to those in the zonal F_p shown in Fig. 4. This result indicates that differences in the zonal F_p between the 4D and 2D experiments can be accounted for to a substantial degree by changes in the k between the 4D and 2D.

The time rate of change of k along rays is determined by the five forcing terms (Eqs. A13 and A17). Among the five terms, the two zonal shear forcing terms [$-k/\left(a\cos \mathit{\varphi }\right)\partial U/\partial \mathit{\lambda }$ and $-l/\left(a\cos \mathit{\varphi }\right)\partial V/\partial \mathit{\lambda }$] and curvature term (lc_gλtan ϕ∕a) are predominant in most cases. The other forcing terms in M_λ∕h_λ (Eq. A17) that depend on the stability (N²) and vertical density variation (α²) are usually 2–3 orders of magnitude smaller (see Fig. S5).

Figure 6 Latitude–height cross sections of total and three major forcing terms – the zonal shear terms of the large-scale zonal and meridional winds (U and V) and the curvature term – of the zonal wave number for (top row) OGW and (bottom row) NOGW_SA11 in the 4D experiment at 00:00 UTC on 20 January 2009. The zonal shear terms of U and V are $-k/\left(a\cos \mathit{\varphi }\right)\partial U/\partial \mathit{\lambda }$ and $-l/\left(a\cos \mathit{\varphi }\right)\partial V/\partial \mathit{\lambda }$, respectively. The curvature term is given by lc_gλtan ϕ∕a. Hatched areas indicate regions where the paired and two-tailed t test for the 4D and 2D experiments gives p values larger than 0.05 (i.e., no statistical significance at the level of 0.05).

Figure 6 shows latitude–height distributions of the total and three major forcing terms in the k equation for OGWs and NOGW_SA11's in the 4D experiment at 00:00 UTC on 20 January 2009. Note that the forcing terms in the k equation in the 2D are all zero. For NOGWs, results for SA11 are presented since the other NOGW schemes give roughly similar results. It is clear that the magnitude of the curvature term is as large as the two zonal shear forcing terms in the midlatitudes as well as the polar regions, which supports the importance of the curvature term presented by Hasha et al. (2008). In the UTLS, the total forcing term for NOGWs is generally negative above the tropospheric jet cores (60–40^∘ S and 20–50^∘ N) where the negative k is predominant and enhanced in the 4D, and the negative total forcing is mainly due to the zonal shear term of zonal winds and curvature effects. For OGWs around 70^∘ S, the enhancement of the positive k in the stratosphere is attributed to curvature effects. For NOGWs, curvature effects are predominant over the other two major forcing terms in the SH stratosphere where winds are steady, but they are a little smaller than the two zonal shear forcing terms in the NH where the polar vortex varies rapidly in space and time.

The structure of the three major forcing terms in the k equation (Fig. 6) is different from that of k (Fig. 5), except for the NH USLM. This difference may occur in relation to space and time propagation of GWs, since certain k's substantially changed in some regions may not be changed a lot as GWs propagate to the other regions. In fact, positive k of NOGWs with eastward F_p in the UTLS of the SH mid- to high-latitude regions is increased by positive zonal shear terms for the zonal wind and curvature terms (Fig. S5). The increase in positive k enhances eastward F_p in the SH UTLS. However, the eastward F_p of NOGWs is reduced in the SH middle atmosphere, even though total forcing in the k equation for NOGWs with eastward F_p is positive (Fig. S5). It is seen that the GWs with enhanced eastward F_p would be dissipated through the saturation as they propagate northward toward the SH midlatitude USLM. This dissipation results in reduction in eastward F_p and positive k in the SH mesosphere (Figs. 4–5). In contrast, distributions of k and major forcing terms (Fig. 6) are correlated with each other in the NH USLM, which indicates that the structure of k in the NH USLM can be locally generated by forcing terms around z=40 km.

Figure 7 Latitude–height cross sections of zonal mean zonal wind (ZU) and ensemble means of zonally averaged meridional wave numbers (l) for OGW and NOGW_SA11 in the 4D and 2D experiments at 00:00 UTC on (top row) 15 and (bottom row) 20 January 2009. Contour interval for zonal mean zonal wind is 10 m s⁻¹, and negative values are plotted in broken lines. Hatched areas on the meridional wave numbers indicate regions where the paired and two-tailed t test for the 4D and 2D experiments gives p values larger than 0.05 (i.e., no statistical significance at the level of 0.05).

Figure 7 shows latitude–height distributions of zonal mean zonal wind and ensemble means of zonally averaged meridional wave numbers (l) for OGWs and NOGW_SA11's in the 4D and 2D experiments at 00:00 UTC on 15 and 20 January 2009. It is found that the structural difference of l between the 4D and 2D is substantial for both OGWs and NOGWs and is more significant than the difference in structure of k. This result may be related to larger meridional variations of the large-scale flow than its zonal variations, given that the time rates of changes of l and k are determined by meridional and zonal variations of the large-scale flow, respectively.

In the SH, large positive l of OGWs appears around 70^∘ S in the 4D compared with the 2D on both of the two dates, and these are due to the positive meridional shear forcing term related to zonal wind [$-(k/a)\partial U/\partial \mathit{\varphi }$], where k>0, and $\partial U/\partial \mathit{\varphi }<\mathrm{0}$ (see Fig. S6). In the NH USLM, positive (negative) l of OGWs in the 4D appears roughly south (north) of the eastward jet axis on 15 January but significantly enhanced positive l is predominant on 20 January. The signs of l across the jet on 15 January mean northward (southward) propagation of OGWs relative to large-scale meridional winds. Therefore, OGWs in the 4D can in fact be converged along the jet axis in the NH USLM, as long as the meridional variations of meridional winds are not significant. On 20 January, as the jet is moved towards 80^∘ N, it is seen in the 4D that positive l of OGWs and poleward propagation of OGWs become predominant south of the jet axis. In the 2D, the structure of l of OGWs in the NH does not seem to exhibit significant responses to spatiotemporal variations of the large-scale flow.

The structure of l for NOGWs is coherent with the large-scale flow structure in the 4D, but it seems more or less random in the 2D. In the SH, l of NOGWs in the 4D is positive overall in the high-latitude stratosphere above z=10 km and in the midlatitude middle atmosphere above z=30 km, and c_gλ is negative overall due to the westward winds. This explains why the curvature term (lc_gλtan ϕ∕a) in the k equation is positive for NOGWs in the SH (ϕ<0) middle atmosphere, as shown in Fig. 6h. The positive l of NOGWs in the SH is mainly due to the forcing term related to the meridional shear of zonal winds and curvature term ($-k{c}_{\mathrm{g}\mathit{\lambda }}\tan \mathit{\varphi }/a$) in the l equation (see Fig. S6). Positive k above z=40 km (Fig. 5) and changes in the sign of $\partial U/\partial \mathit{\varphi }$ around the axis of the westward winds ($\partial U/\partial \mathit{\varphi }<\mathrm{0}$ south of the axis and $\partial U/\partial \mathit{\varphi }>\mathrm{0}$ north of the axis) make the meridional shear forcing term in the l equation positive in the SH polar USLM and negative in the SH midlatitude USLM. In the NH USLM, as in OGWs, the structure of the l on 15 January may indicate latitudinal convergence of NOGW rays along the axis of the stratospheric jet. On 20 January, positive l is predominant in the NH USLM, and NOGWs south of 80^∘ N propagate northward.

Figure 8 Same as Fig. 7 except for the zonally averaged number of GW packets.

Figure 8 shows latitude–height distributions of zonal mean zonal wind and ensemble means of the zonally averaged number of GW packets for OGWs and NOGW_SA11's in the 4D and 2D experiments at 00:00 UTC on 15 and 20 January 2009. The number of GW packets corresponds to the ray counter used for averaging in the gridding method. Comparison between the 4D and 2D indicates that there is convergence of GW packets in the 4D along the jet axis in the NH USLM and in the NH polar stratosphere on 15 January for both OGWs and NOGWs. In the latitude–height region of 40–70^∘ N and 40–90 km on 15 January, the zonally averaged number of rays for OGWs (NOGWs) is about 0.9 (1.9) in the 4D, and it is 1.7 (1.3) times larger than in the 2D. In this region, GWs with westward F_p outnumbers GWs with eastward F_p, and therefore the ray convergence may increase the westward F_p. In the NH, mid- to high-latitude lower thermosphere (30–90^∘ N and 120–140 km), the zonally averaged number of rays for OGWs (NOGWs) is about 2.5 (1.9) times larger in the 4D on 15 January. This difference is mainly due to OGWs and NOGWs with eastward F_p (not shown).

As the eastward jet in the stratosphere moves towards the North Pole on 20 January, the number of GW packets in the 4D increases in the NH polar mid- to upper stratosphere (50–80^∘ N and 30–50 km). In this region, the zonally averaged number of rays for OGWs (NOGWs) is about 0.9 (2.3) in the 4D, and it is 1.7 (1.1) times larger than in the 2D. As on 15 January, GWs with westward F_p are predominant in the NH stratosphere, resulting in some enhancement of westward F_p. In the NH lower thermosphere (30–90^∘ N and 120–140 km), the zonally averaged number of rays for both OGWs and NOGWs is about 1.8 times larger in the 4D. This increase is mostly attributed to OGWs and NOGWs with eastward F_p, as on 15 January.

GWs generally propagate more to the thermosphere in the 4D. Even though the eastward winds are still large in the NH middle atmosphere on 15 and 20 January, both OGWs and NOGWs with eastward F_p (i.e., k> 0 and $c_{\mathrm{p}\mathit{\lambda }}-U>\mathrm{0}$, where c_pλ is the ground-based zonal phase speed) can propagate better to the thermosphere in the 4D since they are less dissipated in the middle atmosphere. There seems to be a tendency in the 4D for GWs to elude critical-level filtering or nonlinear saturation better in the lower atmosphere. This tendency may be attributed to a larger degree of freedom in propagation associated with change in wave number directions due to refraction and curvature effects that can occur before either filtering or saturation is initiated as GWs approach critical levels. Similar results can also be found for NOGWs in the SH. More NOGW packets in the 4D are found in the SH middle atmosphere, and they are related to reduced restriction in the propagation of GWs with eastward F_p towards the middle atmosphere. Also, as GWs propagate better toward the upper atmosphere, GW packets (with eastward F_p) in the 4D look vertically more spread near z=90 km in the SH where the zonal mean zonal wind is reversed. This spread in the 4D compared to the 2D implies that GWs in the 4D may better avoid filtering without being trapped in a narrow vertical layer close to critical levels.

Convergence (or focusing) of GW packets may have some effects on distribution of the GW pseudomomentum fluxes. Song and Chun (2008) emphasized this effect as an important mechanism that accounts for differences between ray-based and columnar GWPs. Discussion regarding Fig. 8 indicates that convergence or divergence (spread of GW packets) effects in the 4D would be smaller than a factor of 2 in terms of the magnitude of F_p. These convergence and divergence effects are, however, likely relatively small compared to impacts due to change in horizontal wave numbers shown in Figs. 4–5, although the refraction effects would not entirely lead to local change in mean flows (see Sect. 5.2).

The spatial distribution of GW packets in the 4D is generally more contiguous in the latitudinal direction than that in the 2D in which latitudinal discontinuity is clear. This difference indicates more than improvement in the smoothness of the GW pseudomomentum fluxes in the 4D. Since McLandress and McFarlane (1993), various modeling studies have suggested that PW momentum forcing tends to be compensated by parameterized GW momentum forcing. Cohen et al. (2013) noted that the GW momentum forcing on small latitudinal scales can induce instability that generates PWs to be involved in the compensation. Given that the latitudinal discontinuity of the pseudomomentum fluxes in the 2D is unrealistic compared to the 4D, the discontinuity due to columnar GWPs can possibly generate spurious PWs in models. Therefore, it cannot be said that compensation between PW and GW forcing always occurs for the right reasons in models with columnar GWPs. For further understanding of the compensation, the 4D formulation beyond columnar GWPs can be required because change in the GW pseudomomentum due to the horizontal refraction can induce change in mixing of the mean potential vorticity (PV) around an aggregate of the refracted GWs. This GW-induced PV mixing can produce changes in PW-induced PV mixing in the surf zone related to PW breaking, resulting in compensation between PW and GW effects (Cohen et al., 2014).

Figure 9 Longitude–latitude cross sections of (a) zonal wind (U), (b) meridional wind (V), and ensemble means of (c–d) zonal pseudomomentum fluxes (F_p), (e–f) zonal wave numbers, and (g–h) vertical wave numbers for OGW at z=38 km in the 4D and 2D experiments at 00:00 UTC on 20 January 2009. OGW F_p is multiplied by the efficiency factor (0.125). Contour interval for zonal and meridional winds is 20 m s⁻¹, and negative values are plotted in broken lines. Hatched areas indicate regions where the paired and two-tailed t test for the 4D and 2D experiments gives p values larger than 0.05 (i.e., no statistical significance at the level of 0.05).

5.2 Horizontal distributions of GW characteristics

Figure 9 shows longitude–latitude distributions of zonal and meridional winds and ensemble averages of zonal pseudomomentum fluxes (F_p) and zonal and vertical wave numbers (k and m) for OGWs at z=38 km in the 4D and 2D experiments at 00:00 UTC on 20 January 2009. Zonal and meridional winds demonstrate that PWs of zonal wave number 2 (ZWN2) are significantly enhanced in the NH stratosphere. The PWs are accompanied by large spatial gradients of horizontal winds in the NH mid- to high-latitude regions. As shown in the previous section, zonal F_p and k can be significantly changed in association with zonal gradients of horizontal winds.

The horizontal structure of zonal F_p is significantly different between the 4D and 2D. First of all, zonal F_p of OGWs in the 4D is widespread in the NH mid- to high-latitude regions, and it is significantly enhanced in some particular areas where large zonal gradients of horizontal winds appear. In the 4D, nonzero zonal F_p is newly found over western Europe, the north of northern Europe, the north of the Kamchatka Peninsula, the west of North America, and over Greenland. This regional difference in nonzero zonal F_p indicates the effects of horizontal propagation. In polar regions and west of North America, zonal F_p of OGWs is directed eastward in the 4D. This eastward F_p is related to positive k induced by refraction or curvature effects. In the 2D, eastward F_p of OGWs appears only over Greenland where F_p is eastward from launch levels.

The increase in the magnitudes of zonal F_p and k appears together in narrow areas between southward and northward winds (80^∘ E and 80–140^∘ W around 30–70^∘ N), between westward and eastward winds (20–60^∘ E and 30^∘ N), and along the polar eastward jets. This enhancement is spatially correlated with the two zonal shear forcing and curvature terms in the k equation (not shown), and therefore it is thought of as being locally induced by the wind shear and curvature terms. In these narrow areas, the magnitudes of k and m become O(10⁻⁴–10⁻³ rad m⁻¹) and O(10⁻²–10⁻¹ rad m⁻¹), respectively. That is, zonal and vertical wavelengths can be as small as O(1–10 km) and O(10–100 m), respectively. In some areas between southward and northward winds (60^∘ E and 130^∘ W), m is positive (i.e., c_gz<0), positive k is large, and as a result, westward F_p is enhanced. This situation shows that large increases in F_p can occur locally as GWs propagate downward, experiencing substantial horizontal refraction.

Figure 10 Same as Fig. 9 except for NOGW_SA11.

Horizontal distributions of zonal F_p, k, and m for NOGW_SA11's at z=38 km at 00:00 UTC on 20 January 2009 (Fig. 10) can also be accounted for in similar ways as in the case of OGWs shown in Fig. 9. Statistically significant differences in zonal F_p between the 4D and 2D are mostly found in the NH mid- to high-latitude regions. Similar to OGWs, zonal F_p of NOGWs is significantly enhanced in regions where the zonal gradients of horizontal wind components and curvature effects are large (not shown). Enhancement of westward F_p near 60^∘ E (north of Canada) is largely due to the zonal shear term of zonal (meridional) winds. Along narrow regions where large zonal gradients of meridional winds appear (Fig. 10b), positive k is substantially increased, m is positive (i.e., c_gz<0), and as a result, zonal F_p is directed westward. In these regions, westward F_p becomes significantly large, and zonal and vertical wavelengths of NOGWs can also be O(1–10 km) and O(10–100 m), respectively, owing to horizontal refraction related to the zonally varying meridional wind.

Figure 11 Longitude–latitude cross sections of (a) horizontal wind (U), ensemble means of (b) the horizontal group velocity (c_gh), (c) the horizontal intrinsic group velocity (c_gh−U), (d) the vertical component of the group velocity (c_gz), and (e) ratio of intrinsic frequency ($\widehat{\mathit{\omega }}$) to Coriolis parameter ($\left|f\right|$) for OGW at z=38 km in the 4D experiment at 00:00 UTC on 20 January 2009. Hatched areas indicate regions where the paired and two-tailed t test for the 4D and 2D experiments gives p values larger than 0.05 (i.e., no statistical significance at the level of 0.05). For hatching over horizontal vector fields (b–c), the mean value of p values for the zonal and meridional components is used.

Figure 11 shows longitude–latitude distributions of horizontal wind velocity (U) and ensemble averages of the horizontal group velocity (c_gh), horizontal intrinsic group velocity (c_gh−U), vertical component of the group velocity (c_gz), and ratio of intrinsic frequency to Coriolis parameter ($\widehat{\mathit{\omega }}/\left|f\right|$) for OGWs at z=38 km in the 4D experiment at 00:00 UTC on 20 January 2009. Distributions of the horizontal group velocity and horizontal wind velocity look similar to each other, but in fact they exhibit clear differences as shown in the intrinsic group velocity.

Over the Tibetan Plateau and Rocky Mountains, the horizontal intrinsic group velocity is not close to zero. Therefore, in these regions, OGWs propagate relative to the horizontal winds. Meanwhile, in the meridionally elongated regions near 60^∘ E, 130^∘ W, and 100–140^∘ E; in the zonally elongated regions around the longitude of 180^∘ between 60 and 80^∘ N; and in isolated regions near (60^∘ N, 160^∘ W), (50^∘ N, 90^∘ W), and (70^∘ N, 60^∘ W), the magnitude of the horizontal intrinsic group velocity is roughly close to zero (see white areas in Fig. 11c). Therefore, OGW packets roughly move at the large-scale horizontal wind velocity in these areas. That is, OGW packets behave like tracers advected by the horizontal winds. The stratospheric jets (Fig. 11a) meander around these regions, but the jet cores are somewhat displaced with respect to the regions except for the meridionally elongated regions near 60^∘ E and 130^∘ W. Vertical group velocity components (c_gz) are at least 1 order of magnitude smaller than the magnitude of the horizontal intrinsic group velocity [O(1–10 m s⁻¹)] in these elongated and isolated regions. In regions near (60^∘ N, 160^∘ W), (70^∘ N, 60^∘ W), and (60^∘ N, 60^∘ E) among the regions of the small horizontal intrinsic group velocity, c_gz is particularly close to 0. Relatively large positive c_gz of O(1–10 m s⁻¹) are found near 120 and 40^∘ W, 20–40^∘ E, and 80–100^∘ E. These regions are roughly found in the peripheries of the regions of the small horizontal intrinsic group velocity.

Equation (13) shows the magnitude of the 3D intrinsic group velocity as a function of wave numbers and intrinsic frequency under the Boussinesq approximation (m²≫α²) as follows:

$$\begin{array}{}\text{(13)}& \left|{\widehat{\mathit{c}}}_{\mathrm{g}}\right|={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }}}\sqrt{{\displaystyle \frac{(N^{\mathrm{2}}-{\widehat{\mathit{\omega }}}^{\mathrm{2}})({\widehat{\mathit{\omega }}}^{\mathrm{2}}-f^{\mathrm{2}})}{{\widehat{\mathit{\omega }}}^{\mathrm{2}}}}},\end{array}$$

where κ is the magnitude of the 3D wave number vector ($\mathit{\kappa }=\sqrt{k^{\mathrm{2}}+l^{\mathrm{2}}+m^{\mathrm{2}}}$). From Eq. (13), it is clear that the magnitude of the intrinsic group velocity approaches zero in case κ significantly increases.

It is already seen that the magnitude of zonal and vertical wave numbers (Fig. 9) are substantially increased in the narrow and elongated regions around 60^∘ E and 130^∘ W and along the axis of the polar eastward jets. In these regions, (c_gλ, c_gϕ) ≈ (U, V) as shown in Fig. 11. For small $\left|{\widehat{\mathit{c}}}_{\mathrm{g}}\right|$, intrinsic frequency need not be necessarily small (Eq. 13), although the ratio $\widehat{\mathit{\omega }}/\left|f\right|$ approaches the limiting value ($\sqrt{\mathrm{2}}$) as the magnitude of the wave number approaches infinity (Bühler, 2014). The ratios, in fact, have a broad range of values in the regions of the small $\left|{\widehat{\mathit{c}}}_{\mathrm{g}}\right|$. The ratios are quite small (1–2) over the west of Greenland and the west of Alaska, but they are quite large (10–20) along 130^∘ W east of Alaska and the north of western Russia.

The significant enhancement of zonal F_p (=kc_gzA) in the elongated areas along 60^∘ E and along 80^∘ N around the longitude of 180^∘ implies an increase in zonal GW pseudomomentum (kA) since c_gz's are not particularly increased. Increases in the zonal pseudomomentum in these regions where the intrinsic group velocity is quite small indicates that there is a possibility of occurrence of the wave capture (Bühler and McIntyre, 2005) in the elongated areas during the evolution of the SSW event. When $\left|{\widehat{\mathit{c}}}_{\mathrm{g}}\right|$ is small in the highly strained large-scale flow, GW packets behave like tracers, and their shape is also substantially stretched, which leads to an increase in wave numbers and thus pseudomomentum. For GW packets in slowly varying mean flows, the change in the pseudomomentum should be balanced by a change in a quantity called impulse, defined by the mean flow vortex structure away from GW packets, unless there are external forces (Bühler, 2014). Hence, the enhanced pseudomomentum of GWs captured by distorted mean flows can cause a change in vortical motions far from GW packets. However, the surged pseudomomentum in the narrow areas may not lead to significant local change in mean flows even when the GWs are dissipated, as described by Bühler (2014).

As shown in Figs. 4–5, zonally averaged eastward F_p in the NH UMLT is enhanced in the 4D, and it is related to positive k. At z=92 km, differences in zonal F_p of NOGWs between the 4D and 2D appear over broad longitudes in the NH polar region north of 60^∘ N (Fig. S7). In the 4D, eastward F_p in the NH polar region is correlated with positive k, and m's are negative (c_gz>0). The positive k is not locally induced by zonal shear or curvature terms near z=92 km (not shown), and it seems to be gradually acquired as NOGWs propagate upward through the middle atmosphere.

Figure 12 Time–height cross sections of (a) zonal wind (U), and (b–e) ensemble means of zonal pseudomomentum fluxes (F_p) averaged over 30–90^∘ N for OGW and NOGW_SA11 in the (top row) 4D and (bottom row) 2D experiments. OGW F_p is multiplied by the efficiency factor (0.125). Contour interval for zonal winds is 10 m s⁻¹, and negative values are plotted in broken lines. Hatched areas over the zonal F_p indicates regions where the paired and two-tailed t test for the 4D and 2D experiments give p values larger than 0.05 (i.e., no statistical significance at the level of 0.05).

5.3 Time variations of pseudomomentum fluxes

Figure 12 shows time–height cross sections of zonal mean zonal wind and ensemble averages of zonal F_p of OGWs and NOGW_SA11's averaged between 30 and 90^∘ N in the 4D and 2D experiments from 8 January to 2 February 2009. Time–height distributions of zonal F_p are also quite different between the 4D and 2D. First, westward F_p in the 4D is significantly enhanced in the upper stratosphere from 13 January, about 10 d before the onset date (24 January). Increase of westward F_p is larger in OGWs than in NOGWs. Second, there is larger eastward F_p of OGWs and NOGWs in the middle stratosphere (z=30–40 km) and lower thermosphere (z=120–160 km) in the 4D from the early stage of the SSW evolution. Third, eastward F_p is substantially enhanced in the upper mesosphere a few days earlier than the onset date in the 4D. Fourth, eastward F_p of OGWs is increased in the 4D in the USLM around z=40 km in the recovery phase after the onset.

Enhanced westward F_p in the upper stratosphere before onset is clearly related to the surge of pseudomomentum associated with GWs captured by vortical mean flows. The other enhancements in the magnitude of zonal F_p in the 4D, however, do not seem to be related to the wave capture, and they are likely due to GWs in the 4D that propagate upward better, avoiding filtering or saturation in the lower atmosphere. In the 2D, significant eastward F_p in the lower thermosphere is not found before the onset, which may indicate that vertical propagation is quite restrictive in the 2D. Eastward F_p of NOGWs in the USLM near the onset occurs in both the 4D and 2D, but eastward F_p of OGWs appears only in the 4D. The eastward F_p is increased around z=40 km below the recovered eastward jets a few days after the onset date in the 4D for both OGWs and NOGWs. For OGWs, this enhanced eastward F_p is induced by upward-propagating OGWs with eastward F_p from source layers, since westward winds prevail from the ground to z=40–50 km for several days after the onset. This enhanced eastward F_p, if it exists, may induce more rapid recovery of the stratospheric jets, accelerating the downward movement of the ESs.

Close inspection of Fig. 12 for the early stage of the SSW evolution indicates that the differences between the 4D and 2D begin from the middle stratosphere around 10–11 January 2009. Albers and Birner (2014) showed that this early period can be important in the development of the SSW, especially in terms of the interaction between PWs and zonally asymmetric GW forcing. In fact, results of the 4D experiments in this study demonstrate a possibility of strong interactions between GWs and distorted vortex (or PWs) in the early period, while such a possibility seems to be low in the 2D. In the 4D, it can also be said that a possibility of a positive feedback exists between GW pseudomomentum fluxes and vortex evolution in that the large-scale flow can finally evolve into a structure of ZWN2.

Figure 13 Relative vorticity at 5 hPa and zonal pseudomomentum fluxes (F_p) for OGW at z=38 km at 00:00 UTC on (a, e) 11 (b, f), 15 (c, g) 17, and (d, h) 19 January in 2009 for the (top row) 4D and (bottom row) 2D experiments. OGW F_p is multiplied by the efficiency factor (0.125). Contour interval for relative vorticity is $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹, and negative values are plotted in broken lines. Hatched areas over the zonal F_p indicate regions where the paired and two-tailed t test for the 4D and 2D experiments give p values larger than 0.05 (i.e., no statistical significance at the level of 0.05). Latitudinal grids are plotted for every 10 from 20^∘ N.

Figure 13 shows time evolutions of relative vorticity at 5 hPa and zonal pseudomomentum fluxes (F_p) for OGWs at z=38 km at 00:00 UTC from 11 to 19 January in 2009 for the 4D and 2D experiments. On 11 January, relative vorticity exhibits a slightly elliptical shape, and PW activity looks relatively weak. In fact, there is almost no PW activity in the NH stratosphere from 1 to 10 January, as shown in Albers and Birner (2014). Although PW activity is weak on 11 January, the zonal F_p of OGWs already exhibits the longitudinal structure of ZWN2. The ZWN2 structure is much more enhanced in the 4D in association with the eastward F_p of OGWs in wide areas over the northern Pacific and northern Atlantic oceans. The OGW F_p in the 2D is confined only over the mountainous regions, but such a restriction of the geographical structure of the OGW F_p does not exist in the 4D. A similar zonally asymmetric structure is found in the meridional F_p of OGWs at z=38 km (see Fig. 14), although there are some phase shifts in the longitudinal direction compared to the zonal F_p.

Figure 14 Same as Fig. 13 except for meridional pseudomomentum fluxes (F_p).

It is clear from Figs. 13 and 14 that time evolution of the ZWN2 structure of the F_p in the 4D correlates spatially and temporally with the evolution of the vorticity. The eastward F_p over the Atlantic Ocean on 11 January moves eastward along the edge of the positive relative vorticity, and its structure is distorted on 15 January as the mean vortex is distorted. The eastward F_p over the northern Pacific Ocean gradually disappears, moving eastward on 15 January. On 17 January, the structure of the eastward F_p near the longitude of 0^∘ (the prime meridian) is more distorted and becomes elongated in the meridional direction along the vortex edge between 0 and 30^∘ E. In the narrow elongated vortex region, the southward F_p is substantially enhanced, while the northward F_p is increased over Asia and North America. On 19 January, the eastward F_p disappears, while the westward F_p is enhanced overall in and around the vortex. At the time, the northward F_p becomes predominant around the North Pole.

On 11 January, between 40 and 60^∘ N outside of the near-circular vortex (i.e., in the regions of near-zero or weakly negative relative vorticity), the eastward F_p is found near 30^∘ W and 150^∘ E in the 4D, and the relatively weak westward F_p appears in the remaining longitudinal sectors (Fig. 13a). In order to understand the distribution of the OGW F_p and possible feedback between OGWs and vortex (or PWs) in the early period of the SSW evolution, we consider a simplified ray-tracing equation for a zonally symmetric vortex. Relative vorticity is presumed to be positive within the vortex and negative outside of the vortex. The time rate of change in the meridional wave number due to horizontal wind shear can be locally approximated in the midlatitude regions by the following:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}l}{\mathrm{d}t}}=-{\displaystyle \frac{\partial U}{\partial y}}k\approx \mathit{\zeta }k,\end{array}$$

where k and l are the zonal and meridional wave numbers, respectively, ζ is the relative vorticity, and U is the zonal wind.

For an initially given aggregate of OGWs with the eastward F_p (i.e., k>0 for upward propagation) and zero l near 30^∘ W and 150^∘ E between 40 and 60^∘ N outside of the zonally symmetric vortex, dl∕dt becomes negative (away from the North Pole) because ζ< 0 outside the vortex, and therefore the pseudomomentum vectors F_p point in the southeastern direction. The generation of the southward F_p near 30^∘ W and 150^∘ E associated with $\mathrm{d}l/\mathrm{d}t<\mathrm{0}$ is seen in Fig. 14. In case this F_p is dissipated, according to the conservation rule of the pseudomomentum and impulse (Bühler, 2014), the OGW aggregate with the southeastward F_p near 30^∘ W can induce the impulse defined by the positive vorticity over northern Europe or western Russia and the negative vorticity over the Atlantic Ocean west of Africa. The aggregate of OGWs with the southeastward F_p near 150^∘ W can generate the impulse defined by the positive vorticity over Alaska and the negative vorticity over the Pacific Ocean.

The OGWs with the southeastward F_p in fact disappear at z=48 km, and change in the OGWs with the northwestward F_p between z=38 and 48 km is not clear (see Figs. S8–S9). The vertical change in the southeastward OGW F_p can generate the positive vorticity over northern Europe or western Russia and over Alaska. This GW-induced vorticity structure can advect the preexisting positive vorticity eastward and stretch the vortex towards Russia and North America across the North Pole. This deformation of the vortex can result in enhancement of PWs with ZWN2. Indeed, the positive vortex is enhanced over Russia, the negative vortex develops west of Africa, and the positive vortex is elongated on 15 January along the meridian across 90^∘ W and 90^∘ E.

On 15 January, an aggregate of OGWs with the large southeastward F_p exists along the prime meridian at z=38 km (Figs. 13–14). These pseudomomentum fluxes almost disappear at z=48 km (Figs. S8–S9). Therefore, this change in the southeastward pseudomomentum between z=38 and 48 km yields the impulse defined by the positive vorticity west of the prime meridian and the negative vorticity east of the prime meridian. This vortex structure induced by GWs on 15 January is consistent with the vortex structure on 17 January. Contrast between the positive and negative vorticities near 30^∘ N and 0–30^∘ E becomes larger on 17 January, and the southeastward F_p is more enhanced along the vortex boundaries. This result suggests a possibility of positive feedback between GW-induced vortices and mean flow vortices (or between GW momentum forcing and PWs) in the stratosphere in the early period of the SSW evolution.

6 Summary and discussion

Effects of realistic propagation of parameterized GWs on GW pseudomomentum fluxes are investigated using a global ray-tracing model for the 2009 SSW event. Two kinds (4D – x–z and t; 2D – z and t) of ray simulations are carried out to understand propagation effects for 20 ensemble members of OGWs and NOGWs for the time period of 25 d from 8 January to 2 February 2009. In each ensemble member of OGWs and NOGWs, a single GW packet is launched at a horizontal grid point, and properties (wavelength, phase speed, propagation direction, and Reynolds stress) of GW packets are randomly chosen from a precomputed set of parameters made based on previous GWP studies for OGWs and NOGWs.

The global ray-tracing model used in this study is composed of two parts, namely ray-tracing and amplitude equations. Ray-tracing equations are formulated considering the curvature effects on the spherical Earth, and they compute the trajectory of GW packets and refraction due to spatiotemporal variations of the large-scale flow. The time evolution of the vertical flux of GW action flux is computed using the amplitude equation. In the amplitude equation, ray-tube effects associated with the geometry of neighboring rays are considered by evaluating group velocity components of GW packets at grid points. For dissipative processes, nonlinear saturation and molecular viscosity are computed along ray trajectories. These dissipations only act on the action flux without affecting GW propagation.

In realistic 4D propagation, horizontal refractions related to large-scale wind shear and curvature effects are essential compared to spatial gradients of thermodynamic large-scale properties such as stability and density. Latitude–height structure of the zonal pseudomomentum fluxes (F_p) is similar overall to that of the zonal wave numbers except for the NH UMLT region. This structural agreement indicates that most GWs propagate upward in the ray simulations in this study. The magnitude of zonal F_p, however, is locally quite different between the 4D and 2D experiments. In the 4D, westward F_p is enhanced in the UTLS in both hemispheres and in the NH upper stratosphere, and eastward F_p's are reduced in the SH USLM. In the NH UMLT, the sign of zonal F_p is reversed between the 4D and 2D. It is seen that eastward F_p in the NH UMLT for the 4D is due to GWs that can propagate upward, better avoiding critical-level filtering and saturation in the lower atmosphere. As GW packets are refracted, GW packets can be converged around the axis of the stratospheric jet. Locally increased numbers of GW packets can have some effects on the zonal F_p by a factor of 2 or less in terms of magnitude. Ray-tube effects are present in the NH upper stratosphere where planetary-scale wave activity is large, but they are not significant in the other regions. Latitudinal distribution of the number of GW packets exhibits discontinuity in the 2D, and this discontinuity can induce instability that generates PWs. Given that the discontinuity in the 2D is unrealistic compared to the 4D, the 4D formulation can help minimize spurious generation of PWs due to columnar GWPs.

In the NH upper stratosphere, westward F_p is significantly enhanced in the 4D experiment along narrow and elongated areas in the mid- to high-latitude regions where spatial variations of the large-scale winds are substantial in association with large planetary wave activity. The significant enhancement in zonal F_p is mainly due to the horizontal refraction related to the horizontal wind shear and curvature effects. In the elongated regions, the magnitude of the horizontal intrinsic group velocity is quite small, which means that GWs travel by roughly following the large-scale winds. This result indicates that the wave capture phenomena may occur along the meandering eastward jets during the evolution of the SSW. For GW packets in slowly varying mean flows, change in GW pseudomomentum due to refraction and curvature terms is balanced by change in the impulse defined by the structure of nearby mean flow vortices. The enhanced F_p of captured GWs may not directly affect local mean flows where GW packets are located, even if dissipative wave–mean interaction is involved.

Enhancement of GW F_p in the 4D experiment begins about 10 d before the SSW onset date, and it remains for several days even after the onset in the recovery phase. Significant increase in the westward F_p related to the wave capture starts 10 d before the onset date in the USLM, and enhancement of the eastward F_p in the middle stratosphere and the lower thermosphere also begins from the early stage of the SSW evolution. In the 2D experiment, vertical propagation is quite restrictive, and significant F_p is not found in the lower thermosphere before the SSW onset. In the mesosphere, eastward F_p is substantially enhanced a few days before the onset in the 4D for both OGWs and NOGWs, but in the 2D relatively weak eastward F_p is induced near the onset for NOGWs alone. In the recovery phase, eastward F_p is also enhanced around z=40 km in the 4D, which implies that recovery of the stratospheric jets would possibly be accelerated when realistic propagation is considered.

In the early stage of the 2009 SSW evolution, it is interesting that the spatial distribution of the OGW F_p exhibits a much clearer ZWN2 structure in the 4D than in the 2D. The clearer ZWN2 structure is associated with enhancement of the eastward OGW F_p in wide areas over the Atlantic and Pacific oceans. In the 4D experiment where restriction such as GW propagation in the vertical direction alone is removed, the OGW activity can be found over the oceans. This enhanced ZWM2 structure of the GW F_p in the 4D shows a possibility of strong interaction and positive feedback between GWs and the vortical mean flow. Dissipation of the enhanced eastward OGW F_p can advect or stretch the preexisting vortex, and the eastward OGW F_p can possibly be further enhanced owing to the modified vortex structure. Interactions between GWs and vortex in the stratosphere can be far more active in the 4D than in the 2D. Given that substantial sensitivity of the SSW to stratospheric conditions has been found in global models with conventional columnar GWPs (e.g., de la Camara et al., 2017), the increased degree of interaction may affect the sensitivity. However, at the moment it is hard to assess how the more active interactions between GWs and vortex (or PWs) in the 4D can modify the sensitivity of the SSW evolution to small perturbations in GW fields or vortex in the stratosphere. Sensitivity might increase in proportion to the enhanced activity of the interaction, or it might not as in the case for which total wave driving remains unchanged as the resolved wave forcing compensates for parameterized wave effects (e.g., Cohen et al., 2013).

Interpretation of results shown in this paper may depend on the gridding method designed to generate gridded model outputs. In this method, the spatial size of GW packets is assumed to be as large as horizontal and vertical grid spacings used in this study. This implicit assumption may lead to overestimation of the magnitude of the GW F_p enhanced by the horizontal refraction since severe horizontal refraction may stretch GW packets significantly (anisotropically). In this case of substantial deformation of GW packets, the packets may not occupy grid spacings entirely. Physically, the size and shape of GW packets are also important since they may affect how GWs interact with mean flows. As Bühler (2014) described, as GW packets occupy more and more spaces in the longitudinal direction, they can influence the large-scale flows where the packets are located more locally. GW packets confined to limited areas may affect the ambient flows in more nonlocal ways. In order to properly consider the size of GW packets in space and time, one may need information about how much of the GW fields are steadily generated from sources (e.g., a generation timescale Δt_g is large (small) for steady (intermittent) sources) and how much of the GW fields occupy horizontal and vertical spaces from the sources (e.g., $\left|c_{\mathrm{gh}}\right|\mathrm{\Delta }{t}_{\mathrm{g}}$ and $\left|c_{\mathrm{gz}}\right|\mathrm{\Delta }{t}_{\mathrm{g}}$).

In the present study we have not discussed how GW momentum forcing can be estimated from the ray-simulation results. As described above, consideration of realistic propagation of localized GW packets in the slowly varying large-scale flows requires GWPs to compute influences of GWs in more nonlocal ways in space and time, which violates the basic assumptions of current modeling frameworks. In SSW cases considered in this study, large-scale flows can vary rapidly in space and time, and the nonlocal approach may particularly be more important since GWs can change vortex structure located around the GWs. However, at this point, it is not straightforward to present how to estimate the nonlocal influences of GWs in a clear way. In order to consider the nonlocality in models, one might either somehow extend columnar GWPs or explicitly implement ray-tracing formulations. One way or another, further theoretical developments of GW processes seem to be necessary as long as physically based methods with minimal ad hoc treatments are preferred.

Appendix A: Detailed model description

A1 Derivation of ray-tracing equations

Local time change of ω is obtained by taking the time derivative of Eq. (3) and can be written using Eqs. (1)–(3) as follows:

$$\begin{array}{}\text{(A1)}& \partial \mathit{\omega }/\partial t=-\partial \mathrm{\Omega }/\partial \mathit{k}\cdot {\mathrm{\nabla }}_r\mathit{\omega }+\left(\partial \mathrm{\Omega }/\partial {\mathrm{\Lambda }}_n\right)\partial {\mathrm{\Lambda }}_n/\partial t,\end{array}$$

where $\partial \mathrm{\Omega }/\partial \mathit{k}$ corresponds to the group velocity c_g.

From the definition of d∕dt ($=\partial /\partial t+{\mathit{c}}_{\mathrm{g}}\cdot {\mathrm{\nabla }}_r$), it is clear that Eq. (A1) is the same as Eq. (5). By definition, ${\mathit{c}}_{\mathrm{g}}=\mathrm{d}\mathit{r}/\mathrm{d}t$ (where r=re_r). This is proved by substituting r=re_r into dr∕dt and by using $\partial \mathit{r}/\partial t=\mathrm{0}$, $\partial \mathit{r}/\partial \mathit{\lambda }=r\cos \mathit{\varphi }{\mathit{e}}_{\mathit{\lambda }}$, $\partial \mathit{r}/\partial \mathit{\varphi }=r{\mathit{e}}_{\mathit{\varphi }}$, and $\partial \mathit{r}/\partial r={\mathit{e}}_{\mathrm{r}}$ (i.e., $\mathrm{d}\mathit{r}=r\cos \mathit{\varphi }\mathrm{d}\mathit{\lambda }{\mathit{e}}_{\mathit{\lambda }}+r\mathrm{d}\mathit{\varphi }{\mathit{e}}_{\mathit{\varphi }}+\mathrm{d}r{\mathit{e}}_{\mathrm{r}}$). As a result, a trajectory of a wave packet is described as follows:

$$\begin{array}{}\text{(A2)}& \left(r\cos \mathit{\varphi }\mathrm{d}\mathit{\lambda }/\mathrm{d}t,r\mathrm{d}\mathit{\varphi }/\mathrm{d}t,\mathrm{d}r/\mathrm{d}t\right)=\left(c_{\mathrm{g}\mathit{\lambda }},c_{\mathrm{g}\mathit{\varphi }},c_{\mathrm{g}r}\right).\end{array}$$

Local time change of k is obtained by taking the time derivative of Eq. (2) and is written using Eqs. (1)–(3) as follows:

$$\begin{array}{}\text{(A3)}& \partial \mathit{k}/\partial t=-\partial \mathrm{\Omega }/\partial k\cdot {\mathrm{\nabla }}_{r}k-\partial \mathrm{\Omega }/\partial {\mathrm{\Lambda }}_n{\mathrm{\nabla }}_r{\mathrm{\Lambda }}_n,\end{array}$$

where $\partial \mathrm{\Omega }/\partial \mathit{k}\cdot {\mathrm{\nabla }}_r\mathit{k}$ is expressed using the summation index as $\partial \mathrm{\Omega }/\partial k_i{\mathrm{\nabla }}_r{k}_i$ because the two k's contract with each other.

Since $k_i={\mathit{e}}_i\cdot \mathit{k}$, $\partial \mathrm{\Omega }/\partial k_i{\mathrm{\nabla }}_r{k}_i=\partial \mathrm{\Omega }/\partial k_i{\mathrm{\nabla }}_r({\mathit{e}}_i\cdot \mathit{k})$, Eq. (A3) becomes the following:

$$\begin{array}{}\text{(A4)}& \mathrm{d}\mathit{k}/\mathrm{d}t=-c_{\mathrm{g}i}\left({\mathrm{\nabla }}_r{\mathit{e}}_i\right)\cdot \mathit{k}-\partial \mathrm{\Omega }/\partial {\mathrm{\Lambda }}_n{\mathrm{\nabla }}_r{\mathrm{\Lambda }}_n.\end{array}$$

Here, c_gi(∇_re_i)⋅k should be 0 for invariance with respect to choice of coordinate system (Sect. 2.2). Consequently, Eq. (A4) is reduced to the equation for k in Eq. (4).

The constraint $c_{\mathrm{g}i}\left({\mathrm{\nabla }}_r{\mathit{e}}_i\right)\cdot k=\mathrm{0}$ indicates that the following two relations should always be satisfied in the following sphere:

$$\begin{array}{}\text{(A5)}& k{c}_{\mathrm{g}\mathit{\varphi }}\tan \mathit{\varphi }+m{c}_{\mathrm{g}\mathit{\lambda }}=l{c}_{\mathrm{g}\mathit{\lambda }}\tan \mathit{\varphi }+k{c}_{\mathrm{g}r},\end{array}$$

and

$$\begin{array}{}\text{(A6)}& l{c}_{\mathrm{g}r}=m{c}_{\mathrm{g}\mathit{\varphi }}.\end{array}$$

Note that these relations are derived from spatial variations of the basis vectors (i.e., ∇_re_i).

Substituting $\mathit{k}=k{\mathit{e}}_{\mathit{\lambda }}+l{\mathit{e}}_{\mathit{\varphi }}+m{\mathit{e}}_{\mathrm{r}}$ into Eq. (A4), where c_gi $\left({\mathrm{\nabla }}_r{\mathit{e}}_i\right)\cdot \mathit{k}=\mathrm{0}$], gives component forms of Eq. (A4) as follows:

$$\begin{array}{}\text{(A7)}& {\displaystyle }\mathrm{d}k/\mathrm{d}t=-\left({\mathrm{\Omega }}_n/h_{\mathit{\lambda }}\right)\partial {\mathrm{\Lambda }}_n/\partial \mathit{\lambda }+C_{k\mathrm{1}}+C_{k\mathrm{2}},\text{(A8)}& {\displaystyle }\mathrm{d}l/\mathrm{d}t=-\left({\mathrm{\Omega }}_n/h_{\mathit{\varphi }}\right)\partial {\mathrm{\Lambda }}_n/\partial \mathit{\varphi }+C_{l\mathrm{1}}+C_{l\mathrm{2}},\end{array}$$

and

$$\begin{array}{}\text{(A9)}& \mathrm{d}m/\mathrm{d}t=-\left({\mathrm{\Omega }}_n/h_{\mathrm{r}}\right)\partial {\mathrm{\Lambda }}_n/\partial r+C_{m\mathrm{1}}+C_{m\mathrm{2}},\end{array}$$

where ${\mathrm{\Omega }}_n=\partial \mathrm{\Omega }/\partial {\mathrm{\Lambda }}_n$, h_λ=rcos ϕ, h_ϕ=r, and h_r=1. Terms denoted by C represent the following curvature effects: $C_{k\mathrm{1}}=\left(l{c}_{\mathrm{g}\mathit{\lambda }}\tan \mathit{\varphi }\right)/r$, $C_{k\mathrm{2}}=-m{c}_{\mathrm{g}\mathit{\lambda }}/r$, $C_{l\mathrm{1}}=-\left(k{c}_{\mathrm{g}\mathit{\lambda }}\tan \mathit{\varphi }\right)/r$, $C_{l\mathrm{2}}=-m{c}_{\mathrm{g}\mathit{\varphi }}/r$, $C_{m\mathrm{1}}=k{c}_{\mathrm{g}\mathit{\lambda }}/r$, and $C_{m\mathrm{2}}=l{c}_{\mathrm{g}\mathit{\varphi }}/r$.

From Eqs. (A7)–(A9), it can be shown that the magnitude of a 3D wave number vector is invariant with respect to the Earth's curvature by multiplying Eq. (A7) by k, Eq. (A8) by l, and Eq. (A9) by m and by adding all these results together.

In the model, Eqs. (4) and (5) are approximated for the shallow atmosphere, and for the dispersion relation Eq. (6), they can be written in a component form as follows:

$$\begin{array}{}\text{(A10)}& {\displaystyle }\mathrm{d}\mathit{\lambda }/\mathrm{d}t=\left[U+k{N}_{\widehat{\mathit{\omega }}}^{\mathrm{2}}/\left(\widehat{\mathit{\omega }}{\mathit{\sigma }}^{\mathrm{2}}\right)\right]/h_{\mathit{\lambda }}=c_{\mathrm{g}\mathit{\lambda }}/h_{\mathit{\lambda }},\text{(A11)}& {\displaystyle }\mathrm{d}\mathit{\varphi }/\mathrm{d}t=\left[V+l{N}_{\widehat{\mathit{\omega }}}^{\mathrm{2}}/\left(\widehat{\mathit{\omega }}{\mathit{\sigma }}^{\mathrm{2}}\right)\right]/h_{\mathit{\varphi }}=c_{\mathrm{g}\mathit{\varphi }}/h_{\mathit{\varphi }},\text{(A12)}& {\displaystyle }\mathrm{d}z/\mathrm{d}t=-m{\widehat{\mathit{\omega }}}_{\mathrm{f}}^{\mathrm{2}}/\left(\widehat{\mathit{\omega }}{\mathit{\sigma }}^{\mathrm{2}}\right)=c_{\mathrm{g}z},\text{(A13)}& {\displaystyle }\mathrm{d}k/\mathrm{d}t=-\left(k{U}_{\mathit{\lambda }}+l{V}_{\mathit{\lambda }}+M_{\mathit{\lambda }}\right)/h_{\mathit{\lambda }}+C_{k\mathrm{1}},\text{(A14)}& {\displaystyle }\mathrm{d}l/\mathrm{d}t=-\left(k{U}_{\mathit{\varphi }}+l{V}_{\mathit{\varphi }}+M_{\mathit{\varphi }}\right)/h_{\mathit{\varphi }}+C_{l\mathrm{1}},\text{(A15)}& {\displaystyle }\mathrm{d}m/\mathrm{d}t=-\left(k{U}_z+l{V}_z+M_z\right),\end{array}$$

and

$$\begin{array}{}\text{(A16)}& \mathrm{d}\mathit{\omega }/\mathrm{d}t=+\left(k{U}_t+l{V}_t+M_t\right),\end{array}$$

where $N_{\widehat{\mathit{\omega }}}^{\mathrm{2}}=N^{\mathrm{2}}-{\widehat{\mathit{\omega }}}^{\mathrm{2}}$; ${\widehat{\mathit{\omega }}}_{\mathrm{f}}^{\mathrm{2}}={\widehat{\mathit{\omega }}}^{\mathrm{2}}-f^{\mathrm{2}}$; ${\mathit{\sigma }}^{\mathrm{2}}=k^{\mathrm{2}}+l^{\mathrm{2}}+m^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}$, h_λ=acos ϕ, and h_ϕ=a; U_λ, U_ϕ, U_z, and U_t (V_λ; and V_ϕ, V_z, and V_t) denote the partial derivatives of U (V) with respect to λ, ϕ, z, and t, respectively. $C_{k\mathrm{1}}=\left(l{c}_{\mathrm{g}\mathit{\lambda }}\tan \mathit{\varphi }\right)/a$ and $C_{l\mathrm{1}}=-\left(k{c}_{\mathrm{g}\mathit{\lambda }}\tan \mathit{\varphi }\right)/a$.

In Eqs. (A13)–(A16), terms starting with M, with effects of background medium properties other than U and V, are given by the following:

$$\begin{array}{}\text{(A17)}& {\displaystyle }{M}_{\mathit{\lambda }}=\mathrm{1}/\left(\mathrm{2}\widehat{\mathit{\omega }}{\mathit{\sigma }}^{\mathrm{2}}\right)\left(k_{\mathrm{h}}^{\mathrm{2}}{N}_{\mathit{\lambda }}^{\mathrm{2}}-{\widehat{\mathit{\omega }}}_{\mathrm{f}}^{\mathrm{2}}{\mathit{\alpha }}_{\mathit{\lambda }}^{\mathrm{2}}\right),\text{(A18)}& {\displaystyle }{M}_{\mathit{\varphi }}=\mathrm{1}/\left(\mathrm{2}\widehat{\mathit{\omega }}{\mathit{\sigma }}^{\mathrm{2}}\right)\left(k_{\mathrm{h}}^{\mathrm{2}}{N}_{\mathit{\varphi }}^{\mathrm{2}}+m_{\mathit{\alpha }}^{\mathrm{2}}{f}_{\mathit{\varphi }}^{\mathrm{2}}-{\widehat{\mathit{\omega }}}_{\mathrm{f}}^{\mathrm{2}}{\mathit{\alpha }}_{\mathit{\varphi }}^{\mathrm{2}}\right),\text{(A19)}& {\displaystyle }{M}_z=\mathrm{1}/\left(\mathrm{2}\widehat{\mathit{\omega }}{\mathit{\sigma }}^{\mathrm{2}}\right)\left(k_{\mathrm{h}}^{\mathrm{2}}{N}_z^{\mathrm{2}}-{\widehat{\mathit{\omega }}}_{\mathrm{f}}^{\mathrm{2}}{\mathit{\alpha }}_z^{\mathrm{2}}\right),\end{array}$$

and

$$\begin{array}{}\text{(A20)}& M_t=\mathrm{1}/\left(\mathrm{2}\widehat{\mathit{\omega }}{\mathit{\sigma }}^{\mathrm{2}}\right)\left(k_{\mathrm{h}}^{\mathrm{2}}{N}_t^{\mathrm{2}}-{\widehat{\mathit{\omega }}}_{\mathrm{f}}^{\mathrm{2}}{\mathit{\alpha }}_t^{\mathrm{2}}\right),\end{array}$$

where $N_{\mathit{\lambda }}^{\mathrm{2}}$, $N_{\mathit{\varphi }}^{\mathrm{2}}$, $N_z^{\mathrm{2}}$, and $N_t^{\mathrm{2}}$ (${\mathit{\alpha }}_{\mathit{\lambda }}^{\mathrm{2}}$, ${\mathit{\alpha }}_{\mathit{\varphi }}^{\mathrm{2}}$, ${\mathit{\alpha }}_z^{\mathrm{2}}$, and ${\mathit{\alpha }}_t^{\mathrm{2}}$) denote the partial derivatives of N² (α²) with respect to λ, ϕ, z, and t, respectively; $f_{\mathit{\varphi }}^{\mathrm{2}}=\partial f^{\mathrm{2}}/\partial \mathit{\varphi }$; $k_{\mathrm{h}}^{\mathrm{2}}=k^{\mathrm{2}}+l^{\mathrm{2}}$; and $m_{\mathit{\alpha }}^{\mathrm{2}}=m^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}$.

Under the shallow atmosphere approximation (Phillips, 1966) where curvature terms related to vertical movements are ignored, there is no relation corresponding to Eq. (A6), and Eq. (A5) is reduced to the following:

$$\begin{array}{}\text{(A21)}& k{c}_{\mathrm{g}\mathit{\varphi }}\tan \mathit{\varphi }=l{c}_{\mathrm{g}\mathit{\lambda }}\tan \mathit{\varphi }.\end{array}$$

Using Eqs. (A13) and (A14), it can be proved that the magnitude of the horizontal wave number vector is invariant with respect to curvature effects as in Eqs. (A7)–(A9).

A2 Effects of viscosity and diffusivity on GWs

Viscous damping and thermal diffusion terms for GWs can be obtained by linearizing the viscosity term (derived from the symmetric stress tensor) in the Navier–Stokes equation and the diffusion term in thermodynamic energy equation (see Kundu, 1990; Vadas and Fritts, 2005) as follows:

$$\begin{array}{}\text{(A22)}& \mathit{\nu }\left[{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{v}}^{\prime }+(\mathrm{1}/\mathrm{3})\mathrm{\nabla }\left(\mathrm{\nabla }\cdot {\mathit{v}}^{\prime }\right)\right],\end{array}$$

and

$$\begin{array}{}\text{(A23)}& \left(\mathit{\nu }/\mathrm{Pr}\right)\left(\mathrm{1}/\stackrel{\mathrm{‾}}{T}\right){\mathrm{\nabla }}^{\mathrm{2}}{T}^{\prime },\end{array}$$

where v^′ [= (u^′, v^′, w^′)] is the 3D perturbation wind vector; u^′, v^′, and w^′ are the zonal, meridional, and vertical perturbation wind components, respectively; T^′ is the temperature perturbation; $\stackrel{\mathrm{‾}}{T}$ is the background temperature; ν is the kinematic viscosity; and Pr is the Prandtl number.

In the viscosity terms in Eq. (A22), $(\mathit{\nu }/\mathrm{3})\mathrm{\nabla }\left(\mathrm{\nabla }\cdot {\mathit{v}}^{\prime }\right)$ is ignored by assuming that GW vertical wavelengths are much smaller than 4πH (Vadas and Fritts, 2005), where H is the density scale height. Diffusivity term Eq. (A23) is reduced to $(\mathit{\nu }/\mathrm{Pr}){\mathrm{\nabla }}^{\mathrm{2}}{b}^{\prime }$ (where $b^{\prime }={\mathit{\theta }}^{\prime }/\stackrel{\mathrm{‾}}{\mathit{\theta }}$, and θ^′, and $\stackrel{\mathrm{‾}}{\mathit{\theta }}$ are the perturbation and background potential temperatures, respectively), by neglecting pressure perturbations and spatiotemporal variations of background variables compared to GW phase variations.

Viscous damping and thermal diffusion may affect propagation of GWs through modification of the dispersion relation (Vadas and Fritts, 2005) and amplitudes, but in this model the effects on amplitudes are only considered. In order to obtain a closed expression for τ_dis in Eq. (7), following Marks and Eckermann (1995), the approximated damping terms [ν∇²v^′ and $(\mathit{\nu }/\mathrm{Pr}){\mathrm{\nabla }}^{\mathrm{2}}{b}^{\prime }$] are modified, albeit somewhat arbitrarily, to a density-weighted form:

$$\begin{array}{}\text{(A24)}& K\left[{\partial }^{\mathrm{2}}{\mathit{\chi }}^{\prime }/\partial x^{\mathrm{2}}+{\partial }^{\mathrm{2}}{\mathit{\chi }}^{\prime }/\partial y^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^{-\mathrm{1}}\partial /\partial z\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}\partial {\mathit{\chi }}^{\prime }/\partial z\right)\right],\end{array}$$

where χ^′ is u^′, v^′, w^′, or b^′; K is either the kinematic viscosity (ν) or the thermal diffusivity (ν∕Pr).

After substituting plane wave solutions such as ${\mathit{\chi }}^{\prime }=e^{z/\left(\mathrm{2}H\right)}\widehat{\mathit{\chi }}{e}^{i(kx+ly+mz-\mathit{\omega }t)}$ into Eq. (A24), derivation of equations for GW energy and action averaged over phases gives the right-hand side of Eq. (7) the following:

$$\begin{array}{}\text{(A25)}& -\mathrm{2}\mathit{\nu }{\mathit{\sigma }}^{\mathrm{2}}\left(X+Y_{\mathrm{pr}}\right)/\left(X+Y\right)A,\end{array}$$

where $X=({\widehat{\mathit{\omega }}}^{\mathrm{2}}+f^{\mathrm{2}})$ $(k^{\mathrm{2}}+l^{\mathrm{2}})/$ $({\widehat{\mathit{\omega }}}^{\mathrm{2}}-f^{\mathrm{2}}{)}^{\mathrm{2}}$, Y_pr= $({\widehat{\mathit{\omega }}}^{\mathrm{2}}+{\mathrm{Pr}}^{-\mathrm{1}}{N}^{\mathrm{2}})$ $(m^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}})/$ $(N^{\mathrm{2}}-{\widehat{\mathit{\omega }}}^{\mathrm{2}}{)}^{\mathrm{2}}$, and $Y=({\widehat{\mathit{\omega }}}^{\mathrm{2}}+N^{\mathrm{2}})(m^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}})/$ $(N^{\mathrm{2}}-{\widehat{\mathit{\omega }}}^{\mathrm{2}}{)}^{\mathrm{2}}$.

Therefore, τ_dis becomes the following:

$$\begin{array}{}\text{(A26)}& {\mathit{\tau }}_{\mathrm{dis}}=\mathrm{1}/\left(\mathrm{2}\mathit{\nu }{\mathit{\sigma }}^{\mathrm{2}}\right)\left(X+Y\right)/\left(X+Y_{\mathrm{pr}}\right)\end{array}$$

when Pr=1, Y_pr=Y, and Eq. (A26) is reduced to Eq. (12).

A3 Details of numerical implementation

The LSODA solver employs sub-time stepping within each δt. The sub-time step is determined so that the maximum norm of relative errors can be less than 1. The relative error (e_r) of each solution (y) is defined by solver-estimated error (e) divided by a weight (w) ($e_{\mathrm{r}}=e/w$), where $w=t_{\mathrm{r}}\left|y\right|+t_{\mathrm{a}}$, and t_r and t_a are relative and absolute tolerances specified for each y, respectively. For λ, ϕ, and z (k, l, m, and ω), t_r and t_a are specified as 10⁻³ and 10⁻⁶ (10⁻⁶ and 10⁻⁹), respectively. Some sensitivity tests on thresholds are carried out, but threshold values smaller than those specified above do not give significantly different results. One example of the sensitivity tests can be found in the Supplement (Fig. S3).

In the gridding method, the horizontal projection of a 3D ray trajectory during δt is assumed to be represented by a great circle path, the shortest path between two points on sphere. For a given initial location (λ_i, ϕ_i, and z_i), time integration of the ray-tracing equations gives a final location (λ_f, ϕ_f, and z_f) after δt. Spherical arc lengths (d) from the final horizontal position to the centers (λ_c and ϕ_c) of eight horizontal grid cells adjacent to the initial horizontal position are computed using the following:

$$\begin{array}{}\text{(A27)}& d={\cos }^{-\mathrm{1}}\left(\sin {\mathit{\varphi }}_{\mathrm{c}}\sin {\mathit{\varphi }}_{\mathrm{f}}+\cos {\mathit{\varphi }}_{\mathrm{c}}\cos {\mathit{\varphi }}_{\mathrm{f}}\cos \mathit{\delta }\mathit{\lambda }\right),\end{array}$$

where $\mathit{\delta }\mathit{\lambda }=|{\mathit{\lambda }}_{\mathrm{c}}-{\mathit{\lambda }}_{\mathrm{f}}|$. Among the eight cell-center locations (λ_c, ϕ_c), one cell that gives minimum d is chosen, and then an identical procedure is repeated for eight neighboring horizontal grid cells around the chosen cell until a grid cell that contains the final horizontal position is approached. Determination of contiguous 3D grid cells between (λ_i, ϕ_i, and z_i) and (λ_f, ϕ_f, and z_f) is completed considering how many vertical grid cells the ray moves through while it passes through the chosen horizontal grid cells.

Using this gridding method, three components of group velocity are stored at the vertices of chosen grid cells between initial and final positions. In addition, various ray properties such as k, l, m, ω, $\widehat{\mathit{\omega }}$, F_A, and F_p, including forcing terms of the ray-tracing equations, are stored at the same grid vertices to generate gridded model outputs.

In the model, rays are eliminated when some criteria are satisfied after time integration for δt as follows: (i) when rays move out of the model atmosphere through top and bottom boundaries, (ii) when rays are 3 d old, (iii) when the magnitude of the pseudomomentum flux ($\left|k_{\mathrm{h}}{c}_{\mathrm{gz}}A\right|$) is less than 10⁻¹⁰, or (iv) when time integration results are numerically invalid. In the present model, rays are not eliminated owing to WKB criteria based on the finding (Sartelet, 2003) that ray theory can work remarkably well in spite of the local breakdown of scale separation between GWs and large-scale flow. For rays to be eliminated, the gridding procedure is not carried out, and thus those rays do not affect τ_def and gridded outputs.