Fluid Mechanics

8.8. Internal Waves in a Continuously Stratified Fluid

Waves may also exist in the interior of a pool, reservoir, lake, or ocean when the fluid’s density in a quiescent state is a continuous function of the vertical coordinate z. The equations of motion for internal waves in such a stratified medium presented here are simplifications of the Boussinesq set specified at the end of Section 4.9. The Boussinesq approximation treats the density as constant, except in the vertical momentum equation. For simplicity, we shall also assume that: 1) the wave motion is effectively inviscid because the velocity gradients are small and the Reynolds number is large, 2) the wave amplitudes are small enough so that the nonlinear advection terms can be neglected, and 3) the frequency of wave motion is much larger than the Coriolis frequency so it does not affect the wave motion. Effects of the earth’s rotation are considered in Chapter 13. The Boussinesq set then simplifies to:

DρDt=0,∂u∂x+∂v∂y+∂w∂z=0,

(4.9, 4.10)

∂u∂t=−1ρ0∂p∂x,∂v∂t=−1ρ0∂p∂y,and∂w∂t=−1ρ0∂p∂z−ρgρ0,

(8.119, 8.120, 8.121)

where ρ₀ is a constant reference density. Here, (4.9) expresses constancy of fluid-particle density while (4.10) is the condition for incompressible flow. If temperature is the only cause of density change, then Dρ/Dt = 0 follows from the heat equation in the non-diffusive form DT/Dt = 0 and a temperature-only equation of state, in the form δρ/ρ = −αδT, where α is the coefficient of thermal expansion. If the density changes are due to changes in the concentration S of a constituent (e.g., salinity in the ocean or water vapor in the atmosphere), then Dρ/Dt = 0 follows from DS/Dt = 0 (the non-diffusive form of the constituent conservation equation) and a concentration-only equation of state, ρ = ρ(S), in the form of δρ/ρ = βδS, where β is the coefficient describing how the density changes due to concentration of the constituent. In both cases, the principle underlying Dρ/Dt = 0 is an equation of state that does not include pressure. In terms of common usage, this equation is frequently called the density equation, as opposed to the continuity equation (4.10).

The five equations (4.9), (4.10), and (8.119) through (8.121) contain five unknowns (u, v, w, p, ρ). Before considering wave motions, first define the quiescent density ρ¯(z)

and pressure p¯(z)

profiles in the medium as those that satisfy a hydrostatic balance:

0=−1ρ0∂p¯∂z−ρ¯gρ0.

(8.122)

When the motion develops, the pressure and density will change relative to their quiescent values:

p=p¯(z)+p′,andρ=ρ¯(z)+ρ′.

(8.123)

The density equation (4.9) then becomes:

∂∂t(ρ¯+ρ′)+u∂∂x(ρ¯+ρ′)+v∂∂y(ρ¯+ρ′)+w∂∂z(ρ¯+ρ′)=0.

(8.124)

Here, ∂ρ¯/∂t=∂ρ¯/∂x=∂ρ¯/∂y=0

, and the nonlinear terms (namely, u∂ρ′/∂x, v∂ρ′/∂y, and w∂ρ′/∂z) are also negligible for small-amplitude motions. The linear part of the fourth term, wdρ¯/dz

, must be retained, so the linearized version of (4.9) is:

∂ρ′∂t+w∂ρ¯∂z=0,

(8.125)

which states that the density perturbation at a point is generated only by the vertical advection of the background density distribution. We now introduce the Brunt–Väisälä frequency, or buoyancy frequency:

N2=−gρ0∂ρ¯∂z.

(8.126)

This is (1.35) when the adiabatic density gradient is zero. As described in Section 1.10, N(z) has units of rad./s and is the oscillation frequency of a vertically displaced fluid particle released from rest in the absence of fluid friction. Using (8.122) and (8.126) in (8.119) through (8.121) and (8.125) produces:

∂u∂t=−1ρ0∂p′∂x,∂v∂t=−1ρ0∂p′∂y,and∂w∂t=−1ρ0∂p′∂z−ρ′gρ0,

(8.127, 8.128, 8.129)

and∂ρ′∂t−N2ρ0gw=0.

(8.130)

Comparing (8.119) through (8.121) and (8.127) through (8.129), we see that the only difference is the replacement of the total density ρ and pressure p with the perturbation density ρ' and pressure p'.

The full set of equations for linear wave motion in a stratified fluid are (4.10) and (8.127) through (8.130), where ρ may be a function of temperature T and concentration S of a constituent, but not of pressure. At first this does not seem to be a good assumption. The compressibility effects in the atmosphere are certainly not negligible; even in the ocean the density changes due to the huge changes in the background pressure are as much as 4%, which is ≈10 times the density changes due to the variations of the salinity and temperature. The effects of compressibility, however, can be handled within the Boussinesq approximation if we regard ρ¯

in the definition of N as the background potential density, that is, the density distribution from which the adiabatic changes of density, due to the changes of pressure, have been subtracted out. The concept of potential density is explained in Chapter 1. Oceanographers account for compressibility effects by converting all their density measurements to the standard atmospheric pressure; thus, when they report variations in density (what they call “sigma tee”) they are generally reporting variations due only to changes in temperature and salinity.

A useful condensation of the above equations involving only w can be obtained by taking the time derivative of (4.10) and using the horizontal momentum equations (8.127) and (8.128) to eliminate u and v. The result is:

1ρ0∇H2p′=∂2w∂z∂t,

(8.131)

where ∇H2=∂2/∂x2+∂2/∂y2

is the horizontal Laplacian operator. Elimination of ρ′ from (8.129) and (8.130) gives:

1ρ0∂2p′∂t∂z=−∂2w∂t2−N2w.

(8.132)

Finally, p′ can be eliminated by taking ∇H2

of (8.132), and inserting the result in (8.131) to find:

∂2∂t∂z(∂2w∂t∂z)=−∇H2(∂2w∂t2+N2w),or∂2∂t2(∇2w)+N2∇H2w=0,

(8.133)

where ∇2=∇H2+∂2/∂z2

is the three-dimensional Laplacian operator. This equation for the vertical velocity w can be used to derive the dispersion relation for internal gravity waves.

Internal Waves in a Stratified Fluid

The situation embodied in (8.133) is fundamentally different from that of interface waves because there is no obvious direction of propagation. For interface waves constrained to follow a horizontal surface with the x-axis chosen along the direction of wave propagation, a dispersion relation ω(k) was obtained that is independent of the wave direction. Furthermore, wave crests and wave groups propagate in the same direction, although at different speeds. However, in the current situation, the fluid is continuously stratified and internal waves might propagate in any direction and at any angle to the vertical. In such a case the direction of the wave number vector K = (k, l, m) becomes important and the dispersion relationship is anisotropic and depends on the wave number components:

ω=ω(k,l,m)=ω(K).

(8.134)

Consequently, the wave number, phase velocity, and group velocity are no longer scalars and the prototype sinusoidal wave (8.2) must be replaced with its three-dimensional extension (8.5). However, (8.134) must still be isotropic in k and l, the wave number components in the two horizontal directions.

The propagation of internal waves is a baroclinic process, in which the surfaces of constant pressure do not coincide with the surfaces of constant density. It was shown in Section 5.2, in connection with Kelvin’s circulation theorem, that baroclinic processes generate vorticity. Internal waves in a continuously stratified fluid are therefore rotational. Waves at a density interface constitute a limiting case in which all the vorticity is concentrated in the form of a velocity discontinuity at the interface. The Laplace equation can therefore be used to describe the flow field within each layer. However, internal waves in a continuously stratified fluid cannot be described by the Laplace equation.

To reveal the structure of the situation described by (8.133) and (8.134), consider the complex version of (8.5) for the vertical velocity:

w=w0ei(kx+ly+mz−ωt)=w0ei(K·x−ωt)

(8.135)

in a fluid medium having a constant buoyancy frequency. Substituting (8.135) into (8.133) with constant N leads to the dispersion relation:

ω2=k2+l2k2+l2+m2N2.

(8.136)

For simplicity choose the x-z plane so it contains K and l = 0. No generality is lost through this choice because the medium is horizontally isotropic, but k now represents the entire horizontal wave number and (8.136) can be written:

ω=kNk2+m2=kNK.

(8.137)

This is the dispersion relation for internal gravity waves and can also be written as:

ω=Ncosθ,

(8.138)

where θ = tan^–1(m/k) is the angle between the phase velocity vector c (and therefore K) and the horizontal direction (Figure 8.29). Interestingly, (8.138) states that the frequency of an internal wave in a stratified fluid depends only on the direction of the wave number vector and not on its magnitude. This is in sharp contrast with surface and interfacial gravity waves, for which frequency depends only on the magnitude. In addition, the wave frequency lies in the range 0 < ω < N, and this indicates that N is the maximum possible frequency of internal waves in a stratified fluid.

Before further investigation of the dispersion relation, consider particle motion in an incompressible internal wave. For consistency with (8.135), the horizontal fluid velocity is written as:

u=u0ei(kx+ly+mz−ωt),

(8.139)

plus two similar expressions for v and w. Differentiating produces:

∂u/∂x=iku0ei(kx+ly+mz−ωt)=iku,

Thus, (4.10) then requires that ku + lv + mw = 0, that is:

K·u=0,

(8.140)

showing that particle motion is perpendicular to the wave number vector (Figure 8.29). Note that only two conditions have been used to derive this result, namely the incompressible continuity equation and trigonometric behavior in all spatial directions. As such, the result is valid for many other wave systems that meet these two conditions. These waves are called shear waves (or transverse waves) because the fluid moves parallel to the constant phase lines. Surface or interfacial gravity waves do not have this property because the field varies exponentially in the vertical.

We can now interpret θ in the dispersion relation (8.138) as the angle between the particle motion and the vertical direction (Figure 8.29). The maximum frequency ω = N occurs when θ = 0, that is, when the particles move up and down vertically. This case corresponds to m = 0 (see (8.137)), showing that the motion is independent of the z-coordinate. The resulting motion consists of a series of vertical columns, all oscillating at the buoyancy frequency N, with the flow field varying in the horizontal direction only.

At the opposite extreme we have ω = 0 when θ = π/2, that is, when the particle motion is completely horizontal. In this limit our internal wave solution (8.137) would seem to require k = 0, that is, horizontal independence of the motion. However, such a conclusion is not valid; pure horizontal motion is not a limiting case of internal waves, and it is necessary to examine the basic equations to draw any conclusion for this case. An examination of the governing set, (4.10) and (8.127) through (8.130), shows that a possible steady solution is w = p′ = ρ′ = 0, with u and v and any functions of x and y satisfying the two-dimensional incompressible continuity equation:

∂u/∂x+∂v/∂y=0.

(7.2)

The z-dependence of u and v is arbitrary. The motion is therefore two dimensional in the horizontal plane, with the motion in the various horizontal planes decoupled from each other. This is why clouds in the upper atmosphere seem to move in flat horizontal sheets, as often observed in airplane flights (Gill, 1982). For a similar reason a cloud pattern pierced by a mountain peak sometimes shows Karman vortex streets, a two-dimensional flow feature; see the striking photograph in Figure 10.21. A restriction of strong stratification is necessary for such almost horizontal flows, because (8.130) suggests that the vertical motion is small if N is large.

The foregoing discussion leads to the interesting phenomenon of blocking in a strongly stratified fluid. Consider a two-dimensional body placed in such a fluid, with its axis horizontal (Figure 8.30). The two dimensionality of the body requires ∂v/∂y = 0, so that the continuity equation (7.2) reduces to ∂u/∂x = 0. A horizontal layer of fluid ahead of the body, bounded by tangents above and below it, is therefore blocked and held motionless. (For photographic evidence see Figure 3.18 in the book by Turner (1973).) This happens because the strong stratification suppresses the w field and prevents the fluid from going below or over the body.

Dispersion of Internal Waves in a Stratified Fluid

The dispersion relationship (8.137) for linear internal waves with constant buoyancy frequency contains a few genuine surprises that challenge our imaginations and violate the intuition acquired by observing surface or interface waves. One of these surprises involves the phase, c, and group, c_g, velocity vectors. In multiple dimensions, these are defined by:

c=(ω/K)eKandcg=ex(∂ω/∂k)+ey(∂ω/∂l)+ez(∂ω/∂m),

(8.8, 8.141)

where e_K = K/K. For interface waves c and c_g are in the same direction, although their magnitudes can be different. For internal waves, (8.137), (8.8), and (8.141) can be used to determine:

c=ωK2(kex+mez),andcg=NmK3(mex−kez).

(8.142, 8.143)

Forming the dot product of these two equations produces:

cg·c=0!

(8.144)

Thus, the phase and group velocity vectors are perpendicular as shown on Figure 8.29. Equations (8.142) and (8.143) do place the horizontal components of c and c_g in the same direction, but their vertical components are equal and opposite. In fact, c and c_g form two sides of a right triangle whose hypotenuse is horizontal (Figure 8.31). Consequently, the phase velocity has an upward component when the group velocity has a downward component, and vice versa. Equations (8.140) and (8.144) are consistent because c and K are parallel and c_g and u are parallel. The fact that c and c_g are perpendicular, and have opposite vertical components, is illustrated in Figure 8.32. It shows that the phase lines are propagating toward the left and upward, whereas the wave group is propagating to the left and downward. Wave crests are constantly appearing at one edge of the group, propagating through the group, and vanishing at the other edge.

The group velocity here has the usual significance of being the velocity of propagation of energy of a certain sinusoidal component. Suppose a source is oscillating at frequency ω. Then its energy will only be transmitted outward along four beams oriented at an angle θ with the vertical, where cosθ = ω/N. This has been verified in a laboratory experiment (Figure 8.33). The source in this case was a vertically oscillating cylinder with its axis perpendicular to the plane of paper. The frequency was ω < N. The light and dark lines in the photograph are lines of constant density, made visible by an optical technique. The experiment showed that the energy radiated along four beams that became more vertical as the frequency was increased, which agrees with cosθ = ω/N.

These results were obtained by assuming that N is depth independent, an assumption that may seem unrealistic at first. Figure 13.2 shows N vs. depth for the deep ocean, and N < 0.01 everywhere, but N is largest between ∼200 m and ∼2 km. These results can be considered locally valid if N varies slowly over the vertical wavelength 2π/m of the motion. The so-called WKB approximation for internal waves, in which such a slow variation of N(z) is not neglected, is discussed in Chapter 13.

Energy Considerations for Internal Waves in a Stratified Fluid

The energy carried by an internal wave travels in the direction and at the speed of the group velocity. To show this is the case, construct a mechanical energy equation from (8.127) through (8.129) by multiplying the first equation by ρ₀u, the second by ρ₀v, the third by ρ₀w, and summing the results to find:

∂∂t[ρ02(u2+v2+w2)]+gρ′w+∇·(p′u)=0.

(8.145)

Here, the continuity equation has been used to write u(∂p′/∂x) + v(∂p′/∂y) + w(∂p′/∂z) = ∇·(p′u)

, which represents the net work done by pressure forces. Another interpretation is that ∇·(p′u)

is the divergence of the energy flux p′u, which must change the wave energy at a point. As the first term in (8.145) is the rate of change of kinetic energy, we can anticipate that the second term gρ′w must be the rate of change of potential energy, E_p. This is consistent with the energy principle derived in Chapter 4 (see (4.56)), except that ρ′ and p′ replace ρ and p because we have subtracted the mean quiescent state here. Using the density equation (8.130), the rate of change of potential energy can be written as:

∂Ep∂t=gρ′w=∂∂t[g2ρ′22ρ0N2],

(8.146)

which shows that the potential energy per unit volume must be the positive quantity E_p = g²ρ′²/2ρ₀N². The potential energy can also be expressed in terms of the displacement ζ of a fluid particle, given by w = ∂ζ/∂t. Using (8.130), we can write:

∂ρ′∂t=N2ρ0g∂ζ∂t,whichrequiresthatρ′=N2ρ0ζg.

(8.147)

The potential energy per unit volume is therefore:

Ep=g2ρ′22ρ0N2=12N2ρ0ζ2.

(8.148)

This expression is consistent with our previous result from (8.96) for two infinitely deep fluids, for which the average potential energy of the entire water column per unit horizontal area was shown to be:

14(ρ2−ρ1)ga2,

(8.149)

where the interface displacement is of the form ζ = acos(kx − ωt) and (ρ₂ − ρ₁) is the density discontinuity. To see the consistency, we shall symbolically represent the buoyancy frequency of a density discontinuity at z = 0 as:

N2=−gρ0dρ¯dz=gρ0(ρ2−ρ1)δ(z),

(8.150)

where δ(z) is the Dirac delta function (see Appendix B.4). (As with other relations involving the delta function, (8.150) is valid in the integral sense, that is, the integral (across z = 0) of the last two terms is equal because ∫δ(z)dz = 1.) Using (8.150), a vertical integral of (8.148), coupled with horizontal averaging over a wavelength, gives the expression (8.149). Note that for surface or interfacial waves, E_k and E_p represent kinetic and potential energies of the entire water column, per unit horizontal area. In a continuously stratified fluid, they represent energies per unit volume.

We shall now demonstrate that the average kinetic and potential energies are equal for internal wave motion. Assume periodic solutions:

[u,w,p′,ρ′]=[uˆ,wˆ,pˆ,ρˆ]ei(kx+ly+mz−ωt).

(8.151)

Then all variables can be expressed in terms of w:

p′=−ωmρ0k2wˆei(kx+mz−ωt),ρ′=iN2ρ0ωgwˆei(kx+mz−ωt),andu=−mkwˆei(kx+mz−ωt),

(8.152)

where p′ is derived from (8.131), ρ′ from (8.130), and u from (8.127). The average kinetic energy per unit volume is therefore:

Ek=12ρ0(u2+w2)¯=14ρ0(m2k2+1)wˆ2,

(8.153)

where we have taken real parts of the various expressions in (8.151) before computing quadratic quantities and used the fact that the average of cos²() over a wavelength is 1/2. The average potential energy per unit volume is:

Ep=g2ρ′22ρ0N2=N2ρ04ω2wˆ2,

(8.154)

where we have used ρ′2¯=wˆ2N4ρ02/2ω2g2

, found from (8.152) after taking its real part. Use of the dispersion relation ω² = k²N²/(k² + m²) shows that:

Ek=Ep,

(8.155)

which is a general result for small oscillations of a conservative system without Coriolis forces. The total wave energy is:

E=Ek+Ep=12ρ0(m2k2+1)wˆ2.

(8.156)

Last, we shall show that c_g times the wave energy equals the energy flux. The average energy flux F across a unit area can be found from (8.151):

F=p′u¯=exp′u¯+ezp′w¯=ρ0ωmwˆ22k2(exmk−ez).

(8.157)

Using equations (8.143) and (8.156), group velocity times wave energy is:

cgE=NmK3(mex−kez)[12ρ0(m2k2+1)wˆ2],

which reduces to (8.157) on using the dispersion relation (8.137), so it follows that:

F=cgE.

(8.158)

This result also holds for surface or interfacial gravity waves. However, in that case F represents the flux per unit width perpendicular to the propagation direction (integrated over the entire depth), and E represents the energy per unit horizontal area. In (8.157), on the other hand, F is the flux per unit area, and E is the energy per unit volume.

Example 8.8

If the wave's frequency is ω ≤ N and the phase speed vector makes an angle θ with the horizontal direction as shown in Figure 8.29, what are the fluid particle trajectories in the x-z plane as an internal wave with vertical velocity amplitude w_o passes?

Solution

In this case, the horizontal and vertical velocities are:

u=−mkwocos(kx+mz−ωt)andw=wocos(kx+mz−ωt).

as follows from (8.151) and (8.152) after real parts are taken. Fluid particle paths are given by x and z components of (3.8):

dxpdt=−mkwocos(kxp+mzp−ωt)anddzpdt=wocos(kxp+mzp−ωt),

where the subscript p indicates the fluid particle's coordinates. The ratio of these two equations leads to: dz_p/dx_p = –k/m, or kx_p + mz_p = constant = kx_o + mz_o, if the fluid particle of interest is at x_o and z_o at t = 0, then:

k(xp−xo)=−m(zp−zo)=mwoω(cos(kxo+mzo−ωt)−cos(kxo+mzo)).

Thus, this fluid particle oscillates along a line segment centered on (x_o,z_o) with horizontal amplitude mw_o/kω = (w_o/ω)tanθ, and vertical amplitude w_o/ω. The orientation of this line segment is perpendicular to the wave number vector K = (k, 0, m). These fluid particle trajectories are different from the circles or ellipses found in the water column below ordinary water-surface waves; see (8.35).