## Abstract

A time-dependent generalization of a Fourier-ray method is presented and tested for fast numerical computation of high-resolution nonhydrostatic mountain-wave fields. The method is used to model mountain waves from Jan Mayen on 25 January 2000, a period when wavelike cloud banding was observed long distances downstream of the island by the Advanced Very High Resolution Radiometer Version 3 (AVHRR-3). Surface weather patterns show intensifying surface geostrophic winds over the island at 1200 UTC caused by rapid eastward passage of a compact low pressure system. The 1200 UTC wind profiles over the island increase with height to a jet maximum of ∼60–70 m s^{−1}, yielding Scorer parameters that indicate vertical trapping of any short wavelength mountain waves. Separate Fourier-ray solutions were computed using high-resolution Jan Mayen orography and 1200 UTC vertical profiles of winds and temperatures over the island from a radiosonde sounding and an analysis system. The radiosonde-based simulations produce a purely diverging trapped wave solution that reproduces the salient features in the AVHRR-3 imagery. Differences in simulated wave patterns governed by the radiosonde and analysis profiles are explained in terms of resonant modes and are corroborated by spatial ray-group trajectories computed for wavenumbers along the resonant mode curves. Output from a nonlinear Lipps–Hemler orographic flow model also compares well with the Fourier-ray solution horizontally. Differences in vertical cross sections are ascribed to the Fourier-ray model’s current omission of tunneling of trapped wave energy through evanescent layers.

## 1. Introduction

When suitable environmental conditions exist, flow over mountains generates quasi-stationary gravity waves that can propagate obliquely away from the parent orography. These waves can have important effects on other atmospheric processes. Air parcels advected through these waves experience rapidly oscillating adiabatic cooling and heating that can have strong net influences on cloud physics and associated chemistry in both the troposphere and stratosphere (e.g., Jensen et al. 1998; Fueglistaler et al. 2003). These influences sometimes take the form of spectacular wave-banded cloud displays that are visible from space (e.g., Fritz 1965; Gjevik and Marthinsen 1978; Burroughs and Larson 1979; Sharman and Wurtele 1983; Mitchell et al. 1990; Worthington 2001). Mountain waves grow in amplitude with altitude and can break, generating drag forces that affect the synoptic-scale circulation and turbulence that mixes chemical species (e.g., Kim et al. 2003).

In addition to these atmospheric effects, mountain waves present hazards to aviation. Severe structural damage and injuries to passengers and crew can occur when aircraft fly through severe clear-air turbulence (CAT) produced by mountain-wave breaking (e.g., Bacmeister et al. 1994; Ralph et al. 1997; de Villiers and van Heerden 2001). Furthermore, sudden changes in flight altitude caused by mountain waves, either through wave-induced CAT and/or short-wavelength vertical wave displacements, are also an important issue given recent enaction of reduced vertical separation minima (RVSM) for commercial air traffic. The large-amplitude mountain-wave event over Colorado documented by Lilly (1978) provides a vivid illustration of both of these hazards, causing a commerical airliner to drop ∼4000 ft in a little over 1 min while simultaneously undergoing severe airframe buffeting due to wave-induced CAT.

Thus, mountain wave events are important to forecast, yet the spatial scales of these waves and the subwavelength instabilities that lead to breaking, drag, and turbulence generation are generally too short for operational numerical weather prediction (NWP) models to resolve fully (e.g., Benoit et al. 2002; Smith 2004). For example, subgrid-scale gravity wave drag must be parameterized in NWP and climate models (e.g., Kim et al. 2003). Until operational NWP model resolutions improve, alternative forecasting algorithms for mountain waves must be developed.

Ongoing efforts in this area have been undertaken at the Naval Research Laboratory for over a decade in developing the MountainWave Forecast Model (MWFM; Eckermann et al. 2004, 2006). Instead of simulating the fully nonlinear discretized Navier–Stokes equations, as in a traditional NWP model, the MWFM approach uses ray methods to simulate the generation, propagation and breakdown of mountain waves within the large-scale environment specified by operational NWP model output. Since ray solutions are easy to interpret, computationally fast, and valid within a broad variety of atmospheric environments, they are attractive as potential forecasting algorithms.

The first version of the MWFM (MWFM-1) tested these ideas using a hydrostatic two-dimensional spatial-ray formulation (Bacmeister et al. 1994). The MWFM-2 was extended to use three-dimensional spatial-ray equations governed by a nonhydrostatic dispersion relation including rotation (Eckermann and Preusse 1999). Both models employ an idealized ridge decomposition of the earth’s topography to define the major terrain features relevant for generation of waves near the surface. MWFM forecasts have been used for over a decade now to direct National Aeronautics and Space Administration (NASA) research aircraft away from hazardous mountain-wave-induced turbulence and into regions where nonbreaking waves generate cloud, tasks for which it has repeatedly shown skill (e.g., Bacmeister et al. 1994; Eckermann et al. 2004, 2006). MWFM-2 hindcasts of stratospheric wave amplitudes have been validated globally using new data from stratospheric research satellites (Eckermann and Preusse 1999; Jiang et al. 2002, 2004), as well as regionally using various suborbital measurements (e.g., Hertzog et al. 2002; Eckermann et al. 2006). MWFM forecasts of stratospheric mountain-wave turbulence were utilized extensively by the U.S. Air Force during Operations Enduring Freedom and Iraqi Freedom (Eckermann 2002), and since 2004 MWFM-2 has been run operationally at the Air Force Weather Agency. MWFM-2 hindcasts have also been used in a variety of research applications, such as the role of mountain waves in accelerating ozone loss chemistry in the Arctic winter stratosphere (e.g., Carslaw et al. 1999; Pierce et al. 2003; Svendsen et al. 2005; Mann et al. 2005).

Thus, ray methods are now an accepted approach to forecasting mountain waves (Eckermann et al. 2004). Yet the current MWFM ray algorithms still contain a number of significant simplifications and shortcomings. First, they use a one-dimensional height-dependent form for the wave action equation that does not include the effects of horizontal geometrical spreading of the rays on the wave amplitude evolution (e.g., Shutts 1998). Second, they use a spatial formulation for the ray solutions, which leads to caustic singularities that are not practical to correct (e.g., Broutman et al. 2001, 2002, 2004). Third, they use idealized ridge databases for their source functions that do not capture the full spectrum of waves radiated by flow over realistic topography.

To fully explore and exploit the capabilities of ray methods for mountain-wave forecasting, over the past several years we have progressively developed improved ray algorithms that reduce or eliminate these and other weaknesses in the current MWFM ray equations (Broutman et al. 2001, 2002, 2003, 2004, 2006). This has led to a new ray algorithm that we refer to here as the Fourier-ray method, because it involves Fourier-synthesized (rather than spatially synthesized) ray solutions. The method has recently been reviewed inter alia by Broutman et al. (2004) and the version of it that we use here is described in section 2a.

This method shows promise as a potential next-generation dynamical ray core for the MWFM, since the idealized solutions derived to date alleviate or eliminate all of the aforementioned weaknesses of the corresponding spatial ray solutions. Specifically, the algorithm incorporates arbritrary topographic forcing at the lower boundary, models both trapped and free-propagating waves, includes horizontal geometric spreading in its wave action solutions, and systematically corrects the caustic singularities found in spatialray solutions. To date, however, we have applied these new algorithms only to idealized mountain-wave problems that use analytical functions for the obstacle shape and background wind profile.

Here we take our next step in assessing the potential of the Fourier-ray method for forecasting by extending and applying these algorithms to a more realistic case study. The major extensions here are incorporation of realistic topography and real background wind and temperature profiles from both radiosonde observations and a meteorological analysis issued by a data assimilation system. We concentrate here on nondissipative wave solutions, deferring developments in estimating locations of wave-induced CAT to later studies.

Our modeling focuses on Jan Mayen (71°N, 8.4°W), a small island in the North Atlantic whose topography is dominated to the north by Mount Beerenberg, a quasi-circular volcanic mountain of width ∼5–10 km and a peak elevation of ∼2270 m (see Fig. 1). We seek to “hindcast” wavelike patterns observed over the island in satellite cloud imagery on 25 January 2000. Previous analysis and linear modeling of some earlier satellite cloud images over Jan Mayen (Gjevik and Marthinsen 1978; Simard and Peltier 1982) associated cloud banding here with three-dimensional trapped lee waves forced by flow over Mount Beerenberg that radiated into much larger atmospheric volumes downstream. Such waves present an excellent test case for our new algorithm’s vertical reflection and high-resolution forecasting capabilities, since Jan Mayen’s topography would not be adequately resolved by current operational NWP systems. Our results here are used to benchmark the initial performance of this new ray code and to target areas for further development for future forecasting applications.

## 2. Models

### a. The Fourier-ray model

Our primary tool is a flexible numerical implementation of a Fourier-ray method. The name “Fourier ray” refers to Fourier-synthesized ray solutions: that is, the ray solutions are computed in a Fourier domain and then superimposed by inverse Fourier transform to give a spatial solution. The method is described in Broutman et al. (2002, 2003, 2006). During earlier stages of its development, different names were used for the same basic method (e.g., Maslov’s method, a simplified Fourier method). The reasons for the name changes are discussed in Broutman et al. (2006).

The ray solution in the Fourier domain is expressed as a function of *k*, *l*, *z*, *t*, where *k* and *l* are the horizontal wavenumbers, *z* is height, and *t* is time. The advantage of solving in the Fourier domain is that *k* and *l* are constant along the ray in a horizontally uniform background, as assumed here. Along with a simple ray-based treatment of the time dependence of the solutions (Broutman et al. 2006), the problem is effectively reduced to a one-dimensional calculation in *z*. This is a major simplification for numerical ray tracing, especially for the wave amplitude calculations and the correction of caustics.

Broutman et al. (2002, 2003) derived steady-state Fourier-ray solutions for hydrostatic and nonhydrostatic waves, respectively. Wave transience was incorporated into the Fourier-ray method by Broutman et al. (2006). We use that transient formulation here to aid direct comparisons with output from nonlinear numerical models at finite times. Besides being helpful for estimating setup times for the wave field, the transient formulation has two computational advantages. First, there are no resonant singularities in the transient solution for the trapped waves, as there are in the steady-state solution. (The steady-state singularities can be removed by adding a small imaginary component to the wavenumber or frequency, but this artificially damps wave amplitudes. The Fourier-ray solutions presented here are nondissipative.) Second, the transient wave field is more limited in spatial extent than the steady-state solution, which is the longtime limit of the transient solution. The transient solution can thus be represented with fewer computational grid points.

We consider linear three-dimensional mountain waves radiated from realistic topography into arbitrary vertical profiles of winds and stability. The mountain waves are stationary, with zero frequency (*ω* = 0) relative to the ground. For a horizontal background wind **U** = (*U*, *V*, 0), the intrinsic frequency is *ω̂* = *ω* − **k** · **U** = −*kU* − *lV*, where **k** = (*k*, *l*, *m*) is the wavenumber vector. The background wind and buoyancy frequency *N* are assumed to depend on height but not on horizontal position or time.

The ray tracing performed here is based on a gravity wave dispersion relation of the form

where *k _{h}* = (

*k*

^{2}+

*l*

^{2})

^{1/2}. We ignore the effects of the earth’s rotation and compressibility in (1) for simplicity: they are included in other versions of the code.

We define *w̃*(*k*, *l*, *z*, *t*) to be the ray solution in the Fourier domain for the vertical velocity associated with the mountain waves. The corresponding spatial solution *w*(*x*, *y*, *z*, *t*) is obtained by inverse Fourier transform:

The Fourier-ray solution for *w̃*(*k*, *l*, *z*, *t*) is given in the appendix, based on derivations in Broutman et al. (2006). We emphasize that *w*(*x*, *y*, *z*, *t*) is not the same as the spatial-ray solution for the vertical velocity, which is obtained from the stationary phase limit of the inverse Fourier transform in (2). The distinction is important because the stationary phase approximation, but not *w*, breaks down at caustics in the spatial domain (see, e.g., Shutts 1998; Broutman et al. 2001).

This Fourier-ray code ingests arbitrary vertical profiles of wind and temperature, which are interpolated onto a vertical grid sufficiently fine to evaluate the following phase integrals accurately:

where *z _{t}* is the height of the turning point for a given (

*k*,

*l*). The first integral is used for vertically propagating waves, the second integral is used for vertically trapped waves at heights below their turning point

*z*, and the third integral is used for vertically trapped waves at heights above their turning point (see the appendix). An absolute value appears in the third integral because

_{t}*m*is imaginary above

*z*. These integrals are computed here numerically by the trapezoidal rule, using a vertical grid size Δ

_{t}*z*= 0.25 km. Comparisons with solutions using smaller Δ

*z*showed that this grid size was adequate.

The code uses a finite Fourier-series approximation to the inverse Fourier transform (2). In the cases presented here, 512 wavenumber values were chosen for *k* and *l*, corresponding to a spatial grid of 1024 by 1024 in *x* and *y*, and a horizontal grid spacing of Δ*x* = Δ*y* = 1 km.

### b. Nonlinear orographic flow model

To assess more rigorously the accuracy of the Fourier-ray solutions to follow, we also performed some companion simulations using a fully nonlinear numerical model. The model in question solves the incompressible nonlinear Navier–Stokes equations over topography using the Lipps–Hemler anelastic approximation (Lipps and Hemler 1982; Nance and Durran 1994).

The model uses terrain-following coordinates (Gal-Chen and Somerville 1975), a centered second-order spatial discretization, and a modified leapfrog time-stepping scheme (Rees 1988). It was run here on a staggered Arakawa C grid of 256 points in each horizontal direction and 32 points in the vertical, with a grid spacing of 1 km in all directions and a time step of 4 s. A free-slip (frictionless) lower boundary condition was used. To absorb waves at heights ≳20 km, we imposed Rayleigh damping with a rate coefficient that increased linearly from zero at 20 km altitude to 4 10^{−3} s^{−1} at the 32-km upper boundary. To minimize reflections from the lateral boundaries, a radiative scheme was used (Miranda and James 1992) as well as some linearly increasing viscous damping within a 10-km-wide region at each side boundary. Away from these upper and side boundaries, the only diabatic terms were a Richardson-number-dependent first-order turbulent closure scheme (Lilly 1962) and a small amount of fourth-order diffusion to suppress grid-scale noise (hyperviscosity coefficient of 10^{7} m^{4} s^{−1}).

### c. Topography

Topographic elevations *h*(*x*, *y*) for Jan Mayen were obtained from the United States Geological Survey (USGS) GTOPO30 database. This global dataset of digital terrain elevation has 30 arc s resolution, corresponding roughly to 1-km spatial resolution, which we interpolated linearly to our horizontal grid spacing of 1 × 1 km^{2} using Delaunay triangulation and full spherical to Cartesian transformations. This model topography for Jan Mayen is plotted in Fig. 2a and can be used in the 1 km × 1 km Fourier-ray model runs.

For the 1 km × 1 km Lipps–Hemler model runs, however, this 1 km × 1 km topography must be smoothed to reduce the potential for unphysical four-point gridpoint noise in the simulations (see, e.g., Davies and Brown 2001). We achieve this using a two-dimensional five-point running average, and the results are shown in Fig. 2b. In the Fourier-ray model runs shown here we also use this smoothed topography to facilitate more direct comparisons with the Lipps–Hemler model output.

## 3. Observational data

### a. Advanced Very High Resolution Radiometer Version 3

The Advanced Very High Resolution Radiometer Version 3 (AVHRR-3) is a cross-scanning passive radiometer, first deployed on the National Oceanic and Atmospheric Administration (*NOAA-15*) polar-orbiting satellite. It has six channels: three thermal infrared (IR) channels and three other channels in the visible and near IR. The system acquires radiances with a 20.32-cm-diameter telescope that scans cross track using a continuously rotating mirror. The telescope’s instantaneous field of view (IFOV) is 0.0745°, corresponding to surface footprint diameters across track and along track of 1.1 km × 1.1 km at nadir and 6.2 km × 2.3 km at the far off-nadir positions. For each scan, 2048 samples are acquired at off-nadir angles between ±55.37° cross track of the subsatellite point, corresponding to a total horizontal swath distance at the surface of ∼2900 km. For further technical details, see section 4.1.1 of Kidder and Vonder Haar (1995) and section 3.1 and appendix J.1 of Goodrum et al. (2000).

This intrinsic data resolution, known as local area coverage (LAC), is too dense to be continuously stored and telemetered to ground stations. Thus, LAC data can only be received at certain scheduled times. Continuous monitoring at all other times is provided by global area coverage (GAC) data, which are onboard averages of a reduced subset of the raw LAC data as described in section 10.2.1 of Kidder and Vonder Haar (1995) and section 3.1.3.2 of Goodrum et al. (2000). Surface footprint diameters for GAC data are ∼4 × 3.3 km^{2} at nadir.

Data from each channel are converted on board into 10-bit binary earth scene counts *C _{E}* prior to being telemetered to ground stations. Earth scene radiances

*R*are derived using the nonlinear radiance correction formula,

_{E}where *a*_{0}, *a*_{1}, and *a*_{2} are regression coefficients derived from prelaunch calibration data (Walton et al. 1998).

### b. Radiosonde data

A meteorological station on Jan Mayen launches radiosondes at 0000 and 1200 UTC each day: see Fig. 1 of Gjevik and Marthinsen (1978) for its precise location on the island. Here we use winds and temperatures acquired from the 1200 UTC sounding on 25 January 2000 in our model simulations, which we reinterpolate onto a high vertical resolution grid, vertically smooth using a running average of 2-km width, then reinterpolate again onto the models’ vertical grids. While we assume here that these data approximate the true vertical profile over the island at this time, the actual balloon trajectory is oblique: given a typical ascent velocity of 5 m s^{−1} (e.g., Lane et al. 2000), passive advection calculations based on radiosonde winds for this day place the balloon ∼100 km downstream by the time it reached 10-km altitude.

### c. Meteorological analyses

To provide some cross validation for the model runs, we also use profiles of winds and temperatures at the closest grid point to Jan Mayen from the 12-hourly Met Office (UKMO) 3.75° × 2.5° analyses (Swinbank and O’Neill 1994). To define surface weather patterns, we use 6-hourly analyzed mean sea level pressures (MSLP) from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) 2.5° × 2.5° reanalyses (Kalnay et al. 1996).

## 4. Jan Mayen wave clouds of 25 January 2000

### a. AVHRR-3 observations

Figure 3 plots a chronological sequence of selected AVHRR-3 channel-5 (∼11.5–12.5 *μ*m) radiances on 25 January 2000. The data have been reinterpolated to Cartesian coordinates (*x* being east–west), with the island of Jan Mayen centered at *x* = *y* = 0 on each map and its coastline plotted in each panel. We use only the thermal IR data since Jan Mayen is in polar night at this time of year.

At ∼1140 UTC, Fig. 3a shows that the island was obscured by cloud decks. Around 3 h later at 1500 UTC these clouds have been advected to the east to reveal evidence of mountain-wave banding of lower-level clouds: at this time we have both conventional GAC data (Fig. 3b) and high-resolution LAC data (Fig. 3c). The LAC image shows cloud banding superficially consistent with a purely diverging three-dimensional ship mountain-wave pattern (Sharman and Wurtele 1983) emanating from Jan Mayen, though only one “arm” of this V-shaped pattern is fully visible, the other still being partially obscured by overlying cloud decks.

Around 90 min later at 1630 UTC, GAC radiances in Fig. 3d again show a diverging wavelike cloud pattern. Though this is a lower-resolution image than the LAC image in Fig. 3c, it provides our least-obscured image of the full wave pattern. In fact this wave pattern extends farther downstream in this image than might be suggested by the limited domain plotted in Fig. 3d.

Figure 3e plots another AVHRR-3 GAC image acquired roughly 2 h later at ∼1820 UTC. Upper-level cloud decks have moved in and partially obscured the wave pattern. Nonetheless, we still see evidence of the lower (southward) arm of the diverging wave pattern in the bottom half of the image.

Figure 3f and a sequence of later images (not shown) do not show evidence of banded clouds downstream of Jan Mayen. While some, like Fig. 3f, are heavily impacted by obscuring cloud layers, the weight of evidence from these images suggests that the banded wave clouds have disappeared at these later times.

Thus, from the AVHRR-3 channel-5 IR imagery in Fig. 3, we deduce an approximate 6-h duration of this apparent Jan Mayen wave event, lasting from ∼1200 to 1800 UTC. The role of high-cloud decks in obscuring the wave banding patterns also indicates that the banded wave clouds occur in the lower or midtroposphere.

### b. Synoptic surface meteorology

Figure 4 plots 6-hourly NCEP–NCAR reanalyzed MSLP in a broad region centered over Jan Mayen from 0600 UTC on 25 January to 0000 UTC on 26 January. We see rapid eastward passage of a compact low pressure cell across the domain. As the core of this low passed to the north of Jan Mayen at ∼1200 UTC, Fig. 4b shows that it produced enhanced surface geostrophic westerlies across Mount Beerenberg. These surface westerlies persist over Jan Mayen at 1800 UTC, though the isobars in Fig. 4c reveal that the core of the low has now passed well to the east of the island, and that surface westerly flow across Mount Beerenberg is about to weaken and be replaced by weaker northerly flow as the low moves farther eastward. This is confirmed by the isobars at 0000 UTC in Fig. 4d.

Thus, the appearance of cloud banding from 1200 to 1800 UTC in the AVHRR-3 radiances in Fig. 3 correlates with a period of enhanced surface westerlies across Mount Beerenberg (and hence enhanced potential for mountain-wave forcing) accompanying the eastward passage of a polar low to the north of Jan Mayen.

### c. Meteorological profiles over Jan Mayen at 1200 UTC

Figure 5 plots the three-dimensional horizontal wind profile over Jan Mayen from the routine 1200 UTC radiosonde ascent on 25 January (solid curve), as well as the nearest gridbox 1200 UTC profile from the UKMO analysis. The radiosonde profiles verify the presence of strong near-westerly surface-level winds inferred from Fig. 4b, with speeds ∼30 m s^{−1}. These wind speeds increase with altitude to a jet maximum of ∼70 m s^{−1} near 10-km altitude. This maximum is nearer 55 m s^{−1} in the UKMO profile.

To produce mountain-wave activity so far downstream of Jan Mayen in Fig. 3, it is likely that these waves were vertically trapped (e.g., Simard and Peltier 1982). Strongly sheared wind profiles like those in Fig. 5 are one way to produce the vertical wave reflection required for trapping. To assess this, we compute profiles of the Scorer parameter. The standard Scorer parameter is a two-dimensional term in which the horizontal wavenumber vectors and the wind vector are coaligned. For three-dimensional problems, a range of angles for the horizontal wavenumber vector exist and the wind vector can rotate with altitude, and thus no unique Scorer parameter exists. Some insights can be gained, however, by studying a generalized form for three-dimensional problems,

where *N*(*z*) is the background buoyancy frequency profile, *Ũ*(*z*) = *U*_{tot}(*z*) cos[*α*(*z*) − *ϕ*] is the profile of the component of the wind vector *U*_{tot}(cos *α*, sin *α*) projected along a given horizontal wavenumber vector’s direction *ϕ* = arctan (*l*/*k*), *Ũ _{zz}* is the curvature in this projected component wind profile, and

*U*

_{tot}(

*z*) = [

*U*

^{2}(

*z*) +

*V*

^{2}(

*z*)]

^{1/2}. An equivalent expression is utilized by Vosper and Worthington (2002) and Eckermann et al. (2006). From the profiles in Fig. 5, specifically the hodographic projection at the surface, we see that the horizontal wind vector is directed at an approximately constant angle of

*α*= 20° north of due east through most of the troposphere.

Figures 6a,b plot component profiles of wind speed and curvature, respectively, for both the radiosonde and UKMO profiles on 25 January at 1200 UTC, using a choice for the direction of the horizontal wavenumber of *ϕ* = *α* = 20° (i.e., wave vectors aligned parallel to the wind vector). Figure 6c plots the buoyancy frequency profiles. The thin black curve in Fig. 6a shows the raw radiosonde wind profile, while the thick black curve shows that same profile after high-resolution resampling and then smoothing using a sliding 2-km vertical averaging window to remove the small-scale structure. This filter is applied to all the profiles (both radiosonde and UKMO) that enter into the Scorer parameter calculation (5), since even slight vertical structure in the wind and temperature profiles gets greatly enhanced during the numerical gradient and curvature calculations in (5), yielding large spurious oscillations in ℓ^{2}(*z*, *ϕ*) (see, e.g., Danielsen and Bleck 1970; Shutts 1992; Smith 2004).

The thick curves in Fig. 6d show Scorer parameters ℓ^{2}(*z*, 20°) calculated by ignoring the curvature term in (5), since the first term *N* ^{2}(*z*)/*Ũ*^{2}(*z*) is generally the most significant. Furthermore, curvature terms are omitted at present from our Fourier-ray model’s wave equations, so this simplified term is roughly equivalent to this model’s Scorer parameter. (The Scorer parameter and the present version of the Fourier-ray model omit Coriolis terms, but their influence on these short-wavelength trapped wave solutions is minimal.) The corresponding dotted curves show ℓ^{2}(*z*, 20°) with the curvature profiles in Fig. 6b retained in the calculations. They produce fairly minimal modifications to the simplified *N* ^{2}(*z*)/*Ũ*^{2}(*z*) profile. Tests that include additional wind shear and density-scale height terms in a generalized Scorer parameter derived from fully compressible wave equations (Nance 1997) also yielded minimal changes to these simplified *N* ^{2}(*z*)/*Ũ*^{2}(*z*) profiles.

Figure 6d shows that the wind and stability profiles over Jan Mayen on 25 January at 1200 UTC yield a deep minimum in ℓ^{2}(*z*, *ϕ*) at altitudes ∼5–10 km. Turning points *z _{t}* for a given mountain wave of total horizontal wavelength

*λ*= 2

_{h}*π*/

*k*occur where ℓ

_{h}^{2}(

*z*,

*ϕ*) −

*k*

^{2}

_{h}= 0. The dashed curves in Fig. 6d plot ℓ

^{2}(

*z*, 20°) −

*k*

^{2}

_{h}for

*λ*= 30 km. We see that the radiosonde-derived Scorer parameter predicts reflection of this horizontal wavelength at

_{h}*z*≈ 5 km, whereas the UKMO profile predicts no reflection (due to its weaker peak winds), though it comes very close to yielding reflection at ∼8 km. Since waves of comparable or shorter horizontal wavelengths are likely to be generated by flow over the horizontally narrow obstacle presented by Mount Beerenberg (see Fig. 1), we see that Jan Mayen mountain waves probably reflected vertically in the midtroposphere at 1200 UTC. Since we have a stable boundary layer (Fig. 6c) containing fairly strong surface flow (Fig. 6a), these reflected waves will reflect again at the surface, rather than being absorbed (Smith et al. 2002), yielding a vertically trapped wave that can penetrate long distances downstream within this waveguide, as observed in Fig. 3.

_{t}The surface wind speed *U*_{0} is ≈30 m s^{−1} (Fig. 6a) and surface buoyancy frequency *N*_{0} is ≈0.015 rad s^{−1} (Fig. 6c). Since the maximum height of the Jan Mayen topography *h _{m}* is ≈2 km (Fig. 2), this yields Fr

^{−1}≈ 1, where Fr =

*U*

_{0}/

*N*

_{0}

*h*is the surface Froude number. This suggests a near-maximum amplitude for the forced wave, with the surface flow passing over (rather than around) the mountain and translating all of the mountain pressure drag into a gravity wave response. Since Fr

_{m}^{−1}is close to unity, some weak gravity wave breaking and/or associated lee-vortex nonlinearity may occur in a thin layer near the surface.

## 5. Model results

### a. Fourier-ray solutions

We begin by studying how horizontal cross sections of our Fourier-ray model solutions compare with the AVHRR-3 imagery in Fig. 3. Our best (least obscured) AVHRR-3 image of the wavelike cloud banding in Fig. 3d is replotted in Fig. 7a as our observational reference. Below it, we plot two Fourier-ray solutions for vertical velocity at *z* = 3 km and *t* = 4 h, as calculated using the Jan Mayen radiosonde profiles (Fig. 7b) and the UKMO profiles (Fig. 7c).

Both simulations show similarities with the AVHRR-3 image in Fig. 7a. The major response in each case is a V-shaped trapped mountain-wave pattern extending long distances downstream at angles to the island similar to those observed in Fig. 7a. There are, however, some noticeable differences between the two solutions.

Sharman and Wurtele (1983) showed that analytical solutions for trapped mountain waves from three-dimensional symmetric obstacles have properties and sensitivities that resemble classical Kelvin solutions for ocean surface ship wakes (Reed and Milgram 2002). Like ship wakes, trapped mountain waves can exhibit both transverse waves, which appear downstream of the obstacle as linear phase lines orthogonal to the flow direction, and diverging waves, which appear as curved phase lines in two V-shaped wedges of activity either side of the downstream flow axis.

The Fourier-ray model simulation based on the UKMO profiles in Fig. 7c contains both diverging and transverse waves, but the solution based on the Jan Mayen radiosonde profile in Fig. 7b contains only diverging waves. The AVHRR-3 images in Figs. 3 and 7a show clear evidence of diverging waves but no evidence of transverse waves in the cloud banding: see especially the high-resolution LAC image in Fig. 3c. This indicates that the (presumably) more accurate representation of the meteorological profiles over Jan Mayen from the radiosonde ascent is important for an accurate simulation of the observed disturbance. We will have more to say in section 6b about why the transverse waves are present for the UKMO profile simulations but not for the radiosonde profile simulations.

Other aspects of the radiosonde-based Fourier-ray solution also agree better with the imagery than does the UKMO-based solution. For example, the radiosonde-based solution shows greater “flaring” of the diverging waves into progressively wider wedge angles with increasing distance downstream. By wedge angle, we mean the angle relative to a line passing from the center of the topography through the approximate center of the downstream wave field. There is also a clear horizontal wavelength dependence to the diverging waves, with the shortest wavelengths confined to the smaller wedge angles and the longer wavelengths extending into larger wedge angles. A similar wavelength dependence is seen in idealized numerical model solutions (see, e.g., Fig. 9a of Sharman and Wurtele 2004).

Since the solutions based on the radiosonde profile agree better with the AVHRR-3 observations, we focus on them for further analysis.

### b. Nonlinear numerical model simulation

To assess more rigorously the accuracy of the Fourier-ray solutions in Fig. 7, we performed companion runs using a nonlinear numerical model containing identical topography, as described in sections 2b and 2c. Figure 8 plots this model’s vertical velocity field at *z* = 3 km after 4 h using the radiosonde profiles. The line sloping upward from left to right indicates the horizontal coordinates of a vertical cross section that is shown in Fig. 9b.

The similarities between the numerical model solution in Fig. 8 and the corresponding Fourier-ray solution in Fig. 7b are striking. Both solutions reproduce a purely diverging trapped wave solution with similar wavelengths and phase orientation. Both also include very short wavelengths on the inside of the upper arm of the V-shaped pattern, located just below the line indicating the position of the vertical cross section, with the phases nearly parallel to that line.

Figure 9 plots height-varying cross sections of the vertical velocities along the line in Fig. 8 for the Fourier-ray model (Fig. 9a) and the nonlinear numerical model (Fig. 9b). Both show a disturbance trapped between the ground and ∼8 km altitude with similar horizontal wavelength structure between 0 and 8 km. Differences at higher altitudes and downwind are discussed below.

## 6. Discussion

Our Fourier-ray simulations of mountain waves forced by observed flow over Jan Mayen on 25 January 2000 produced vertical velocity fields *w*(*x*, *y*, *z*, *t*) at *z* = 3 km and *t* = 4 h that reproduced the salient features of high-resolution banded IR cloud imagery acquired from space by AVHRR-3. For further validation, we also compared those model results with output from companion runs using a nonlinear numerical model, again revealing good overall agreement.

On a 2.8-GHz single processor workstation, these Fourier-ray solutions took ∼5 min to compute in FORTRAN 90 on a 1024 × 1024 horizontal grid at 40 levels from the ground to 20 km. The corresponding nonlinear numerical model runs at 256 × 256 × 32 with a 4-s time step on the same workstation took 38 h to yield the *t* = 4 h solution. Since the Lipps–Hemler code we use has not been aggressively optimized, we performed identical runs on the same workstation using the Weather Research and Forecasting (WRF) Advanced Research (ARW) model (Skamarock et al. 2005). These ARW runs took 21 h to yield the *t* = 4 h solution using a 4-s time step, and 8.5 h using a 10-s time step.

Some differences among the various model solutions were noted. We investigate them further here.

### a. Vertical wave structure

The vertical cross sections of the Fourier-ray and numerical model solutions in Fig. 9, while similar in some respects, also show some differences. Trapped waves in the nonlinear numerical model solution decay in amplitude more rapidly with downstream distance than the Fourier-ray solution, while at stratospheric altitudes the nonlinear numerical model solution contains more energy in freely propagating (phase tilted) mountain waves than the Fourier-ray solution.

Both differences seem to be related to the tunnelling of trapped waves (Sutherland and Yewchuk 2004) through an evanescent layer surrounding the wind jet at heights ∼7–10 km, which then emerge into a free-propagating region above ∼10 km (see Fig. 6d). This “leakiness” of the tropospheric duct leads to greater decay with downwind distance of the vertically trapped tropospheric wave amplitudes in Fig. 9b. Similar tunnelling of trapped mountain waves into the stratosphere was noted by Eckermann et al. (2006) in a Lipps–Hemler model forecast for this same date over northern Scandinavia, a region ∼30° to the east of this Jan Mayen wave event.

Our Fourier-ray model does not yet account for tunnelling, which requires a two turning-point analysis (Bender and Orszag 1978). One turning point is at the lower edge of the evanescent region, and the other is at the upper edge. Our single turning-point Fourier-ray solution breaks down near and above the upper turning point. Here, as a temporary measure, we have set the wave amplitude to zero whenever the waves approach the second turning point. Specifically, we set the wave amplitude to zero whenever *z* > *z*_{tp} + 0.9Δ*z*_{tp}, where *z*_{tp} is the height of the lower turning point and Δ*z*_{tp} is the height difference between the upper and lower turning points. This was done only for waves that have two turning points. The solutions for the vertically propagating waves and for waves with a single turning point were left unchanged.

To include the effects of tunnelling in the Fourier-ray model, we would need either to match the ray solution to separate Airy function solutions around each of the two turning points, or to introduce a uniformly valid two turning-point solution in terms of parabolic cylinder functions [sometimes called Weber functions: see Kravtsov and Orlov (1999)]. We hope to test such solutions in future work. We have, however, estimated how much attenuation should occur through the evanescent layer in the current problem by evaluating exp(−∫^{b}_{a} |*m*| *dz*) for each *k*, *l* value that results in two turning points. Here *a*, *b* are the heights of the lower and upper turning points, respectively, and the exponential is a measure of the wave attenuation in the evanescent region between the turning points. We found that important parts of the spectrum had attenuation factors greater than 0.1, suggesting that their tunnelled wave amplitudes would be noticeable above the wind jet maximum. For example, for the *k*, *l* value considered in calculating the Scorer parameter of Fig. 6d, the value of the exponential is about 0.6.

### b. Resonant modes

While the radiosonde profile produced a purely diverging trapped wave solution in Fig. 7b, the UKMO profile produced a solution with both diverging and transverse components (Fig. 7c).

To investigate this further, we calculated the resonant modes for each of these solutions. These are the singularities of the steady-state (*t* → ∞) trapped wave solution (A11), which involve the Airy function Ai(*r*). The Airy-function argument *r*, defined in (A8), is positive in the evanescent region above the turning point, negative in the propagation region below the turning point, and zero at the turning point. Resonances occur whenever Ai(*r*_{0}) = 0, with *r*_{0} = *r*(*z* = 0). This occurs approximately at *r*_{0} = −2.34, −4.09, and −5.52, for modes 1, 2, and 3, respectively.

Contours of *r*_{0} as a function of the horizontal wavenumbers *k*, *l* are plotted in Fig. 10 for the radiosonde profiles (Fig. 10a) and for the UMKO profiles (Fig. 10b). The bold curves indicate the resonant values where Ai(*r*_{0}) = 0. A gap in, or termination of, a particular *r*_{0} contour occurs where the local *k*, *l* values do not correspond to a trapped wave. Generally, the *k*, *l* values to the upper right of the plotted contours correspond to propagating waves without turning points, while the *k*, *l* values to the lower left of the plotted contours correspond to evanescent waves whose intrinsic frequency is greater than *N* at the ground.

For the radiosonde profiles, there is a gap in the mode-1 resonance curve, centered around *k*/*k*_{0} = −0.5. For the UKMO case, the mode-1 curve is continuous. This seems to be the reason why transverse waves occur for the UKMO profile solution in Fig. 7c but not for the radiosonde profile solution in Fig. 7b.

Spatial-ray tracing illustrates the point. Figure 11 plots horizontal group trajectories of some rays with (*k*, *l*) values sampled from the corresponding mode-1 resonance curves in Fig. 10. These trajectories were computed by numerically integrating the spatial ray equations (Lighthill 1978)

out to 4 h, the same time at which the corresponding Fourier-ray solutions in Figs. 7b,c were evaluated. Here **x** is the position vector of the ray and **c*** _{g}* is its ground-based group velocity vector.

Figure 11 shows that the presence (UKMO case) or absence (radiosonde case) of mode-1 rays in the region directly downwind of the mountain corresponds to the presence or absence of transverse waves in the Fourier-ray solutions in Fig. 7. For the radiosonde profile, there are rays with wavenumbers (*k*, *l*) lying off of this mode-1 curve that propagate downwind. However, these rays are either vertically propagating (as opposed to vertically trapped), and thus unable to produce a transverse wave train at long downwind distances, or they are vertically trapped but nonresonant. Resonance is important because these are the waves that constructively interfere in the solution (A3) to produce the largest downstream amplitudes. We also checked the rays for the mode-2 resonances in the radiosonde data, and again we found an absence of rays in this same region.

Note that the spatial-ray tracing also provides other consistency checks on our Fourier-ray solutions. For example, Fig. 11 shows rays governed by the radiosonde profile “flaring” more in the diverging wings and propagating farther downstream after 4 h than those rays governed by the UKMO profile, consistent with the Fourier-ray solutions in Fig. 7.

The physical interpretation of these results is as follows. Transverse waves are absent in the radiosonde-based solution because its mode-1 resonant wavenumber magnitudes *k _{h}* in Fig. 10 are smaller than those of the UKMO-based solution. The smaller radiosonde-based mode-1

*k*values do not satisfy the reflection criterion ℓ

_{h}^{2}(

*z*,

*ϕ*) −

*k*

^{2}

_{h}= 0 for the transverse wavenumber alignments

*ϕ*, and so these rays propagate freely in the vertical, leaving a purely diverging trapped wave pattern. Conversely, the larger

*k*values on the UKMO-based resonant mode-1 curve satisfy ℓ

_{h}^{2}(

*z*,

*ϕ*) −

*k*

^{2}

_{h}= 0 for all wavenumber alignments, resulting in a trapped wave solution with both diverging and transverse components. The smaller

*k*(longer

_{h}*λ*) resonant values of the radiosonde-based solution compared to the UKMO-based solution are evident in the wave patterns in Figs. 7b,c.

_{h}## 7. Summary

We have used the Fourier-ray method to simulate mountain waves that were observed by high-resolution IR satellite cloud imagery and that extended long distances downwind of the island of Jan Mayen on 25 January 2000. This was a first step in applying the Fourier-ray method to a real-world problem. The results highlight some of the desirable features of the Fourier-ray method, such as the ability to produce very high resolution solutions quickly for realistic topography and arbitrary wind and stability profiles.

Comparison with a companion simulation using a nonlinear numerical model of flow over orography revealed similar wave patterns but also indicated the possible importance (see Fig. 9) of wave tunnelling through the wind jet centered at a height of about 10 km. Wave tunnelling (i.e., the penetration of waves through a localized evanescent region) is a two turning-point effect, which is not included in the present formulation of the Fourier-ray model. The analysis of two turning-point problems is standard (Bender and Orszag 1978), and we hope to add this capability to the Fourier-ray model in the future.

The speed and accuracy of this algorithm, coupled with its use of arbitrary surface topography and vertical meteorological profiles, make it promising as a mountain-wave forecasting tool, particularly as a potential next-generation dynamical core for the MWFM. Future work will apply the method to more complex real-world problems, such as extended topography and profiles based on operational NWP fields (see Eckermann et al. 2004).

## Acknowledgments

This research was supported by The Office of Naval Research through NRL’s base 6.2 research program, and by the NASA Geospace Sciences Program. Funding was also provided by NSF Grants ATM-0435789 (DB) and ATM-0448888 (DB and JL). Comments from Dr. James Rottman and the reviewers are appreciated.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**.**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### The Transient Fourier-Ray Solution

We consider mountain waves that begin to emerge from the mountain at *t* = 0 and are continuously generated thereafter. We assume that all rays have the full wave amplitude of the longtime steady-state solution. The resulting time-dependent ray solutions were used to produce the results in the present study and are given below. They were derived in Broutman et al. (2006), who expressed them in terms of the vertical displacement *η̃*(*k*, *l*, *z*, *t*). This is related to the vertical velocity (used in the present study) by *w̃* = −*iω̂η̃*, where *ω̂* is the intrinsic frequency.

The ray solution is the sum of solutions for the vertically propagating waves and the vertically trapped waves:

For the vertically propagating waves

For the vertically trapped waves

These are Eqs. (5) and (22), respectively, of Broutman et al. (2006), with the latter equation assuming equal numbers of upgoing and downgoing waves. The terms in these equations are defined below.

The factor *F* accounts for the time dependence in the vertically propagating waves and is defined by

The Heaviside function *H*(*ξ*) is zero for *ξ* < 0 and unity for *ξ* > 0. The integral in the argument of *H* is the height at time *t* of the initial ray (for that *k*, *l*) generated at the mountain at *t* = 0.

The topography *h*(*x*, *y*) has a Fourier transform *h̃*(*k*, *l*). The factor *G*, defined such that *G*|*η̃*|^{2} is the vertical flux of wave action, is

where *ρ* is the background density and *c _{g}*

_{3}is the vertical group velocity. Here

*G*

_{0}is

*G*at

*z*= 0.

The wave phase is

The sign convention here is *m* ≥ 0, with a negative sign inserted in (A6) to produce downward phase propagation and upward group propagation.

For the vertically trapped waves, Ai is the Airy function, and

The phases *ϕ*_{1} and *ϕ*_{2} are for a ray incident upon and reflected from the caustic, respectively.

The time dependence of the trapped wave solution (A3) comes in through *S _{M}*, defined by

where

The counting variable *M* indicates the number of pairs of incident and reflected rays. Initially *M* = 0 with *S*_{0} defined to be zero. Each time the initial ray makes a complete round trip from the ground to the turning point and back to the ground, *M* increases by 1. We compute *M* as a function of time using Eq. (33) of Broutman et al. (2006).

For the vertically trapped waves, the longtime steady-state solution is

where *r*_{0} is *r* at *z* = 0. This is Eq. (A9) of Broutman et al. (2003). It is used here in section 6b for the calculation of resonant modes, which occur where Ai(*r*_{0}) = 0.

## Footnotes

*Corresponding author address:* Stephen D. Eckermann, E. O. Hulburt Center for Space Research, Naval Research Laboratory, Code 7646, Washington, DC 20375. Email: stephen.eckermann@nrl.navy.mil