# Internal gravity waves in geophysical fluids

Chantal Staquet

► To cite this version:

Chantal Staquet. Internal gravity waves in geophysical fluids. Environmental Stratified Flows, Springer, pp.75-131, 2005, 10.1007/3-211-38078-7_2. hal-00264751

HAL Id: hal-00264751

https://hal.science/hal-00264751v1

Submitted on 8 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Chapter 2:
Internal gravity waves in geophysical fluids

Chantal Staquet

Laboratoire des Ecoulements Géophysiques et Industriels,
Université Joseph Fourier, Grenoble, France
Chantal.Staquet@hmg.inpg.fr
**Abstract** In this second Chapter, we recall some basic notions on stably-stratified flows before focusing on internal gravity wave dynamics. In Section 1, we illustrate the occurrence of stably-stratified flows in nature and, in Section 2, we derive the Boussinesq approximation. Except for a very brief presentation of the Kelvin-Helmholtz instability in Section 3, we start discussing about internal gravity waves from this section on. The linear properties of the wave are discussed in this Section. The mechanisms that can lead the wave field to breaking are addressed in Sections 4 to 6: parametric and buoyancy-induced instabilities (Section 4), interaction with a shear flow (Section 5), interaction with a sloping boundary (Section 6). We briefly discuss about the statistical properties of the breaking wave field in Section 7 and introduce some notions on mixing in Section 8.

# 1 Stably-stratified fluids: why are they important

## 1.1 Stably-stratified fluids are ubiquitous

An incompressible fluid is stably-stratified if it displays a vertical density gradient with negative sign (or a vertical temperature gradient with positive sign): the fluid becomes less dense, or warmer, as the altitude increases. Such fluids are ubiquitous in nature: the oceanic water masses, the stratosphere (which is the part of the atmosphere comprised between ≈ 10 and 50 kms) are stably-stratified. This is exemplified in the oceanic case in Figure 1: the temperature profile measured in the Mindanao trench during the 1929-1930 Snellius expedition in eastern Indonesia is displayed as a function of depth. The gradient is positive down to 3500m but compressibility effects become important below this depth leading to a slight warming of the water masses. Such a change of sign in the temperature gradient therefore does not imply that the fluid is unstably stratified but simply results from compression of the water columns. The concept of potential temperature, which is discussed in detail in Chapter 4 of the lecture notes, is introduced to account for compressibility effects. Indeed, the vertical gradient of the potential temperature field is now everywhere positive.

In this example, the huge depth of the temperature record, close to 10 km, along with the sparsity of the sampling, implies that the overall (large scale) behaviour of the temperature profile is captured. If one considers much smaller scales however, the vertical potential temperature gradient is not always positive. An example is provided in the stratosphere in Figure 2, which displays the potential temperature as a function of altitude. While the gradient is positive in overall, negative gradients exist locally, which are due to turbulent motions. In the stratosphere, the temperature and velocity fluctuations are mostly due to internal gravity waves so that this turbulence results from the breaking of the waves, as we shall further discuss it.

Stably-stratified media are also found in the interior of stars. In the Sun for instance, a stably stratified region exists below the convective zone, in the so-called radiative zone. As we shall show it, internal gravity waves can propagate in a stably-stratified medium. It has been proposed by Schatzman (1996) that the superposition of these waves results in a net transport of mass toward the interior of the Sun, which would account for the under-abundance of Lithium observed in this star: the Lithium would be transported
<figure><img src="image_2.png" /></figure>

**Figure 1.** Temperature and potential temperature profiles from data in the Mindanao trench collected during the 1929-1930 Snellius expedition around Indonesia.

toward the core where it would be burned.

Stably-stratified media may also be used within an industrial, or environmental, con-

<figure><img src="image_3.png" /></figure>

**Figure 2.** High resolution (20cm) temperature profile, just above the tropopause. Boxes marked with letters (a...g) contain examples of sheets. Two close-ups of sheets "b" and "e" are displayed with the three profiles measured on the main gondola. Ongoing mixing is clearly present in the sheet "b" (from Dalaudier et al., 1994).
<figure><img src="image_4.png" /><figcaption>Figure 3. A solar pond (from www.solarpond.utep.edu).</figcaption></figure>

text. The solar pond (Figure 3) is a nice illustration. The idea is to capture heat by a stable stratification. The principle is very simple: the bottom of a pond is made dark and very salty, in order for the bottom water, when heated by the sun (more precisely, when heated by the dark bottom) to remain heavier than the layer above. The temperature of the bottom water may reach 70 to 100 degrees and is used, for instance, to warm water inside a pipe that would cross the pond in its lowest part.

Stably-stratified fluids are also found in the cooling circuits of nuclear reactors.

## 1.2 Stably-stratified fluids give rise to phenomena with environmental implications

The solar pond just discussed provides one example that stably-stratified fluids may give rise to phenomena with environmental implications. Another common phenomenon, which has a negative impact on the environment, is the well-known situation of a thermal inversion in the atmosphere of a valley in winter. The process is again simple (f.i. Stull, 1988): the ground cools when no longer heated by the sun and, when its temperature becomes lower than that of the neighbouring air, the air cools in turn in transfering heat to the ground. A very cold layer of air thus exists at the bottom of the valley, which is surmounted by warmer air, leading to a stable (potential) temperature gradient. This stable situation may persist for the whole day in winter. If cities are located at the bottom of the valley, the pollutants emitted by road traffic and factories remain trapped in the stably-stratified layer of air, leading to high concentrations of pollutants. Two examples are provided in Figure 4, in the Grenoble valley and the Chamonix valley.

## 1.3 Quantification of the stability: the Brunt-Väisälä frequency

The Brunt-Väisälä frequency is the frequency of a fluid particle displaced adiabatically from its equilibrium position along the vertical. The demonstration is very standard. Let
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) The valley of Grenoble in winter, when a thermal inversion exists all day long. (b) Photograph of the Chamonix valley, with the Green Needle (Aiguille Verte) in background. Water and pollutant particles are again trapped in a stable thermal layer at the bottom of the valley.</figcaption></figure>

us do a thought experiment and consider one fluid particle of unit volume, at equilibrium position $z$ in a stably-stratified fluid with density profile $\bar{\rho}(z)$. When the particle is displaced adiabatically by a vertical distance $\zeta$, the local density at the new position is $\bar{\rho}(z + \zeta)$. The fluid particle is thus subjected to two forces, its weight $g\bar{\rho}(z)$ and the Archimedean force $-g\bar{\rho}(z + \zeta)$. Applying Newton's law to the fluid particle results in the equation $\ddot{\zeta} + N^2\zeta = 0$, where $N^2 = -\frac{g}{\bar{\rho}(z)}\frac{d\bar{\rho}}{dz}$ is the Brunt-Väisälä, or buoyancy, frequency. Stability implies that $N^2$ is positive.

A standard approximation to describe stably-stratified incompressible fluid motions is the Boussinesq approximation, described in Section 2.2; in this case, the density and the temperature profiles are linearly related. Therefore, the Brunt-Väisälä frequency may also be defined as $N^2 = \frac{g}{\bar{T}(z)} \frac{dT}{dz}$. When the fluid is compressible, the potential temperature is used instead of the temperature and the expression of the Brunt-Väisälä frequency is the same as the previous one with the temperature profile being simply replaced by the potential temperature profile.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Atmospheric lee waves generated by a wind blowing over two islands.</figcaption></figure>

## 1.4 What are the generic motions?

As announced above, the generic fluid motions in a stably-stratified fluid are internal gravity waves. Indeed, the sum of the weight and of the Archimedean force acting on a fluid particle displaced from equilibrium position is (at first order in $\zeta$) $-g\zeta \frac{d\bar{\rho}}{dz}$, which is a restoring force. This is the restoring buoyancy force. A restoring force in a continuous medium generates wave motions. We shall show in Section 3 that linear waves are indeed a solution of the equations of motions. Before doing so, we provide two illustrations of internal gravity waves in the atmosphere (Figure 5) and in the ocean (Figure 6).

In the atmosphere, the main mechanisms for gravity wave generation is the wind blowing over topography. The generated waves are called lee waves and are made visible by clouds. Figure 5 thus displays the cloud cover over the sea when the blowing wind interacts with two islands. The island located on the left of the picture is too high for the air advected by the wind to go over the island and the advected air simply flows around the island. In this case, the fluid particles are not lifted and no wave is generated. Instead, a quasi-two-dimensional pattern (in a horizontal plane) of vortices are produced, of the Karmann street type. By contrast, the island located on the right of the picture is low enough for the advected air to go over the island and a quite different pattern is observed: in this case, waves are generated, with horizontal wave length of the order of the island size.

In the ocean, the two main mechanisms for the generation of internal waves are the wind and the tide. The wind forcing is indirect: the wind mixes the upper layer of the ocean (over a depth of order 100m) and motions at the bottom of this layer perturbs the stably-stratified fluid below which generates internal gravity waves. By contrast, the tide generates internal gravity waves directly, through its interaction with topography (seamounts, continental slopes, etc.) Internal gravity waves generated in the neighbourhood of a continental slope by the tide are displayed in Figure 6. The waves are in the interior of the water masses but, because their amplitude is large and they propagate not far from the surface, the velocity field they induce perturbs the free surface which makes them visible using a radar on board of a satellite.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Internal waves (propagating below the free surface) generated by the interaction of the tide with a continental slope. This is a radar image acquired by the Synthetic Radar Aperture aboard of a satellite.</figcaption></figure>

## 1.5 What about unstably stratified flows?

Unstably stratified flows are also very common in natural media. In this situation, no waves are generated since the buoyancy force is not a restoring one. On the opposite, this force lifts the fluid particle further away from the position it has been brought. The commonest example occurs in the atmospheric boundary layer: the ground is heated by the Sun during the day, which heats the air (the direct heating of air by the Sun is much less important, especially if the air is dry). The vertical gradient of the potential temperature profile is therefore negative in this case. Such a situation is highly unstable and gives rise to large scale turbulent motions, which mix the fluid. Therefore, if pollutants are emitted at the ground level, the local concentrations are generally smaller than in a stable situation, but may be found very far from their emission region. For instance, in summer in the Chamonix valley, high concentrations of heavy metals (such as Mercury) have been found at 4300m altitude, at Dôme du Goûter (Veyssyre et al., 2000).

# 2 Basic equations

## 2.1 Governing equations

Let $(x, y, z)$ be a cartesian coordinate system in an absolute reference frame where the fluid velocity is $\vec{u} = (u, v, w)$. The fluid motions in this reference frame are governed by the compressible Navier-Stokes equations. These equations are composed of

* the momentum equation, which is, when expressed per unit mass:

$$ \frac{D\vec{u}}{Dt} = -\nabla p + \vec{g} - 2\vec{\Omega} \wedge \vec{u} + \frac{\mu}{\rho} \nabla^2 \vec{u}, \quad (2.1) $$
where $\frac{D}{Dt} = \frac{\partial}{\partial t} + \vec{u} \cdot \nabla$ refers to the material derivative, $p$ is the pressure field, $-2\vec{\Omega} \wedge \vec{u}$ is the Coriolis acceleration due to earth rotation and $\mu$ is the dynamical viscosity;

* the mass conservation equation:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{u}) = 0, \qquad (2.2) $$

where $\rho$ is the density field;

* and a dynamical equation for the internal energy $e$, stemming from the first law of thermodynamics applied to a system made of a unit mass of volume $v = 1/\rho$:

$$ \frac{De}{Dt} = Q - p \frac{Dv}{Dt}. \qquad (2.3) $$

$Q$ is the amount of heat added to the system $v$ and $-p \frac{Dv}{Dt}$ is the work of pressure forces on the system.

## 2.2 The Boussinesq approximation

Water is a very weakly compressible fluid so that oceanic fluid motions may be assumed to be incompressible over the typical vertical scale of the motions. The fact that density, that is the mass of the fluid system, still varies in space and time may be incompatible with the assumption that the fluid velocity is incompressible. A consistent, and very well-known, approximation that satisfies the former requirement while assuming that the fluid is close to incompressibility is the Boussinesq approximation. We now describe in detail this approximation (see also Cushman-Roisin, 1994, p. 37).

In the following, we set $\rho = \rho_0 + \tilde{\rho}$ where $\rho_0$ is some constant reference density and $\tilde{\rho}$ refers to the deviation field about $\rho_0$. In the ocean for instance, $\rho_0$ is the mean value of the density, equal to $1028 \text{ kg/m}^3$ (at atmospheric pressure) and $\tilde{\rho}/\rho_0 \simeq 10^{-3}$. We also introduce reference values for the temperature and the salinity: $T_0 = 10^\circ \text{C}$ and $S_0 = 34 \text{ g/kg}$.

**How does the mass conservation equation become $\nabla \cdot \vec{u} = 0$?** Let us find in which sense the mass conservation equation can be approximated by $\nabla \cdot \vec{u} = 0$.

We introduce the small parameter $\epsilon$ such that $\tilde{\rho}/\rho_0 = \epsilon \rho^{(1)}/\rho_0$, with $\rho^{(1)}/\rho_0 = O(1)$. We also decompose $\vec{u}$ into $\vec{u}^{(0)} + \epsilon \vec{u}^{(1)}$; hence, $\vec{u}^{(0)}$ is the velocity field in the limit $\epsilon \to 0$, that is, in the limit of infinitely small density fluctuations. Introducing this decomposition of the fields into the mass conservation equation, one gets

$$ \rho_0 \nabla \cdot \vec{u}^{(0)} + \epsilon \left( \frac{D}{Dt} \rho^{(1)} + \rho_0 \nabla \cdot \vec{u}^{(1)} + \rho^{(1)} \nabla \cdot \vec{u} \right) = 0. \qquad (2.4) $$

At order zero in $\epsilon$, the mass conservation reduces to $\nabla \cdot \vec{u}^{(0)}$. This implies that the momentum equations for an incompressible fluid are also valid in this limit (that is, with $\vec{u} \simeq \vec{u}^{(0)}$). This also implies that the gravity acceleration $g$ should be large enough for the product $g\tilde{\rho}/\rho_0$ to stay finite.
**How does the energy equation become the usual heat equation?** For any constant volume transformation, the internal energy $e$ is equal to $C_v T$ where $C_v$ is the heat capacity per unit mass and $T$ is the temperature. Equation (2.3) can thus be rewritten as

$$C_v \frac{DT}{Dt} = Q + \frac{p}{\rho^2} \frac{D\rho}{Dt}. \qquad (2.5)$$

Since the state variable $T$ has been introduced, we need an equation of state. In the limit of very small density deviation about a constant reference density, the fluid may be assumed to be incompressible, as just seen. Hence, $\rho$ only weakly depends on pressure. Density fluctuations arise from temperature and salinity changes about the background state so that, at first approximation, a linear equation of state may be assumed

$$\rho = \rho_0[1 - \alpha(T - T_0) + \beta(S - S_0)], \qquad (2.6)$$

where $\alpha$ is the thermal expansion coefficient and $\beta$ is the haline contraction coefficient.

Assuming that the heat added to the fluid volume $v$ results from diffusion process, $Q$ can be related to $T$ using Fourier law: $\rho Q = \kappa \nabla^2 T$. Replacing $Q$ by this expression in equation (2.5) and using $\frac{D\rho}{Dt} = -\rho \nabla \cdot \vec{u}$, one gets ($k$ is the thermal conductivity of the fluid)

$$\rho C_v \frac{DT}{Dt} + p \nabla \cdot \vec{u} = k \nabla^2 T. \qquad (2.7)$$

In the limit $\tilde{\rho} \ll \rho_0$, that is, $\vec{u} \simeq \vec{u}^{(0)}$ one eventually recover the usual heat equation

$$\frac{DT}{Dt} = \frac{k}{\rho_0 C_v} \nabla^2 T. \qquad (2.8)$$

The coefficient $\frac{k}{\rho_0 C_v}$ is the thermal diffusivity in the fluid, usually denoted $\kappa$. If the density fluctuation $\tilde{\rho}$ mainly arises from temperature fluctuations, the above equations become, using (2.6):

$$\frac{D\tilde{\rho}}{Dt} = \kappa \nabla^2 \tilde{\rho}. \qquad (2.9)$$

**And what about the momentum equation?** We use the assumptions $\tilde{\rho} \ll \rho_0$ and $\rho \simeq \rho_0$. In the momentum equations, $\rho$ comes into play in two terms of quite different physical meaning: (i) in the acceleration and Coriolis terms, $\rho D\vec{u}/Dt$ and $-2\rho\vec{\Omega} \wedge \vec{u}$, and (ii) in the gravity term $g\rho$. While $\rho$ can be replaced by $\rho_0$ in the former two terms, it certainly cannot be replaced also in the latter term because the fluid would be homogeneous otherwise, with constant density $\rho_0$. Therefore

* $\rho D\vec{u}/Dt \rightarrow \rho_0 D\vec{u}/Dt$; $-2\rho\vec{\Omega} \wedge \vec{u} \rightarrow -2\rho_0\vec{\Omega} \wedge \vec{u}$.

* $-\nabla p + \rho\vec{g} \rightarrow -\nabla\tilde{p} + \tilde{\rho}\vec{g}$, decomposing the dynamical pressure $p$ into $P_0$ (hydrostatic pressure) plus $\tilde{p}$ and using the hydrostatic balance of the rest state $dP_0/dz = -\rho_0 g$.

The equations of motions become

$$\rho_0 \frac{D\vec{u}}{dt} = -\nabla\tilde{p} + \tilde{\rho} \vec{g} - 2\rho_0 \vec{\Omega} \wedge \vec{u} + \mu \nabla^2 \vec{u}. \qquad (2.10)$$
The Boussinesq equations are equations (2.9) and (2.10) above along with the incomp-
ressibility condition $\nabla \cdot \vec{u} = 0$. We recall that this approximation is strictly valid in the
limit $\vec{u} = \vec{u}^{(0)}$ and $\tilde{\rho} \ll \rho_0$.

## 2.3 Non-dimensional parameters

Models of geophysical (or astrophysical) flows are often designed in the laboratory. The idea is to model not the whole system but only a part of it, whereof dynamics are to be studied. For instance, the formation and dynamics of the Great Red Spot of Jupiter has been modelled in a laboratory experiment using two counter-rotating disks, in order to reproduce the zonal shear that prevails on this planet (as in any strongly rotating planet) and controls the dynamics. The scales are very different in the natural and laboratory systems but comparison can be performed by considering the non-dimensional equations of motions. These non-dimensional equations involve a few non-dimensional parameters and the application of the laboratory results to the natural system can be done if the non-dimensional parameters are the same. The other interest of working with non-dimensional equations, especially when natural systems are to be studied, is that one deals with values of order 1. Such a non-dimensionalization is actually necessary when the equations of motions are solved on a computer.

We now list the main non-dimensional parameters that govern the dynamics of a
Boussinesq flow (that is, a flow whose dynamics can be described by the Navier-Stokes
equations in the Boussinesq approximation) subjected to the Coriolis force.

Let $U$ be the scale of the velocity field, $L$ be a typical length scale of the fluid motions and $\tau$ the time scale. Three non-dimensional parameters compare the nonlinear time scale $L/U$ to the time scale $\tau$, depending upon the physical process whereof $\tau$ is the time scale: (i) if $\tau$ is the time scale of viscous effects, then $\tau = L^2/\nu$ and $Re = \frac{\tau}{L/U} = \frac{UL}{\nu}$ is the Reynolds number; (ii) if $\tau$ is the time scale of stratification effects, then $\tau = N^{-1}$ and $Fr = \frac{\tau}{L/U} = \frac{U}{NL}$ is the Froude number; (iii) if $\tau$ is the time scale of rotation (or Coriolis) effects, then $\tau = (2\Omega)^{-1}$ and $Ro = \frac{\tau}{L/U} = \frac{U}{2\Omega L}$ is the Rossby number. Finally, the Prandtl number compares the viscosity to the thermal diffusivity: $Pr = \frac{\nu}{\kappa}$.

A few other non-dimensional parameters can be inferred from these four ones. These
are for instance the Ekman number, which compares viscous and Coriolis effects: $Ek = \frac{\nu}{2\Omega L^2} = ReRo$; and the ratio of rotation to stratification effects (also referred to as the
Prandtl ratio): $\frac{N}{2\Omega} = \frac{Ro}{Fr}$.

## 2.4 Linearization about a basic state

Since we are interested in motions of small amplitude relative to a basic state, it is
useful to decompose the fields into a part associated with the basic state and a part
associated with the small amplitude motions. Thus we write: $\vec{u} = \vec{U} + \vec{u}'$, $\tilde{\rho} = R +$
$\rho'$, $\tilde{p} = P + p'$ where $\vec{U}$ and $R$ may depend upon the spatial coordinates but (this is
our assumption) do not depend upon time, R and P are in hydrostatic balance, and the amplitude of $\vec{u}'$, $\rho'$ and $p'$ is much smaller than that of U, R and P respectively. Introducing this decomposition into the Boussinesq equations and linearizing about the basic state yield the linearized Boussinesq equations. We shall consider two specific cases.

• The basic state is a fluid at rest in hydrostatic balance.

In this case, $U = 0$ and, therefore, the amplitude of $\vec{u}'$ is much smaller than 1. The linearized Boussinesq equations become in this case ($\nu$ is the kinematic viscosity):

$$ \frac{\partial \vec{u}'}{\partial t} = -\frac{1}{\rho_0} \nabla p' - g \frac{\rho'}{\rho_0} \vec{i}_z - 2\vec{\Omega} \wedge \vec{u}' + \nu \nabla^2 \vec{u}' \quad (2.11) $$

$$ \frac{\partial \rho'}{\partial t} = -w' \frac{dR}{dz} + \kappa \nabla^2 \rho' \quad (2.12) $$

$$ \nabla \cdot \vec{u}' = 0. \quad (2.13) $$

• The basic state is a parallel shear flow $U(z)$ with density field $R(z)$ in hydrostatic balance.

The linearized Boussinesq equations of motions become in this case

$$ \frac{\partial u'}{\partial t} + U \frac{\partial u'}{\partial x} = -w' \frac{dU}{dz} - \frac{1}{\rho_0} \frac{\partial p'}{\partial x} + 2\Omega v' + \nu \nabla^2 u' \quad (2.14) $$

$$ \frac{\partial v'}{\partial t} + U \frac{\partial v'}{\partial x} = -\frac{1}{\rho_0} \frac{\partial p'}{\partial y} - 2\Omega u' + \nu \nabla^2 v' \quad (2.15) $$

$$ \frac{\partial w'}{\partial t} + U \frac{\partial w'}{\partial x} = -\frac{1}{\rho_0} \frac{\partial p'}{\partial z} - g \frac{\rho'}{\rho_0} + \nu \nabla^2 w' \quad (2.16) $$

$$ \frac{\partial \rho'}{\partial t} + U \frac{\partial \rho'}{\partial x} = -w' \frac{dR}{dz} + \kappa \nabla^2 \rho' \quad (2.17) $$

$$ \nabla \cdot \vec{u}' = 0. \quad (2.18) $$

# 3 Solutions of the linearized Boussinesq equations: internal gravity waves, Kelvin-Helmholtz instability

## 3.1 Generalities on linear waves

A linear wave in a homogeneous and steady medium<fnref n="1" /> is characterized by a frequency $\Omega$ and a wave vector $\vec{k}$. For a plane wave in an infinite medium, the motion is of the form

$$ A e^{i(\vec{k} \cdot \vec{x} - \Omega t)}, \quad (3.1) $$

<footnote n="1">¹The properties of the medium remain unchanged by translation of the reference frame and do not depend upon time.</footnote>
where A is the wave amplitude. The properties of the wave are completely set by the dispersion relation, namely a relation between $\Omega$ and $\vec{k}$:

$$ \Omega = F(\vec{k}), \qquad (3.2) $$

where the function F also involves the fluid parameters upon which the wave propagation depends (N, g, etc.) The dispersion relation does not depend upon space and time because of the homogeneity and steadiness assumptions.

Since the dispersion relation characterizes the wave field, its derivation is an essential step in the wave study. One can distinguish two methods for this purpose. In the classical method, one first needs to derive the linearized equations of motion. Since the medium is homogeneous, the coefficients of these equations are constant. The general form of the wave solution is then substituted in these equations, from where the dispersion relation is obtained.

A general method that gives access directly to the dispersion relation, without using the equations of motion, can be used instead (Lighthill, 1978). Because it is simpler than the classical method and very powerful, we shall rely on this second method to compute the dispersion relation.

## 3.2 A general method to compute the dispersion relation

The frequency $\Omega$ of an oscillating system, whether consisting of a finite number of degrees of freedom (like the pendulum and or chain of springs) or evolving in a continuous medium, can be computed directly in three different ways (Lighthill, 1978):

1. From the restoring force. The square of $\Omega$ is then obtained as follows:

$$ \Omega^2 = \frac{\text{restoring force}}{\text{displacement x mass}}; \qquad (3.3) $$

2. From energetic considerations. $\Omega^2$ is then defined as:

$$ \Omega^2 = \frac{\text{generalized stiffness}}{\text{generalized inertia}}, \qquad (3.4) $$

where the *generalized stiffness* is the coefficient of the displacement squared in the expression of the potential energy and the *generalized inertia* is the coefficient of the temporal derivative of the displacement squared in the expression of the kinetic energy;

3. From dimensional analysis.

Before applying this general method to internal gravity waves, let us show how it works for a simple discrete system. For this purpose, we consider the harmonic oscillator, whose frequency is $\Omega = \sqrt{\frac{K}{m}}$, where K is the stiffness of the spring and m the mass of the object attached to the spring.

In case 1., the restoring force is $-K.x$ where x is the displacement of the mass from equilibrium position. Hence, $\Omega^2$ should be equal to $|K.x|$ divided by $|x|$ times m, which yields the actual expression for $\Omega$. In case 2., the potential energy is $0.5Kx^2$ so that
the “generalized” stiffness is $K$ and the kinetic energy is $0.5m\dot{x}^2$ so that the generalized inertia is $m$. Finally, in case 3., let us assume that $\Omega$ is a function of $K$ and $m$, of the form $K^\alpha m^\beta$. By writing that the unit of $\Omega$, $m$ and $K$ are rad.s$^{-1}$, kg (also Newton.m$^{-1}$.s$^2$) and Newton.m$^{-1}$ respectively, one gets $\alpha = -\beta = 1/2$.

### 3.3 Internal gravity waves in an infinite medium

**The dispersion relation from general principles.** In the following, we assume that the stably-stratified medium is homogeneous, which implies that $N$ is constant. We shall assume that the medium is steady throughout the paper.

To compute $\Omega$, we shall use the second expression provided above (case 2.). We first need to compute the potential energy to infer the generalized stiffness. The potential energy is the work done by the restoring force from the current position to the equilibrium position. We showed in section 1.4 that the restoring force is $g\frac{d\bar{\rho}}{dz}\zeta$ so that

$$E_p = \int_{\zeta}^{0} g \frac{d\bar{\rho}}{dz} zdz = -N^2 \rho_0 \int_{\zeta}^{0} zdz = 0.5 \rho_0 N^2 \zeta^2.$$ The generalized stiffness is therefore $\rho_0 N^2$. Let us now compute the kinetic energy to get the generalized inertia. As we show it in the next section, the motion induced by a monochromatic internal gravity wave occurs in a two-dimensional vertical plane, made by $\vec{k}$ and $\vec{g}$ (this is the propagation plane of the wave). We assume that this plane coincides with the $(x, z)$ plane. Hence, the kinetic energy is $E_k = 0.5\rho_0(u^2 + w^2)$. From continuity equation: $k_x u + k_z w = 0$, implying that $u = -\tan\theta w$, where $\theta$ is the angle of the wave vector with the horizontal.

Since, $w = \dot{\zeta}$, one gets $E_k = 0.5\rho_0(1 + \tan^2\theta)\dot{\zeta}^2 = 0.5\rho_0\frac{\dot{\zeta}^2}{\cos^2\theta}$. The generalized inertia is therefore $\frac{\rho_0}{\cos^2\theta}$.

Hence, the frequency (and dispersion relation) of a monochromatic internal gravity waves in a homogeneous medium is

$$\Omega^2 = N^2 \cos^2 \theta. \qquad (3.5)$$

It follows that internal gravity waves are *dispersive waves*, because the phase velocity $c = \Omega/|\vec{k}|$ depends upon $\vec{k}$, with an *anisotropic* dispersion relation, because $\Omega$ depends upon the angle that $\vec{k}$ makes with the horizontal.

**Main properties of internal gravity waves.** The properties of internal gravity waves are described in several textbooks (Lighthill, 1978; Leblond and Mysak, 1978; Cushman-Roisin, 1994) and we refer to these books for a detailed account. We make a brief summary of the important properties which we illustrate with laboratory experiments.

We cannot avoid discussing the beautiful experiment of Mowbray and Rarity (1967), which was designed to check the dispersion relation (3.5). A horizontal cylinder oscillates vertically at a constant frequency $\Omega$ in a constant $N$ fluid. The vertical motion of the cylinder perturbs the flow and relation (3.5) show that, if $\Omega^2 \le N^2$, waves will be produced. Linear waves in a non dissipative medium carry the momentum and the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Laboratory experiments performed by Mowbray and Rarity (1967). Visualization is made in a vertical plane and the dark vertical bar is the cylinder support.</figcaption></figure>

energy of their source without material transport (f.i. Andrews et al., 1987). The energy is transported at the group velocity $\partial\Omega/\partial\vec{k}$. For internal gravity waves, the group velocity is perpendicular to $\vec{k}$. Hence the energy should be carried away by the waves along directions making an angle $\theta$ with the vertical. This peculiar pattern is displayed in Figure 7 and is referred to as the St-Andrews cross.

Perhaps the most important property of internal gravity waves in a homogeneous stratified medium is that only a time scale is imposed, namely $N^{-1}$. There is no scale selection by the stratified medium<fnref n="2" />. It follows that in laboratory experiments of internal gravity waves forced by a paddle, the size of the paddle (or, more generally, of the forcing mechanism) selects the scale (f.i. McEwan, 1971). In the experiment by Mowbray and Rarity discussed above, the diameter of the cylinder selects the width of the energetic region. It may thus be expected that in experiments where the forcing mechanism at frequency $\Omega$ has no length scale, all wave vectors which make an angle $\cos^{-1}(\Omega/N)$ with the horizontal will be excited, according to the dispersion relation (3.5). In the laboratory experiments of McEwan and Robinson (1975) and Benielli and Sommeria (1998), internal gravity waves were forced by parametric instability: a tank stably-stratified with salted water oscillates vertically. In a frame of reference attached to the tank, the gravity appears to oscillate which triggers a parametric instability which generates in turn internal gravity waves. Though the forcing mechanism does not impose any length scale, it eventually appears that the size of the tank provides the scale selection: only the largest wavelength was excited.

When Coriolis effects are important, the waves are named “inertia-gravity waves”, with dispersion relation $\Omega^2 = N^2\cos^2\theta + f^2\sin^2\theta$ where $f$ is twice the projection of the angular vector of the Earth onto the local vertical. In a laboratory experiment, $f$ is

<footnote n="2">²By contrast, in the stratified shear flow problem for instance, the basic flow is made of a velocity profile of the form $U \tanh(z/\delta)$, which imposes both a velocity scale ($U$) and a length scale ($\delta$), so that a time scale ($\delta/U$) is imposed as well.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** (a) Laboratory experiment of a vertically oscillating sphere in a stably-stratified rotating medium with constant *N* and *f* (half a domain is shown due to symmetry). (b) Verification of the dispersion relation (3.6) for different values of *N* and *f*. Note that *θ* denotes here the angle that the wave vector makes with the vertical and that the intrinsic wave frequency is denoted as *ω* (from Peacock and Tabaei, 2005).</figcaption></figure>

twice the angular frequency of the rotating tank. It follows that Ω now becomes double-bounded: *f*² ≤ Ω² ≤ *N*². The dispersion relation can also be written as

$$ \cos^2\theta = \frac{\Omega^2 - f^2}{N^2 - f^2} \qquad (3.6) $$

so that, again, one angle should be selected in a rotating stably-stratified medium by a forcing mechanism with frequency Ω. The laboratory experiment by Peacock and Tabaei (2005) nicely illustrates this relation (Figure 8).

## 3.4 The Kelvin-Helmholtz instability in a stably-stratified fluid

**Evidence for the Kelvin-Helmholtz instability in natural media.** The Kelvin-Helmholtz instability is also described by the linearized Boussinesq equations, when the basic flow is made of a velocity profile of the hyperbolic (or error function) type in a stably-stratified medium.

The Kelvin-Helmholtz instability develops at the interface between two fluids flowing horizontally at different velocities. This configuration is strongly unstable when the two fluids have the same density: a perturbation amplifies as soon as the Reynolds number associated with the shear flow is of order 1 or larger. As the perturbation amplifies, nonlinear effects come into play which make the instability saturate. The nonlinear development of the instability is manifested as Kelvin-Helmholtz vortices. Examples for this instability are ubiquitous in geophysical flows, whether in the atmosphere (Figure 9) or in the ocean (Figure 10). In the former case, the instability results from a vertically sheared wind; in latter case, the instability is due to tidal motion over a sill.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Kelvin-Helmholtz instability visualized by clouds on top of a mountain.</figcaption></figure>

**The Richardson number.** A shear flow may exist in a constant N stably-stratified medium but a configuration most encountered in real flows is the one when the two fluids have also different densities (the lower fluid being of course heavier), leading to a vertical density profile of the hyperbolic tangent type, like the velocity profile. In this case, the density difference between the two fluids comes into play and the Froude number is no longer the most appropriate parameter to account for stratification effects. A bulk (or global or overall) Richardson number is rather introduced, whose definition comes naturally when the Boussinesq equations are made non-dimensional. The Richardson number is a measure of the stabilizing effect of buoyancy relative to the destabilizing effect of the shear.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Kelvin-Helmholtz instability due to tidal motion over a sill (from Farmer and Armi, 1999).</figcaption></figure>
Let us start from the Boussinesq equations linearized about a basic state consisting of the stratified shear flow $(U(z), R(z))$ (see the final set of equations at the end of section 2). We nondimensionalize these equations by choosing the following scales: $U$, the velocity scale, is half the velocity difference across the shear flow; $\delta$, the length scale, is the thickness of the shear flow; $\Delta\rho$, the density scale, is half the density difference across the shear flow. $T$, the time scale, is equal to $\delta/U$. It follows that the pressure scale is $\rho_0U^2$.

The non-dimensional equations are (using the same notation for dimensional and dimensionless variables):

$$ \frac{\partial u'}{\partial t} + U \frac{\partial u'}{\partial x} = -w' \frac{dU}{dz} - \frac{\partial p'}{\partial x} + \frac{1}{Re} \nabla^2 u' \quad (3.7) $$

$$ \frac{\partial v'}{\partial t} + U \frac{\partial v'}{\partial x} = - \frac{\partial p'}{\partial y} + \frac{1}{Re} \nabla^2 v' \quad (3.8) $$

$$ \frac{\partial w'}{\partial t} + U \frac{\partial w'}{\partial x} = - \frac{\partial p'}{\partial z} - J \rho' + \frac{1}{Re} \nabla^2 w' \quad (3.9) $$

$$ \frac{\partial \rho'}{\partial t} + U \frac{\partial \rho'}{\partial x} = - \frac{dR}{dz} w' + \frac{1}{Pr Re} \nabla^2 \rho' \quad (3.10) $$

$$ \nabla \cdot \vec{u}' = 0. \quad (3.11) $$

These equations involve three non-dimensional parameters: the Reynolds number, $Re = \frac{U\delta}{\nu}$, the Prandtl number and the bulk Richardson number $J = g \frac{\Delta\rho L}{\rho_0 U^2}$.

**The Richardson number criterion.** When the two fluids have different densities such that the lower fluid is heavier, the instability cannot develop whatever the density difference: as fluid particles raise because of the instability, kinetic energy is converted into potential energy so that, intuitively, a too strong density difference prevents the instability from developing. Chandrashekar (1961) derived an heuristic criterion for local instability based upon energetic considerations, which was demonstrated rigorously by Howard (1961) (see also Drazin and Reid, 1981, p. 326). The energetic considerations compare kinetic and potential energies, that is nonlinear and stratification effects, which relates to the Richardson number. Since the bulk Richardson number does not provide any information about the local stability of the shear flow, a local Richardson number was introduced for this purpose:

$$ Ri(z) = - \frac{g \frac{dR}{dz}}{\rho_0 \left( \frac{dU}{dz} \right)^2}. \quad (3.12) $$

The stability criterion is the following: a necessary condition for instability is $Ri(z) < 1/4$ somewhere in the fluid. It has been shown by Hazel (1972) that this condition is necessary and sufficient for a velocity and a density profile of the hyperbolic tangent type.
# 4 Instability of internal gravity waves in a constant $N$ fluid at rest

## 4.1 Resonant interaction theory

Resonant interactions in the fluid dynamical context were discovered independently by Phillips (1960) and Hasselmann (1962) for surface gravity waves and progressively extended to internal gravity waves by considering, in particular, the interaction of two surface waves with an internal gravity wave (Ball, 1964); the case of an internal wave triad was eventually considered by Phillips (1966) and Hasselmann (1967). Extensive work was next performed in the oceanographic context (see the reviews by Phillips, 1981; Muller et al., 1986; Staquet and Sommeria, 2002). It was shown in particular that energy transfers in the spectrum of large scale, low frequency internal gravity waves in the ocean (the Garrett-Munk spectrum) mostly consist of resonant interactions (McComas and Bretherton, 1977). In the atmosphere, the decrease of density with altitude amplifies the waves, so that resonant interactions are only rarely expected (a significant fraction of the wave energy is still backscattered downward, see f.i. Barat and Cot, 1992). Perhaps one exception is the Antartic atmospheric boundary layer, which is surmounted by a very cold layer of air in which the waves may be trapped and interact resonantly (Rees et al., 2000).

We briefly summarize this theory in the following, which we illustrate by results from Koudella and Staquet (2005).

**Wave steepness.** In the following, we consider a monochromatic internal gravity wave ($\vec{k}, \Omega$). We show below that fluid motions may be assumed to occur in the propagation plane ($\vec{k}, \vec{g}$) of the wave if the wave amplitude is infinitely small. Henceforth we assume two-dimensional fluid dynamics in the vertical ($x, z$) propagation plane. Let

$$ (u, w) = (A_u, A_w) \exp [ I(\vec{k} \cdot \vec{x} - \Omega t) ] + \text{c.c.} \quad (4.1) $$

denote the velocity field induced by the wave, where $\Omega = F(\vec{k})$ obeys the dispersion relation (3.5).

The normalized amplitude, or steepness, of the wave is defined as:

$$ s = \frac{u_{max}}{c_x}, \quad (4.2) $$

where $c_x$ is the phase velocity in the $x$-direction. As illustrated in Figure 11, this definition implies that isopycnals (i.e. constant density surfaces) are nowhere overturned when $s < 1$ (the wave is also said statically stable) and that the isopycnals are locally overturned when $s > 1$ (statically unstable wave). The case $s = 1$ implies that vertical isopycnals have locally formed. Since $u_{max} = A_u$ and $c_x = \frac{\Omega}{k_x} = \frac{N}{|\vec{k}|}$, the steepness is also defined as

$s = \frac{A_u |\vec{k}|}{N}$ for a monochromatic wave. Hence, taking $A_u$ as a velocity scale for the wave, $|\vec{k}|$ as a length scale and $N^{-1}$ as a time scale, the steepness is simply the wave Froude
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Illustration of the physical meaning of the steepness *s* of a monochromatic wave, defined by (4.2). A constant density line is drawn for different values of *s*.</figcaption></figure>

number. Let us show that fluid motions may be assumed to occur in the propagation plane ($\vec{k}, \vec{g}$) of the wave for an infinitely small wave amplitude. In three dimensions, the equations of motion for the vorticity $\vec{\xi} = \nabla \times \vec{u}$ are expressed by:

$$ \frac{\partial \vec{\xi}}{\partial t} + \vec{u} \cdot \nabla \vec{\xi} = \vec{\xi} \cdot \nabla \vec{u} - \nabla \times \frac{g}{\rho_0} \vec{\rho} \vec{e}_z + \nu \nabla^2 \vec{\xi}. \quad (4.3) $$

The magnitude of the advective term $\vec{u} \cdot \nabla \vec{\xi}$ or of the vortex stretching/tilting term $\vec{\xi} \cdot \nabla \vec{u}$ relative to the magnitude of the baroclinic term $\nabla \times \frac{g}{\rho_0} \vec{\rho} \vec{e}_z$ scales as the square of the wave steepness. This shows that, when $s^2 \ll 1$, the dynamics are not only weakly nonlinear but the fluid motions are mostly confined in a vertical plane if the vertical vorticity is initially zero. The linear stability analysis of a monochromatic internal gravity wave to three-dimensional perturbations by Klostermeyer (1991), later discussed in this paper, actually showed that this result holds as long as $s \lesssim 0.7$.

**Resonance relations.** The resonant interaction theory describes the interaction among a triad of internal gravity waves of infinitely small steepness ($s \ll 1$). Consider a triad of waves ($\vec{k}_0, \Omega_0$), ($\vec{k}_1, \Omega_1$) and ($\vec{k}_2, \Omega_2$) with $\Omega_i = F(\vec{k}_i)$, $0 \le i \le 2$. Since the waves are nonlinearly interacting,

$$ \vec{k}_0 + \vec{k}_1 + \vec{k}_2 = 0. \quad (4.4) $$

This relation is a spatial resonance relation. If, moreover,

$$ \Omega_0 + \Omega_1 + \Omega_2 = 0, \quad (4.5) $$

with $\Omega_i$ of any sign, efficient interactions can occur among the wave triad. This relation is a temporal resonance relation. If the temporal resonance relation is not satisfied, the triad interaction produces harmonic perturbation, at frequency $\Omega_0 + \Omega_1 + \Omega_2$, without cumulative effects on wave amplitudes. These selective interactions are named 'resonant interactions'. They are selective because they imply both a spatial and a temporal resonance relation. It was shown by Hasselmann (1967) that "sum interactions" $|\Omega_1| + |\Omega_2| = |\Omega_0|$ are unstable. In the limit of small-scale secondary waves, that is $|\vec{k}_1|, |\vec{k}_2| \gg |\vec{k}_0|$, the instability is of the parametric type: $\Omega_1 \simeq \Omega_2 \simeq -\Omega_0/2$. Hasselmann also showed that "difference interactions" $|\Omega_1| - |\Omega_2| = |\Omega_0|$ are stable.

The energy and the quasi-momentum of the waves ($\vec{P}_i = \frac{E_i}{\Omega_i} \vec{k}_i$ for wave *i*) are conserved within the triad: $E_0 + E_1 + E_2 = \text{constant}$ and $\vec{P}_0 + \vec{P}_1 + \vec{P}_2 = \text{constant}$. But the
wave action, defined for each wave by $\frac{E_i}{\Omega_i}$, is not conserved: $\frac{E_0}{\Omega_0} + \frac{E_1}{\Omega_1} + \frac{E_2}{\Omega_2} \neq$ constant (Muller et al., 1986).

**Resonant interaction equations.** Resonant interactions among a triad occur over a slow time scale, inversely proportional to the wave steepness $s$ (we recall that the theory is valid for $s \ll 1$). A slow time scale $t_1$ is therefore introduced, in addition to the (rapid) time scale $t$ ($\simeq N^{-1}$) over which the waves evolve: $t_1 = \frac{t}{s}$. To describe the interaction, the amplitude and the phase of each wave is assumed to vary over the slow time scale $t_1$. A two-time scale asymptotic expansion of the velocity and density fields is performed and, at first order in $s$, the evolution equations of the amplitudes $A_i$ are obtained (f.i. McEwan and Plumb, 1977):

$$ \dot{A}_0 + T_0 A_0 = S \Omega_0 \overline{A}_1 \overline{A}_2, \quad (4.6) $$

$$ \dot{A}_1 + T_1 A_1 = S \Omega_1 \overline{A}_0 \overline{A}_2, \quad (4.7) $$

$$ \dot{A}_2 + T_2 A_2 = S \Omega_2 \overline{A}_0 \overline{A}_1, \quad (4.8) $$

where $T_i$ is the viscous damping rate, $S$ is the interaction coefficient, which depends upon the wave vectors of the triad, and the bar denotes the complex conjugate.

## 4.2 Stability of a primary wave of infinitely small steepness

**Theoretical results.** The amplitude equations (4.6)-(4.8) can be used to investigate the stability of a monochromatic internal gravity wave, termed the primary wave. Let $(\vec{k}_0, \Omega_0)$ be the primary wave and assume that $A_1, A_2 \ll A_0$. Equations (4.6)-(4.8) may then be linearized about the basic state of the primary wave. An evolution equation for $A_1$ and $A_2$ is obtained (neglecting viscous effects for simplicity):

$$ \ddot{A}_1 = S^2 \Omega_1 \Omega_2 |A_0|^2 A_1. \quad (4.9) $$

It follows immediately that if $\Omega_1\Omega_2 > 0$, exponential growth of the secondary waves is possible, with growth rate equal to $S|A_0|\sqrt{\Omega_1\Omega_2}$. This condition was shown by Hasselmann (1967) to be equivalent to sum interactions, that is, the triad is unstable. The growth rate of the secondary waves is plotted as a function of $\vec{k}_1$ in Figure 12 for the run we consider below, which we shall refer to as run R1 for simplicity. For very large wavenumbers, resonance occurs through parametric instability and the growth rate saturates: there is no scale selection in the inviscid case, arbitrary small scales being excited by the instability. When viscous effects are taken into account, $A_0$ must be larger than a threshold amplitude for amplification to occur. The smallest scales are damped by viscosity, whatever its (finite) value, so that the viscosity provides the scale selection. In the case displayed in Figure 12, a wavenumber close to 5 is selected.

**Numerical results.** The theoretical results are illustrated by studying numerically the stability of a primary wave of large scale relative to the numerical domain (Koudella and Staquet, 2005). The Boussinesq equations are solved in a vertical plane for this purpose.
emerge as time elapses. The perturbation vorticity appears to be organized as bands of alternate sign, which make an angle $\theta_1$ such that $N\cos\theta_1 \approx \omega_0/2$. Also, the sign of a given band changes sign after one primary wave period and recovers its sign after two periods. These observations are consistent with the perturbation being parametrically excited by the primary wave.

The fact that the bands are associated with parametrically excited secondary waves is confirmed in Figure 14b. This calculation differs from the one displayed in frame (a) by the value of the viscosity, which has been decreased by a factor 100. The same bands are visible, with the same inclination (because the excited waves have a frequency close to $\Omega_0/2$), but the wavelength is much smaller because of the smaller scale selection by the lower viscosity.

**Energetics of the parametric instability.** In this section, we analyse the mechanism through which the perturbation grows. As usual in this type of study, we consider the energetics of the perturbation. The novel point here is that the spatio-temporal dependency of the forcing mechanisms of the perturbation is described (Koudella and Staquet, 2005). We start with a rotation of the reference frame (f.i. Mied, 1976) so that $Z$ is aligned with the primary wave vector $\vec{k}_0$ in the new $(0, X, Z)$ frame. We decompose the velocity, density and pressure fields into a part associated with the primary wave and a part associated with the perturbation. The Boussinesq equations linearized about the primary wave state are, for the kinetic and potential energy of the perturbation (Lombard

<figure><img src="image_13.png" /><figcaption>**Figure 13.** (a) Vorticity $\zeta_0 = dU_0/dZ$ and density fields of a monochromatic wave. Note the quadrature relationship between the two fields. A few isopycnals are shown. The dark (light) inclined bands are the regions where $dU_0/dZ$ is negative (positive). (b) Same as (a), with regions where either the kinetic energy or the potential energy of the perturbation is forced. RSS (ISS): reduced (increased) static stability region (from Koudella and Staquet, 2005).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 14.** (a) Constant contours of the perturbation vorticity for run 1. (b) Same as (a) except that the viscosity has been decreased by a factor 100 (from Koudella and Staquet, 2005).</figcaption></figure>

and Riley, 1996):

$$ \frac{D}{Dt} \left[ \frac{1}{2} (U'^2 + W'^2) \right] = -\nabla(\vec{U}' \cdot P') - BF - U'W' \frac{dU_0}{dZ} \quad (4.10) $$

$$ \frac{D}{Dt} \left[ \frac{1}{2} \frac{\rho'^2}{N^2} \right] = BF - \frac{1}{N^2} \rho' W' \frac{dR_0}{dZ}. \quad (4.11) $$

where $(U', W', P', \rho')$ are the perturbation fields in the rotated $(0, X, Z)$ frame and BF is the buoyancy flux. In this frame, the primary wave velocity is $(U_0, 0)$; we denote its density field by $R_0$.

These equations immediately show that the kinetic energy of the perturbation is forced by the term: $-U'W' \frac{dU_0}{dZ}$, which involves the shear of the primary wave (as for any unstable shear flow) and the potential energy is forced by the term: $\frac{1}{N^2} \rho' W' \frac{dR_0}{dZ}$, which involves the vertical gradient of the primary wave density field.

Let us consider the forcing of the kinetic energy. If we assume that the perturbation is made of a superposition of linear internal gravity waves, the correlation $U'W'$ is one-signed. It is not difficult to show here that $U'W' \ge 0$ in the rotated reference frame: the velocity field of the perturbation is aligned with the vorticity bands displayed in Figure 14 so that $U'$ has the same sign as $W'$. On the other hand, $dU_0/dZ$ oscillates in time and space. At a given position, this implies that the forcing term $-U'W' \frac{dU_0}{dZ}$ is positive, i.e. the perturbation is forced, when $dU_0/dZ$ is negative; but the forcing would be negative when $dU_0/dZ$ becomes positive at that position, implying that the perturbation would give back the acquired energy. How is a non zero net forcing possible? The point is that the perturbation will not give its kinetic energy back if the correlation term $U'W'$ is zero, or very small, during the half primary wave period where $dU_0/dZ$ is positive. This
will be the case if the perturbation energy is of potential form during this half-period. It follows that optimum energy forcing occurs when the perturbation energy is of kinetic form during the half-period when $dU_0/dZ < 0$ and of potential form during the other half-period when $dU_0/dZ > 0$.

Since $dU_0/dZ$ is of the form $\exp(I\Phi_0)$, the energy of the perturbation should be of this form too. It follows that the velocity and density fields of the perturbation should behave as $(U', W', \sigma') \sim \exp(I\frac{\Phi_0}{2})$. In other words, the instability must be of the parametric type for optimal forcing to occur.

What about the potential energy of the perturbation? It can be shown that its forcing term $\frac{1}{N^2}\rho'W'\frac{dR_0}{dZ}$ is always positive. Since $dR_0/dZ$ oscillates at the primary wave period, this implies that the perturbation potential energy is forced twice per primary wave period: at the locations where $dR_0/dZ$ is maximum, where the vertical density gradient experienced by the perturbation is decreased, and at locations where $dR_0/dZ$ is minimum, where the vertical density gradient experienced by the perturbation is increased. The former region is referred to as RSS (reduced static stability regions) and the latter as ISS (increased static stability regions). By contrast, as shown above, the kinetic energy is forced once, in regions where the primary wave shear $dU_0/dZ$ is negative, and this is where vorticity bands are visible in Figure 14. Due to the continuous conversion of kinetic energy into potential energy, it can be shown that potential energy is converted into kinetic energy in a ISS region, the reverse occurring in a RSS region. As a consequence, the perturbation potential energy always increases in RSS regions.

A summary of the spatial occurrence of these forcing mechanisms is displayed in Figure 13b.

As mentioned above, the perturbation growth is controlled by a few parametrically resonant triads, which may be assumed to be independent as they grow. Hence, the energetics of the perturbation may be reduced to that of a single resonant triad and approximate expressions of the perturbation energy and of energy transfer terms between the primary wave and the perturbation may be computed from this assumption. It can be shown in particular that the phase average over $2\pi$ of the kinetic energy transfer term $(-U'W'\frac{dU_0}{dZ})$ to the potential energy transfer term $(\frac{1}{N^2}\rho'W'\frac{dR_0}{dZ})$ is equal to $\cos(\theta_1-\theta_0)$, that is, $\approx 0.91$ for run R1. Hence, potential energy is transferred at a higher rate than kinetic energy in an unstable primary wave of small amplitude. This is actually a very general result in stratified turbulence, where it is found that potential energy transfer toward small scales is higher than kinetic energy transfers. The analogy can be pursued in the present case since parametric instability excites small scale waves. Two arguments can be provided to account for this behaviour : the incompressibility constraint, which applies to the velocity but not to the density field (Lesieur, 1997) and the subset of two-dimensional interactions which prevent kinetic energy from being transferred toward small scales (Holloway and Ramsden, 1988).

What does happen as the perturbation amplifies? The amplifying perturbation vorticity field in $dU_0/dZ < 0$ regions may become unstable to a (large scale) Kelvin-Helmholtz instability. On the other hand, the perturbation potential energy continuously increases in RSS regions, eventually leading to unstable fluid layers; hence, a (small-scale)
buoyancy-induced instability may also occur there. However, the unstable fluid layers are stabilized by the shear of the primary wave so that, in a two-dimensional computation, a Kelvin-Helmhotz instability eventually occurs. In a three-dimensional medium by contrast, the buoyancy-induced instability can develop in the plane perpendicular to the primary wave propagation plane because there is no shear there. Hence, in a three-dimensional fluid, both instabilities compete. Koudella and Staquet (2005) showed that the buoyancy-induced instability grows faster than the Kelvin-Helmholtz instability, unless viscous effects are strong enough, thereby confirming earlier heuristic analysis (Munk, 1981). It follows that an unstable primary wave of small steepness eventually breaks down through a three-dimensional buoyancy-induced instability (if viscous effects are low enough).

## 4.3 Stability of a primary wave of finite steepness

**Theoretical results.** When the primary wave steepness is smaller than 1 but not infinitesimal, the use of a linear stability analysis against two-dimensional perturbations allows one to go beyond the resonant interaction theory. The early analyses neglected Coriolis effects. Drazin (1977) and Mied (1976) independently considered the case of a perfect fluid while Klostermeyer (1983) extended their results to the case of a real fluid. The fundamental conclusion of this two-dimensional analysis is that a primary wave is always unstable, whatever the associated Richardson number (that is whatever the intensity of stratification $N$ relative to the wave-induced shear). The instability is again of the parametric type in the case of small scale perturbation. In the limit of a vanishing primary wave steepness, Drazin (1977) showed that the results of the linear stability analysis coincide with those of the resonant interaction theory.

By contrast, when the primary wave amplitude is so strong that isopycnals are overturned ($s > 1$), the analysis must be conducted in a three-dimensional framework. A linear stability analysis of this basic state to two- and three-dimensional perturbations was conducted by Klostermeyer (1991) ; only one orientation of the primary wave vector and of the plane of the three-dimensional perturbation was considered (this plane was perpendicular to the primary wave propagation plane ($\vec{k}, \vec{g}$)). Klostermeyer (1991) found that the most amplified perturbation is three-dimensional and grows at small scale through a parametric instability mechanism. This pioneering work was complemented by Lombard and Riley (1996) and Sonmor and Klaassen (1996). Largest growth rates were found for perturbations at an angle oblique to the primary wave propagation plane but, as Lombard and Riley (1996) first realized, the features of the three-dimensional instability (wave number and amplitude) strongly depend upon the primary wave propagation angle (see their Figure 6).

**Numerical results.** These results are illustrated by three computations, referred to as runs R2, R3 and R4 (Koudella, 1999). The first one (run R2) is analogous to run R1 above, except that the initial primary wave now evolves in three dimensions and is perturbed by a three-dimensional noise of small amplitude (Figure 15a). Such a small steepness wave is unstable to parametric instability in its propagation plane implying that the flow remains two-dimensional during the instability growth (Figure 15b). The
secondary resonant waves amplify most in RSS regions, as we showed it above. This leads to overturned isopycnals ($s$ is locally greater than 1), which are unstable. A buoyancy-induced instability develops in the plane perpendicular to the propagation plane of the primary wave and the flow becomes three-dimensional (Figure 15c). Like the parametric instability, the buoyancy-induced instability amplifies arbitrary small scale waves so that the viscosity provides the scale selection mechanism. In Figure 15c, a wavelength equal to $\simeq 1/7$ of the domain size along the $y$-direction is visible. Because the instability has already reached a nonlinear stage in this frame, pairs of counter-rotating vortices (one pair per wavelength) are actually visible in the Figure. These small energetic scales are eventually dissipated by viscosity and the flow relaminarizes, because the primary wave is not forced.

It is remarkable that the validity of the resonant interaction theory extends much beyond its theoretical domain of validity, as we shall show it. We present in Figure 16 a run analogous to run R2 except that the initial primary wave steepness is twice larger ($s = 0.724$, run R3). An in run R2, the isopycnals locally amplify up to overturning while remaining two-dimensional: the primary wave becomes unstable through parametric instability during this stage (first two frames of the Figure). The overturned isopycnals become unstable through a buoyancy-induced instability (third frame). Because the viscosity is twice higher in this run than in run R2, a larger wavelength is selected along the $y$-direction, equal to 1/3 of the domain size. The whole flow eventually becomes turbulent and dissipates.

Details of the formation of the buoyancy-induced instability in this run are presented in Figure 17. Snapshots of the density field are displayed in a plane perpendicular to the propagation plane of the primary wave. Three wavelengths are thus visible in the second frame (instead of $\simeq 7$ for run R2). The nonlinear development of the instability gives rise to coherent structures with a mushroom form (fifth frame), which yield a pair of vortices of alternate sign. A zoom of the formation of the mushroom structures is displayed on the third row.

The formation of mushroom structures as a result of a buoyancy-induced instability has also been observed in laboratory experiments. Figure 18a stems from the experiment of Benielli and Sommeria (1998), already referred to in Section 3.3: it displays a visualization in a plane perpendicular to the propagation plane of a large scale internal gravity wave excited through parametric instability; two mushroom structures are visible. Figure 18b stems from the laboratory experiment of a stably-stratified shear layer performed by Schowalter et al. (1994).

The above results show that overturned isopycnals are spontaneously unstable to a three-dimensional buoyancy-induced instability. This is illustrated in Figure 19 through the temporal evolution of three isopycnal surfaces for an internal gravity wave of initial steepness larger than 1 (run R4). In this case, the viscosity is 50% larger than in run R3 and two wavelengths only appear only the $y$-direction. As previously noted, the whole wave field eventually breaks down through this instability. Such large amplitude internal gravity waves are encountered in the atmosphere, due to the rarefaction of air (lowering of density) as the wave propagates upward. Other mechanisms leading to a local amplification of the wave are discussed in the next section.
<figure><img src="image_15.png" /><figcaption>**Figure 15.** Breaking of an internal gravity wave with $s = 0.36$ (resolution $256^3$). Six constant density surfaces are plotted at successive times: (a): $t = 16.1$ BVP, (b): $t = 38.2$ BVP, (c): $t = 42$ BVP, (d): $t = 44$ BVP. In frames (a) and (b), a two-dimensional secondary instability develops, leading to local overturnings of the isopycnals; in frame (c), a three-dimensional instability develops upon overturned isopycnals, and in (d), local breaking occurs. BVP is the Brunt-Väisälä period $2\pi/N$ (from Koudella, 1999).</figcaption></figure>

**Influence of Coriolis effects on the instability of a primary wave.** When Coriolis effects are taken into account, the linear stability analysis of Dunkerton (1997) (see also references therein) and numerical experiments of Lelong and Dunkerton (1998) show that the instability mechanism depends upon the intrinsic frequency of the primary
wave. When $\Omega$ is close to $f$, the velocity field induced by the wave is nearly steady (relative to the perturbation growth rate) and nearly horizontal (with a sinusoidal vertical dependency) so that this velocity field is close to a parallel steady shear flow: it may thus bear a Kelvin-Helmholtz instability. Such an instability is found to occur whatever the primary wave amplitude. By contrast, when $\Omega$ is larger than $f$, the behaviour depends upon the primary wave amplitude. Thus, resonant interactions occur for a vanishingly small amplitude (Miyazaki and Adachi, 1998). When $s$ has a finite value, even larger than 1, a shear instability prevails as long as $\Omega$ is not too large. When $s > 1$, this instability is more and more modified by convection as $\Omega$ increases and, in the limit of no rotation effect ($\Omega \gg f$), it becomes a three-dimensional convective instability.

<figure><img src="image_16.png" /><figcaption>**Figure 16.** Breaking of an internal gravity wave with initial steepness twice larger than in Figure 15 (resolution 256³). Three isopycnals are displayed at successive times. First row, from left to right, in unit of the Brunt-Väisälä period (BVP) $2\pi/N$: $t = 18.9, 20, 21.25$. Second row, from left to right: $t = 21.7, 22.4, 23.6$ (from Koudella, 1999).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 17.** Vertical cross-sections of the density field for the primary wave displayed in Figure 16, at successive times. The last two frames are a zoom of the buoyancy induced instability. Visualizations are made in the plane perpendicular to the primary wave propagation plane (from Koudella, 1999).</figcaption></figure>

# 5 Interaction of internal gravity waves with an ambient shear flow

## 5.1 Geophysical context

In geophysical flows, internal gravity waves generally interact with an ambient shear flow: for instance, the shear flow is the wind in the atmosphere (Bretherton, 1966) or a current in the ocean (Badulin et al., 1990). The former case is very common : the blowing of a wind over a mountain range generates lee waves which next interact with the wind. Recent research in the Sun has revealed that the differential rotation that exists in this star creates a very strong shear at the basis of the convective zone, which must interact with the waves that are produced there (Talon et al., 2002). How does this shear flow interact with the waves is unknown. An analogous interaction occurs when the waves encounter a temperature (or a density) front. Actually, when Coriolis effects
<figure><img src="image_18.png" /></figure>

**Figure 18.** Laboratory experiments in which a buoyancy-induced instability develops (see text for details). Left frame: from Benielli and Sommeria (1998); right frame: from Schowalter et al. (1994).

come into play, such a front creates a vertical shear via the thermal wind balance so that the wave interact with both a density front and a shear. And when there is no ambient (large scale) shear or density front, as in the deep ocean, the wave field interacts with the shear flow it induces.

Why is the wave-vortex interaction an important geophysical problem? A very simple example is provided by lee waves. Since lee waves are produced by the wind, energy is extracted from the wind to generate these waves. As a result, the wind speed decreases over a mountain range because of this generation process, as if the mountain were exerting a drag force on the wind. This drag force has been introduced (as a subgrid-scale model) for the alpine range in the ECMWF<fnref n="3" /> weather forecast model, thereby improving the weather predictions for this region (Lott, 1995). Another example may be provided by the permeability of the polar vortex. The polar vortex forms over the pole in the winter hemisphere and is a quasi-two-dimensional stratospheric vortex. It extends from the bottom of the stratosphere (≈ 8 kms at the pole) up to an altitude of 70 kms or so, while its horizontal extent is of a few thousand kms. Such a quasi-two-dimensional vortex is known to be a barrier to transport: laboratory experiments by Heijst et al. (1991) for instance show that dye injected in an isolated quasi-two-dimensional vortex remains capped inside the vortex and do not mix with the fluid outside the vortex. Atmospheric measurements at mid-latitudes still led to the suspiscion that air possibly leaks outside the polar vortex and a mechanism was seeked for to explain this leakage. McIntyre (1995) suggested that the interaction of inertia-gravity waves (which are internal gravity waves subjected to rotation) with the vortex edge could result in an irreversible transport of mass across the vortex edge. We finally report about a subtle mechanism of wave-vortex interactions in the ocean. Klein et al. (2003) showed that meso-scale eddies (which are

<footnote n="3">³European Center for Meteorological Weather Forecast.</footnote>
<figure><img src="image_19.png" /></figure>

**Figure 19.** Breaking of an internal gravity wave with overturned isopycnals at $t = 0$ ($s = 1.086$, resolution $256^3$). Three constant density surfaces are plotted at successive times. First row, from left to right, in unit of the Brunt-Väisälä period (BVP): $t = 9.7$, 10.6, 11.5. Second row, from left to right: $t = 11.8$, 12.3, 14.1 (from Koudella, 1999).

responsible for the transport of heat toward the poles), through their effects on near-inertial ($\Omega \approx f$) oscillations, can induce a spatially heterogeneous diffusion that depends on the vorticity sign. In turn, the heterogeneous diffusion breaks the cyclone-anticyclone symmetry on short time scales (compared to the diffusive time scales) and favors the emergence of cyclonic structures.

## 5.2 Ray theory

In the following, we focus upon the interaction in an infinite medium of internal gravity waves with an ambient large scale horizontal shear flow $\vec{U}$, and, when rotation is present, with a density front as well.

We consider a monochromatic wave of intrinsic frequency $\Omega$ and wave vector $\vec{k}$ (or a wave packet with main intrinsic frequency and wave vector $\Omega$ and $\vec{k}$). The intrinsic frequency is the frequency measured in a frame of reference attached to the ambient flow. The simplest, and main, effect of the interaction on the internal wave is the change of
the intrinsic wave frequency through the Doppler effect:

$$ \omega_0 = \Omega + \vec{k} \cdot \vec{U}, \qquad (5.1) $$

$\omega_0$ being the frequency of the source that emits the waves. This source is supposed to be fixed in the frame of reference relative to which $\vec{U}$ is measured. $\Omega$ changes as the wave propagates in the changing velocity field, and may approach its lower or upper bound. In this case, further propagation is no longer possible and, in the linear limit, the wave may either be trapped or reflected. Two academic situations are usually considered. Within an atmospheric context, the mean flow is a horizontal wind with a vertical shear $U(z)\vec{i}_x$; within an oceanic context, the mean flow is a horizontal current with a horizontal shear $U(y)\vec{i}_y$.

The change in $\Omega$ as the wave packet propagates into the current is most easily predicted within the WKB approximation (see Olbers, 1981, for a very clear presentation of this approximation). This approximation relies upon the assumption that the properties of the fluid medium that affect the wave propagation ($\vec{U}$ and $N$, in the present case) vary slowly in time and space relative to the wave intrinsic frequency and wavelength respectively. Hence, the medium may be assumed to be uniform and steady over a length scale of order $|\vec{k}|^{-1}$ and over a time of order $\Omega^{-1}$. Under this assumption, the evolution of the wave vector is known along a ray (defined as $d\vec{x}/dt = \vec{c}_g + \vec{U}$) and is driven by the gradients of the ambient velocity and buoyancy fields. In this paper, we shall consider that the fluid medium is steady (i.e. $\partial U/\partial t = 0$, $\partial N/\partial t = 0$). In this case, the absolute frequency is constant along a ray: $d\omega_0/dt = 0$, where $d/dt = \partial/\partial t + (\vec{c}_g + \vec{U}) \cdot \nabla$ denotes the material derivative following a ray. The equations governing the refraction of the wave vector along a ray, known as the ray equations, are:

$$ \frac{dk_i}{dt} = -\frac{\partial \Omega}{\partial N} \frac{\partial N}{\partial x_i} - k_j \frac{\partial U_j}{\partial x_i}, \qquad (5.2) $$

that is, in the present case, with $U = (U(y, z), 0, 0)$ and $N(y, z)$

$$ \frac{dk_x}{dt} = 0 \qquad (5.3) $$

$$ \frac{dk_y}{dt} = -\frac{\partial \Omega}{\partial N} \frac{\partial N}{\partial y} - k_x \frac{\partial U}{\partial y}, \qquad (5.4) $$

$$ \frac{dk_z}{dt} = -\frac{\partial \Omega}{\partial N} \frac{\partial N}{\partial z} - k_x \frac{\partial U}{\partial z}. \qquad (5.5) $$

As for the relative frequency, it changes along a ray according to the equation

$$ \frac{d\Omega}{dt} = -c_{gi}k_x \frac{dU}{dx_i}, \qquad (5.6) $$

where $c_{gi}$ is the component of the group velocity along the direction of inhomogeneity $x_i$.

Changes of the wave amplitude are inferred from the conservation of wave action. For any slowly varying background, the action $A = E/\Omega$, where $E$ is the wave energy,
<figure><img src="image_20.png" /></figure>

**Figure 20.** Sketch of a wave packet approaching a critical level $z = z_c$, within the WKB approximation. Two rays are shown in the vertical $(x, z)$ plane, with $\vec{c}_g + U(z)\vec{i}_x$ being the absolute group velocity. The intrinsic frequency of the wave packet, $\omega_0 - U(z_c)k_x$, from the Doppler relation, vanishes at the critical level.

satisfies the conservation equation (Bretherton, 1968):

$$ \frac{\partial A}{\partial t} + \nabla \cdot [(\vec{c}_g + \vec{U})A] = 0. \qquad (5.7) $$

Equation (5.7) implies that the action contained in a small volume $\delta V$ moving with the absolute group velocity is conserved, that is

$$ \frac{d(A\delta V)}{dt} = 0. \qquad (5.8) $$

The form of the WKB theory we use is the approximation of geometrical optics but, for simplicity, the terminology *WKB approximation* will be employed hereafter.

## 5.3 Horizontal mean flow with a vertical shear $U(z)\vec{i}_x$

This situation was first investigated without rotation effect, when a wave packet propagates upwards in a vertical shear flow such that its intrinsic frequency decreases. If there exists an altitude at which $\Omega$ vanishes, the wave cannot propagate beyond this altitude, known as the critical level. Hence the critical level acts as a wave filter. This level is also defined as the altitude at which the component of the phase velocity of the wave along the wind direction equals the wind velocity (Figure 20).

The first theoretical approach to this problem was performed by Bretherton (1966), who considered the Boussinesq equations in a non-rotating frame, linearized about a basic state defined by the wind and a linear stratification, under the WKB approximation. Note that this approximation implies that the Richardson number ($Ri = N^2/(dU/dz)^2$) is much larger than 1 everywhere (Olbers, 1981). For a wave packet approaching the critical level, Bretherton (1966) showed that the vertical component of the group velocity $c_{gz}$ decreases as $\eta^2$, where $\eta$ is the distance of the wave packet to the critical level, while $k_z$ increases as $\eta^{-1}$. From the dispersion relation (3.5), $\Omega \to 0$ (since $k_x$ and $k_y$ are constant). Even if these results are correct as the critical level is approached, the WKB approximation diverges in the immediate vicinity of that level. The theory predicts
that this level is reached in an infinite time (so that there should be no transmitted component) and that the wave-induced energy increases without bound, as $\eta^{-1}$.

Relaxing the WKB approximation, but still considering the linearized equations of motions and assuming that the Richardson number is everywhere greater than 0.25, Booker and Bretherton (1967) showed that, for a monochromatic wave, momentum is actually transferred from the wave to the mean flow at the critical level, except for a weak transmitted component whose energy is the incident wave energy reduced by the factor $\exp[-2\pi(Ri_c - 1/4)^{1/2}]$, where $Ri_c$ is the value of the Richardson number at the critical level. The transfer of momentum to the mean flow is manifested as a discontinuity in the vertical flux of wave-induced horizontal momentum $\overline{\rho u' w'}(z)$ at the critical level (where the bar denotes an average over a horizontal wavelength); this flux is otherwise constant with altitude in the absence of critical level (Eliassen and Palm, 1961). No component is reflected in this linear limit.

For a high enough wave amplitude, the incoming wave-induced energy is higher than that absorbed by the shear flow and, consequently, the energy density increases in the neighbourhood of the trapping level. Breaking eventually occurs, thereby partly dissipating the accumulated energy (since momentum is conserved however, momentum will be locally deposited by the breaking waves). This process is well-known to occur in the atmosphere. However, as noted by McIntyre (2000), there is a subtle effect which implies that, unless viscous effects are too strong, the linear theory always fails after some time: as the wave approaches the critical level, its intrinsic phase velocity decreases faster than the wave-induced velocity along the wind direction does because of absorption so that their ratio, which is one measure of the wave steepness, increases. Hence, breaking may occur as well.

This situation was studied in detail through several experimental works (f.i. Koop, 1981; Koop and McGee, 1986) and is exemplified in Figure 21a for a weak amplitude wave and in Figure 21b for a large amplitude wave. Koop and McGee (1986) found that when breaking occurs, it does so through a two-dimensional Kelvin-Helmholtz instability when no critical level exists and through a three-dimensional buoyancy-induced instability when there is a critical level. The former behaviour can be explained as follows: as the waves propagate in the shear flow, the intrinsic frequency decreases so that wave-induced motions are close to horizontal, slow and of small steepness; such motions are therefore prone to Kelvin-Helmholtz instability. The latter behaviour is consistent with the general analysis of Lelong and Dunkerton (1998), who show that a gravity wave of low frequency and high steepness breaks through a buoyancy-induced instability. In numerical experiments, the breaking problem was first addressed in a two-dimensional vertical plane. The idea was that, as $k_z \to 0$, a strong wave-induced vertical shear develops (which actually dominates the mean flow shear), so that destabilisation should occur via a Kelvin-Helmholtz instability. The two-dimensional simulations realized by Winters and D'Asaro (1989) actually revealed an unexpected behaviour once the strong wave shear has built up, consisting of a long regime of overturned isopycnals without breaking. This led Winters and Riley (1992) to design a model for a gravity wave packet in the neighbourhood of the critical level, as a basic state for a linear stability analysis. They found that the most unstable modes are both of the Kelvin-Helmholtz type (occurring in the propagation plane of the wave packet) and of the convective type (occurring in a
<figure><img src="image_21.png" /><figcaption>**Figure 21.** Laboratory experiment of an internal gravity wave interacting with a vertical shear flow. Both the shear flow and the wave field are created by a moving topography. (a) The wave amplitude is small enough for complete absorption of the wave field by the shear flow to occur; (b) for a larger wave amplitude, the wave breaks below the critical level (from Koop and McGee, 1986).</figcaption></figure>

plane perpendicular to it). Three-dimensional numerical simulations (Fritts et al., 1994; Winters and D'Asaro, 1994) later confirmed that the dominant instability is a three-dimensional convective instability, which the two-dimensional simulations were unable to reproduce. Only Winters and D'Asaro (1994) studied the influence of the breaking of a wave packet upon the fluid medium, for $Ri = 25$ and for three values of the wave packet initial amplitude. The main result is as follows: the energy carried by the wave packet is first (in time) transferred to the mean flow, the more so when the initial wave amplitude is smaller, in agreement with linear theory. A reflected component is next produced, whose energy content grows with the initial wave amplitude (no component seems to be transmitted, whatever the amplitude). Consequently, as noted by Winters and D'Asaro (1994), very little energy is left for fluid mixing and for potential vorticity changes.

The work of Booker and Bretherton (1967) was extended to rotating flows by Jones (1967) and by Wurtele et al. (1996). Using a linearized approach, Jones (1967) showed that, in addition to the classical critical level (characterized by $\Omega = 0$), two additional singular levels exist, corresponding to $\Omega = \pm f$. His paper contains another important result: in a rotating fluid, the vertical flux of wave-induced horizontal momentum is no longer conserved (away from the critical level). It should be replaced by the vertical flux of wave-induced angular momentum for this conservation property to hold again (away from the singular levels). Wurtele et al. (1996) further showed that the wave becomes evanescent in the neighbourhood of the critical level (since the wave can only propagate for $\Omega^2 > f^2$) so that the only effective singular levels are those where $\Omega^2 = f^2$. Wurtele et al. (1996) also investigated numerically the situation of a time developing wave, being emitted from a source at $t = 0$. In this case, the singularity develops in time as well as at early times, the propagating wave crosses the $\Omega = f$ (for instance) level and decays
in the evanescent region. As time elapses and the wave reaches a steady state (i.e. its amplitude becomes constant), the singularity develops but non linear effects develop as well. As a result, the wave breaks in the neighbourhood of the singular level and a reflected component is emitted; the wave would be absorbed at the singular level in a linear regime. These results are thus analogous to those found by Booker and Bretherton (1967) at a critical level in a non rotating fluid. Wurtele et al. (1996) also showed that the behaviour is notably different for a continuous spectrum of frequencies (such as a lee wave generated by a flow over an arbitrary topography), in which case no singular behaviour is encountered.

In the discussion above, the monochromatic wave propagates in the shear flow so that its intrinsic frequency decreases. When the wave propagates in a shear flow so that its intrinsic frequency increases, it reflects onto the horizontal plane where $\Omega = N$, in the linear limit. Jones (1968) showed that overreflection occurs when the local Richardson number is smaller than 0.25 somewhere, that is, the shear flow is potentially unstable; the wave extracts energy from that flow when reflecting. When nonlinear effects come into play (due to the high initial amplitude of the wave for instance) and the background shear flow is stable, Sutherland (2001) showed that a horizontally periodic wave packet permanently deposits momentum to the mean flow at altitudes close to and below the reflecting level; when the wave packet is horizontally compact, a substantial part of the finite amplitude wave packet energy may be transmitted.

## 5.4 Horizontal mean flow with a horizontal shear $U(y)\vec{i}_x$

**Theoretical results.** The behaviour of an internal gravity wave packet in a horizontally sheared current $U(y)\vec{i}_x$ was studied thirty years ago within an oceanographic context, using the linear theory in the WKB approximation (Ivanov and Morozov, 1974; Olbers, 1981; Basovich and Tsimring, 1984; Badulin et al., 1985; Badulin and Shrira, 1993). Coriolis effects are ignored in all studies except in Olbers' work. When the wave packet enters into the current and propagates against it so that its intrinsic frequency increases (Eq. (5.6) with $x_i = y$), it cannot propagate beyond the position $y_t$ at which $\Omega = N$. The mean flow being barotropic, this position actually is a vertical plane if $N$ is constant; this plane is hereafter referred to as the trapping plane (Figure 22).

Since the properties of the medium in which the wave propagates vary only with $y$, $k_x$ and $k_z$ remain unchanged (as well as $\omega_0$ as already noted, because the medium is steady). From the linear dispersion relation (3.5), one easily infers that $k_y$ goes to infinity as $\Omega \to N$. More precisely, it can be shown from WKB theory that $k_y \sim \eta^{-1/2}$ as $\eta \to 0$, where $\eta$ refers again to the distance of the wave packet to the trapping plane; moreover, $c_{gy} \sim \eta^{3/2}$ and $c_{gx}$, $c_{gz} \to 0$ as well (Staquet and Huerre, 2002). The latter property implies that the wave packet slows down in the neighbourhood of the trapping plane so that its energy density locally increases (in other words, the wave-induced energy accumulates in the neighbourhood of this plane). However, the WKB theory also predicts that this energy tends to infinity ($E \sim \eta^{-3/2}$) and that the trapping plane is reached in an infinite time. As in the critical level situation, the two latter results are unphysical and stem from the failure of the WKB theory in the immediate vicinity of the trapping plane. Olbers (1981) actually noted that, if $E \sim \eta^{-\mu}$, the asymptotic behaviour of the wave is
<figure><img src="image_22.png" /><figcaption>**Figure 22.** Sketch of a wave packet approaching a trapping plane, within the WKB approximation, and propagating against the current ($\vec{k} \cdot \vec{U} = k_x U < 0$). Two rays are shown in the horizontal ($x, y$) plane, with $\vec{c}_{gh} + U(y)\vec{i}_x$ being the component of the absolute group velocity in that plane. The intrinsic frequency of the wave packet, $\omega_0 - U(y_t)k_x$, from the Doppler relation, is equal to the Brünt-Väisälä frequency at the trapping plane (adapted from Staquet and Huerre, 2002).</figcaption></figure>

regular or singular depending upon whether $\mu < 1$ or $\mu \ge 1$. In the former case, the wave reaches the plane in a finite time, and this plane is for instance a reflexion plane. In the latter case, the wave reaches the plane in an infinite time; very strong gradients of the wave-induced Reynolds stress form in the close neighbourhood of the plane, which yield momentum exchange with the mean flow. Information upon the actual wave behaviour when $\mu \ge 1$ can be obtained by solving the linearized equations of motions. Ivanov and Morozov (1974) thus found that the total wave-induced energy may indeed increase, as opposed to the critical level situation studied by Booker and Bretherton (1967). In this situation, momentum is transferred from the shear flow to the wave, so that the potential for wave breaking exists. The three-dimensional numerical study of Staquet and Huerre (2002), whereof a summary is provided below, shows that an inertia-gravity wave packet may indeed break in the neighbourhood of a trapping plane.

When the wave propagates along the current such that its intrinsic frequency decreases, $k_y$ decreases as well. If $k_y$ decays down to zero, $\Omega$ reaches a minimum value $\Omega_{min}$, obtained by setting $k_y$ to 0 in the dispersion relation. The behaviour of the wave as $\Omega \to \Omega_{min}$ may be guessed by using WKB theory. Note that the theory becomes less and less valid as the plane gets closer since the wavelength along the y-direction increases. Assuming that the theory remains valid, $k_y \sim \eta^{1/2}$, $c_{gy} \sim \eta^{1/2}$ and $E \sim \eta^{-1/2}$ as $\eta \to 0$. Here, according to Olbers' (1981) argument, $\mu = 1/2$ so that the $\Omega = \Omega_{min}$ plane is a reflexion plane for the wave.

Ollers et al. (2003) solved numerically the equation for the amplitude of a hydrostatic internal gravity wave emitted away from the shear flow and propagating toward it. The wave behaviour (transmitted, reflected, over-reflected) depends upon the stability of the shear flow. When the latter flow is inertially stable, the wave is always reflected, with a possible transmitted component. When the shear flow is inertially unstable by contrast,
over-reflection is possible, with a reflection coefficient up to 3.25 (for the case considered in the paper). One may wonder whether a high enough reflexion coefficient would not lead to the instability of the reflected wave. Note the analogy between this behaviour and that found by Jones (1968) for an unstable vertical shear flow.

**Numerical results.** The behaviour of a wave packet in a rotating, constant-$N$ medium propagating into a barotropic shear flow $U(y)\dot{\bar{i}}_x$ has been investigated numerically by Staquet and Huerre (2002). The shear flow consists of a horizontal shear layer (with a tanh profile) while the wave packet is a plane wave whose amplitude is modulated by a gaussian function along the $y$-direction. The parameters of the wave (its wave vector), of the shear flow (its maximum amplitude) and of the medium (the Brunt-Väisälä and Coriolis frequencies) are chosen so that (i) the wave intrinsic frequency increases as the wave propagates in the shear flow and (ii) a trapping plane exists.

The wave behaviour is illustrated in Figure 23 for a cyclonic shear flow (that is, the vorticity of the shear flow is of the same sign as the Coriolis frequency). Constant contours of the density field are displayed at successive times in a vertical $(y, z)$ plane. At $t=0$ (Figure 23a), the wave packet is hardly visible due to its small steepness ($s=0.26$). Because the horizontal shear flow does not displace the isopycnals, it is not visible either. Since $N$ is constant, the trapping plane is a vertical plane and its intersection with the $(y, z)$ plane is marked with a vertical line in the Figure. The wave packet exhibits two major changes as it propagates toward the shear flow: the isopycnals steepen and the wave amplitude increases. The former effect is accounted for by noting that, as $k_y$ increases, the incompressibility condition reduces to $k_yu_y + k_zu_z \approx 0$. Hence, the vector $(k_y, k_z)$ is perpendicular to the phase lines in the $(y, z)$ plane. Since $k_y \to +\infty$ while $k_z$ remains constant, the phase lines steepen. The second effect results from the trapping of the wave. The local increase of the wave amplitude makes the wave packet break (Figure 23e) and small-scale motions are produced. The latter motions are quickly dissipated however because the primary wave packet is not forced. As well, the shear flow is hardly modified by the momentum deposit that occurs during wave breaking.

The stage of the flow that follows breaking dramatically changes when the shear flow is anticyclonic (Figure 24). Indeed, in this situation, the shear flow is subjected to an inertial instability, through which small-scale motions are most amplified. The point is that the small-scale motions resulting from wave breaking act as a perturbation to the shear flow, which triggers the instability. The medium is therefore considerably modified by the breaking of the wave in this situation, because it initiates the inertial instability of the (very energetic) shear flow. The latter instability results in momentum and mass transport: the shear of the background flow is weakened and a passive scalar is transported across the trapping plane, namely across the shear flow.

## 5.5 Horizontal mean flow with both a horizontal and a vertical shear

What does happen if a wave packet in a stably-stratified rotating medium interacts with a (thermal wind) balanced shear flow $\vec{U}(y, z)$ involving both a horizontal and a vertical shear? As discussed above, the situation is not simple. For instance, when the intrinsic frequency $\Omega$ increases and approaches $N$, the wave packet should be reflected
<figure><img src="image_23.png" /><figcaption>**Figure 23.** Interaction of a wave packet (confined along the *y*-direction only) with a barotropic cyclonic shear layer. Constant contours of the density field are displayed at successive times, expressed in Brunt-Väisälä periods; (a) *t* = 0, (b) *t* = 1.9, (c) *t* = 3.8, (d) *t* = 4.6, (e) *t* = 5.6, (f) *t* = 11.4. The vertical line in frame (a) marks the trapping plane.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 24.** Interaction of a wave packet (confined along the *y*-direction only) with a barotropic anticyclonic shear layer. Constant contours of the density field are displayed at successive times, expressed in Brunt-Väisälä periods (a) *t* = 0, (b) *t* = 1.9, (c) *t* = 3.8, (d) *t* = 4.6, (e) *t* = 4.8, (f) *t* = 9.4. The vertical line in frame (a) marks the trapping plane (from Staquet and Huerre, 2002).</figcaption></figure>
by the vertical shear $\partial U/\partial z$ but trapped by the horizontal shear $\partial U/\partial y$. Also, since the shear flow satisfies the thermal wind balance, a buoyancy field $B(y, z)$ exists as well such that
$$f \frac{\partial U}{\partial z} = -\frac{\partial B}{\partial y}.$$

We have investigated this situation in a simple context, when a wave packet confined both in the $y$ and $z$ directions interacts with a horizontal shear layer with a sinusoidal vertical dependency (Edwards and Staquet, 2005). We chose the parameters such that $\Omega$ increases because of the horizontal shear of the background flow. We first explored the parameter range by solving the ray equations and performed three-dimensional direct numerical simulations (DNS) to investigate the influence of nonlinear effects on the behaviour of the wave packet.

A general behaviour is observed in the WKB theory, which is illustrated in Figure 25. The ray trajectories are displayed in frame (a) and the intrinsic frequency $\Omega$ is plotted versus time in frame (b). In frame (a), the wave packet at initial time is represented by a set of twenty rays aligned along two perpendicular segments. The $\Omega = N$ surface, plotted with a thick dashed line, displays two important locations when the interaction with the wave packet is considered: (i) where the local radius of curvature is minimum, corresponding to a maximum value of $\partial U/\partial y$ (location 1) and where the local radius of curvature is maximum, corresponding to a minimum value of $\partial U/\partial y$ (location 2). In the former case, the $\Omega = N$ surface is nearly vertical and is a trapping surface; in the latter case, the surface is nearly horizontal and is reflecting. All rays propagate toward the $\Omega = N$ surface and reach it, either in the neighbourhood of location 1 or of location 2. In the former case (location 1), the rays are trapped at the surface and propagate along it downward, with a nearly vertical group velocity, toward location 2. Note that the group velocity has strongly decreased when location 1 is reached. All rays sooner or later reach location 2 and reflect there. The rays then propagate in the interior of the shear flow within a wave guide made by the $\Omega = N$ surface. Frame (a) also displays grey points, at which the steepness of the wave packet exceeds 1. This suggests that breaking may occur there, resulting in irreversible mass and momentum transport. Note however that WKB theory is no longer valid at the trapping plane (and at a reflecting surface) and that, most importantly as we shall see, molecular effects have been ignored in this analysis.

DNS results are displayed in Figure 26 through constant contours of the fluctuating buoyancy field $b'$. It should be stressed that the steepness of the wave packet is twice smaller than in the baroclinic cases and that the varying horizontal shear ($\partial U/\partial y$) is not larger. Hence, the wave-shear interaction (which scales like $s^2\partial U/\partial y$) is smaller by a factor 4 at least. The interaction between the shear flow and the wave packet is therefore weak and one may consider that the buoyancy fluctuations displayed in Figure 26 solely belong to the wave field. The DNS behaviour is close to the WKB prediction up to the time the wave packet reaches the trapping surface at location 2 (at $t \simeq 176 \simeq 28$ BVP). This is attested in frame b) where the intrinsic frequency predicted by WKB theory and the numerical simulation are compared. Molecular effects deeply change the subsequent wave packet behaviour in the DNS, for two reasons. First, $k_y$ increases as the packet approches the trapping surface, implying that small scales along the $y$-direction are produced. These small scales are of course very sensitive to molecular effects. The
<figure><img src="image_25.png" /></figure>

**Figure 25.** WKB predictions. (a) Trajectories of rays at $t = 127.3$ BVP (Brunt-Väisälä period) in a vertical $(y, z)$ plane. The rays start from the initial wave packet location. The thick dashed line marks the intersection of the $\Omega = N$ trapping surface with the $(y, z)$ plane and dotted lines represent contours of the shear flow velocity $U(y, z)$. Light grey circles are plotted at each time the steepness along a given ray exceeds the value of 1. (b) Temporal evolution of the intrinsic frequency $\Omega$ of the wave packet normalized by the Coriolis frequency $f$. Dashed line: WKB prediction for a central ray; solid line: DNS result at the packet centre; dotted line: $N/f$ (upper bound of $\Omega/N$). Results are compared over the duration of the DNS ($t=240$, that is 38.2 BVP) (from Edwards and Staquet, 2005).

second reason, connected to the first, is that the packet slows down as it approaches the trapping plane, which makes it also prone to molecular effects. As a consequence, the nearly steady, small-scale packet is dissipated locally and does not penetrate into the wave guide.

The general conclusion that can be drawn from this study is that a single wave packet is unlikely to modify its environment. In geophysical flows however, waves are most often generated by a permanent source (like the interaction of the tide with the topography in the ocean) or, at least, are emitted during a long time with respect to their intrinsic period. Since, in our study, the wave packet slows down as it approches the trapping surface, a continuous emission of such packets would result in their superposition in the neighbourhood of the reflecting surface, and possibly in breaking. Higher resolution simulations should also be conducted to reduce the influence of molecular effects on the wave packet behaviour.

# 6 Interaction of internal gravity waves with a boundary

## 6.1 Focusing by wall reflexion

In this section, we show that internal gravity waves can locally amplify by a simple geometrical effect first discovered by Phillips (1966), which may result in breaking and
<figure><img src="image_26.png" /><figcaption>**Figure 26.** DNS results. Constant contours of the fluctuating buoyancy field $b'$ are plotted in a vertical $(y, z)$ plane at successive times (in unit of Brunt-Väisälä period): (a) $t=0$; (b) $t=8.9$, (c) $t=13.4$, (d) $t=28$. In frame (a), dotted lines represent contours of the shear flow velocity $U(y, z)$. The surface $\Omega = N$ is displayed with a solid line on all frames (from Edwards and Staquet, 2005).</figcaption></figure>

turbulence. Coriolis effects are ignored in this paragraph.

When an internal gravity wave reflects onto a sloping boundary, the frequency is conserved upon reflection, as for any reflecting wave in a steady medium. The dispersion relation (3.5) thus implies a peculiar property (Phillips, 1966): the angle that the wave vector makes *with the horizontal* is preserved as well upon reflexion<fnref n="4" />. This geometrical property implies that a ray tube (made of curves everywhere tangent to the group velocity) reflecting on a sloping boundary may be focused (see Figure 27). Since the flux of energy across the tube section is preserved during reflexion, the energy flux increases after reflexion. A particularly interesting situation occurs when the angle that the inci-

<footnote n="4">⁴The Snell-Descartes law valid for non dispersive waves is recovered only when the slope is horizontal or vertical.</footnote>
<figure><img src="image_27.png" /><figcaption>**Figure 27.** Reflexion of an internal wave on a sloping wall near the critical angle. The angle $\theta$ of propagation with respect to the vertical is conserved, so the reflexion angle $\epsilon_r$ is much smaller than the incidence angle $\epsilon_i$, resulting in wave focusing (from Staquet and Sommeria, 2002).</figcaption></figure>

dent ray makes with the horizontal is equal to the angle of the slope ("critical angle"). In this case, the linear theory predicts a reflected wave of infinite amplitude, infinitesimal wavelength and zero group velocity. The Reynolds number is actually conserved after reflexion (the wavelength decreasing in inverse proportion to the wave-induced velocity) so that, if the incident Reynolds number is large, turbulent effects are expected to occur in the reflected ray tube. This situation has been investigated theoretically using a weakly nonlinear approach (Dauxois and Young, 1999; Tabaei et al., 2005), experimentally (Ivey et al., 2000; McPhee-Shaw and Kunze, 2002; Peacock and Weidman, 2005) and with numerical simulations (Slinn and Riley, 1998; Zikanov and Slinn, 2001). Dauxois and Young (1999) predicted the generation of a harmonic component, which was recently confirmed by the laboratory experiment of Peacock and Weidman (2005), and the formation of an array of counter-rotating vortices along the wall. Slinn and Riley (1998) showed that a strongly turbulent boundary layer can form along the wall with efficient mixing property and, according to the laboratory experiment of McPhee-Shaw and Kunze (2002), this results in irreversible mass transport along the isopycnals and therefore away from the boundary, in so-called nepheloid layers (Figure 28) (see Staquet and Sommeria, 2002, for further discussion).

In a closed domain in which one wall (at least) is inclined to the vertical, the solutions of the linear equation are not longer regular solutions (seiches) but have a fractal structure. In this situation, the incident ray tubes focus at each reflexion on the sloping wall. It has been shown by Maas and Lam (1995) that the focused tubes actually converge toward an "attractor", namely toward a localized region in space (Figure 29a). A laboratory experiment has been performed to check the validity of these theoretical ideas (Maas et al., 1997): the occurence of the attractor is clearly visible in Figure 29b. The occurence of turbulence and mixing along the attractor is still an open question.
<figure><img src="image_28.png" /></figure>

**Figure 28.** Laboratory experiments of an internal gravity wave reflecting on an inclined wall. Fluoresceine dye at the wall is transported along isopycnals, away from the wall, as a result of the breaking of the wave (from McPhee-Shaw and Kunze, 2002).

## 6.2 The internal tide

The internal tide is the field of internal gravity waves generated in the ocean by the interaction of the barotropic tide (the common tide) with the bathymetry. Munk and

<figure><img src="image_29.png" /></figure>

**Figure 29.** (a) Sketch of the attractor: the successive reflexions of an internal wave in a closed domain with one inclined wall converge toward a limited region, referred to as an attractor. (b) Laboratory experiment of the attractor in a tank stably-stratified with salted water. Horizontal lines of dye have been introduced at initial time, whose subsequent motions attest of the internal wave formation. To make the attractor more visible, the image of the flow at initial time has been substracted to that at a time where the attractor has formed. The localized presence of dye lines attest of localized fluid motions (from Maas et al., 1997).
Wunsch (1998) suggested that the breaking and mixing of the internal tide could fundamentally contribute to the global thermal equilibrium of the ocean: the cold bottom water masses in the thermohaline circulation indeed raise toward the surface through mixing; and mixing in the abyssal ocean, away from the boundaries, mostly occurs through internal wave breaking. Several questions therefore stem from this argument: what is the rate of conversion from the barotropic tide into the baroclinic tide? How and where does the internal tide lose energy? In particular, where does breaking occur? Since the conversion of energy is largest at continental slopes, we considered the generation of the internal tide at a continental slope.

The numerical modelling of this problem is not simple because the internal tide dynamics are nonlinear, non-hydrostatic, three-dimensional and a bathymetry is involved. Moreover, a continental slope typically extends from a depth of about 200m down to the abyssal plain, at 4000m. We are currently studying this problem with the code developed by John Marshall at MIT (Marshall et al., 1997) with the investigation of breaking and mixing as the main objective.

Figure 30a displays the internal tide velocity field for a linear hydrostatic computation performed by Gerkema (2002). This simpler context allows one to focus on the generation and propagation of the internal tide. The Brunt-Väisälä frequency is constant in this situation and equal to $2 \times 10^{-3}$ rad/s, the Coriolis frequency is that of mid-latitudes ($10^{-4}$ rad/s) and the topography has a simple shape. Figure 30a has two remarkable features: the internal tide is generated at the shelf break and manifests itself as a single ray reflecting at the bottom and at the surface. Simple qualitative arguments can be provided to account for these features. Because the tide must jump over the slope, its velocity field acquires a vertical velocity over the slope which is largest on top of the slope. Let us consider a frame of reference attached to the barotropic tide. In this frame the topography has an oscillatory motion at the tidal frequency. This is analogous to the motion of an oscillating body in a stably-stratified medium already discussed. An internal gravity wave field is thus generated in the region of largest vertical velocity (if the oscillating frequency is comprised between the local values of N and f), whose energy propagates along a cross with angle set by the dispersion relation. In the present case, only parts of the cross can form because the body (the topography) is not spatially bounded, which propagate toward the abyssal plain and toward the shelf.

The inclusion of nonlinear non-hydrostatic effects strongly modify this regular pattern (Figure 30b). As discussed by Gerkema et al. (2005), an internal wave beam in an infinite medium is a solution of the fully nonlinear Boussinesq equations (see also Tabaei and Akylas, 2003). Hence, nonlinear effects are expected to occur where the beam meets a boundary. Figure 30b shows indeed that small-scale structures are observed in the generation region and where the beam meets the abyssal plain. One can explain the former observation by noticing that the beam is tangent to the slope in the generation region and is therefore aligned with the barotropic flow, so that forcing by the latter flow is optimum. In the reflexion region, higher harmonics are generated as a result of the superposition of the incident and reflected beams. As discussed in Section 6.1, the generation of a second harmonic, predicted theoretically by Dauxois and Young (1999), has been observed in a laboratory experiment of Peacock and Weidman (2005) as a result of the reflexion of an internal gravity beam onto a sloping boundary. No higher harmonics
<figure><img src="image_30.png" /><figcaption>**Figure 30.** (a) Horizontal velocity field at 10 tidal periods, using a linear hydrostatic numerical code (from Gerkema, 2002). (b) Same as (a), using the nonlinear non-hydrostatic MIT code.</figcaption></figure>

are observed however in the experiment, as opposed to the present simulation, because the internal wave beam is not forced.

It follows that higher frequencies are generated once the beam reflects, and this is attested by frequency spectra of the kinetic and potential energy (Figure 31). The slope spectrum is close to -2 and is made of quasi-linear waves, suggesting that the spectra may be of the Garrett-Munk type (see next Section).

A practical application of this study is that two well-distinct regions appear to be turbulent, as further discussed in Gerkema et al. (2005), so that any parameterization of the internal tide may be sought for independently in each region. Of course, this point should be complemented by simulations in a more realistic context.

## 7 Statistical properties of “wave turbulence”

As we recall it in Section 3.3, internal gravity waves transport the momentum and energy of their source, without any transport of matter as long as their dynamics are linear and non-dissipative (see, f.i., Andrews et al., 1987). The degree of nonlinearity depends upon the length scales involved and is quantified by the wave steepness (defined in Section 4.1).

In the ocean, $s \ll 1$ for vertical wavelengths comprised between a few tens and a few hundreds of meters and the wave dynamics are weakly nonlinear<fnref n="5" /> (Garrett and Munk, 1979). Velocity and temperature spectra measured in the ocean at these scales remarkably display the same dependency on the vertical wavenumber and frequency

<footnote n="5">⁵At such scales, Coriolis effects are important and the Rossby number $Ro = A_u |\vec{k}|/f$ should actually be used in place of the Froude number. These large scales are characterized by $Ro \ll 1$. This implies that $Fr \ll 1$, since $Fr = Ro(f/N)$ and $f/N \ll 1$ in geophysical fluids. The Froude number becomes the relevant parameter at smaller scales.</footnote>
whatever the measurement location, as long as this location is far away from sources and sinks of the waves. These spectral features are very well accounted for by a semi-empirical model proposed by Garrett and Munk (1979), based upon the assumption that the waves are linear and influenced by the Coriolis force. The so-called Garrett-Munk spectrum depends separately upon $k_z$ and upon the wave intrinsic frequency $\Omega$ and scales as $k_z^{-2}$ for high vertical wave numbers and as $\Omega^{-2}$ for high frequencies $\Omega$.

For vertical wavelengths between a few meters and a few tens of meters, the wave dynamics become nonlinear ($Fr < 1$) but the spectral features still display a universal dependency, which is distinct from the Garrett-Munk spectrum (Holloway, 1980); remarkably, the kinetic and (available) potential energy spectra display the same law and level as a function of the vertical wave number $k_z$, whatever the flow conditions: $E(k) \approx 0.2N^2k_z^{-3}$. These nonlinear waves break when their Froude number becomes greater than 1. Breaking produces motions close to the dissipation scale, at which fluid mixing occurs. It is important to realize that the waves irreversibly modify their environment when they break, through a transport of mass both along the vertical direction (i.e. through fluid mixing$^6$) and along the horizontal direction (as a mean flow is induced through deposition of momentum) (McIntyre and Norton, 1990).

The ubiquitous occurrence of the latter spectrum is illustrated in Figure 32, through spectra of the available potential energy. In frame (a), data collected from oceanic measurements of the vertical temperature gradient are displayed as a function of $k_z$. Note that the spectrum is nearly flat for small $k_z$ (the exponent is 0.1), corresponding to the Garrett-Munk spectrum. Frame (b) displays very high resolution measurements of the horizontal fluctuating velocity field in the stratosphere. The spectrum is of the

<figure><img src="image_31.png" /><figcaption>**Figure 31.** Frequency kinetic and potential energy spectra of the internal tide field (a) before reflexion of the wave beam; (b) after reflexion of the wave beam. The spectra are averaged over the vertical direction (from Gerkema et al., 2005).</figcaption></figure>

<footnote>$^6$Because of the existence of the restoring force of buoyancy along the vertical direction, no net vertical transport of mass can occur in a stably stratified fluid without mixing.</footnote>
Kolmogoroff type ($\sim k^{-5/3}$) for a vertical wavelength smaller than $\simeq 10$m. In frame (c), the spectrum stems from a laboratory experiment of mixing of a density interface by an oscillating grid (Hannoun and List, 1988). For a certain parameter regime, internal gravity waves develop at the interface and break. In frame (d), the result from a two-dimensional (in a vertical plane) direct numerical simulation of a large scale breaking internal gravity wave in a constant $N$ stratified fluid at rest is presented (Bouruet-Aubertot et al., 1996). The spectra are averaged over 14 Brunt-Väisälä periods (BVP) around the breaking time ($t \simeq 35$ BVP). An identical $k_z^{-3}$ spectrum is obtained when the computation is performed in three dimensions (Carnevale et al., 2001).

# 8 Mixing by breaking internal gravity waves

The importance of mixing stems from the fact that small-scale processes, such as internal gravity waves, have an impact upon their environment only through irreversible processes, such as mixing. Hence, in large scale circulation models, these small-scale motions are usually represented through their mixing properties, via a turbulent diffusivity. Internal gravity waves deposit momentum when they break so that this effect should also be taken into account through a (isopycnal) transport coefficient.

It is not possible to ignore mixing from these small scale motions. As we already stressed it, in the ocean for instance, the cold bottom water masses raise through mixing and it is very likely that the thermohaline circulation, with typical time scale of one thousand years, would stop without mixing, with typical time scale of one minute. As well, in the atmosphere, the climate dynamics strongly depend upon mixing processes, which arise in convective clouds (Bony et al., 1995).

In this section, we only address basic concepts of mixing. A more extensive presentation along with measurements of mixing in a few laboratory and numerical experiments and in geophysical fluids can be found in Staquet (2004).

## 8.1 Basic concepts

Motions in a stably-stratified fluid are associated with kinetic and available potential energy (APE). The APE is partly reversibly converted into kinetic energy through the buoyancy flux, and partly irreversibly converted into background potential energy (BPE) through mixing (Lorenz, 1955; Thorpe, 1977; Winters et al., 1995). The BPE is the potential energy of the background density field $\rho_b(z)$, namely the density field that would exist in the fluid if it were brought to rest. The BPE is therefore the minimum potential energy of the fluid. The physical mechanism to account for the irreversible conversion from APE into BPE is simple: when fluid motions occur, fluid particles of different density are brought into contact and mix through molecular effects. This point can be discussed further. Mixing between two particles depends upon the time the particles are in contact. In other words, the Reynolds number of the fluid motion that advects the particles comes into play: if fluid particles are transported by large scale motions, with large Reynolds number, the diffusive time will be large relative to the advective time and mixing hardly occurs (but stirring does). By contrast, small scale motions will lead to more mixing. These ideas are illustrated in Staquet (2000) for the
<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 32.** Evidence of $k_z^{-3}$ potential energy spectra $E_p(k)$ in the presence of internal gravity wave breaking. a) Schematic spectra of vertical temperature gradients $k_z^2 E_p(k)$ in the oceanic thermocline (from Gregg, 1987). b) Stratospheric temperature spectra measured from instrumented balloons, showing a $k_z^{-3}$ range with a Kolmogoroff $k_z^{-5/3}$ tail at high $k_z$ (provided by F. Dalaudier). c) Normalized internal wave spectra, computed from fluoresceine concentration profiles in a laboratory experiment of mixing of a density interface by an oscillating grid. $\Delta b$ is the buoyancy jump across the density interface, $l_0$ is an integral scale of turbulence and $R_j$ is a Richardson number defined by $\Delta b l_0 / u_0^2$, with $u_0$ being the rms of the horizontal velocity measured in a homogeneous fluid at the same distance from the oscillating grid as the density interface (from Hannoun and List, 1988). d) Two-dimensional direct numerical simulations of a large scale breaking internal gravity wave ($k_y$ is the vertical wavenumber) (from Bouruet-Aubertot et al., 1996).</figcaption></figure>
stably-stratified shear layer.

At the scale of the fluid system, mixing is quantified by the rate of change of the BPE. This rate of change is always positive, a consequence of the second law of thermodynamics. This rate is also the dissipation rate of the APE and will be referred to as the mixing rate $\epsilon_{APE}$. The efficiency of mixing is characterized by the ratio of this rate of change to the amount of energy brought to the system per time unit. If a statistically steady state is reached, the amount of energy brought to the system per time unit is equal to the rate of energy dissipated per time unit in the system, namely $\epsilon_{APE} + \epsilon_{KE}$. The efficiency of mixing is therefore $\frac{\epsilon_{APE}}{\epsilon_{APE} + \epsilon_{KE}}$ and is referred to as the flux Richardson number in the oceanographic community (see Toole, 1998, for a discussion of mixing in the oceanographic context). The crucial point is therefore to compute $\epsilon_{APE}$.

## 8.2 How to compute the mixing rate $\epsilon_{APE}$?

An exact expression for the APE is presented in Holliday and McIntyre (1981), in terms of an expansion that holds for general $N(z)$. For constant $N$ only the first term of the expansion remains, giving $APE = -\frac{g}{2d\bar{\rho}/dz} \langle \rho'^2 \rangle$, where $\langle \rangle$ denotes a volume average and $\bar{\rho}(z)$ is the horizontally averaged density profile. Note that for a constant $N$, the scale of the vertical gradient of the density profile is theoretically infinite so that this case corresponds to the linear limit. From a practical point of view, this expression of the APE is also used when the displacement of the fluid particles is small with respect to the scale of the vertical density gradient $d\bar{\rho}/dz$. In this case, from the linearized Boussinesq equations, the dissipation rate of APE is simply $-\kappa \frac{g}{d\bar{\rho}/dz} \langle |\nabla \rho'|^2 \rangle$.

This method cannot be used when the particle displacements are not small with respect to the scale of the density gradient (but the Boussinesq approximation still holds). The stably-stratified shear layer is one example. In this case, Winters et al. (1995) proposed an exact analytical expression for $\epsilon_{APE}$, based upon the background density profile. Hence, one has to compute this profile and two methods have been proposed for this purpose. The first method was proposed by Thorpe (1977) and extended to two and three dimensions by Winters et al. (1995). The background density profile is obtained by sorting adiabatically the density field so as to get a stable profile (note that the energy released by the sorting process is the APE). From a practical point of view, one sorts the values of the density field at each grid point so that a stable profile is reached; each plane of the numerical domain (supposed to be parallelepipedic), starting from the lowest altitude, is then filled with the heaviest particles and the density is averaged over each plane to get a monotonic density profile. Tseng and Ferziger (2001) proposed another method to access the background density profile, which is more efficient numerically than the sorting method and relies upon the probability density function of the density field.

As said above, these theoretical approaches have been developed to compute the mixing efficiency. Mixing is also quantified by a (diapycnal) turbulent diffusivity $K_\rho$, defined as $\epsilon_{APE}/g = -K_\rho \frac{d\rho_b}{dz}$. In the linear limit, $\bar{\rho} \approx \rho_b$ for all times and one recovers
the expression widely used in oceanography: $K_ρ = -κ \frac{\langle |\nabla ρ'|^2 \rangle}{(d\bar{ρ}/dz)^2}$.

## 9 Conclusion

We tried in these Notes to broadly describe the very rich dynamics of internal gravity waves. The very study of those waves started in the mid-sixties, when Ball (1964) realized that two surface waves could resonantly interact with an internal gravity wave. Much has been learned about these waves, especially from a theoretical point of view, although the dispersion relation of inertia-gravity waves has just been verified experimentally (Peacock and Tabaei, 2005). And much remains to be done when natural media are considered. The sources of internal gravity waves in the atmosphere, other than the wind blowing over orography, are still the subject of active research. For instance, convective clouds at low latitudes generate such waves, which may interact in turn with the clouds. As well, the analysis of atmospheric data (FASTEX campaign) revealed that energetic large-scale low-frequency waves, with velocity comprised between 8 and 10 m/s, are emitted by the tropospheric jet (Plougonven et al., 2003). The important role of internal gravity waves in mixing the abyssal water masses has been conjectured but not proved. And the role of these waves in the dynamics of the radiative zone of the Sun is totally unknown. From a modelling point of view, weakly nonlinear theories have been developed to account for the statistical properties of those waves (see Staquet and Sommeria, 2002). But, from a practical point of view, no reliable parameterization of the waves exists, which would relate the transport coefficients they induce to the dynamical parameters of the large scale motions which create them. Such parameterizations are crucially needed in general circulation models of the atmosphere and the ocean.

## 10 Acknowledgments

The author thanks V. Armenio and S. Sarkar for their invitation to the CISM Summer School on "Environmental stratified flows" in Udine, for very fruitful discussions and for providing the opportunity to write these lecture notes.

## Bibliography

D. G. Andrews, J. R. Holton, and C. B. Leovy. *Middle atmosphere dynamics*. Academic Press, 1987.

S. I. Badulin and V. I. Shrira. On the irreversibility of internal-wave dynamics due to wave trapping by mean flow inhomogeneities. part 1. local analysis. *J. Fluid Mech.*, 251:21-53, 1993.

S. I. Badulin, V. I. Shrira, and L. S. Tsimring. The trapping and vertical focusing of internal waves in a pycnocline due to the horizontal inhomogeneities of density and currents. *J. Fluid Mech.*, 158:199-218, 1985.

S. I. Badulin, V. M. Vasilenko, and N. N. Golenko. Transformation of internal waves in the equatorial lomonosov current. *Atmospheric and Oceanic Physics.*, 26(2):110-117, 1990.
K. Ball. Energy transfer between external and internal gravity waves. *J. Fluid Mech.*, 19:465-478, 1964.

J. Barat and C. Cot. Wind shear rotary spectra in the atmosphere. *Geophys. Res. Letters.*, 19(2):103-106, 1992.

A. Y. Basovich and L. S. Tsimring. Internal waves in a horizontally inhomogeneous flow. *J. Fluid Mech.*, 142:233-49, 1984.

D. Benielli and J. Sommeria. Excitation and breaking of internal gravity waves by parametric instability. *J. Fluid Mech.*, 374:117-144, 1998.

S. Bony, J.-P. Duvel, and H. Le Treut. Observed dependence of the water vapor and clear-sky greenhouse effect on sea-surface temperature: comparison with climate warming experiments. *Climate Dynamics*, 11:307-320, 1995.

J. R. Booker and F. P. Bretherton. The critical layer for internal gravity waves in a shear flow. *J. Fluid Mech.*, 27:513-539, 1967.

P. Bouruet-Aubertot, J. Sommeria, and C. Staquet. Stratified turbulence produced by internal wave breaking: two-dimensional numerical experiments. *Dyn. Atmos. Oceans*, 23:357, 1996.

F. P. Bretherton. The propagation of groups of internal gravity waves in a shear flow. *Quart. J. Roy. Met. Soc.*, 92:466-480, 1966.

F. P. Bretherton. Propagation in slowly varying waveguides. *Proc. Roy. Soc. A*, 302: 555-576, 1968.

G. F. Carnevale, M. Briscolini, and P. Orlandi. Buoyancy- to inertial-range transition in forced stratified turbulence. *J. Fluid Mech.*, 427:205-239, 2001.

S. Chandrashekar. *Hydrodynamic and hydromagnetic stability*. Oxford: Clarendon press, 1961.

B. Cushman-Roisin. *Introduction to Geophysical Fluid Dynamics*. Prentice-Hall, 1994.

F. Dalaudier, C. Sidi, M. Crochet, and J. Vernin. Direct evidence of "sheets" in the atmospheric temperaturefield. *J. Atmosph. Sciences*, 51:237-248, 1994.

T. Dauxois and W. R. Young. Near-critical reflection of internal waves. *J. Fluid Mech.*, 390:271-295, 1999.

P. G. Drazin and W. H. Reid. *Hydrodynamic Stability*. Cambridge University Press, Cambridge, 1981.

P.G. Drazin. On the parametric instability of an internal gravity wave. *Proc. R. Soc. Lond. A*, 356:411, 1977.

T. J. Dunkerton. Shear instability of internal inertia-gravity waves. *J. Atmos. Sciences*, 54:1628, 1997.

N. R. E. Edwards and C. Staquet. Focusing of an inertia-gravity wave packet by a baroclinic shear flow. *Dynamics of Atmospheres and Oceans*, In press, 2005.

A. Eliassen and E. Palm. On the transfer of energy in stationary mountain waves. *Geofys. Publ.*, 22(3):1-23, 1961.

D. M. Farmer and L. Armi. Stratified flow over topography: the role of small scale entrainment and mixing in flow establishment. *Proc. Roy. Soc., London A*, 455:3221-3258, 1999.
D.C. Fritts, J.R. Isler, and O. Andreassen. Gravity wave breaking in two and three dimensions. 2. Three-dimensional evolution and instability structure. *J. Geophys. Res.*, 99(D4):8109–8123, 1994.

C. Garrett and W. Munk. Internal waves in the ocean. *Ann. Rev. Fluid Mech.*, 11: 339–369, 1979.

T. Gerkema. Application of an internal-tide generation model to baroclinic spring-neap cycles. *J. Geophys. Res.*, 107(C9):3124, 2002.

T. Gerkema, C. Staquet, and P. Bouruet-Aubertot. Nonlinear effects in internal-tide beams, and mixing. *Ocean modelling*, Accepted for publication, 2005.

M. Gregg. Diapycnal mixing in the thermocline: a review. *J. Geophys. Res.*, 92(C5): 5249–5286, 1987.

I. A. Hannoun and E. J. List. Turbulent mixing at a shear-free density interface. *J. Fluid Mech.*, 189:211–234, 1988.

K. Hasselmann. On the non-linear energy transfer in a gravity-wave spectrum. Part 1: General theory. *J. Fluid Mech.*, 12:481–500, 1962.

K. Hasselmann. A criterion for nonlinear wave stability. *J. Fluid Mech.*, 30:737, 1967.

P. Hazel. Numerical studies of the stability of inviscid stratified shear flows. *J. Fluid Mech.*, 51:39–61, 1972.

G. J. F. van Heijst, R. C. Kloosterziel, and C. W. M. Williams. Laboratory experiments on the tripolar vortex in a rotating fluid. *J. Fluid Mech.*, 225:301–331, 1991.

D Holliday and M. E. McIntyre. On potential energy density in an incompressible stratified fluid. *J. Fluid Mech.*, 107:221–225, 1981.

G. Holloway. Oceanic internal waves are not weak waves. *J. Phys. Oceanogr.*, 10:906–914, 1980.

G. Holloway and D. Ramsden. Theories of internal wave interaction and stably stratified turbulence: testing against numerical simulations. In J. C. Nihoul and B.M. Jamard, editors, *Small-scale Turbulence and Mixing in the Ocean*, 'Aha Huliko'a Hawaiian Winter Workshop, pages 363–377. Elsevier, New-York, 1988.

L. N. Howard. Note on a paper of john w. miles. *J. Fluid Mech.*, 10:509–512, 1961.

Y. A. Ivanov and Y. G. Morozov. Deformation of internal gravity waves by a stream with horizontal shear. *Okeanologie*, 14:457–461, 1974.

G. N. Ivey, K. B. Winters, and I. P. D. De Silva. Turbulent mixing in a sloping benthic boundary layer energized by internal waves. *J. Fluid Mech.*, 418:59–76, 2000.

W. J. Jones. Propagation of internal gravity waves in fluids with shear flow and rotation. *J. Fluid Mech.*, 30:439–449, 1967.

W. J. Jones. Reflexion and stability of waves in stably stratified fluids with shear flow: a numerical study. *J. Fluid Mech.*, 34:609–624, 1968.

P. Klein, B. L. Hua, and X. Carton. Emergence of cyclonic structures due to the interaction between near-inertial oscillations and mesoscale eddies. *Quart. J. Roy. Met. Soc.*, 129:1–20, 2003.

J. Klostermeyer. Parametric instabilities of internal gravity waves in boussinesq fluids with large reynolds numbers. *Geophys. Astrophys. Fluid. Dyn.*, 26:85, 1983.

J. Klostermeyer. Two and three-dimensional parametric instabilities in finite amplitude internal gravity waves. *Geophys. Astrophys. Fluid Dyn.*, 61:1, 1991.
C. G. Koop. A preliminary investigation of the interaction of internal gravity waves with a steady shearing motion. *J. Fluid Mech.*, 113:347-386, 1981.

C. G. Koop and B. McGee. Measurements of internal gravity waves in a continuously stratified shear flow. *J. Fluid Mech.*, 172:453-480, 1986.

C. R. Koudella. *Ondes internes de gravité : instabilités, déferlement et vorticité potentielle*. PhD thesis, École Normale Supérieure de Lyon, France, 1999.

C. R. K. Koudella and C. Staquet. Instability mechanisms of a two-dimensional progressive internal gravity wave. *J. Fluid Mech.*, submitted, 2005.

P. H. Leblond and L. A. Mysak. *Waves in the Ocean*. Number 20 in Elsevier Oceanography Series. Elsevier Scientific Publishing Company, 1978.

M. P. Lelong and T. J. Dunkerton. Inertia-gravity wavesbreaking in three dimensions: 2. convectively unstable waves. *J. Atmos. Sci.*, 55(15):2489, 1998.

M. Lesieur. *Turbulence in Fluids*. Kluwer Academic Publishers, 3 edition, 1997.

M.J. Lighthill. *Waves in Fluids*. Cambridge University Press, 1978.

P.N. Lombard and J.J. Riley. Instability and breakdown of internal gravity waves. I. linear stability analysis. *Phys. Fluids*, 8(12):3271, 1996.

E.N. Lorenz. Available potential energy and the maintenance of the general circulation. *Tellus*, 7:157-167, 1955.

F. Lott. Comparison between the orographic response of the ecmwf model and the pyrex 1990 data. *Quarterly Journal of the Royal Meteorological Society*, 121:1323-1348, 1995.

L. R. M. Maas, D. Benielli, J. Sommeria, and F.-P. A. Lam. Observation of an internal wave attractor in a confined, stably stratified fluid. *Nature*, 388:557, 1997.

L. R. M. Maas and F.-P. A. Lam. Geometric focusing of internal waves. *J. Fluid Mech.*, 300:1-41, 1995.

J. Marshall, A. Adcroft, C. Hill, L. Perelman, and C. Heisey. A finite-volume, incompressible navier-stokes model for studies of the ocean on parallel computers. *J. Geophys Res.*, 102(C3):5753-5766, 1997.

C. H. McComas and F. P. Bretherton. Resonant interactions of oceanic internal waves. *J. Geophys. Res.*, 82:1397-1411, 1977.

A.D. McEwan. Degeneration of resonantly-excited standing internal gravity waves. *J. Fluid Mech.*, 50:431-448, 1971.

A.D. McEwan and R.A. Plumb. Off-resonant amplification of finite internal wave packets. *Dyn. Atmos. Ocean*, 2:83-105, 1977.

A.D. McEwan and R. M. Robinson. Parametric instability of internal gravity waves. *J. Fluid Mech.*, 67(4):667-687, 1975.

M. E. McIntyre. The stratospheric polar vortex and sub-vortex : fluid dynamics and mid-latitude ozone loss. *Phil. Trans. Roy. Soc. London*, 352:227-240, 1995.

M. E. McIntyre. On global-scale atmospheric circulations. In G. K. Batchelor, H. K. Mofatt, and M. G. Worster, editors, *Perspectives in Fluid Dynamics*, pages 557-624. Cambridge University Press, 2000.

M. E. McIntyre and W. A. Norton. Dissipative wave-mean interactions and the transport of vorticity or potential vorticity. *J. Fluid Mech.*, 212:403-435, 1990.
E.E. McPhee-Shaw and E. Kunze. Boundary-layer intrusions from a sloping bottom: A mechanism for generating intermediate nepheloid layers. *J. Geophys. Res.*, 107: DOI:10.1029/2001JC000801, 2002.

R.P. Mied. The occurrence of parametric instabilities in finite-amplitude internal gravity waves. *J. Fluid Mech.*, 78:763-784, 1976.

T. Miyazaki and K. Adachi. Short-wavelength instabilities of waves in rotating stratified fluids. *Phys. Fluids*, 10(12):3168-3177, 1998.

D. E. Mowbray and B. S. H. Rarity. A theoretical and experimental investigation of the phase configuration of internal waves of small amplitude in a density stratified liquid. *J. Fluid Mech.*, 28:1-16, 1967.

W. Munk and C. Wunsch. Abyssal recipes ii: energetics of tidal and wind mixing. *Deep-Sea Research I*, 45:1977, 1998.

W. H. Munk. A survey of internal waves and small scale processes. In B. A. Warren and C. Wunsch, editors, *Evolution of Physical Oceanography*, pages 264-291. MIT Press, Cambridge, Mass., 1981.

P. Muller, G. Holloway, F. Henyey, and N. Pomprey. Nonlinear interactions among internal gravity waves. *Rev. Geophys.*, 24(3):493-536, 1986.

D. J. Olbers. The propagation of internal waves in a geostrophic current. *J. Phys. Oceanogr.*, 11:1224-1233, 1981.

M. C. Ollers, L. P. J. Kamp, F. Lott, P. F. J. van Velthoven, H. M. Kelder, and F. W. Sluijter. Propagation properties of inertia-gravity waves through a barotropic shear layer and application to the antarctic polar vortex. *Quartely Journal of the Royal Meteorological Society*, 129:2495-2911, 2003.

T. Peacock and A. Tabaei. Visualization of nonlinear effects in reflecting internal wave beams. *Phys. Fluids*, Accepted for publication, 2005.

T. Peacock and P. Weidman. The effect of rotation on conical wave beams in a stratified fluid. *Expts. Fluids*, To appear, 2005.

O. M. Phillips. On the dynamics of unsteady gravity waves of finite amplitude, i. *J. Fluid Mech.*, 9:193-217, 1960.

O. M. Phillips. Wave interactions - the evolution of an idea. *J. Fluid Mech.*, 106:215-227, 1981.

O.M. Phillips. *The dynamics of the upper ocean*. Cambridge University Press, 1966.

R. Plougonven, H. Teitelbaum, and V. Zeitlin. Inertia-gravity wave generation by the tropospheric mid-latitude jet as given by the fastex radiosoundings. *Journal of Geophysical Research*, 108(D21):4686, 2003.

J. M. Rees, J. C. W. Denholm-Price, J. C. King, and P. S. Anderson. A climatological study of internal gravity waves in the atmospheric boundary layer overlying the brunt ice shelf, antarctica. *Journal of the Atmospheric Sciences*, 57:511-526, 2000.

E. Schatzman. Diffusion process produced by random internal waves. *J. Fluid Mech.*, 322:355-382, 1996.

D. G. Schowalter, J. C. Lasheras, and C. W. van Atta. A study of streamwise vortex structures in a stratified shear layer. *J. Fluid Mech.*, 281:247-291, 1994.

D. N. Slinn and J. J. Riley. Turbulent dynamics of a critically reflecting internal gravity wave. *Theor. Comp. Fluid Dyn.*, 11:281-303, 1998.
L. J. Sonmor and G. P. Klaassen. Higher-order resonant instabilities of internal gravity waves. *J.Fluid Mechanics*, 324:1–23, 1996.

C. Staquet. Mixing in a stably stratified shear layer: two- and three-dimensional numerical experiments. *Fluid Dyn. Research*, 27:367–404, 2000.

C. Staquet. Gravity and inertia-gravity internal waves : breaking processes and induced mixing. *Surveys in Geophysics*, 25(3-4):281–314, 2004.

C. Staquet and G. Huerre. On transport across a barotropic shear flow by breaking inertia-gravity waves. *Physics of Fluids*, 14(6):1993–2006, 2002.

C. Staquet and J. Sommeria. Mixing by breaking internal gravity waves: from instabilities to turbulence. *Annual Reviews Fluid Mech*, 34:559–593, 2002.

R. B. Stull. *Boundary layer meteorology*. Kluwer Academic Publisher, 1988.

B. R. Sutherland. Finite-amplitude internal wave-packet dispersion and breaking. *J. Fluid Mech.*, 429:343–380, 2001.

A. Tabaei and T. R. Akylas. Nonlinear internal gravity wave beams. *J. Fluid Mech.*, 482:141–161, 2003.

A. Tabaei, T. R. Akylas, and K. G. Lamb. Nonlinear effects in reflecting and colliding internal wave beams. *J. Fluid Mech.*, 526:217–243, 2005.

S. Talon, P. Kumar, and J.-P. Zahn. Angular momentum extraction by gravity waves in the sun. *Astrophys. J. Letters*, 574:L175–L178, 2002.

S. A. Thorpe. Turbulence and mixing in a scottish loch. *Phil. Trans. R. Soc. Lond. A*, 286:125–181, 1977.

J. Toole. Turbulent mixing in the ocean: intensity, causes and consequences. In E. Chas-signet and J. Verron, editors, *Ocean modelling and parameterization*, Proceedings of Les Houches Winter School, pages 171–190. Kluwer Academic Publishers, 1998.

Y. H. Tseng and J. H. Ferziger. Mixing and available potential energy in stratified flows. *Physics Fluids*, 33(2):1281–1293, 2001.

A. Veysseyre, K. Moutard, C. Ferrari, K. Van de Velde, C. Barbante, G. Cozzi, G. Capodaglio, and C. Boutron. Heavy metals in fresh snow collected at different altitudes in the chamonix and maurienne valleys, french alps: initial results. *Atmo-spheric Environment*, 35(2):415–425, 2000.

K. B. Winters and E. A. D'Asaro. Two-dimensional instability of finite amplitude internal gravity wave packets near a critical level. *J. Geophys. Res.*, 94(C9):12,709–19, 1989.

K. B. Winters and E. A. D'Asaro. Three-dimensional wave breaking near a critical level. *J. Fluid Mech.*, 272:255–284, 1994.

K. B. Winters, P. N. Lombard, J. J. Riley, and E. A. D'Asaro. Available potential energy and mixing in density-stratified fluids. *J. Fluid Mech.*, 289:115–128, 1995.

K. B. Winters and J. J. Riley. Instability of internal waves near a critical level. *Dyn. Atmos. Oceans*, 16:249–278, 1992.

M. G. Wurtele, A. Datta, and R. D. Sharman. The propagation of gravity-inertia waves and lee waves under a critical level. *J. Atmos. Sciences*, 53(11):1505, 1996.

O. Zikanov and D. N. Slinn. Along-slope current generation by obliquely incident internal waves. *J. Fluid Mech.*, 445:235–261, 2001.