15 NOVEMBER 1999

SATOMURA AND SATO

3847

Secondary Generation of Gravity Waves Associated with the Breaking ofMountain Waves
TAKEHIKO SATOMURA AND KAORU SATO
Division of Geophysics, Graduate School of Science, Kyoto University, Kyoto, Japan(Manuscript received 24 June 1997, in final form 21 January 1999)
ABSTRACTThe generation of small-scale gravity waves associated with the breaking of mountain waves in the stratosphere
has been simulated within a fully compressible, nonhydrostatic, two-dimensional numerical model. The model
includes the mesosphere and offers sufficiently high spatial resolution to characterize the breaking and generationof gravity waves in the stratosphere. The mean flow is initialized by CIRA86 at 408N in February. A bell-shaped
mountain with 30-km half-width and 1.5-km height is located in the bottom of the modeled domain.
The primary wave forced by the mountain propagates into the stratosphere with amplitude increasing withheight and subsequently breaks in the lower stratosphere. After the primary wave breaking, significant wave
activity is simulated in the stratosphere. These secondary waves are generated at both the upstream and thedownstream edges of the breaking zone. Analysis is mostly focused on downstream small-scale waves. The
horizontal and vertical wavelengths of the secondary gravity waves are 3–8 km and 3–20 km, respectively, andphase velocities are 21.5 to 14 m s21 . The amplitudes of the vertical velocity component of the secondary
waves are 0.1–0.2 m s21 at altitudes of about 20 km. Theoretical consideration and model simulation suggestthat the winter stratosphere can be a wave duct for small-scale gravity waves once generated there as in thepresent simulation, because zonal winds are minimized in the winter stratosphere between the height regions of
the subtropical jet and the mesospheric jet. It is also suggested that both the convective instability and an
instability related with normal-modes act separately in different areas to generate the secondary gravity waves.

1. Introduction
The amplitude of internal gravity waves propagating
upward increases with height owing to the decreasing
air density. The gravity waves break in high altitudes
because statically unstable areas are locally produced,
if no critical level is met before the wave reaches there.
Today it is believed that the momentum deposition associated with the breaking of gravity waves is a significant forcing of the mean flow and causes the observed weak mean zonal winds near the mesopause (see,e.g., Holton 1982; Matsuno 1982) and the lower stratosphere (e.g., Tanaka and Yamanaka 1985; Palmer et al.
1986; McFarlane 1987). Further, gravity waves are also
considered to play an important role in the quasi-biennial oscillation observed in the equatorial stratosphere[Dunkerton (1997a) for a review].
What, then, is the source of gravity waves in the
middle atmosphere and how do they propagate from thesource? From observational studies using aircraft (e.g.,
Fritts and Nastrom 1992) or mesosphere, stratosphere,

Corresponding author address: Dr. T. Satomura, Division of Geophysics, Graduate School of Science, Kyoto University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan.
E-mail: satomura@kugi.kyoto-u.ac.jp

q 1999 American Meteorological Society

and troposphere radar (e.g., Sato 1994), the dominantsource of gravity waves is believed to be topography in
midlatitudes. To understand topographic gravity waves,
extensive field experiments have been completed overthe Rocky Mountains (e.g., Lilly and Kennedy 1973),
around the Alps (Kuettner 1986), and the Pyrenees(Bougeault et al. 1993). In general, however, observations by sondes, radars, or lidars are fixed in horizontal
position and the observational tracks of aircraft offer
sparse coverage in the vertical. Owing to these constraints it is difficult to study stationary waves such as
mountain waves in detail by using the observational data
only. Therefore, studies using numerical models andobservational studies are complementary to each other.
While the detailed propagation processes of mountainwaves from the ground to the middle atmosphere are
not yet clearly understood, studies to introduce the effects of orographic gravity waves into the middle atmospheric general circulation model have been carriedout since the 1980s. Lindzen (1981) hypothesized thatthe amplitudes of gravity waves would be ‘‘saturated’’owing to local instability acting to prevent further
growth of amplitudes with height. The obliterated wave
momentum due to the gravity wave saturation is assumed to be deposited in background flow. Based on
this assumption, the momentum deposition through the3848

JOURNAL OF THE ATMOSPHERIC SCIENCES

breaking of stationary waves has been parameterized inmany global and regional numerical models. These attempts have succeeded in reducing the westerly bias ofmodel results in the lower stratosphere and the troposphere (e.g., Palmer et al. 1986; McFarlane 1987; Iwasaki 1989). Recently, nonstationary gravity waves havealso been considered important for precise modeling of
the middle atmospheric circulation and new hypotheses
and parameterizations have been proposed (e.g., Hines
1997a,b). These studies of gravity wave drag parameterizations have focused on momentum deposited bywave breaking; they ignored detailed processes of wave
breaking and did not trace propagation of large-amplitude waves from the troposphere to the middle atmosphere for their practical purposes.
The breaking of gravity waves propagating from below has been studied by two-dimensional numericalmodels. Bacmeister and Schoeberl (1989) simulated thebreaking of mountain waves in a deep (0–85 km) atmosphere using anelastic equations with constant basic
horizontal wind and Brunt–Väisälä frequency. Theirsimulation indicated that a nonlinear wave–wave interaction produced large-amplitude gravity waves emitted
downward from the breaking zone of the mountainwave. The downward gravity waves have almost thesame horizontal wavelengths as the mountain wave.
These downward gravity waves greatly reduced the momentum flux due to mountain waves below the overturning level. Walterscheid and Schubert (1990) studied
the breaking processes of a long (300-km wavelength)monochromatic gravity wave in an isothermal middle
atmosphere at rest. They used a compressible numericalmodel and concluded that the amplitude growth of the
wave with height that is supposed to occur due to density
decrease is substantially reduced because of early breakdown and conversion of kinetic energy into potential
energy. Prusa et al. (1996) recently employed an anelastic numerical model and simulated the propagationand the breaking of gravity waves excited by an unsteady Gaussian obstacle located at an altitude of 15 km(the bottom of the model domain). It was shown that amonochromatic wave dominates the wave field at highaltitudes because of the dispersion of gravity waves during vertical propagation and that the breaking levels
lowered in time from z ø 100 km to 60 km. They
initialized the whole model domain by constant horizontal wind and Brunt–Väisälä frequency as did Bacmeister and Schoeberl (1989). Owing to the simplified
basic state and the small obstacle located at a high altitude of 15 km, the gravity waves propagate from thestratosphere to the mesosphere without breaking. The
processes of propagation and the breaking of mountainwaves under realistic atmospheric conditions have notbeen examined in depth.
The purpose of the present work is to simulate breaking of mountain waves under realistic atmospheric conditions. It is expected that gravity waves are generated
through the process of mountain wave breaking. The

VOLUME 56

characteristics of the secondary gravity waves would bedetermined by the breaking processes and backgroundflow. The excitation of internal gravity waves by instability of parallel flows was studied theoretically and also
numerically (e.g., McIntyre and Weissman 1978; Fritts
1982; Sutherland and Peltier 1992; Sutherland et al.
1994; Dunkerton 1997b). The basic conditions of thesestudies are, however, confined to simple parallel flows.
Dunkerton and Robins (1992) found a new mode ofsecondary instability ‘‘radiating’’ from an overturning
region of primary monochromatic gravity waves. The
amplitude of this radiating mode decreases so rapidly
from the overturning zone that its vertical extent is about
one vertical wavelength of the primary gravity wave orless.
To simulate propagation of mountain waves througha large vertical extent from the surface to the mesopause
correctly, a fully compressible numerical model is employed. The horizontal and the vertical resolutions aremaintained sufficiently high to resolve the breaking andassociated processes.
2. Model description and experimental setup
The numerical model we use is based on the twodimensional fully compressible model developed by Satomura (1989). While this model is able to generate
terrain-following coordinates numerically, it incorporates terrain at the lower boundary through the use ofa hybrid coordinates of the Cartesian coordinate and aconventional coordinate transformation given by GalChen and Somerville (1975) asz (z 2 h)z* 5 1for z # z1
z1 2 hfor z . z1 ,
z* 5 z(1)where h(x) is the lower boundary height of the integration domain and z 1 5 15 km. The compressible equations of motion, the continuity equation, and the thermodynamic equation aredvu
1 cp (Q 1 u)=p 5 kg 1 Fv ,dt
Qd(P 1 p)R (P 1 p)
R (P 1 p)
1 d=·v 5 dF,dtcycy (Q 1 u) u

(2)
(3)

and
du
5 Fu ,
dt

(4)

where
P(z) 1 p (x, z, t) 5

12 ,
p
p0

R d /cp

dP(z)
g
52
,
dz
cpQ(z)

15 NOVEMBER 1999

SATOMURA AND SATO

and v 5 (u, w). The right-hand sides of (2)–(4) are thecontributions from turbulent mixing, for which we employ a conventional first-order closure formulation (Lilly 1962; Satomura 1989).
The model uses a second-order finite-difference approximation centered both in time and in space. To gain
computational efficiency, both a time splitting method
similar to Klemp and Wilhelmson (1978) and an implicit
method in the vertical direction are applied to soundwave terms in the compressible equations. A fourthorder artificial viscous term is added to Eqs. (2) and (4)to damp small-scale computational modes.
The 512-km-wide integration domain is covered bygrids with 500-m spacing in the x direction. In the vertical direction, the domain height is 97 km and is divided
into 600 layers whose depth increases smoothly withheight from 25 m near the ground to 500 m in themesosphere. In the stratosphere and most of the middle
and the upper troposphere, which are focused on in thisstudy, the vertical grid interval is fixed at 100 m in orderto keep sufficient resolution.
At the lateral boundaries, an open boundary condition(Miller and Thorpe 1981) based on Orlanski (1976) is
employed. The upper boundary is a rigid wall and a
thick ‘‘sponge’’ layer is put above a height of 50 km
to avoid reflection of vertically propagating waves at
the upper boundary. The lower boundary is a free-slip
and thermally insulated rigid wall. A bell-shaped mountaina 2 h0,a 1 (x 2 x 0 ) 2

h(x) 52(5)where h 0 5 1.5 km and a 5 30 km, is located at thecenter (x 0 5 256 km) of the lower boundary.
The initial wind and temperature fields are horizontally uniform. The vertical profiles of the initial wind
and temperature are shown in Fig. 1, which are takenfrom the COSPAR International Reference Atmosphere
86 at 408N in February and interpolated to the model
grid points by cubic-spline functions. In this profile, thetropopause and the stratopause are located at altitudes
of about 10 and 47 km, respectively. The wind has local
maxima of 29 m s21 at 13-km height corresponding to
the subtropical jet and 68 m s21 at 61-km height (themesospheric jet). The vertical profile of the critical horizontal wavelength

lc [

2p
,
l

(6)

where l is the Scorer parameter
l2 5

N2
1 d2u
2
,
u2
u dz 2(7)

and N the Brunt–Väisälä frequency, is also shown in
Fig. 1. The critical horizontal wavelength has a local
maximum in the upper troposphere. The linear theory
indicates that mountain waves whose horizontal wave-

3849

FIG. 1. Vertical profiles of initial horizontal wind (solid), temperature (dotted), and the critical wavelength (dashed) taken from
CIRA86 at 408N on February.

lengths l x are shorter than l c are vertically evanescent(e.g., Smith 1977) and reflected downward. It should
be noted that gravity waves with l x , 10 km generated
by mountains cannot propagate upward into the stratosphere in the present basic state.
The simulation was run for 10 h and diagnostic field
variables were saved at 5-min intervals between t 5 6
h and 10 h for the analysis. A contour plot of horizontally averaged momentum flux, u9w9, where u9 and w9are wave components and the overbar indicates horizontal average, is shown in Fig. 2. This height–timediagram of averaged momentum flux indicates that themountain waves were well established up to the stratosphere after 6-h integration and had reached a practicallysteady state at 10 h.
3. Results
Snapshots at 5 h, 7 h, 9 h, and 10 h around the mountain are shown in Fig. 3. Blue and red contours represent
downward and upward vertical velocities, respectively,
and the black lines are contours of potential temperature.
The primary mountain wave propagates upward fromthe mountain to the stratosphere. The contours of potential temperature in Figs. 3b–d clearly show that thewave breaking occurs in the region of x 5 250–300 km
and z 5 21–24 km. While the initial wind is weak atthese heights (see Fig. 1), it is still more than 10 m s21 ;no critical level exists in the initial wind there. Thus itis probable that the increase of wave amplitude due tothe decrease of air density with height results in this3850

JOURNAL OF THE ATMOSPHERIC SCIENCES

VOLUME 56

horizontally. A close-up view of filtered vertical velocity
at t 5 9 h near the tail of the breaking zone is shownin Fig. 4. We can find wave activity both above (regionA) and below (region B) the breaking zone indicatedby a border line between the regions A and B. The
maximum amplitudes of the vertical velocities of the
waves are 0.1–0.2 m s21 . The waves in region A transmit
energy upward and those in region B are evanescent asshown in section 3b. Waves in both regions A and B
seem to originate from a thin sheet along their mutual
border. We can also identify waves transmitting energy
downward (region C) and upward (region D) from asmall area around (x, z) ø (313 km, 20 km).
a. Phase velocity and wavelength

FIG. 2. Vertical momentum flux, u9w9, as a function of z and t.
The contour interval is 2.5 3 10 4 N m21 .

wave breaking in the lower stratosphere. The group velocity of linear steady gravity waves in the lower stratosphere is more than several m s21 for l x less than several
tens of kilometers. Therefore, Fig. 3 shows that a partof the wave energy of the primary mountain wave istransmitted upward through the breaking zone.
Significant wave activity besides the mountain wave
are found in two regions in Figs. 3c and 3d: downstreamof the mountain throughout the troposphere and stratosphere, and over the mountain in the stratosphere. Waves
over the mountain have longer wavelengths than thedownstream waves. The phase lines of the waves overthe mountain were nearly parallel to the mountain wave
phase lines near the breaking zone, which is clear inclose-up views (not shown). The constant phase linesof the downstream waves are not parallel to the mountain wave phases. The difference in phases and wavelengths suggests that the wave excitation mechanisms
are different between the two regions. We will focus on
the downstream short waves in the following sections,because they are easily separated from the primary
mountain waves.
The short waves seem to propagate both upward (x5 310–330 km, z 5 20–40 km in Fig. 3c and x 5 320–
340 km, z 5 20–40 km in Fig. 3d) into the stratosphereand downward (x 5 310–320 km, z 5 5–18 km in Figs.
3c and 3d) into the troposphere from the breaking zone.
The waves propagating upward have shorter wavelengths than those propagating downward. There areother small-scale and small-amplitude waves around thebreaking zone. We will examine these waves in detailbelow.
To analyze the structure of the short waves, a highpass filter whose cutoff wavelength is 10 km was applied

Figure 5 shows an x–t diagram of filtered verticalvelocities at z 5 27 km (regions A and D) and z 5 15
km (regions B and C). It should be noted that waves
have small but nonzero horizontal phase velocities. Further, waves in region A have a distinct negative phasevelocity. The wave parameters and environmental variables are summarized in Table 1 where horizontal andvertical wavelengths, l x and l z , are estimated by assuming plane waves in each region. The mean values
and variations of l x and l z given in Table 1 are calculated from 10 samples around the selected points foru and N (see below for point selection) in each region.
The lengths of each sample are 40 km and 18 km forthe estimation of horizontal and vertical wavelength,
respectively. The vertical wavelength of linear gravity
wave with horizontal wavelength of l x can be also calculated from the linear dispersion relation of internal
gravity waves:

[1

2 1 2

N
l zcp 5 2p
u 2 cp

2

2p
2
lx

2

]

1
2
4H 2

21/2

,

(8)

where H is the density-scale height. We selected points
for the estimation of u and N as listed in Table 1, where
spatial variations of u and N in both x and z directions
are rather small and neglected the effect of d 2 u/dz 2 . The
density scale height H decreases slightly from 6.6 km
at z 5 15 km to 6.2 km at z 5 27 km.
The theoretically estimated vertical wavelength, l zcp ,corresponds well with the observed wavelength l z , indicating that waves in regions A, C, and D are likelyinternal gravity waves. These waves have similar horizontal wavelengths of between 3 and 8 km and horizontal phase velocities within 21.5 to 12 m s21 , whereas vertical wavelengths vary from 3 km in region D to
30 km in C.
2
For waves in region B, however, l zcp
is negative and,thus, waves are evanescent (or external) in the vertical.
In fact, the amplitudes of waves in region B clearlydecrease before the waves reach z 5 15 km (Fig. 4),
although the altitude of source area of waves in region
B is almost the same as that in region C.

15 NOVEMBER 1999

SATOMURA AND SATO

3851

FIG. 3. Vertical velocity, w, in m s21 and potential temperature at (a) 5 h, (b) 7 h, (c) 9 h, and (d) 10 h of simulation time. The blue andred lines are contours of negative and positive vertical velocity, respectively. The contour intervals of w are 60.1, 60.2, 60.5, and 61 m
s21 . The black lines are contours of potential temperature with an interval of 20 K. The mountain top is located at x 5 256 km.

b. Momentum flux
The waves found in the regions A–D have the horizontal phase velocities less than the background horizontal wind velocity. If internal gravity waves have negative phase speed with respect to the mean flow as thewaves in the four regions, the momentum flux ru9w9
should be negative for upward energy-propagatingwaves and positive for downward propagating waves,
where dashed variables indicate wave components.

Figure 6 shows the time series of the momentum flux
ru0w0 starting from t 5 8 h 30 min with a 30-min
interval, where u0 and w0 are filtered velocities. It isclear that upward propagating waves in regions A andD have negative momentum fluxes and downward propagating waves in region C have positive momentum
flux: the signs of the momentum fluxes are consistent
with the internal gravity wave characteristics mentioned
at the beginning of this section. The maximum absolutevalues of these fluxes are greater than 5 3 1023 N m22 .

3852

JOURNAL OF THE ATMOSPHERIC SCIENCES

VOLUME 56

FIG. 4. Filtered vertical velocity at t 5 9 h near the downstream
tail of the breaking zone. The contour interval is 0.04 m s 21 . Dashedcontours represent negative values. Thick dash–dotted lines separate
regions A, B, C, and D.

Evanescent waves in region B, however, again show
different characteristics. Here the momentum flux alternates between positive and negative signs with almost
the same amplitude. The horizontally averaged momentum flux in region B is nearly equal to zero.
Upward and downward propagating waves are reflected and the sign of the momentum flux is changed
at z ø 35 km in region D and z ø 15 km in region C,
respectively. This reflection of waves will be discussed
in the next section.
c. Spectral characteristics
Before we start discussion on the waves found aroundthe breaking zone of the mountain wave in detail, it isworth confirming that waves in the regions A–D are notpropagating from the troposphere but are generated in
the stratosphere. The power spectrum of vertical velocity at t 5 10 h and the curve of the upper height limit
where wave packets started at the ground at t 5 0 canreach at t 5 10 h are shown in Fig. 7. The upper height
limit is calculated by integrating the group velocity ofeach wavenumber estimated from the linear dispersion
relation using the initial horizontal velocity and potential temperature profiles. Waves with horizontal wavelengths shorter than 10 km cannot propagate from theground above z 5 15 km where the critical wavelength
is maximized (see Fig. 1). The upper limit of propagation for waves with longer wavelengths is lower because of their smaller vertical group velocities. The
power spectra are calculated only for z . 15 km
because the model coordinates are curvilinear below 15
km and hence the interpretation of the Fourier transform
within such coordinates is difficult there. The line of
the upper limit coincides well with the contour of 1025
kg s22, particularly for z . 20 km in the wavelength
range longer than about 10 km corresponding to the
maximum critical wavelength at z 5 15 km. Thus theline of upper limit could be used as an indicator of the

FIG. 5. Horizontal-time section diagram of filtered vertical velocity
at (a) z 5 27 km and (b) z 5 15 km. The contour interval is (a) 0.02
m s21 and (b) 0.01 m s21 . Dashed contours represent negative values.
Thick solid lines show horizontal phase velocity c p 5 2 m s21 . Dash–dotted lines show (a) c p 5 21.5 m s21 and (b) c p 5 4 m s21 .

upper-limit altitude reached by the primary mountainwave.
Substantial power penetrates into the ‘‘inhibited’’ areaon the shortwave side (l x , 10 km and z . 15 km) in
Fig. 7, corresponding to waves whose wavelengths areshorter than 10 km observed around the breaking zone.
Therefore, the short waves around the breaking zone
TABLE 1. Wave parameters and environmental variables at 9 h.
Here lx is horizontal wavelength, u is unfiltered horizontal velocity,
lzcp is vertical wavelength calculated from the dispersion relation, —is external mode, cp is horizontal phase velocity, and lz is vertical
wavelength.
N
Region lx (km) cp (m s21) lz (km) lzcp (km) u (m s21) (31022 s21)
A
B
C
D

3.3–3.7
5.0–5.7
5.5–8.5
5.0–1221.5
4
2
2

2.3–9.0 2.8–7.6 5.3–9.3a——
26–30b
12–33 8.0–41 17–26c
2.5–2.8 1.4–3.3 6.4–12d

At (x, z) 5 (290–310 km, 24–25 km).
At (x, z) 5 (310–330 km, 17–19 km).
c
At (x, z) 5 (310–330 km, 17–19 km).
d
At (x, z) 5 (290–310 km, 24–25 km).
ab

1.9–2.2
2.1–2.3
2.0–2.1
1.9–2.2

15 NOVEMBER 1999

SATOMURA AND SATO

are not those propagating from the bottom of the model,
but are, in fact, those generated in the stratosphere. In
this sense, it is appropriate to call these waves ‘‘secondary gravity waves’’ since they are generated by processes associated with the breaking of the primary
mountain waves.
Note, also, that the thermal and wind structure of thewinter stratosphere forms a wave duct for secondary
gravity waves. As shown in a region of x 5 320–350
km and z 5 25–35 km in Figs. 6c and 6d, the energyof the secondary waves cannot propagate into the mesosphere but is reflected at z 5 35–40 km since thewavelengths of secondary waves are less than 10 km
and the critical wavelength l c surpasses 10 km once
more in the upper stratosphere (Fig. 1). Owing to theirnonzero phase velocities, the secondary gravity waves
in regions C and D will propagate away downstream
within the stratospheric duct from the breaking zone ofthe primary mountain wave. Thus, these secondary
gravity waves will possibly disturb the stratospheric atmosphere even over the ocean or over flat surfaces.
4. Discussion
a. Generation mechanisms
In our model, gravity waves with different wavelengths and phase velocities and vertically evanescentdisturbances were generated in different areas. To look
into the generation mechanisms of each region definedin Fig. 4, greatly magnified views of the Brunt–Väisälä
frequency, horizontal velocity, and log 2Ri, where Ri is
the Richardson number, near the downstream tail of the
breaking zone at 9 h, are shown in Fig. 8 with filteredvertical velocity. At 9 h, the mountain waves havereached steady development for several hours as shownin Fig. 2.
It is clear that gravity waves in region A originate in
the middle or upper part of the wave-breaking zone of
the mountain wave where the vertical stratification isunstable (Fig. 8a). Strong convection is expected to prevail under such unstable conditions. Given an unstablelayer of depth 1 km, the linear theory of convective
instability between two rigid stress-free plates predicts
that the horizontal wavelength of the fastest growing
mode is 2.8 km (e.g., Chandrasekhar 1961); this value
is similar to that of the gravity waves observed in region
A. The horizontally averaged heat flux associated withthe filtered components is clearly positive from z 5 20.1
km to 21.2 km corresponding to the breaking zone (notshown).
Convective motions are considered to generate gravity waves by both mechanical oscillator and obstaclerelated mechanisms (e.g., Clark et al. 1986; Fovell etal. 1992). Mechanical oscillators would generate gravity
waves that propagate to right and left, whereas gravity
waves generated by the obstacle effect propagate withthe same horizontal phase velocity as the convective

3853

movement. The phase velocity of the gravity waves inregion A (21.5 m s21 , see Fig. 5a) is almost the sameas the horizontal wind velocity in the overturning zone(ø21 m s21 , see Fig. 8b). Considering that the horizontal wind speed gives a good estimation of the advection speed of the embedded convection, the obstacle
effect may be the generation mechanism of the secondary waves in region A.
The generation of secondary gravity waves in regions
B–D may be explained in terms of the instability studied
by Sutherland and Peltier (1992, referred to as SP92
hereafter) and Dunkerton (1997b). SP92 found unstable
but propagating modes by a linear stability analysis ofhorizontally uniform flow. They used a vertical profile
of horizontal velocity and Brunt–Väisälä frequency extracted from one wavelength of Long’s solution (Long1953) above a critical position in the lee of the obstacle.
The origin of disturbances in region B is located ina stable zone just beneath the breaking point. This stablezone overlaps a strong shear zone (Fig. 8), which is 2
km thick at x 5 305 km and in which the horizontal
wind decreases from 30 m s21 at z 5 19 km to a slowly
reversed wind at z 5 21 km. The wavenumber of thetemporally fastest growing mode (6.546 in dimensionless units, see SP92) corresponds to about 7.5 km if the
thickness and strength of shear are 2 km and 1.5 3 1022
s21 , respectively, as in the strong shear zone in region
B. The most unstable mode analyzed by SP92 has amaximum amplitude just below the neutral buoyancy
layer. Therefore the characteristics of the structure andwavelength of the most unstable mode analyzed by SP92
also agrees well with disturbances in region B.
The phase velocity of the most unstable mode of SP92
is about 0.2 m s21 , while the phase velocity of waves
in region B is 4 m s21 . It should be noticed that thelevels of the minimum horizontal velocity and the neutral stratification coincide in SP92 whereas the primary
mountain wave in our simulation developed into a nonlinear state much further than the condition of SP92. In
fact, looking into Figs. 8a and 8b, the neutral stratification is formed 500 m below the level of the minimum
horizontal velocity at 9 h. If we assume that the phasevelocity of the most unstable mode is almost stationary
relative to the horizontal velocity at the level of neutral
stratification, it then becomes about 3 m s21 . This value
corresponds well with the phase velocity of the disturbances in region B.
In region B, the local critical wavelength1 l c 5 6–7.5
km at z 5 18–19 km owing to the stronger horizontal
wind in region B. The horizontal wavelength of waves
is l x ø 5.4 km in region B. Therefore the propagation
condition l x . l c is not satisfied and waves are vertically evanescent in region B.

1
Here l c is calculated by Eqs. (6) and (7), which are not exact forwaves with nonzero phase velocities. The error is, however, only
about 10% because the background wind is strong (ø30 m s21 ).

3854

JOURNAL OF THE ATMOSPHERIC SCIENCES

VOLUME 56

FIG. 6. Horizontal-height section of momentum flux ru0w0 (N m22 ) at (a) 8 h 30 min, (b) 9 h, (c) 9 h 30 min, and (d) 10 h. The
momentum flux is calculated using filtered wind velocities.

The source area of waves in region C seems to be(x, z) ø (312 km, 21 km) at 9 h, near the downstream
end of the breaking zone. The thickness of the vertical
shear zone at x 5 312 km is about 3 km (from u 5 28
m s21 at z 5 18 km to u 5 1 m s21 at z 5 21 km). The
Brunt–Väisälä frequency decreases to almost zero at z
5 21 km where u 5 1 m s21 (Figs. 8a and 8b). The
most unstable mode analyzed by SP92 for this background state has about a 10-km wavelength and about

1 m s21 phase velocity. These values agree well withthe parameters of the gravity waves in region C and theinstability possibly causes the gravity waves in this region. The horizontal length of the source area, however,extends over about one wavelength only. Thus, it may
be difficult to apply the analysis of SP92, which assumed a horizontally uniform wind and temperature.
It is possible that the gravity wave generation in region C is controlled by spatial structure of the unstable15 NOVEMBER 1999

SATOMURA AND SATO

3855

FIG. 7. Power spectrum of vertical velocity weighted by air density
at t 5 10 h. The contours are 1023 , 1025 , 1027 kg s22 . The shaded
area represents values .1025 kg s22 . The thick solid curve indicatesthe upper limit of the height where wave packets starting from theground at t 5 0 h can reach at 10 h.

region. Scinocca and Peltier (1993) observed a largeamplitude stationary disturbance similar to the most unstable mode of SP92 in their nonlinear numerical simulation of downslope windstorms. The wavelength ofthe disturbance agrees with the most unstable mode ofSP92 and the source area of the disturbance is locatedat the downstream edge of the breaking zone: the characteristics of the large-amplitude stationary disturbanceobserved by Scinocca and Peltier (1993) are very similar
to those of the waves in region C.
Figure 9 shows close-up views of x–t diagrams offiltered vertical velocities near the downstream end of
the breaking zone. The left edge of the wave packet inregion C (Fig. 9b), which is considered to be its source
area, moves downstream from (x, t) ø (302 km, 7.4 h)to (x, t) ø (310 km, 8.4 h). This movement corresponds
closely with the stretch of the mechanically unstable (Ri, 0.25) layer in both their locations and their speeds.
As seen clearly in Fig. 9b, the phase velocity of gravity
waves coincides with the stretching speed of the unstable layer. Since this unstable layer is associated with theprimary wave breaking, the correspondence between the
wave characteristics and the stretch of the unstable layer
suggests that the gravity waves in region C are generated
by a similar mechanism as shown by Scinocca and Peltier (1993), that is, through a nonlinear process at theedge of the breaking zone.
The local critical wavelength l c in region C is 5.5–
7.0 km. The horizontal wavelength of the secondary
wave l x ø 7.6 km in region C. Therefore the propa-

FIG. 8. Close-up of (a) Brunt–Väisälä frequency, N; (b) horizontal
velocity, u; and (c) log 2Ri, where Ri is the Richardson number, at t
5 9 h near the tail of the breaking zone. The contour interval in (a)is 0.0025 s21 , (b) 2 m s21 , and (c) 1. A contour line of u 5 21 m
s21 is added in (b). Shaded areas indicate (a) N , 0, (b) u , 0, and(c) Ri , 0.25. Filtered vertical velocity is also plotted in blue (negative) and red (positive) contours for 60.05 m s21 , 60.1 m s21 , and60.2 m s21 . Regions are indicated by letters A–D.

gation condition l x . l c is satisfied there and gravity
waves can propagate far from the source area.
Gravity waves in region D appear from (x, z) ø (314km, 21 km) at t 5 9 h, which is also near the downstreamedge of the breaking zone. Judged from the fact thatthe horizontal wavelength and phase velocity of gravity
waves in region D are very similar to those in region
C, it is inferred that gravity waves in region D are generated by the same generation mechanism in region C.
It is worth noting that the leftmost wave in region
D appears from almost the same location x ø 312 km

3856

JOURNAL OF THE ATMOSPHERIC SCIENCES

VOLUME 56

gravity waves in regions C and D be due to the difference of the background wind field.
b. Reflection of generated waves
If the critical wavelength l c increases slowly in the
direction of the wave propagation, linear theory predicts
that plain stationary gravity waves are reflected at aheight where l c 5 l x . Assuming that gravity waves arealmost stationary in regions C and D, the heights wherel c 5 l x are about 31 and 16 km, respectively (Fig. 1
and Table 1). In fact the gravity waves are reflected atthese altitudes, as shown in Fig. 6. Departure from zero
phase velocity of about 2 m s21 modifies the reflection
altitudes only slightly to z ø 35 km for region C andz ø 15 km for region D. The wave reflection in theupper stratosphere is also observed in the power spectra(Fig. 7) where strong and weak powers in w appear
alternately in the vertical. The winter stratosphere acts
as a wave duct for small-scale gravity waves and suchtrapped gravity waves can be generated in association
with breaking of large-scale gravity waves propagatingfrom the lower troposphere as shown in this study.
c. Observational evidence

FIG. 9. Horizontal-time section diagram of filtered vertical velocity(thin lines) at (a) z 5 22 km and (b) z 5 19 km. The contour intervalis (a) 0.02 m s21 and (b) 0.04 m s21 . Dashed contours represent
negative values. Areas where the minimum Ri between z 5 19 km
and 22 km is less than 0.25 are shaded. Thick solid contours in (a)denote averaged intrinsic phase velocities u 2 c p 5 0, 1, 2, 3, and4 m s21 .

for t 5 8–8.5 h (Fig. 9a), whereas gravity waves inregion C appear along the edge of the breaking zone
propagating downstream (Fig. 9b). The waves in region
D are observed within an area where the intrinsic phasevelocity is rather larger than the horizontal wind. Thicksolid curves in Fig. 9a show contours of averaged intrinsic phase velocitiy u 2 c p in the height range of 19–
22 km:
u 2 cp 5

E

z522 km(u 2 cp ) dz,

z519 kmwhere c p 5 2 m s21 as in Table 1. The value of u 2 c p
is smaller than 1 m s21 for x , 310 km. Although the
WKB approach is not applicable to the region where
the wind varies with similar scale of gravity waves, it
is natural to consider that the existence of the area where
the intrinsic phase velocity becomes nearly equal to zeroprevents secondary gravity waves from propagating
freely though it to region D. Therefore it is deduced
that the difference in the timing of appearance between

The secondary gravity waves simulated in this study
have short horizontal wavelengths l x 5 3–8 km, long
vertical wavelengths l z 5 3–20 km, and small phasevelocities c p 5 21.5 to 14 m s21 . Among the observational studies reporting strong vertical wind fluctuations in the middle atmosphere, Sato and Hirota (1988)
first succeeded in determining directly horizontal structures of gravity waves with the MU radar (Shigaraki,
Japan; 358N, 1368E) by using its multibeam capability
to observe the winter lower stratosphere. The determined wave parameters were l x , 5–30 km, l z ø 5km, and c p 5 1–5 m s21 in the height range of 12–22
km; these parameters are very similar to the secondary
waves found in this study.
A subtropical westerly jet (ø80 m s21 ) was locatedat an altitude of around 12 km during the observation
period of Sato and Hirota (1988). The small-scale gravity waves as observed by them cannot propagate fromthe ground through the jet to the lower stratosphere. In
the winter, a strong monsoon wind (west-north-westerlywind) blows over Japan and mountain waves develop
frequently. On the windward side of the MU radar in
winter, there are mountains of about 1 km in height.
Therefore, it is reasonable to consider that the gravity
waves detected by Sato and Hirota (1988) were the secondary gravity waves generated in association withbreaking of mountain waves as observed in our simulation. Sato (1990) showed a high correlation between
the activity of vertical wind disturbances observed by
the MU radar in winter and the horizontal wind at 90
kPa (ø1 km) at Yonago, which is about 250 km westnorthwest of the radar site and on the upstream side of15 NOVEMBER 1999

3857

SATOMURA AND SATO

the mountain. This fact also suggests that mountainwaves control the small-scale wave activity in the downstream area.
5. Conclusions
Using a two-dimensional nonhydrostatic model, we
have simulated the generation of new secondary gravity
waves associated with the breaking of primary mountainwaves in the lower stratosphere. The wave parametersof the secondary waves are 3–8 and 2–20 km for horizontal and vertical wavelengths, respectively, and 21.5
to 14 m s21 for phase velocities, which well satisfy thetheoretical dispersion relation of internal gravity waves.
The amplitudes of the vertical velocities of the secondary waves are 0.1–0.2 m s21 . These secondary gravity
waves propagate horizontally as well as vertically fromthe primary wave-breaking area because of nonzero
phase velocities. These waves are reflected verticallywhere vertical wavelengths become infinite and aretrapped in the stratosphere, because the vertical structureof the environmental wind and temperature in the winterlower stratosphere works as a duct for gravity waves
with such short horizontal wavelengths. The horizontal
wavelengths of the secondary gravity waves studied inthis paper are much shorter than the primary mountainwave (ø180 km). Hence, these waves are clearly different from the downward propagating waves found by
Bacmeister and Shoeberl (1989), whose horizontal
wavelength was almost the same as that of the primary
mountain wave.
Comparing the wave parameters with unstable modes
in linear and nonlinear analyses made in the previous
studies, three kinds of generation mechanisms are suggested for the secondary internal gravity waves: 1) an
obstacle effect due to convective motions in the shear
flow near the breaking zone (region A); 2) an unstablenormal mode of slowly growing breaking mountainwave (SP92; Scinocca and Peltier 1993) (region B); 3)
a mode growing nonlinearly at the edge of the breaking
mountain wave (regions C and D). Mountain wave simulations over narrower mountains (a 5 10 km and 5km) have also been performed (not shown). Secondary
gravity waves similar to those reported in this paper are
also generated in these cases. This fact suggests that thegeneration mechanisms of small-scale gravity waves as
discussed in this study act over a wide range of breaking
mountain wave parameters.
Recent numerical simulations indicate that disturbances with three-dimensional structures appear in welldeveloped K–H billows (e.g., Fritts et al. 1996). Since
the mechanisms of the secondary wave generation simulated in this paper are not the K–H instability, it isuncertain whether the three-dimensionality is essential
to generate the secondary gravity waves associated withmountain wave breaking in the real atmosphere. However, it is necessary to confirm our results by threedimensional models in the future. Time-dependent en-

vironmental flow is possibly of more relevance to mountain wave generation (e.g., Lott and Teitelbaum 1993)
and consequently to the generation mechanisms of thesecondary gravity waves. Thus, studies incorporating
the time evolution of the background flow should bealso made in order to fully understand the dynamical
effects of mountains on the atmosphere.
Acknowledgments. We acknowledge valuable discussions with Dr. David C. Fritts. The numerical experiment was done on the KDK system at Radio Atmospheric Science Center and on the Data Processing Center at Kyoto University. This research was supported by
the Center for Climate System Research, University ofTokyo; a Grant-in-Aid for Scientific Research B-208458146 of the Ministry of Education, Science andCulture, Japan (TS); and the Nissan Science Foundation(KS). Figures were drawn by the GFD-DENNOU Library.
REFERENCES
Bacmeister, J. T., and M. R. Schoeberl, 1989: Breakdown of vertically
propagating two-dimensional gravity waves forced by orography. J. Atmos. Sci., 46, 2109–2134.
Bougeault, P., and Coauthors, 1993: The atmospheric momentum budget over a major mountain range: First results of the PYREX
field program. Ann. Geophys., 11, 395–418.
Chandrasekhar, S., 1961: Hydrodynamic and Hydromagnetic Stability. Clarendon Press, 652 pp.
Clark, T. L., T. Hauf, and J. P. Kuettner, 1986: Convectively forced
internal gravity waves: Results from two-dimensional numericalexperiments. Quart. J. Roy. Meteor. Soc., 112, 899–925.
Dunkerton, T. J., 1997a: The role of gravity waves in the quasibiennial oscillation. J. Geophys. Res., 102, 26 053–26 076.
, 1997b: Shear instability of internal inertia-gravity waves. J.
Atmos. Sci., 54, 1628–1641.
, and R. E. Robins, 1992: Radiating and nonradiating modes ofsecondary instability in a gravity-wave critical layer. J. Atmos.
Sci., 49, 2546–2559.
Fovell, R., D. Durran, and J. R. Holton, 1992: Numerical simulationsof convectively generated stratospheric gravity waves. J. Atmos.
Sci., 49, 1427–1442.
Fritts, D. C., 1982: Shear excitation of atmospheric gravity waves.
J. Atmos. Sci., 39, 1936–1952.
, and G. D. Nastrom, 1992: Sources of mesoscale variability ofgravity waves. Part II: Frontal, convective, and jet stream excitation. J. Atmos. Sci., 49, 111–147.
, T. L. Palmer, Ø. Andreassen, and I. Lie, 1996: Evolution andbreakdown of Kelvin–Helmholtz billows in stratified compressible flows. Part I: Comparison of two- and three-dimensionalflows. J. Atmos. Sci., 53, 3173–3191.
Gal-Chen, T., and R. C. J. Somerville, 1975: On the use of a coordinate transformation for the solution of the Navier–Stokes equations. J. Comput. Phys., 17, 209–228.
Hines, C. O., 1997a: Doppler-spread parameterization of gravitywave momentum deposition in the middle atmosphere. Part 1:
Basic formulation J. Atmos. Solar-Terr. Phys., 59, 371–386.
, 1997b: Doppler-spread parameterization of gravity-wave momentum deposition in the middle atmosphere. Part 2: Broad andquasi monochromatic spectra, and implementation. J. Atmos. Solar-Terr. Phys., 59, 387–400.
Holton, J. R., 1982: The role of gravity wave induced drag and diffusion in the momentum budget of the mesopshere. J. Atmos.
Sci., 39, 791–799.
Iwasaki, T., 1989: A parameterization scheme of orographic gravity

3858

JOURNAL OF THE ATMOSPHERIC SCIENCES

wave drag with two different vertical partitionings. J. Meteor.
Soc. Japan, 67, 11–27.
Klemp, J. B., and R. B. Wilhelmson, 1978: The simulation of threedimensional convective storm dynamics. J. Atmos. Sci., 35,
1070–1096.
Kuettner, J. P., 1986: The Aim and Conduct of ALPEX. Vol. 27, GARP
Publ. Ser., WMO, 3–14 pp.
Lilly, D. K., 1962: On the numerical simulation of buoyant convection. Tellus, 14, 148–172.
, and P. J. Kennedy, 1973: Observations of a stationary mountainwave and its associated momentum flux and energy dissipation.
J. Atmos. Sci., 30, 1135–1152.
Lindzen, R. S., 1981: Turbulence and stress owing to gravity wave
and tidal breakdown. J. Geophys. Res., 86, 9707–9714.
Long, R. R., 1953: Some aspects of the flow of stratified fluids. Part
I: A theoretical investigation. Tellus, 5, 42–58.
Lott, F., and H. Teitelbaum, 1993: Linear unsteady mountain waves.
Tellus, 451, 201–220.
Matsuno, T., 1982: A quasi one-dimensional model of the middle
atmosphere circulation interacting with internal gravity waves.
J. Meteor. Soc. Japan, 60, 215–226.
McFarlane, N. A., 1987: The effect of orographically excited gravity
wave drag on the general circulation of the lower stratosphereand troposphere. J. Atmos. Sci., 44, 1775–1800.
McIntyre, M. E., and M. A. Weissman, 1978: On radiating instabilities
and resonant overreflection. J. Atmos. Sci., 35, 1190–1196.
Miller, M. J., and A. J. Thorpe, 1981: Radiation conditions for thelateral boundaries of limited-area numerical models. Quart. J.
Roy. Meteor. Soc., 107, 615–628.
Orlanski, I., 1976: A simple boundary condition for unbounded hyperbolic flows. J. Comput. Phys., 21, 251–269.
Palmer, T. N., G. J. Shutts, and R. Swinbank, 1986: Alleviation of asystematic westerly bias in general circulation and numerical

VOLUME 56

weather prediction models through an orographic gravity wave
drag parametrization. Quart. J. Roy. Meteor. Soc., 112, 1001–
1039.
Prusa, J. M., P. K. Smolarkiewicz, and R. R. Garcia, 1996: Propagation and breaking at high altitudes of gravity waves excitedby tropospheric forcing. J. Atmos. Sci., 53, 2186–2216.
Sato, K., 1990: Vertical wind disturbances in the troposphere andlower stratosphere observed by the MU radar. J. Atmos. Sci., 47,
2803–2817.
, 1994: A statistical study of the structure, saturation and sourcesof inertio-gravity waves in the lower stratosphere observed with
the MU radar. J. Atmos. Terr. Phys., 56, 755–774.
, and I. Hirota, 1988: Small-scale gravity waves in the lower
stratosphere revealed by the MU radar multi-beam observation.
J. Meteor. Soc. Japan, 66, 987–999.
Satomura, T., 1989: Compressible flow simulations on numericallygenerated grids. J. Meteor. Soc. Japan, 67, 473–482.
Scinocca, J. F., and W. R. Peltier, 1993: The instability of Long’s
stationary solution and the evolution toward severe downslope
windstorm flow. Part I: Nested grid numerical simulations. J.
Atmos. Sci., 50, 2245–2263.
Smith, R. B., 1977: The influence of mountains on the atmosphere.
Advances in Geophsyics, Vol. 21, Academic Press, 87–230.
Sutherland, B. R., and W. R. Peltier, 1992: The stability of stratified
jets. Geophys. Astrophys. Fluid Dyn., 66, 101–131.
, C. P. Caulfield, and W. R. Peltier, 1994: Internal gravity wave
generation and hydrodynamic instability. J. Atmos. Sci., 51,
3261–3280.
Tanaka, H., and M. D. Yamanaka, 1985: Atmospheric circulation in
the lower stratosphere induced by the mesoscale mountain wave
breakdown. J. Meteor. Soc. Japan, 63, 1047–1054.
Walterscheid, R. L., and G. Schubert, 1990: Nonlinear evolution ofan upward propagating gravity wave: Overturning, convection,transience and turbulence. J. Atmos. Sci., 47, 101–125.