## Abstract

This work examines the influence of horizontal propagation of three-dimensional (3D) mountain waves on the wave momentum flux (WMF) within finite domains (e.g., the grid cell of general circulation models). Under the Wentzel–Kramers–Brillouin (WKB) approximation, analytical solutions are derived for hydrostatic nonrotating mountain waves using the Gaussian beam approximation (GBA), which incorporates both the wind vertical curvature effect and the height variation of stratification. The GBA solutions are validated against numerical simulations conducted using the Advanced Regional Prediction System (ARPS). In the situation of idealized terrain, wind, and stratification, the WMF obtained from the GBA shows a good agreement with the numerical simulation. The effect of wind curvature in enhancing the WMF is captured, although the WKB-based GBA solution tends to overestimate the WMF, especially at small Richardson numbers of order unity. For realistic terrain and/or atmospheric conditions, there are some biases between the WKB GBA and simulated WMFs, arising from, for example, the missing physics of wave reflection. Nonetheless, the decreasing trend of finite-domain WMF with height, because of the horizontal propagation of 3D mountain waves, can be represented fairly well. Using the GBA, a new scheme is proposed to parameterize the orographic gravity wave drag (OGWD) in numerical models. Comparison with the traditional OGWD parameterization scheme reveals that the GBA-based scheme tends to produce OGWD at higher altitudes, as the horizontal propagation of mountain waves can reduce the wave amplitude and thus inhibit wave breaking.

## 1. Introduction

Orographic gravity waves (OGWs), or mountain waves, are generated as stably stratified airflow goes over mountains. The breaking of OGWs can deposit the wave momentum into the mean flow, resulting in the well-known orographic gravity wave drag (OGWD). The deposition of the wave momentum plays an important role in the momentum budget of atmospheric general circulation (Holton 1982). Because OGWs are, in general, subgrid-scale physical processes, the OGWD needs to be parameterized in general circulation models (GCMs; Kim et al. 2003). Indeed, the parameterization of OGWD is now a common feature of contemporary numerical models, which can help alleviate the systematic wind and temperature biases in the model (e.g., Palmer et al. 1986; McFarlane 1987).

Current OGWD parameterization exclusively makes use of the single-column assumption: namely, OGWs propagate vertically within the model grid cell where the subgrid-scale orography (SSO) is present (Kim and Arakawa 1995; Lott and Miller 1997; Gregory et al. 1998; Scinocca and McFarlane 2000; Webster et al. 2003; Kim and Doyle 2005). This assumption is appropriate in the case of two-dimensional (2D), hydrostatic, and nonrotating waves that have zero horizontal group velocity. In contrast, three-dimensional (3D) OGWs can experience significant lateral propagation, thus violating the single-column assumption. For instance, Smith (1980) and Shutts (1998) studied the structure of 3D mountain waves generated in environments with zero and directional wind shear over a circular bell-shaped mountain, respectively. (Directional wind shear means that the ambient wind direction changes with height.) In both cases, the wave fields radiate away from the mountain, forming a bow-shaped pattern. In the directional wind case, the waves are advected downstream for an infinitely long distance on meeting their selective critical levels (Broad 1995; Shutts 1995).

The horizontal propagation of 3D mountain waves has not yet been considered in existing OGWD parameterization schemes. Recently, Eckermann et al. (2015a,b, hereafter E15a,b) for the first time examined the effect of “horizontal geometrical spreading.” Using the Fourier-ray method (Broutman et al. 2002), E15a addressed 3D mountain waves generated in unidirectional flow over elliptical obstacles. The peak wave amplitude was shown to decrease monotonically with height. The accumulated reduction of wave amplitude was comparable to, or even larger than, that caused by the refraction of vertical wavenumber. E15b simplified the Fourier-ray solution and derived two asymptotic solutions for the vertical displacement and steepness that were relevant to the parameterization of OGWD. The derived analytical solutions, which depend only on the terrain’s elliptical aspect ratio and a normalized height coordinate involving wind and stability variations, agree very well with the exact solution.

E15a and E15b mainly focused on the reduction of wave amplitude with height due to horizontal geometrical spreading, which leads to a higher wave breaking level as compared to the 2D wave case. E15b also noted that 3D OGWs could propagate out of the domain where they were generated (e.g., a GCM grid cell), especially at high altitudes. Since the waves outside the domain of interest (i.e., grid cell) do not exert a direct influence on the *local* mean flow, their momentum should be discarded from the budget of domain of interest. Therefore, the lateral propagation of 3D waves has twofold influence. For one thing, it reduces the wave amplitude with height and consequently increases the altitude of wave breaking. For another, it decreases the amount of wave momentum flux (WMF) deposited into the *local* mean flow. This second effect, to the authors’ knowledge, has not yet been quantitatively studied in detail. In principle, quantifying the percentage of 3D waves remaining within a finite domain is nontrivial. It changes with the domain size, mountain shape, and location relative to the domain, environmental wind, and stratification, among other factors. The only way is to explicitly compute the WMF within the domain. However, analytical solutions only exist in very limited cases: for example, constant wind over elliptical terrain (Phillips 1984). Asymptotic solutions that can accurately represent the OGWs are thus needed.

Using the ray-tracing method, Shutts (1998) derived a far-field stationary-phase solution for hydrostatic, nonrotating mountain waves in a turning wind. This solution breaks down in the region directly over the terrain where multiple rays intersect [i.e., caustics (Lighthill 1978)], however. To improve the ray solution, Broutman et al. (2002) suggested using Maslov’s method that expressed the ray solution in the spectral domain rather than in the spatial domain. Although the Fourier-ray solution can avoid caustics in the spatial domain, it may suffer from caustics in the spectral domain. Pulido and Rodas (2011, hereafter PR11) proposed another approach: that is, the Gaussian beam approximation (GBA), which is a higher-order ray approximation. For each ray, the GBA considers not only the contribution of this ray but also a beam of rays around it to the wave field. Thus, the ray solution is well defined even at caustics. In addition, the GBA can be readily applied to mountain waves forced by complex topography via the superposition of multiple Gaussians.

The performance of GBA is examined for a few idealized wind profiles in PR11, which gives encouraging results. However, the effect of wind vertical curvature is not considered. Teixeira et al. (2004, hereafter TMV04) developed an analytical model of OGWs using the second-order Wentzel–Kramers–Brillouin (WKB) approximation. It shows that the surface WMF increases as the Richardson number decreases for the wind turning with height while maintaining its magnitude. This drag enhancement due to the wind vertical curvature appears to play an important role in affecting the global atmospheric torque (Miranda et al. 2009). Whether the GBA can capture the effect of wind vertical curvature will be studied in this work. The behavior of GBA under the circumstance of variable buoyancy frequency will also be examined, while only constant buoyancy frequency is considered in PR11.

The main goal of this work is to establish a new OGWD parameterization scheme using the GBA method, which is able to account for the horizontal propagation and the horizontal geometrical spreading of mountain waves. The rest of this paper is organized as follows. In section 2, analytical solutions are derived for hydrostatic nonrotating mountain waves using the GBA. Section 3 describes the terrain, wind, and stratification, as well as the setup of numerical experiments. Section 4 validates the GBA solutions against the numerical simulation for a number of cases, including both idealized and realistic terrain and atmospheric conditions. In section 5, a new OGWD parameterization scheme is proposed using the GBA and compared to the traditional scheme. Section 6 summarizes the paper with discussions.

## 2. Theory

For hydrostatic mountain waves generated in adiabatic, inviscid, and Boussinesq flow, with Earth’s rotation neglected, the vertical wavenumber is given by

Here, *N* is the Brunt–Väisälä frequency (or buoyancy frequency), *K* is the magnitude of horizontal wavenumber **K** = (*k*, *l*), , **V**(*z*) = [*U*(*z*), *V*(*z*)] is the horizontally uniform mean wind, and . (The subscript *z* denotes partial derivative with respect to *z* unless otherwise stated.) Unlike PR11, the wind curvature effect **V**_{zz} is considered. Note that *m* is singular at [i.e., **V**(*z*) = 0 or **V**(*z*)⊥**K**]. In the former case, it is the well-known total critical level where all wave components are attenuated (Booker and Bretherton 1967). The latter is termed the selective critical level, because only the waves normal to the mean flow are absorbed (Broad 1995; Shutts 1995). Readers are referred to Xu et al. (2012, hereafter XWX12) for more details about the selective critical-level filtering.

Under the WKB approximation, the vertical velocity of mountain waves can be solved as [cf. Eq. (25) of TMV04]

where is the Fourier transform of vertical velocity *w*(*x*, *y*, *z*). Linear waves satisfy the free-slip lower boundary condition: that is,

with being the Fourier transform of terrain height *h*(*x*, *y*) and = **V**_{0} · **K**. On substituting Eq. (3) into Eq. (2), one can readily obtain

By definition of *w* = *dη*/*dt*, the vertical displacement *η* in the spectral space is

The corresponding vertical displacement in the physical domain can be recovered by virtue of inverse Fourier transform: that is,

where is the wave phase; is the Heaviside function equal to zero (unity) for the negative (positive) argument, which is to account for the selective critical level filtering by directional wind shear (XWX12). In the absence of wind curvature, the above equation is reduced to

The height variation of buoyancy frequency is involved, consistent with Eq. (26) of Broutman et al. (2002) and Eq. (13) of E15a. In the situation of constant buoyancy frequency, Eq. (7) can be simplified to Eq. (4) of PR11: that is,

Under the GBA, the amplitude of ray **K**_{c} = (*k*_{c}, *l*_{c}) is assumed to be represented locally as a Gaussian-shaped envelope centered at **K**_{c}:

where is the Gaussian amplitude determined by the terrain spectrum; *σ*_{k} and *σ*_{l} are the Gaussian width in *k* and *l* directions. Hereinafter, the subscript *c* indicates evaluation at the central wavenumber **K**_{c}. Using Taylor’s expansion, the phase function can be expressed as

with higher-order terms omitted. On substituting the above two equations into Eq. (6), the integral can be solved exactly [cf. Eq. (21) in PR11]:

where the matrices and are given by

and the determinant of is

It is the *real* part of Eq. (11) that denotes the vertical displacement in the physical domain, whereas its *magnitude* represents the wave amplitude, similar to Eq. (28) in E15a. Equation (11) is the vertical displacement arising from the wave centered at **K**_{c}. The total wave field forced by an isolated mountain is often composed of a number of waves, which requires superposition of multiple Gaussians. An algorithm conserving the wave energy was proposed in PR11 (see their Fig. 7), with the amplitude of each Gaussian given by

where Δ*k* and Δ*l* are the spectral resolutions, and (*σ*_{k}, *σ*_{l}) = *G*_{c}(Δ*k*, Δ*l*) with *G*_{c} being a tunable parameter.

The GBA solutions for other wave quantities [e.g., horizontal velocity **v**_{h} = (*u*, *υ*)] can be derived from the polarization relation of plane waves. Once these variables are obtained, the WMF within a finite domain is [cf. Eq. (66) in E15b]

where *M* and *N* are the number of grid points within the domain, Δ*x* and Δ*y* are the grid spacings, and *ρ*_{0} is the reference density of Boussinesq flow.

## 3. Methodology

### a. Terrain

Two mountains are examined in this study (see Table 1). The first is the idealized 3D circular bell-shaped mountain (BELL; Fig. 1a)

where *h*_{m} is the mountain amplitude and *r*_{a} is the half-width. For linear and hydrostatic waves, we choose *h*_{m} = 10 m and *r*_{a} = 50 km, respectively. The second is the Hainan Island (HNI; Fig. 1b) located in the South China Sea in the area 18°–20°N, 109°–111°E. This realistic mountain is about 250 km long and 200 km wide, consisting of several isolated hills of up to 1800 m in height. Since we are only interested in linear mountain waves, the Hainan Island height is scaled to 10 m.

The terrain spectrum, which is required by the GBA, is obtained by virtue of the fast Fourier transform (FFT) software FFTPACK (Swarztrauber 1982). Given a finite domain of *M* × *N* points, the spectral resolutions are Δ*k* = 2*π*/(*M*Δ*x*) and Δ*l* = 2*π*/(*N*Δ*y*), respectively. The FFT can give a number of discrete wavenumbers that are integer multiples of Δ*k* and Δ*l*. However, the GBA solutions utilize waves at half grid points in the spectral space: that is, *k*_{c} = (0.5 + *ii*)Δ*k* (with *ii* = 0, 1, 2, …), and *l*_{c} = (0.5 + *jj*)Δ*l* (with *jj* = 0, 1, 2, …), to avoid the zero wavenumber. The Fourier spectra of these half-grid wavenumbers are calculated using Eq. (14).

### b. Environment wind and buoyancy frequency

To exclude the influence of selective critical-level absorption, only unidirectional wind profiles are considered. The first wind profile (wCST) is constant, with *U* = 10 m s^{−1} and *V* = 0. The second wind profile (wCUR) is defined as

with *z*_{c} being the total critical level at which the wind vanishes. This wind profile was used in Teixeira and Miranda (2004) to study the effect of wind curvature. It features a nonuniform vertical shear of . The Richardson number at the total critical level is thus , which has an important influence on the WMF. The last one is the monthly mean wind speed at Hainan Island in July 2014 (wHN7; Fig. 2a), derived from the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim; Dee et al. 2011).

We also make use of two different potential temperature profiles. In the first case (bCST), the atmosphere is of constant *N* = 0.01 s^{−1}, with surface potential temperature *θ*_{sfc} = 288 K and surface pressure *P*_{sfc} = 1000 hPa. The second (bHN7) is derived from ERA-Interim in the same way as the Hainan Island wind profile. Seen from Fig. 2b, the buoyancy frequency is relatively low in the troposphere below about *z* = 12 km, and it increases rapidly with height in the lower stratosphere, almost doubling its value at *z* =15 km.

### c. Setup of numerical experiments

The GBA solution is validated against numerical simulations using the Advanced Regional Prediction System (ARPS): that is, a fully compressible, nonhydrostatic model developed at the Center for Analysis and Prediction of Storms (CAPS), the University of Oklahoma (Xue et al. 2000). This numerical model has been well tested in Xue et al. (2000) and Xu et al. (2013) for both linear and nonlinear mountain waves.

A series of numerical experiments is conducted with different configurations of terrain and base-state wind and buoyancy frequency (see Table 2). The model domain has 203 × 203 points in the horizontal with a horizontal grid spacing of 5 km and a vertical grid spacing of 100 m. The mountain is placed at the domain center. A free-slip boundary condition is applied at the bottom, along with a radiation condition at the four side boundaries. For the upper boundary condition, a Rayleigh damping layer is used in the upper portion of the model. A fourth-order horizontal computational mixing is adopted to damp small-scale numerical noises. The Coriolis parameter is set to zero for irrotational gravity waves, and all other physical parameterizations are turned off. The model is integrated for 50 h such that the simulated waves reach a steady state.

## 4. Comparison of GBA with numerical simulations

The GBA solutions are obtained by superposition of 1600 Gaussians (unless otherwise stated), with *k*_{min} = 0.5Δ*k* and *l*_{min} = 0.5Δ*l.* The tunable parameter is set to *G _{c}* = 0.55, a little larger than suggested in PR11. To make the non-Boussinesq model result comparable with the GBA solution, the simulated waves are multiplied by a factor of [

*ρ*(0)/

*ρ*(

*z*)]

^{−1/2}to eliminate the height variation of base-state density.

### a. Constant wind

In the case of constant wind past a circular bell-shaped mountain (ExpBCST), there are exact solutions [cf. Eq. (20) in Smith (1980)]. Figure 3 shows the perturbed vertical velocity at *z* = 4 km (left) and *z* = 8 km (right) from the exact solution (Figs. 3a,b), the ARPS model (Figs. 3c,d), and the GBA solution (Figs. 3e,f), respectively. Consistent with Smith (1980), the wave fields exhibit a bow-shaped pattern that widens with height, propagating laterally and downstream for a few hundred kilometers. Both the model (Figs. 3c,d) and GBA wave fields (Figs. 3e,f) agree well with the exact solution (Figs. 3a,b).

Of particular interest is the WMF, which is most relevant in the parameterization of OGWD. From Fig. 4, the WMF derived from the GBA is in good agreement with the numerical simulation, despite both being a little smaller than the exact solution. The WMF shows a decreasing trend with height. At *z* = 10 km, it is only about 80% of that at *z* = 1 km. This is quite different from the WMF in an infinite domain, which is constant with height. The height decrease of WMF is attributed to the lateral propagation of waves out of the domain. It is thus more notable for a smaller domain. For example, at *z* = 10 km, the WMF within the small box (Figs. 3a,b) is decreased to less than 50% of that at *z* = 1 km (Fig. 4).

### b. Wind profile with vertical curvature

Since exact solutions do not exist in this case, we only compare the GBA solution to the numerical simulation. Two wind profiles wCUR1 and wCUR2 are considered, which have the same surface wind of *U*_{0} = 20 m s^{−1} but different critical-level heights of *z*_{c} = 8000 and 4000 m (Table 1). The corresponding Richardson numbers at the critical level are Ri_{c} = 4 and 1, respectively.

The general patterns of the GBA wave fields are very similar to their numerical counterparts (not shown). Figure 5 presents the vertical distribution of the GBA and simulated WMFs for the two wind profiles. Note that the simulated WMF decreases to zero at the critical level where the mountain waves become nonlinear and finally break. However, the GBA WMF is derived from linear wave theory such that the two WMFs are not comparable near the critical line. As a result, we will only focus on the WMF below. In fact, the waves from the GBA also break at the critical level, as will be shown in section 5. In the case of Ri_{c} = 4, the GBA-based WMF shows broad agreement with the model prediction, which is only about 2.9% larger than its numerical analog (Fig. 5a). In contrast, the GBA overestimates the WMF by about 10.4% when the Richardson number is decreased to Ri_{c} = 1 (Fig. 5b). Nonetheless, the decreasing trend of the WMF with height, which is due to the lateral spread of 3D mountain waves, is well captured by the GBA solution in both cases.

The relatively large discrepancy between the GBA and numerical model at small Richardson numbers is because the WKB solution [Eq. (2)] only has first-order accuracy. As studied in TMV04, extending the WKB solution to second order can better correct the leading-order WMF caused by the wind vertical shear and curvature. According to TMV04’s second-order WKB solution [see their Eq. (50)], the surface WMF for the wind profile [Eq. (17)] is

where is the leading-order WMF. The increase of WMF is inversely proportional to the critical level Richardson number. Figure 6 displays the normalized WMF as a function of Ri_{c}. The enhancement of WMF predicted by the GBA is about twice the second-order WKB solution. This is consistent with the behavior of the first-order asymptotic solution found in TMV04 (see their Fig. 7). The difference between the GBA and the second-order WKB solutions is smaller than 2% for about Ri > 5. This explains why in PR11 the GBA solution agrees well with the “exact solution” (obtained by solving the wave equation using a fourth-order Runge–Kutta algorithm in the wavenumber space), even though the wind curvature is not considered. The Richardson number in their study is Ri = 10.

### c. Realistic wind and buoyancy frequency

Figure 7a shows the vertical distribution of the WMF generated in the realistic wind and buoyancy frequency over a circular bell-shaped mountain (ExpBHN7a). There is good consistency between the GBA-derived and the model-simulated WMFs at altitudes above 8 km. Yet the two WMF profiles diverge significantly in the mid- to lower troposphere. The GBA solution tends to underestimate the WMF by up to 18.6% at *z* = 1 km. As mentioned above, the wind vertical curvature can help enhance the WMF. However, the Richardson number of the realistic wind is basically greater than 10 (not shown), which would result in little WMF enhancement. Thus, it is seemingly ascribed to the realistic buoyancy frequency used. Seen from Fig. 2b, the buoyancy frequency varies radically in the upper troposphere and lower stratosphere (near *z* = 13 km). Wave reflection can take place at the interface of layers of different stratifications (Teixeira et al. 2005; Jiang et al. 2014). Interaction between the incident and reflected waves can amplify the WMF under favorable conditions (e.g., Klemp and Lilly 1975). However, these downward reflected waves are not taken into account in the GBA solution, which is based on the WKB approximation.

The influence of the realistic buoyancy frequency on the WMF can be verified by a sensitivity experiment ExpBHN7b, which is the same as ExpBHN7a but uses the constant wind wCST. In this case, the numerical model also produces a greater WMF than the GBA, especially in the troposphere (Fig. 7b). The largest difference between the model and GBA WMFs is ~13.3% at *z* = 1 km. Note that the difference between the model and GBA WMFs is reduced above *z* = 13 km. It is likely because parts of the waves are reflected downward in the model. If there were a significant amount of wave reflection, the simulated WMF may become smaller than the WKB-GBA WMF above the reflection level.

### d. Realistic terrain

For mountain waves generated in constant wind over the realistic Hainan Island terrain (ExpHCST), Figs. 8a and 8b display the perturbed vertical velocity at *z* = 5 km from the model and GBA, respectively. While the GBA can capture the overall wave pattern, the detailed structure is not as coherent as in the idealized mountain case (Fig. 3). For instance, the model produces small-scale wave activities near the mountain peak in western Hainan Island, which are absent in the GBA wave field. To better capture the finescale structure (i.e., high-wavenumber part of the spectrum), more Gaussians are needed (see Fig. 8c, which uses 3600 Gaussians), and vice versa (Fig. 8d, which uses 400 Gaussians).

Figure 9 shows the vertical distribution of the WMF. For the *x-*component WMF, similar height-decreasing trends are found in both the model and GBA. But the GBA WMF (using 1600 Gaussians) is smaller than its numerical counterpart by about 9% at *z* = 1 km. This difference is mainly due to the lack of small-scale waves in the GBA wave field. Increasing (decreasing) the number of Gaussians results in an increase (decrease) of the WMF. Moreover, different from the WMF produced by the circular bell-shaped mountain, the *y*-component WMF is nonzero, which is attributed to the anisotropy of Hainan Island. As noticed in Phillips (1984), a lateral force will be created when the upstream wind is not parallel to the terrain symmetry axis.

## 5. Parameterization of OGWD using the GBA

For the parameterization of OGWD, there are in general two issues that need to be addressed: that is, at which altitude the waves break, and how much WMF is absorbed by the mean flow.

For the first one, mountain waves are assumed to break at the height where the steepness of vertical displacement (i.e., ) exceeds a threshold *α*, which equals to unity for the simplest case of convective overturning (Fritts 1984). On the other hand, Palmer et al. (1986) suggested the use of “minimum Richardson number” in their parameterization scheme (i.e., the smallest Richardson number under the influence of gravity waves):

where is the mean-flow Richardson number. Wave breaking takes place when falls below a critical value Ri_{c}, which is often set to 0.25, corresponding to the onset of Kelvin–Helmholtz instability. Evidently, this wave-modulated Richardson number is closely related to the wave amplitude. Figure 10 shows the wave amplitude in the cases of wCST and wCUR1 over a circular bell-shaped mountain, obtained from the GBA solution. (Note that the growth of wave amplitude due to decreasing density is not considered.) In the constant wind case, the maximum wave amplitude is located near the mountaintop and decreases with height, owing to the horizontal geometrical spreading (Figs. 10a,b). The peak wave amplitude at *z* = 5 km accounts for about 33% of that at the surface, in quantitative agreement with the Fourier-ray solution of E15a (see their Fig. 1). In the latter case, the peak wave amplitude is also concentrated near the mountaintop (Fig. 10c). Moving aloft, the wave amplitude first decreases with height as a result of horizontal geometrical spreading. On approaching the critical level at *z* = 8 km, the wave amplitude increases rapidly (Fig. 10d), which is attributed to the accumulation of wave action (Shutts 1998). The increasing wave amplitude leads to a wave-modulated Richardson number less than 0.25 at about *z* = 7.95 km (Fig. 11), implying wave breaking there.

For the second issue, the wave saturation hypothesis (Lindzen 1981) is often adopted, which reduces the amplitude of breaking waves to a saturation wave amplitude *η*_{sat}. In current OGWD parameterization schemes, the WMF deposited into the mean flow is equal to (see page 2333 of E15a)

where *τ*_{ref} is the WMF at the reference level *z*_{ref}, and *η*_{a} is the peak wave amplitude. However, as shown above, the WMF within a finite domain decreases with height due to horizontal propagation. Therefore, the amount of WMF absorbed by the local mean flow at the wave breaking level *z*_{b} should be

### a. New OGWD parameterization scheme

A simple OGWD parameterization scheme is proposed for linear OGWs, with the wave horizontal propagation considered. For nonlinear low-level wave breaking and flow blocking, they can be treated similarly as in previous schemes (e.g., Kim and Arakawa 1995; Lott and Miller 1997; Webster et al. 2003). The vertical profile of the WMF is calculated as follows:

Based on a predefined critical Froude number, a reference level

*z*_{ref}is determined according to the low-level wind and stratification and SSO elevation [see Eq. (4) in Lott and Miller (1997)]. The SSO*above z*_{ref}is Fourier transformed to obtain the terrain spectrum.The

*U*and*V*winds at model levels are fitted by cubic spline interpolation (or other methods) to gain a smooth wind profile, which is of a continuous second derivative used by the vertical wavenumber.Compute the maximum wave amplitude

*η*_{a}and the wave-modulated Richardson number from the first model level above*z*_{ref}. Once wave breaking is activated at a level*z*_{b}where is smaller than Ri_{c}, compute the finite-domain WMF(*τ**z*_{b}) and deposit the wave momentum into the mean flow according to Eq. (21).The

*residual*WMF is passed to the next model level until meeting another wave breaking level or reaching the model top.

For the practical use of OGWD parameterization, the efficiency of computation is of paramount importance. To quantify the effect of horizontal propagation, one needs to calculate the wave fields of *u*, *υ*, and *w* at the breaking level. This is time consuming, especially when the model grid cell is divided into a large number of subgrids. A coarser SSO (i.e., fewer subgrid-scale grid points) is thus preferred. Furthermore, on determining the wave breaking level, one may only need to compute the wave amplitude over the largest SSO, because the peak wave amplitude tends to stay near the terrain summit (E15a).

The accuracy of the GBA solution is sensitive to the spectral resolution (Δ*k* and Δ*l*) determined by the domain size. We have examined the WMF in constant wind over a circular bell-shaped mountain as an example, for three different domain sizes (i.e., D1 = 1000 × 1000 km^{2}, D2 = 500 × 500 km^{2}, and D3 = 250 × 250 km^{2}). The WMF within D3, which is of the coarsest spectral resolution, is underestimated most notably (not shown). To increase the spectral resolution, the “zero padding” technique can be used when performing the Fourier transform of SSO (i.e., adding a number of zeros beyond the side boundaries). As such, the SSO spectrum is sampled at finer resolution, with its distribution unchanged.

It is noteworthy that the waves propagating out of the model grid cell are discarded. Thus, the new parameterization scheme may fail in the case of a high-resolution model, given the significant fictitious loss of wave momentum.

### b. Comparison with traditional OGWD parameterization scheme

The above OGWD parameterization scheme is first tested for OGWs generated in the wind profile wCUR1 over a circular bell-shaped mountain with *h*_{m} = 500 m and *r*_{a} = 50 km. Figure 12a presents the vertical variation of the wave amplitude from the GBA solution, with the height decay of density considered. For comparison, the wave amplitude ignoring the horizontal propagation of 3D mountain waves is also shown (termed as 2D wave amplitude), which is obtained from the conservation of wave action [cf. Eq. (13) in E15a]:

The second approximate equality holds when the wind curvature effect is neglected. Here, is the wave amplitude at the surface, which equals to the mountain height for linear waves. The threshold of is used for wave breaking such that the saturation wave amplitude is [cf. Eqs. (11) and (12) in Palmer et al. (1986)]

As shown in Fig. 12a, the 2D wave amplitude keeps increasing with height and exceeds the saturation wave amplitude at about *z* = 5.2 km. According to Eq. (22), the amplification of wave amplitude is attributed to the decreasing wind speed. In contrast, the GBA wave amplitude exhibits a decrease with height below about *z* = 5 km such that the waves saturate at a higher altitude of about 6.7 km. Upon saturation, mountain waves break and deposit the excessive wave momentum into the mean flow, leading to an OGWD: that is,

Figure 12b shows the mean-flow acceleration caused by breaking waves. Compared to the traditional scheme, the deposition of wave momentum produces a larger OGWD at a higher altitude in the new scheme.

The proposed new scheme is also tested using the realistic wind and buoyancy frequency and the Hainan Island rescaled to 500 m. In this case, there is no critical level. The GBA wave amplitude is always smaller than the saturation wave amplitude (Fig. 13), indicating no wave breaking and OGWD (at least below *z* = 30 km). The waves may break at a higher altitude (such as in the mesosphere), given the salient decrease of density there. For the 2D wave amplitude, it is greater than the saturation wave amplitude between about *z* = 2.5 and 8.5 km, as well as above about 16 km. The growth in the mid- to lower troposphere is mainly attributed to the decreasing wind speed and buoyancy frequency (Fig. 2), whereas the height decay of density plays a more vital role in producing the large amplitude in the stratosphere. The traditional scheme mistakenly deposits the wave momentum in the lower stratosphere, which has been found in previous studies (e.g., Milton and Wilson 1996). Clearly, this excessive stratospheric OGWD can be alleviated by the new scheme based on the GBA.

## 6. Summary

The effect of horizontal propagation of three-dimensional (3D) mountain waves on the wave momentum flux (WMF) and orographic gravity wave drag (OGWD) is investigated using the Gaussian beam approximation (GBA). For one thing, it reduces the wave amplitude with height and hence increases the wave breaking altitude (E15a,b). For another, it decreases the amount of WMF that can be attenuated by the local mean flow within finite domains (such as a GCM grid cell). Both of these two effects are missing in current OGWD parameterization schemes.

Under the WKB approximation and GBA, analytical expressions are derived for hydrostatic nonrotating mountain waves, with the height variation of stratification and wind vertical curvature considered. The behaviors of mountain waves and WMF are investigated for both idealized circular bell-shaped mountain and the realistic Hainan Island. Three unidirectional wind profiles are examined, including a constant wind, an idealized wind of vertical curvature, and a realistic wind derived from ERA-Interim data over Hainan Island. Furthermore, two potential temperature profiles are examined, with one of constant buoyancy frequency and the other also derived from the ERA-Interim dataset. The GBA solutions are validated against numerical simulations performed using the ARPS model.

Because of horizontal propagation of 3D mountain waves, the WMF decreases with height within a finite domain, even in the absence of wave breaking and/or other dissipative processes. The GBA solution can well capture the decay of finite-domain WMF with height, especially for idealized terrain, wind, and buoyancy frequency. The effect of wind curvature on improving the WMF is also captured, which cannot be omitted at small Richardson numbers of order unity. In comparison with TMV04’s second-order WKB solution, the first-order WKB approximation used in our GBA solution tends to overpredict the WMF, but these two solutions agree very well at Richardson numbers greater than 5. For mountain waves generated in realistic wind and buoyancy frequency, the WMF is underestimated as compared to its numerical counterpart, since the WKB-based GBA solution cannot represent wave reflection at the atmospheric internal boundaries (e.g., tropopause). In the case of mountain waves forced by Hainan Island, which contains several isolated hills, more Gaussians are needed in the GBA solution to better represent the small-scale waves and the WMF.

According to the GBA, a new OGWD parameterization scheme is proposed and compared with the traditional scheme that does not take into account the horizontal propagation of 3D mountain waves. In the case of idealized wind (with a critical level) and constant buoyancy frequency, the wave amplitude keeps increasing with height in the traditional scheme. However, the wave amplitude obtained from the GBA solution first decreases with height, then increases rapidly on approaching the critical level. Consequently, the waves break at a higher altitude and produce a greater OGWD than in the traditional scheme. In another case of realistic wind (without critical level) and buoyancy frequency, the GBA wave amplitude is largely reduced at high altitudes because of the horizontal geometrical spreading. As a result, there is no wave breaking in the lower stratosphere, which can help alleviate the excessive stratospheric OGWD found when using the traditional OGWD parameterization scheme.

The horizontal propagation of mountain waves may produce a nonlocal OGWD that influences the mean flow far away from the mountain. For instance, gravity waves triggered in the southern Andes and Antarctic Peninsula were found to contribute remarkably to the stratospheric OGWD over the Southern Ocean near 60°S through meridional propagation (e.g., Jiang et al. 2014, and references therein). In order to represent the propagating waves, there have been a few studies in the past (e.g., Marks and Eckermann 1995; Hasha et al. 2008; Song and Chun 2008) that have used the standard ray-tracing theory. This issue is not considered in the currently proposed OGWD parameterization scheme but will be revisited using the GBA.

The GBA solution derived can also address the selective critical-level filtering caused by directional wind shear (PR11). In accordance with XWX12, the absorption of 3D mountain waves at selective critical levels will result in a body force normal to the mean flow, which is not included in any existing OGWD parameterization scheme. In this situation, it is the 3D divergence of horizontal momentum flux [i.e., and ] that is important in forcing the local mean flow (e.g., Miyahara 2006). This will be studied in the future.

## Acknowledgments

We appreciate Dr. Manuel Pulido for his help with the Gaussian beam approximation. We are also very grateful to the anonymous reviewers for their constructive comments and suggestions. This study is jointly supported by National Fundamental Research 973 Program of China (2015CB452801 and 2013CB430103) and National Science Foundation of China (Grant 41505046).

## REFERENCES

*Waves in Fluids.*Cambridge University Press, 504 pp.

*Parallel Computations*, G. Rodrigue, Ed., Academic Press, 51–83.

## Footnotes

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).