# On transport across a barotropic shear flow by breaking inertia-gravity waves

Chantal Staquet, Gaël Huerre

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

## To cite this version:

Chantal Staquet, Gaël Huerre. On transport across a barotropic shear flow by breaking inertia-gravity waves. Physics of Fluids, 2002, 14, pp.1993-2006. 10.1063/1.1476303. hal-00264604

HAL Id: hal-00264604

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

Submitted on 5 Feb 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On transport across a barotropic shear flow by breaking inertia-gravity waves

C. Staquet<sup>a)</sup> and G. Huerre

LEGI, B.P. 53, 38041 Grenoble Cedex 9, France

(Received 22 January 2001; accepted 14 March 2002; published 6 May 2002)

We investigate the interaction of an inertia-gravity wave packet with a barotropic dynamical barrier, by means of three-dimensional direct numerical simulations of the Boussinesq equations in a rotating reference frame. A simple model is used for the barrier, namely a barotropic shear layer with a horizontal shear. We show that the increase of the intrinsic frequency $\Omega$ of the wave packet by Doppler shifting may result in the breaking of the wave packet, in the neighborhood of a trapping plane where $\Omega=N$ ($N$ is the Brunt–Väisälä frequency of the fluid and is constant in the present study). The breaking process triggers a turbulent diffusive transport of mass and momentum across the barrier, with a turbulent Prandtl number of order 1. © 2002 American Institute of Physics.
[DOI: 10.1063/1.1476303]

## I. INTRODUCTION

Dynamical barriers are ubiquitous large-scale features in geophysical flows. Fronts of temperature or the edges of vortices, such as meddies in the ocean or the polar vortex in the stratosphere, are examples of dynamical barriers. They are characterized by a strong local gradient of temperature (for fronts) or of potential vorticity (for vortex edges). These barriers, as they are entitled, are supposedly impermeable to transport by quasi-two-dimensional turbulence but this impermeability has recently been put into question (e.g., McIntyre¹) from observational evidence. The existence of a transport of material across a dynamical barrier is a crucial issue in geophysical fluid dynamics because, for instance, of the possible release of chemical reactants or pollutants originally trapped or isolated by the barrier.

In the present paper, we investigate the possible transport across a dynamical barrier by inertia-gravity waves, following a conjecture proposed by McIntyre¹ within the context of the polar vortex. We shall make use of a very simplified model for the barrier and for the waves so that our work must be considered as a study of the fundamentals of the wave-barrier interaction, motivated by the geophysical considerations discussed above.

The simple model we consider consists in a barotropic current, that is, a vertically uniform horizontal flow with a horizontal shear, toward which an inertia-gravity wave packet propagates, the fluid being uniformly stably stratified (the Brunt–Väisälä frequency $N$ is therefore constant). The interaction of internal gravity waves with such a current has been studied in the literature within an oceanic context, using a linear analysis² and a Wentzel–Kramers–Brillouin (WKB) approximation.³⁻⁶ Our study will thus rely on these theoretical works.

These works display strong analogies with those related

to the vertical propagation of internal gravity waves in a horizontal wind with a vertical shear, a topic extensively studied within an atmospheric context (starting with the work of Bretherton⁷ and Booker and Bretherton;⁸ see Baik et al.⁹ for recent development). The main results of the latter works are as follows. When an internal gravity wave propagates upwards along a vertically sheared horizontal wind whose velocity increases with height, the wave may be trapped near a level at which its horizontal phase speed matches the wind velocity. This is the so-called critical level. As the wave propagates toward that level, its intrinsic frequency $\Omega$ (i.e., the frequency measured by an observer moving with the wind) decreases by Doppler shifting and reaches a minimum value at this level (this value is zero if the Coriolis force is ignored). Hence, the wave cannot propagate further up. Using the WKB approximation, Bretherton⁷ first showed that, for a wave packet, the vertical component of the wave vector increases as $\eta^{-1}$, where $\eta$ is the distance of the wave packet from the critical level (located at $\eta=0$) along the vertical direction and the vertical component of the group velocity decreases as $\eta^2$ [Fig. 1(a)]. However, the WKB approximation also predicts that the wave induced horizontal velocity $u'$ goes to infinity ($u' \sim \eta^{-1/2}$ as $\eta \to 0$) and the critical level is never reached (the time required to reach that level scales like $\eta^{-1}$). This unphysical behavior is due to the WKB approximation to becoming invalid at the critical level. Using the linear Boussinesq approximation, without WKB approximation, Booker and Bretherton⁸ showed that the wave energy actually decays, being totally absorbed into the mean flow (which is therefore accelerated) at the critical level. Viscous and thermal dissipation are neglected in this approach. As pointed out by McIntyre,¹⁰ the horizontal phase speed of the wave packet $c_x$ decays faster than the wave induced horizontal velocity $u'$ as $t \to +\infty$ so that, in spite of this absorption, the ratio $u'/c_x$, which is one measure of the wave steepness, increases. As a consequence, the linear theory becomes invalid for large time and breaking

<footnote>a) Author to whom correspondence should be addressed. Electronic mail: chantal.staquet@hmg.inpg.fr</footnote>
<figure><img src="image_3.png" /></figure>

FIG. 1. Interaction of an internal gravity wave of absolute frequency $\omega_0$ and wave vector $\vec{k}$ with a horizontal flow $U_b \vec{t}_x$ with either a vertical shear (a) or a horizontal shear (b), in a stably stratified fluid with constant Brunt–Väisälä frequency $N$ ($f=0$). The shear is taken constant for simplicity. In (a), the interaction process is two dimensional and occurs in the vertical ($\vec{k}, \vec{g}$) plane. A critical level is located at $z=z_c$, at which the intrinsic frequency of the wave $\Omega$ vanishes. As the critical layer is approached, $c_{gz} \to 0$, $c_{gx} \to U_b$, $k_x$ remains constant while $k_z$ increases without limit. In (b), the interaction process is three dimensional. It is displayed in a horizontal $(x,y)$ plane; $\vec{k}_h$ and $\vec{c}_{gh}$, respectively, denote the restriction to this plane of the wave vector and of the intrinsic group velocity of the wave. A trapping plane is located at $y=y_t$, at which $\Omega$ is equal to $N$. As the trapping plane is approached, $c_{gy} \to 0$, $c_{gx} \to U_b$, $k_x$ remains constant while $k_y$ increases without limit.

may occur before the critical level is reached. Using three-dimensional numerical simulations, Winters and D'Asaro¹¹ showed that the occurrence of breaking depends upon the initial steepness of the wave packet. For a large enough initial steepness, the wave does break. The details of the breaking process have been investigated by Dörnbrack¹² and Fritts et al.:¹³ breaking first occurs through a convective instability. This result is consistent with the numerical finding by Lelong and Dunkerton¹⁴ that a low frequency internal gravity wave with a high vertical shear and a strong steepness, which are the kinematic properties of a wave packet near a critical level, breaks through convective instability. The behavior of internal gravity waves near a critical level has also been investigated using laboratory experiments by Thorpe¹⁵ and Koop and McGee¹⁶ (see Staquet and Sommeria¹⁷ for a review).

When the horizontal mean flow displays a horizontal shear (such as an oceanic current) and increases as the wave propagates into the current, a trapping process may also occur when the wave propagates against the current [Fig. 1(b)]. Indeed, in this case, the wave intrinsic frequency $\Omega$ increases by Doppler shifting. If there exists a location in the flow where the intrinsic frequency reaches the Brunt–Väisälä frequency $N$, the wave cannot propagate beyond that location (since $\Omega \le N$ for internal gravity waves, as we recall it in the next section). When $N$ is constant, this location is a vertical plane, hereafter referred to as a trapping plane. Let $\eta$ be again the distance of a wave packet from the trapping plane along the direction of inhomogeneity, which is now horizontal. The WKB approximation dictates that, as $\eta$ goes to zero, the component of the wave vector along that direction increases again without bound, as $\eta^{-1/2}$, and the group velocity decreases as $\eta^{3/2}$ (see Ref. 3). Like for the critical level situation, the wave energy increases without bound ($\sim \eta^{-3/2}$) and the trapping plane is never reached, due to the invalidity of the WKB approximation at that plane. Hints about the actual wave behavior near the trapping plane may be obtained from the linear analysis about the mean shear state performed by Ivanov and Morozov.² The numerical resolution of the eigenvalue problem thus obtained yields that the wave amplitude may increase in this case, unlike the

critical level situation. The source of energy is of course provided by the current. As pointed out by Badulin et al.,⁵ the wave packet slows down as it propagates into the current (because the intrinsic group velocity goes to zero), which provides an additional contribution to the local increase of the wave energy. As a result, the wave may break before the trapping plane location. This situation, despite its close analogy to the critical level problem, has never been addressed numerically and the purpose of the present paper is to investigate whether breaking may indeed occur as a result of the interaction with the current and to analyze the resulting transport properties across that current. Badulin et al.⁵ also show that horizontal inhomogeneities of density, namely fronts, result in the same trapping process. The effect of the Coriolis force does not fundamentally change these conclusions.³

When the velocity of the horizontal flow increases as the wave propagates into that flow, reflection occurs when the wave propagates *against* the flow in the vertical shear case (e.g., a wind) and *along* the flow in the horizontal shear case (e.g., a current). This corresponds to an increase of the wave intrinsic frequency in the former case (up to $N$) and to a decrease in the latter case (down to a minimum frequency). This reflection situation will not be considered here (see Sutherland¹⁸,¹⁹ for a theoretical and numerical study of internal gravity wave packets reflecting in a vertically sheared horizontal flow).

In the present paper, the interaction of an inertia-gravity wave packet with a parallel shear flow in a horizontal plane (the barrier) is investigated using three-dimensional direct numerical simulations of the Boussinesq equations in a rotating reference frame. An equation for the transport of a passive scalar is coupled to the equations so as to quantify the possible resulting transport process across the barrier. In the next section, we specify the physical and numerical models we use for the study of the interaction. Section III reports on the breaking process and the resulting transport properties are analyzed in Sec. IV. Conclusions are drawn in the final section.
## II. A SIMPLIFIED MODEL OF THE WAVE-BARRIER INTERACTION

### A. Physical model

We consider a stably stratified fluid medium in hydrostatic equilibrium, within a rotating frame of reference with Coriolis frequency $f$. Let $(x,y,z)$ be a coordinate system in this reference frame, with $z$ pointing upwards along the rotation axis. $\vec{x}=(x,y,z)$ is thus the Eulerian position of a fluid particle in this frame. The unit vectors along the axis of coordinates are denoted $(\vec{i}_x, \vec{i}_y, \vec{i}_z)$. Let $\rho_0+\rho_{\text{lin}}(z)$ be the density field of the hydrostatic equilibrium state, where $\rho_{\text{lin}} \ll \rho_0$ for any $z$ and $\rho_0$ is a constant reference density. The linear density profile $\rho_{\text{lin}}(z)$ sets the Brunt-Väisälä frequency

$$N^2 = - \frac{g}{\rho_0} \frac{d\rho_{\text{lin}}}{dz}, \quad (1)$$

where $g$ is the acceleration of gravity. $N$ is therefore constant in the present study.

The barrier is a horizontal shear flow, whose shear is also horizontal,

$$\vec{U}_b(y) = U_b(y) \vec{i}_x. \quad (2)$$

This flow is in geostrophic equilibrium since the Coriolis force $-f\vec{i}_z \wedge \vec{U}_b$ is balanced by a pressure gradient along the $y$ direction. The barrier is barotropic for simplicity. Indeed, any change of $U_b$ along the vertical direction would yield a horizontal inhomogeneity of the density field, according to the thermal wind balance (e.g., Gill,²⁰ p. 217). Hence, both a density front and a shear flow would be present in the flow and refract the wave. In the present paper, we focus our attention about a single mechanism for the refraction of the wave, so that a barotropic shear flow is considered. The use of a shear flow represents either a horizontal current or the edge of a vortex whose horizontal extent is much larger than the horizontal wavelength of the wave. Two different shear flows have been considered, a shear layer and a plane jet; the same qualitative behavior, with respect to the breaking and transport properties, has been obtained for either flow and results for the shear layer will be presented in this paper. The initial velocity profile we chose for the latter flow is

$$U_b(y) = U_0 \tanh \left( \frac{y - y_1}{a} \right), \quad (3)$$

where $a$ and $U_0$ denote, respectively, half the (so-called) vorticity thickness of the shear layer and half the velocity difference across the shear layer. Note that this shear flow is anticyclonic: the vertical component of its vorticity $-dU_b/dy(<0)$ is of opposite sign to the Coriolis frequency $f$.

An inertia-gravity wave packet of the main wave vector $\vec{k}=(k_x,k_y,k_z)$ is introduced as the initial condition. The wave vector $\vec{k}$ sets the intrinsic frequency $\Omega$ of the wave packet via the dispersion relation

$$\Omega^2 = N^2 \left( \frac{k_h}{k} \right)^2 + f^2 \left( \frac{k_z}{k} \right)^2, \quad (4)$$

<figure><img src="image_4.png" /><figcaption>FIG. 2. Sketch of the initial condition of the simulations: an inertia-gravity wave packet is introduced at $t=0$ about the $y=y_0$ plane and propagates in a horizontal shear layer with no vertical dependency centered about the $y=y_1$ plane.</figcaption></figure>

where $k_h = \sqrt{k_x^2 + k_y^2}$ is the horizontal wave number and $k = |\vec{k}|$. Relation (4) implies that $f \le \Omega \le N$. The group velocity of the wave packet in the rotating reference frame (that is, the absolute group velocity) is $\vec{c}_g + \vec{U}_b$, where $\vec{c}_g = \partial\Omega/\partial\vec{k}$ is the intrinsic group velocity (e.g., Lighthill,²¹ p. 327). In the present case where the shear flow is directed along the horizontal $x$ direction, the group velocity is expressed by $(\partial\Omega/\partial k_x + U_b, \partial\Omega/\partial k_y, \partial\Omega/\partial k_z)$. We also recall that the energy of an internal wave packet is transported at the group velocity along trajectories referred to as rays; a ray path is thus defined by $d\vec{x}/dt = \vec{c}_g + \vec{U}_b$.

The energy of the wave packet is confined about the $y=y_0$ vertical plane (with $y_0<y_1$) at initial time. The initial velocity and density fields induced by the wave packet are thus given by (e.g., Olbers,³ Gill,²⁰ p. 263)

$$\begin{aligned} \dot{\vec{u}}'(\vec{x}, t=0) ={}& A e^{-(y-y_0)^2} e^{i\vec{k} \cdot \vec{x}} \left[ -\frac{k_z}{k_h^2} \left( k_x + i k_y \frac{f}{\Omega} \right), \right. \\ & \left. -\frac{k_z}{k_h^2} \left( k_y - i k_x \frac{f}{\Omega} \right), 1 \right], \end{aligned} \quad (5)$$

$$\rho'(\vec{x}, t=0) = i \frac{N^2}{\Omega} w'(\vec{x}, t=0). \quad (6)$$

$A$ is the initial amplitude of the vertical velocity component $w'$. This simple physical model of the wave-barrier interaction is sketched in Fig. 2.

The wave packet being introduced in a moving fluid (the shear layer), its frequency $\omega_0$ in the rotating frame of reference (the absolute frequency) differs from its intrinsic frequency. $\omega_0$ is defined by the Doppler relation

$$\omega_0 = \Omega + \vec{k} \cdot \vec{U}_b = \Omega + k_x U_b. \quad (7)$$

Because the background shear flow is not homogeneous along the $y$ direction, the properties of the wave packet change as it propagates into this flow. Those changes are easily predicted if the wavelength along the $y$ direction is much smaller than the scale over which $U_b$ varies. When this assumption holds, a WKB approximation can be used which predicts that, at lowest orders, the dispersion relation remains
locally valid but the wave is refracted along rays by the
spatial variability of the background flow. (A very clear in-
troduction to the WKB approximation for internal gravity
waves is provided by Olbers.³) Since U<sub>b</sub> depends only upon
y, the k<sub>x</sub> and k<sub>z</sub> components of the wave packet do not
change along a ray. By contrast, k<sub>y</sub> changes according to

$$
\frac{dk_y}{dt} = -k_x \frac{dU_b}{dy}, \qquad (8)
$$

where $d/dt = \partial/\partial t + (\vec{c}_g + \vec{U}_b) \cdot \nabla$ denotes the rate of change along a ray in the rotating reference frame. Equation (8) implies that $k_y$ increases along the ray when $k_x$ is negative, since $dU_b/dy$ is positive from definition (3). (This choice for the sign of $k_x$ will be made in the next section.) Note that the increase of $k_y$ as the wave propagates toward the barrier implies that the wavelength of the wave packet along this direction becomes increasingly smaller with respect to the characteristic scale of $dU_b/dy$, which enforces the validity of the WKB approximation (but the latter approximation fails in the immediate vicinity of the trapping plane, as already noted in the Introduction). Since $U_b$ does not depend upon time, the absolute frequency $\omega_0$ stays constant along a ray ($d\omega_0/dt=0$). Hence, using the Doppler relation (7), the intrinsic frequency changes along a ray according to

$$
\frac{d\Omega}{dt} = -k_x c_{gy} \frac{dU_b}{dy}, \quad (9)
$$

where $c_{gy}$ is the y component of the intrinsic group velocity.
Since the wave packet propagates toward the barrier, $c_{gy}$ is
positive. Thus, the intrinsic frequency increases along a ray
when $k_x$ is negative. Like $\Omega$, the wave energy $E$ is not con-
served along a ray (it would be if the fluid were at rest or
uniformly translating) but the wave action $A=E/\Omega$ is con-
served instead of the energy$^7$

$$
\frac{\partial A}{\partial t} + \nabla \cdot [(\vec{c}_g + \vec{U}_b)A] = 0, \quad (10)
$$

for a nondissipative and nondiffusive medium.

As noted in the Introduction, the trapping location is a
vertical y=y<sub>t</sub> plane obtained by setting Ω to N in the Dop-
pler relation (7). Using ω<sub>0</sub>=Ω(y<sub>0</sub>) + k<sub>x</sub>U<sub>b</sub>(y<sub>0</sub>), one gets an
expression for the location of the trapping plane as a function
of the initial parameters of the wave

$$
U_b(y_t) = \frac{\Omega(y_0) - N}{k_x} + U_b(y_0). \qquad (11)
$$

In the present case where $U_b$ is monotonically varying,
$y_t$ is uniquely defined by expression (11). This expression
shows that $y_t$ can be located anywhere along the y axis.

B. Mathematical and numerical models

The fluid motions are described by the Navier–Stokes
equations in the Boussinesq approximation within the rotat-
ing reference frame introduced in the preceding section. Let
$\vec{u}=(u,v,w)$ be the components of the velocity field in this
reference frame and $\vec{\rho}(x,y,z,t)$, the deviation of the density
field from the hydrostatic density field $\rho_0+\rho_{\text{lin}}(z)$. As is

customary, we scale $\vec{\rho}$ as an acceleration using the variable $\rho' = (g/\rho_0) \vec{\rho}$ in order for $N$ to come into play in the equations of motion. These equations thus are (written per unit mass)

$$
\frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla) \vec{u} = - \nabla P - \rho' \vec{i}_z - f \vec{i}_z \wedge \vec{u} + \nu \nabla^2 \vec{u}, \quad (12)
$$

$$
\frac{\partial \rho'}{\partial t} + (\vec{u} \cdot \nabla) \rho' = N^2 w + \frac{\nu}{\text{Pr}} \nabla^2 \rho', \quad (13)
$$

$$
\frac{\partial S}{\partial t} + (\vec{u} \cdot \nabla) S = \frac{\nu}{\text{Pe}} \nabla^2 S, \quad (14)
$$

$$
\nabla \cdot \vec{u} = 0. \tag{15}
$$

In Eq. (12), the term *P* stands for the dynamical pressure divided by *ρ*₀ and *ν* is the kinematic viscosity. In Eq. (13), *Pr* is the Prandtl number (i.e., *ν*/*Pr* defines the molecular diffusivity for density changes). Equation (14) describes the evolution of the passive scalar *S*(*x*,*t*); the equation involves the Peclet number *Pe* defined such that *ν*/*Pe* is the molecular diffusivity of the passive scalar.

The initial condition is provided by the physical model
described in Sec. II A: at *t*=0, an inertia-gravity wave packet
localized about the *y*=*y*₀ plane propagates in a horizontal
shear layer centered about the *y*=*y*₁ plane. Periodic bound-
ary conditions are imposed along the *x* and *z* directions and
free slip ones are used along the inhomogeneous *y* direction.

Equations (12)–(15) are solved numerically in a rectangular box. The spatial derivatives are computed by a pseudospectral method in all three directions. Along the y direction, periodicity is recovered by using a classical technique:²² symmetry properties of the fields are assumed about the y boundaries so that the fields become periodic when the domain is doubled by symmetry about these boundaries. Time advancement is performed by a third order Adams–Bashforth scheme. Note that no subgrid scale modelling is used: the terms responsible for the diffusion of momentum, density, and passive scalar [$\nu\nabla^2\vec{u}$, ($\nu$/Pr)$\nabla^2\rho'$ and ($\nu$/Pe)$\nabla^2S$, respectively] involve a constant molecular viscosity $\nu$, whose value is specified below.

C. Physical and numerical parameters of the runs

For simplicity, the physical parameters *a* (half the vorticity thickness of the shear layer), *A* (the initial amplitude of the vertical velocity induced by the wave) and *N* are set to 1, 0.1 and 1, respectively. The other parameters of the problem ($\vec{k}$ and $y_0$ for the wave; $U_0$ and $y_1$ for the barrier; *f*, *v*, *Pr*, and *Pe*) as well as the size of the computational domain are dictated by two sets of constraints. A first set consists of physical constraints: (i) the horizontal shear layer should remain stable against the Kelvin–Helmholtz instability during the propagation and evolution of the wave; (ii) there should exist a trapping plane somewhere in the flow; (iii) the horizontal wavelength of the wave should be larger, by a factor of 10 or so, than the thickness of the shear flow (this ratio is dictated by geophysical application). The second set consists of numerical constraints: (iv) the trapping plane should be located within the computational domain; (v) the wave
<table><caption>TABLE I. Description of the runs. The resolution along the x, y, and z directions, respectively, is indicated. $U_0$ is half the velocity difference across the shear layer, $\nu$ the molecular viscosity, and Pe is the Peclet number associated with the passive scalar. $y_t$ is the position of the trapping plane, defined by Eq. (11) and $y_s$ is the position along the y axis of the maximum concentration of the passive scalar [defined by Eq. (22)]. The other parameters of the runs are described in Sec. II C.</caption><thead><tr><th>Run</th><th>Resolution</th><th>$U_0$</th><th>$\nu$</th><th>Pe</th><th>$y_t$</th><th>$y_s$</th></tr></thead><tbody><tr><td>1a</td><td>(256,257,128)</td><td>1</td><td>$10^{-4}$</td><td>10</td><td>0.55</td><td>1</td></tr><tr><td>1b</td><td>(128,129,64)</td><td>1</td><td>$10^{-4}$</td><td>10</td><td>0.55</td><td>1</td></tr><tr><td>2</td><td>(256,257,128)</td><td>1</td><td>$10^{-4}$</td><td>1</td><td>0.55</td><td>1</td></tr><tr><td>3</td><td>(256,257,128)</td><td>1</td><td>$10^{-3}$</td><td>10</td><td>0.55</td><td>1</td></tr><tr><td>4a</td><td>(256,257,128)</td><td>1</td><td>$10^{-4}$</td><td>10</td><td>0.55</td><td>1.65</td></tr><tr><td>4b</td><td>(128,129,64)</td><td>1</td><td>$10^{-4}$</td><td>10</td><td>0.55</td><td>1.65</td></tr><tr><td>5</td><td>(128,129,64)</td><td>0.5</td><td>$10^{-4}$</td><td>10</td><td>1.15</td><td>1</td></tr><tr><td>6</td><td>(128,129,64)</td><td>2</td><td>$10^{-4}$</td><td>10</td><td>0.12</td><td>1</td></tr><tr><td>7</td><td>(128,129,64)</td><td>3</td><td>$10^{-4}$</td><td>10</td><td>-0.107</td><td>-0.107</td></tr><tr><td>8</td><td>(128,129,64)</td><td>4</td><td>$10^{-4}$</td><td>10</td><td>-0.26</td><td>-0.26</td></tr></tbody></table>

packet should quickly reach the neighborhood of the trapping plane, in order to save computer time; (vi) the $N/f$ ratio should not exceed a value of order 10 or so because a very small vertical wavelength relative to the horizontal one should then be employed, thus imposing a high resolution along the vertical direction; (vii) for a given resolution, the molecular diffusion coefficients should be high enough for the smallest scales along the inhomogeneous y direction to be properly solved (we recall that the wavelength along the y direction decreases as the wave packet propagates toward the barrier).

Table I gives the values chosen in view of these constraints. In all runs, the wave packet is introduced at $t=0$ about location $y_0=-2$ and the main wave vector is $\vec{k} = (-\pi/3, \pi/3, 4\pi/3)$; the barrier is located at $y_1=1$ and $U_0$ is equal to 1 except in runs 5–8 where it is varied from 0.5 to 4. In all runs, $\nu=10^{-4}$ (except in run 3), the Coriolis parameter $f=0.23$ (so that $N/f=4.35$) and Pr=1. The Peclet number is equal to 10 (except in run 2). The size of the computational domain is $(L_x, L_y, L_z)=(12,12,6)$. Preliminary computations at resolution (128,129,64) have been performed in order to explore the parameter space. Several computations have next been carried out at the higher resolution of (256,257,128). The interaction between the wave and the barrier displays the same features from one run to the other for a given resolution and the same dynamical properties are recovered when the resolution is increased. Results for run 1a are thus presented in the paper and the influence of a varying viscosity, Peclet number and shear $dU_b/dy$ will be discussed in connection with these results.

### III. THE BREAKING PROCESS

#### A. Overall behavior of the wave packet

The overall dynamics of the wave packet in a vertical $(y,z)$ plane are illustrated in Figs. 3 and 4 through contours at successive times of the total density field $-N^2z+\rho'$. Figure 3 displays the density field before breaking occurs and Fig. 4 illustrates the turbulent regime that follows breaking. The shear flow is not visible in these figures because, having no component in the $(y,z)$ plane, it does not displace the isopycnals. All properties of the wave in this plane are thus

intrinsic properties. Contours of the entrophy $Z$ are plotted in Fig. 5 at successive times in the same $(y,z)$ plane to illustrate the wave-barrier interaction (both the shear flow and the wave packet being rotational); the maximum value of $Z$ over the numerical domain is displayed versus time in Fig. 6. $Z$ is defined by $0.5|\xi|^2$ where $\xi=\nabla \wedge \vec{u}$ is the vorticity field of the flow.

The changes in the wave kinematics described in the preceding section are apparent in Fig. 3. As the wave packet propagates, the wavelength along the y direction decreases (that is, $k_y$ increases) while the wavelength along the vertical direction remains unchanged (that is, $k_z$ is unchanged). As a consequence, the direction of the phase lines tilts toward the vertical direction. Indeed, a vector perpendicular to the phase lines in this plane is $(k_y, k_z)$ ($k_x$ being constant in time), which tends toward the horizontal direction as $k_y \to +\infty$. Another feature is that the steepness of the isopycnals increases from Figs. 3(a) to 3(c), the isopycnals being locally vertical in Fig. 3(c). They have overturned in Fig. 3(d), implying that the wave packet is prone to breaking. Note that the breaking process, which may be defined as the generation of small dissipative scales, has not occurred yet: the isopycnals are still continuous lines in Fig. 3(d). Figure 3 also shows that the wave packet propagates beyond the theoretical location of the trapping plane [marked by a straight vertical line in Fig. 3(a)]. This point will be further discussed in the next section.

The generation of small dissipative scales occurs at $t \approx 30$, as displayed in Fig. 4(a). This event is manifested as a first peak in the entrophy curve (Fig. 6). It occurs at $t=30$ so that this time will be taken as the breaking time. A turbulent regime next sets in. The turbulent region broadens about the trapping plane, through horizontal intrusions of turbulent fluid within the shear flow [Figs. 4(b) and 5(c)], and intensifies [Figs. 4(c) and 5(d)], up to the time viscous dissipation at small scale overcomes turbulent production: at $t=58.8$, as seen in Fig. 6, the entrophy has reached a maximum value (over space and time) and decays afterward.

It should be noted that breaking does not always occur,
<figure><img src="image_5.png" /></figure>

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

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

<figure><img src="image_8.png" /><figcaption>FIG. 3. Contours of the total density field $-N^2z + \rho'$ in a vertical (y,z) plane. (a) $t=0$; the vertical line marks the theoretical location of the trapping plane; (b) $t=12$; (c) $t=24$; (d) $t=28.8$.</figcaption></figure>

even if there exists a trapping plane in the flow. In run 3 where the viscosity is 10 times larger than in run 1a, the phase lines tilt toward the vertical as the wave propagates toward that plane but the steepness of the isopycnals hardly

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

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

<figure><img src="image_11.png" /><figcaption>FIG. 4. Same as Fig. 3 at later times: (a) $t=30$; (b) $t=37.8$; (c) $t=58.8$.</figcaption></figure>

increases. The wave packet dissipates without breaking in the neighborhood of the trapping plane before $t=30$.

## B. Energetics

Insight into the energetics of the wave packet can be gained by linearizing the equations of motions about the shear flow. Let $(u' + \bar{U}, v', w')$ be the decomposition of the velocity field into a fluctuating part $\vec{u}'$ and a mean part $\bar{U}\vec{t}_x$, where the overline denotes a spatial average over $x$ and $z$. (Note that $v' = v$ and $w' = w$ but the primes will be retained for these two components, for consistency with $u'$.) As long as the amplitude of the wave packet remains small with respect to $U_0$, $\vec{u}'$ and $\bar{U}\vec{t}_x$ may be identified with the velocity field induced by the wave packet and of the horizontal shear flow, respectively. The linearized equations for the fluctuat-
<figure><img src="image_12.png" /><figcaption>FIG. 5. (Color) Contours of the entrophy in the same vertical (y,z) plane as in Figs. 3 and 4. (a) t=0; (b) t=24; (c) t=37.8; (d) t=58.8.</figcaption></figure>

ing kinetic energy components (⟨u'²⟩/2, ⟨v'²⟩/2, ⟨w'²⟩/2) and available potential energy [defined by ⟨ρ'²⟩/(2N²), within a good approximation²³] averaged over the numerical domain are (per unit mass, ignoring dissipative effects):

<figure><img src="image_13.png" /><figcaption>FIG. 6. Maximum values of the entrophy Z over the numerical domain versus time. Z is normalized by its initial maximum value.</figcaption></figure>

$$ \frac{1}{2} \frac{d\langle u'^2 \rangle}{dt} = \left\langle p' \frac{\partial u'}{\partial x} \right\rangle - \left\langle u' v' \frac{d\bar{U}}{dy} \right\rangle + f\langle u' v' \rangle, \quad (16) $$

$$ \frac{1}{2} \frac{d\langle v'^2 \rangle}{dt} = \left\langle p' \frac{\partial v'}{\partial y} \right\rangle - f\langle u' v' \rangle, \quad (17) $$

$$ \frac{1}{2} \frac{d\langle w'^2 \rangle}{dt} = \left\langle p' \frac{\partial w'}{\partial y} \right\rangle - \langle \rho' w' \rangle, \quad (18) $$

$$ \frac{1}{2N^2} \frac{d\langle \rho'^2 \rangle}{dt} = \langle \rho' w' \rangle. \quad (19) $$

$d/dt$ now denotes the Eulerian time derivative and $p'$ is the fluctuating pressure. The term $\langle u'v' \rangle$ is the horizontal Reynolds stress averaged over the numerical domain, $-\langle u'v' d\bar{U}/dy \rangle$ is responsible for energy exchange between the wave packet and the mean flow and $\langle \rho'w' \rangle$, the buoyancy flux, reversibly transfers vertical kinetic energy into available potential energy.

In order to track the location of the wave packet as time evolves, we plot in Fig. 7(a) the shear of the unperturbed mean flow $dU_b/dy$ at the y location, denoted $y_w$, the x and z averaged wave kinetic energy is maximum. The shear of the actual mean flow $d\bar{U}/dy$ at location $y_w$ is also displayed in the figure. $y_w$ itself is drawn as a function of time in Fig. 7(b). Temporal evolutions are displayed up to the time breaking occurs. Figure 7(a) shows that the wave packet enters the region of strong shear at $t \approx 10$. The shear of the unperturbed mean flow $dU_b/dy(y_w)$ strongly increases from this time on, thus the wave packet continues to propagate into the mean flow. From $t \approx 20$, $dU_b/dy(y_w)$ saturates at a value close to 0.8, just below its value at the trapping plane (equal to 0.819). This means that the propagation of the wave packet has nearly stopped in the neighborhood of this plane. This is confirmed by Fig. 7(b): $y_w$ is nearly constant from $t \approx 20$ and close to, but smaller than, $y_t$ (equal to 0.55). The
<figure><img src="image_14.png" /><figcaption>FIG. 7. (a) Horizontal shear of the unperturbed shear layer $dU_b/dy$ (solid line) and of the actual mean flow $d\bar{U}/dy$ (dashed line) at the y location, denoted $y_w$, the x and z averaged wave kinetic energy is maximum versus time. The value of $dU_b/dy$ at the theoretical trapping plane is indicated by a dashed-dotted line. (b) Temporal evolution of $y_w$. The trapping plane position $y_t$ is indicated with a dashed-dotted line. (c) Horizontal Reynolds stress of the wave packet times the Coriolis frequency $f(u'v')$ (plotted with a minus sign to yield a positive quantity) (dashed line) and rate of change of the energy of the wave packet due to energy exchange with the mean flow $-\langle u'v' d\bar{U}/dy \rangle$ (solid line); both quantities are averaged over the numerical domain and plotted versus time.</figcaption></figure>

shear of the actual mean flow $d\bar{U}/dy(y_w)$ follows that of the undisturbed mean flow up to $t \approx 24$. From this time on, $d\bar{U}/dy(y_w)$ decays. This decay may be due to finite amplitude effects (such as a self-induced mean flow) resulting from the growth of the wave packet energy. The local decay of the actual mean shear results in a decay of the intrinsic frequency of the wave packet [by Eq. (9)], so that the packet can propagate further down. This is visible in Fig. 7(b): the decrease of $d\bar{U}/dy(y_w)$ at $t=25$ makes $y_w$ reach the trapping plane. The latter plane is eventually crossed when $d\bar{U}/dy(y_w)$ further decays at $t=30$. This effect was also noted in Sec. III A.

The domain averaged horizontal Reynolds stress times the Coriolis frequency $f(u'v')$ is plotted in Fig. 7(c) (with a minus sign) along with the energy exchange term

$\langle -u'v' d\bar{U}/dy \rangle$. An estimate of the former term may be obtained, using WKB approximation, from expressions (5) for $u'$ and $v'$ without the Gaussian modulation (we recall that $u'$ and $v'$ are normalized by $w'$). This yields

$$ \langle u'v' \rangle = \frac{s^2}{2} k_x \frac{k_y}{(k_h^2)^2} (\Omega^2 - f^2) \sim s^2 \frac{k_x}{k_y^3} (\Omega^2 - f^2) \quad \text{as } k_y \to +\infty. \quad (20) $$

$\langle u'v' \rangle$ is thus negative in the present case, $k_x$ being negative (also, $k_y>0$ and $\Omega>f$). Since $d\bar{U}/dy$ varies much slower than $u'v'$ in the neighborhood of the trapping plane, $-\langle u'v' d\bar{U}/dy \rangle$ may be approximated by $\langle -u'v' \rangle d\bar{U}/dy$. Because $d\bar{U}/dy$ is positive, the latter term is positive: the kinetic energy of the wave field is fed by the mean shear
<figure><img src="image_15.png" /><figcaption>FIG. 8. Temporal evolution of the volume-averaged fluctuating kinetic energy (per unit mass) $\langle u'^2 \rangle/2$ (bottom solid line), $\langle v'^2 \rangle/2$ (dashed line), and $\langle w'^2 \rangle/2$ (dotted line); the volume-averaged potential energy $\langle p'^2/(2N^2) \rangle$ is also shown (dashed-dotted line). The upper solid line displays the volume averaged kinetic energy of the mean shear $\langle \bar{U}^2 \rangle/2$, divided by a factor 60 to fit the range of the vertical coordinate of the fluctuating components.</figcaption></figure>

flow. This is attested in Fig. 7(c): $-f\langle u'v' \rangle$ is positive and decays with time (relative to $\langle w'^2 \rangle$), in agreement with the estimate (20) when $k_y$ increases (since $k_x$ is constant and $s$ remains of order 1). As for the energy exchange term, Fig. 7(c) shows that it nearly vanishes at initial time because the wave is located out of the shear region; it displays a (weak) positive maximum value at $t \approx 10$, when the wave packet has penetrated the strong shear region and stays about constant up to $t \approx 25$, implying that the decay of $-f\langle u'v' \rangle$ is partly made up for by the increase of the shear. The estimate (20) also shows that the extraction of energy from the mean flow is reduced by the Coriolis force, owing to the $\Omega^2 - f^2$ factor.

The components of the fluctuating kinetic energy and of the energy of the mean flow are plotted in Fig. 8 along with the available potential energy of the wave (the latter curve is plotted up to $t=30$ since the notion of wave becomes meaningless when breaking occurs). Equation (16) shows that the energy extracted from the mean flow is fed into the component of the wave energy parallel to the mean flow, as physically expected. This energy is transferred into the $\langle v'^2 \rangle$ component by the Coriolis term $f\langle u'v' \rangle$ and into the $\langle w'^2 \rangle$ (and $\langle v'^2 \rangle$ as well) components by the pressure-velocity correlation term. As a result, $\langle u'^2 \rangle$ constantly decreases during the propagation of the wave toward the barrier. By contrast, the feeding of the $\langle v'^2 \rangle$ component through the Coriolis term results in a nearly constant level of this component up to $t \approx 20$, despite the progressive orientation of the wave motions toward a vertical $(x,z)$ plane. As expected, most of the energy gained by the wave packet is fed into vertical motions: the vertical kinetic energy component increases by a factor of 10 between $t=0$ and $t=26$. Figure 8 also shows that the total wave energy (kinetic energy plus available potential energy) increases by a factor of 1.8 between $t=0$ and $t=30$. This increase cannot be attributed to the accumulation

of the energy of the wave packet near the trapping plane. Indeed, the accumulation will lead to a local increase of the wave energy density but will leave unchanged the volume averaged energy.

From the breaking time $t=30$, the dynamics of the wave packet change dramatically (Fig. 8): the $\langle u'^2 \rangle$ component suddenly increases and, by the Coriolis force, the $\langle v'^2 \rangle$ component increases as well. The $\langle w'^2 \rangle$ component displays the same behavior, with some delay. Figure 8 shows that the energy of the mean shear flow appreciably decays during this turbulent stage: the small scale turbulence is thus sustained by the mean shear. We found indeed that the correlated growth of $\langle u'^2 \rangle$ and $\langle v'^2 \rangle$ gives a strong increase in $-\langle u'v' d\bar{U}/dy \rangle$, which in turn feeds the small scale turbulence $(-\langle u'v' d\bar{U}/dy \rangle)$ increases by a factor of 100 between its minimum value at $t \approx 24$ and its maximum value at $t \approx 50$. Note that the energy exchange term can no longer be interpreted as an energy transfer from the mean shear flow to the wave (since the notion of wave has become meaningless) but, because of the clear scale separation that exists during this regime, from a large scale shear flow to small scale turbulence. Such a regime is triggered by the breaking process.

When $U_0$ is increased, all other parameters being unchanged, $\langle d\bar{U}/dy \rangle$ increases as well so that the location of the trapping plane gets closer to the initial position of the wave packet (see the value of $y_t$ in the table for runs 6 to 8). As a consequence, breaking occurs earlier than in run 1a.

### C. Instability and breaking

The maxima of the components of the fluctuating vorticity $\bar{\zeta}' = \nabla \wedge \vec{u}'$ are plotted versus time in Fig. 9. For $t \le 30$, the main result is the growth at a constant rate of the streamwise component $\zeta_x'$ while the other components remain at a nearly constant level. The growth of $\zeta_x'$ is easily accounted for by the tilting of the fluid motion toward a vertical plane while the wave number $k_y$ increases, as the packet propagates toward the barrier: the dominant velocity gradient is then $\partial w'/\partial y$, which contributes to (and controls) $\zeta_x'$. $\zeta_x'$ in turn controls the dissipation rate of kinetic energy for $t \le 30$: Fig. 6 shows indeed that the maximum value of $Z$ behaves identically to that of $\zeta_x'$ up to the breaking time.

As shown from Figs. 3(d) and 4(a), the instability through which the wave breaks is a convective instability (the term “buoyancy induced” instability, employed by Schowalter et al.²⁴ would actually be more appropriate). This instability results in the formation of vortices of alternate sign, with axis aligned along the $x$ direction. These vortices are visualized in Fig. 10 by a cross section of the streamwise vorticity in the $(y,z)$ plane. An noted in the Introduction, such a convective instability is known to develop from a large amplitude wave of low frequency but the result is also valid when the wave is of high frequency (see Fritts et al.¹³ for numerical evidence and Benielli and Sommeria²⁵ for experimental evidence). The breaking process eventually leads to a turbulent flow with all vorticity components having the same amplitude (Fig. 9).

The turbulent regime does not result from a Kelvin-
<figure><img src="image_16.png" /><figcaption>FIG. 9. Maximum values over the numerical domain $D$ of the fluctuating vorticity components: (a) $\max_D \zeta_x'$; (b) $\max_D \zeta_y'$; (c) $\max_D \zeta_z'$.</figcaption></figure>

Helmholtz instability of the shear layer. The flow at $t=0$ is solely made of the wave packet and of the shear layer so that no background noise has perturbed the shear layer when the wave packet reaches its highest vorticity region. Of course,

<figure><img src="image_17.png" /><figcaption>FIG. 10. (Color) Contours of the streamwise vorticity component $\zeta_x = (\nabla \times \vec{u}) \cdot \vec{v}_x$ at $t=31.8$ in the same $(y,z)$ plane as in Figs. 3 and 4.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 11. (Color) Contours of the vertical vorticity component $\zeta_z = (\nabla \times \vec{u}) \cdot \vec{v}_z$ in a horizontal $(x,y)$ plane. (a) $t=36$; (b) $t=40.8$.</figcaption></figure>
the wave packet may play the role of a perturbation of finite amplitude for the shear layer. Two factors actually stabilize the layer against this two-dimensional instability. The first one is that the streamwise extent of the shear layer has been chosen to be a little lower (in view of the constraints listed in Sec. II C) than the most amplified wavelength, in order to lower the growth rate of the Kelvin–Helmholtz instability; this yields $L_x=12$ instead of $L_x \approx 14.5$ (for a hyperbolic tangent velocity profile, see Michalke$^{26}$). The second factor, which is much more efficient at stabilizing the shear flow against a Kelvin–Helmholtz instability, stems from the small scale turbulence itself: the velocity profile of the shear layer is broadened by the turbulence through momentum transport, as we show in the next section. The layer thickness $\delta(t)$ increases by a factor of 2 so that the ratio $L_x/\delta$ is greatly reduced.

By contrast, the shear layer is subjected to a (three-dimensional) inertial instability, which is triggered by the small scale turbulence stemming from the breaking of the wave packet. This instability next enhances this small scale turbulence. Indeed, a sufficient condition for a plane parallel shear flow to become locally unstable through an inertial instability is$^{27}$

$$f\left(f - \frac{dU_b}{dy}\right) < 0, \qquad (21)$$

which is fully satisfied with the present set of parameters ($f=0.23$ and $dU_b/dy=0.8$ in the neighborhood of the trapping plane). Preliminary calculations with an inertially stable shear layer (obtained by reversing the sign of $dU_b/dy$, and of $k_x$ as well to keep $-k_x dU_b/dy > 0$) show that the breaking process is similar but the resulting small scale turbulence quickly decays because it is no longer sustained by the mean shear flow.$^{28}$

The stability of the shear layer with respect to the two-dimensional Kelvin–Helmholtz instability is attested in Fig. 11, where contours of the vertical vorticity component $\zeta_z = (\nabla \times \vec{u}) \cdot \vec{t}_z$ are displayed at the beginning of the breaking process in a horizontal plane. We found that these contours are identical to those of $\zeta'_z$ at the same times, implying that the vertical vorticity is dominated by its fluctuating part. Very small structures of negative and positive sign are visible. A Kelvin–Helmholtz instability would have led to a large scale structure along the streamwise direction (Corcos and Sherman$^{29}$). The small scale structures belong to the wave field because they propagate along the streamwise direction when the vertical plane is moved vertically. Their small scale along the y direction is a manifestation of the decrease of the wave vector component along this direction as the wave propagates toward the barrier. The scale of a structure of a given sign along the y direction, $l_h$ say, may be approximately predicted by stating that $\zeta'_x$ is mostly contributed by $\partial w'/\partial y$. Let us assume that $w'$ is of the form $W \cos(k_y y - Nt)$ just before breaking. Thus, $l_h \approx \pi/k_y$. The breaking event being local, maximum values of $\zeta'_x$ and $w'$ over the numerical domain $D$ may be used to estimate $l_h$ so that $l_h \approx \pi \max_D w'/\max_D \zeta'_x$. From Fig. 9(a), $\max_D \zeta'_x \approx 11$ at $t=40$. We computed the maximum values over the numeri-

<figure><img src="image_19.png" /><figcaption>FIG. 12. Distribution of the passive scalar averaged over x and z as a function of y for five successive times: solid lines are used for $t=0$ and $t=84$, and dashed-dotted lines for intermediate times ($t=42, 54, 66$).</figcaption></figure>

cal domain of $w'$ (not shown) and found that $\max_D w' \approx 0.32$ at the same time. This yields $l_h \approx 0.09$, which is in good agreement with the value obtained from Fig. 11, equal to $\approx 0.12$.

## IV. TRANSPORT PROPERTIES

A passive scalar located on the opposite side of the trapping plane with respect to the wave packet is added to the initial condition in order to quantify the transport triggered by the breaking process. The scalar field $S$ has a Gaussian distribution about a vertical plane $y=y_s$,

$$S(x, y, z, t=0) = A_s \exp^{-(y-y_s/\sigma_0)^2}. \qquad (22)$$

The value of $y_s$ was varied from 0.65 to 1.65 without any apparent change in the behavior of the scalar field as we shall show it. We thus describe the results for the case $y_s=1$, that is when the Gaussian distribution is centered about the same value as the vorticity of the undisturbed shear layer. $\sigma_0$ is equal to 1.3, in units of the half-thickness of the shear layer $a$.

The temporal evolution of the scalar distribution averaged over x and z is plotted in Fig. 12. The distribution remains nearly Gaussian as time elapses. Diffusion starts at a different time on the left-hand side than on the right-hand side so that, at a given time, each side of the distribution was fitted by a Gaussian function with a different standard deviation. These standard deviations are referred to as $\sigma_l$ and $\sigma_r$, respectively, and their average value as $\sigma$.

The variance $\sigma^2$ is plotted as a function of time in Fig. 13(a), while the right and left variances $\sigma_r^2$ and $\sigma_l^2$ are displayed in Figs. 13(b) and 13(c), respectively. $\sigma^2$ remains constant and equal to $\sigma_0^2$ as long as breaking has not occurred. When breaking occurs, $\sigma^2$ starts departing from its initial value and displays a nearly linear evolution from $t=35$. A change in the slope occurs at $t \approx 55$, which coincides with the decay of the small scale turbulence. An identical behavior is found for run 2, in which $Pe=1$ and in run 4a, in
<figure><img src="image_20.png" /><figcaption>FIG. 13. (a) Variance of the passive scalar distribution $\sigma^2$ as a function of time. *, run 1a; o, run 2; +, run 4a. (b) Variance of the right-hand side of the passive scalar distribution, symmetrized about the instantaneous maximum value and fitted by a Gaussian function, as a function of time. (c) Same as (b), for the left-hand side of the distribution.</figcaption></figure>

which $y_s=1.65$. In the latter run, the change of slope is delayed because of the more distant position of $y_s$ from the breaking region. Overall, $\sigma^2$ grows up to a value close to 5 and $\sigma_l^2$ up to 5.5, which is nearly three times larger than the shear layer thickness: this means that transport occurs across the barrier. The linear growth of $\sigma^2$ attests that the transport by the small scale turbulence amounts to a diffusive process. The diffusion coefficient is defined by $\kappa_t = 0.25 d\sigma^2/dt$ and from Fig. 13(a), a value of 0.026 is obtained for $\kappa_t$ for $30 \le t \le 55$. In order to compute a Peclet number associated with the transport process, we normalize $\kappa_t$ by a characteristic length scale and a characteristic velocity scale associated with the breaking wave packet along the direction of transport $y$. We choose $l_h$ for the length scale (equal to 0.12 at $t=40$) and the maximum value of the $v'$ component for the velocity scale (equal to 0.4 at the same time). The Peclet number $l_h v/\kappa_t$ thus obtained is equal to 0.5, which is of order one: the transport induced by the small scale turbulence is indeed of turbulent nature. This conclusion is confirmed by the fact that this turbulent Peclet number remains unchanged when the molecular Peclet number is decreased by a factor of 10. (We also found that $\kappa_t$ increases linearly with $U_0$ but we did not compute the dependence of the Peclet number on $U_0$.)

The small scale turbulence also transports momentum. This is attested in Fig. 14 from the temporal evolution of the mean velocity profile of the shear layer $\bar{U}(y,t)$. $\bar{U}(y,t)$ remains unchanged as long as breaking has not occurred and smoothes out for $t>30$ (significant fluctuations about the mean profile are however observed at $t=48$). Like for the passive scalar, the diffusive nature of the transport implies that the velocity profile should evolve self-similarly, being of the form $f(y/\delta(t))$. If the velocity profile were of the error function type, $\delta^2$ should be equal to $\pi v_t + C$, where $v_t$ is a turbulent viscosity and $C$ a constant value. A hyperbolic tangent function being very close to an error function, we shall assume that this scaling holds in the present case. $\delta^2$ is therefore plotted as a function of time in Fig. 14(b). It remains constant up to $t \approx 30$ and starts increasing from this time. The growth is not linear but can reasonably be approximated by a linear law, displayed with a dotted line in the figure. This linear law yields a value of 0.025 for $v_t$. Hence, the turbulent Prandtl number $v_t/\kappa_t$ associated with the transport by the small scale turbulence is nearly equal to 1.

## V. CONCLUSION

We have shown that an inertia-gravity wave packet may break in the neighborhood of a barotropic dynamical barrier and trigger a turbulent diffusive transport of mass and momentum across the barrier. Three-dimensional numerical simulations of the Boussinesq equations in a rotating reference frame were performed for this purpose. The dynamical barrier consists of a horizontally sheared horizontal shear layer, toward which an inertia-gravity wave packet propagates. The basic mechanism of the interaction is the increase of the intrinsic frequency of the wave packet $\Omega$ by Doppler shifting, which leads to the formation of a trapping plane at the location where $\Omega=N$. As this plane is approached, en-
<figure><img src="image_21.png" /><figcaption>FIG. 14. (a) Mean velocity profile averaged over $x$ and $z$, $\bar{U}(y,t)$, at successive times. The curves are plotted every 12 time units, from $t=0$ (stepest solid line) to $t=84$ (smoothest solid line), a dashed-dotted line being used for intermediate times. (b) Squared thickness $\delta^2(t)$ of the mean velocity profile, fitted by the function $\tanh(y-y_1)/\delta(t)$, as a function of time. The slope of the dotted line, equal to 0.08, is the mean value of $d(\delta^2)/d(t)$ for $33 \le t \le 57$.</figcaption></figure>

ergy is transferred from the mean flow to the wave packet, the wavelength along the direction of inhomogeneity decreases and the phase surfaces tilt toward a vertical plane; the wave steepness increases so that breaking may occur. The breaking process triggers an inertial instability of the shear layer (but no Kelvin-Helmholtz instability develops). As a result, a small scale turbulent regime sets in, which the mean flow sustains. Enhanced transport properties are found during this regime, with a turbulent Prandtl number of order one. Preliminary computations with an inertially stable shear layer show that the breaking process remains unchanged. But the small scale turbulence stemming from breaking is no longer sustained by the mean flow and, therefore, it quickly decays. The resulting transport properties are weakened and are currently quantified.²⁸

Though this study was motivated by geophysical considerations, the simple configuration we have addressed makes it premature to draw any implication to realistic cases. We still note that gravity wave frequency spectra with peaks at a frequency close to $N$ have been observed in the ocean,³⁰ a phenomenon that may be accounted for by the interaction of the waves with a horizontally sheared current [D'Asaro (private communication)]. As well, Jacobitz and Sarkar³¹ found that turbulence production in a stably stratified fluid is strongly increased when the flow involves a horizontally sheared horizontal mean flow, a numerical result which may be explained by the occurrence of wave breaking in the neighborhood of trapping locations.

The next step is to investigate the interaction of inertia-gravity waves with a baroclinic barrier, and this is the purpose of current research.

## ACKNOWLEDGMENTS

This work benefitted from discussions with Jan-Bert Flor, Claude Sidi, and Neil Edwards. Support was provided by the PNCA research program of CNRS/INSU. Computer time was allocated by the Computer Center of CNRS (ID-RIS) under Contracts Nos. 990580 and 000580.

<footnote n="1">¹M.E. McIntyre, “The stratospheric polar vortex and sub-vortex: fluid dynamics and mid-latitude ozone loss,” Philos. Trans. R. Soc. London **352**, 227 (1995).</footnote>

<footnote n="2">²Y.A. Ivanov and Y.G. Morozov, “Deformation of internal gravity waves by a stream with horizontal shear,” Oceanology **14**, 457 (1974).</footnote>

<footnote n="3">³D.J. Olbers, “The propagation of internal waves in a geostrophic current,” J. Phys. Oceanogr. **11**, 1224 (1981).</footnote>

<footnote n="4">⁴A.Y. Basovich and L.S. Tsimring, “Internal waves in a horizontally inhomogeneous flow,” J. Fluid Mech. **142**, 233 (1984).</footnote>

<footnote n="5">⁵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 (1985).</footnote>

<footnote n="6">⁶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 (1993).</footnote>

<footnote n="7">⁷F.P. Bretherton, “The propagation of groups of internal gravity waves in a shear flow,” Q. J. R. Meteorol. Soc. **92**, 466 (1966).</footnote>

<footnote n="8">⁸J.R. Booker and F.P. Bretherton, “The critical layer for internal gravity waves in a shear flow,” J. Fluid Mech. **27**, 513 (1967).</footnote>

<footnote n="9">⁹J.B. Baik, H.S. Hwang, and H.Y. Chun, “Transient critical-level effect for internal gravity waves in a stably stratified shear flow with thermal forcing,” Phys. Fluids **11**, 238 (1999).</footnote>

<footnote n="10">¹⁰M.E. McIntyre, “On global-scale atmospheric circulations,” in *Perspectives in Fluid Dynamics*, edited by G.K. Batchelor, H.K. Mofatt, and M.G. Worster (Cambridge University Press, Cambridge, UK, 2000), p. 557.</footnote>

<footnote n="11">¹¹K.B. Winters and E.A. D'Asaro, “Three-dimensional wave instability near a critical level,” J. Fluid Mech. **272**, 255 (1994).</footnote>

<footnote n="12">¹²A. Dörnbrack, “Turbulent mixing by breaking gravity waves,” J. Fluid Mech. **375**, 113 (1998).</footnote>

<footnote n="13">¹³D.C. Fritts, S. Arendt, and O. Andreassen, “Vorticity dynamics in a breaking internal gravity wave. Part 2. Vortex interactions and transition to turbulence,” J. Fluid Mech. **367**, 65 (1998).</footnote>

<footnote n="14">¹⁴M.-P. Lelong and T.J. Dunkerton, “Inertia-gravity wave breaking in three dimensions. Part II: Convectively unstable waves,” J. Atmos. Sci. **55**, 15 (1998); **55**, 2489 (1998).</footnote>

<footnote n="15">¹⁵S.A. Thorpe, “An experimental study of critical layers,” J. Fluid Mech. **103**, 321 (1981).</footnote>

<footnote n="16">¹⁶C.G. Koop and B. McGee, “Measurements of internal gravity waves in a continuously stratified shear flow,” J. Fluid Mech. **172**, 453 (1986).</footnote>

<footnote n="17">¹⁷C. Staquet and J. Sommeria, “Internal gravity waves, from instabilities to turbulence,” Annu. Rev. Fluid Mech. **34**, 539 (2002).</footnote>
<footnote n="18">¹⁸B.R. Sutherland, “Propagation and reflection of large amplitude internal gravity waves,” Phys. Fluids ¹¹, 1081 (1999).</footnote>

<footnote n="19">¹⁹B.R. Sutherland, “Internal wave reflection in uniform shear,” Q. J. R. Meteorol. Soc. ¹²⁶, 3255 (2000).</footnote>

<footnote n="20">²⁰A.E. Gill, *Atmosphere-Ocean Dynamics* (Cambridge University Press, Cambridge, UK, 1985).</footnote>

<footnote n="21">²¹J. Lighthill, *Waves in Fluids* (Cambridge University Press, Cambridge, UK, 1978).</footnote>

<footnote n="22">²²S.A. Orszag, “Numerical simulation of incompressible flows within simple boundaries. I. Galerkin (spectral) representation,” Stud. Appl. Math. ⁵⁰, 293 (1971).</footnote>

<footnote n="23">²³D. Holliday and M.E. McIntyre, “On potential energy density in an incompressible stratified fluid,” J. Fluid Mech. ¹⁰⁷, 221 (1981).</footnote>

<footnote n="24">²⁴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. ²⁸¹, 247 (1994).</footnote>

<footnote n="25">²⁵D. Benielli and J. Sommeria, “Excitation and breaking of internal gravity waves by parametric instability,” J. Fluid Mech. ³⁷⁴, 117 (1998).</footnote>

<footnote n="26">²⁶A. Michalke, “On the inviscid instability of the hyperbolic-tangent velocity profile,” J. Fluid Mech. ¹⁹, 543 (1964).</footnote>

<footnote n="27">²⁷T.J. Pedley, “On the stability of viscous flow in a rapidly rotating pipe,” J. Fluid Mech. ³⁵, 97 (1969).</footnote>

<footnote n="28">²⁸N.R. Edwards and C. Staquet, “Interactions of an inertia-gravity wave packet with a barotropic shear layer” (in preparation).</footnote>

<footnote n="29">²⁹G.M. Corcos and F.S. Sherman, “The mixing layer: deterministic models of a turbulent flow. Part 1: Introduction and the two-dimensional flow,” J. Fluid Mech. ¹³⁹, 29 (1984).</footnote>

<footnote n="30">³⁰E.A. D'Asaro and R.C. Lien, “The wave-turbulence transition in stratified flows,” J. Phys. Oceanogr. ³⁰, 1669 (2000).</footnote>

<footnote n="31">³¹F.G. Jacobitz and S. Sarkar, “The effect of nonvertical shear on turbulence in a stably stratified medium,” Phys. Fluids ¹⁰, 1158 (1998).</footnote>