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

Abstract. 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.


from k i = e i · k = e i · (ke λ + le φ + me r ), but it must become zero 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).

Dispersion relation
Dispersion function Ω is required to compute the ray-tracing equations, and is given by wave dispersion relation.
The large-scale flow variables (U , V , N 2 , and α 2 ) and f 2 correspond to Λ n s in Eq.
(3), and they are assumed to vary slowly in space and time with respect to GW phases.

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): where A is the GW action density [= E/ω (> 0)]; E is the phase-averaged GW energy per unit volume; τ dis (> 0) is the wave 20 dissipation timescale (see Appendix A2 for details).
In the wave-mean interaction, the vertical flux of IGW horizontal pseudomomentum (F p = c gz p h = c gz k h A, where p 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 horizontally varying medium, but |F p | does not vary owing to curvature terms, as |k h | is invariant with 5 respect to the curvature (Sect. 2.2).

Dissipation mechanisms
For dissipation of GWs, two separate processes are employed: Nonlinear wave saturation and molecular diffusion.
Nonlinear saturation is computed by forcing |F p | not to exceed values for saturated GWs in GW-induced turbulence fields (Lindzen, 1981). The saturation flux (F p,sat ) formulated under the mid-frequency approximation (Senf and Achatz, 2011, SA11 15 hereafter) is employed and can be written as  In the Northern Hemisphere (NH), the polar-night jet is already 10 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: 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 15 polar upper mesosphere.
For ray simulations, G2S data at hourly interval are spatially smoothed using the vertical 1-2-1 smoother and horizontal moving averaging on spherical surface. Horizontal averaging is done using variables within the area of a spherical cap centered at every lat-long 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 20 spherical cap's surface is equivalent to the area of a circle with a radius of 300 km [π(300 km) 2 ] on a flat surface. Smoothed G2S data are regridded at 2.5 • ×2.5 • horizontal resolution for use in ray simulations.

GW ensembles
In this study, orographic and nonorographic GWs are separately considered. Properties of orographic GWs are given by the M87 scheme. Nonorographic GWs are specified based on SA11 and Warner and McIntyre (1996)  For nonorographic GWs (NGWs), 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 (Ê) in WM96 is given by a separable function of m,ω, and the azimuth angle (ϕ) [i.e., The pseudomomentum flux spectrumρc gzÊ (m,ω, ϕ)k/ω can be written as a function of k, ω, and ϕ by multiplying the spectrum by the Jacobian factor (J = m/|k|). This spectrum is discretized to obtain 14-wave schemes 10 through numerical integrations around appropriately chosen ks and ωs for 14 sets of ϕ and c p as in SA11 in a quiescent atmosphere with a specified stability (0.01 rad s −1 ). Two 14-wave WM96 schemes are obtained by using two different values (1 and 5/3) of p in the spectrum B(ω) given by B 0 (p)ω −p . NGWs 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 NGW schemes. In 15 these schemes, horizontal propagation directions (ϕs) and ground-based phase speeds (c p s) 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 = 5/3 (WM96a), the range of λ h is broader (309-782 km) compared with SA11, and in WM96 with p = 1 (WM96b), the range is much broader (128-942 km).
Each GW has identical amount of Reynolds stress in the three schemes. For this, the stresses for GWs with c p > 20 m s −1 20 are reduced in SA11, and integration intervals of k and ω for the spectra in WM96 are appropriately adjusted. As a result, NGWs with c p s smaller (larger) than 20 m s −1 have Reynolds stress of the order of 10 −3 (10 −5 ) N m −2 . Reynolds stresses exhibit slightly more westward flux, but they are almost isotropic. The total magnitudes of Reynolds stresses used in this study are 3.6×10 −3 N m −2 in the eastward, northward and southward directions and 4.1×10 −3 N m −2 in the westward direction.
These magnitudes are comparable to the total momentum flux of 3.75×10 −3 N m −2 in each Cardinal direction at 450 hPa that 25 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.

Generation of GW ensembles
In ray simulations, 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 30 for differences between ray simulations.
The OGW scheme (M87) launches single OGW at a horizontal grid point, but NGW schemes usually specify multiple GWs.
Hence, GW ensembles are separately generated for OGWs and NGWs. 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 NGWs, single GW is selected with uniform chance from 14 discrete waves. GW ensembles are precomputed using a random-number generator for reproducibility. teristic structure is more clearly revealed in its ensemble average (Fig. 3f). The sign of zonal NGW F p s is generally opposite to that of the tropospheric zonal flows.

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

20
For each GW ensemble member, two kinds of simulations are carried out: Four-dimensional (4D) (x-z, t) and two-dimensional (2D) (z, t) experiments. In 4D experiments, GW rays propagate horizontally as well as vertically in spatially varying background media, but in 2D, rays are allowed to propagate only in the vertical direction. In both 4D and 2D, GW rays propagate through time-varying flows, and therefore modulations of the observed frequencies of GWs occur. In 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, 25 3-h averaged gridded outputs are generated. GW rays are launched every 3 h, and 3-day old rays are eliminated. These launch interval and ray life time are chosen considering computational time, the time scale 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). Additional 4D OGW and NGW 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 15 significant differences between nonzero and zero τ def s are found in relatively narrow regions around the NH jet core (75 • N and 40 km at 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 s when τ def = 0 is reduced by less than half compared with when τ def = 0. These differences are localized compared with the differences between the 4D and 2D experiments, which may indicate that ray-tube effects be relatively limited. However, it 20 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 single GW ray is launched at a grid point. the sign of zonal-mean ks is overall the same as that of the zonal F p s shown in Fig. 4, except for some regions in the NH USLM. This similarity in the signs of ks and zonal F p s 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 ks at z = 30-50 km are related to the GWs with negative c gz and westward F p (see Sect. 5.2 for details).

30
Statistically significant differences of zonal-mean ks between the 4D and 2D are also found in similar regions to those of zonal-mean zonal F p s shown in Fig. 4. In the NH polar UMLT, the sign of ks is reversed between the 4D and 2D, and the magnitude of positive ks in the 4D is 1.2-6.3 times larger than that of negative ks in the 2D for both OGWs and NGWs. In the UTLS, negative ks of NGWs are 1.4-2.4 times larger in the 4D. In the SH mesosphere, positive ks of NGWs in the mid-latitude regions are reduced by about half in the 4D, and positive ks of OGWs around 70 • S are about 1-3 times enhanced in the 4D.

35
These changes in ks in the 4D with respect to the 2D are roughly similar to those in the zonal F p s 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 ks between the 4D and 2D.
Time rate of change of ks along rays is determined by the five forcing terms [Eqs. (A13) and (A17)]. Among the 5 terms, the two zonal shear forcing terms [−k/(a cos φ)∂U/∂λ and −l/(a cos φ)∂V /∂λ] and curvature term (lc gλ tan φ/a) are predomi-5 nant in most cases. Thermodynamic forcing terms are usually 2-3 order of magnitude smaller (see Fig. S5).  Structure of ls for NGWs is coherent with the large-scale flow structure in the 4D, but it seems more or less random in the 2D. In the SH, ls of NGWs in the 4D are overall positive in the high-latitude stratosphere above z = 10 km and in the mid-latitude middle atmosphere above z = 30 km, and c gλ is overall negative due to the westward winds. This explains why  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 NGWs. In the latitude-height region of 40 • N-70 • N and 40-90 km on 15 January, zonally averaged number of rays for OGWs (NGWs) 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 outnumber GWs with eastward F p , and therefore the ray convergence may increase the 25 westward F p . In the NH mid-to high-latitude lower thermosphere (30 • N-90 • N and 120-140 km), zonally averaged number of rays for OGWs (NGWs) is about 2.5 (1.9) times larger in the 4D on 15 January. This difference is due mainly to OGWs and NGWs 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 • N-80 • N and 30-50 km). In this region, zonally averaged number of 30 rays for OGWs (NGWs) is about 0.9 (2.3) in the 4D, and it is 1.7 (1.1) times larger than in the 2D. As in 15 January, GWs with westward F p are predominant in the NH stratosphere, resulting in some enhancement of westward F p s. In the NH lower thermosphere (30 • N-90 • N and 120-140 km), zonally averaged number of rays for both OGWs and NGWs is about 1.8 times larger in the 4D. This increase is mostly attributed to OGWs and NGWs with eastward F p s, as in 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 NGWs with eastward F p (i.e., k > 0 and c pλ − U > 0, where c pλ is the ground-based zonal phase speed) can propagate better to the thermosphere in the 4D, being less dissipated in the middle atmosphere. There seems a tendency in the 4D for GWs to better elude critical-level filtering or nonlinear saturation in the lower atmosphere. This tendency may be attributed to larger degree of freedom in propagation, associated with change in 5 wavenumber 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 NGWs in the SH. More NGW 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 10 near z = 90 km in the SH where the zonal-mean zonal wind is reversed. This spread in the 4D compared with 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 account for differences between ray-based and columnar GWPs. Discussion regarding Fig. 8 indicates that convergence or divergence (spread of GW packets) effects in the 15 4D would be smaller than a factor of 2 in terms of the magnitude of F p s. These convergence and divergence effects are, however, likely relatively small compared with impacts due to change in horizontal wavenumbers shown in Figs. 4-5, although the refraction effects would not entirely lead to local change in mean flows (see Sect. 5.2). winds in the NH mid-to high-latitude regions. As shown in the previous section, zonal F p s and ks can be significantly changed in association with zonal gradients of horizontal winds.

25
Horizontal structure of zonal F p s is significantly different between the 4D and 2D. First of all, zonal F p s of OGWs in the 4D are widespread in the NH mid-to high-latitude regions, and they are significantly enhanced in some particular areas where large zonal gradients of horizontal winds appear. In the 4D, nonzero zonal F p s are newly found over the Western Europe, north of Northern Europe, north of the Kamchatka peninsula, west of North America, and over Greenland. This regional difference in nonzero zonal F p s indicates effects of horizontal propagation. In polar regions and west of North America, zonal F p s of 30 OGWs become eastward in the 4D. These eastward F p s are related to positive ks induced by refraction or curvature effects.
In the 2D, eastward F p s of OGWs appear only over Greenland and some regions along the Antarctic coastlines where F p s are eastward from launch levels.
Increase in the magnitudes of zonal F p s and ks appears together in narrow areas between southward and northward winds  Over the Tibetan plateau and North America in which the eastward jet cores are located, both components of the horizontal intrinsic group velocity are not close to zero. Therefore, in these regions, OGWs propagate relative to the horizontal winds.

35
Equation (13) shows the magnitude of the 3D intrinsic group velocity as a function of wavenumbers and intrinsic frequency under the Boussinesq approximation (m 2 α 2 ): where κ is the magnitude of the 3D wavenumber vector (κ = √ k 2 + l 2 + m 2 ). From Eq. (13), it is clear that the magnitude of the intrinsic group velocity approaches zero in case κ significantly increases.

5
It is already seen that the magnitude of zonal and vertical wavenumbers (Fig. 9)  In realistic 4D propagation, horizontal refractions related to large-scale wind shear and curvature effects are essential compared with spatial gradients of thermodynamic large-scale properties such as stability and density. Latitude-height structure of the zonal pseudomomentum fluxes (F p s) is overall similar to that of the zonal wavenumbers except for the NH UMLT region.

30
This structural agreement indicates that most of GWs propagate upward in the ray simulations in this study. The magnitude of zonal F p s, however, is locally quite different between the 4D and 2D experiments. In the 4D, westward F p s are 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 s are reversed between the 4D and 2D. It is seen that eastward F p s in the NH UMLT for the 4D are 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 number of GW packets can have some effects on the zonal F p s in a factor of 2 or less in terms of magnitude. Ray-tube effects are present in the NH upper stratosphere where planetary-scale wave activities are large, but they are not significant in the other regions. Enhancement of GW F p s in the 4D experiment begins about 10 days before the SSW onset date, and it remains several days 15 even after the onset in the recovery phase. Significant increase of the westward F p s related to the wave capture starts 10 days before the onset date in the USLM, and enhancement of the eastward F p in the lower thermosphere also begins from early stage of the SSW evolution. In the 2D experiment, vertical propagation is quite restrictive, and significant F p are not found in the lower thermosphere before the SSW onset. In the mesosphere, eastward F p s are substantially enhanced a few days before the onset in the 4D for both OGWs and NGWs, but relatively weak eastward F p s are induced near the onset in the 2D for 20 NGWs alone. In the recovery phase, eastward F p s are 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.
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 s enhanced by the horizontal 25 refraction, since severe horizontal refraction may stretch significantly anisotropically GW packets. In this case of substantial deformation of GW packets, the packets may not occupy entirely grid spacings. 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 more locally the large-scale flows where the packets are located. GW packets confined in limited areas may affect the ambient flows in more nonlocal ways. In order to 30 consider properly size of GW packets in space and time, one may need information about how much GW fields are steadily generated from sources [e.g., A generation time scale ∆t g is large (small) for steady (intermittent) sources] and how much GW fields occupy horizontal and vertical spaces from the sources (e.g., |c gh |∆t g and |c gz |∆t g ).
In the present study we have not discussed about 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 35 flows requires for 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 as 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 in a clear way how to estimate the nonlocal influences of GWs. In order to consider the nonlocality in models, one might either somehow extend columnar GWPs or 5 explicitly impliment ray-tracing formulations. One way or the other, further theoretical developments on GW processes seem to be necessary, as long as physically-based methods with minimal ad-hoc treatments are preferred. .
where ∂Ω/∂k corresponds to the group velocity c g .
Local time change of k is obtained by taking the time derivative of Eq.
(2) and is written using Eqs.
(1)-(3) as where ∂Ω/∂k · ∇ r k is expressed using summation index as ∂Ω/∂k i ∇ r k i because the two ks contract with each other.
Since k i = e i · k, ∂Ω/∂k i ∇ r k i = ∂Ω/∂k i ∇ r (e i · k), and thus Eq. (A3) becomes Here, c gi (∇ r e i ) · k should be zero 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 gi (∇ r e i ) · k = 0 indicates that the following two relations should always be satisfied on sphere: kc gφ tan φ + mc gλ = lc gλ tan φ + kc gr , and lc gr = mc gφ . (A6) 10 Note that these relations are derived from spatial variations of the basis vectors (i.e., ∇ r e i ).

10
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 Using Eqs. (A13) and (A14), it can be proved that the magnitude of horizontal wavenumber 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 Navier-Stokes equation and the diffusion term in thermodynamic energy equation (see Kundu, 1990;Vadas and Fritts, 2005) as follows: 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;T is the background temperature; ν is the kinematic viscosity; Pr is the Prandtl number.
In the viscosity terms Eq. (A22), (ν/3)∇ (∇ · v ) is ignored by assuming that GW vertical wavelengths are much smaller 25 than 4πH (Vadas and Fritts, 2005), where H is the density scale height. Diffusivity term Eq. (A23) is reduced to (ν/Pr)∇ 2 b (where b = θ /θ, θ andθ 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 dispersion relation (Vadas and Fritts, 2005) as well as amplitudes, but in this model effects on amplitudes are only considered. In order to obtain a 30 closed expression for τ dis in Eq. (7), following Marks and Eckermann (1995), the approximated damping terms [ν∇ 2 v and (ν/Pr)∇ 2 b ] are modified, albeit somewhat arbitrarily, to a density-weighted form: where χ is u , v , w , or b ; K is either the kinematic viscosity (ν) or the thermal diffusivity (ν/Pr).

A3 Details of numerical implementation
The LSODA solver employs subtime stepping within each δt. Sub-timestep 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 an weight 15 (w) (e r = e/w), where w = t r |y| + t 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 −3 , and 10 −6 (10 −6 , and 10 −9 ), respectively. Some sensitivity tests on thresholds are carried out, but threshold values smaller than 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-20 circle path, the shortest path between two points on sphere. For an initial location (λ i , φ i , z i ), time integration of the ray tracing equations gives a final location (λ f , φ f , z f ) after δt. Spherical arc lengths (ds) from the final horizontal position to the centers (λ c , φ c ) of 8 horizontal grid cells adjacent to the initial horizontal position are computed using where δλ = |λ c − λ f |. Among the 8 cell-center locations (λ c , φ c ), one cell that gives minimum d is chosen, and then identical 25 procedure is repeated for 8 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 , z i ) and (λ f , φ f , z f ) is completed considering how many vertical grid cells the ray move through while it pass 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, ω,ω, 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: (i) when rays move out of the model atmosphere through top and bottom boundaries, (ii) when rays are 3-day old, (iii) when magnitude of the pseudomo-