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

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 V of the wave packet by Doppler shifting may result in the breaking of the wave packet, in the neighborhood of a trapping plane where V 5 N ~ N is the Brunt–Va¨isa¨la¨ frequency of the ﬂuid 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. @

On transport across a barotropic shear flow by breaking inertia-gravity waves C. Staquet a) and G. Huerre 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 ⍀ of the wave packet by Doppler shifting may result in the breaking of the wave packet, in the neighborhood of a trapping plane where ⍀ϭ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 1 ͒ 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 1 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 2 and a Wentzel-Kramers-Brillouin ͑WKB͒ approximation. [3][4][5][6] 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 7 and Booker and Bretherton; 8 see Baik et al. 9 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 ⍀ ͑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 7 first showed that, for a wave packet, the vertical component of the wave vector increases as Ϫ1 , where is the distance of the wave packet from the critical level ͑located at ϭ0͒ along the vertical direction and the vertical component of the group velocity decreases as 2 ͓Fig. 1͑a͔͒. However, the WKB approximation also predicts that the wave induced horizontal velocity uЈ goes to infinity ͑uЈϳ Ϫ1/2 as →0͒ and the critical level is never reached ͑the time required to reach that level scales like Ϫ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 8 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, 10 the horizontal phase speed of the wave packet c x decays faster than the wave induced horizontal velocity uЈ as t→ϩϱ 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 a͒ Author to whom correspondence should be addressed. Electronic mail: chantal.staquet@hmg.inpg.fr may occur before the critical level is reached. Using threedimensional numerical simulations, Winters and D'Asaro 11 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 12 and Fritts et al.: 13 breaking first occurs through a convective instability. This result is consistent with the numerical finding by Lelong and Dunkerton 14 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 15 and Koop and McGee 16 ͑see Staquet and Sommeria 17 for a re-view͒. 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 ⍀ 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 ⍀р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 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 goes to zero, the component of the wave vector along that direction increases again without bound, as Ϫ1/2 , and the group velocity decreases as 3/2 ͑see Ref. 3͒. Like for the critical level situation, the wave energy increases without bound (ϳ Ϫ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. 2 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., 5 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. 5 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. 3 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 18,19 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. The shear is taken constant for simplicity. In ͑a͒, the interaction process is two dimensional and occurs in the vertical (k ជ ,g ជ ) plane. A critical level is located at zϭz c , at which the intrinsic frequency of the wave ⍀ vanishes. As the critical layer is approached, c gz →0, c gx →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; k ជ h and c ជ g h , 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 ⍀ is equal to N. As the trapping plane is approached, c gy →0, c gx →U b , k x remains constant while k y increases without limit.

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. 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 ( i ជ x , i ជ y , i ជ z ). Let 0 ϩ lin (z) be the density field of the hydrostatic equilibrium state, where lin Ӷ 0 for any z and 0 is a constant reference density. The linear density profile lin (z) sets the Brunt-Väisälä frequency 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, This flow is in geostrophic equilibrium since the Coriolis force Ϫ f ı ជ z ٙ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, 20 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

͑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 k ជ ϭ(k x ,k y ,k z ) is introduced as the initial condition. The wave vector k ជ sets the intrinsic frequency ⍀ of the wave packet via the dispersion relation where k h ϭͱk x 2 ϩk y 2 is the horizontal wave number and k ϭ͉k ជ ͉. Relation ͑4͒ implies that f р⍀рN. The group velocity of the wave packet in the rotating reference frame ͑that is, the absolute group velocity͒ is c ជ g ϩU ជ b , where c ជ g ϭ‫ץ‬⍀/‫ץ‬k ជ is the intrinsic group velocity ͑e.g., Lighthill, 21 p. 327͒. In the present case where the shear flow is directed along the horizontal x direction, the group velocity is expressed by (‫ץ‬⍀/‫ץ‬k x ϩU b ,‫ץ‬⍀/‫ץ‬k y ,‫ץ‬⍀/‫ץ‬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 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, 3 Gill,20 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 0 in the rotating frame of reference ͑the absolute frequency͒ differs from its intrinsic frequency. 0 is defined by the Doppler relation 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 introduction to the WKB approximation for internal gravity waves is provided by Olbers. 3 ͒ Since U b depends only upon y, the k x and k z components of the wave packet do not change along a ray. By contrast, k y changes according to where d/dtϭ‫ץ/ץ‬tϩ(c ជ g ϩU ជ b )•ٌ 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 0 stays constant along a ray (d 0 /dtϭ0). Hence, using the Doppler relation ͑7͒, the intrinsic frequency changes along a ray according to 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 ⍀, the wave energy E is not conserved along a ray ͑it would be if the fluid were at rest or uniformly translating͒ but the wave action AϭE/⍀ is conserved instead of the energy 7 for a nondissipative and nondiffusive medium. As noted in the Introduction, the trapping location is a vertical yϭy t plane obtained by setting ⍀ to N in the Doppler relation ͑7͒. Using 0 ϭ⍀(y 0 )ϩk x U b (y 0 ), one gets an expression for the location of the trapping plane as a function of the initial parameters of the wave 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 rotating reference frame introduced in the preceding section. Let u ជ ϭ(u,v,w) be the components of the velocity field in this reference frame and (x,y,z,t), the deviation of the density field from the hydrostatic density field 0 ϩ lin (z). As is customary, we scale as an acceleration using the variable Јϭ (g/ 0 ) in order for N to come into play in the equations of motion. These equations thus are ͑written per unit mass͒ In Eq. ͑12͒, the term P stands for the dynamical pressure divided by 0 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 0 plane propagates in a horizontal shear layer centered about the yϭy 1 plane. Periodic boundary 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: 22 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 ٌ͓ 2 u ជ , (/Pr) ٌ 2 Ј and (/Pe) ٌ 2 S, respectively͔ involve a constant molecular viscosity , 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 ͑k ជ and y 0 for the wave; U 0 and y 1 for the barrier; f , , 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 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 k ជ ϭ(Ϫ/3,/3,4/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, ϭ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.

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 2 zϩЈ. 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 →ϩϱ. 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 Ӎ30, as displayed in Fig. 4͑a͒. This event is manifested as a first peak in the enstrophy 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 enstrophy has reached a maximum value ͑over space and time͒ and decays afterward.
It should be noted that breaking does not always occur, 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, 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. 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 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ЈϩŪ ,vЈ,wЈ) be the decomposition of the velocity field into a fluctuating part u ជ Ј and a mean part Ū ı ជ 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 , u ជ Ј and Ū ı ជ 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- ing kinetic energy components ͑͗uЈ 2 ͘/2, ͗vЈ 2 ͘/2, ͗wЈ 2 ͘/2͒ and available potential energy ͓defined by ͗Ј 2 ͘/(2N 2 ), within a good approximation 23 ͔ averaged over the numerical domain are ͑per unit mass, ignoring dissipative effects͒: d/dt now denotes the Eulerian time derivative and pЈ is the fluctuating pressure. The term ͗uЈvЈ͘ is the horizontal Reynolds stress averaged over the numerical domain, Ϫ͗uЈvЈdŪ /dy͘ is responsible for energy exchange between the wave packet and the mean flow and ͗ЈwЈ͘, 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Ū /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Ӎ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Ӎ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 Ӎ20 and close to, but smaller than, y t ͑equal to 0.55͒. The shear of the actual mean flow dŪ /dy(y w ) follows that of the undisturbed mean flow up to tӍ24. From this time on, dŪ /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Ū /dy(y w ) at tϭ25 makes y w reach the trapping plane. The latter plane is eventually crossed when dŪ /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͒ ͗uЈvЈ͘ is thus negative in the present case, k x being negative ͑also, k y Ͼ0 and ⍀Ͼ f ͒. Since dŪ /dy varies much slower than uЈvЈ in the neighborhood of the trapping plane, Ϫ͗uЈvЈ dŪ /dy͘ may be approximated by ͗ϪuЈvЈ͘dŪ /dy.
Because dŪ /dy is positive, the latter term is positive: the kinetic energy of the wave field is fed by the mean shear flow. This is attested in Fig. 7͑c͒: Ϫ f ͗uЈvЈ͘ is positive and decays with time ͑relative to ͗wЈ 2 ͘), 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Ӎ10, when the wave packet has penetrated the strong shear region and stays about constant up to tӍ25, implying that the decay of Ϫ f ͗uЈvЈ͘ 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 ⍀ 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 ͗vЈ 2 ͘ component by the Coriolis term f ͗uЈvЈ͘ and into the ͗wЈ 2 ͘ ͑and ͗vЈ 2 ͘ as well͒ components by the pressure-velocity correlation term. As a result, ͗uЈ 2 ͘ constantly decreases during the propagation of the wave toward the barrier. By contrast, the feeding of the ͗vЈ 2 ͘ component through the Coriolis term results in a nearly constant level of this component up to t Ӎ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 ͗uЈ 2 ͘ component suddenly increases and, by the Coriolis force, the ͗vЈ 2 ͘ component increases as well. The ͗wЈ 2 ͘ 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 ͗uЈ 2 ͘ and ͗vЈ 2 ͘ gives a strong increase in Ϫ͗uЈvЈ dŪ /dy͘, which in turn feeds the small scale turbulence (Ϫ͗uЈvЈ dŪ /dy͘ increases by a factor of 100 between its minimum value at tӍ24 and its maximum value at t Ӎ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, ͗dŪ /dy͘ 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 ជ Јϭٌٙu ជ Ј are plotted versus time in Fig. 9. For tр30, the main result is the growth at a constant rate of the streamwise component x Ј while the other components remain at a nearly constant level. The growth of 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 ‫ץ‬wЈ/‫ץ‬y, which contributes to ͑and controls͒ x Ј . x Ј in turn controls the dissipation rate of kinetic energy for t р30: Fig. 6 shows indeed that the maximum value of Z behaves identically to that of 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. 24 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. 13 for numerical evidence and Benielli and Sommeria 25 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- 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 Ӎ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 ␦(t) increases by a factor of 2 so that the ratio L x /␦ 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 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 twodimensional Kelvin-Helmholtz instability is attested in Fig.  11, where contours of the vertical vorticity component z ϭ(ٌٙu ជ)•ı ជ z are displayed at the beginning of the breaking process in a horizontal plane. We found that these contours are identical to those of 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 x Ј is mostly contributed by ‫ץ‬wЈ/‫ץ‬y. Let us assume that wЈ is of the form W cos(k y yϪNt) just before breaking. Thus, l h Ӎ/k y . The breaking event being local, maximum values of x Ј and wЈ over the numerical domain D may be used to estimate l h so that l h Ӎ max D wЈ/max D x Ј . From Fig. 9͑a͒, max D x ЈӍ11 at tϭ40. We computed the maximum values over the numeri-cal domain of wЈ ͑not shown͒ and found that max D wЈ Ӎ0.32 at the same time. This yields l h Ӎ0.09, which is in good agreement with the value obtained from Fig. 11, equal to Ӎ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 / 0 ) 2 . ͑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. 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 l and r , respectively, and their average value as .
The variance 2 is plotted as a function of time in Fig.  13͑a͒, while the right and left variances r 2 and l 2 are displayed in Figs. 13͑b͒ and 13͑c͒, respectively. 2 remains constant and equal to 0 2 as long as breaking has not occurred. When breaking occurs, 2 starts departing from its initial value and displays a nearly linear evolution from t Ӎ35. A change in the slope occurs at tӍ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 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, 2 grows up to a value close to 5 and 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 2 attests that the transport by the small scale turbulence amounts to a diffusive process. The diffusion coefficient is defined by t ϭ0.25 d 2 /dt and from Fig. 13͑a͒, a value of 0.026 is obtained for t for 30 рtр55. In order to compute a Peclet number associated with the transport process, we normalize 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/ 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 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 Ū (y,t). Ū (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/␦(t)). If the velocity profile were of the error function type, ␦ 2 should be equal to t tϩC, where 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. ␦ 2 is therefore plotted as a function of time in Fig. 14͑b͒. It remains constant up to tӍ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 t . Hence, the turbulent Prandtl number t / 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 ⍀ by Doppler shifting, which leads to the formation of a trapping plane at the location where ⍀ϭN. As this plane is approached, en- FIG. 13. ͑a͒ Variance of the passive scalar distribution 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. 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. 28 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, 30 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 31 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 inertiagravity waves with a baroclinic barrier, and this is the purpose of current research.
FIG. 14. ͑a͒ Mean velocity profile averaged over x and z, Ū (y,t), at successive times. The curves are plotted every 12 time units, from tϭ0 ͑steepest solid line͒ to tϭ84 ͑smoothest solid line͒, a dashed-dotted line being used for intermediate times. ͑b͒ Squared thickness ␦ 2 (t) of the mean velocity profile, fitted by the function tanh((yϪy 1 )/␦(t)), as a function of time. The slope of the dotted line, equal to 0.08, is the mean value of d(␦ 2 )/d(t) for 33рtр57.