## Abstract

Operational gravity wave parameterization schemes in GCMs are columnar; that is, they are based on a ray-tracing model for gravity wave propagation that neglects horizontal propagation as well as refraction by horizontally inhomogeneous basic flows. Despite the enormous conceptual and numerical simplifications that these approximations provide, it has never been clearly established whether horizontal propagation and refraction are indeed negligible for atmospheric climate dynamics. In this study, a three-dimensional ray-tracing scheme for internal gravity waves that allows wave refraction and horizontal propagation in spherical geometry is formulated. Various issues to do with three-dimensional wave dynamics and wave–mean interactions are discussed, and then the scheme is applied to offline computations using GCM data and launch spectra provided by an operational columnar gravity wave parameterization scheme for topographic waves. This allows for side-by-side testing and evaluation of momentum fluxes in the new scheme against those of the parameterization scheme. In particular, the wave-induced vertical flux of angular momentum is computed and compared with the predictions of the columnar parameterization scheme. Consistent with a scaling argument, significant changes in the angular momentum flux due to three-dimensional refraction and horizontal propagation are confined to waves near the inertial frequency.

## 1. Introduction

In columnar ray-tracing schemes both the horizontal location as well as the horizontal wavenumber vector are constant along a ray. In addition, for a steady nondissipative wave train the net vertical flux of horizontal pseudomomentum is constant along a ray tube. Because the vertical flux of horizontal pseudomomentum equals the wave-induced vertical flux of mean horizontal momentum (e.g., Andrews and McIntyre 1978, section 5.2), this leads to the standard paradigm of columnar wave–mean interaction, which is the “pseudomomentum rule.” According to this rule persistent, cumulative forcing of the mean flow occurs only where waves dissipate and then the mean force is equal to the rate of dissipation of the horizontal pseudomomentum.

Columnar parameterization schemes for gravity waves (or other processes) offer enormous conceptual and numerical simplifications, the latter being of particular importance for the current generation of highly parallelized GCMs. Nevertheless, it has never been clearly established whether horizontal propagation and refraction are, indeed, negligible for atmospheric climate dynamics. This is an interesting question not only for parameterization efforts but also for understanding the dynamics of explicit gravity waves that are increasingly produced and observed in high-resolution GCMs (e.g., Plougonven and Snyder 2005). Such explicitly resolved waves may also affect data assimilation procedures.

The present study makes a step toward estimating how three-dimensional refraction of internal gravity waves affects the global vertical transport of angular momentum. To this end, we built an offline ray-tracing model in spherical geometry and applied it to output from the Third Generation Atmospheric General Circulation Model (AGCM3) of the Canadian Centre for Climate Modeling and Analysis.

The AGCM3 includes a columnar gravity wave parameterization scheme for topographic waves and multiple runs of the ray tracer were computed using the same topographic waves at the lower boundary as the parameterization scheme. In addition, further runs with artificially lowered intrinsic frequencies were computed so as to assess the sensitivity of the angular momentum transport to the details of the wave spectrum. Overall, the scientific aims of this study were twofold: to clarify the fundamental issues to do with three-dimensional refraction and wave–mean interactions in atmospheric models and to produce a clear, quantitative comparison between a three-dimensional ray-tracing model and a columnar parameterization scheme in the context of the angular momentum budget of the atmosphere.

We briefly summarize the scientific background of three-dimensional ray tracing, which was already described analytically in Jones (1969). Subsequent numerical studies noted that inclusion of three-dimensional effects in the middle atmosphere could lead to significant changes in the wave dynamics at least for some wave rays, especially those with reduced vertical group velocities (e.g., Dunkerton and Butchart 1984; Dunkerton 1984). Meridional propagation alone leads to well-known changes in the wave structure due to the latitude dependence of the Coriolis parameter, and such changes are important for the interpretation of the observed latitudinal structure of wave energy spectra (Bühler 2003). However, such latitude-dependent effects would not change the vertical flux of zonal angular momentum associated with a wave train. On the other hand, it is well known that horizontal refraction by zonally asymmetric mean fields would alter this flux (the standard argument is presented in section 2).

What is less well known is how such horizontal refraction affects the wave–mean interactions. Recently, this problem was considered theoretically by Bühler and McIntyre (2003, 2005, hereafter BM05), who demonstrated that the refractive changes in the wave-induced momentum fluxes go hand in hand with persistent changes in the horizontally inhomogeneous mean flow that caused the refraction. This extends the pioneering earlier work of Bretherton (1969). In particular, for simple idealized scenarios a new conservation law for the sum of horizontal pseudomomentum and the impulse of the layerwise distribution of potential vorticity was derived in BM05. For instance, this conservation law allows understanding how total horizontal momentum (or angular momentum in the spherical case) is conserved when three-dimensional refraction is present. Persistent, cumulative mean-flow changes due to refraction are indeed possible, and they fall outside the columnar pseudomomentum rule because they do not rely on wave dissipation in any obvious manner. These considerations motivated the present study.

The plan for the paper is as follows: The basic theory for three-dimensional wave dynamics is laid out in section 2, and the numerical ray-tracing model is described in section 3. The model was carefully tested for consistency with the spherical coordinates and in columnar mode it was calibrated against the topographic gravity wave scheme used in AGCM3. This close comparison gave us confidence that the ray tracing scheme produces realistic fluxes of angular momentum. The numerical experiments and results are discussed in sections 4 and 5.

## 2. Three-dimensional gravity wave dynamics

For simplicity, we first work in Cartesian geometry and later describe the adaptation to spherical geometry in section 2c.

### a. Ray tracing and refraction

The ray-tracing equations describe the linear evolution of a slowly varying wave train containing small-scale waves that propagate relative to a large-scale basic flow (e.g., Whitham 1974; Lighthill 1978). This is the natural setting for unresolved, subgrid-scale internal waves. Thus, all wave fields are taken to be proportional to the real part of *a*(**x**, *t*) exp[*i*Θ(**x**, *t*)], where *a* ≪ 1 is the slowly varying wave amplitude and Θ is the rapidly varying wave phase whose derivatives; that is, the wavenumber vector **k** = (*k*, *l*, *m*) = **∇**Θ and the absolute frequency *ω* = −Θ* _{t}* are again slowly varying. A suitable dispersion relation

*ω*= Ω(

**k**,

**x**,

*t*) then links

**k**and

*ω*and this serves as a first-order PDE for Θ whose characteristics are the ray-tracing equations

which form a Hamiltonian system of ODEs for the canonical variables (**x**, **k**) based on the Hamiltonian function Ω. Physically, *d*/*dt* = ∂/∂*t* + **c*** _{g}* ·

**∇**is the time derivative following a ray and

**c**

*=*

_{g}**∇**Ω is the absolute group velocity. As is well known, the spatial and temporal symmetries of the basic flow induce conservation laws along the rays; that is, if Ω(

_{k}**k**,

**x**,

*t*) is independent of one of (

*x*,

*y*,

*z*,

*t*), then the corresponding member of (

*k*,

*l*,

*m*,

*ω*) is constant along a ray. In a columnar parameterization scheme the ray is constrained not to travel horizontally, and therefore for consistency the horizontal gradient of background fields must be ignored. In this context Ω(

**k**,

**x**,

*t*) is formally independent of

*x*and

*y*and, therefore, (

*k*,

*l*) are both constant along a ray.

In the absence of forcing and dissipation the wave train amplitude along ray tubes formed by nonintersecting rays is governed by the conservation law for the *O*(*a*^{2}) wave action density *A*; that is,

Here *A* = *E*/ω̂, where *E* is the phase-averaged disturbance energy density per unit volume (Bretherton and Garrett 1968) and ω̂ is the intrinsic frequency. The wave action is conserved by virtue of the assumed scale separation between the small-scale waves and the large-scale basic flow.

Of particular interest is the pseudomomentum vector **p** = **k***A*, which arises naturally in wave–mean interaction theory. For instance, **p** appears naturally in averaged expressions for the potential vorticity (Bühler and McIntyre 1998, 2005) and, as noted above, the vertical flux of horizontal pseudomomentum **p*** _{H}* = (

*k*,

*l*)

*A*, say, equals the wave-induced vertical flux of mean horizontal momentum (e.g., Andrews and McIntyre 1978, section 5.2). In columnar parameterization schemes the horizontal components of pseudomomentum are conserved by construction because

*dk*/

*dt*=

*dl*/

*dt*= 0 in such schemes and, therefore, horizontal pseudomomentum conservation is inherited from wave action conservation. This leads to the pseudomomentum rule discussed in section 1.

On the other hand, by combining (1) and (2) it is clear in principle that a component of **p** is conserved only if the corresponding component of **k** is invariant along a ray, that is, only if the basic flow has a symmetry in the corresponding spatial direction. In particular, zonal pseudomomentum *p* = *kA* is conserved only if the basic flow is zonally symmetric. This elementary fact implies that the wave-induced vertical flux of mean zonal momentum along a ray tube is *not constant* in principle, even if the wave field is steady and nondissipating. This does not contradict the usual nonacceleration conditions based on the Eliassen–Palm flux theorems because these conditions are always derived under the assumption that the basic flow is zonally symmetric (e.g., Andrews et al. 1987). Indeed, there is no contradiction with momentum conservation either because more detailed studies of the wave–mean interactions that accompany such nondissipative changes in pseudomomentum illustrate that the changes in the pseudomomentum flux are precisely compensated for by wave-induced changes in the basic flow (Bühler and McIntyre 2003; BM05).

A specific example of how pseudomomentum can be created or destroyed by a three-dimensional basic flow is wave refraction by a shear flow. With a basic velocity field **U**(**x**, *t*) = (*U*, *V*, *W*), the dispersion function Ω has the generic Doppler-shifting form

where ω̂ is the intrinsic frequency, and the group velocity **c**_{g} = **∇**_{k}ω̂ + **U** now consists of advection by the basic flow plus intrinsic propagation relative to it. For simplicity, we restrict ω̂ to be a function of **k** only, but the general case will be treated in section 3. In this case, the symmetries of Ω are those of **U,** and we have the explicit refraction equations

Clearly, any dependence of **U** on the zonal coordinate *x* (say, due to synoptic-scale vortices, meandering jet streams, or stratospheric Rossby waves) induces concomitant changes in *k* and therefore in the vertical flux of zonal momentum.

Actually, slightly more can be predicted from (4) because refraction affects **k** (and also **p**) in a manner partially analogous to the evolution of the gradient of a passive tracer advected by **U**. Specifically, for a passive tracer *ϕ* we have

where *D*/*Dt* is the material derivative based on **U**. Thus, the evolution of **k** and **∇***ϕ* differ only due to advection by the intrinsic group velocity, which is the difference between advection along group-velocity rays in (4) and along material trajectories of the basic flow in (5). This partial analogy can be exploited in a number of ways (e.g., Badulin and Shrira 1993; BM05). In particular, by aligning **k** with a right eigenvector of **∇U** it is possible to sustain unbounded exponential growth in time of all three components of **k** (Jones 1969). This leads to a peculiar three-dimensional version of a critical layer in which the local amplitude also grows exponentially while the intrinsic group velocity, which for gravity waves is inversely proportional to |**k**| at fixed ω̂, decays exponentially to zero. In this scenario the wave pattern is eventually passively advected by the basic flow; the transition to this passive advection regime and the concomitant wave–mean interactions have been termed “wave capture” in BM05. Such dynamics is of course absent in columnar parameterization schemes, although some evidence for the occurrence of this kind of effect is available from direct numerical simulation of gravity waves with high-resolution GCMs, as reported in Plougonven and Snyder (2005).

Two simple scaling arguments demarcate the intrinsic frequency range for which three-dimensional refraction can be expected to be important. This will be the case if three-dimensional refraction can induce significant relative changes in the horizontal wavenumbers **k*** _{H}* = (

*k*,

*l*), and corresponding changes in the wave-induced vertical fluxes of horizontal momentum along the ray tube, before the ray encounters a critical layer or is otherwise dissipated. We assume that the basic flow is described by a layerwise horizontal velocity field such that |

**U**| = |(

*U*,

*V*, 0)| ∼ 𝒰 and

*L*and

*H*are horizontal and vertical length scales such that |

**∇**

_{H}**U**| ∼ 𝒰/

*L*and |

**U**

*| ∼ 𝒰/*

_{z}*H*. Then, for a ray propagating through a depth of atmosphere Δ

*r*during a residence time Δ

*t*, we have

Using ω̂ defined by the dispersion relation for gravity waves, given in Eq. (10) below, and *f* defined as the magnitude of the inertial frequency, the vertical group velocity scales as

so that as ω̂ → *f*, *dr*/*dt* → 0. This implies a longer residence time Δ*t* at fixed Δ*r* for low frequency waves. On the other hand, the relative importance of changes in horizontal and vertical wavenumber can be estimated from (1b) as

Under quasigeostrophic scaling the typical aspect ratio of the basic flow is the Prandtl ratio *H*/*L* = | *f* |/*N*, and hence horizontal refraction is comparable to vertical refraction when |*m*|/|**k*** _{H}*| ≥ |

*f*|/

*N*. Again, this scaling supports the importance of three-dimensional effects for low frequency waves. Making use again of the anelastic dispersion relationship for internal gravity waves, the ratio (7) is unity if ω̂ = 2 |

*f*|. The ratio tends to zero for ω̂ ≫ |

*f*|, suggesting that vertical refraction dominates horizontal refraction in the rest of the frequency spectrum. Both of these scaling arguments suggest that three-dimensional refraction is most important for the near-inertial part of the atmospheric wave spectrum. This dependency of horizontal refraction and propagation on intrinsic frequency has been noted in numerous ray tracing studies, such as Dunkerton (1984), Eckermann (1992), and Hertzog et al. (2001). The robust nature of the dependence is further supported by the experimental results presented in section 4.

Additionally, allowing horizontal propagation can have a large impact whenever the absolute horizontal group velocity of a wave is large. In such cases, the path of a three-dimensionally propagating ray quickly diverges from that of its columnar analog, and it moves through a different environment. Even if horizontal refraction does not play a large role, the final destination of the wave’s pseudomomentum flux could be significantly changed by three-dimensional effects. Topographic waves in two dimensions, for which columnar propagation schemes were originally developed, have absolute horizontal group velocities near zero, so this effect should be of little importance for them. But it could play a significant role for many waves in the nonorographic spectrum and for certain three-dimensional topographic waves (e.g., Smith 1980; Shutts 1998).

### b. Local amplitude prediction

The local wave amplitude is important because wave dissipation due to nonlinear wave breaking must be modeled in parameterization schemes and incipient wave breaking is sensitive to local wave amplitudes. In a columnar scheme (2) is applied to a steady and horizontally homogeneous wave train, in which case it reduces to *Aw _{g}* =

*A*

_{0}

*w*

_{g}_{0}in which

*w*is the vertical group velocity and the index denotes the known conditions at the lower boundary. From this equation

_{g}*A*, and, therefore the local wave amplitude

*a*, can be inferred using the usual polarization relations. One-dimensional caustics correspond to

*w*= 0 and these occur at reflection and critical layers where ray-tracing predicts infinite amplitudes. Now, in a three-dimensional setting, it is simplest to retain the assumption of a steady wave train and then (2) reduces to

_{g}**∇**· (

*A*

**c**

*) = 0. In principle, this equation can be solved along a bundle of rays forming a ray tube with constant net wave-action flux along the cross section of the tube.*

_{g}The kinematic ray theory of Hayes (1970) achieves this by adding differential equations for the spatial gradients of **k** to the ray tracing set and using these quantities to compute **∇** · **c*** _{g}* along a single ray. The local amplitudes along the ray may then be calculated, provided the gradient of

**k**is specified at the lower boundary. This requirement can pose a difficulty in many practical applications of ray tracing where the boundary conditions may not be specified with sufficient resolution to constrain

**∇k**.

Furthermore, even if **∇k** is well known at the boundary, the use of the Hayes ray equations is complicated by the ubiquitous appearance of caustics. At caustics neighboring rays intersect, the ray tube area vanishes locally, ray theory breaks down, and the predicted amplitude is infinite. Though it is possible to formulate the ray-tracing equations in a way that prevents singularities at caustics (e.g., Broutman et al. 2004, 2001), there is no simple way of computing the correct amplitude near a caustic in a general three-dimensional basic state. As is well known from more elaborate ray theory, the correct wave amplitude at such a caustic is not infinite and can be computed, in principle, by resolving the singularity of ray theory with a more detailed local theory. In idealized geometries, the solution involves Airy functions across the caustic (Lighthill 1978). A series of studies (Broutman 1984, 1986; Eckermann and Marks 1996; Sonmor and Klaassen 2000; Walterscheid 2000) has used the Hayes equations to investigate the formation of caustics in idealized background flows. All noted that the location and severity of amplitude amplification at caustics is extremely sensitive to the boundary conditions on **∇k** and the details of the background flow. The Hayes amplitude equations are therefore inappropriate for applications in which these details are poorly constrained, as they are in the present case. Therefore, after some experimentation with Hayes’s extended equations, we restricted our study to the crude approximation in which the local wave action density is computed by assuming that the cross-sectional area of the ray tubes is constant (e.g., Marks and Eckermann 1995).

### c. Spherical geometry

The ray-tracing equations (1)–(2) apply in spherical geometry provided that due care is taken when adapting the relevant vector expressions to the nonconstant basis vectors of the spherical coordinate system. Our resulting ray-tracing equations differ from the equations used in the Gravity Wave Regional or Global Ray Tracer (GROGRAT) model of Marks and Eckermann (1995) by several metric terms that become relevant if rays travel substantial horizontal distances. Our derivation is detailed in the model description in section 3 and the appendix.

As for the wave–mean interaction problem in spherical geometry, the computation of the detailed structure of a wave field and the corresponding local response of the mean flow is a subtle and difficult problem. We hope to report some interesting progress in this area in the near future, but in the present study we concern ourselves only with the computation of the aggregate vertical transport of angular momentum around the earth’s rotation axis by an ensemble of ray tubes. The net angular momentum transport by an individual ray tube can be equated to the cross-sectional integral of the vertical flux of angular pseudomomentum, which has density *kAr* cos *θ* (Jones 1969). Here *k* is the local zonal wavenumber, *r* is the distance to the center of the earth, and *θ* the latitude such that *r* cos *θ* is the distance from the rotation axis. The net vertical flux is

where the integral is over the horizontal cross section of the ray tube with area element *dS*. The second expression involves the cross-sectionally integrated wave action flux

and the useful approximate equality holds if *k* is essentially constant over the horizontal extent of the ray tube. For a steady nondissipating wave train *F* is constant along a ray tube and, therefore, *P* ∝ *kr* cos*θ* in this case. Again, it follows from (1) that for a zonally symmetric basic state the net flux of angular pseudomomentum *P* is constant, and this shows that the combination *kr* cos*θ* must be constant along a ray in a zonally symmetric atmosphere. This is a good test for a ray-tracing scheme in spherical geometry because (*k*, *r*, *θ*) all vary individually in this case. Finally, if the basic flow is not zonally symmetric, then *P* will not be constant along the tube and, therefore, the wave-induced flux of angular momentum will be variable in altitude even without wave dissipation. This is analogous to the nonconstant flux of zonal momentum in the Cartesian geometry discussed in section 2a.

## 3. Model description

A global ray-tracing model in spherical coordinates has been developed in Matlab. It computes the propagation, refraction, and dissipation of ray tubes through steady, three- dimensionally varying stratification and horizontal winds **U** = (*U*, *V*, 0) in a compressible deep atmosphere. The nonessential restriction to steady basic flows is common to all parameterization schemes. The model uses the dispersion relationship for inertia–gravity waves in a rotating stratified compressible atmosphere, which is

Here *N* is the buoyancy frequency, *f* the local Coriolis parameter, and *α* = −0.5*ρ*^{−1}∂*ρ*/∂*z*. This general dispersion relationship, also used in the GROGRAT ray tracing model of Marks and Eckermann (1995), is accurate over the entire frequency range | *f* | < ω̂ < *N*. The sign convention ω̂ > 0 has been chosen, which implies that *m* < 0 corresponds to gravity waves with upward propagating group velocity, *k* > 0 denotes waves with positive zonal intrinsic group velocity—and therefore with prograde zonal momentum flux—and *l* > 0 denotes northward intrinsic group velocity.

### a. Spherical ray tracing

As noted before, finding the correct ray-tracing equations in spherical coordinates is a slightly subtle matter. For instance, the wavenumber vector in these coordinates is **k** = *k*** λ̂** +

*l*

**Θ̂**+

*m*

**r̂**, where

**,**

*λ̂***Θ̂**, and

**r̂**are the usual zonal, meridional, and radial unit vectors. The components of

**k**in spherical coordinates evolve not only due to changes in

**k**but also due to the spatial variation of the coordinate frame. Furthermore, as described in Francis (1972), the curvature of the earth introduces a gentle refraction of

**k**that becomes significant when rays propagate horizontally on global scales. The appendix demonstrates how to compute the correct ray-tracing equations for an arbitrary dispersion relation in curvilinear coordinates. The equations resulting from (10) are

where

To facilitate comparison with columnar gravity wave parameterization schemes, the ray tracer can also be run in columnar mode in which

The last two terms in each of (14), (15), (16) were absent in the GROGRAT model of Marks and Eckermann (1995). These terms are necessary, for instance, to ensure that *kr* cos*θ* (and therefore the angular momentum flux) is constant along a ray propagating in a quiescent basic flow. As an example, consider a quiescent, nonrotating, and spherically symmetric atmosphere. Ray solutions in this environment must respect the atmosphere’s spherical symmetry since without rotation there is no difference between the zonal and meridional directions. This means that the horizontal projection of rays follows great circles, and this is true if these terms are included. On the other hand, if these terms are absent, a ray launched eastward in the midlatitudes will propagate along a line of constant latitude, violating spherical symmetry. Of course, these differences only matter if rays travel substantial horizontal distances or if the rays stray close to the poles, as noted in Dunkerton and Butchart (1984).

### b. Numerical implementation

To initialize the model, the steady background fields (*U, V, N*, *α*) are given on a regular (*λ, θ, r*) grid. Rays are initialized by specifying the launch location, horizontal wavenumber, and intrinsic frequency ω̂ > 0. The vertical wavenumber *m* < 0 is then computed from (10) consistent with upward vertical propagation. The local wave amplitude in terms of vertical particle displacements and the cross-sectionally integrated wave action flux *F* are also specified. Though the ray-tracing equations are valid for time-dependent basic flows, the model currently does not allow for them. In this setting, *t* is not a real physical time but merely a parameterization of distance along a ray tube: it measures the amount of time it would take for the steady wave field to develop from undisturbed initial conditions due to forcing at the lower boundary.

The gridded background fields are extended to the ray location by local interpolation in (*λ*, *θ, r*) using the tricubic scheme of Lekien and Marsden (2005). The scheme is local in the sense that it requires only data from the corners of the grid box containing the ray and smooth in the sense that the function values and all first partial derivatives agree on the grid-cell boundaries. The function values and a set of seven partial and mixed partial derivatives must be prescribed at each grid point. These derivatives are computed by second-order finite differences. The grid must be regularly spaced in (*λ, θ, r*) for the interpolation to be smooth, so datasets that use pressure levels and/or Gaussian grids are regridded to meet this requirement.

The ray-tracing equations are integrated using the built-in Matlab function ode45. This ODE solver uses an explicit Runge–Kutta (4, 5) formula, the Dormand–Prince pair, and an adaptive step size algorithm to keep relative error within prescribed tolerances. For the experiments reported in this paper, the relative error tolerance is 10^{−3}. In addition, we found it essential to include redundant equations for the background fields as dependent variables in the ray-tracing system in order for the adaptive time stepping to work consistently. Otherwise, it would not be the case that *dω*/*dt* = 0 for steady background fields, for instance. The GROGRAT model of Marks and Eckermann (1995) accomplishes this by including an equation for *dω*/*dt*, but we preferred to leave conservation of *ω* as an independent check on the validity of the solutions.

Ray integrations are terminated when a ray leaves the model domain or the wave action flux falls below a prescribed tolerance. The wave action flux tolerance used in the following experiments was 1% of the total wave action flux launched, divided by the number of rays.

Some assumption was also required to model the behavior of rays undergoing vertical reflection. As discussed in section 2b, the detailed calculation of the wave field in the neighborhood of a caustic is a subtle and poorly constrained problem. In particular, the amplitude depends on the phase of the incoming and outgoing waves. Some ray tracing codes, such as Marks and Eckermann (1995), terminate rays upon vertical reflection because of a local breakdown of the WKB approximation. They were primarily concerned with the characteristics of waves that entered the stratosphere and not vertically trapped tropospheric waves, so the latter were removed using this WKB criterion. This procedure is equivalent to assuming that the angular momentum flux of the upward traveling wave is perfectly canceled by the flux of downward flux of the reflected wave, which should be accurate in the absence of dissipation when horizontal propagation and refraction are neglected. However, this cancellation is not to be expected in the present case, when **k*** _{H}* evolves along the ray. It is well known (e.g., Lighthill 1978) that ray theory can be “healed” by the consideration of extra terms, which become significant only in the neighborhood of a reflection. The solutions obtained by those methods partially justify the continuation of the ray-tracing equations through reflection points, so vertical reflection does not automatically trigger ray termination in our model.

### c. Nonlinear wave saturation

Wave saturation schemes are heuristic methods by which nonlinear wave breaking can be modeled within a linear ray-tracing computation. Numerous saturation schemes exist in the literature and, to facilitate comparison of our model to the topographic wave parameterization of AGCM3, its saturation scheme has been used here in a modified form that applies throughout the entire frequency range | *f* | < ω̂ < *N*. Full details of this scheme can be found in McFarlane (1987). Basically, it is assumed that max|*u*′_{z}|/*N* ≤ *F*_{c}, where *u*′ is the component of the disturbance horizontal velocity in the direction of the horizontal wavenumber vector **k*** _{H}*, the maximum is over a wave cycle, and the threshold value of the Froude number

*F*

_{c}= ½. This can be converted into a saturation threshold for the vertical particle displacement

*ζ*′ = Re[

*ζ̂*

*ρ*

_{0}/

*ρ*exp(

*i*

*θ*)] by using the standard polarization relations. This yields

for the density-adjusted amplitude *ζ̂*. The vertical action flux density is related to *ζ̂* by

therefore, (18) implies a local saturation threshold on this flux.

However, the model only carries the cross-sectionally integrated flux *F* of (9) along the ray tube, so, as discussed in section 2b, at this step we must make the crude assumption that the cross-sectional area is constant along a ray tube. This means that *F* ∝ *Aw _{g}* and, therefore, (18) implies a threshold directly on

*F*. Now, at each time step along the ray the value of

*ζ̂*required to hold the wave action flux

*F*constant is computed using Eq. (19). The saturation amplitude

*ζ̂*

_{sat}is then computed using Eq. (18) and, if the local amplitude satisfies

*ζ̂*<

*ζ̂*

_{sat}, then no dissipation occurs and

*F*maintains its previous value. However, if the amplitude exceeds the saturation amplitude, then the saturation

*ζ̂*=

*ζ̂*

_{sat}is applied and a new, reduced value of

*F*is computed. It is clear that

*F*can only decrease along the ray. On the other hand, the flux of angular pseudomomentum

*P*=

*Fkr*cos

*θ*may increase or decrease along a ray, depending on the refraction by the basic flow.

We exploit a side effect to allow for vertical wave reflection without saturation, despite the fact that ray theory predicts infinite amplitude when *m* → 0. The point is that both *ζ̂* and *ζ̂*_{sat} diverge when *m* → 0 but the ratio *ζ̂*^{2}/*ζ̂*^{2}_{sat} ∝ |*m*|(*m*^{2} + *α*^{2}) goes to zero. In other words, the numerical scheme continues the ray through the reflection process without activating the saturation scheme.

## 4. Experiments and results

AGCM3 provided the background fields on a T47 Gaussian grid with 31 vertical levels. The model domain extends from the ground to 40 500-m altitude. Because the ray tracer’s interpolation scheme requires a regular spatial grid, the raw GCM data were regridded using two-dimensional cubic spline interpolation in the horizontal on each level, followed by one-dimensional cubic spline interpolation along each vertical column. The *α*^{2} field was computed before regridding using cubic splines along each vertical column to compute the vertical derivative of density.

AGCM3 uses the topographic gravity wave parameterization scheme of Scinocca and Mc-Farlane (2000). The scheme employs two waves at launch to represent anisotropy in the topographic wave field. These waves are launched from each surface grid point over land with zero absolute frequency at each model time step. Properties of this scheme employed for the present study are derived from instantaneous output of the prameterization every six model hours. Included in this output is the launch value of the vertical displacement amplitude *ζ̂*_{0}.

The ray tracer was used for offline computations based on the output of the AGCM3 model. Specifically, we compute the global wave–induced angular momentum transport for a realistic topographic launch spectrum with and without three-dimensional propagation. For topographic waves with *ω* = 0, the intrinsic frequency at the launch level is

For typical topographic wavelengths, ω̂ falls in the intermediate range | *f* | ≪ ω̂ ≪ *N*, for which it is expected that three-dimensional refraction effects have a weak impact. This exercise is intended to demonstrate that, although the assumption of horizontal homogeneity in the parameterization of topographic waves is incorrect in principle, it has a weak impact on the global wave-induced vertical transport of angular momentum.

It is also interesting to test the dependence of three-dimensional refraction effects on intrinsic frequency. From the viewpoint of testing their impact on angular momentum transport, it would be ideal to conduct a similar comparison of a physically realistic nonorographic parameterization with fluxes computed from three-dimensional ray tracing solutions. The authors hope to undertake such a project in the future. However, AGCM3 does not include a nonorographic parameterization scheme.

Thus, as a simple test of the frequency dependence of angular momentum transport, an artificial, ad hoc, low frequency spectrum was created by remapping the intrinsic frequencies of the AGCM3 launch spectrum from the interval | *f* | < ω̂ < *N* onto the interval | *f* | < ω̂ < 3| *f* |, while holding *k*, *l*, *ζ̂*, *F* and the launch Froude number constant. Specifically, we use the power-law mapping

While this remapping is unphysical, as a thought experiment it has the advantage of redistributing the frequency content of the topographic spectrum while holding constant all other properties that could potentially impact angular momentum transport (i.e., horizontal wavenumber, the quantity and spatial distribution of vertical wave action flux, and launch stability).

Another important difference between the topographic launch spectrum and the remapped spectrum is that the topographic waves satisfy Eq. (20) so that *ω* = 0, whereas the remapped waves do not. We note in passing that a third computation was performed with a second spectrum in which the launch values of *k* and *l* were rescaled to maintain *ω* = 0. The results obtained were qualitatively similar to those of the first remapping, so this indicates a certain robustness of the effects we describe here.

### a. Ray examples

Examining individual ray solutions provides a qualitative illustration of the dependence of three-dimensional propagation effects on intrinsic frequency. As mentioned above, previous ray tracing studies, such as Dunkerton (1984), have described similar solutions, but we include these here for purposes of illustration.

Figures 1 and 2 show the detailed evolution of two waves with altitude in the topographic and low frequency cases. Both figures show the vertical evolution of *k,* the zonal velocity *U* at the location of the ray, and the cross-sectionally integrated vertical angular pseudomomentum flux *P* from (8) for both three-dimensional and columnar propagation. The topographic case is shown in Fig. 1. Here, the relative change in *k* is small over the depth of the atmosphere. Furthermore, because the rays do not propagate far horizontally, the zonal velocity profile experienced by the three-dimensional ray is almost identical to that experienced by the columnar ray. Not surprisingly, both topographic rays have avoided saturation in both cases.

The profiles for the low frequency rays launched from the same site, shown in Fig. 2, show a large impact from three-dimensional refraction. There is a large relative change in *k* over the depth of the atmosphere. Because the three-dimensional rays propagate significant horizontal distances, the zonal velocity profile that they experience differs significantly from the columnar case. The changed environment leads to the deposition of wave momentum flux at a significantly different altitude.

### b. Comparison with GCM parameterization scheme

Though individual solutions of the ray-tracing equations show that three-dimensional propagation can lead to significantly different outcomes for individual ray tubes compared with columnar propagation, by themselves they give no information about how the new effects impact global, aggregate forcing of mean winds by internal gravity waves.

To investigate this question, the ray tracer was used to compute the trajectory of every wave in the AGCM3 parameterization launch spectrum. Though the waves are launched from just above ground level in the GCM, they are launched from the GCM model level nearest to 4000 m in these computations. Thus, the experiments with 3D propagation are equivalent to switching on 3D effects once waves have propagated above 4000 m. The GCM parameterization includes detailed consideration of near-ground dynamics such as blocking and downslope windstorm flow. Launching from the middle of the troposphere avoids consideration of these boundary layer issues. The ray tracing solutions of all launched rays are added to construct vertical profiles of globally integrated angular pseudomomentum flux, which is equivalent to the flux of mean angular momentum.

### c. Results for topographic waves

Figure 3 shows the January mean globally integrated vertical angular pseudomomentum flux profiles computed by the ray tracer with and without three-dimensional refraction. The daily variability of the profile at three reference levels is indicated by error bars spanning an interval containing 90% of the variability from time step to time step. To validate the ray-tracing scheme, the same flux profile computed from the output of the AGCM3 parameterization scheme is also shown. The AGCM3 scheme is equivalent to the ray-tracing model in columnar propagation mode, except for differences in the representation of the background fields arising from the use of different interpolation methods and the modification of the saturation scheme to handle rotation and nonhydrostatic waves.

It is evident that the introduction of three-dimensional refraction does not have a strong impact on the global vertical transport of angular pseudomomentum by topographic gravity waves. Both profiles reach a minimum flux of approximately −2.0 × 10^{19} kg m^{2} s^{−2} near the tropopause, rapidly increasing to approximately −3.2 × 10^{18} kg m^{2} s^{−2} at 20 000 m, then gently increasing to approximately −9 × 10^{17} kg m^{2} s^{−2} at the top of the model domain. This profile shape corresponds to a strong retrograde forcing of the lower stratosphere, relatively little forcing of the middle stratosphere, and approximately 5% of the launched retrograde angular pseudomomentum flux reaching the upper stratosphere or mesosphere. The relative difference between the two profiles averages 4% through most of the domain and increases to 13% at 40 000 m. The angular pseudomomentum flux escaping the top of the model domain is more negative for three-dimensional propagation than for columnar propagation, corresponding to a slightly stronger retrograde forcing of the upper atmosphere. A globally integrated vertical angular pseudomomentum flux *O*(2 × 10^{19} kg m^{2} s^{−2}) is equivalent to an average vertical pseudomomentum flux density of 8 × 10^{−3} kg m^{−1} s^{−2}.

The qualitative shape of the vertical flux profile can be explained by considering the mean vertical profiles of zonal velocity and density, shown in Fig. 4b. The launch spectrum is divided between waves carrying prograde angular pseudomomentum flux, for which *k* < 0, and those carrying retrograde angular momentum flux, for which *k* < 0. Figure 4a shows how the vertical flux profile is decomposed into a contribution from retrograde waves and one from prograde waves.

Broadly speaking, prograde waves are generated by winds with a westerly component near the ground, and prograde waves are generated by winds with an easterly component. Two processes combine to dissipate the majority of both prograde and retrograde waves below 40 000 m. First, rays can be dissipated at critical surfaces on which ω̂ = −**U** · **k** → | *f* |. Second, the oscillation amplitude increases like *ρ*^{−1/2} to conserve action flux as density decays. Whether the ray encounters a critical surface or not, increasing amplitudes eventually lead to saturation of the wave field by nonlinear instabilities.

As the waves propagate upward, prograde waves are typically annihilated at critical surfaces as the zonal wind shifts from easterly to westerly leading into the jet. In Fig. 4, one can see that the contribution from prograde waves is almost completely suppressed by 13 000 m, where the westerly jet reaches its maximum. On the other hand, the positive shear environment suppresses breaking for retrograde waves, and there is only weak decay of their contribution over the same interval. Above the jet, the negative shear environment begins to Doppler shift retrograde waves back toward lower intrinsic phase velocities, making them more susceptible to saturation as the atmospheric density continues to decay. Between 15 000 and 20 000 m, the majority of the pseudomomentum flux in retrograde waves is dissipated by saturation processes. The remaining flux is associated with waves with amplitudes sufficiently low to escape significant saturation between 20 000 and 40 000 m.

### d. Results for low-frequency waves

As expected, the vertical angular pseudomomentum flux profiles computed for the low-frequency spectrum show a much greater impact from three-dimensional propagation and refraction. Figure 5 shows the angular pseudomomentum flux profiles for the same month computed with the hypothetical low-frequency launch spectrum. The profile for three-dimensional rays becomes more negative than the columnar profile immediately above the launch level, and retains this bias through the entire depth of the model domain. The three-dimensional profile is 41% more negative on average and transmits 81% more retrograde angular pseudomomentum flux to the upper atmosphere. The daily variance of the flux profile is increased in the three-dimensional case by an average of 30%. The increase in variance is primarily due to increased variability of the angular momentum flux above 20 000 m.

A heuristic argument to explain why three-dimensional refraction should lead to increased negative angular pseudomomentum flux can be based on the analogy between *k* and the gradient of a passive tracer discussed below (5). Haynes and Anglade (1997) demonstrated that passive advection in typical stratospheric flows has a robust statistical tendency to increase the magnitude of the tracer gradient in a fashion compatible with the usual results from random advection theory. As this result did not rely on the physics underlying the advecting velocity field, and this is the only sense in which wave vector refraction differs from tracer gradient evolution, it seems likely that this cumulative gradient stretching is inherited by *k* in the ray-tracing problem. Though individual rays may experience a decrease in |*k*|, the statistical tendency for |*k*| to increase is apparent. Indeed, a similar cascade to shorter horizontal wavelengths has been observed in a study of internal waves in the ocean propagating through a background drawn from the Garrett–Munk spectrum (see Flatté et al. 1985; Henyey et al. 1986). Thus, as prograde waves are filtered out at low altitudes, it can be expected that retrograde waves will have more negative *k* on average when three-dimensional refraction is allowed.

Figures 6 and 7 indicate that increased retrograde flux *P* = *kFr* cos*θ* in the three-dimensional case is due to the tendency for *k* to become more negative combined with a tendency for three-dimensional rays to transmit more wave action flux. In Fig. 6a, the value of *k* for the three-dimensional ray at the top of the model domain has been plotted against its launch value *k*_{0} for all rays that escape to the upper atmosphere in both three-dimensional and columnar mode. The median relative change (*k* − *k*_{0})/*k*_{0} is −17%. This tendency is more pronounced when we consider only the strongest rays. Figure 6b shows rays with the largest action flux at the top of the model domain, which explain 90% of the flux transmitted to the upper atmosphere. For these rays, the median relative change of *k* is −32%.

Meridional propagation contributes very little to the change in angular pseudomomentum flux. The resulting change in *r* cos *θ* has a mean of less than 5%.

Figure 7 compares the wave action flux at the top of the model domain for individual rays with columnar and three-dimensional propagation. The flux increases 52% with three-dimensional propagation, which together with the mean decrease of *k* contributes to the increase of retrograde momentum flux to the upper atmosphere. Figures 6 and 7 demonstrate that the tendencies just described are statistical and that there is great variation in the evolution of these quantities between individual rays.

## 5. Concluding comments

Consistent with the scaling argument presented in section 2, we have found that significant changes in angular momentum fluxes due to three-dimensional refraction and horizontal propagation are more or less confined to low frequency waves. This result indicates that columnar parameterization schemes for topographic waves should produce results that are little affected by the neglect of three-dimensional effects, which is encouraging. There is, however, the caveat that we used data at the fairly low horizontal resolution of T47. It remains to be seen what the impact of increased horizontal resolution would be and, also, how three-dimensional effects would change the parameterization of nontopographic waves, especially in the low frequency range. We hope to follow up on these questions in the near future.

## Acknowledgments

The authors thank Stephen Eckermann for providing the source code to GROGRAT, which was an extremely useful resource. AH’s work is supported by a National Science Foundation Graduate Research Fellowship and a National Defense Science and Engineering Graduate Fellowship. OB gratefully acknowledges financial support for this work under the National Science Foundation Grants OCE-0324934 and DMS-0604519. We are grateful for the comments of two anonymous referees that significantly improved an earlier version of this paper.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**(

**,**(

**,**

**,**

**,**(

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**(

### APPENDIX

#### Ray Tracing on a Sphere

If the wavelengths of internal gravity waves are such that the curvature of the earth is negligible on the scale of a wavelength, then the dispersion relation *ω* = Ω(*k*, *x*) computed in Cartesian coordinates may be retained, simply substituting *k* = **k** · ** λ̂**,

*l*=

**k**·

*and*

**θ̂***m*=

**k**·

**r̂**. The components of

**k**in the underlying Cartesian coordinates will be denoted in index notation

*k*. In general, the dispersion relation may be written in Cartesian coordinates

_{i}As the unit vectors of the spherical coordinate system are slowly varying on the scale of a wavelength, one can use this dispersion relationship to obtain ray-tracing equations. This approach was used in Francis (1972) and further justified in Dong and Yeh (1993), and we use it now to demonstrate how the use of spherical coordinates introduces metric terms to the refraction of the wavenumber for a general dispersion relationship of the form in Eq. (A1).

We then have

Therefore,

and likewise *c _{gλ}* = Ω

*,*

_{k}*c*= Ω

_{gθ}*.*

_{l}However, simply taking the partial derivative of the dispersion relationship with respect to the spatial variables will not yield the correct refraction equations. We find

By the chain rule,

This equation combined with the expression for **c*** _{g}* and Eq. (A3) implies that

An equation of this form would be obtained for any curvilinear coordinate system. The metric terms in parentheses are a property of the coordinate system. For spherical coordinates,

In this notation, the vector to the left contracts with the gradient. Simplifying, we find

Similar computations give

## Footnotes

*Corresponding author address:* Alexander Hasha, Courant Institute of Mathematical Sciences, Room 1108, 251 Mercer Street, New York, NY 10012. Email: hasha@cims.nyu.edu