INTERNAL GRAVITY WAVES: From Instabilities to Turbulence

■ Abstract We review the mechanisms of steepening and breaking for internal gravity waves in a continuous density stratification. After discussing the instability of a plane wave of arbitrary amplitude in an infinite medium at rest, we consider the steepening effects of wave reflection on a sloping boundary and propagation in a shear flow. The final process of breaking into small-scale turbulence is then presented. The influence of those processes upon the fluid medium by mean flow changes is discussed. The specific properties of wave turbulence, induced by wave-wave interactions and breaking, are illustrated by comparative studies of oceanic and atmospheric observations, as well as laboratory and numerical experiments. We then review the different attempts at a statistical description of internal gravity wave fields, whether weakly or strongly interacting.


INTRODUCTION
Internal (gravity) waves propagate in any fluid stably stratified with density.The stratification can be due to vertical temperature or concentration gradients (e.g., salinity in the oceans) with similar dynamical properties.In both cases, fluid particles oscillate owing to the restoring force of buoyancy.
Internal waves are observed in the atmosphere, at scales ranging from meters to kilometers, and dominate fluctuations in the stratosphere.The most common sources of internal gravity waves are the blowing of a wind over topography (see, e.g., Wurtele et al. 1996 for a review) and cumulus convective clouds.Other sources, about which much less is known, are the adjustment toward equilibrium of jets and of unsteady vortical flows (e.g., as observed by Vaughan & Worthington 2000) and the sudden formation of turbulent patches by a dynamical instability.Internal gravity waves can be revealed by cloud patterns (see Figure 1) and are measured by direct probing or remote sensing, such as radar and lidar.Star scintillation has been used to analyze internal waves in the stratosphere of Earth, and also on other planets, thanks to the observation of star occultation from space (Sicardi  et al. 1999).Internal waves are also expected to occur in stellar interior, being excited by perturbations from the outer convective zones (Schatzman 1996).
Internal waves in the ocean's interior are as common as waves at the sea surface.They have been reported since the beginning of the twentieth century as fluctuations of velocity and temperature, forming background "noise," which often dominates the mean currents.In some cases, the SAR (synthetic aperture radar), in being sensitive to the roughness of the sea surface, is able to provide "footprints" of oceanic internal waves (whose associated flow perturbs the sea surface).Internal waves are excited by the interaction of tide with topography, by wind stress fluctuations, and by various other processes that are still poorly documented.Owing to cascades of nonlinear interactions, it is often not possible to track back the origin of the observed waves.
Internal waves propagating along a density interface behave similarly to surface waves.We restrict this review to the quite different case of a continuously stratified fluid in which propagation occurs in a slanting direction.The restoring buoyancy force acting on a vertically displaced fluid particle is characterized by the Brunt-Väisälä (or buoyancy) frequency N defined (in an incompressible fluid) by where g denotes the acceleration of gravity, z the vertical (upward) coordinate, and ρ b (z) the density of the background state, assumed in hydrostatic balance, with horizontal isopycnals (constant density surfaces).Compressibility effects may be ignored in the ocean but not in the atmosphere (in the latter case, the vertical density gradient must be replaced by its difference from the adiabatic gradient).
We often consider a uniformly stratified fluid, i.e., with constant N, which is a good approximation when considering short wavelengths in comparison with the scale of density changes.The internal wave period ranges from the minimum value 2π/N , typically a fraction of an hour in the atmosphere or oceans, to much longer periods, typically a day, for nearly vertical wavevectors.In the latter case the effect of the Coriolis force is important, and internal waves are modified into inertia-gravity waves.We focus the review on pure internal gravity waves.
At wavelengths shorter than the atmospheric scale height, internal gravity waves are clearly distinguished from sound waves by their transverse (divergenceless) velocity field and much slower propagation.Their vertical velocity w obeys the following equation (obtained by eliminating all components but w from the Boussinesq equations linearized about the hydrostatic balanced state): where ∇ 2 H denotes the Laplacian operator in the horizontal plane.From Equation 2, the well-known dispersion relationship of internal gravity waves (in a medium at rest) is inferred as ω = ±N cosθ. (3) Equation 3 (first determined by Lord Rayleigh in 1883) relates the wave frequency ω to the wave vector k, whose angle with the horizontal is denoted θ .The velocity oscillations are perpendicular to k in the vertical plane.In the presence of a Coriolis force, Equation 3 is modified into ω 2 = N 2 cos 2 θ + f 2 sin 2 θ , where f is the Coriolis parameter and the polarization is elliptical instead of rectilinear in a vertical plane.The phase velocity ω/k is directed along k, whereas the group velocity c g = ∇ k ω is perpendicular to k, as the frequency ω does not depend upon the magnitude of the wavevector.Therefore the energy propagation and ray path (defined by ẋ = c g ) are remarkably perpendicular to k.We refer to Lighthill (1978), Gill (1985), or LeBlond & Mysak (1978) for classical reviews of the elementary internal wave properties.Equation 2 also admits the trivial solution w = 0: This corresponds to a purely horizontal motion, whose frequency is therefore zero.This motion is referred to as a vortex mode if it possesses vertical vorticity, but it may also contain a horizontal mean flow.It exists for each wavevector, in addition to the two wave modes (with frequency ±ω) characterized by the dispersion relation (Equation 3).Assuming each mode is nondivergent, this wave-vortex (or poloidal-toroidal) decomposition of the velocity field u can be written as u = ∇ × ψi z + (∇ h φ + wi z ) + ū(z, t). (4) The vortex mode ∇ × ψi z contains all the vertical vorticity, the wave mode ∇ h φ + wi z contains all the vertical velocity of the flow, and the last term is the horizontal mean flow.This decomposition (Equation 4) is mathematically possible for any incompressible velocity field, but its physical interpretation in terms of waves and vortices is valid only in the linear regime.Extensions to nonlinear regimes are reviewed by Riley & Lelong (2000).
Internal gravity waves can thus be considered as short-time and short-scale perturbations in natural media.For this reason, they are filtered out from weather forecast models, and their effect upon the largest scales must be modeled through parameterizations.The need for a suitable modeling is assessed by the global role of internal gravity waves at transporting momentum from the lower to the upper atmosphere (e.g., Holton et al. 1995, McIntyre 1995).In the ocean, the main role of internal gravity waves is a vertical (diapycnal) mixing of heat and solutants occuring during wave breaking.In all these processes, internal gravity waves transport energy and momentum from their sources, which they release during breaking events, thus inducing mean flow changes [linear nondissipative waves do not induce mean flow changes, as stated by the classical nonacceleration theorem (e.g., Andrews et al. 1987)].
Like the Stokes drift for surface gravity waves, a mean flow can be directly produced by nonlinear effects but the resulting changes are weak and reversible: The mean flow disappears after the wave emission has stopped (except when the mean flow itself propagates as Rossby waves; see Bühler & McIntyre 1998).A cumulative, and significant, effect on the mean flow (such as a transport of potential vorticity) occurs only in the presence of irreversible processes, namely, viscosity or diffusion (McIntyre & Norton 1990).In most cases of geophysical interest, these dissipative effects are quite weak on the wave itself, but they act at the end of a turbulent cascade resulting from wave instability and breaking.Therefore the cumulative exchange of momentum is indirectly controlled by nonlinear effects, through the enhanced dissipation.A similar comment applies for the vertical transport of buoyancy (temperature or salinity).
Although it transports momentum from its source to the dissipation region, a single plane wave has no momentum density, in the sense that the mean fluctuating velocity is zero.However, another related conserved quantity can be defined, the wave pseudomomentum p.For a plane wave in a slowly varying medium, this quantity has the expression p = (E/ω)k, like momentum in quantum mechanics or for electromagnetic waves.E is the wave energy density (per unit volume) and the ratio A = E/ω defines the wave action.Mean flow changes due to wave nonlinearity or dissipation amount to a force F acting on the background, equal to minus the rate of change of p.These conservation laws have been generalized to large amplitude waves by Andrews & McIntyre (1978) using a Lagrangian formalism.
The wave nonlinearity is quantified by the Froude number Fr, a ratio between inertial and buoyancy effects where U is a typical oscillating velocity and L ≡ 1/|k| is the wavelength.If one sets for U the extremum horizontal velocity u max induced by the wave, the Froude number becomes equal to the wave "steepness" where c x is the intrinsic (relative to the fluid) phase speed of the wave along the horizontal direction.Waves with s < 1 are statically stable, whereas for s > 1, the isopycnals are overturned and the waves are statically unstable.We expect breaking to occur in the latter case, as for usual surface waves.Breaking processes are reviewed in Section 3. We do not address the properties and implications of mixing by breaking gravity waves, owing to page limitations (see, e.g., Toole 1998, Staquet et al. 2001).
Most internal gravity waves are generated with small steepness, so they are far from the direct breaking condition s > 1.However, this condition can be realized locally in a random superposition of linear waves through constructive interferences.Also, in the atmosphere, the steepening of upward propagating waves simply occurs because of air rarefaction with altitude.Various other processes (discussed in Section 2) can progressively increase the steepness of a monochromatic wave until breaking occurs.A pure plane wave in uniform N is always unstable, whatever its amplitude (in the inviscid diffusionless case): It feeds secondary waves by resonant wave interactions, which therefore steepen.Reflection on a sloping boundary and the interaction with a mean shear flow provide other, purely linear, mechanisms for wave steepening.Breaking results in the excitation of many waves, even in the case of a single initial wave, requiring a statistical description.This is still more true in the oceanic and atmospheric cases, owing to the superposition of many wave sources.Field data and laboratory and numerical experiments are presented in Section 4, and the statistical models of this wave turbulence are discussed in Section 5.

Resonant Interaction Theory
The theory of resonant wave interactions was first introduced in fluid dynamics by Phillips in 1960, in the context of surface waves, and later extended to internal gravity waves.The books by Phillips (1977), LeBlond &Mysak (1978), andCraik (1985) and the review paper by Phillips (1981) provide excellent references to this subject.
Internal waves in a uniform stratification interact through triads of Fourier modes (unlike surface gravity waves).The dynamical variables are either Lagrangian [the displacement ξ = (ξ x , ξ y , ξ z ) of the fluid particles from their unperturbed reference position] or Eulerian (velocity u and density fluctuation ρ ).In the limit of weak amplitude, the two sets of variables are simply related [ ξ = u and ρ = −ρ b (N 2 /g)ξ z ], but their relation becomes very complex at larger amplitudes.
The triad interaction is the first term in an expansion in the wave amplitudes (the Froude number Fr or equivalently the wave steepness).The field variables oscillate with the linear frequency ω on the rapid time scale N −1 , but wave interactions produce modulations of the amplitude and phase on the slow time scale N −1 /Fr.This wave interaction is efficient only when a condition of temporal resonance is satisfied, in addition to the spatial resonance, so that the three wavevectors, k 0,1,2 , and their associated frequencies, ω 0,1,2 , satisfy We follow the convention that both positive and negative frequencies are possible to avoid resonance conditions with negative signs (reversing both ω and k in the complex oscillation α exp(ik.r− iωt) does not change the real wave).When the temporal resonance relation (Equation 8) is not satisfied, the triad interaction produces harmonic perturbations, at frequency ω 0 + ω 1 + ω 2 , without cumulative effects on wave amplitudes.The complex amplitudes α i of each wave of the resonant triad satisfy the dynamical system: where * denotes the complex conjugate and S is an interaction coefficient, depending only on the three wave vectors (see, e.g., McEwan et al. 1972 for explicit calculations).The dot refers to the temporal derivative with the slow time scale.Each wave has an energy density E i = α i α * i , an action A i = E i /ω i , and a pseudomomentum p i = A i k i .The energy E 0 + E 1 + E 2 is conserved by the interacting triad, as well as the pseudo-momentum vector p 0 + p 1 + p 2 (by contrast, the wave action generally is not conserved by triad interaction).The system (Equation 9) can be solved analytically (in terms of elliptic functions), yielding a periodic exchange of energy among the wavevectors.However many triads simultaneously interact, and the corresponding terms must be added, resulting in complex behavior, often described by statistical models, as discussed in Section 5.
General triad interactions, whether resonant or not, may be obtained with the Eulerian variables, as the exact spatial Fourier transform of the quadratic advective terms in the Boussinesq equations (see Ripa 1981).A dynamical equation analogous to Equation 9 is derived for the Fourier amplitudes a i = α i exp(−iωt), without hypothesis on wave amplitude.Because of the different form of the interaction coefficients in this case, the vertical component of the pseudo-momentum is no longer conserved in the nonresonant interaction (but the horizontal components and the total energy are).For weak amplitudes, the only efficient triad interactions are the resonant ones, and the dynamical system (Equation 9) is recovered.The Eulerian variables allow one to study the resonant interactions of the waves with the vortex modes (see Riley & Lelong 2000), which is not possible with the Lagrangian variables (the displacement is never small in a vortex mode as a result of the absence of restoring force).
An important case of resonant triad interaction arises in the problem of stability of a single wave.The interaction with each possible triad involving the primary wave k 0 is then described by the system (Equation 9), with small perturbations, α 1,2 |α 0 |, and the possibility of exponential growth is sought.Hasselmann (1967a) proved by quite general arguments that this is only possible for sum inter- this can be seen directly on the system (Equation 9)].This implies that growth is only possible for perturbing waves with frequency lower (in absolute value) than the primary wave (within the approximation of resonant triad interactions).For internal gravity waves, it appears that the maximum growth rate occurs for perturbations of half the primary wave frequency ω 1 =ω 2 =−ω 0 /2, in the limit of high wavenumbers k 1 ∼ k 2 k 0 .The instability then takes the form of a parametric instability (like that occurring for a pendulum with oscillating suspension).Interestingly it provides a mechanism of direct energy transfer toward small scales, without a turbulent cascade process (Figure 2a).This point was already noted by Pierrehumbert (1986) in the context of the instability of an elliptic vortex, another form of parametric instability.
The stability problem of a single wave is quite relevant for laboratory studies.Davis & Acrivos (1967) first observed this process for gravity waves propagating at a thin density interface, and Martin et al. (1972) did so for waves propagating in a long channel uniformly stably stratified with salted water.The destabilization of standing internal gravity waves in a closed container can be treated like traveling waves.It was studied by Thorpe (1968), McEwan (1971), Orlanski (1972), McEwan & Plumb (1977), and Benielli & Sommeria (1998).In the latter experiments, the instability was observed to grow as a localized wave packet, represented in Figure 2b.
Several numerical studies of internal waves in a vertical plane, either propagating or standing, were carried out during 1970-1980 (a comprehensive review of which is provided by Orlanski & Cerasoli 1981).These studies were mostly performed within an oceanic context and considered a random field of internal waves.The subharmonic parametric instability has been revisited only recently with high resolution numerical simulations, focusing on the spatial structure of the instability and the physical mechanisms of amplification.In Bouruet-Aubertot et al. (1995) and C.R. Koudella & C. Staquet (submitted), a monochromatic low-amplitude high-frequency (ω 0 ) primary wave is considered in a uniformly stratified fluid.It is either standing (in the former paper) or propagating (in the latter).The instability  Staquet (submitted).The vorticity of the growing perturbation is displayed using grey levels (the primary wave field has been filtered out).(b) Case of a standing wave produced in a laboratory experiment by Benielli & Sommeria (1998).The growing wave packet is visualized by the deformation of horizontal dye strips.organizes as localized vorticity bands of alternate sign (see Figure 2a), making a mean angle θ 1 with the horizontal such that cosθ 1 = 1 2 ω 0 /N , consistent with the condition of subharmonic instability.This is in agreement with the laboratory experiments.The selection of the possible perturbation modes is influenced by the quantization set by the domain size.This condition becomes less important for high perturbation wavenumbers, and the selection is then mostly controlled by viscosity.
The theory of resonant interactions has been little applied to the atmosphere, likely because a significant part of the waves reaches large amplitude owing to air rarefaction with altitude: Breaking occurs through instabilities that are more rapid than those stemming from resonant triad interactions, as discussed in the next section.There is no clear evidence of organized triad interactions from oceanic observations, as observed with pure waves in laboratory experiments so that, although the theory was mostly developed in the oceanic context, it was generally used in its statistical version, as discussed in Section 5.

Linear stability theory
The resonant interaction theory is only valid for a small wave amplitude (steepness s).The linear stability of a primary plane wave with arbitrary amplitude in an inviscid fluid has been analyzed by the method of normal modes and Floquet theory, suited to the time periodicity of the primary wave.Drazin (1977) considered mathematically the case of two-dimensional perturbations upon a small amplitude primary wave: Instability then occurs through second-order interactions as well as weaker third-order interactions (of the form k 1 + k 2 + nk 0 = 0 and ω 1 + ω 2 + nω 0 = 0, where n + 1 is the order and the 0 index refers to the primary wave), thus recovering predictions from the resonant interaction theory.An important conclusion of Drazin's work is that an internal gravity wave is always unstable, whatever the value of its local Richardson number.It follows that Miles-Howard's criterion about the stability of parallel steady shear flows does not apply to an internal gravity wave (which locally induces a parallel but oscillating shear flow).The systematic investigation of the two-dimensional stability problem was initiated by Mied (1976) and complemented by Klostermeyer (1982Klostermeyer ( , 1991)).These authors numerically solved the eigenvalue Floquet problem for a statically stable primary wave (s ≤ 0.7).The main conclusion is that the subharmonic parametric instability predicted by resonant interaction theory for small-scale perturbations remains valid when s is increased.Klostermeyer (1991) also investigated the stability of a primary wave against three-dimensional perturbations, in a plane normal to the primary wave plane.He considered only one direction of the wave vector (θ = −72 • ), but his conclusions were confirmed by a more systematic exploration of the primary wave and perturbation parameters by Lombard & Riley (1996a) and Sonmor & Klaassen (1997).Klostermeyer (1991) showed that the primary wave is most unstable to twodimensional, small-scale perturbations when its steepness s < 1.When s > 1, he found that three-dimensional perturbations grow fastest, the instability being either of the harmonic or the subharmonic parametric type; in addition, a purely growing mode of zero frequency is obtained, which is a vortex mode.Lombard & Riley (1996a) actually showed that a two-dimensional instability with comparable growth rate may also be present, depending on the propagation angle of the primary wave.The underlying structure of the three-dimensional instability appears when s → 0: It consists of resonant interactions of order higher than 3. Lombard & Riley (1996a) showed that the wavenumber range of the fastest growing, threedimensional modes is broad-banded and has no upper bound (in the inviscid limit), which reflects the underlying resonant interaction process.Two-dimensional instability modes have been detected in the atmosphere by sensitive Doppler radars (Klostermeyer 1990), but no three-dimensional modes have been identified; however, the structures of streaks in noctilucent clouds displayed in Figure 1 might be explained in terms of the three-dimensional instability of an internal gravity wave (Fritts et al. 1993).Lombard & Riley (1996a) and Sonmor & Klaassen (1997) investigated the energy sources of the perturbation: These are the shear and the buoyancy gradients of the primary wave, which feed the kinetic energy and the potential energy of the perturbation, respectively.Sonmor & Klaassen (1997) thus showed that the two-dimensional small-scale parametric instability is primarily buoyancy driven.This result implies that potential energy is transferred toward the small-scale perturbation at a higher rate than kinetic energy, a behavior also observed for stratified turbulence by Holloway & Ramsden (1988).This point is further discussed in Section 5.4.In general, both forcing processes occur whatever s is so that, according to Lombard & Riley (1996a), no classification into a shear-driven or a buoyancy-driven instability seems possible.Sonmor & Klaassen (1997) made a step forward in proposing a classification of the fastest growing linear instabilities as a function of the primary wave steepness and propagating angle (see their Figure 22).

Focusing by Wall Reflection
Internal waves obey simple, but unusual, reflection laws at a rigid sloping boundary.In optics or acoustics, the incident and reflected wave rays make the same angle with respect to the normal direction to the reflecting surface, whereas for internal waves they make the same angle with respect to the gravity direction.This is a direct consequence of the peculiar dispersion relation (Equation 3): Because the frequency is preserved by reflection on a fixed wall, the angle θ must indeed be preserved as well.The ray tube sketched in Figure 3 is then focused into a narrower ray tube, leading to a concentration of the energy flux (of course the reverse process of defocusing reflection is equally possible, but it does not lead to the interesting steepening effects).Near the critical angle, for which the reflected ray is along the wall, the linear theory predicts a reflected wave of infinite amplitude, infinitesimal wavelength, and zero group velocity, leading to the trapping of oncoming wave energy in the boundary region (as first discussed by Phillips 1966).
Dauxois & Young (1999) nicely specified, using a matched asymptotic analysis, how the linear theory breaks down near the critical angle.The wall region is analyzed as a boundary layer.As the group velocity of the reflected wave tends to zero, the time for wave interation with the wall diverges accordingly, so that a time dependent model is considered: The incident wave amplitude is turned on at some time origin.A first nonlinear effect is the possible emission of a wave Figure 3 Reflection of an internal gravity wave on a sloping wall near the critical angle.harmonic.Concerning the wall layer itself, an array of counter-rotating vortices is formed, whose transverse width (normal to the wall) decreases with time t in t −2 at the critical angle, while transverse oscillations build up.Density overturning and wavebreaking are, however, expected to occur before this, owing to the wave steepness amplification in t 2 .(The implication of breaking on the fluid is discussed in Section 3.3.)In contrast to these inviscid results, viscosity and diffusion stabilize a steady solution even at the critical angle (Kistovich & Chashechkin 1995).
Another possibility of strong linear wave amplification occurs by successive reflections in a closed domain, as first proposed by Maas & Lam (1995).In a domain with a sloping boundary, rays exponentially approach a closed circuitan attractor-after several reflections.The corresponding standing waves [timeperiodic solutions of the wave equation (Equation 2)] have a fractal structure along the attractor, as shown in Figure 4. Smooth eigenmodes are recovered by introducing a weak viscosity.They tend to an inviscid fractal solution in the limit of small viscosity, as calculated by Rieutord et al. (2001).These authors consider a generalization to inertia-gravity waves in a spherical envelope, describing oscillations in a rapidly rotating star.Maas et al. (1997) experimentally observed the existence of such a wave attractor.A strong wave steepness is obtained along the attractor, producing some turbulence.This could have applications for mixing in closed oceanic basins or lakes, but it requires that waves keep their coherence for a few turns around the attractor.This will depend on the effectiveness of scattering by nonlinear wave interactions (see Section 5).

Figure 4
Standing wave inviscid solution in a basin with sloping side wall (from Maas et al. 1997).The streamfunction, represented in grey levels, has a fractal structure.

Interaction with a Mean Shear Flow
An internal gravity wave may also amplify by interaction with the mean shear flow in which it propagates.In the atmosphere, for instance, a horizontal wind with a vertical shear acts as a barrier to vertically propagating waves if the horizontal phase velocity matches the horizontal velocity of the wind at some altitude.This altitude is a critical level.If propagating upward, the wave is trapped below that level and may break there, thereby irreversibly transferring energy toward smaller scales.The wind thus plays a filtering role in upward propagation.Henyey et al. (1986) showed, in the ocean, that small-scale waves encounter a critical level within a few inertial periods when propagating through a realistic shear, mostly as a result of the internal wave field.Hence, interaction with a mean shear is one of the mechanisms through which internal waves provide a link from large to small scales.The reverse link should not be forgotten: Waves transport momentum over long ranges, as illustrated by McIntyre (2000).
The mean flow with which the waves interact may be external to the waves, as is a wind or an oceanic current, or induced by them, because of nonlinear or dissipative effects.In the former case, theoretical understanding has first been gained using WKB (or eikonal) approximation (e.g., Olbers 1981, LeBlond & Mysak 1978).The medium in which the wave propagates, characterized by the Brunt-Väisälä frequency N (x, t) and a local background mean flow ū(x, t), is assumed to be weakly inhomogeneous and weakly unsteady, in the sense that it varies on a much larger scale than the wavelength and the wave period.The dispersion relation is thus assumed to hold locally, whereas the waves feel the mean gradients (of N and ū) through propagation and refraction only.The wave dynamics is then described in the absolute reference frame (relative to which ū is measured) along ray paths (ẋ = c g + ū) by the classical ray equations k ≡ ∂k/∂t + (c g + ū)∇k = −∇ ω.The absolute frequency ω is defined by where the relative frequency ω is given by the dispersion relation (Equation 3).If the fluid medium is at rest and nondissipative, the intrinsic wave energy density E is conserved along ray paths.By contrast, in the presence of a background mean shear, the wave action density A = E/ω is conserved instead of the energy (Bretherton 1966, Bretherton & Garrett 1969): When dissipative or of finite amplitude, the waves may induce mean flow changes, and Grimshaw (1974) introduced a transport equation for the slowly varying mean flow coupled to that for the wave action density (also Andrews & McIntyre 1978).VERTICALLY VARYING MEAN FLOW The case of a steady horizontal mean flow with a vertical shear ū(z)i x has been intensively studied since the theoretical work by Bretherton (1966), who used a linear inviscid WKB analysis.If the wave propagates in such a way that k x ū increases, then the intrinsic frequency ω decreases from Equation 10(because ω remains constant in a steady medium) and the vertical wavenumber k z increases from Equation 3.This is a common situation in the atmosphere, as the waves generally propagate upward and the wind increases with altitude.If a critical altitude z c exists at which ω vanishes (in the absence of Coriolis effect), Bretherton (1966) showed that, for a wave packet, , where c gz is the vertical component of the group velocity (Figure 5a).Also, the distance |z − z c | varies as t −1 so that an infinite time is needed for this level to be reached.This unphysical behavior results from the linear inviscid WKB solution being singular at the critical level.
Relaxing the WKB approximation, Booker & Bretherton (1967) showed that if the local Richardson number associated with the mean flow Ri is everywhere greater than 1/4, then the wave energy density is attenuated by a factor exp [−2π(Ri c − 1/4) 1/2 ] as the wave passes through the critical level.Hence, no reflected and transmitted components are produced for Ri 1/4, implying that the wave is fully absorbed and irreversibly transfers all its energy to the mean flow.In spite of this absorption, the steepness diverges (because c x decreases faster than u ) and in practice, breaking occurs before the critical level is reached (McIntyre 2000).
The mathematical singularity of the linear inviscid solutions can be removed by taking into account viscous or nonlinear effects near the critical level, in a region called the critical layer.Hazel (1967) studied the viscous critical layer.Brown & Stewartson (1982) investigated the effect of nonlinearity (see also Maslowe 1986).They found no transmitted component, but a reflected wave is produced by finite amplitude effects, implying that the wave absorption by the mean flow is no longer complete.This result is confirmed by the three-dimensional numerical experiments of Winters & D'Asaro (1994), as specified in Section 3.3.
Mean flow changes induced by finite amplitude effects may in turn modify the wave properties (such as its intrinsic frequency), a process Fritts & Dunkerton (1984) called self-acceleration.By this process the wave can propagate well beyond its initial critical level [an analogous behavior is obtained for a wave packet propagating toward a reflecting level, at which ω = N (Sutherland 2000)].Winters & Riley (1992) investigated the stability of an internal gravity wave near a critical level.A simple model for the primary wave was inferred from the classical predictions of the WKB theory (without the Coriolis effect).Two limiting cases of amplified perturbations were found, with comparable growth rate: a two-dimensional Kelvin-Helmholtz instability in the plane of the wave and a buoyancy-induced instability, commonly referred to as a convective instability, in the transverse plane.The three-dimensional nature of the latter is consistent with Deardorff's (1965) conclusion that convection is damped in the plane of a weak shear flow.Any unstable perturbation in general should be a combination of these limiting cases, which is consistent with the conclusions drawn by Klostermeyer (1991) and Lombard & Riley (1996a) for a large amplitude primary wave in the absence of mean shear.
Finally, we note that a word of caution is exerted in some papers about the importance of wave breaking near critical levels.As stressed by Badulin & Shrira (1993), this may be a consequence of the problem idealization, which involves a unidirectional parallel shear flow.Broad (1999) and Shutts (1998) found that, if the mean wind changes direction with height, the waves need not be trapped below a critical level and, indeed, may not break at all.Also, as the intrinsic frequency decreases, the wave dynamics becomes influenced by Earth rotation and by low frequency variations in the mean current.
HORIZONTALLY VARYING MEAN FLOW The interaction of internal gravity waves with a horizontal shear flow has received much attention from Russian scientists, within an oceanic context, using linear inviscid WKB theory.Interest in this topic has grown recently, motivated by the permeability of dynamical barriers (such as the edge of the polar vortex or oceanic currents).The interaction of an internal gravity wave with a horizontal shear flow may lead to a "trapping level" analogous to the critical level of a vertical shear.Ivanov & Morozov (1974) and Basovich & Tsimring (1984) first considered a steady barotropic horizontal parallel flow, e.g., ū(y)i x , in a uniformly stratified fluid (the problem is thus inherently three dimensional).The approach of the trapping level is characterized by a decay to zero of both the intrinsic group velocity and the wavelength along the direction of inhomogeneity y.Because the vertical wavenumber component is conserved in this process, the wavevector tends to the horizontal y direction (Figure 5b).Therefore the intrinsic frequency tends to N at the trapping level, instead of zero for the critical level in a vertical shear.
Furthermore the momentum of the mean flow is irreversibly transferred to the wave, contributing to wave amplification.This problem has been addressed by Staquet & Huerre (2001) using direct numerical simulations of an inertia-gravity wave: Breaking may indeed be observed in the neighborhood of a trapping plane.Olbers (1981) and Badulin et al. (1985) considered the case of a baroclinic horizontal parallel shear flow in a stably stratified fluid.The shear flow is associated with a horizontal density gradient by the thermal wind balance (in addition to the vertical stratification).Olbers (1981) showed that the trapping mechanism is not stable when the Brunt-Väisälä frequency is constant: Even a weak vertical shear triggers reflection when ω approaches N. The trapping process becomes stable when a vertical inhomogeneity of N exists.In this case, the wave amplitude increases algebraically, and its energy appears to focus at the depth of maximum value of N. Evidence for such a trapping phenomenon in the Lomonosov current is provided in Badulin et al. (1990).This work was further extended by Badulin & Shrira (1993) to the propagation of guided waves in an arbitrary inhomogeneous background current field, again within the frame of the WKB approximation.Preliminary conclusions confirm the tendency for irreversible wave transformation to small scales.

INTERACTION WITH THE WAVE-INDUCED MEAN FLOW
The self-acceleration process put forward by Fritts & Dunkerton (1984) was further studied by Sutherland (2001), who considered a wave packet of finite amplitude propagating in a twodimensional medium at rest.On the basis of a study performed by Grimshaw (1977), Sutherland proposed that the wave packet may resonantly interact with its induced mean flow, considered as a wave with infinite horizontal wavelength, and become unstable.The point is that the critical wave amplitude for the occurrence of this instability appears to be smaller than for overturning and, a fortiori, for convective instability: This resonant interaction mechanism may therefore be relevant, in agreement with two-dimensional numerical simulations performed in parallel.

General Remarks
Wave breaking is defined as the production of turbulence and irreversible energy dissipation.As intuitively expected, it corresponds to the occurrence of isopycnal overturning, defined in the introduction by the condition s ≡ u max /c x > 1 for the wave steepness s.However overturning does not always imply breaking, but it is a necessary condition for breaking (e.g., Thorpe 1994b).An equivalent definition of s (for a monochromatic wave) is s = A ξ k z , where A ξ denotes the vertical displacement amplitude (e.g., Thorpe 1999a): This shows that a strong steepness can be associated with a very weak displacement amplitude of the wave.An example is provided by secondary waves resonantly interacting with a (weak steepness) primary wave.The velocity amplitude of the secondary waves cannot exceed that of the initial primary wave because they draw energy from it, but their steepness can.The secondary wave steepness thus grows through the interaction, via a parametric instability; when their steepness reaches a critical value of order one, the secondary waves becomes unstable themselves and break.It should be noted that the primary wave steepness does not grow at all during this process and remains very small.The breaking of the secondary waves eventually dissipates the energy of this primary wave so that the latter wave is sometimes said to break.It would be more appropriate to say that it loses its coherence; what break in this process are the secondary waves fed by the primary wave through resonant interactions.

Breaking Processes
We restrict ourselves in this section to a deterministic description of the breaking process so that the idealized case of a monochromatic wave is considered.The breaking process thus results from the growth of an instability upon a basic state made of a wave with a steepness larger than one.[The influence of the Coriolis force profoundly modifies this view (see Lelong & Dunkerton 1998, Thorpe 1999a).]This instability is three dimensional and is mainly driven by buoyancy forces, as discussed in Section 2.2.Lombard & Riley (1996b) numerically found that a monochromatic wave propagating in a quiescent fluid breaks through a combination of shear and convective instability.For a higher numerical resolution (256 3 ), Koudella (1999) concluded that the instability leading to breaking is dominated by unstable buoyancy effects, even for small initial steepness of the primary wave (s 0.3).Afanasyev & Peltier (2001), when studying the breaking processes of a growing amplitude lee wave, found as well that breaking occurs through a convective instability once the wave front has overturned.The formation of small scale convective structures in the numerical experiments of Koudella (1999) is illustrated in Figure 6a.In a stably stratified fluid, those mushroom-like structures indicate that a dominant source of potential energy is available for the perturbation.Indeed, they are also observed in the stratified shear layer, inside the Kelvin-Helmholtz vortices, as the result of statically unstable regions that form as the vortices roll up (Figure 6c).Note that alternate mechanisms may lead to such structures, as they are a common feature of three-dimensional unstable shear flows.For instance, they are found in the unstratified shear layer, in between Kelvin-Helmholtz vortices, as a result of the intense strain (Schowalter et al. 1994 and references therein).
These conclusions remain unchanged in the presence of a background shear: The pioneering direct numerical simulations of Winters & D'Asaro (1994), performed at the moderate resolution of 32 × 32 × 200, conclude again that both shear and unstable buoyancy drive the instability, leading to breaking with no dominant energy source.At a higher resolution, an internal wave train appears to break through a convective instability (which next triggers a dynamical shear instability).The latter results have been obtained for internal gravity waves forced by a topography (Dörnbrack 1998) and by a localized body force (Andreassen et al. 1998), using large-eddy simulations.These above three articles concluded that the initial processes leading to the steepening of the wave are two dimensional [in agreement with the finding by Klostermeyer (1999)] but that the breaking process itself is genuinally three dimensional.Fritts et al. (1994) suggested that two-dimensional computations with an appropriate parameterization of the breaking process may actually provide a satisfactory description of wave-shear interaction.
Laboratory experiments display a somewhat more complex view when the breaking process is considered.The laboratory experiments performed by McEwan (1983), Thorpe (1994c), and Benielli & Sommeria (1998) (Figure 6b) provide evidence of convective instability and mixing.McEwan (1983) showed that breaking occurs through very localized and intermittent small scale unstable events, a crucial property when mixing or statistical modeling is to be determined.The same experiment was reproduced by Taylor (1992), but as reported by Thorpe (1994b), a variety of structures and mechanisms leading to small-scale energization was found instead.Koop & McGee (1986) studied breaking in critical layer experiments.They inferred a simple model from the equations for the wave action density and induced mean flow derived by Grimshaw (1974), which they compared with their carefully instrumented experiments.The model provides accurate predictions of convective instability near the critical level and of Kelvin-Helmholtz instability when no critical level exists.
The interaction of a wave with a sloping boundary, whether inferred from numerical or laboratory experiments, does not allow for a simple description of the instability mechanism, because of the solid boundary.Dauxois & Young (1999) still argue that the instability should be convective (because overturning occurs before the Miles-Howard stability condition is violated), in agreement with Taylor's (1993) experimental measurements.Slinn & Riley (2001) actually showed that the details of breaking depend upon the slope of the bottom.

Induced Mean Flow Changes
The major motivation in wave breaking studies, beyond the knowledge of the breaking process itself, lies in the resulting changes to the background.Koudella (1999) computed the mean flow induced by the breaking of a progressive statically stable primary wave (s = 0.54 and s = 0.72).He found that, at the end of the computation, the mean flow kinetic energy is as high as the total kinetic energy in the wave part of the flow [using decomposition (Equation 4)].
When a wave packet interacts with a critical level, Winters & D'Asaro (1994) showed that the first effect (in time) of the wave upon the fluid is the transfer of momentum to the mean flow (Ri = 25 in these simulations).This effect also appears to be the largest energy sink; it increases as the initial wave steepness decreases, as expected from linear theory.The second most important energetic process is wave reflection, which increases with s.As for the transmitted component, it remains smaller than a few percents of the initial wave energy, whatever s is.As a consequence, little energy is left for energy dissipation, for mixing, and for the production of potential vorticity.
The interaction with a sloping boundary, discussed in Section 2.3, leads to wave breaking, forming a bore of mixed fluid (if the bottom slope is high enough, 30 • , according to Slinn & Riley 2001), which propagates along the wall; this eventually initiates horizontal intrusions in the stratified outer region, as observed in several laboratory experiments (e.g., Ivey et al. 2000).The numerical simulations by Slinn & Riley (1998, 2001) also carefully document the energy budget at critical angle.According to the (two-dimensional) numerical simulations by Javam et al. (1999), the wave reflectivity decreases as the angle of incidence approaches the critical value, owing to enhanced energy dissipation.
With an oblique incidence (tilted with respect to the normal vertical plane), wave absorption must generate alongslope currents owing to the release of the momentum transported by the wave, as observed with surface waves.Thorpe (1999b) provided some theoretical estimates of this effect, and Dunkerton et al. (1998) have performed laboratory experiments.They suggested the interesting possibility that even a horizontally isotropic incident wave field could spontaneously produce a transverse flow as an instability process.A numerical study carried out by Zikanov & Slinn (2001) shows that the wave does not penetrate down to the boundary, owing to its interaction with the induced alongslope current at a critical level.As a result, small scale mixing processes are produced away from the boundary.
In the oceanic case, with a random superposition of incident waves, we expect that the waves approaching the critical frequency will be selectively amplified near a sloping bottom.Measurements by Eriksen (1998) support this idea, and various observations suggest that mixing and sediment suspension are indeed enhanced by this effect on the continental shelf; see the introductions of Slinn & Riley (1998) and Dauxois & Young (1999) for further references.

Ocean
We refer to Munk (1981) for an excellent introduction to internal waves in the ocean and to Müller et al. (1986) and Olbers (1983) for more advanced reviews.The ocean is generally stably stratified in the so-called thermocline, below the upper mixed layer that is 100 m deep on average.Although coherent oscillations are excited by tides in coastal regions (baroclinic tide), internal wave activity mostly occurs as random fluctuations.Their energy and statistical properties show remarkable similarities at any time in most parts of the ocean.Various spectral measurements have been summarized in the model of Garrett & Munk (1979) (GM), which has been empirically improved in successive versions (see Munk 1981 for the latest).Strong deviations from this form have been observed only in special places under polar caps, near the equator, near strong currents and sloping bottom.In the latter case, the GM spectrum was found to relax back to the universal form only a few hundred meters off the bottom (see Eriksen 1998).The GM spectrum summarizes measured temporal and frequency spectra into a single energy spectrum E zω (k z , ω) depending both on vertical wavenumber k z and frequency ω and assumed to be shaped by second order nonlinear interactions.A superposition of linear waves is assumed so that any spectral quantity can be deduced: For instance, the vertical energy spectrum E(k z ) is given by the integration of E zω (k z , ω) over ω.The GM spectrum scales as k −2 z for high vertical wave numbers k z and in ω −2 for high frequencies ω.It is dominated by the lowest vertical wavenumbers and the nearinertial frequencies (ω close to the Coriolis frequency f ).The corresponding wave vector is close to vertical, with velocity dominantly horizontal.The total kinetic energy is three times the total potential energy, expressing the dominance of nearinertial waves (equipartition occurs for pure gravity waves).
The Froude number at the dominant vertical wavenumber k z ∼ 10 −2 m −1 is u k z /N ∼ 0.1 (with u ∼ 7 cm/s) justifying a linear wave model, but it increases for higher wave numbers.The typical velocity at scale k −1 z can be estimated from general dimensional arguments as u(k z ) ∼ √ k z E(k z ), and the corresponding Froude number It increases in k 1/2 z in the GM spectrum, and Fr(k z ) ∼ 1 for k z ∼ 1 m −1 .Beyond this value, a k −3 z energy spectrum, or equivalently k −1 z spectrum of shear and temperature gradients, is observed over one decade (see Figure 7a).Remarkably a universal spectral level is observed, for both kinetic (E c ) and potential (E p ) energy (expressed per unit mass): Parameters α and β remain close to 0.2 (Gregg 1987).This spectral range corresponds to a constant Froude number among the scales, by Equation 12.It can be viewed as the result of a "saturation" process: Beyond this value, energy is quickly dissipated by breaking events.Unlike the GM spectrum, this spectral range corresponds to genuinely nonlinear, or "compliant" waves, meaning that propagation is dominated by advection owing to larger waves.Finally, at centimeter scales, sporadic patches of turbulence are observed, providing the last step of a global cascade degrading the energy-containing waves.The time of energy decay has been estimated at approximately 100 days, partly explaining the remarkable wave persistence in the presence of variable sources.

Atmosphere
Internal waves in the atmosphere are studied by rising balloons, providing highresolution vertical profiles of temperature and velocity (e.g., Dalaudier et al. 1994).Remote sensing, using radars or lidars, yields a more global coverage but with indirect data processing.Although internal waves can be observed as coherent oscillations above mountains, they appear in most cases as random motion, for which a statistical description is required.In interpreting the data, a fundamental question is whether the observed randomness reflects the source variability or is the consequence of the complex wave dynamics itself, involving scattering, refraction, and breaking.Probably all effects are influent.A given location in the stratosphere indeed receives waves emitted over thousands of square kilometers near the ground.Thus the influence of many independent sources is superposed, with the possibility of multiple scattering along the propagation path.Although the waves dominantly propagate upward, a significant fraction of their energy is backscattered downward (Barat & Cot 1992).
Many measurements yield the spectral law (Equation 13), in the range of vertical scales from a few kilometers (increasing with altitude) to about 100 m, see Figure 7b.The prefactors (α and β) seem however less constant than in the ocean.They decrease with altitude and change with time (within an order of magnitude), suggesting less universal "saturation" mechanisms.Furthermore, tendencies to smaller spectral slopes (−2.5) are observed in the upper stratosphere.In the horizontal direction, k −3 x spectra have been observed for scales smaller than a few kilometers, but with smaller slopes at lower k x (Bacmeister et al. 1996).Horizontal spectra in k −2 x and time spectra in ω −2 have been measured simultaneously with k −3 z vertical spectra by Hostetler & Gardner (1994).
As in the ocean the wave energy is dissipated during breaking events.In these local overturning regions, typically less than 100 m high, Kolmogoroff spectra in 2/3 k −5/3 are observed.Their measurement allows the evaluation of the energy dissipation rate (see, e.g., Alisse & Sidi 2000).Similar Kolmogoroff spectra, depending on the dissipation rate of temperature variance, are obtained for the temperature fluctuations.These dissipation rates are of high interest, as they are related to the vertical transport properties of the wave field.

Laboratory and Numerical Experiments
Whereas laboratory experiments have addressed various elementary wave processes, such as resonant interaction and instabilities, few studies have been concerned with the statistical aspects.Stratified grid turbulence has been extensively studied, but the Froude number is quite high in this case.The turbulence obtained at low Froude numbers, from wave breaking, has quite distinct properties, but it is not clear to what extent universal regimes, which disregard the details of the forcing conditions, can be obtained.Time spectra from a fixed local density probe have been measured by Hannoun & List (1988) and Benielli & Sommeria (1996).The spectra, now depending upon the frequency ω, display a ω −3 law, and when properly scaled, superpose whatever the flow conditions (Figure 7c).
It should be noted that this ω −3 frequency range is beyond the buoyancy frequency N, so that it cannot be related to weakly linear waves.Benielli & Sommeria (1996) interpreted the observed ω −3 spectra as spatial spectra advected by the waves with vertical velocity U, using a Taylor's hypothesis k z = ω/U .For a variety of experimental conditions, they deduce the same spectral law (Equation 13) for E p with α = 0.2.We expect from these results that the minimum wavenumber of this compliant wave spectrum is N /U (reached for ω = N ).Direct spatial measurements would be useful to check these ideas.
Early numerical simulations of breaking internal gravity waves did not reach clear conclusions about the spectral level (Orlanski & Cerasoli 1981), owing to insufficient resolution.Using 512 2 direct numerical simulations of a standing internal gravity wave, of moderate amplitude, Bouruet-Aubertot et al. (1996) found an inertial range displaying clear k −3 z spectra after wave breaking.The spectra lock onto the law (Equation 13) with constant spectral levels α and β, persisting for many primary wave periods, in spite of a global decay of the energy.This provides support for an interpretation in terms of wave saturation.Carnevale et al. (2001) have recently obtained similar results in three-dimensional computations using subgrid scale modeling.The results display a tendency to a k −5/3 Kolmogoroff cascade beyond the k −3 range (Figure 7d), in remarkable agreement with atmospheric data.The analysis of the buoyancy flux u z θ in the simulations by Ramsden & Holloway (1992) indicates that the potential energy feeds the kinetic energy at small scales.Bouruet-Aubertot et al. (1996) confirmed this result and found a power law in k 0.75 for the correlation of u z and θ , in agreement with theoretical estimates by Holloway (1988).
In summary, all these experiments and simulations indicate that the k −3 z range with constant saturated levels α ∼ β can result from a single wave breaking.A more global spectral analysis, involving the horizontal direction and time, is needed, allowing in particular a direct evaluation of the departure from the linear dispersion relation.More advanced statistical analysis would also be needed to characterize the strong intermittency of this wave turbulence.The recent analysis of atmospheric data by Alisse & Sidi (2000) leads the way in this field.
The GM spectrum corresponds to a regime with weaker wave interactions.In numerical simulations by Winters & D'Asaro (1997), a GM spectrum was found to be steadily maintained by a random large scale forcing.The corresponding rate of energy dissipation is in good agreement with oceanic data, providing a further test of the theoretical models (Figure 8).This opens wide perspectives for the direct numerical modeling of oceanic waves.

STATISTICAL MODELS OF WAVE TURBULENCE
We have shown in the previous section that universal wave spectra are observed in most of the atmosphere and the ocean, in the presence of a variety of wave sources.It is therefore reasonable to seek an interpretation in terms of a local statistical steady state resulting from complex wave interactions.Then the details of the source do not matter.

Linear Models and Saturated Energy Spectra
The description of linear wave propagation in various stratification and shear flow conditions is a first step for parameterization in circulation models.In a linear analysis, the wave spectrum corresponds to the forcing spectrum, multiplied by some linear response function characterizing the wave dynamics.
In the mid ocean, far from coastal effects and intense currents, the main source of internal waves seems to be the wind, acting through different mechanisms (see Olbers 1983 for a review).Rubenstein (1994) has considered the linear response of internal waves to random propagating wind stress perturbations, represented by an empirical fit to observations.He has obtained some similarities with the GM spectrum and was able to explain the observed dependency of the dominant vertical modes with the thickness of the upper mixed layer, at the bottom of which the waves are emitted.This dependency results in a buffer effect, partly explaining the observed constancy of internal wave energy.The results of this model, however, depend on little-known parameters describing the wind stress statistics and also on an empirical damping coefficient.This damping is in fact an implicit representation of the cascade process caused by nonlinear wave interactions.
In the atmosphere, the simplest assumption is that the waves propagate upward from their source, according to an eikonal approximation.Wave breaking is then represented by an ad hoc dissipation model, maintaining the Froude number smaller than 1.Such saturated wave models have been adapted to various atmospheric conditions (see Fritts 1984, Warner & McIntyre 1996).

Weakly Interacting Waves
Kinetic theories of waves aim at finding closed evolution equations for the energy spectra.The case of weakly interacting internal waves was first derived by Hasselmann (1966Hasselmann ( , 1967b) ) and is reviewed in detail by Müller et al. (1986).It is obtained as a first-order perturbation to freely propagating waves and relies on general methods of field theoretical physics, first developed to describe thermal fluctuations in crystals.The general approach for waves in continuous media is discussed in the book by Zakharov et al. (1992), but the standard Hamiltonian description used in this book does not apply directly to internal waves (see Caillol & Zeitlin 2000).
The starting point is a set of free waves with random phases, which is justified for weak wave interactions.Then wave packets from far distant origins, hence with uncorrelated phase, interact at a given location (see Hasselmann 1967b).Dispersive effects also favor random phase, as the Fourier components of a single wave packet become increasingly spread in phase with time.The existence of coherent structures, such as solitons or shocks, is excluded by the theory.
The probability distribution of the wave amplitudes a i is Gaussian in the limit of weak interactions.This can be roughly understood as a consequence of the central limit theorem, considering each spatial Fourier amplitude a i (the index i represents a wave vector) as the sum of contributions from many distant independent oscillations.The correlation relations a i a * j = a i a * i δ i j and the closure (exact for a Gaussian distribution) directly result from the random-phase assumption (here the brackets denote an ensemble average over many realizations).Furthermore, the triple correlations (and all the odd moments) vanish.However, this property cannot be exact, as the energy transfers among the waves are associated with these triple correlations.
Similarly, the dynamical equation for the triple correlations depends on fourthorder correlations.The closure problem is solved thanks to Equation 14, as if the field were Gaussian.As a result, the triple correlations a i a j a k oscillate in time at frequency ω i + ω j + ω k and amplitude (ω i + ω j + ω k ) −1 , and a net effect on the energy equation only occurs for the resonant triad interactions, for which ω i + ω j + ω k = 0.The resulting kinetic equation for the energies a i a * i has therefore been called the resonant interaction approximation (RIA).It formally resembles Boltzmann's equations for a perfect gas.The wave energy density a i a * i , depending on the wave pseudo-momentum p i , replaces the distribution function in particle momentum, and the effect of scattering is analogous to binary collisions, conserving both momentum and energy.
The earlier approaches, by Hasselmann (1966) and followers, use Lagrangian variables.More recently Caillol & Zeitlin (2000) have used Eulerian variables, allowing an introdution of interactions with the vortex modes, but they found no interaction at lowest order.This confirms results obtained directly on triad interactions, without statistical averaging (Riley & Lelong 2000).The wave-vortex interactions can only "catalyze" wave-wave interactions, which should favor randomness of the wave field.In the absence of vortex modes, the Lagrangian and Eulerian variables are linearly related in the limit of small amplitudes, as mentioned in Section 2.1, so they should give the same kinetic equations.
The RIA theory has provided the first quantitative understanding of the statistics of oceanic internal waves (Olbers 1976).The GM spectrum is viewed as a steady energy cascade from large vertical scales to small vertical scales.An external forcing at large vertical scale and a small scale dissipation are therefore assumed.A mechanism of removal of near inertial wave energy also needs to be invoked, to balance an excessive build up of energy predicted by RIA.The existence of a k −2 z vertical wavenumber spectral range with ω −2 frequency spectra are well explained by the theory, using analytical arguments (McComas & Müller 1981).Caillol & Zeitlin (2000) also found k −2 z vertical spectra, but associated with a ω −3/2 frequency spectrum.They have obtained this result as self-similar solutions (without Coriolis effects), in the spirit of the Kolmogoroff theory, whereas McComas & Müller (1981) consider a finite frequency domain, bounded by the Coriolis frequency f, which may explain the difference.Numerical simulations of the RIA equations by McComas & Müller (1981) indicate a time evolution converging to the k −2 z ω −2 spectra, in agreement with the GM spectrum.This would exclude the k −2 z ω −3/2 spectra, but improved computations would be useful, with a careful comparison of the different RIA approaches.
RIA calculations furthermore indicate that the GM spectrum quickly relaxes to its steady state after various distortions, explaining the universality of this spectrum.Idealized distortions were actually introduced, which prevents quantitative comparisons with observations, for instance in the case of a reflection on a sloping bottom.Müller (1976) related the random wave momentum transport ("eddy" viscosity) to relaxation in the presence of a shear but found a strong overestimation in the vertical direction (see Müller et al. 1986).Refined models would be useful.
A major outcome of the RIA is the finding by McComas & Bretherton (1977) that among the many resonant wave interactions, three classes dominate the dynamics in the GM spectrum.
Elastic scattering denotes the backscatter of a downward-propagating highfrequency wave into an upward-propagating wave (or the reverse) by interaction with a low-frequency near-inertial wave.There is little energy exchange in this process, and the vertical wavenumber is unchanged in modulus but reversed in direction.By this mechanism the upward or downward origin of the waves is quickly forgotten, explaining the symmetry of the GM spectrum with respect to the propagation direction.This phenomenon is probably active as well in the atmosphere, but the down-going waves decay and the up-going waves are systematically amplified by air rarefaction, so the latter are always dominant.
Parametric subharmonic instability is the decay of a large-scale wave into two smaller-scale waves of approximately half the frequency.This process was described as a deterministic mechanism in Section 2.1.It accumulates energy at the lowest available frequency f, contributing to the near-inertial spectral peak.
Induced diffusion denotes the scattering of a high wavenumber high-frequency wave into another nearby wave vector by interaction with a large-scale low frequency wave.This process, obtained as a particular limit of the energy transfer integral, can be viewed more physically as a wave packet propagating through a random shear, with scales much larger than the packet size.The wave packet then experiences random perturbations of its wave vector, prescribing a random walk.An equation of diffusion for the high-frequency action spectrum A(k) = E(k)/ω(k) results.Because the shear depends mostly on the vertical coordinate z, the diffusion is limited to the z direction, where the diffusion coefficient D depends on the shear spectrum at the vertical wavenumber ( f /ω)k z .The action A(k) diffuses toward larger and larger k z along lines of constant horizontal wavenumber k H , therefore transferring action to lower frequencies ω ∼ N k H /k z .The total wave action A(k)dk z is conserved in this diffusion process, in agreement with an eikonal transport (see Section 2.4), so that the energy (A(k)/ω)dk z is not conserved: It interacts with the large-scale low-frequency background.
The dynamics of this low-frequency wave background is controlled by the parametric subharmonic instability.Adding the two effects at any given vertical wave vector k z , McComas & Müller (1981) found a constant-flux energy cascade in vertical wavenumber.This downscale cascade corresponds to a k −2 z vertical energy spectrum, as observed in the GM spectrum, from a forcing scale k −1 0 10 −2 m −1 .McComas & Müller (1981) computed the total dissipation rate for the GM spectrum.Assuming that 1/5 of this energy is used for vertical mixing, they deduced a wave diffusivity of 10 −5 m 2 s −1 , in agreement with various estimates.The prediction, however, scales as the square of the little-known forcing scale k −1 0 .The limitations of the RIA have been discussed in detail by Holloway (1980) and Müller et al. (1986).Wave interactions necessarily lead to highly nonlinear (steep) waves for which the approximation fails, but the range of validity is difficult to evaluate.Carnevale & Frederiksen (1983) pointed out that the RIA artificially conserves the vertical pseudo-momentum in contrast with higher order approximations (see Section 2.1).This might spuriously constrain the RIA statistics.Ray tracing (eikonal) models (discussed next) keep the scale separation hypothesis of the induced diffusion mechanism but without limitations on the interaction strength.

Eikonal Theories
With induced diffusion as an important contribution to the cascade in the GM spectrum, Henyey & Pomphrey (1983) kept the scale separation hypothesis but relaxed the weak-interaction approximation.They therefore used an eikonal approximation, describing propagation along rays in a random medium.The fluctuations of N have negligible effects and the shear is the only influent quantity (Henyey et al. 1986).The authors generated the shear by a random superposition of 50 (threedimensional) propagating waves, a Monte Carlo sampling of the GM spectrum.They numerically solved the classical ray equations for wave packet propagation in this background field.Whereas the horizontal wave vector does not change much, as expected, the vertical wavenumber k z grows in time by rapid excursions, like at the approach of critical levels (see Section 2.4).This strong interaction is in contrast with the progressive diffusion in the k z space described by the RIA (Equation 15).Henyey et al. (1986) have provided an empirical model of the Monte Carlo results, giving an expression for the energy dissipation.It has the same form as for the RIA theory but with a smaller coefficient, by a factor of 3 to 4. Ocean measurements as well as the numerical computations of Winters & D'Asaro (1997) indicate a result just in between these two estimates (Figure 8).Sun & Kunze (1999) have recently revisited the model of Henyey et al. (1986) and found a significant effect of the vertical divergence ∂u z /∂z in addition to the shear ∂u H /∂z, which was only considered previously.Their results are closer to the estimates of the energy dissipation rate in the GM spectrum and could be generalized to the parameterization of wave effects in other conditions.Hines (1991) applied a similar idea of "random Doppler spreading" by horizontal wind to the atmosphere, leading to k −3 z spectra.Souprayen et al. (2001) have similarly obtained a diffusion effect in vertical wavenumber, as in Equation 15.This diffusion does not result, however, from a resonant interaction hypothesis but from the assumption of a Gaussian shear with white noise derivatives, which prevent the trapping in critical layers.This contrasts with the realistic shear used in oceanic studies.Another difference from the oceanic studies is the existence of a systematic upward vertical propagation from a ground source, occurring simultaneously with the diffusion toward large vertical wave numbers.This model remarkably leads to a k −3 energy spectrum whose level is in N 2 but with a coefficient depending on the forcing conditions, unlike what would be expected for saturated spectra.Note that a wavenumber cutoff has to be introduced to ensure energy dissipation at small scale (a rough model of wave breaking) but it does not control the k −3 level, unlike in models of saturated spectra.Other features of the atmospheric spectra are also remarkably reproduced, and the relative simplicity of the calculation allows applications to parameterize the internal wave effects in the presence of realistic sources.The results rely, however, on a closure hypothesis, assuming that the shear is produced by the wave field, which is in principle inconsistent with the scale separation required for the validity of ray tracing.This is a general limitation of the eikonal approaches.

Models with Strong Wave Interactions
Quasi-normal closures were developed in the 1960s and 1970s, starting with the direct interaction approximation (DIA) model of Kraichnan, and they successfully describe the dynamics of isotropic turbulence, when intermittency effects can be neglected (see, e.g., Lesieur 1997).Extensions to anisotropic problems provide remarkable results, but with a tedious computational effort (see Godeferd et al. 2001 for a recent review).Although such models are designed for strongly nonlinear turbulence, they naturally include linear effects and should match with the weakly nonlinear theories, as discussed by Holloway (1980) and Müller et al. (1986).
The Eulerian approach seems appropriate, as large-particle displacements are expected, although Allen & Joseph (1989) made a valuable attempt at describing the wave field as a statistical equilibrium in the Lagrangian coordinates.This is probably unrealistic, but the authors provide interesting relations between Eulerian and Lagrangian variables for strong amplitudes.
The quasi-normal closures start with the same assumption (Equation 14) as for weakly nonlinear waves.However, the resulting equations lead to spurious negative spectra in ordinary turbulence.In an attempt to solve this problem, an eddy damping of the triple correlations is introduced.The resulting equation for the triple correlations V (3) = a i (k 0 )a j (k 1 )a k (k 2 ) is symbolically written as V (3) + LV (3) +µV (3) = V (2) V (2) , where L is the linear operator representing the free-wave propagation and µ is the rate of eddy damping.In the eddy damping quasi-normal model (EDQNM) discussed here, this damping is assumed to occur on a short timescale in comparison with the evolution of the spectra V (2) .In the RIA, a small damping of the triple correlation is also introduced (to damp the "memory" of the interactions, thus introducing the required time irreversibility), but the limit µ → 0 is taken, resulting in the selection of resonant triads.Therefore the RIA is the weak-amplitude limit of the EDQNM.
The application of EDQNM to the stratified case was first used by Carnevale & Frederiksen (1983) in the vertical plane.For simplicity, they further assumed the absence of correlations u z ρ (which is also satisfied in the RIA owing to the wave character of the individual modes).This is, however, a severe limitation, as this correlation (buoyancy flux) represents the mixing effects, which is of primary interest.Holloway (1988) introduced this correlation and found that the advective transfer to smaller scales is more efficient for potential energy than for kinetic energy.Potential energy thus accumulating at small scale is transferred to kinetic energy by the buoyancy flux u z ρ .Therefore, although the spectral integral of this flux is positive, expressing a global (down-gradient) mixing effect, it is positive only at large scale and negative (counter-gradient) at small scales: The small scales tend to restratify the background.This point of view is quite opposite to an earlier theory by Lumley (1964), who introduced buoyancy effects as a perturbation to a usual Kolmogoroff cascade, leading to a progressive change from the k −5/3 range to a k −3 range at smaller wave numbers.Although this theory may correctly describe stratified turbulence generated by shear effects, it does not match with the generation by a wave field at larger scales.
The extension of EDQNM to the three-dimensional case was performed by Godeferd & Cambon (1994), following the general procedure of Cambon (1989), using the wave-vortex decomposition (Equation 4).Staquet & Godeferd (1998) applied the model to the collapse of strong turbulence, initially isotropic, making precise comparisons with direct numerical simulations and laboratory experiments.The link with the weakly nonlinear wave regimes remains to be made.In particular the possibility of obtaining a k −3 range is an open question.Intermittency is quite strong in this regime, which may rule out quasi-normal theories.Powerful theoretical methods have been recently developed to describe how intermittency is generated in a passive field transported by a random Gaussian velocity field (see Shraiman & Siggia 2000 for a review).Similarly the large-scale weakly nonlinear waves randomly strain the "compliant" waves of the k −3 range, which could possibly be described by similar theoretical approaches.

CONCLUSIONS
Internal gravity waves offer a fascinating variety of flow phenomena.The linear propagation is already a rich subject, with various focusing and steepening effects.The final instability, leading to breaking itself, seems more universal, with a similar overturning event occurring in various cases.These processes are now fairly well known from theory, laboratory, or numerical experiments.In natural media, individual processes are however difficult to identify, and a statistical description is more appropriate.It is a necessity for parameterizing the mean effects of the internal waves in the atmospheric or oceanic circulation.
There is a plethora of statistical models for internal gravity waves in the ocean or the atmosphere, with some duplications between the two domains.Testing such theories on field measurements is difficult owing to the sparsity of the data and the lack of knowledge of the forcing conditions.Reproduction of natural internal wave spectra in well-controlled laboratory or numerical experiments is therefore a necessity.The reproduction of k −3 spectra in both laboratory and numerical experiments as well as of the GM spectrum in numerical simulations are encouraging steps in this direction.We still need to select good standard configurations for testing genuine statistical properties of the wave interactions, independent of the details of the forcing conditions.The classical grid turbulence quickly collapses in a stratified medium, and it seems always quite different from the "wave turbulence" sustained by internal gravity waves.
Many theoretical questions have been left unsolved.Even the classical weakinteraction theory is not so clearly understood, although it has been successful in interpreting the GM spectrum.The link with stronger nonlinear regimes at smaller scales could be made thanks to recent developments in anisotropic twopoint closures.It is not clear whether k −3 spectra can be explained by such theories.
Intermittency is much higher than in usual turbulence, which may rule out such closure models.Recent progress in "scalar turbulence," describing the cascade of a transported quantity under random large-scale straining, could provide tools for a rigorous approach to this problem.
However the available methods could already be useful for improving the parameterization of internal waves effects on global circulation.The difficulty is to match, in a consistent closure scheme, different models valid in widely different ranges of temporal and spatial scales.The generation process generally produces large-scale linear modes, as resonant scattering is important at intermediate scales, whereas eikonal approaches seem successful in describing smaller scales.We hope that this review will help in connecting these different approaches.

Figure 1
Figure 1 View of noctilucent clouds from Kustavi, Finland (61 • N, 21 • E) on 22 July 1989 showing characteristic bands and streak structures.In this case, bands are separated by ∼50 km and streaks by ∼3 to 5 km (from Fritts et al. 1993; photograph by Pakka Parviainen).

Figure 2
Figure 2 The instability of an internal gravity wave by parametric subharmonic instability.(a) Plane progressive wave obtained in the numerical experiments by C.R. Koudella & C.Staquet (submitted).The vorticity of the growing perturbation is displayed using grey levels (the primary wave field has been filtered out).(b) Case of a standing wave produced in a laboratory experiment byBenielli & Sommeria (1998).The growing wave packet is visualized by the deformation of horizontal dye strips.

Figure 5
Figure 5 Interaction of an internal gravity wave with a horizontal shear flow ūi x with either a vertical shear (a) or a horizontal shear (b), in a stably stratified fluid with constant Brunt-Väisälä frequency N (the Coriolis parameter is set to zero).The sketch is drawn in a vertical (x, z) plane in (a) and in a horizontal (x, y) plane in (b).Rays are plotted in the absolute reference frame (relative to which ω is measured) so that c g + ū(y)i x is the absolute group velocity of the wave.

Figure 7
Figure 7 Evidence of k −3 potential energy spectra E p (k) in the presence of internal gravity wave breaking.(a) Schematic spectra of vertical temperature gradients k 2 E p (k) in the oceanic thermocline (from Gregg 1987).(b) Stratospheric temperature spectra measured from instrumented balloons, showing a k −3 z range with a Kolmogoroff k z −5/3 tail at high k z (provided by F. Dalaudier).(c) Laboratory experiments of a breaking standing internal gravity wave, from Benielli & Sommeria 1996.(d ) Three-dimensional large eddy simulations of a forced gravity wave, from Carnevale et al. (2001); an analytical fit in k −3 at small k, ending in the Kolmogoroff spectrum k −5/3 at large k, is drawn.

Figure 8
Figure 8 Energy dissipation rate in the wave spectrum of Garrett & Munk (1979), as a function of the mean energy density (from Winters & D'Asaro 1997).Results from direct numerical computations are compared with the eikonal model of Henyey et al. (1986) (solid line), the calculations of McComas & Muller (1981) using resonant interaction approximation (dashed line), and oceanic measurements, from Gregg (1987), (grey line).

Figure 6
Figure 6 Evidence of convective instabilities resulting from internal gravity wave breaking.(a) Direct numerical simulation of a progressive monochromatic internal gravity wave of initial steepness s = 0.3.Six isopycnals are shown shortly after breaking has occurred (from Koudella 1999).(b) Convective mushrooms visualized by dye bands for a standing wave in a laboratory experiment (from Benielli & Sommeria 1998).(c) Same as (b) for a stably stratified shear layer (from Schowalter et al. 1994).