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:

and

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:

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:

and

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:

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:

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:

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:

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:

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:

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:

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.

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 2Angular 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^∘).

Download

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 3Longitude–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.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 4Latitude–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.

Download

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 5Latitude–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).

Download

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 6Latitude–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).

Download

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 7Latitude–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).

Download

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 8Same as Fig. 7 except for the zonally averaged number of GW packets.

Download

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 9Longitude–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 10Same 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 11Longitude–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:

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 12Time–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).

Download

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 13Relative 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 14Same 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:

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.